- added get slip rate helper function

- added vacancy generation rate
- cleaned up some bugs

code/constitutive.f90

@@ -42,7 +42,8 @@ module constitutive
  
  private :: &
    constitutive_hooke_TandItsTangent, &
-   constitutive_getAccumulatedSlip
+   constitutive_getAccumulatedSlip, &
+   constitutive_getSlipRate
  
 contains
 
@@ -784,6 +785,10 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
    tick, tock, & 
    tickrate, &
    maxticks
+ real(pReal), dimension(:), allocatable :: &
+   accumulatedSlip
+ integer(pInt) :: &
+   nSlip
  
  if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) &
    call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
@@ -808,7 +813,7 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
    case (LOCAL_DAMAGE_brittle_ID)
      call damage_brittle_dotState(ipc, ip, el)
    case (LOCAL_DAMAGE_ductile_ID)
-     call damage_ductile_dotState(Lp, ipc, ip, el)
+     call damage_ductile_dotState(ipc, ip, el)
    case (LOCAL_DAMAGE_gurson_ID)
      call damage_gurson_dotState(Lp, ipc, ip, el)
  end select
@@ -820,7 +825,9 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
 
  select case (phase_vacancy(material_phase(ipc,ip,el)))
    case (LOCAL_VACANCY_generation_ID)
-     call vacancy_generation_dotState(Tstar_v, Lp, ipc, ip, el)
+     call constitutive_getAccumulatedSlip(nSlip,accumulatedSlip,FpArray(1:3,1:3,ipc,ip,el),ipc,ip,el)
+     call vacancy_generation_dotState(nSlip,accumulatedSlip,Tstar_v,constitutive_getTemperature(ipc,ip,el), &
+                                      ipc, ip, el)
  end select
 
  if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) then
@@ -1217,7 +1224,7 @@ subroutine constitutive_getAccumulatedSlip(nSlip,accumulatedSlip,Fp,ipc, ip, el)
    pInt
  use math, only: &
    math_mul33xx33, &
-   math_norm33, &
+   math_equivStrain33, &
    math_I3
  use material, only: &
    phase_plasticity, &
@@ -1260,7 +1267,7 @@ subroutine constitutive_getAccumulatedSlip(nSlip,accumulatedSlip,Fp,ipc, ip, el)
    case (PLASTICITY_J2_ID)
      nSlip = 1_pInt
      allocate(accumulatedSlip(nSlip))
-     accumulatedSlip(1) = math_norm33(math_mul33xx33(transpose(Fp),Fp) - math_I3)/2.0_pReal
+     accumulatedSlip(1) = math_equivStrain33((math_mul33xx33(transpose(Fp),Fp) - math_I3)/2.0_pReal)
    case (PLASTICITY_PHENOPOWERLAW_ID)
      call constitutive_phenopowerlaw_getAccumulatedSlip(nSlip,accumulatedSlip,ipc, ip, el)
    case (PLASTICITY_DISLOTWIN_ID)
@@ -1275,6 +1282,73 @@ subroutine constitutive_getAccumulatedSlip(nSlip,accumulatedSlip,Fp,ipc, ip, el)
 
 end subroutine constitutive_getAccumulatedSlip
 
+!--------------------------------------------------------------------------------------------------
+!> @brief returns accumulated slip rates on each system defined
+!--------------------------------------------------------------------------------------------------
+subroutine constitutive_getSlipRate(nSlip,slipRate,Lp,ipc, ip, el)
+ use prec, only: &
+   pReal, &
+   pInt
+ use math, only: &
+   math_mul33xx33, &
+   math_equivStrain33, &
+   math_I3
+ use material, only: &
+   phase_plasticity, &
+   material_phase, &
+   PLASTICITY_none_ID, &
+   PLASTICITY_j2_ID, &
+   PLASTICITY_phenopowerlaw_ID, &
+   PLASTICITY_dislotwin_ID, &
+   PLASTICITY_dislokmc_ID, &
+   PLASTICITY_titanmod_ID, &
+   PLASTICITY_nonlocal_ID
+ use constitutive_phenopowerlaw, only: &
+   constitutive_phenopowerlaw_getSlipRate
+ use constitutive_dislotwin, only: &
+   constitutive_dislotwin_getSlipRate
+ use constitutive_dislokmc, only: &
+   constitutive_dislokmc_getSlipRate
+ use constitutive_titanmod, only: &
+   constitutive_titanmod_getSlipRate
+ use constitutive_nonlocal, only: &
+   constitutive_nonlocal_getSlipRate
+
+ implicit none
+ 
+ real(pReal), dimension(:), allocatable :: &
+   slipRate
+ integer(pInt) :: &
+   nSlip
+ real(pReal),  intent(in), dimension(3,3) :: &
+   Lp                                                                                               !< plastic velocity gradient
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ 
+ select case (phase_plasticity(material_phase(ipc,ip,el)))
+   case (PLASTICITY_none_ID)
+     nSlip = 0_pInt
+     allocate(slipRate(nSlip))
+   case (PLASTICITY_J2_ID)
+     nSlip = 1_pInt
+     allocate(slipRate(nSlip))
+     slipRate(1) = math_equivStrain33(Lp)
+   case (PLASTICITY_PHENOPOWERLAW_ID)
+     call constitutive_phenopowerlaw_getSlipRate(nSlip,slipRate,ipc, ip, el)
+   case (PLASTICITY_DISLOTWIN_ID)
+     call constitutive_dislotwin_getSlipRate(nSlip,slipRate,ipc, ip, el)
+   case (PLASTICITY_DISLOKMC_ID)
+     call constitutive_dislokmc_getSlipRate(nSlip,slipRate,ipc, ip, el)
+   case (PLASTICITY_TITANMOD_ID)
+     call constitutive_titanmod_getSlipRate(nSlip,slipRate,ipc, ip, el)
+   case (PLASTICITY_NONLOCAL_ID)
+     call constitutive_nonlocal_getSlipRate(nSlip,slipRate,ipc, ip, el)
+ end select
+
+end subroutine constitutive_getSlipRate
+
 !--------------------------------------------------------------------------------------------------
 !> @brief returns array of constitutive results
 !--------------------------------------------------------------------------------------------------

code/constitutive_dislokmc.f90

@@ -153,6 +153,7 @@ module constitutive_dislokmc
    constitutive_dislokmc_LpAndItsTangent, &
    constitutive_dislokmc_dotState, &
    constitutive_dislokmc_getAccumulatedSlip, &
+   constitutive_dislokmc_getSlipRate, &
    constitutive_dislokmc_postResults
  private :: &
    constitutive_dislokmc_stateInit, &
@@ -1674,6 +1675,7 @@ subroutine constitutive_dislokmc_dotState(Tstar_v,Temperature,ipc,ip,el)
  enddo
  
 end subroutine constitutive_dislokmc_dotState
+
 !--------------------------------------------------------------------------------------------------
 !> @brief returns accumulated slip
 !--------------------------------------------------------------------------------------------------
@@ -1719,6 +1721,50 @@ subroutine constitutive_dislokmc_getAccumulatedSlip(nSlip,accumulatedSlip,ipc, i
    
 end subroutine constitutive_dislokmc_getAccumulatedSlip
 
+!--------------------------------------------------------------------------------------------------
+!> @brief returns accumulated slip
+!--------------------------------------------------------------------------------------------------
+subroutine constitutive_dislokmc_getSlipRate(nSlip,slipRate,ipc, ip, el)
+ use lattice, only: &
+   lattice_maxNslipFamily
+ use material, only: &
+   mappingConstitutive, &
+   plasticState, &
+   phase_plasticityInstance
+
+   implicit none
+ 
+ real(pReal), dimension(:), allocatable :: &
+   slipRate
+ integer(pInt) :: &
+   nSlip
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ integer(pInt) :: &
+   offset, &
+   phase, &
+   instance, &
+   offset_accshear_slip, &
+   f, j, i
+
+ offset = mappingConstitutive(1,ipc,ip,el)
+ phase = mappingConstitutive(2,ipc,ip,el)
+ instance = phase_plasticityInstance(phase)
+ nSlip = constitutive_dislokmc_totalNslip(instance)
+ allocate(slipRate(nSlip))
+ offset_accshear_slip = 2_pInt*nSlip
+
+ j = 0_pInt
+ do f = 1_pInt,lattice_maxNslipFamily                                 ! loop over all slip families
+   do i = 1_pInt,constitutive_dislokmc_Nslip(f,instance)             ! process each (active) slip system in family
+      j = j+1_pInt
+      slipRate(j) = plasticState(phase)%dotState(offset_accshear_slip+j,offset)
+   enddo
+ enddo 
+   
+end subroutine constitutive_dislokmc_getSlipRate
 
  
 !--------------------------------------------------------------------------------------------------

code/constitutive_dislotwin.f90

@@ -167,6 +167,7 @@ module constitutive_dislotwin
    constitutive_dislotwin_LpAndItsTangent, &
    constitutive_dislotwin_dotState, &
    constitutive_dislotwin_getAccumulatedSlip, &
+   constitutive_dislotwin_getSlipRate, &
    constitutive_dislotwin_postResults
  private :: &
    constitutive_dislotwin_stateInit, &
@@ -1942,6 +1943,52 @@ subroutine constitutive_dislotwin_getAccumulatedSlip(nSlip,accumulatedSlip,ipc,
 end subroutine constitutive_dislotwin_getAccumulatedSlip
 
  
+!--------------------------------------------------------------------------------------------------
+!> @brief returns accumulated slip rate
+!--------------------------------------------------------------------------------------------------
+subroutine constitutive_dislotwin_getSlipRate(nSlip,slipRate,ipc, ip, el)
+ use lattice, only: &
+   lattice_maxNslipFamily
+ use material, only: &
+   mappingConstitutive, &
+   plasticState, &
+   phase_plasticityInstance
+
+   implicit none
+ 
+ real(pReal), dimension(:), allocatable :: &
+   slipRate
+ integer(pInt) :: &
+   nSlip
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ integer(pInt) :: &
+   offset, &
+   phase, &
+   instance, &
+   offset_accshear_slip, &
+   f, j, i
+
+ offset = mappingConstitutive(1,ipc,ip,el)
+ phase = mappingConstitutive(2,ipc,ip,el)
+ instance = phase_plasticityInstance(phase)
+ nSlip = constitutive_dislotwin_totalNslip(instance)
+ allocate(slipRate(nSlip))
+ offset_accshear_slip = 2_pInt*nSlip
+
+ j = 0_pInt
+ do f = 1_pInt,lattice_maxNslipFamily                                 ! loop over all slip families
+   do i = 1_pInt,constitutive_dislotwin_Nslip(f,instance)             ! process each (active) slip system in family
+      j = j+1_pInt
+      slipRate(j) = plasticState(phase)%dotState(offset_accshear_slip+j,offset)
+   enddo
+ enddo 
+   
+end subroutine constitutive_dislotwin_getSlipRate
+
+ 
 !--------------------------------------------------------------------------------------------------
 !> @brief return array of constitutive results
 !--------------------------------------------------------------------------------------------------

code/constitutive_nonlocal.f90

@@ -251,6 +251,7 @@ module constitutive_nonlocal
  constitutive_nonlocal_deltaState, &
  constitutive_nonlocal_updateCompatibility, &
  constitutive_nonlocal_getAccumulatedSlip, &
+ constitutive_nonlocal_getSlipRate, &
  constitutive_nonlocal_postResults
  
  private :: &
@@ -3579,11 +3580,49 @@ subroutine constitutive_nonlocal_getAccumulatedSlip(nSlip,accumulatedSlip,ipc, i
  nSlip = totalNslip(instance)
  allocate(accumulatedSlip(nSlip))
  forall (s = 1:nSlip) &
-   accumulatedSlip(s) = plasticState(phase)%dotState(iGamma(s,instance),offset)
+   accumulatedSlip(s) = plasticState(phase)%state(iGamma(s,instance),offset)
    
 end subroutine constitutive_nonlocal_getAccumulatedSlip
 
  
+!--------------------------------------------------------------------------------------------------
+!> @brief returns accumulated slip rate
+!--------------------------------------------------------------------------------------------------
+subroutine constitutive_nonlocal_getSlipRate(nSlip,slipRate,ipc, ip, el)
+ use lattice, only: &
+   lattice_maxNslipFamily
+ use material, only: &
+   mappingConstitutive, &
+   plasticState, &
+   phase_plasticityInstance
+
+   implicit none
+ 
+ real(pReal), dimension(:), allocatable :: &
+   slipRate
+ integer(pInt) :: &
+   nSlip
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ integer(pInt) :: &
+   offset, &
+   phase, &
+   instance, &
+   s
+
+ offset = mappingConstitutive(1,ipc,ip,el)
+ phase = mappingConstitutive(2,ipc,ip,el)
+ instance = phase_plasticityInstance(phase) 
+ nSlip = totalNslip(instance)
+ allocate(slipRate(nSlip))
+ forall (s = 1:nSlip) &
+   slipRate(s) = plasticState(phase)%dotState(iGamma(s,instance),offset)
+   
+end subroutine constitutive_nonlocal_getSlipRate
+
+ 
 !--------------------------------------------------------------------------------------------------
 !> @brief return array of constitutive results
 !--------------------------------------------------------------------------------------------------

code/constitutive_phenopowerlaw.f90

@@ -90,6 +90,7 @@ module constitutive_phenopowerlaw
    constitutive_phenopowerlaw_LpAndItsTangent, &
    constitutive_phenopowerlaw_dotState, &
    constitutive_phenopowerlaw_getAccumulatedSlip, &
+   constitutive_phenopowerlaw_getSlipRate, &
    constitutive_phenopowerlaw_postResults
  private :: &
    constitutive_phenopowerlaw_aTolState, &
@@ -1039,6 +1040,54 @@ subroutine constitutive_phenopowerlaw_getAccumulatedSlip(nSlip,accumulatedSlip,i
 end subroutine constitutive_phenopowerlaw_getAccumulatedSlip
 
  
+!--------------------------------------------------------------------------------------------------
+!> @brief returns accumulated slip rate
+!--------------------------------------------------------------------------------------------------
+subroutine constitutive_phenopowerlaw_getSlipRate(nSlip,slipRate,ipc, ip, el)
+ use lattice, only: &
+   lattice_maxNslipFamily
+ use material, only: &
+   mappingConstitutive, &
+   plasticState, &
+   phase_plasticityInstance
+
+   implicit none
+ 
+ real(pReal), dimension(:), allocatable :: &
+   slipRate
+ integer(pInt) :: &
+   nSlip
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ integer(pInt) :: &
+   offset, &
+   phase, &
+   instance, &
+   offset_accshear_slip, &
+   nTwin, &
+   f, j, i
+
+ offset = mappingConstitutive(1,ipc,ip,el)
+ phase = mappingConstitutive(2,ipc,ip,el)
+ instance = phase_plasticityInstance(phase)
+ nSlip = constitutive_phenopowerlaw_totalNslip(instance)
+ nTwin = constitutive_phenopowerlaw_totalNtwin(instance)
+ offset_accshear_slip = nSlip + nTwin + 2_pInt
+ 
+ allocate(slipRate(nSlip))
+ j = 0_pInt
+ slipFamiliesLoop: do f = 1_pInt,lattice_maxNslipFamily
+   do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance)                                       ! process each (active) slip system in family
+     j = j+1_pInt
+     slipRate(j) = plasticState(phase)%dotState(offset_accshear_slip+j,offset)
+   enddo
+ enddo slipFamiliesLoop
+   
+end subroutine constitutive_phenopowerlaw_getSlipRate
+
+ 
 !--------------------------------------------------------------------------------------------------
 !> @brief return array of constitutive results
 !--------------------------------------------------------------------------------------------------

code/constitutive_titanmod.f90

@@ -177,6 +177,7 @@ module constitutive_titanmod
    constitutive_titanmod_LpAndItsTangent, &
    constitutive_titanmod_dotState, &
    constitutive_titanmod_getAccumulatedSlip, &
+   constitutive_titanmod_getSlipRate, &
    constitutive_titanmod_postResults, &
    constitutive_titanmod_homogenizedC
 
@@ -1822,6 +1823,50 @@ subroutine constitutive_titanmod_getAccumulatedSlip(nSlip,accumulatedSlip,ipc, i
 end subroutine constitutive_titanmod_getAccumulatedSlip
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief returns accumulated slip
+!--------------------------------------------------------------------------------------------------
+subroutine constitutive_titanmod_getSlipRate(nSlip,slipRate,ipc, ip, el)
+ use lattice, only: &
+   lattice_maxNslipFamily
+ use material, only: &
+   mappingConstitutive, &
+   plasticState, &
+   phase_plasticityInstance
+
+   implicit none
+ 
+ real(pReal), dimension(:), allocatable :: &
+   slipRate
+ integer(pInt) :: &
+   nSlip
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ integer(pInt) :: &
+   offset, &
+   phase, &
+   instance, &
+   offset_accshear_slip, &
+   f, j, i
+
+ offset = mappingConstitutive(1,ipc,ip,el)
+ phase = mappingConstitutive(2,ipc,ip,el)
+ instance = phase_plasticityInstance(phase)
+ nSlip = constitutive_titanmod_totalNslip(instance)
+ allocate(slipRate(nSlip))
+ offset_accshear_slip = 2_pInt*nSlip
+
+ j = 0_pInt
+ do f = 1_pInt,lattice_maxNslipFamily                                 ! loop over all slip families
+   do i = 1_pInt,constitutive_titanmod_Nslip(f,instance)             ! process each (active) slip system in family
+      j = j+1_pInt
+      slipRate(j) = plasticState(phase)%dotState(offset_accshear_slip+j,offset)
+   enddo
+ enddo 
+   
+end subroutine constitutive_titanmod_getSlipRate
 
 
 !--------------------------------------------------------------------------------------------------

code/damage_ductile.f90

@@ -259,7 +259,7 @@ end subroutine damage_ductile_aTolState
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
-subroutine damage_ductile_dotState(Lp, ipc, ip, el)
+subroutine damage_ductile_dotState(ipc, ip, el)
  use material, only: &
    mappingConstitutive, &
    damageState
@@ -269,8 +269,6 @@ subroutine damage_ductile_dotState(Lp, ipc, ip, el)
    lattice_DamageMobility
 
  implicit none
- real(pReal), intent(in), dimension(3,3) :: &
-   Lp
  integer(pInt), intent(in) :: &
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point

code/vacancy_generation.f90

@@ -25,7 +25,14 @@ module vacancy_generation
    vacancy_generation_Noutput                                                                   !< number of outputs per instance of this damage 
 
  real(pReal),                         dimension(:),           allocatable,         public :: &
-   vacancy_generation_aTol
+   vacancy_generation_aTol, &
+   vacancy_generation_freq, &
+   vacancy_generation_energy, &
+   vacancy_generation_C1, &
+   vacancy_generation_C2
+
+ real(pReal),                                                 parameter,           private :: &
+   kB = 1.38e-23_pReal                                                                              !< Boltzmann constant in J/Kelvin
 
  enum, bind(c) 
    enumerator :: undefined_ID, &
@@ -118,6 +125,10 @@ subroutine vacancy_generation_init(fileUnit)
  allocate(vacancy_generation_outputID(maxval(phase_Noutput),maxNinstance),      source=undefined_ID)
  allocate(vacancy_generation_Noutput(maxNinstance),                             source=0_pInt) 
  allocate(vacancy_generation_aTol(maxNinstance),                                source=0.0_pReal) 
+ allocate(vacancy_generation_freq(maxNinstance),                                source=0.0_pReal) 
+ allocate(vacancy_generation_energy(maxNinstance),                              source=0.0_pReal) 
+ allocate(vacancy_generation_C1(maxNinstance),                                  source=0.0_pReal) 
+ allocate(vacancy_generation_C2(maxNinstance),                                  source=0.0_pReal) 
 
  rewind(fileUnit)
  phase = 0_pInt
@@ -155,6 +166,18 @@ subroutine vacancy_generation_init(fileUnit)
        case ('atol_vacancyGeneration')
          vacancy_generation_aTol(instance) = IO_floatValue(line,positions,2_pInt)
 
+       case ('vacancy_frequency')
+         vacancy_generation_freq(instance) = IO_floatValue(line,positions,2_pInt)
+
+       case ('vacancy_energy')
+         vacancy_generation_energy(instance) = IO_floatValue(line,positions,2_pInt)
+
+       case ('vacancy_C1')
+         vacancy_generation_C1(instance) = IO_floatValue(line,positions,2_pInt)
+
+       case ('vacancy_C2')
+         vacancy_generation_C2(instance) = IO_floatValue(line,positions,2_pInt)
+
      end select
    endif; endif
  enddo parsingFile
@@ -250,7 +273,7 @@ end subroutine vacancy_generation_aTolState
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
-subroutine vacancy_generation_dotState(Tstar_v, Lp, ipc, ip, el)
+subroutine vacancy_generation_dotState(nSlip, accumulatedSlip, Tstar_v, Temperature, ipc, ip, el)
  use lattice, only: &
    lattice_massDensity, &
    lattice_specificHeat
@@ -259,28 +282,36 @@ subroutine vacancy_generation_dotState(Tstar_v, Lp, ipc, ip, el)
    phase_vacancyInstance, &
    vacancyState
  use math, only: &
-   math_Mandel6to33
+   math_Mandel6to33, &
+   math_trace33
 
  implicit none
  integer(pInt), intent(in) :: &
+   nSlip, &
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point
    el                                                                                               !< element
+ real(pReal), dimension(nSlip), intent(in) :: &
+   accumulatedSlip
  real(pReal),  intent(in), dimension(6) :: &
    Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor (Mandel)
- real(pReal),  intent(in), dimension(3,3) :: &
+ real(pReal),  intent(in) :: &
-   Lp
+   Temperature                                                                                          !< 2nd Piola Kirchhoff stress tensor (Mandel)
+ real(pReal) :: &
+   pressure                                                                                          !< 2nd Piola Kirchhoff stress tensor (Mandel)
  integer(pInt) :: &
    instance, phase, constituent 
 
  phase = mappingConstitutive(2,ipc,ip,el)
  constituent = mappingConstitutive(1,ipc,ip,el)
  instance = phase_vacancyInstance(phase)
+ pressure = math_trace33(math_Mandel6to33(Tstar_v))
  
  vacancyState(phase)%dotState(1,constituent) = &
-    0.95_pReal &
+    vacancy_generation_freq(instance)* &
-  * sum(abs(math_Mandel6to33(Tstar_v)*Lp)) &
+    (1.0_pReal + vacancy_generation_C2(instance)*sum(accumulatedSlip))* &     
-  / (lattice_massDensity(phase)*lattice_specificHeat(phase))
+    exp(-(vacancy_generation_energy(instance) - vacancy_generation_C2(instance)*pressure)/ &
+         (kB*Temperature))
   
 end subroutine vacancy_generation_dotState