get_SHC_R – Wannier90

public subroutine get_SHC_R(dis_manifold, kmesh_info, kpt_latt, print_output, pw90_oper_read, pw90_spin_hall, SH_R, SHR_R, SR_R, v_matrix, eigval, scissors_shift, wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, num_wann, num_valence_bands, have_disentangled, seedname, stdout, timer, error, comm)

Uses

- w90_postw90_types
- w90_types
Help

Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.

Compute several matrices for spin Hall conductivity SR_R = <0n|sigma_{x,y,z}.(r-R)alpha|Rm> SHR_R = <0n|sigma.H.(r-R)alpha|Rm> SH_R = <0n|sigma.H|Rm>

Arguments

Type	Intent	Attributes		Name
type(dis_manifold_type),	intent(in)		::	dis_manifold
type(kmesh_info_type),	intent(in)		::	kmesh_info
real(kind=dp),	intent(in)		::	kpt_latt(:,:)
type(print_output_type),	intent(in)		::	print_output
type(pw90_oper_read_type),	intent(in)		::	pw90_oper_read
type(pw90_spin_hall_type),	intent(in)		::	pw90_spin_hall
complex(kind=dp),	intent(inout),	allocatable	::	SH_R(:,:,:,:)
complex(kind=dp),	intent(inout),	allocatable	::	SHR_R(:,:,:,:,:)
complex(kind=dp),	intent(inout),	allocatable	::	SR_R(:,:,:,:,:)
complex(kind=dp),	intent(in)		::	v_matrix(:,:,:)
real(kind=dp),	intent(in)		::	eigval(:,:)
real(kind=dp),	intent(in)		::	scissors_shift
type(wigner_seitz_type),	intent(in)		::	wigner_seitz
type(ws_distance_type),	intent(in)		::	ws_distance
type(ws_region_type),	intent(in)		::	ws_region
integer,	intent(in)		::	num_bands
integer,	intent(in)		::	num_kpts
integer,	intent(in)		::	num_wann
integer,	intent(in)		::	num_valence_bands
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

Called by

Help

Source Code

  subroutine get_SHC_R(dis_manifold, kmesh_info, kpt_latt, print_output, pw90_oper_read, &
                       pw90_spin_hall, SH_R, SHR_R, SR_R, v_matrix, eigval, scissors_shift, &
                       wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, num_wann, &
                       num_valence_bands, have_disentangled, seedname, stdout, timer, error, comm)
    !================================================
    !
    !! Compute several matrices for spin Hall conductivity
    !! SR_R  = <0n|sigma_{x,y,z}.(r-R)_alpha|Rm>
    !! SHR_R = <0n|sigma_{x,y,z}.H.(r-R)_alpha|Rm>
    !! SH_R  = <0n|sigma_{x,y,z}.H|Rm>
    !
    !================================================

    use w90_postw90_types, only: pw90_oper_read_type, pw90_spin_hall_type, wigner_seitz_type
    use w90_types, only: dis_manifold_type, kmesh_info_type, ws_distance_type, ws_region_type, &
                         print_output_type, timer_list_type

    implicit none

    ! arguments
    type(dis_manifold_type), intent(in) :: dis_manifold
    type(kmesh_info_type), intent(in) :: kmesh_info
    type(pw90_oper_read_type), intent(in) :: pw90_oper_read
    type(print_output_type), intent(in) :: print_output
    type(pw90_spin_hall_type), intent(in) :: pw90_spin_hall
    type(wigner_seitz_type), intent(in) :: wigner_seitz
    type(ws_distance_type), intent(in) :: ws_distance
    type(ws_region_type), intent(in) :: ws_region
    type(timer_list_type), intent(inout) :: timer
    type(w90_comm_type), intent(in) :: comm
    type(w90_error_type), allocatable, intent(out) :: error

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

    real(kind=dp), intent(in) :: eigval(:, :)
    real(kind=dp), intent(in) :: scissors_shift
    real(kind=dp), intent(in) :: kpt_latt(:, :)

    complex(kind=dp), intent(in) :: v_matrix(:, :, :)
    complex(kind=dp), allocatable, intent(inout) :: SR_R(:, :, :, :, :) ! <0n|sigma_x,y,z.(r-R)_alpha|Rm>
    complex(kind=dp), allocatable, intent(inout) :: SHR_R(:, :, :, :, :) ! <0n|sigma_x,y,z.H.(r-R)_alpha|Rm>
    complex(kind=dp), allocatable, intent(inout) :: SH_R(:, :, :, :) ! <0n|sigma_x,y,z.H|Rm>

    complex(kind=dp), allocatable :: SR_R_temp(:, :, :, :, :)
    complex(kind=dp), allocatable :: SHR_R_temp(:, :, :, :, :)
    complex(kind=dp), allocatable :: SH_R_temp(:, :, :, :)

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

    ! local variables
    complex(kind=dp), allocatable :: SR_q(:, :, :, :, :)
    complex(kind=dp), allocatable :: SHR_q(:, :, :, :, :)
    complex(kind=dp), allocatable :: SH_q(:, :, :, :)

    complex(kind=dp), allocatable :: S_o(:, :)
    complex(kind=dp), allocatable :: spn_o(:, :, :, :), spn_temp(:, :)
    complex(kind=dp), allocatable :: H_o(:, :, :)
    complex(kind=dp), allocatable :: SH_o(:, :, :, :)
    complex(kind=dp)              :: SM_o(num_bands, num_bands, 3)
    complex(kind=dp)              :: SHM_o(num_bands, num_bands, 3)

    complex(kind=dp)              :: SS_q(num_wann, num_wann, 3)
    complex(kind=dp)              :: SM_q(num_wann, num_wann, 3)
    complex(kind=dp)              :: SHM_q(num_wann, num_wann, 3)

    real(kind=dp)                 :: s_real, s_img
    integer                       :: spn_in, counter, ierr, s, is

    integer                       :: n, m, &
                                     ik, ik2, ik_prev, nn, inn, nnl, nnm, nnn, &
                                     idir, ncount, nn_count, mmn_in, &
                                     nb_tmp, nkp_tmp, nntot_tmp
    integer, allocatable          :: num_states(:)
    real(kind=dp)                 :: m_real, m_imag
    logical                       :: nn_found
    character(len=60)             :: header
    logical :: on_root = .false.

    if (mpirank(comm) == 0) on_root = .true.

    if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
      call io_stopwatch_start('get_oper: get_SHC_R', timer)

    allocate (SR_R_temp(num_wann, num_wann, wigner_seitz%nrpts, 3, 3))
    allocate (SHR_R_temp(num_wann, num_wann, wigner_seitz%nrpts, 3, 3))
    allocate (SH_R_temp(num_wann, num_wann, wigner_seitz%nrpts, 3))

    if (.not. allocated(SR_R)) then
      allocate (SR_R(num_wann, num_wann, wigner_seitz%nrpts_pw90, 3, 3))
    else
      if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
        call io_stopwatch_stop('get_oper: get_SHC_R', timer)
      return
    end if
    if (.not. allocated(SHR_R)) then
      allocate (SHR_R(num_wann, num_wann, wigner_seitz%nrpts_pw90, 3, 3))
    else
      if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
        call io_stopwatch_stop('get_oper: get_SHC_R', timer)
      return
    end if
    if (.not. allocated(SH_R)) then
      allocate (SH_R(num_wann, num_wann, wigner_seitz%nrpts_pw90, 3))
    else
      if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
        call io_stopwatch_stop('get_oper: get_SHC_R', timer)
      return
    end if

    ! start copying from get_SS_R, Junfeng Qiao
    ! read spn file
    if (on_root) then

      allocate (spn_o(num_bands, num_bands, num_kpts, 3))

      allocate (num_states(num_kpts))
      do ik = 1, num_kpts
        if (have_disentangled) then
          num_states(ik) = dis_manifold%ndimwin(ik)
        else
          num_states(ik) = num_wann
        end if
      end do

      ! Read from .spn file the original spin matrices <psi_nk|sigma_i|psi_mk>
      ! (sigma_i = Pauli matrix) between ab initio eigenstates
      !
      if (pw90_oper_read%spn_formatted) then
        open (newunit=spn_in, file=trim(seedname)//'.spn', form='formatted', &
              status='old', err=109)
        write (stdout, '(/a)', advance='no') &
          ' Reading spin matrices from '//trim(seedname)//'.spn in get_SHC_R : '
        read (spn_in, *, err=110, end=110) header
        write (stdout, '(a)') trim(header)
        read (spn_in, *, err=110, end=110) nb_tmp, nkp_tmp
      else
        open (newunit=spn_in, file=trim(seedname)//'.spn', form='unformatted', &
              status='old', err=109)
        write (stdout, '(/a)', advance='no') &
          ' Reading spin matrices from '//trim(seedname)//'.spn in get_SHC_R : '
        read (spn_in, err=110, end=110) header
        write (stdout, '(a)') trim(header)
        read (spn_in, err=110, end=110) nb_tmp, nkp_tmp
      end if
      if (nb_tmp .ne. num_bands) then
        call set_error_fatal(error, trim(seedname)//'.spn has wrong number of bands', comm)
        return
      end if
      if (nkp_tmp .ne. num_kpts) then
        call set_error_fatal(error, trim(seedname)//'.spn has wrong number of k-points', comm)
        return
      end if
      if (pw90_oper_read%spn_formatted) then
        do ik = 1, num_kpts
          do m = 1, num_bands
            do n = 1, m
              read (spn_in, *, err=110, end=110) s_real, s_img
              spn_o(n, m, ik, 1) = cmplx(s_real, s_img, dp)
              read (spn_in, *, err=110, end=110) s_real, s_img
              spn_o(n, m, ik, 2) = cmplx(s_real, s_img, dp)
              read (spn_in, *, err=110, end=110) s_real, s_img
              spn_o(n, m, ik, 3) = cmplx(s_real, s_img, dp)
              ! Read upper-triangular part, now build the rest
              spn_o(m, n, ik, 1) = conjg(spn_o(n, m, ik, 1))
              spn_o(m, n, ik, 2) = conjg(spn_o(n, m, ik, 2))
              spn_o(m, n, ik, 3) = conjg(spn_o(n, m, ik, 3))
            end do
          end do
        end do
      else
        allocate (spn_temp(3, (num_bands*(num_bands + 1))/2), stat=ierr)
        if (ierr /= 0) then
          call set_error_alloc(error, 'Error in allocating spm_temp in get_SHC_R', comm)
          return
        end if
        do ik = 1, num_kpts
          read (spn_in) ((spn_temp(s, m), s=1, 3), m=1, (num_bands*(num_bands + 1))/2)
          counter = 0
          do m = 1, num_bands
            do n = 1, m
              counter = counter + 1
              spn_o(n, m, ik, 1) = spn_temp(1, counter)
              spn_o(m, n, ik, 1) = conjg(spn_temp(1, counter))
              spn_o(n, m, ik, 2) = spn_temp(2, counter)
              spn_o(m, n, ik, 2) = conjg(spn_temp(2, counter))
              spn_o(n, m, ik, 3) = spn_temp(3, counter)
              spn_o(m, n, ik, 3) = conjg(spn_temp(3, counter))
            end do
          end do
        end do
        deallocate (spn_temp, stat=ierr)
        if (ierr /= 0) then
          call set_error_dealloc(error, 'Error in deallocating spm_temp in get_SHC_R', comm)
          return
        end if
      end if

      close (spn_in)

    end if !on_root
    ! end copying from get_SS_R, Junfeng Qiao

    ! start copying from get_HH_R, Junfeng Qiao
    ! Note this is different from get_HH_R, at here we need the
    ! original Hamiltonian to construct SHR_R, SH_R.
    if (on_root) then
      allocate (H_o(num_bands, num_bands, num_kpts))
      H_o = cmplx_0
      do ik = 1, num_kpts
        do m = 1, num_bands
          H_o(m, m, ik) = eigval(m, ik)
        end do
        ! scissors shift applied to the original Hamiltonian
        if (num_valence_bands > 0 .and. abs(scissors_shift) > 1.0e-7_dp) then
          do m = num_valence_bands + 1, num_bands
            H_o(m, m, ik) = H_o(m, m, ik) + scissors_shift
          end do
        else if (pw90_spin_hall%bandshift) then
          do m = pw90_spin_hall%bandshift_firstband, num_bands
            H_o(m, m, ik) = H_o(m, m, ik) + pw90_spin_hall%bandshift_energyshift
          end do
        end if
      end do
    end if !on_root
    ! end copying from get_HH_R, Junfeng Qiao

    ! start copying from get_AA_R, Junfeng Qiao
    ! read mmn file
    !
    if (on_root) then

      allocate (SR_q(num_wann, num_wann, num_kpts, 3, 3))
      allocate (SHR_q(num_wann, num_wann, num_kpts, 3, 3))
      allocate (SH_q(num_wann, num_wann, num_kpts, 3))
      allocate (S_o(num_bands, num_bands))

      open (newunit=mmn_in, file=trim(seedname)//'.mmn', &
            form='formatted', status='old', action='read', err=101)
      write (stdout, '(/a)', advance='no') &
        ' Reading overlaps from '//trim(seedname)//'.mmn in get_SHC_R   : '
      ! Read the comment line (header)
      read (mmn_in, '(a)', err=102, end=102) header
      write (stdout, '(a)') trim(header)
      ! Read the number of bands, k-points and nearest neighbours
      read (mmn_in, *, err=102, end=102) nb_tmp, nkp_tmp, nntot_tmp
      ! Checks
      if (nb_tmp .ne. num_bands) then
        call set_error_fatal(error, trim(seedname)//'.mmn has wrong number of bands', comm)
        return
      end if
      if (nkp_tmp .ne. num_kpts) then
        call set_error_fatal(error, trim(seedname)//'.mmn has wrong number of k-points', comm)
        return
      end if
      if (nntot_tmp .ne. kmesh_info%nntot) then
        call set_error_fatal(error, trim(seedname)//'.mmn has wrong number of nearest neighbours', comm)
        return
      end if

      SR_q = cmplx_0
      SHR_q = cmplx_0
      SH_q = cmplx_0
      ik_prev = 0

      ! QZYZ18 Eq.(48)
      allocate (SH_o(num_bands, num_bands, num_kpts, 3))
      SH_o = cmplx_0
      do ik = 1, num_kpts
        do is = 1, 3
          SH_o(:, :, ik, is) = matmul(spn_o(:, :, ik, is), H_o(:, :, ik))

          call get_gauge_overlap_matrix(num_bands, num_wann, eigval, v_matrix, dis_manifold, &
                                        ik, num_states(ik), ik, num_states(ik), &
                                        SH_o(:, :, ik, is), have_disentangled, SH_q(:, :, ik, is))
        end do
      end do

      ! Composite loop over k-points ik (outer loop) and neighbors ik2 (inner)
      do ncount = 1, num_kpts*kmesh_info%nntot
        !
        !Read from .mmn file the original overlap matrix
        ! S_o=<u_ik|u_ik2> between ab initio eigenstates
        !
        read (mmn_in, *, err=102, end=102) ik, ik2, nnl, nnm, nnn
        do n = 1, num_bands
          do m = 1, num_bands
            read (mmn_in, *, err=102, end=102) m_real, m_imag
            S_o(m, n) = cmplx(m_real, m_imag, kind=dp)
          end do
        end do
        !debug
        !OK
        !if(ik.ne.ik_prev .and.ik_prev.ne.0) then
        !   if(nn_count.ne.nntot)&
        !        write(stdout,*) 'something wrong in get_AA_R!'
        !endif
        !enddebug
        if (ik .ne. ik_prev) nn_count = 0
        nn = 0
        nn_found = .false.
        do inn = 1, kmesh_info%nntot
          if ((ik2 .eq. kmesh_info%nnlist(ik, inn)) .and. &
              (nnl .eq. kmesh_info%nncell(1, ik, inn)) .and. &
              (nnm .eq. kmesh_info%nncell(2, ik, inn)) .and. &
              (nnn .eq. kmesh_info%nncell(3, ik, inn))) then
            if (.not. nn_found) then
              nn_found = .true.
              nn = inn
            else
              call set_error_fatal(error, 'Error reading '//trim(seedname)//'.mmn.&
                   & More than one matching nearest neighbour found', comm)
              return
            end if
          end if
        end do
        if (nn .eq. 0) then
          write (stdout, '(/a,i8,2i5,i4,2x,3i3)') ' Error reading '//trim(seedname)//'.mmn:', &
            ncount, ik, ik2, nn, nnl, nnm, nnn
          call set_error_fatal(error, 'Neighbour not found', comm)
          return
        end if
        nn_count = nn_count + 1 !Check: can also be place after nn=inn (?)

        SM_o = cmplx_0
        SHM_o = cmplx_0
        SS_q = cmplx_0
        SM_q = cmplx_0
        SHM_q = cmplx_0
        do is = 1, 3
          ! QZYZ18 Eq.(50)
          SM_o(:, :, is) = matmul(spn_o(:, :, ik, is), S_o(:, :))
          ! QZYZ18 Eq.(51)
          SHM_o(:, :, is) = matmul(SH_o(:, :, ik, is), S_o(:, :))

          ! Transform to projected subspace, Wannier gauge
          !
          ! QZYZ18 Eq.(50)
          call get_gauge_overlap_matrix(num_bands, num_wann, eigval, v_matrix, dis_manifold, &
                                        ik, num_states(ik), ik, num_states(ik), &
                                        spn_o(:, :, ik, is), have_disentangled, SS_q(:, :, is))
          ! QZYZ18 Eq.(50)
          call get_gauge_overlap_matrix(num_bands, num_wann, eigval, v_matrix, dis_manifold, &
                                        ik, num_states(ik), kmesh_info%nnlist(ik, nn), &
                                        num_states(kmesh_info%nnlist(ik, nn)), SM_o(:, :, is), &
                                        have_disentangled, SM_q(:, :, is))
          ! QZYZ18 Eq.(51)
          call get_gauge_overlap_matrix(num_bands, num_wann, eigval, v_matrix, dis_manifold, &
                                        ik, num_states(ik), kmesh_info%nnlist(ik, nn), &
                                        num_states(kmesh_info%nnlist(ik, nn)), SHM_o(:, :, is), &
                                        have_disentangled, SHM_q(:, :, is))

          ! Assuming all neighbors of a given point are read in sequence!
          !
          do idir = 1, 3
            ! QZYZ18 Eq.(50)
            SR_q(:, :, ik, is, idir) = SR_q(:, :, ik, is, idir) &
                                       + kmesh_info%wb(nn)*kmesh_info%bk(idir, nn, ik) &
                                       *(SM_q(:, :, is) - SS_q(:, :, is))
            ! QZYZ18 Eq.(51)
            SHR_q(:, :, ik, is, idir) = SHR_q(:, :, ik, is, idir) &
                                        + kmesh_info%wb(nn)*kmesh_info%bk(idir, nn, ik) &
                                        *(SHM_q(:, :, is) - SH_q(:, :, ik, is))
          end do
        end do

        ik_prev = ik
      end do !ncount

      close (mmn_in)

      do is = 1, 3
        ! QZYZ18 Eq.(46)
        call fourier_q_to_R(num_kpts, wigner_seitz%nrpts, wigner_seitz%irvec, kpt_latt, &
                            SH_q(:, :, :, is), SH_R_temp(:, :, :, is))
        do idir = 1, 3
          ! QZYZ18 Eq.(44)
          call fourier_q_to_R(num_kpts, wigner_seitz%nrpts, wigner_seitz%irvec, kpt_latt, &
                              SR_q(:, :, :, is, idir), SR_R_temp(:, :, :, is, idir))
          ! QZYZ18 Eq.(45)
          call fourier_q_to_R(num_kpts, wigner_seitz%nrpts, wigner_seitz%irvec, kpt_latt, &
                              SHR_q(:, :, :, is, idir), SHR_R_temp(:, :, :, is, idir))
        end do
      end do

      do is = 1, 3
        ! QZYZ18 Eq.(46)
        call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, &
                                   SH_R_temp(:, :, :, is), SH_R(:, :, :, is))
        do idir = 1, 3
          ! QZYZ18 Eq.(44)
          call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, &
                                     SR_R_temp(:, :, :, is, idir), SR_R(:, :, :, is, idir))
          ! QZYZ18 Eq.(45)
          call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, &
                                     SHR_R_temp(:, :, :, is, idir), SHR_R(:, :, :, is, idir))
        end do
      end do

      SR_R = cmplx_i*SR_R
      SHR_R = cmplx_i*SHR_R

    end if !on_root

    call comms_bcast(SH_R(1, 1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts_pw90*3, error, comm)
    if (allocated(error)) return
    call comms_bcast(SR_R(1, 1, 1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts_pw90*3*3, error, comm)
    if (allocated(error)) return
    call comms_bcast(SHR_R(1, 1, 1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts_pw90*3*3, error, comm)
    if (allocated(error)) return

    ! end copying from get_AA_R, Junfeng Qiao

    if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
      call io_stopwatch_stop('get_oper: get_SHC_R', timer)
    return

    deallocate (SH_R_temp)
    deallocate (SR_R_temp)
    deallocate (SHR_R_temp)

101 call set_error_file(error, 'Error: Problem opening input file '//trim(seedname)//'.mmn', comm)
    return
102 call set_error_file(error, 'Error: Problem reading input file '//trim(seedname)//'.mmn', comm)
    return
109 call set_error_file(error, 'Error: Problem opening input file '//trim(seedname)//'.spn', comm)
    return
110 call set_error_file(error, 'Error: Problem reading input file '//trim(seedname)//'.spn', comm)
    return

  end subroutine get_SHC_R

get_SHC_R Subroutine

Contents

Source Code

Uses

Arguments

Calls

Called by

Source Code