berry_get_imfgh_klist Subroutine

public subroutine berry_get_imfgh_klist(dis_manifold, fermi_energy_list, kpt_latt, ws_region, print_output, wannier_data, ws_distance, wigner_seitz, AA_R, BB_R, CC_R, HH_R, u_matrix, v_matrix, eigval, kpt, real_lattice, scissors_shift, mp_grid, fermi_n, num_bands, num_kpts, num_wann, num_valence_bands, effective_model, have_disentangled, seedname, stdout, timer, error, comm, imf_k_list, img_k_list, imh_k_list, occ, ladpt)

Uses

  • proc~~berry_get_imfgh_klist~~UsesGraph proc~berry_get_imfgh_klist berry_get_imfgh_klist module~w90_comms w90_comms proc~berry_get_imfgh_klist->module~w90_comms module~w90_constants w90_constants proc~berry_get_imfgh_klist->module~w90_constants module~w90_postw90_common w90_postw90_common proc~berry_get_imfgh_klist->module~w90_postw90_common module~w90_postw90_types w90_postw90_types proc~berry_get_imfgh_klist->module~w90_postw90_types module~w90_types w90_types proc~berry_get_imfgh_klist->module~w90_types module~w90_utility w90_utility proc~berry_get_imfgh_klist->module~w90_utility module~w90_wan_ham w90_wan_ham proc~berry_get_imfgh_klist->module~w90_wan_ham module~w90_comms->module~w90_constants module~w90_error_base w90_error_base module~w90_comms->module~w90_error_base module~w90_postw90_common->module~w90_constants module~w90_error w90_error module~w90_postw90_common->module~w90_error module~w90_postw90_types->module~w90_comms module~w90_postw90_types->module~w90_constants module~w90_types->module~w90_constants module~w90_utility->module~w90_comms module~w90_utility->module~w90_constants module~w90_wan_ham->module~w90_constants module~w90_wan_ham->module~w90_error module~w90_error->module~w90_comms module~w90_error->module~w90_error_base

Calculates the three quantities needed for the orbital magnetization:

  • -2Im[f(k)] [Eq.33 CTVR06, Eq.6 LVTS12]
  • -2Im[g(k)] [Eq.34 CTVR06, Eq.7 LVTS12]
  • -2Im[h(k)] [Eq.35 CTVR06, Eq.8 LVTS12] They are calculated together (to reduce the number of Fourier calls) for a list of Fermi energies, and stored in axial-vector form.

Arguments

Type IntentOptional Attributes Name
type(dis_manifold_type), intent(in) :: dis_manifold
real(kind=dp), intent(in), allocatable :: fermi_energy_list(:)
real(kind=dp), intent(in) :: kpt_latt(:,:)
type(ws_region_type), intent(in) :: ws_region
type(print_output_type), intent(in) :: print_output
type(wannier_data_type), intent(in) :: wannier_data
type(ws_distance_type), intent(inout) :: ws_distance
type(wigner_seitz_type), intent(inout) :: wigner_seitz
complex(kind=dp), intent(inout), allocatable :: AA_R(:,:,:,:)
complex(kind=dp), intent(inout), allocatable :: BB_R(:,:,:,:)
complex(kind=dp), intent(inout), allocatable :: CC_R(:,:,:,:,:)
complex(kind=dp), intent(inout), allocatable :: HH_R(:,:,:)
complex(kind=dp), intent(in) :: u_matrix(:,:,:)
complex(kind=dp), intent(in) :: v_matrix(:,:,:)
real(kind=dp), intent(in) :: eigval(:,:)
real(kind=dp), intent(in) :: kpt(3)
real(kind=dp), intent(in) :: real_lattice(3,3)
real(kind=dp), intent(in) :: scissors_shift
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: fermi_n
integer, intent(in) :: num_bands
integer, intent(in) :: num_kpts
integer, intent(in) :: num_wann
integer, intent(in) :: num_valence_bands
logical, intent(in) :: effective_model
logical, intent(in) :: have_disentangled
character(len=50), intent(in) :: seedname
integer, intent(in) :: stdout
type(timer_list_type), intent(inout) :: timer
type(w90_error_type), intent(out), allocatable :: error
type(w90_comm_type), intent(in) :: comm
real(kind=dp), intent(out), optional :: imf_k_list(:,:,:)
real(kind=dp), intent(out), optional :: img_k_list(:,:,:)
real(kind=dp), intent(out), optional :: imh_k_list(:,:,:)
real(kind=dp), intent(in), optional :: occ(:)
logical, intent(in), optional :: ladpt(:)

Calls

proc~~berry_get_imfgh_klist~~CallsGraph proc~berry_get_imfgh_klist berry_get_imfgh_klist proc~pw90common_fourier_r_to_k pw90common_fourier_R_to_k proc~berry_get_imfgh_klist->proc~pw90common_fourier_r_to_k proc~pw90common_fourier_r_to_k_vec pw90common_fourier_R_to_k_vec proc~berry_get_imfgh_klist->proc~pw90common_fourier_r_to_k_vec proc~utility_im_tr_prod utility_im_tr_prod proc~berry_get_imfgh_klist->proc~utility_im_tr_prod proc~utility_re_tr_prod utility_re_tr_prod proc~berry_get_imfgh_klist->proc~utility_re_tr_prod proc~utility_zgemm_new utility_zgemm_new proc~berry_get_imfgh_klist->proc~utility_zgemm_new proc~wham_get_eig_uu_hh_jjlist wham_get_eig_UU_HH_JJlist proc~berry_get_imfgh_klist->proc~wham_get_eig_uu_hh_jjlist proc~wham_get_occ_mat_list wham_get_occ_mat_list proc~berry_get_imfgh_klist->proc~wham_get_occ_mat_list zgemm zgemm proc~utility_zgemm_new->zgemm proc~get_hh_r get_HH_R proc~wham_get_eig_uu_hh_jjlist->proc~get_hh_r proc~pw90common_fourier_r_to_k_new pw90common_fourier_R_to_k_new proc~wham_get_eig_uu_hh_jjlist->proc~pw90common_fourier_r_to_k_new proc~utility_diagonalize utility_diagonalize proc~wham_get_eig_uu_hh_jjlist->proc~utility_diagonalize proc~wham_get_jjp_jjm_list wham_get_JJp_JJm_list proc~wham_get_eig_uu_hh_jjlist->proc~wham_get_jjp_jjm_list proc~pw90common_get_occ pw90common_get_occ proc~wham_get_occ_mat_list->proc~pw90common_get_occ proc~set_error_input set_error_input proc~wham_get_occ_mat_list->proc~set_error_input proc~get_hh_r->proc~set_error_input interface~comms_bcast comms_bcast proc~get_hh_r->interface~comms_bcast proc~fourier_q_to_r fourier_q_to_R proc~get_hh_r->proc~fourier_q_to_r proc~get_win_min get_win_min proc~get_hh_r->proc~get_win_min proc~io_stopwatch_start io_stopwatch_start proc~get_hh_r->proc~io_stopwatch_start proc~io_stopwatch_stop io_stopwatch_stop proc~get_hh_r->proc~io_stopwatch_stop proc~mpirank mpirank proc~get_hh_r->proc~mpirank proc~operator_wigner_setup operator_wigner_setup proc~get_hh_r->proc~operator_wigner_setup proc~set_error_fatal set_error_fatal proc~get_hh_r->proc~set_error_fatal proc~set_error_file set_error_file proc~get_hh_r->proc~set_error_file 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~utility_diagonalize->proc~set_error_fatal zhpevx zhpevx proc~utility_diagonalize->zhpevx proc~utility_rotate_new utility_rotate_new proc~wham_get_jjp_jjm_list->proc~utility_rotate_new proc~comms_bcast_char comms_bcast_char interface~comms_bcast->proc~comms_bcast_char proc~comms_bcast_cmplx comms_bcast_cmplx interface~comms_bcast->proc~comms_bcast_cmplx proc~comms_bcast_int comms_bcast_int interface~comms_bcast->proc~comms_bcast_int proc~comms_bcast_logical comms_bcast_logical interface~comms_bcast->proc~comms_bcast_logical proc~comms_bcast_real comms_bcast_real interface~comms_bcast->proc~comms_bcast_real proc~set_error_fatal->proc~comms_sync_error proc~set_error_fatal->proc~set_base_error proc~set_error_file->proc~comms_sync_error proc~set_error_file->proc~set_base_error proc~utility_rotate_new->proc~utility_zgemm_new proc~comms_bcast_char->proc~comms_sync_error proc~comms_no_sync_bcast_char comms_no_sync_bcast_char proc~comms_bcast_char->proc~comms_no_sync_bcast_char proc~comms_bcast_cmplx->proc~comms_sync_error proc~comms_no_sync_bcast_cmplx comms_no_sync_bcast_cmplx proc~comms_bcast_cmplx->proc~comms_no_sync_bcast_cmplx proc~comms_bcast_int->proc~comms_sync_error proc~comms_no_sync_bcast_int comms_no_sync_bcast_int proc~comms_bcast_int->proc~comms_no_sync_bcast_int proc~comms_bcast_logical->proc~comms_sync_error proc~comms_no_sync_bcast_logical comms_no_sync_bcast_logical proc~comms_bcast_logical->proc~comms_no_sync_bcast_logical proc~comms_bcast_real->proc~comms_sync_error proc~comms_no_sync_bcast_real comms_no_sync_bcast_real proc~comms_bcast_real->proc~comms_no_sync_bcast_real

Called by

proc~~berry_get_imfgh_klist~~CalledByGraph proc~berry_get_imfgh_klist berry_get_imfgh_klist proc~berry_get_imf_klist berry_get_imf_klist proc~berry_get_imf_klist->proc~berry_get_imfgh_klist proc~berry_main berry_main proc~berry_main->proc~berry_get_imfgh_klist proc~berry_main->proc~berry_get_imf_klist proc~gyrotropic_get_k_list gyrotropic_get_k_list proc~gyrotropic_get_k_list->proc~berry_get_imfgh_klist proc~gyrotropic_get_k_list->proc~berry_get_imf_klist proc~k_path k_path proc~k_path->proc~berry_get_imfgh_klist proc~k_path->proc~berry_get_imf_klist proc~k_slice k_slice proc~k_slice->proc~berry_get_imfgh_klist proc~k_slice->proc~berry_get_imf_klist proc~gyrotropic_main gyrotropic_main proc~gyrotropic_main->proc~gyrotropic_get_k_list program~postw90 postw90 program~postw90->proc~berry_main program~postw90->proc~k_path program~postw90->proc~k_slice program~postw90->proc~gyrotropic_main

Source Code

  subroutine berry_get_imfgh_klist(dis_manifold, fermi_energy_list, kpt_latt, ws_region, print_output, &
                                   wannier_data, ws_distance, wigner_seitz, AA_R, BB_R, CC_R, HH_R, &
                                   u_matrix, v_matrix, eigval, kpt, real_lattice, &
                                   scissors_shift, mp_grid, fermi_n, num_bands, num_kpts, &
                                   num_wann, num_valence_bands, effective_model, &
                                   have_disentangled, seedname, stdout, timer, error, comm, &
                                   imf_k_list, img_k_list, imh_k_list, occ, ladpt)
    !================================================!
    !
    !! Calculates the three quantities needed for the orbital
    !! magnetization:
    !!
    !! * -2Im[f(k)] [Eq.33 CTVR06, Eq.6 LVTS12]
    !! * -2Im[g(k)] [Eq.34 CTVR06, Eq.7 LVTS12]
    !! * -2Im[h(k)] [Eq.35 CTVR06, Eq.8 LVTS12]
    !! They are calculated together (to reduce the number of
    !! Fourier calls) for a list of Fermi energies, and stored
    !! in axial-vector form.
    !
    ! The two optional output parameters 'imh_k_list' and
    ! 'img_k_list' are only calculated if both of them are
    ! present.
    !
    !================================================!

    use w90_comms, only: w90_comm_type, mpirank
    use w90_constants, only: dp, cmplx_i
    use w90_types, only: print_output_type, wannier_data_type, &
                         dis_manifold_type, kmesh_info_type, ws_region_type, ws_distance_type, timer_list_type
    use w90_postw90_common, only: pw90common_fourier_R_to_k_vec, pw90common_fourier_R_to_k
    use w90_postw90_types, only: wigner_seitz_type
    use w90_utility, only: utility_re_tr_prod, utility_im_tr_prod, utility_zgemm_new
    use w90_wan_ham, only: wham_get_eig_UU_HH_JJlist, wham_get_occ_mat_list

    implicit none

    ! arguments
    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(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, fermi_n
    integer, intent(in) :: mp_grid(3)
    integer, intent(in) :: stdout

    real(kind=dp), intent(in) :: kpt(3)
    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(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) :: BB_R(:, :, :, :) ! <0|H(r-R)|R>
    complex(kind=dp), allocatable, intent(inout) :: CC_R(:, :, :, :, :) ! <0|r_alpha.H(r-R)_beta|R>

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

    real(kind=dp), intent(out), optional :: imf_k_list(:, :, :)
    real(kind=dp), intent(out), optional :: img_k_list(:, :, :)
    real(kind=dp), intent(out), optional :: imh_k_list(:, :, :)
    real(kind=dp), intent(in), optional :: occ(:)
    logical, intent(in) :: effective_model
    logical, intent(in), optional :: ladpt(:)

    ! local variables
    complex(kind=dp), allocatable :: HH(:, :)
    complex(kind=dp), allocatable :: UU(:, :)
    complex(kind=dp), allocatable :: f_list(:, :, :)
    complex(kind=dp), allocatable :: g_list(:, :, :)
    complex(kind=dp), allocatable :: AA(:, :, :)
    complex(kind=dp), allocatable :: BB(:, :, :)
    complex(kind=dp), allocatable :: CC(:, :, :, :)
    complex(kind=dp), allocatable :: OOmega(:, :, :)
    complex(kind=dp), allocatable :: JJp_list(:, :, :, :)
    complex(kind=dp), allocatable :: JJm_list(:, :, :, :)
    ! Temporary space for matrix products
    complex(kind=dp), allocatable :: tmp(:, :, :)

    real(kind=dp) :: eig(num_wann)
    real(kind=dp) :: s

    integer :: i, j, ife, nfermi_loc

    logical :: todo(fermi_n)

    if (present(occ)) then
      nfermi_loc = 1
    else
      nfermi_loc = fermi_n
    end if

    if (present(ladpt)) then
      todo = ladpt
    else
      todo = .true.
    end if

    allocate (HH(num_wann, num_wann))
    allocate (UU(num_wann, num_wann))
    allocate (f_list(num_wann, num_wann, nfermi_loc))
    allocate (g_list(num_wann, num_wann, nfermi_loc))
    allocate (JJp_list(num_wann, num_wann, nfermi_loc, 3))
    allocate (JJm_list(num_wann, num_wann, nfermi_loc, 3))
    allocate (AA(num_wann, num_wann, 3))
    allocate (OOmega(num_wann, num_wann, 3))

    ! Gather W-gauge matrix objects
    !

    if (present(occ)) then
      call wham_get_eig_UU_HH_JJlist(dis_manifold, fermi_energy_list, kpt_latt, ws_region, &
                                     print_output, wannier_data, ws_distance, wigner_seitz, HH, &
                                     HH_R, JJm_list, JJp_list, u_matrix, UU, v_matrix, 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, occ=occ)
      if (allocated(error)) return
      call wham_get_occ_mat_list(fermi_energy_list, f_list, g_list, UU, num_wann, error, comm, &
                                 occ=occ)
      if (allocated(error)) return

    else
      call wham_get_eig_UU_HH_JJlist(dis_manifold, fermi_energy_list, kpt_latt, ws_region, &
                                     print_output, wannier_data, ws_distance, wigner_seitz, HH, &
                                     HH_R, JJm_list, JJp_list, u_matrix, UU, v_matrix, 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 wham_get_occ_mat_list(fermi_energy_list, f_list, g_list, UU, num_wann, error, comm, &
                                 eig=eig)
      if (allocated(error)) return

    end if

    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, OO_pseudo=OOmega)
    if (allocated(error)) return

    if (present(imf_k_list)) then
      ! Trace formula for -2Im[f], Eq.(51) LVTS12
      !
      do ife = 1, nfermi_loc
        if (todo(ife)) then
          do i = 1, 3
            !
            ! J0 term (Omega_bar term of WYSV06)
            imf_k_list(1, i, ife) = &
              utility_re_tr_prod(f_list(:, :, ife), OOmega(:, :, i))
            !
            ! J1 term (DA term of WYSV06)
            imf_k_list(2, i, ife) = -2.0_dp* &
                                    ( &
                                    utility_im_tr_prod(AA(:, :, alpha_A(i)), JJp_list(:, :, ife, beta_A(i))) &
                                    + utility_im_tr_prod(JJm_list(:, :, ife, alpha_A(i)), AA(:, :, beta_A(i))) &
                                    )
            !
            ! J2 term (DD of WYSV06)
            imf_k_list(3, i, ife) = -2.0_dp* &
                                    utility_im_tr_prod(JJm_list(:, :, ife, alpha_A(i)), JJp_list(:, :, ife, beta_A(i)))
          end do
        end if
      end do
    end if

    if (present(img_k_list)) img_k_list = 0.0_dp
    if (present(imh_k_list)) imh_k_list = 0.0_dp

    if (present(img_k_list) .and. present(imh_k_list)) then
      allocate (BB(num_wann, num_wann, 3))
      allocate (CC(num_wann, num_wann, 3, 3))

      allocate (tmp(num_wann, num_wann, 5))
      ! tmp(:,:,1:3) ... not dependent on inner loop variables
      ! tmp(:,:,1) ..... HH . AA(:,:,alpha_A(i))
      ! tmp(:,:,2) ..... LLambda_ij [Eq. (37) LVTS12] expressed as a pseudovector
      ! tmp(:,:,3) ..... HH . OOmega(:,:,i)
      ! tmp(:,:,4:5) ... working matrices for matrix products of inner loop

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

      do j = 1, 3
        do i = 1, j
          call pw90common_fourier_R_to_k(ws_region, wannier_data, ws_distance, wigner_seitz, &
                                         CC(:, :, i, j), CC_R(:, :, :, i, j), kpt, real_lattice, &
                                         mp_grid, 0, num_wann, error, comm)
          if (allocated(error)) return

          CC(:, :, j, i) = conjg(transpose(CC(:, :, i, j)))
        end do
      end do

      ! Trace formula for -2Im[g], Eq.(66) LVTS12
      ! Trace formula for -2Im[h], Eq.(56) LVTS12
      !
      do i = 1, 3
        call utility_zgemm_new(HH, AA(:, :, alpha_A(i)), tmp(:, :, 1))
        call utility_zgemm_new(HH, OOmega(:, :, i), tmp(:, :, 3))
        !
        ! LLambda_ij [Eq. (37) LVTS12] expressed as a pseudovector
        tmp(:, :, 2) = cmplx_i*(CC(:, :, alpha_A(i), beta_A(i)) &
                                - conjg(transpose(CC(:, :, alpha_A(i), beta_A(i)))))

        do ife = 1, nfermi_loc
          !
          ! J0 terms for -2Im[g] and -2Im[h]
          !
          ! tmp(:,:,5) = HH . AA(:,:,alpha_A(i)) . f_list(:,:,ife) . AA(:,:,beta_A(i))
          call utility_zgemm_new(tmp(:, :, 1), f_list(:, :, ife), tmp(:, :, 4))
          call utility_zgemm_new(tmp(:, :, 4), AA(:, :, beta_A(i)), tmp(:, :, 5))

          s = 2.0_dp*utility_im_tr_prod(f_list(:, :, ife), tmp(:, :, 5))
          img_k_list(1, i, ife) = utility_re_tr_prod(f_list(:, :, ife), tmp(:, :, 2)) - s
          imh_k_list(1, i, ife) = utility_re_tr_prod(f_list(:, :, ife), tmp(:, :, 3)) + s

          !
          ! J1 terms for -2Im[g] and -2Im[h]
          !
          ! tmp(:,:,1) = HH . AA(:,:,alpha_A(i))
          ! tmp(:,:,4) = HH . JJm_list(:,:,ife,alpha_A(i))
          call utility_zgemm_new(HH, JJm_list(:, :, ife, alpha_A(i)), tmp(:, :, 4))

          img_k_list(2, i, ife) = -2.0_dp* &
                                  ( &
                                  utility_im_tr_prod(JJm_list(:, :, ife, alpha_A(i)), BB(:, :, beta_A(i))) &
                                  - utility_im_tr_prod(JJm_list(:, :, ife, beta_A(i)), BB(:, :, alpha_A(i))) &
                                  )
          imh_k_list(2, i, ife) = -2.0_dp* &
                                  ( &
                                  utility_im_tr_prod(tmp(:, :, 1), JJp_list(:, :, ife, beta_A(i))) &
                                  + utility_im_tr_prod(tmp(:, :, 4), AA(:, :, beta_A(i))) &
                                  )

          !
          ! J2 terms for -2Im[g] and -2Im[h]
          !
          ! tmp(:,:,4) = JJm_list(:,:,ife,alpha_A(i)) . HH
          ! tmp(:,:,5) = HH . JJm_list(:,:,ife,alpha_A(i))
          call utility_zgemm_new(JJm_list(:, :, ife, alpha_A(i)), HH, tmp(:, :, 4))
          call utility_zgemm_new(HH, JJm_list(:, :, ife, alpha_A(i)), tmp(:, :, 5))

          img_k_list(3, i, ife) = -2.0_dp* &
                                  utility_im_tr_prod(tmp(:, :, 4), JJp_list(:, :, ife, beta_A(i)))
          imh_k_list(3, i, ife) = -2.0_dp* &
                                  utility_im_tr_prod(tmp(:, :, 5), JJp_list(:, :, ife, beta_A(i)))
        end do
      end do
      deallocate (tmp)
    end if

  end subroutine berry_get_imfgh_klist