get_AA_R Subroutine

public subroutine get_AA_R()

Uses

  • proc~~get_aa_r~~UsesGraph proc~get_aa_r get_AA_R module~w90_postw90_common w90_postw90_common proc~get_aa_r->module~w90_postw90_common module~w90_constants w90_constants proc~get_aa_r->module~w90_constants module~w90_comms w90_comms proc~get_aa_r->module~w90_comms module~w90_parameters w90_parameters proc~get_aa_r->module~w90_parameters module~w90_io w90_io proc~get_aa_r->module~w90_io module~w90_postw90_common->module~w90_constants module~w90_postw90_common->module~w90_comms module~w90_comms->module~w90_constants module~w90_comms->module~w90_io module~w90_parameters->module~w90_constants module~w90_parameters->module~w90_io module~w90_io->module~w90_constants

AA_a(R) = <0|r_a|R> is the Fourier transform of the Berrry connection AA_a(k) = i (a=x,y,z)

Arguments

None

Calls

proc~~get_aa_r~~CallsGraph proc~get_aa_r get_AA_R proc~io_error io_error proc~get_aa_r->proc~io_error proc~io_file_unit io_file_unit proc~get_aa_r->proc~io_file_unit

Called by

proc~~get_aa_r~~CalledByGraph proc~get_aa_r get_AA_R proc~k_slice k_slice proc~k_slice->proc~get_aa_r proc~k_path k_path proc~k_path->proc~get_aa_r proc~gyrotropic_main gyrotropic_main proc~gyrotropic_main->proc~get_aa_r proc~berry_main berry_main proc~berry_main->proc~get_aa_r proc~berry_get_sc_klist berry_get_sc_klist proc~berry_main->proc~berry_get_sc_klist proc~wham_get_eig_uu_hh_aa_sc_tb_conv wham_get_eig_UU_HH_AA_sc_TB_conv proc~wham_get_eig_uu_hh_aa_sc_tb_conv->proc~get_aa_r proc~berry_get_sc_klist->proc~wham_get_eig_uu_hh_aa_sc_tb_conv

Contents

Source Code


Source Code

  subroutine get_AA_R
    !==================================================
    !
    !! AA_a(R) = <0|r_a|R> is the Fourier transform
    !! of the Berrry connection AA_a(k) = i<u|del_a u>
    !! (a=x,y,z)
    !
    !==================================================

    use w90_constants, only: dp, cmplx_0, cmplx_i
    use w90_parameters, only: num_kpts, nntot, num_wann, wb, bk, timing_level, &
      num_bands, ndimwin, nnlist, have_disentangled, &
      transl_inv, nncell, effective_model
    use w90_postw90_common, only: nrpts
    use w90_io, only: stdout, io_file_unit, io_error, io_stopwatch, &
      seedname
    use w90_comms, only: on_root, comms_bcast

    complex(kind=dp), allocatable :: AA_q(:, :, :, :)
    complex(kind=dp), allocatable :: AA_q_diag(:, :)
    complex(kind=dp), allocatable :: S_o(:, :)
    complex(kind=dp), allocatable :: S(:, :)
    integer                       :: n, m, i, j, &
                                     ik, ik2, ik_prev, nn, inn, nnl, nnm, nnn, &
                                     idir, ncount, nn_count, mmn_in, &
                                     nb_tmp, nkp_tmp, nntot_tmp, file_unit, &
                                     ir, io, ivdum(3), ivdum_old(3)
    integer, allocatable          :: num_states(:)
    real(kind=dp)                 :: m_real, m_imag, rdum1_real, rdum1_imag, &
                                     rdum2_real, rdum2_imag, rdum3_real, rdum3_imag
    logical                       :: nn_found
    character(len=60)             :: header

    if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_AA_R', 1)

    if (.not. allocated(AA_R)) then
      allocate (AA_R(num_wann, num_wann, nrpts, 3))
    else
      if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_AA_R', 2)
      return
    end if

    ! Real-space position matrix elements read from file
    !
    if (effective_model) then
      if (.not. allocated(HH_R)) call io_error( &
        'Error in get_AA_R: Must read file'//trim(seedname)//'_HH_R.dat first')
      AA_R = cmplx_0
      if (on_root) then
        write (stdout, '(/a)') ' Reading position matrix elements from file ' &
          //trim(seedname)//'_AA_R.dat'
        file_unit = io_file_unit()
        open (file_unit, file=trim(seedname)//'_AA_R.dat', form='formatted', &
              status='old', err=103)
        read (file_unit, *) ! header
        ir = 1
        ivdum_old(:) = 0
        n = 1
        do
          read (file_unit, '(5I5,6F12.6)', iostat=io) &
            ivdum(1:3), j, i, rdum1_real, rdum1_imag, &
            rdum2_real, rdum2_imag, rdum3_real, rdum3_imag
          if (io < 0) exit
          if (i < 1 .or. i > num_wann .or. j < 1 .or. j > num_wann) then
            write (stdout, *) 'num_wann=', num_wann, '  i=', i, '  j=', j
            call io_error('Error in get_AA_R: orbital indices out of bounds')
          endif
          if (n > 1) then
            if (ivdum(1) /= ivdum_old(1) .or. ivdum(2) /= ivdum_old(2) .or. &
                ivdum(3) /= ivdum_old(3)) ir = ir + 1
          endif
          ivdum_old = ivdum
          AA_R(j, i, ir, 1) = AA_R(j, i, ir, 1) + cmplx(rdum1_real, rdum1_imag, kind=dp)
          AA_R(j, i, ir, 2) = AA_R(j, i, ir, 2) + cmplx(rdum2_real, rdum2_imag, kind=dp)
          AA_R(j, i, ir, 3) = AA_R(j, i, ir, 3) + cmplx(rdum3_real, rdum3_imag, kind=dp)
          n = n + 1
        enddo
        close (file_unit)
        ! AA_R may not contain the same number of R-vectors as HH_R
        ! (e.g., if a diagonal representation of the position matrix
        ! elements is used, but it cannot be larger
        if (ir > nrpts) then
          write (stdout, *) 'ir=', ir, '  nrpts=', nrpts
          call io_error('Error in get_AA_R: inconsistent nrpts values')
        endif
      endif
      call comms_bcast(AA_R(1, 1, 1, 1), num_wann*num_wann*nrpts*3)
      if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_AA_R', 2)
      return
    endif

    ! Real-space position matrix elements calculated by Fourier
    ! transforming overlap matrices defined on the ab-initio
    ! reciprocal mesh
    !
    ! Do everything on root, broadcast AA_R at the end (smaller than S_o)
    !
    if (on_root) then

      allocate (AA_q(num_wann, num_wann, num_kpts, 3))
      allocate (AA_q_diag(num_wann, 3))
      allocate (S_o(num_bands, num_bands))
      allocate (S(num_wann, num_wann))

      allocate (num_states(num_kpts))
      do ik = 1, num_kpts
        if (have_disentangled) then
          num_states(ik) = ndimwin(ik)
        else
          num_states(ik) = num_wann
        endif
      enddo

      mmn_in = io_file_unit()
      open (unit=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_AA_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) &
        call io_error(trim(seedname)//'.mmn has wrong number of bands')
      if (nkp_tmp .ne. num_kpts) &
        call io_error(trim(seedname)//'.mmn has wrong number of k-points')
      if (nntot_tmp .ne. nntot) &
        call io_error &
        (trim(seedname)//'.mmn has wrong number of nearest neighbours')

      AA_q = cmplx_0
      ik_prev = 0

      ! Composite loop over k-points ik (outer loop) and neighbors ik2 (inner)
      do ncount = 1, num_kpts*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)
          enddo
        enddo
        !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, nntot
          if ((ik2 .eq. nnlist(ik, inn)) .and. &
              (nnl .eq. nncell(1, ik, inn)) .and. &
              (nnm .eq. nncell(2, ik, inn)) .and. &
              (nnn .eq. nncell(3, ik, inn))) then
            if (.not. nn_found) then
              nn_found = .true.
              nn = inn
            else
              call io_error('Error reading '//trim(seedname)//'.mmn.&
                   & More than one matching nearest neighbour found')
            endif
          endif
        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 io_error('Neighbour not found')
        end if
        nn_count = nn_count + 1 !Check: can also be place after nn=inn (?)

        ! Wannier-gauge overlap matrix S in the projected subspace
        !
        call get_gauge_overlap_matrix( &
          ik, num_states(ik), &
          nnlist(ik, nn), num_states(nnlist(ik, nn)), &
          S_o, S)

        ! Berry connection matrix
        ! Assuming all neighbors of a given point are read in sequence!
        !
        if (transl_inv .and. ik .ne. ik_prev) AA_q_diag(:, :) = cmplx_0
        do idir = 1, 3
          AA_q(:, :, ik, idir) = AA_q(:, :, ik, idir) &
                                 + cmplx_i*wb(nn)*bk(idir, nn, ik)*S(:, :)
          if (transl_inv) then
            !
            ! Rewrite band-diagonal elements a la Eq.(31) of MV97
            !
            do i = 1, num_wann
              AA_q_diag(i, idir) = AA_q_diag(i, idir) &
                                   - wb(nn)*bk(idir, nn, ik)*aimag(log(S(i, i)))
            enddo
          endif
        end do
        ! Assuming all neighbors of a given point are read in sequence!
        if (nn_count == nntot) then !looped over all neighbors
          do idir = 1, 3
            if (transl_inv) then
              do n = 1, num_wann
                AA_q(n, n, ik, idir) = AA_q_diag(n, idir)
              enddo
            endif
            !
            ! Since Eq.(44) WYSV06 does not preserve the Hermiticity of the
            ! Berry potential matrix, take Hermitean part (whether this
            ! makes a difference or not for e.g. the AHC, depends on which
            ! expression is used to evaluate the Berry curvature.
            ! See comments in berry_wanint.F90)
            !
            AA_q(:, :, ik, idir) = &
              0.5_dp*(AA_q(:, :, ik, idir) &
                      + conjg(transpose(AA_q(:, :, ik, idir))))
          enddo
        end if

        ik_prev = ik
      enddo !ncount

      close (mmn_in)

      call fourier_q_to_R(AA_q(:, :, :, 1), AA_R(:, :, :, 1))
      call fourier_q_to_R(AA_q(:, :, :, 2), AA_R(:, :, :, 2))
      call fourier_q_to_R(AA_q(:, :, :, 3), AA_R(:, :, :, 3))

    endif !on_root

    call comms_bcast(AA_R(1, 1, 1, 1), num_wann*num_wann*nrpts*3)

    if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_AA_R', 2)
    return

101 call io_error &
      ('Error: Problem opening input file '//trim(seedname)//'.mmn')
102 call io_error &
      ('Error: Problem reading input file '//trim(seedname)//'.mmn')
103 call io_error('Error in get_AA_R: problem opening file '// &
                  trim(seedname)//'_AA_R.dat')

  end subroutine get_AA_R