TDF_kpt Subroutine

private subroutine TDF_kpt(pw90_boltzwann, ws_region, pw90_spin, wannier_data, ws_distance, wigner_seitz, HH_R, SS_R, deleig_k, eig_k, EnergyArray, kpt, real_lattice, TDF_k, mp_grid, num_wann, num_elec_per_state, spin_decomp, error, comm)

Uses

  • proc~~tdf_kpt~~UsesGraph proc~tdf_kpt TDF_kpt module~w90_comms w90_comms proc~tdf_kpt->module~w90_comms module~w90_constants w90_constants proc~tdf_kpt->module~w90_constants module~w90_postw90_types w90_postw90_types proc~tdf_kpt->module~w90_postw90_types module~w90_spin w90_spin proc~tdf_kpt->module~w90_spin module~w90_types w90_types proc~tdf_kpt->module~w90_types module~w90_utility w90_utility proc~tdf_kpt->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 TDF of a single k point

This routine does not use the adaptive smearing; in fact, for non-zero temperatures one often doesn't even need to smear. It simply uses a standard smearing as defined by the variables pw90_boltzwann_TDF_smr_fixed_en_width and pw90_boltzwann_TDF_smr_index

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

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/cell_volume if we want to calculate 1/(2*pi)^3 * \int_{BZ} d^3 k The only factor that is included INSIDE this routine is the spin degeneracy factor (=2 if spinors is .false., =1 if spinors is .true.) The EnergyArray is assumed to be evenly spaced (and the energy spacing is taken from EnergyArray(2)-EnergyArray(1)) The routine is assuming that EnergyArray has at least two elements. The meaning of the three indices of the TDF_k array is different with respect to those of the dos_k array returned by the dos_get_k routine The TDF_k array must have dimensions 6 * size(EnergyArray) * ndim, where ndim=1 if spin_decomp==false, or ndim=3 if spin_decomp==true. This is not checked.

Arguments

Type IntentOptional Attributes Name
type(pw90_boltzwann_type), intent(in) :: pw90_boltzwann
type(ws_region_type), intent(in) :: ws_region
type(pw90_spin_mod_type), intent(in) :: pw90_spin
type(wannier_data_type), intent(in) :: wannier_data
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(:,:,:,:)
real(kind=dp), intent(in) :: deleig_k(:,:)

array with the band derivatives at the given k point (in eV * angstrom / (2pi) as internally given by the code) already corrected in case of degeneracies, as returned by the wham_get_deleig_a routine
real(kind=dp), intent(in) :: eig_k(:)

array with the eigenvalues at the given k point (in eV)
real(kind=dp), intent(in) :: EnergyArray(:)

array with the energy grid on which to calculate the DOS (in eV) It must have at least two elements
real(kind=dp), intent(in) :: kpt(3)

the three coordinates of the k point vector whose DOS contribution we want to calculate (in relative coordinates)
real(kind=dp), intent(in) :: real_lattice(3,3)
real(kind=dp), intent(out) :: TDF_k(:,:,:)

TDF_k array in which the contribution is stored. Three dimensions: TDF_k(ij, energyidx, spinidx), where: - ij indexes the components of the TDF (symmetric) tensor (1=XX, 2=XY, ...); see the global constants defined in the module - 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
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: num_wann
integer, intent(in) :: num_elec_per_state
logical, intent(in) :: spin_decomp
type(w90_error_type), intent(out), allocatable :: error
type(w90_comm_type), intent(in) :: comm

Calls

proc~~tdf_kpt~~CallsGraph proc~tdf_kpt TDF_kpt proc~spin_get_nk spin_get_nk proc~tdf_kpt->proc~spin_get_nk proc~utility_w0gauss utility_w0gauss proc~tdf_kpt->proc~utility_w0gauss 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~set_error_input set_error_input proc~utility_w0gauss->proc~set_error_input 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~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~~tdf_kpt~~CalledByGraph proc~tdf_kpt TDF_kpt proc~calctdfanddos calcTDFandDOS proc~calctdfanddos->proc~tdf_kpt proc~boltzwann_main boltzwann_main proc~boltzwann_main->proc~calctdfanddos program~postw90 postw90 program~postw90->proc~boltzwann_main

Source Code

  subroutine TDF_kpt(pw90_boltzwann, ws_region, pw90_spin, wannier_data, ws_distance, &
                     wigner_seitz, HH_R, SS_R, deleig_k, eig_k, EnergyArray, kpt, real_lattice, &
                     TDF_k, mp_grid, num_wann, num_elec_per_state, spin_decomp, error, comm)
    !================================================!
    !! This subroutine calculates the contribution to the TDF of a single k point
    !!
    !!  This routine does not use the adaptive smearing; in fact, for non-zero temperatures
    !!       one often doesn't even need to smear. It simply uses a standard smearing as defined by
    !!       the variables pw90_boltzwann_TDF_smr_fixed_en_width and pw90_boltzwann_TDF_smr_index
    !!
    !! still to do: adapt spin_get_nk to read in input the UU rotation matrix
    !!
    !! 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/cell_volume
    !!       if we want to calculate 1/(2*pi)^3 * \int_{BZ} d^3 k
    !! The only factor that is included INSIDE this routine is the spin degeneracy
    !!       factor (=2 if spinors is .false., =1 if spinors is .true.)
    !! The EnergyArray is assumed to be evenly spaced (and the energy spacing
    !!       is taken from EnergyArray(2)-EnergyArray(1))
    !! The routine is assuming that EnergyArray has at least two elements.
    !!  The meaning of the three indices of the TDF_k array is different with respect to
    !!       those of the dos_k array returned by the dos_get_k routine
    !! The TDF_k array must have dimensions 6 * size(EnergyArray) * ndim, where
    !!       ndim=1 if spin_decomp==false, or ndim=3 if spin_decomp==true. This is not checked.
    !!
    !================================================!

    use w90_constants, only: dp, smearing_cutoff, min_smearing_binwidth_ratio
    use w90_utility, only: utility_w0gauss
    use w90_types, only: print_output_type, wannier_data_type, ws_region_type, ws_distance_type
    use w90_postw90_types, only: pw90_boltzwann_type, pw90_spin_mod_type, wigner_seitz_type
    use w90_spin, only: spin_get_nk
    use w90_utility, only: utility_w0gauss
    use w90_comms, only: w90_comm_type

    implicit none

    ! arguments
    type(pw90_boltzwann_type), intent(in) :: pw90_boltzwann
    type(ws_region_type), intent(in) :: ws_region
    type(pw90_spin_mod_type), intent(in) :: pw90_spin
    type(wannier_data_type), intent(in) :: wannier_data
    type(ws_distance_type), intent(inout) :: ws_distance
    type(wigner_seitz_type), intent(in) :: wigner_seitz
    type(w90_error_type), allocatable, intent(out) :: error
    type(w90_comm_type), intent(in) :: comm

    integer, intent(in) :: num_wann
    integer, intent(in) :: mp_grid(3)

    real(kind=dp), intent(in) :: kpt(3)
    !! the three coordinates of the k point vector whose DOS contribution we
    !! want to calculate (in relative coordinates)
    real(kind=dp), intent(in) :: EnergyArray(:)
    !! array with the energy grid on which to calculate the DOS (in eV)
    !! It must have at least two elements
    real(kind=dp), intent(in) :: eig_k(:)
    !! array with the eigenvalues at the given k point (in eV)
    real(kind=dp), intent(in) :: deleig_k(:, :)
    !! array with the band derivatives at the given k point
    !! (in eV * angstrom / (2pi) as internally given by the code)
    !! already corrected in case of degeneracies, as returned by the
    !! wham_get_deleig_a routine
    real(kind=dp), intent(out) :: TDF_k(:, :, :)
    !! TDF_k array in which the contribution is stored. Three dimensions:
    !! TDF_k(ij, energyidx, spinidx), where:
    !! - ij indexes the components of the TDF (symmetric) tensor (1=XX, 2=XY, ...);
    !!   see the global constants defined in the module
    !! - 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
    real(kind=dp), intent(in) :: real_lattice(3, 3)

    complex(kind=dp), allocatable, intent(inout) :: HH_R(:, :, :) !  <0n|r|Rm>
    complex(kind=dp), allocatable, intent(inout) :: SS_R(:, :, :, :) ! <0n|sigma_x,y,z|Rm>

    logical, intent(in) :: spin_decomp
    integer, intent(in) :: num_elec_per_state

    ! local variables
    real(kind=dp) :: smear, 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 :: BandIdx, loop_f, min_f, max_f
    logical :: DoSmearing

    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)

    TDF_k = 0.0_dp
    do BandIdx = 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(BandIdx))/2.0_dp ! |alpha|^2
        ! Contribution to spin-down DOS
        beta_sq = 1.0_dp - alpha_sq ! |beta|^2 = 1 - |alpha|^2
      end if

      ! Do not use an adaptive smearing here, it would require the knowledge of second derivatives
      ! And typically, when working at not too small temperatures, smearing is not needed

      ! Faster optimization: I precalculate the indices
      ! Value of the smearing in eV; default = 0 eV, i.e. no smearing
      smear = pw90_boltzwann%tdf_smearing%fixed_width
      if (smear/binwidth < min_smearing_binwidth_ratio) then
        min_f = max(nint((eig_k(BandIdx) - EnergyArray(1))/ &
                         (EnergyArray(size(EnergyArray)) - EnergyArray(1)) &
                         *real(size(EnergyArray) - 1, kind=dp)) + 1, 1)
        max_f = min(nint((eig_k(BandIdx) - 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(BandIdx) - smearing_cutoff*smear - EnergyArray(1))/ &
                         (EnergyArray(size(EnergyArray)) - EnergyArray(1)) &
                         *real(size(EnergyArray) - 1, kind=dp)) + 1, 1)
        max_f = min(nint((eig_k(BandIdx) + smearing_cutoff*smear - 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
        if (DoSmearing) then
          arg = (EnergyArray(loop_f) - eig_k(BandIdx))/smear
          rdum = utility_w0gauss(arg, pw90_boltzwann%tdf_smearing%type_index, error, comm)/smear
          if (allocated(error)) return
        else
          rdum = 1._dp/(EnergyArray(2) - EnergyArray(1))
        end if
        !
        ! Contribution to total DOS
        !
        TDF_k(XX, loop_f, 1) = TDF_k(XX, loop_f, 1) + rdum* &
                               r_num_elec_per_state*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 1)
        TDF_k(XY, loop_f, 1) = TDF_k(XY, loop_f, 1) + rdum* &
                               r_num_elec_per_state*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 2)
        TDF_k(YY, loop_f, 1) = TDF_k(YY, loop_f, 1) + rdum* &
                               r_num_elec_per_state*deleig_k(BandIdx, 2)*deleig_k(BandIdx, 2)
        TDF_k(XZ, loop_f, 1) = TDF_k(XZ, loop_f, 1) + rdum* &
                               r_num_elec_per_state*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 3)
        TDF_k(YZ, loop_f, 1) = TDF_k(YZ, loop_f, 1) + rdum* &
                               r_num_elec_per_state*deleig_k(BandIdx, 2)*deleig_k(BandIdx, 3)
        TDF_k(ZZ, loop_f, 1) = TDF_k(ZZ, loop_f, 1) + rdum* &
                               r_num_elec_per_state*deleig_k(BandIdx, 3)*deleig_k(BandIdx, 3)

        ! 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
          TDF_k(XX, loop_f, 2) = TDF_k(XX, loop_f, 2) + rdum* &
                                 alpha_sq*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 1)
          TDF_k(XY, loop_f, 2) = TDF_k(XY, loop_f, 2) + rdum* &
                                 alpha_sq*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 2)
          TDF_k(YY, loop_f, 2) = TDF_k(YY, loop_f, 2) + rdum* &
                                 alpha_sq*deleig_k(BandIdx, 2)*deleig_k(BandIdx, 2)
          TDF_k(XZ, loop_f, 2) = TDF_k(XZ, loop_f, 2) + rdum* &
                                 alpha_sq*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 3)
          TDF_k(YZ, loop_f, 2) = TDF_k(YZ, loop_f, 2) + rdum* &
                                 alpha_sq*deleig_k(BandIdx, 2)*deleig_k(BandIdx, 3)
          TDF_k(ZZ, loop_f, 2) = TDF_k(ZZ, loop_f, 2) + rdum* &
                                 alpha_sq*deleig_k(BandIdx, 3)*deleig_k(BandIdx, 3)
          ! Spin-down contribution
          TDF_k(XX, loop_f, 3) = TDF_k(XX, loop_f, 3) + rdum* &
                                 beta_sq*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 1)
          TDF_k(XY, loop_f, 3) = TDF_k(XY, loop_f, 3) + rdum* &
                                 beta_sq*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 2)
          TDF_k(YY, loop_f, 3) = TDF_k(YY, loop_f, 3) + rdum* &
                                 beta_sq*deleig_k(BandIdx, 2)*deleig_k(BandIdx, 2)
          TDF_k(XZ, loop_f, 3) = TDF_k(XZ, loop_f, 3) + rdum* &
                                 beta_sq*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 3)
          TDF_k(YZ, loop_f, 3) = TDF_k(YZ, loop_f, 3) + rdum* &
                                 beta_sq*deleig_k(BandIdx, 2)*deleig_k(BandIdx, 3)
          TDF_k(ZZ, loop_f, 3) = TDF_k(ZZ, loop_f, 3) + rdum* &
                                 beta_sq*deleig_k(BandIdx, 3)*deleig_k(BandIdx, 3)
        end if
      end do
    end do !loop over bands

    ! I multiply it here, since I am assuming a constant relaxation time, independent of the band index
    ! (actually, it is also independent of k)
    TDF_k = TDF_k*pw90_boltzwann%relax_time

  end subroutine TDF_kpt