Updated calculate_Y_OTF to use vector_allreps.

.github/workflows/build_test_fpm.yml

@@ -40,4 +40,5 @@ jobs:
         source fpm_config_caf.sh
         fpm build; fpm install && \
         fpm test test_autodiff test_misc test_periodictable && \
+        fpm test test_vector_allreps && \
         fpm test bte_regression --runner="sh test/3C-SiC/fpm_run_bte_caf.sh"

src/bte.f90

@@ -606,7 +606,7 @@ subroutine dragless_phbte_full(Tdir, self, num, crys, sym, ph, el)
     end if
 
     do it_ph = 1, num%maxiter
-       call iterate_bte_ph(crys%T, num, ph, el, self%ph_rta_rates_ibz, &
+       call iterate_bte_ph(crys%T, num, crys, ph, el, self%ph_rta_rates_ibz, &
             self%ph_field_term_T, self%ph_response_T)
 
        !Calculate phonon transport coefficients
@@ -731,9 +731,9 @@ subroutine dragfull_ephbtes(Tdir, self, num, crys, sym, ph, el)
        !Scheme: for each step of phonon response, fully iterate the electron response.
 
        !Iterate phonon response once          
-       call iterate_bte_ph(crys%T, num, ph, el, self%ph_rta_rates_ibz, &
+       call iterate_bte_ph(crys%T, num, crys, ph, el, self%ph_rta_rates_ibz, &
             self%ph_field_term_T, self%ph_response_T, self%el_response_T)
-       call iterate_bte_ph(crys%T, num, ph, el, self%ph_rta_rates_ibz, &
+       call iterate_bte_ph(crys%T, num, crys, ph, el, self%ph_rta_rates_ibz, &
             self%ph_field_term_E, self%ph_response_E, self%el_response_E)
 
        !Calculate phonon transport coefficients
@@ -1004,13 +1004,15 @@ subroutine calculate_field_term(species, field, nequiv, ibz2fbz_map, &
     end if
   end subroutine calculate_field_term
 
-  subroutine iterate_bte_ph(T, num, ph, el, rta_rates_ibz, &
+  subroutine iterate_bte_ph(T, num, crys, ph, el, rta_rates_ibz, &
        field_term, response_ph, response_el)
     !! Subroutine to iterate the phonon BTE one step.
     !! 
     !! T Temperature in K
-    !! drag Is drag included?
+    !! num Numerics object
+    !! crys Crystal object
     !! ph Phonon object
+    !! el Electron object
     !! rta_rates_ibz Phonon RTA scattering rates
     !! field_term Phonon field coupling term
     !! response_ph Phonon response function
@@ -1019,7 +1021,7 @@ subroutine iterate_bte_ph(T, num, ph, el, rta_rates_ibz, &
     type(phonon), intent(in) :: ph
     type(electron), intent(in) :: el
     type(numerics), intent(in) :: num
-    !logical, intent(in) :: drag
+    type(crystal), intent(in) :: crys
     real(r64), intent(in) :: T, rta_rates_ibz(:,:), field_term(:,:,:)
     real(r64), intent(in), optional :: response_el(:,:,:)
     real(r64), intent(inout) :: response_ph(:,:,:)
@@ -1136,7 +1138,7 @@ subroutine iterate_bte_ph(T, num, ph, el, rta_rates_ibz, &
 
           if(present(response_el)) then
              if(num%Y_OTF) then
-                call calculate_Y_OTF(el, ph, num, istate1, T, Y, istate_el1, istate_el2)
+                call calculate_Y_OTF(el, ph, num, crys, istate1, T, Y, istate_el1, istate_el2)
                 nprocs_phe = size(Y)
              else
                 !Set Y filename

src/interactions.f90

@@ -1817,7 +1817,7 @@ subroutine calculate_eph_interaction_ibzq(wann, crys, el, ph, num, key)
     sync all
   end subroutine calculate_eph_interaction_ibzq
 
-  subroutine calculate_Y_OTF(el, ph, num, istate, T, &
+  subroutine calculate_Y_OTF(el, ph, num, crys, istate, T, &
        Y_istate, istate_el1, istate_el2)
     !! On-the-fly (OTF), serial calcualator of the ph-e transition probability
     !! for all IBZ phonon wave vectors. This subroutine calculates
@@ -1827,6 +1827,7 @@ subroutine calculate_Y_OTF(el, ph, num, istate, T, &
     type(electron), intent(in) :: el
     type(phonon), intent(in) :: ph
     type(numerics), intent(in) :: num
+    type(crystal), intent(in) :: crys
     real(r64), intent(in) :: T
     real(r64), intent(out), allocatable :: Y_istate(:)
     integer(i64), intent(out), allocatable, optional :: istate_el1(:), istate_el2(:)
@@ -1841,6 +1842,7 @@ subroutine calculate_Y_OTF(el, ph, num, istate, T, &
     character(len = 1024) :: filename
     procedure(delta_fn), pointer :: delta_fn_ptr => null()
     logical :: keep_interaction_tally
+    type(vec) :: q_vec, k_vec, kp_vec
 
     !Do I need to keep a tally of the all the interacting states?
     keep_interaction_tally = present(istate_el1) .and. present(istate_el2)
@@ -1875,19 +1877,16 @@ subroutine calculate_Y_OTF(el, ph, num, istate, T, &
     !Energy of phonon
     en_ph = ph%ens(iq_fbz, s)
 
+    !Create phonon wave vector
+    q_vec = vec(iq_fbz, ph%wvmesh, crys%reclattvecs)
+
     !1/(1 + Bose factor) for phonon
     if(en_ph /= 0.0_r64) then
        invboseplus1 = 1.0_r64/(1.0_r64 + Bose(en_ph, T))
     else
        invboseplus1 = 0.0_r64
     end if
 
-    !Initial (IBZ blocks) wave vector (crystal coords.)
-    q = ph%wavevecs(iq_fbz, :)
-
-    !Convert from crystal to 0-based index vector
-    q_indvec = nint(q*ph%wvmesh)
-
     !Load g2_istate from disk for scattering rates calculation
     !Change to data output directory
     call chdir(trim(adjustl(num%g2dir)))
@@ -1911,19 +1910,13 @@ subroutine calculate_Y_OTF(el, ph, num, istate, T, &
     !Run over initial (in-window, FBZ blocks) electron wave vectors
     do ik = 1, el%nwv
        !Initial wave vector (crystal coords.)
-       k = el%wavevecs(ik, :)
-
-       !Convert from crystal to 0-based index vector
-       k_indvec = nint(k*el%wvmesh)
+       k_vec = vec(el%indexlist(ik), el%wvmesh, crys%reclattvecs)
 
        !Find final electron wave vector
-       kp_indvec = modulo(k_indvec + el%mesh_ref_array*q_indvec, el%wvmesh) !0-based index vector
-
-       !Muxed index of kp
-       ikp = mux_vector(kp_indvec, el%wvmesh, 0_i64)
+       kp_vec = vec_add(k_vec, q_vec, el%wvmesh, crys%reclattvecs)
 
        !Check if final electron wave vector is within energy window
-       call binsearch(el%indexlist, ikp, ikp_window)
+       call binsearch(el%indexlist, kp_vec%muxed_index, ikp_window)
        if(ikp_window < 0) cycle
 
        !Run over initial electron bands
@@ -2506,7 +2499,7 @@ subroutine calculate_ph_rta_rates(rta_rates_3ph, rta_rates_phe, num, crys, ph, e
 
           if(present(el)) then
              if(num%Y_OTF) then
-                call calculate_Y_OTF(el, ph, num, istate, crys%T, Y)
+                call calculate_Y_OTF(el, ph, num, crys, istate, crys%T, Y)
              else
                 !Set Y filename
                 filepath_Y = trim(adjustl(num%Ydir))//'/Y.istate'//trim(adjustl(tag))