dos_get_k Subroutine

public subroutine dos_get_k(num_elec_per_state, ws_region, kpt, EnergyArray, eig_k, dos_k, num_wann, wannier_data, real_lattice, mp_grid, pw90_dos, spin_decomp, pw90_spin, ws_distance, wigner_seitz, HH_R, SS_R, smearing, error, comm, levelspacing_k, UU)

Uses

  • proc~~dos_get_k~~UsesGraph proc~dos_get_k dos_get_k module~w90_comms w90_comms proc~dos_get_k->module~w90_comms module~w90_constants w90_constants proc~dos_get_k->module~w90_constants module~w90_postw90_types w90_postw90_types proc~dos_get_k->module~w90_postw90_types module~w90_spin w90_spin proc~dos_get_k->module~w90_spin module~w90_types w90_types proc~dos_get_k->module~w90_types module~w90_utility w90_utility proc~dos_get_k->module~w90_utility module~w90_comms->module~w90_constants module~w90_error_base w90_error_base module~w90_comms->module~w90_error_base module~w90_postw90_types->module~w90_comms module~w90_postw90_types->module~w90_constants module~w90_spin->module~w90_constants module~w90_error w90_error module~w90_spin->module~w90_error module~w90_types->module~w90_constants module~w90_utility->module~w90_comms module~w90_utility->module~w90_constants module~w90_error->module~w90_comms module~w90_error->module~w90_error_base

This subroutine calculates the contribution to the DOS of a single k point

\todo still to do: adapt spin_get_nk to read in input the UU rotation matrix

\note This routine simply provides the dos contribution of a given point. This must be externally summed after proper weighting. The weight factor (for a full BZ sampling with N^3 points) is 1/N^3 if we want the final DOS to be normalized to the total number of electrons. \note The only factor that is included INSIDE this routine is the spin degeneracy factor (=num_elec_per_state variable) \note The EnergyArray is assumed to be evenly spaced (and the energy spacing is taken from EnergyArray(2)-EnergyArray(1)) \note The routine is assuming that EnergyArray has at least two elements. \note The dos_k array must have dimensions size(EnergyArray) * ndim, where ndim=1 if spin_decomp==false, or ndim=3 if spin_decomp==true. This is not checked. \note If smearing/binwidth < min_smearing_binwidth_ratio, no smearing is applied (for that k point)

\param kpt the three coordinates of the k point vector whose DOS contribution we want to calculate (in relative coordinates) \param EnergyArray array with the energy grid on which to calculate the DOS (in eV) It must have at least two elements \param eig_k array with the eigenvalues at the given k point (in eV) \param dos_k array in which the contribution is stored. Three dimensions: dos_k(energyidx, spinidx), where: - energyidx is the index of the energies, corresponding to the one of the EnergyArray array; - spinidx=1 contains the total dos; if if spin_decomp==.true., then spinidx=2 and spinidx=3 contain the spin-up and spin-down contributions to the DOS \param smr_index index that tells the kind of smearing \param smr_fixed_en_width optional parameter with the fixed energy for smearing, in eV. Can be provided only if the levelspacing_k parameter is NOT given \param adpt_smr_fac optional parameter with the factor for the adaptive smearing. Can be provided only if the levelspacing_k parameter IS given \param levelspacing_k optional array with the level spacings, i.e. how much each level changes near a given point of the interpolation mesh, as given by the dos_get_levelspacing() routine If present: adaptive smearing If not present: fixed-energy-width smearing

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: num_elec_per_state
type(ws_region_type), intent(in) :: ws_region
real(kind=dp), intent(in) :: kpt(3)
real(kind=dp), intent(in) :: EnergyArray(:)
real(kind=dp), intent(in) :: eig_k(:)
real(kind=dp), intent(out) :: dos_k(:,:)
integer, intent(in) :: num_wann
type(wannier_data_type), intent(in) :: wannier_data
real(kind=dp), intent(in) :: real_lattice(3,3)
integer, intent(in) :: mp_grid(3)
type(pw90_dos_mod_type), intent(in) :: pw90_dos
logical, intent(in) :: spin_decomp
type(pw90_spin_mod_type), intent(in) :: pw90_spin
type(ws_distance_type), intent(inout) :: ws_distance
type(wigner_seitz_type), intent(in) :: wigner_seitz
complex(kind=dp), intent(inout), allocatable :: HH_R(:,:,:)
complex(kind=dp), intent(inout), allocatable :: SS_R(:,:,:,:)
type(pw90_smearing_type), intent(in) :: smearing
type(w90_error_type), intent(out), allocatable :: error
type(w90_comm_type), intent(in) :: comm
real(kind=dp), intent(in), optional :: levelspacing_k(:)
complex(kind=dp), intent(in), optional :: UU(:,:)

Calls

proc~~dos_get_k~~CallsGraph proc~dos_get_k dos_get_k proc~set_error_input set_error_input proc~dos_get_k->proc~set_error_input proc~spin_get_nk spin_get_nk proc~dos_get_k->proc~spin_get_nk proc~utility_w0gauss utility_w0gauss proc~dos_get_k->proc~utility_w0gauss proc~comms_sync_error comms_sync_error proc~set_error_input->proc~comms_sync_error proc~set_base_error set_base_error proc~set_error_input->proc~set_base_error proc~pw90common_fourier_r_to_k pw90common_fourier_R_to_k proc~spin_get_nk->proc~pw90common_fourier_r_to_k proc~utility_diagonalize utility_diagonalize proc~spin_get_nk->proc~utility_diagonalize proc~utility_rotate_diag utility_rotate_diag proc~spin_get_nk->proc~utility_rotate_diag proc~utility_w0gauss->proc~set_error_input proc~set_error_fatal set_error_fatal proc~utility_diagonalize->proc~set_error_fatal zhpevx zhpevx proc~utility_diagonalize->zhpevx proc~utility_matmul_diag utility_matmul_diag proc~utility_rotate_diag->proc~utility_matmul_diag proc~utility_zgemm_new utility_zgemm_new proc~utility_rotate_diag->proc~utility_zgemm_new proc~set_error_fatal->proc~comms_sync_error proc~set_error_fatal->proc~set_base_error zgemm zgemm proc~utility_zgemm_new->zgemm

Called by

proc~~dos_get_k~~CalledByGraph proc~dos_get_k dos_get_k proc~calctdfanddos calcTDFandDOS proc~calctdfanddos->proc~dos_get_k proc~dos_main dos_main proc~dos_main->proc~dos_get_k proc~boltzwann_main boltzwann_main proc~boltzwann_main->proc~calctdfanddos program~postw90 postw90 program~postw90->proc~dos_main program~postw90->proc~boltzwann_main

Source Code

  subroutine dos_get_k(num_elec_per_state, ws_region, kpt, EnergyArray, eig_k, dos_k, num_wann, &
                       wannier_data, real_lattice, mp_grid, pw90_dos, spin_decomp, &
                       pw90_spin, ws_distance, wigner_seitz, HH_R, SS_R, &
                       smearing, error, comm, levelspacing_k, UU)
    !================================================!
    use w90_constants, only: dp, smearing_cutoff, min_smearing_binwidth_ratio
    use w90_utility, only: utility_w0gauss
    use w90_postw90_types, only: pw90_spin_mod_type, pw90_dos_mod_type, pw90_smearing_type, &
                                 wigner_seitz_type
    use w90_types, only: wannier_data_type, ws_region_type, ws_distance_type
    use w90_spin, only: spin_get_nk
    use w90_utility, only: utility_w0gauss
    use w90_comms, only: w90_comm_type

    ! Arguments
    type(pw90_dos_mod_type), intent(in) :: pw90_dos
    type(pw90_spin_mod_type), intent(in) :: pw90_spin
    type(ws_region_type), intent(in) :: ws_region
    type(wannier_data_type), intent(in) :: wannier_data
    type(wigner_seitz_type), intent(in) :: wigner_seitz
    type(ws_distance_type), intent(inout) :: ws_distance
    type(pw90_smearing_type), intent(in) :: smearing
    type(w90_error_type), allocatable, intent(out) :: error
    type(w90_comm_type), intent(in) :: comm

    integer, intent(in) :: mp_grid(3)
    integer, intent(in) :: num_elec_per_state
    integer, intent(in) :: num_wann
    !integer, intent(in) :: smr_index

    real(kind=dp), intent(in) :: kpt(3)
    real(kind=dp), intent(in) :: eig_k(:)
    real(kind=dp), intent(in) :: EnergyArray(:)
    real(kind=dp), intent(in) :: real_lattice(3, 3)
    real(kind=dp), intent(out) :: dos_k(:, :)
    real(kind=dp), intent(in), optional :: levelspacing_k(:)
    !real(kind=dp), intent(in), optional :: adpt_smr_fac
    !real(kind=dp), intent(in), optional :: adpt_smr_max
    !real(kind=dp), intent(in), optional :: smr_fixed_en_width

    complex(kind=dp), allocatable, intent(inout) :: HH_R(:, :, :)
    complex(kind=dp), allocatable, intent(inout) :: SS_R(:, :, :, :)
    complex(kind=dp), intent(in), optional :: UU(:, :)

    logical, intent(in) :: spin_decomp

    ! local variables
    real(kind=dp) :: eta_smr, arg ! Adaptive smearing
    real(kind=dp) :: rdum, spn_nk(num_wann), alpha_sq, beta_sq
    real(kind=dp) :: binwidth, r_num_elec_per_state
    integer :: i, j, loop_f, min_f, max_f
    logical :: DoSmearing

    if (present(levelspacing_k)) then
      if (.not. smearing%use_adaptive) then
        call set_error_input(error, 'Cannot call doskpt with levelspacing_k and ' &
                             //'without adptative smearing', comm)
        return
      end if
    else
      if (smearing%use_adaptive) then
        call set_error_input(error, 'Cannot call doskpt without levelspacing_k and ' &
                             //'with adptative smearing', comm)
        return
      end if
    end if

    r_num_elec_per_state = real(num_elec_per_state, kind=dp)

    ! Get spin projections for every band
    !
    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

    end if

    binwidth = EnergyArray(2) - EnergyArray(1)

    dos_k = 0.0_dp
    do i = 1, num_wann
      if (spin_decomp) then
        ! Contribution to spin-up DOS of Bloch spinor with component
        ! (alpha,beta) with respect to the chosen quantization axis
        alpha_sq = (1.0_dp + spn_nk(i))/2.0_dp ! |alpha|^2
        ! Contribution to spin-down DOS
        beta_sq = 1.0_dp - alpha_sq ! |beta|^2 = 1 - |alpha|^2
      end if

      if (.not. present(levelspacing_k)) then
        eta_smr = smearing%fixed_width
      else
        ! Eq.(35) YWVS07
        eta_smr = min(levelspacing_k(i)*smearing%adaptive_prefactor, smearing%adaptive_max_width)
!          eta_smr=max(eta_smr,min_smearing_binwidth_ratio) !! No: it would render the next if always false
      end if

      ! Faster optimization: I precalculate the indices
      if (eta_smr/binwidth < min_smearing_binwidth_ratio) then
        min_f = max(nint((eig_k(i) - EnergyArray(1))/ &
                         (EnergyArray(size(EnergyArray)) - EnergyArray(1)) &
                         *real(size(EnergyArray) - 1, kind=dp)) + 1, 1)
        max_f = min(nint((eig_k(i) - EnergyArray(1))/ &
                         (EnergyArray(size(EnergyArray)) - EnergyArray(1)) &
                         *real(size(EnergyArray) - 1, kind=dp)) + 1, size(EnergyArray))
        DoSmearing = .false.
      else
        min_f = max(nint((eig_k(i) - smearing_cutoff*eta_smr - EnergyArray(1))/ &
                         (EnergyArray(size(EnergyArray)) - EnergyArray(1)) &
                         *real(size(EnergyArray) - 1, kind=dp)) + 1, 1)
        max_f = min(nint((eig_k(i) + smearing_cutoff*eta_smr - EnergyArray(1))/ &
                         (EnergyArray(size(EnergyArray)) - EnergyArray(1)) &
                         *real(size(EnergyArray) - 1, kind=dp)) + 1, size(EnergyArray))
        DoSmearing = .true.
      end if

      do loop_f = min_f, max_f
        ! kind of smearing read from input (internal smearing_index variable)
        if (DoSmearing) then
          arg = (EnergyArray(loop_f) - eig_k(i))/eta_smr
          rdum = utility_w0gauss(arg, smearing%type_index, error, comm)/eta_smr
          if (allocated(error)) return
        else
          rdum = 1._dp/(EnergyArray(2) - EnergyArray(1))
        end if

        !
        ! Contribution to total DOS
        !
        if (pw90_dos%num_project == num_wann) then
          !
          ! Total DOS (default): do not loop over j, to save time
          !
          dos_k(loop_f, 1) = dos_k(loop_f, 1) + rdum*r_num_elec_per_state
          ! [GP] I don't put num_elec_per_state here below: if we are
          ! calculating the spin decomposition, we should be doing a
          ! calcultation with spin-orbit, and thus num_elec_per_state=1!
          if (spin_decomp) then
            ! Spin-up contribution
            dos_k(loop_f, 2) = dos_k(loop_f, 2) + rdum*alpha_sq
            ! Spin-down contribution
            dos_k(loop_f, 3) = dos_k(loop_f, 3) + rdum*beta_sq
          end if
        else ! 0<num_dos_project<num_wann
          !
          ! Partial DOS, projected onto the WFs with indices
          ! n=dos_project(1:num_dos_project)
          !
          do j = 1, pw90_dos%num_project
            dos_k(loop_f, 1) = dos_k(loop_f, 1) + rdum*r_num_elec_per_state &
                               *abs(UU(pw90_dos%project(j), i))**2
            if (spin_decomp) then
              ! Spin-up contribution
              dos_k(loop_f, 2) = dos_k(loop_f, 2) &
                                 + rdum*alpha_sq*abs(UU(pw90_dos%project(j), i))**2
              ! Spin-down contribution
              dos_k(loop_f, 3) = dos_k(loop_f, 3) &
                                 + rdum*beta_sq*abs(UU(pw90_dos%project(j), i))**2
            end if
          end do
        end if
      end do !loop_f
    end do !loop over bands

  end subroutine dos_get_k