tran_lcr Subroutine

private subroutine tran_lcr()

Uses

  • proc~~tran_lcr~~UsesGraph proc~tran_lcr tran_lcr module~w90_constants w90_constants proc~tran_lcr->module~w90_constants module~w90_io w90_io proc~tran_lcr->module~w90_io module~w90_parameters w90_parameters proc~tran_lcr->module~w90_parameters module~w90_io->module~w90_constants module~w90_parameters->module~w90_constants module~w90_parameters->module~w90_io

Arguments

None

Calls

proc~~tran_lcr~~CallsGraph proc~tran_lcr tran_lcr proc~io_date io_date proc~tran_lcr->proc~io_date proc~io_error io_error proc~tran_lcr->proc~io_error proc~tran_transfer tran_transfer proc~tran_lcr->proc~tran_transfer zgemm zgemm proc~tran_lcr->zgemm proc~tran_read_htx tran_read_htX proc~tran_lcr->proc~tran_read_htx proc~tran_read_htc tran_read_htC proc~tran_lcr->proc~tran_read_htc proc~io_file_unit io_file_unit proc~tran_lcr->proc~io_file_unit proc~tran_green tran_green proc~tran_lcr->proc~tran_green zgbsv zgbsv proc~tran_lcr->zgbsv proc~tran_read_htxy tran_read_htXY proc~tran_lcr->proc~tran_read_htxy proc~tran_transfer->proc~io_error zgesv zgesv proc~tran_transfer->zgesv zaxpy zaxpy proc~tran_transfer->zaxpy zcopy zcopy proc~tran_transfer->zcopy proc~tran_read_htx->proc~io_file_unit proc~tran_read_htc->proc~io_file_unit proc~tran_green->proc~io_error proc~tran_green->zgemm proc~tran_green->zgesv proc~tran_read_htxy->proc~io_file_unit

Called by

proc~~tran_lcr~~CalledByGraph proc~tran_lcr tran_lcr proc~tran_main tran_main proc~tran_main->proc~tran_lcr program~wannier wannier program~wannier->proc~tran_main proc~wannier_run wannier_run proc~wannier_run->proc~tran_main

Contents

Source Code


Source Code

  subroutine tran_lcr()
    !==================================================================!

    use w90_constants, only: dp, cmplx_0, cmplx_1, cmplx_i, pi
    use w90_io, only: io_error, io_stopwatch, seedname, io_date, &
      stdout, io_file_unit
    use w90_parameters, only: tran_num_ll, tran_num_rr, tran_num_cc, tran_num_lc, &
      tran_num_cr, tran_num_bandc, &
      tran_win_min, tran_win_max, tran_energy_step, &
      tran_use_same_lead, timing_level, tran_read_ht

    implicit none

    integer :: qc_unit, dos_unit
    integer :: ierr
    integer :: KL, KU, KC
    integer :: n_e, n, i, j, k, info
    integer, allocatable :: ipiv(:)
    real(kind=dp) ::  qc, dos
    real(kind=dp) ::  e_scan
    real(kind=dp), allocatable, dimension(:, :) :: hCband
    complex(kind=dp) :: e_scan_cmp
    complex(kind=dp), allocatable, dimension(:, :) :: hLC_cmp, hCR_cmp, &
                                                      totL, tottL, totR, tottR, &
                                                      g_surf_L, g_surf_R, g_C, g_C_inv, &
                                                      gR, gL, sLr, sRr, s1, s2, c1, c2
    character(len=50) :: filename
    character(len=9)  :: cdate, ctime

    if (timing_level > 1) call io_stopwatch('tran: lcr', 1)

    call io_date(cdate, ctime)

    qc_unit = io_file_unit()
    open (qc_unit, file=trim(seedname)//'_qc.dat', status='unknown', &
          form='formatted', action='write')
    write (qc_unit, *) '## written on '//cdate//' at '//ctime ! Date and time

    dos_unit = io_file_unit()
    open (dos_unit, file=trim(seedname)//'_dos.dat', status='unknown', &
          form='formatted', action='write')
    write (dos_unit, *) '## written on '//cdate//' at '//ctime ! Date and time

    KL = max(tran_num_lc, tran_num_cr, tran_num_bandc) - 1
    KU = KL
    KC = max(tran_num_lc, tran_num_cr)

    allocate (hCband(2*KL + KU + 1, tran_num_cc), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating hCband in tran_lcr')
    allocate (hLC_cmp(tran_num_ll, tran_num_lc), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating hLC_cmp in tran_lcr')
    allocate (hCR_cmp(tran_num_cr, tran_num_rr), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating hCR_cmp in tran_lcr')

    !If construct used only when reading matrices from file
    if (tran_read_ht) then
      allocate (hL0(tran_num_ll, tran_num_ll), stat=ierr)
      if (ierr /= 0) call io_error('Error in allocating hL0 in tran_lcr')
      allocate (hL1(tran_num_ll, tran_num_ll), stat=ierr)
      if (ierr /= 0) call io_error('Error in allocating hL1 in tran_lcr')
      allocate (hC(tran_num_cc, tran_num_cc), stat=ierr)
      if (ierr /= 0) call io_error('Error in allocating hC in tran_lcr')
      allocate (hLC(tran_num_ll, tran_num_lc), stat=ierr)
      if (ierr /= 0) call io_error('Error in allocating hLC in tran_lcr')
      allocate (hCR(tran_num_cr, tran_num_rr), stat=ierr)
      if (ierr /= 0) call io_error('Error in allocating hCR in tran_lcr')

      filename = trim(seedname)//'_htL.dat'
      call tran_read_htX(tran_num_ll, hL0, hL1, filename)

      if (.not. tran_use_same_lead) then
        allocate (hR0(tran_num_rr, tran_num_rr), stat=ierr)
        if (ierr /= 0) call io_error('Error in allocating hR0 in tran_lcr')
        allocate (hR1(tran_num_rr, tran_num_rr), stat=ierr)
        if (ierr /= 0) call io_error('Error in allocating hR1 in tran_lcr')
        filename = trim(seedname)//'_htR.dat'
        call tran_read_htX(tran_num_rr, hR0, hR1, filename)
      end if

      filename = trim(seedname)//'_htC.dat'
      call tran_read_htC(tran_num_cc, hC, filename)
      filename = trim(seedname)//'_htLC.dat'
      call tran_read_htXY(tran_num_ll, tran_num_lc, hLC, filename)
      filename = trim(seedname)//'_htCR.dat'
      call tran_read_htXY(tran_num_cr, tran_num_rr, hCR, filename)
    endif

    !  Banded matrix H_C  :  save memory !
    do j = 1, tran_num_cc
      do i = max(1, j - KU), min(tran_num_cc, j + KL)
        hCband(KL + KU + 1 + i - j, j) = hC(i, j)
      end do
    end do
    deallocate (hC, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating hC in tran_lcr')

    !  H_LC : to a complex matrix
    hLC_cmp(:, :) = cmplx(hLC(:, :), kind=dp)
    deallocate (hLC, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating hLC in tran_lcr')

    !  H_CR : to a complex matrix
    hCR_cmp(:, :) = cmplx(hCR(:, :), kind=dp)
    deallocate (hCR, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating hCR in tran_lcr')

    allocate (totL(tran_num_ll, tran_num_ll), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating totL in tran_lcr')
    allocate (tottL(tran_num_ll, tran_num_ll), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating tottL in tran_lcr')
    if (.not. tran_use_same_lead) then
      allocate (totR(tran_num_rr, tran_num_rr), stat=ierr)
      if (ierr /= 0) call io_error('Error in allocating totR in tran_lcr')
      allocate (tottR(tran_num_rr, tran_num_rr), stat=ierr)
      if (ierr /= 0) call io_error('Error in allocating tottR in tran_lcr')
    end if
    allocate (g_surf_L(tran_num_ll, tran_num_ll), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating g_surf_L in tran_lcr')
    allocate (g_surf_R(tran_num_rr, tran_num_rr), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating g_surf_R in tran_lcr')
    allocate (g_C_inv(2*KL + KU + 1, tran_num_cc), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating g_C_inv in tran_lcr')
    allocate (g_C(tran_num_cc, tran_num_cc), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating g_C in tran_lcr')
    allocate (sLr(tran_num_lc, tran_num_lc), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating sLr in tran_lcr')
    allocate (sRr(tran_num_cr, tran_num_cr), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating sRr in tran_lcr')
    allocate (gL(tran_num_lc, tran_num_lc), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating gL in tran_lcr')
    allocate (gR(tran_num_cr, tran_num_cr), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating gR in tran_lcr')
    allocate (c1(tran_num_lc, tran_num_ll), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating c1 in tran_lcr')
    allocate (c2(tran_num_cr, tran_num_rr), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating c2 in tran_lcr')
    allocate (s1(KC, KC), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating s1 in tran_lcr')
    allocate (s2(KC, KC), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating s2 in tran_lcr')
    allocate (ipiv(tran_num_cc), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating ipiv in tran_lcr')

    !  Loop over the energies
    n_e = floor((tran_win_max - tran_win_min)/tran_energy_step) + 1

    write (stdout, '(/1x,a)', advance='no') 'Calculating quantum conductance and density of states...'

    do n = 1, n_e

      e_scan = tran_win_min + real(n - 1, dp)*tran_energy_step

      !    compute conductance according to Fisher and Lee
      !    compute self-energies following Datta

      e_scan_cmp = e_scan + eta

      ! Surface green function for the left lead : g_surf_L
      call tran_transfer(totL, tottL, hL0, hL1, e_scan_cmp, tran_num_ll)
      call tran_green(totL, tottL, hL0, hL1, e_scan, g_surf_L, -1, 1, tran_num_ll)

      ! Self-energy (Sigma_L) : sLr = (hLC_cmp)^+ * g_surf_L * hLC_cmp
      c1 = cmplx_0
      sLr = cmplx_0
      call ZGEMM('C', 'N', tran_num_lc, tran_num_ll, tran_num_ll, cmplx_1, &
                 hLC_cmp, tran_num_ll, g_surf_L, tran_num_ll, cmplx_0, c1, tran_num_lc)
      call ZGEMM('N', 'N', tran_num_lc, tran_num_lc, tran_num_ll, cmplx_1, &
                 c1, tran_num_lc, hLC_cmp, tran_num_ll, cmplx_0, sLr, tran_num_lc)

      ! Surface green function for the right lead : g_surf_R
      if (tran_use_same_lead) then
        call tran_green(totL, tottL, hL0, hL1, e_scan, g_surf_R, 1, 1, tran_num_rr)
      else
        call tran_transfer(totR, tottR, hR0, hR1, e_scan_cmp, tran_num_rr)
        call tran_green(totR, tottR, hR0, hR1, e_scan, g_surf_R, 1, 1, tran_num_rr)
      end if

      ! Self-energy (Sigma_R) : sRr = hCR_cmp * g_surf_R * (hCR_cmp)^+
      c2 = cmplx_0
      sRr = cmplx_0
      call ZGEMM('N', 'N', tran_num_cr, tran_num_rr, tran_num_rr, cmplx_1, &
                 hCR_cmp, tran_num_cr, g_surf_R, tran_num_rr, cmplx_0, c2, tran_num_cr)
      call ZGEMM('N', 'C', tran_num_cr, tran_num_cr, tran_num_rr, cmplx_1, &
                 c2, tran_num_cr, hCR_cmp, tran_num_cr, cmplx_0, sRr, tran_num_cr)

      ! g_C^-1 = -H
      g_C_inv(:, :) = cmplx(-hCband(:, :), kind=dp)

      ! g_C^-1 = -H - Sigma_L^r
      do j = 1, tran_num_lc
        do i = max(1, j - KU), min(tran_num_lc, j + KL)
          g_C_inv(KL + KU + 1 + i - j, j) = g_C_inv(KL + KU + 1 + i - j, j) - sLr(i, j)
        end do
      end do

      ! g_C^-1 = -H - Sigma_L^r - Sigma_R^r
      do j = (tran_num_cc - tran_num_cr) + 1, tran_num_cc
        do i = max((tran_num_cc - tran_num_cr) + 1, j - (tran_num_cr - 1)), min(tran_num_cc, j + (tran_num_cr - 1))
          g_C_inv(KL + KU + 1 + i - j, j) = &
            g_C_inv(KL + KU + 1 + i - j, j) - &
            sRr(i - (tran_num_cc - tran_num_cr), j - (tran_num_cc - tran_num_cr))
        end do
      end do

      ! g_C^-1 = eI - H - Sigma_L^r - Sigma_R^r
      do i = 1, tran_num_cc
        g_C_inv(KL + KU + 1, i) = e_scan + g_C_inv(KL + KU + 1, i)
      end do

      ! invert g_C^-1 => g_C
      g_C = cmplx_0
      do i = 1, tran_num_cc
        g_C(i, i) = cmplx_1
      end do

      call ZGBSV(tran_num_cc, KL, KU, tran_num_cc, g_C_inv, 2*KL + KU + 1, ipiv, g_C, tran_num_cc, info)
      if (info .ne. 0) then
        write (stdout, *) 'ERROR: IN ZGBSV IN tran_lcr, INFO=', info
        call io_error('tran_lcr: problem in ZGBSV')
      end if

      ! Gamma_L = i(Sigma_L^r-Sigma_L^a)
      gL = cmplx_i*(sLr - conjg(transpose(sLr)))

      ! s1 = Gamma_L * g_C^r
      s1 = cmplx_0
      do j = 1, KC
        do i = 1, tran_num_lc
          do k = 1, tran_num_lc
            s1(i, j) = s1(i, j) + gL(i, k)*g_C(k, j + (tran_num_cc - KC))
          end do
        end do
      end do

      ! Gamma_R = i(Sigma_R^r-Sigma_R^a)
      gR = cmplx_i*(sRr - conjg(transpose(sRr)))

      ! s2 = Gamma_R * g_C^a
      s2 = cmplx_0
      do j = 1, KC
        do i = 1, tran_num_cr
          do k = 1, tran_num_cr
            s2(i + (KC - tran_num_cr), j) = s2(i + (KC - tran_num_cr), j) &
                                            + gR(i, k)*conjg(g_C(j, k + (tran_num_cc - tran_num_cr)))
          end do
        end do
      end do

      qc = 0.0_dp
      do i = 1, KC
        do j = 1, KC
          qc = qc + real(s1(i, j)*s2(j, i), dp)
        end do
      end do
      write (qc_unit, '(f15.9,f18.9)') e_scan, qc

      ! compute density of states for the conductor layer

      dos = 0.0_dp
      do i = 1, tran_num_cc
        dos = dos - aimag(g_C(i, i))
      end do
      dos = dos/pi
      write (dos_unit, '(f15.9,f18.9)') e_scan, dos

    end do

    write (stdout, '(a)') ' done'

    close (qc_unit)
    close (dos_unit)

    deallocate (ipiv, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating ipiv in tran_lcr')
    deallocate (s2, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating s2 in tran_lcr')
    deallocate (s1, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating s1 in tran_lcr')
    deallocate (c2, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating c2 in tran_lcr')
    deallocate (c1, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating c1 in tran_lcr')
    deallocate (gR, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating gR in tran_lcr')
    deallocate (gL, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating gL in tran_lcr')
    deallocate (sRr, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating sRr in tran_lcr')
    deallocate (sLr, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating sLr in tran_lcr')
    deallocate (g_C, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating g_C in tran_lcr')
    deallocate (g_C_inv, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating g_C_inv in tran_lcr')
    deallocate (g_surf_R, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating g_surf_R in tran_lcr')
    deallocate (g_surf_L, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating g_surf_L in tran_lcr')
    if (allocated(tottR)) deallocate (tottR, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating tottR in tran_lcr')
    if (allocated(totR)) deallocate (totR, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating totR in tran_lcr')
    deallocate (tottL, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating tottL in tran_lcr')
    deallocate (totL, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating totL in tran_lcr')
    deallocate (hCR_cmp, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating hCR_cmp in tran_lcr')
    deallocate (hLC_cmp, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating hLC_cmp in tran_lcr')
    deallocate (hCband, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating hCband in tran_lcr')

    if (timing_level > 1) call io_stopwatch('tran: lcr', 2)

    return

  end subroutine tran_lcr