Debugged calculate_eph_interaction_ibzk where a wavevector access bug…

… was introduced in the last commit.

app/elphbolt.f90

@@ -154,13 +154,13 @@ program elphbolt
 
      call subtitle("Calculating interactions...")
 
-     !TEST/DUBUG
-     !Calculate RPA dielectric for q over Gamma-Gamma along x over a uniform boson energy mesh
-     call t_event%start_timer('RPA dielectric')
-     call calculate_RPA_dielectric_3d_G0_qpath(el, crys, num)
-     call t_event%end_timer('RPA dielectric')
-     call exit
-     !!
+!!$     !TEST/DUBUG
+!!$     !Calculate RPA dielectric for q over Gamma-Gamma along x over a uniform boson energy mesh
+!!$     call t_event%start_timer('RPA dielectric')
+!!$     call calculate_RPA_dielectric_3d_G0_qpath(el, crys, num)
+!!$     call t_event%end_timer('RPA dielectric')
+!!$     call exit
+!!$     !!
 
      if(num%phdef_Tmat) then
         !Calculate phonon-defect interactions

src/interactions.f90

@@ -2004,7 +2004,7 @@ subroutine calculate_eph_interaction_ibzk(wann, crys, el, ph, num, key)
 
     !Local variables
     integer(i64) :: nstates_irred, istate, m, ik, ik_fbz, n, ikp, s, &
-         iq, start, end, chunk, count, nprocs, num_active_images
+         iq_fine, iq_coarse, start, end, chunk, count, nprocs, num_active_images
     real(r64) :: ph_ens_iq(1, ph%numbands), qlist(1, 3), &
          const, bosefac, fermi_minus_fac, fermi_plus_fac, en_ph, en_el, delta, occup_fac
     real(r64), allocatable :: g2_istate(:), Xplus_istate(:), Xminus_istate(:)
@@ -2071,7 +2071,7 @@ subroutine calculate_eph_interaction_ibzk(wann, crys, el, ph, num, key)
           end if
 
           !Get the muxed index of FBZ wave vector from the IBZ index list
-          ik_fbz = el%indexlist(ik)
+          ik_fbz = el%indexlist_irred(ik)
 
           !Electron energy
           en_el = el%ens_irred(ik, m)
@@ -2116,15 +2116,16 @@ subroutine calculate_eph_interaction_ibzk(wann, crys, el, ph, num, key)
           !Run over final (FBZ blocks) electron wave vectors
           do ikp = 1, el%nwv
              !Create final electron wave vector
-             kp_vec = vec(ikp, el%wvmesh, crys%reclattvecs)
+             kp_vec = vec(el%indexlist(ikp), el%wvmesh, crys%reclattvecs)
 
              !Find interacting phonon wave vector.
              !This is represented with respect to the electronic mesh.
              q_vec = vec_sub(kp_vec, k_vec, el%wvmesh, crys%reclattvecs)
-
+             iq_fine = q_vec%muxed_index
+
              !Note that q, k, and k' are all on the same mesh.
              !However, there are some common q vector on which the
-             !phonon quantities have already been computer. Here, compute
+             !phonon quantities have already been computed. Here, compute
              !only the new k' - k quantities.
              needfinephon = .false.
              if(any(mod(q_vec%int(:), el%mesh_ref_array) /= 0_i64)) then
@@ -2133,8 +2134,9 @@ subroutine calculate_eph_interaction_ibzk(wann, crys, el, ph, num, key)
                 !Calculate the fine mesh phonon.
                 qlist(1, :) = q_vec%frac
                 call wann%ph_wann(crys, 1_i64, qlist, ph_ens_iq, ph_evecs_iq)
-             else !Get the q vector represented in the phonon mesh
+             else !Get the q vector represented in the (coarser) phonon mesh
                 q_vec_coarse = vec_change_grid(q_vec, ph%wvmesh)
+                iq_coarse = q_vec_coarse%muxed_index
              end if
 
              !Run over final electron bands
@@ -2157,7 +2159,7 @@ subroutine calculate_eph_interaction_ibzk(wann, crys, el, ph, num, key)
                       else
                          g2_istate(count) = wann%g2(crys, k_vec%frac, q_vec_coarse%frac, &
                               el%evecs_irred(ik, m, :), el%evecs(ikp, n, :), &
-                              ph%evecs(iq, s, :), ph%ens(iq, s), &
+                              ph%evecs(iq_coarse, s, :), ph%ens(iq_coarse, s), &
                               gkRp_ik, 'ph')
                       end if
                    end if
@@ -2167,7 +2169,7 @@ subroutine calculate_eph_interaction_ibzk(wann, crys, el, ph, num, key)
                       if(needfinephon) then
                          en_ph = ph_ens_iq(1, s)
                       else
-                         en_ph = ph%ens(iq, s)
+                         en_ph = ph%ens(iq_coarse, s)
                       end if
 
                       !Bose and Fermi factors
@@ -2212,9 +2214,9 @@ subroutine calculate_eph_interaction_ibzk(wann, crys, el, ph, num, key)
                       if(needfinephon) then
                          !Write fine phonon index as negative so that the iterator
                          !knows to interpolate phonon quantities at this wave vector.
-                         istate_ph(count) = -mux_state(ph%numbands, s, iq)
+                         istate_ph(count) = -mux_state(ph%numbands, s, iq_fine)
                       else
-                         istate_ph(count) = mux_state(ph%numbands, s, iq)
+                         istate_ph(count) = mux_state(ph%numbands, s, iq_coarse)
                       end if
                    end if
                 end do !s

src/vector.f90

@@ -103,7 +103,7 @@ end function vector_allreps_sub
 
   pure function vector_allreps_change_grid(vin, grid) result(vout)
     !! Change grid.
-    !! Beware, this is only well-defined is the vector is representable
+    !! Beware, this is only well-defined if the vector is representable
     !! in both original and the new grids.
 
     type(vector_allreps), intent(in) :: vin

test/screening_comparison.f90

@@ -18,7 +18,7 @@ program screening_comparison
 
   !concentration and temperature
   real(r64), parameter :: conc = 1.0e18 !cm^-3
-  real(r64), parameter :: T = 0.9
+  real(r64), parameter :: T = 1.0
 
   !real(r64), parameter :: m_eff = 0.267*me !Si
   !real(r64), parameter :: m_eff = 0.2*me !wGaN
@@ -195,46 +195,46 @@ subroutine calculate_Imeps(qmags, ens, chempot, m_eff, eF, kF, beta, Imeps)
     real(r64), intent(in) :: qmags(:), ens(:), chempot, m_eff, eF, kF, beta
     real(r64), allocatable :: Imeps(:, :)
 
-!!$    !Locals
-!!$    integer :: iOmega
-!!$    real(r64) :: u(size(qmags), size(ens))
-!!$    
-!!$    allocate(Imeps(size(qmags), size(ens)))
-!!$    
-!!$    call outer(0.5_r64*kF/qmags, ens/eF, u)
-!!$
-!!$    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
-!!$    end do
-!!$    Imeps(1, :) = 0.0_r64
-!!$    Imeps = (m_eff/me/bohr2nm/kF/eF/beta)*Imeps
-
     !Locals
-    integer :: iOmega, iq
-    real(r64) :: E1(size(qmags), size(ens)), E2(size(qmags), size(ens)), Eq(size(qmags))
-
+    integer :: iOmega
+    real(r64) :: u(size(qmags), size(ens))
+    
     allocate(Imeps(size(qmags), size(ens)))
-
-    Eq = energy_parabolic(qmags, m_eff)
-
-    do iq = 1, size(qmags)
-       E1(iq, :) = (ens - Eq)**2/4.0/Eq
-       E2(iq, :) = (ens + Eq)**2/4.0/Eq
-    end do
+
+    call outer(0.5_r64*kF/qmags, ens/eF, u)
 
     do iOmega = 1, size(ens)
        Imeps(:, iOmega) = &
             real(log(&
-            (1.0_r64 + exp(beta*(chempot - E1(:, iOmega))))/ &
-            (1.0_r64 + exp(beta*(chempot - E2(:, iOmega))))))/ &
-            (qmags(:))**3
+            (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 = 2.0*me/bohr2nm/(beta*hbar*hbar_eVps)*1.0e6*Imeps
+    Imeps(1, :) = 0.0_r64
+    Imeps = (m_eff/me/bohr2nm/kF/eF/beta)*Imeps
+
+!!$    !Locals
+!!$    integer :: iOmega, iq
+!!$    real(r64) :: E1(size(qmags), size(ens)), E2(size(qmags), size(ens)), Eq(size(qmags))
+!!$
+!!$    allocate(Imeps(size(qmags), size(ens)))
+!!$
+!!$    Eq = energy_parabolic(qmags, m_eff)
+!!$
+!!$    do iq = 1, size(qmags)
+!!$       E1(iq, :) = (ens - Eq)**2/4.0/Eq
+!!$       E2(iq, :) = (ens + Eq)**2/4.0/Eq
+!!$    end do    
+!!$
+!!$    do iOmega = 1, size(ens)
+!!$       Imeps(:, iOmega) = &
+!!$            log(&
+!!$            (1.0_r64 + exp(beta*(chempot - E1(:, iOmega))))/ &
+!!$            (1.0_r64 + exp(beta*(chempot - E2(:, iOmega)))))/ &
+!!$            (qmags(:))**3
+!!$    end do
+!!$    !Imeps(1, :) = 0.0_r64
+!!$    Imeps = 2.0*m_eff**2/(me*bohr2nm*beta*hbar**2)*1.0e6*qe*Imeps
   end subroutine calculate_Imeps
 
   subroutine calculate_Reeps(qmags, ens, chempot, m_eff, eF, eplasmon, &
@@ -267,8 +267,8 @@ subroutine calculate_Reeps(qmags, ens, chempot, m_eff, eF, eplasmon, &
     !Here we need an extra factor of ms/me.
 
     !Magic numbers?
-    ngrid = 500
-    ymax = 10.0_r64
+    ngrid = 10000
+    ymax = 20.0_r64
 
     allocate(y(ngrid), I0(ngrid), Reeps(size(qmags), size(ens)))