berry_get_kubo_k Subroutine

  1. berry.F90
  2. w90_berry
  3. berry_get_kubo_k

private subroutine berry_get_kubo_k(kpt, kubo_H_k, kubo_AH_k, jdos_k, kubo_H_k_spn, kubo_AH_k_spn, jdos_k_spn)

Uses

  • proc~~berry_get_kubo_k~~UsesGraph proc~berry_get_kubo_k berry_get_kubo_k module~w90_postw90_common w90_postw90_common proc~berry_get_kubo_k->module~w90_postw90_common module~w90_utility w90_utility proc~berry_get_kubo_k->module~w90_utility module~w90_constants w90_constants proc~berry_get_kubo_k->module~w90_constants module~w90_get_oper w90_get_oper proc~berry_get_kubo_k->module~w90_get_oper module~w90_spin w90_spin proc~berry_get_kubo_k->module~w90_spin module~w90_wan_ham w90_wan_ham proc~berry_get_kubo_k->module~w90_wan_ham module~w90_parameters w90_parameters proc~berry_get_kubo_k->module~w90_parameters module~w90_postw90_common->module~w90_constants module~w90_comms w90_comms module~w90_postw90_common->module~w90_comms module~w90_utility->module~w90_constants module~w90_get_oper->module~w90_constants module~w90_spin->module~w90_constants module~w90_wan_ham->module~w90_constants module~w90_parameters->module~w90_constants module~w90_io w90_io module~w90_parameters->module~w90_io module~w90_comms->module~w90_constants module~w90_comms->module~w90_io module~w90_io->module~w90_constants

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

Arguments

Type IntentOptional AttributesName
real(kind=dp), intent(in) :: kpt(3)
complex(kind=dp), intent(out), dimension(:, :, :):: kubo_H_k
complex(kind=dp), intent(out), dimension(:, :, :):: kubo_AH_k
real(kind=dp), intent(out), dimension(:):: jdos_k
complex(kind=dp), intent(out), optional dimension(:, :, :, :):: kubo_H_k_spn
complex(kind=dp), intent(out), optional dimension(:, :, :, :):: kubo_AH_k_spn
real(kind=dp), intent(out), optional dimension(:, :):: jdos_k_spn

Calls

proc~~berry_get_kubo_k~~CallsGraph proc~berry_get_kubo_k berry_get_kubo_k interface~pw90common_kmesh_spacing pw90common_kmesh_spacing proc~berry_get_kubo_k->interface~pw90common_kmesh_spacing proc~pw90common_fourier_r_to_k_vec pw90common_fourier_R_to_k_vec proc~berry_get_kubo_k->proc~pw90common_fourier_r_to_k_vec proc~wham_get_eig_deleig wham_get_eig_deleig proc~berry_get_kubo_k->proc~wham_get_eig_deleig proc~wham_get_d_h wham_get_D_h proc~berry_get_kubo_k->proc~wham_get_d_h proc~utility_w0gauss utility_w0gauss proc~berry_get_kubo_k->proc~utility_w0gauss proc~utility_rotate utility_rotate proc~berry_get_kubo_k->proc~utility_rotate proc~utility_diagonalize utility_diagonalize proc~berry_get_kubo_k->proc~utility_diagonalize proc~spin_get_nk spin_get_nk proc~berry_get_kubo_k->proc~spin_get_nk proc~kmesh_spacing_mesh kmesh_spacing_mesh interface~pw90common_kmesh_spacing->proc~kmesh_spacing_mesh proc~kmesh_spacing_singleinteger kmesh_spacing_singleinteger interface~pw90common_kmesh_spacing->proc~kmesh_spacing_singleinteger proc~wham_get_eig_deleig->proc~utility_diagonalize proc~get_hh_r get_HH_R proc~wham_get_eig_deleig->proc~get_hh_r proc~pw90common_fourier_r_to_k pw90common_fourier_R_to_k proc~wham_get_eig_deleig->proc~pw90common_fourier_r_to_k proc~wham_get_d_h->proc~utility_rotate proc~io_error io_error proc~utility_diagonalize->proc~io_error proc~spin_get_nk->proc~utility_diagonalize proc~utility_rotate_diag utility_rotate_diag proc~spin_get_nk->proc~utility_rotate_diag proc~spin_get_nk->proc~pw90common_fourier_r_to_k proc~get_hh_r->proc~io_error proc~get_win_min get_win_min proc~get_hh_r->proc~get_win_min proc~fourier_q_to_r fourier_q_to_R proc~get_hh_r->proc~fourier_q_to_r proc~io_file_unit io_file_unit proc~get_hh_r->proc~io_file_unit proc~utility_zgemm_new utility_zgemm_new proc~utility_rotate_diag->proc~utility_zgemm_new proc~utility_matmul_diag utility_matmul_diag proc~utility_rotate_diag->proc~utility_matmul_diag

Called by

proc~~berry_get_kubo_k~~CalledByGraph proc~berry_get_kubo_k berry_get_kubo_k proc~berry_main berry_main proc~berry_main->proc~berry_get_kubo_k

Contents

Source Code

berry_get_kubo_k

Source Code

  subroutine berry_get_kubo_k(kpt, kubo_H_k, kubo_AH_k, jdos_k, &
                              kubo_H_k_spn, kubo_AH_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_utility, only: utility_diagonalize, utility_rotate, utility_w0gauss
    use w90_parameters, only: num_wann, kubo_nfreq, kubo_freq_list, &
      fermi_energy_list, kubo_eigval_max, &
      kubo_adpt_smr, kubo_smr_fixed_en_width, &
      kubo_adpt_smr_max, kubo_adpt_smr_fac, &
      kubo_smr_index, berry_kmesh, spin_decomp
    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_wan_ham, only: wham_get_D_h, wham_get_eig_deleig
    use w90_get_oper, only: HH_R, AA_R
    use w90_spin, only: spin_get_nk

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

    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(:, :, :)

    ! Adaptive smearing
    !
    real(kind=dp)    :: del_eig(num_wann, 3), joint_level_spacing, &
                        eta_smr, Delta_k, arg, vdum(3)

    integer          :: i, j, n, m, ifreq, ispn
    real(kind=dp)    :: eig(num_wann), occ(num_wann), delta, &
                        rfac1, rfac2, occ_prod, spn_nk(num_wann)
    complex(kind=dp) :: cfac, omega

    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 (kubo_adpt_smr) then
      call wham_get_eig_deleig(kpt, eig, del_eig, HH, delHH, UU)
      Delta_k = pw90common_kmesh_spacing(berry_kmesh)
    else
      call pw90common_fourier_R_to_k_new(kpt, HH_R, OO=HH, &
                                         OO_dx=delHH(:, :, 1), &
                                         OO_dy=delHH(:, :, 2), &
                                         OO_dz=delHH(:, :, 3))
      call utility_diagonalize(HH, num_wann, eig, UU)
    endif
    call pw90common_get_occ(eig, occ, fermi_energy_list(1))
    call wham_get_D_h(delHH, UU, eig, D_h)

    call pw90common_fourier_R_to_k_vec(kpt, AA_R, OO_true=AA)
    do i = 1, 3
      AA(:, :, i) = utility_rotate(AA(:, :, i), UU, num_wann)
    enddo
    AA = AA + cmplx_i*D_h ! Eq.(25) WYSV06

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

    kubo_H_k = cmplx_0
    kubo_AH_k = cmplx_0
    jdos_k = 0.0_dp
    if (spin_decomp) then
      call spin_get_nk(kpt, spn_nk)
      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) > kubo_eigval_max .or. eig(n) > 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 (kubo_adpt_smr) 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*kubo_adpt_smr_fac, &
                        kubo_adpt_smr_max)
        else
          eta_smr = kubo_smr_fixed_en_width
        endif
        rfac1 = (occ(m) - occ(n))*(eig(m) - eig(n))
        occ_prod = occ(n)*(1.0_dp - occ(m))
        do ifreq = 1, kubo_nfreq
          !
          ! Complex frequency for the anti-Hermitian conductivity
          !
          if (kubo_adpt_smr) then
            omega = real(kubo_freq_list(ifreq), dp) + cmplx_i*eta_smr
          else
            omega = kubo_freq_list(ifreq)
          endif
          !
          ! 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, kubo_smr_index)/eta_smr
          !
          ! 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)
              endif
            enddo
          enddo
        enddo
      enddo
    enddo

  end subroutine berry_get_kubo_k