berry_get_kubo_k Subroutine

  subroutine berry_get_kubo_k(pw90_berry, dis_manifold, fermi_energy_list, kpt_latt, &
                              pw90_band_deriv_degen, pw90_spin, ws_region, print_output, &
                              wannier_data, ws_distance, wigner_seitz, AA_R, HH_R, kubo_AH_k, &
                              kubo_H_k, SS_R, u_matrix, v_matrix, eigval, kpt, real_lattice, &
                              jdos_k, scissors_shift, mp_grid, num_bands, num_kpts, num_wann, &
                              num_valence_bands, effective_model, have_disentangled, spin_decomp, &
                              seedname, stdout, timer, error, comm, kubo_AH_k_spn, kubo_H_k_spn, &
                              jdos_k_spn)
    !================================================!
    !
    !! Contribution from point k to the complex interband optical
    !! conductivity, separated into Hermitian (H) and anti-Hermitian (AH)
    !! parts. Also returns the joint density of states
    !
    !================================================!

    use w90_constants, only: dp, cmplx_0, cmplx_i, pi
    use w90_comms, only: w90_comm_type
    use w90_utility, only: utility_diagonalize, utility_rotate, utility_w0gauss, &
                           utility_recip_lattice_base
    use w90_types, only: print_output_type, wannier_data_type, &
                         dis_manifold_type, ws_region_type, ws_distance_type, timer_list_type
    use w90_postw90_types, only: pw90_berry_mod_type, pw90_spin_mod_type, &
                                 pw90_band_deriv_degen_type, wigner_seitz_type
    use w90_postw90_common, only: pw90common_get_occ, pw90common_fourier_R_to_k_new, &
                                  pw90common_fourier_R_to_k_vec, pw90common_kmesh_spacing
    use w90_spin, only: spin_get_nk
    use w90_wan_ham, only: wham_get_D_h, wham_get_eig_deleig

    implicit none

    ! arguments
    type(pw90_berry_mod_type), intent(inout) :: pw90_berry
    type(dis_manifold_type), intent(in) :: dis_manifold
    real(kind=dp), allocatable, intent(in) :: fermi_energy_list(:)
    real(kind=dp), intent(in) :: kpt_latt(:, :)
    type(pw90_band_deriv_degen_type), intent(in) :: pw90_band_deriv_degen
    type(pw90_spin_mod_type), intent(in) :: pw90_spin
    type(print_output_type), intent(in) :: print_output
    type(ws_region_type), intent(in) :: ws_region
    type(w90_comm_type), intent(in) :: comm
    type(wannier_data_type), intent(in) :: wannier_data
    type(wigner_seitz_type), intent(inout) :: wigner_seitz
    type(ws_distance_type), intent(inout) :: ws_distance
    type(timer_list_type), intent(inout) :: timer
    type(w90_error_type), allocatable, intent(out) :: error

    integer, intent(in) :: num_wann, num_bands, num_kpts, num_valence_bands
    integer, intent(in) :: mp_grid(3)
    integer, intent(in) :: stdout

    real(kind=dp), intent(in) :: kpt(3)
    real(kind=dp), intent(out) :: jdos_k(:)
    real(kind=dp), intent(in) :: eigval(:, :)
    real(kind=dp), intent(in) :: real_lattice(3, 3)
    real(kind=dp), intent(in) :: scissors_shift

    complex(kind=dp), intent(out) :: kubo_H_k(:, :, :)
    complex(kind=dp), intent(out) :: kubo_AH_k(:, :, :)
    complex(kind=dp), intent(in) :: u_matrix(:, :, :), v_matrix(:, :, :)
    complex(kind=dp), allocatable, intent(inout) :: HH_R(:, :, :) !  <0n|r|Rm>
    complex(kind=dp), allocatable, intent(inout) :: AA_R(:, :, :, :) ! <0n|r|Rm>
    complex(kind=dp), allocatable, intent(inout) :: SS_R(:, :, :, :) ! <0n|sigma_x,y,z|Rm>

    character(len=50), intent(in) :: seedname
    logical, intent(in) :: have_disentangled
    logical, intent(in) :: spin_decomp
    logical, intent(in) :: effective_model

    ! Last three arguments should be present iff spin_decomp=T (but
    ! this is not checked: do it?)
    real(kind=dp), optional, intent(out) :: jdos_k_spn(:, :)
    complex(kind=dp), optional, intent(out) :: kubo_AH_k_spn(:, :, :, :)
    complex(kind=dp), optional, intent(out) :: kubo_H_k_spn(:, :, :, :)

    ! local variables
    complex(kind=dp), allocatable :: HH(:, :)
    complex(kind=dp), allocatable :: delHH(:, :, :)
    complex(kind=dp), allocatable :: UU(:, :)
    complex(kind=dp), allocatable :: D_h(:, :, :)
    complex(kind=dp), allocatable :: AA(:, :, :)

    real(kind=dp) :: recip_lattice(3, 3), volume
    ! Adaptive smearing
    !
    real(kind=dp) :: del_eig(num_wann, 3), joint_level_spacing, eta_smr, Delta_k, arg, vdum(3)
    real(kind=dp) :: eig(num_wann), occ(num_wann), delta, rfac1, rfac2, occ_prod, spn_nk(num_wann)
    complex(kind=dp) :: cfac, omega
    integer :: i, j, n, m, ifreq, ispn

    allocate (HH(num_wann, num_wann))
    allocate (delHH(num_wann, num_wann, 3))
    allocate (UU(num_wann, num_wann))
    allocate (D_h(num_wann, num_wann, 3))
    allocate (AA(num_wann, num_wann, 3))

    if (pw90_berry%kubo_smearing%use_adaptive) then
      call wham_get_eig_deleig(dis_manifold, kpt_latt, pw90_band_deriv_degen, ws_region, print_output, wannier_data, &
                               ws_distance, wigner_seitz, delHH, HH, HH_R, u_matrix, UU, v_matrix, &
                               del_eig, eig, eigval, kpt, real_lattice, scissors_shift, mp_grid, &
                               num_bands, num_kpts, num_wann, num_valence_bands, effective_model, &
                               have_disentangled, seedname, stdout, timer, error, comm)
      if (allocated(error)) return

      call utility_recip_lattice_base(real_lattice, recip_lattice, volume)
      Delta_k = pw90common_kmesh_spacing(pw90_berry%kmesh%mesh, recip_lattice)
    else
      call pw90common_fourier_R_to_k_new(ws_region, wannier_data, ws_distance, wigner_seitz, HH_R, &
                                         kpt, real_lattice, mp_grid, num_wann, error, comm, OO=HH, &
                                         OO_dx=delHH(:, :, 1), OO_dy=delHH(:, :, 2), &
                                         OO_dz=delHH(:, :, 3))
      if (allocated(error)) return

      call utility_diagonalize(HH, num_wann, eig, UU, error, comm)
      if (allocated(error)) return
    end if
    call pw90common_get_occ(fermi_energy_list(1), eig, occ, num_wann)

    call wham_get_D_h(delHH, D_h, UU, eig, num_wann)

    call pw90common_fourier_R_to_k_vec(ws_region, wannier_data, ws_distance, wigner_seitz, AA_R, &
                                       kpt, real_lattice, mp_grid, num_wann, error, comm, &
                                       OO_true=AA)
    if (allocated(error)) return

    do i = 1, 3
      AA(:, :, i) = utility_rotate(AA(:, :, i), UU, num_wann)
    end do
    AA = AA + cmplx_i*D_h ! Eq.(25) WYSV06

    ! Replace imaginary part of frequency with a fixed value
    if (.not. pw90_berry%kubo_smearing%use_adaptive .and. pw90_berry%kubo_smearing%fixed_width /= 0.0_dp) &
      pw90_berry%kubo_freq_list = real(pw90_berry%kubo_freq_list, dp) &
                                  + cmplx_i*pw90_berry%kubo_smearing%fixed_width

    kubo_H_k = cmplx_0
    kubo_AH_k = cmplx_0
    jdos_k = 0.0_dp
    if (spin_decomp) then
      call spin_get_nk(ws_region, pw90_spin, wannier_data, ws_distance, wigner_seitz, HH_R, SS_R, &
                       kpt, real_lattice, spn_nk, mp_grid, num_wann, error, comm)
      if (allocated(error)) return

      kubo_H_k_spn = cmplx_0
      kubo_AH_k_spn = cmplx_0
      jdos_k_spn = 0.0_dp
    end if
    do m = 1, num_wann
      do n = 1, num_wann
        if (n == m) cycle
        if (eig(m) > pw90_berry%kubo_eigval_max .or. eig(n) > pw90_berry%kubo_eigval_max) cycle
        if (spin_decomp) then
          if (spn_nk(n) >= 0 .and. spn_nk(m) >= 0) then
            ispn = 1 ! up --> up transition
          elseif (spn_nk(n) < 0 .and. spn_nk(m) < 0) then
            ispn = 2 ! down --> down
          else
            ispn = 3 ! spin-flip
          end if
        end if
        if (pw90_berry%kubo_smearing%use_adaptive) then
          ! Eq.(35) YWVS07
          vdum(:) = del_eig(m, :) - del_eig(n, :)
          joint_level_spacing = sqrt(dot_product(vdum(:), vdum(:)))*Delta_k
          eta_smr = min(joint_level_spacing*pw90_berry%kubo_smearing%adaptive_prefactor, &
                        pw90_berry%kubo_smearing%adaptive_max_width)
        else
          eta_smr = pw90_berry%kubo_smearing%fixed_width
        end if
        rfac1 = (occ(m) - occ(n))*(eig(m) - eig(n))
        occ_prod = occ(n)*(1.0_dp - occ(m))
        do ifreq = 1, pw90_berry%kubo_nfreq
          !
          ! Complex frequency for the anti-Hermitian conductivity
          !
          if (pw90_berry%kubo_smearing%use_adaptive) then
            omega = real(pw90_berry%kubo_freq_list(ifreq), dp) + cmplx_i*eta_smr
          else
            omega = pw90_berry%kubo_freq_list(ifreq)
          end if
          !
          ! Broadened delta function for the Hermitian conductivity and JDOS
          !
          arg = (eig(m) - eig(n) - real(omega, dp))/eta_smr
          ! If only Hermitean part were computed, could speed up
          ! by inserting here 'if(abs(arg)>10.0_dp) cycle'
          delta = utility_w0gauss(arg, pw90_berry%kubo_smearing%type_index, error, comm)/eta_smr
          if (allocated(error)) return
          !
          ! Lorentzian shape (for testing purposes)
          ! delta=1.0_dp/(1.0_dp+arg*arg)/pi
          ! delta=delta/eta_smr
          !
          jdos_k(ifreq) = jdos_k(ifreq) + occ_prod*delta
          if (spin_decomp) &
            jdos_k_spn(ispn, ifreq) = jdos_k_spn(ispn, ifreq) + occ_prod*delta
          cfac = cmplx_i*rfac1/(eig(m) - eig(n) - omega)
          rfac2 = -pi*rfac1*delta
          do j = 1, 3
            do i = 1, 3
              kubo_H_k(i, j, ifreq) = kubo_H_k(i, j, ifreq) &
                                      + rfac2*AA(n, m, i)*AA(m, n, j)
              kubo_AH_k(i, j, ifreq) = kubo_AH_k(i, j, ifreq) &
                                       + cfac*AA(n, m, i)*AA(m, n, j)
              if (spin_decomp) then
                kubo_H_k_spn(i, j, ispn, ifreq) = &
                  kubo_H_k_spn(i, j, ispn, ifreq) &
                  + rfac2*AA(n, m, i)*AA(m, n, j)
                kubo_AH_k_spn(i, j, ispn, ifreq) = &
                  kubo_AH_k_spn(i, j, ispn, ifreq) &
                  + cfac*AA(n, m, i)*AA(m, n, j)
              end if
            end do
          end do
        end do
      end do
    end do

  end subroutine berry_get_kubo_k