Added specularity factor term to ph-thin-film term.

src/bte.f90

@@ -485,8 +485,9 @@ subroutine dragless_phbte_RTA(Tdir, self, num, crys, sym, ph, el)
          self%ph_rta_rates_bound_ibz + self%ph_rta_rates_4ph_ibz
 
     ! phonon-thin-film
-    call calculate_thinfilm_scatt_rates(ph%prefix, num%phthinfilm, num%phthinfilm_ballistic, crys%thinfilm_height, &
-         crys%thinfilm_normal, ph%vels, ph%indexlist_irred, self%ph_rta_rates_ibz, self%ph_rta_rates_thinfilm_ibz)
+    call calculate_thinfilm_scatt_rates(ph%prefix, num%phthinfilm, num%phthinfilm_ballistic, &
+         crys%specfac, crys%thinfilm_height, crys%thinfilm_normal, &
+         ph%vels, ph%indexlist_irred, self%ph_rta_rates_ibz, self%ph_rta_rates_thinfilm_ibz)
 
     !Matthiessen's rule with thin-film
     self%ph_rta_rates_ibz = self%ph_rta_rates_ibz + self%ph_rta_rates_thinfilm_ibz

src/crystal.f90

@@ -100,6 +100,8 @@ module crystal_module
      !! Characteristic boundary scattering length in mm
      real(r64) :: thinfilm_height
      !! Height of thin-film in mm
+     real(r64) :: specfac
+     !! Specularity factor
      character(1) :: thinfilm_normal
      !! Normal direction of the thin-film: 'x', 'y', or 'z'.
 
@@ -121,7 +123,7 @@ subroutine read_input_and_setup_crystal(self)
     real(r64), allocatable :: masses(:), born(:,:,:), basis(:,:), &
          basis_cart(:,:), subs_perc(:), subs_masses(:), subs_conc(:), &
          dopant_masses(:, :), dopant_conc(:, :)
-    real(r64) :: epsilon(3,3), lattvecs(3,3), T, &
+    real(r64) :: epsilon(3,3), lattvecs(3,3), T, specfac, &
          epsilon0, epsiloninf, subs_mavg, bound_length, thinfilm_height
     character(len=3), allocatable :: elements(:)
     character(len=100) :: name
@@ -133,7 +135,7 @@ subroutine read_input_and_setup_crystal(self)
          polar, born, epsilon, read_epsiloninf, epsilon0, epsiloninf, &
          masses, T, VCA, DIB, twod, subs_masses, subs_conc, bound_length, &
          numdopants_types, dopant_masses, dopant_conc, thinfilm_height, &
-         thinfilm_normal
+         thinfilm_normal, specfac
 
     call subtitle("Setting up crystal...")
 
@@ -187,6 +189,7 @@ subroutine read_input_and_setup_crystal(self)
     bound_length = 1.e12_r64 !mm, practically inifinity
     thinfilm_height = 1.e12_r64 !mm, practically inifinity
     thinfilm_normal = 'z'
+    specfac = 0.0_r64 !Diffusive scattering by default
     numdopants_types = [1, 1]
     dopant_masses = 0.0_r64
     dopant_conc = 0.0_r64
@@ -236,6 +239,7 @@ subroutine read_input_and_setup_crystal(self)
     self%bound_length = bound_length
     self%thinfilm_height = thinfilm_height
     self%thinfilm_normal = thinfilm_normal
+    self%specfac = specfac
     self%numdopants_types = numdopants_types
 
     if(product(numdopants_types) <= 0) then

src/interactions.f90

@@ -3391,13 +3391,14 @@ subroutine calculate_bound_scatt_rates(prefix, finite_crys, length, vels_fbz, &
     call write2file_rank2_real(prefix // '.W_rta_'//prefix//'bound', scatt_rates)
   end subroutine calculate_bound_scatt_rates
 
-  subroutine calculate_thinfilm_scatt_rates(prefix, finite_crys, ballistic_limit, &
+  subroutine calculate_thinfilm_scatt_rates(prefix, finite_crys, ballistic_limit, specfac, &
        height, normal, vels_fbz, indexlist_irred, other_scatt_rates, thin_film_scatt_rates)
     !! Subroutine to calculate the phonon/electron-thin-film scattering rates.
     !!
     !! prefix Type of particle
     !! finite_crys Is the crystal finite?
     !! ballistic_limit Use ballistic limit of Fuchs-Sondheimer theory?
+    !! specfun Specularity factor
     !! height Height of thin-film in mm
     !! normal Normal direction to thin-film
     !! vels Velocities on the FBZ
@@ -3408,45 +3409,47 @@ subroutine calculate_thinfilm_scatt_rates(prefix, finite_crys, ballistic_limit,
     character(len = 2), intent(in) :: prefix
     logical, intent(in) :: finite_crys
     logical, intent(in) :: ballistic_limit
+    real(r64), intent(in) :: specfac
     real(r64), intent(in) :: height
     character(1), intent(in) :: normal
-    real(r64), intent(in) :: vels_fbz(:,:,:)
+    real(r64), intent(in) :: vels_fbz(:, :, :)
     integer(i64), intent(in) :: indexlist_irred(:)
-    real(r64), allocatable, intent(in) :: other_scatt_rates(:,:)
-    real(r64), allocatable, intent(out) :: thin_film_scatt_rates(:,:)
+    real(r64), allocatable, intent(in) :: other_scatt_rates(:, :)
+    real(r64), allocatable, intent(out) :: thin_film_scatt_rates(:, :)
 
     !Local variables
-    integer(i64) :: ik, ib, nk_irred, nb, dir
+    integer(i64) :: ik, ib, nk_irred, nb
+    integer :: dir
     real(r64), allocatable :: Knudsen(:, :), suppression_FS(:, :)
 
     !Number of IBZ wave vectors and bands
     nk_irred = size(indexlist_irred(:))
-    nb = size(vels_fbz(1,:,1))
+    nb = size(vels_fbz(1, :, 1))
 
     !Allocate boundary scattering rates and initialize to infinite crystal values
     allocate(thin_film_scatt_rates(nk_irred, nb))
     thin_film_scatt_rates = 0.0_r64
 
     if(normal == 'x') then
-       dir = 1_i64
+       dir = 1
     else if(normal == 'y') then
-       dir = 2_i64
+       dir = 2
     else if(normal == 'z') then
-       dir = 3_i64
+       dir = 3
     else
        call exit_with_message("Bad thin-film normal direction in calculate_thinfilm_scattrates. Exiting.")
     end if
 
     !Check finiteness of crystal
     if(finite_crys) then
        if(ballistic_limit) then !Large Knudsen number limit
-          do ik = 1, nk_irred
-             do ib = 1, nb
+          do ib = 1, nb
+             do ik = 1, nk_irred
                 thin_film_scatt_rates(ik, ib) = abs(vels_fbz(indexlist_irred(ik), ib, dir)) &
                      /height*1.e-6_r64 !THz
              end do
           end do
-          thin_film_scatt_rates = 2.0_r64*thin_film_scatt_rates
+          thin_film_scatt_rates = 2.0_r64*(1.0_r64 - specfac)/(1.0_r64 + specfac)*thin_film_scatt_rates
        else
           allocate(Knudsen(nk_irred, nb), suppression_FS(nk_irred, nb))
 

src/numerics.f90

@@ -482,6 +482,7 @@ subroutine read_input_and_setup(self, crys)
           if(self%phthinfilm) then
              write(*,"(A,(1E16.8,x),A,A,A)") 'Height for ph-thin-film scattering =', &
                   crys%thinfilm_height, 'mm along the ', crys%thinfilm_normal, ' direction'
+             write(*,"(A,1E16.8)") 'Specularity factor =', crys%specfac
           end if
           write(*, "(A, L)") "Include el-charged impurity interaction: ", self%elchimp
           write(*, "(A, L)") "Include el-boundary interaction: ", self%elbound

test/screening_comparison.f90

@@ -61,7 +61,7 @@ program screening_comparison
   beta = 1.0_r64/kB/crys%T
 
   !Create grid of probe |q|
-  allocate(qmags(el%nwv_irred))
+!!$  allocate(qmags(el%nwv_irred))
 !!$  do ik = 1, el%nwv_irred
 !!$     qmags(ik) = qdist(el%wavevecs_irred(ik, :), crys%reclattvecs)
 !!$  end do
@@ -92,7 +92,7 @@ program screening_comparison
 
   !Create bosonic energy mesh
   numomega = 100
-  call linspace(Omegas, 0.0_r64, 1.0_r64, numomega)
+  call linspace(Omegas, 0.0_r64, 10.0_r64*eF, numomega)
   call write2file_rank1_real("RPA_test_Omegas", Omegas)
 
   !Calculate analytic Im RPA dielectric function
@@ -199,10 +199,10 @@ subroutine calculate_Imeps(qmags, ens, chempot, m_eff, eF, kF, beta, Imeps)
 
     do iOmega = 1, size(ens)
        Imeps(:, iOmega) = &
-            real(log((1.0_r64 + &
-                exp(beta*(chempot - eF*(u(:, iOmega) - qmags(:)/kF/2.0_r64)**2)))/ &
-                (1.0_r64 + exp(beta*(chempot - eF*(u(:, iOmega) + &
-                           qmags(:)/kF/2.0_r64)**2)))))/(qmags(:)/kF)**3
+            real(log(&
+            (1.0_r64 + exp(beta*(chempot - eF*(u(:, iOmega) - qmags(:)/kF/2.0_r64)**2)))/ &
+            (1.0_r64 + exp(beta*(chempot - eF*(u(:, iOmega) + qmags(:)/kF/2.0_r64)**2)))))/ &
+            (qmags(:)/kF)**3
     end do
     Imeps(1, :) = 0.0_r64
     Imeps = (m_eff/me/bohr2nm/kF/eF/beta)*Imeps