subroutine tran_lcr(transport, hC, hCR, hL0, hL1, hLC, hR0, hR1, timing_level, stdout, seedname, &
timer, error, comm)
!================================================!
use w90_constants, only: dp, cmplx_0, cmplx_1, cmplx_i, pi
use w90_io, only: io_stopwatch_start, io_stopwatch_stop, io_date
use w90_wannier90_types, only: transport_type
use w90_error, only: w90_error_type, set_error_alloc, set_error_dealloc, set_error_fatal
use w90_types, only: timer_list_type
implicit none
! arguments
integer, intent(in) :: stdout
integer, intent(in) :: timing_level
real(kind=dp), allocatable, intent(inout) :: hC(:, :)
real(kind=dp), allocatable, intent(inout) :: hCR(:, :)
real(kind=dp), allocatable, intent(inout) :: hL0(:, :)
real(kind=dp), allocatable, intent(inout) :: hL1(:, :)
real(kind=dp), allocatable, intent(inout) :: hLC(:, :)
real(kind=dp), allocatable, intent(inout) :: hR0(:, :)
real(kind=dp), allocatable, intent(inout) :: hR1(:, :)
type(transport_type), intent(in) :: transport
type(timer_list_type), intent(inout) :: timer
type(w90_error_type), allocatable, intent(out) :: error
type(w90_comm_type), intent(in) :: comm
character(len=50), intent(in) :: seedname
! local variables
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 :: hCband(:, :)
complex(kind=dp) :: e_scan_cmp
complex(kind=dp), allocatable :: hLC_cmp(:, :), hCR_cmp(:, :)
complex(kind=dp), allocatable :: totL(:, :), tottL(:, :), totR(:, :), tottR(:, :)
complex(kind=dp), allocatable :: g_surf_L(:, :), g_surf_R(:, :)
complex(kind=dp), allocatable :: g_C(:, :), g_C_inv(:, :)
complex(kind=dp), allocatable :: gR(:, :), gL(:, :)
complex(kind=dp), allocatable :: sLr(:, :), sRr(:, :)
complex(kind=dp), allocatable :: s1(:, :), s2(:, :)
complex(kind=dp), allocatable :: c1(:, :), c2(:, :)
character(len=50) :: filename
character(len=9) :: cdate, ctime
if (timing_level > 1) call io_stopwatch_start('tran: lcr', timer)
call io_date(cdate, ctime)
open (newunit=qc_unit, file=trim(seedname)//'_qc.dat', status='unknown', &
form='formatted', action='write')
write (qc_unit, *) '## written on '//cdate//' at '//ctime ! Date and time
open (newunit=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(transport%num_lc, transport%num_cr, transport%num_bandc) - 1
KU = KL
KC = max(transport%num_lc, transport%num_cr)
allocate (hCband(2*KL + KU + 1, transport%num_cc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hCband in tran_lcr', comm)
return
end if
allocate (hLC_cmp(transport%num_ll, transport%num_lc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hLC_cmp in tran_lcr', comm)
return
end if
allocate (hCR_cmp(transport%num_cr, transport%num_rr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hCR_cmp in tran_lcr', comm)
return
end if
!If construct used only when reading matrices from file
if (transport%read_ht) then
allocate (hL0(transport%num_ll, transport%num_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hL0 in tran_lcr', comm)
return
end if
allocate (hL1(transport%num_ll, transport%num_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hL1 in tran_lcr', comm)
return
end if
allocate (hC(transport%num_cc, transport%num_cc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hC in tran_lcr', comm)
return
end if
allocate (hLC(transport%num_ll, transport%num_lc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hLC in tran_lcr', comm)
return
end if
allocate (hCR(transport%num_cr, transport%num_rr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hCR in tran_lcr', comm)
return
end if
filename = trim(seedname)//'_htL.dat'
call tran_read_htX(transport%num_ll, hL0, hL1, filename, stdout, error, comm)
if (allocated(error)) return
if (.not. transport%use_same_lead) then
allocate (hR0(transport%num_rr, transport%num_rr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hR0 in tran_lcr', comm)
return
end if
allocate (hR1(transport%num_rr, transport%num_rr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hR1 in tran_lcr', comm)
return
end if
filename = trim(seedname)//'_htR.dat'
call tran_read_htX(transport%num_rr, hR0, hR1, filename, stdout, error, comm)
if (allocated(error)) return
end if
filename = trim(seedname)//'_htC.dat'
call tran_read_htC(transport%num_cc, hC, filename, stdout, error, comm)
if (allocated(error)) return
filename = trim(seedname)//'_htLC.dat'
call tran_read_htXY(transport%num_ll, transport%num_lc, hLC, filename, stdout, error, comm)
if (allocated(error)) return
filename = trim(seedname)//'_htCR.dat'
call tran_read_htXY(transport%num_cr, transport%num_rr, hCR, filename, stdout, error, comm)
if (allocated(error)) return
end if
! Banded matrix H_C : save memory !
do j = 1, transport%num_cc
do i = max(1, j - KU), min(transport%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) then
call set_error_dealloc(error, 'Error in deallocating hC in tran_lcr', comm)
return
end if
! H_LC : to a complex matrix
hLC_cmp(:, :) = cmplx(hLC(:, :), kind=dp)
deallocate (hLC, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating hLC in tran_lcr', comm)
return
end if
! H_CR : to a complex matrix
hCR_cmp(:, :) = cmplx(hCR(:, :), kind=dp)
deallocate (hCR, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating hCR in tran_lcr', comm)
return
end if
allocate (totL(transport%num_ll, transport%num_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating totL in tran_lcr', comm)
return
end if
allocate (tottL(transport%num_ll, transport%num_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating tottL in tran_lcr', comm)
return
end if
if (.not. transport%use_same_lead) then
allocate (totR(transport%num_rr, transport%num_rr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating totR in tran_lcr', comm)
return
end if
allocate (tottR(transport%num_rr, transport%num_rr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating tottR in tran_lcr', comm)
return
end if
end if
allocate (g_surf_L(transport%num_ll, transport%num_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating g_surf_L in tran_lcr', comm)
return
end if
allocate (g_surf_R(transport%num_rr, transport%num_rr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating g_surf_R in tran_lcr', comm)
return
end if
allocate (g_C_inv(2*KL + KU + 1, transport%num_cc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating g_C_inv in tran_lcr', comm)
return
end if
allocate (g_C(transport%num_cc, transport%num_cc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating g_C in tran_lcr', comm)
return
end if
allocate (sLr(transport%num_lc, transport%num_lc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating sLr in tran_lcr', comm)
return
end if
allocate (sRr(transport%num_cr, transport%num_cr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating sRr in tran_lcr', comm)
return
end if
allocate (gL(transport%num_lc, transport%num_lc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating gL in tran_lcr', comm)
return
end if
allocate (gR(transport%num_cr, transport%num_cr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating gR in tran_lcr', comm)
return
end if
allocate (c1(transport%num_lc, transport%num_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating c1 in tran_lcr', comm)
return
end if
allocate (c2(transport%num_cr, transport%num_rr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating c2 in tran_lcr', comm)
return
end if
allocate (s1(KC, KC), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating s1 in tran_lcr', comm)
return
end if
allocate (s2(KC, KC), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating s2 in tran_lcr', comm)
return
end if
allocate (ipiv(transport%num_cc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating ipiv in tran_lcr', comm)
return
end if
! Loop over the energies
n_e = floor((transport%win_max - transport%win_min)/transport%energy_step) + 1
write (stdout, '(/1x,a)', advance='no') 'Calculating quantum conductance and &
&density of states...'
do n = 1, n_e
e_scan = transport%win_min + real(n - 1, dp)*transport%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, transport%num_ll, stdout, error, comm)
if (allocated(error)) return
call tran_green(totL, tottL, hL0, hL1, e_scan, g_surf_L, -1, 1, transport%num_ll, stdout, &
error, comm)
if (allocated(error)) return
! Self-energy (Sigma_L) : sLr = (hLC_cmp)^+ * g_surf_L * hLC_cmp
c1 = cmplx_0
sLr = cmplx_0
call ZGEMM('C', 'N', transport%num_lc, transport%num_ll, transport%num_ll, cmplx_1, &
hLC_cmp, transport%num_ll, g_surf_L, transport%num_ll, cmplx_0, c1, transport%num_lc)
call ZGEMM('N', 'N', transport%num_lc, transport%num_lc, transport%num_ll, cmplx_1, &
c1, transport%num_lc, hLC_cmp, transport%num_ll, cmplx_0, sLr, transport%num_lc)
! Surface green function for the right lead : g_surf_R
if (transport%use_same_lead) then
call tran_green(totL, tottL, hL0, hL1, e_scan, g_surf_R, 1, 1, transport%num_rr, stdout, &
error, comm)
if (allocated(error)) return
else
call tran_transfer(totR, tottR, hR0, hR1, e_scan_cmp, transport%num_rr, stdout, error, comm)
if (allocated(error)) return
call tran_green(totR, tottR, hR0, hR1, e_scan, g_surf_R, 1, 1, transport%num_rr, stdout, &
error, comm)
if (allocated(error)) return
end if
! Self-energy (Sigma_R) : sRr = hCR_cmp * g_surf_R * (hCR_cmp)^+
c2 = cmplx_0
sRr = cmplx_0
call ZGEMM('N', 'N', transport%num_cr, transport%num_rr, transport%num_rr, cmplx_1, &
hCR_cmp, transport%num_cr, g_surf_R, transport%num_rr, cmplx_0, c2, transport%num_cr)
call ZGEMM('N', 'C', transport%num_cr, transport%num_cr, transport%num_rr, cmplx_1, &
c2, transport%num_cr, hCR_cmp, transport%num_cr, cmplx_0, sRr, transport%num_cr)
! g_C^-1 = -H
g_C_inv(:, :) = cmplx(-hCband(:, :), kind=dp)
! g_C^-1 = -H - Sigma_L^r
do j = 1, transport%num_lc
do i = max(1, j - KU), min(transport%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 = (transport%num_cc - transport%num_cr) + 1, transport%num_cc
do i = max((transport%num_cc - transport%num_cr) + 1, j - (transport%num_cr - 1)), &
min(transport%num_cc, j + (transport%num_cr - 1))
g_C_inv(KL + KU + 1 + i - j, j) = &
g_C_inv(KL + KU + 1 + i - j, j) - &
sRr(i - (transport%num_cc - transport%num_cr), j - (transport%num_cc - transport%num_cr))
end do
end do
! g_C^-1 = eI - H - Sigma_L^r - Sigma_R^r
do i = 1, transport%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, transport%num_cc
g_C(i, i) = cmplx_1
end do
call zgbsv(transport%num_cc, KL, KU, transport%num_cc, g_C_inv, 2*KL + KU + 1, ipiv, g_C, &
transport%num_cc, info)
if (info .ne. 0) then
write (stdout, *) 'ERROR: IN ZGBSV IN tran_lcr, INFO=', info
call set_error_fatal(error, 'tran_lcr: problem in ZGBSV', comm)
return
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, transport%num_lc
do k = 1, transport%num_lc
s1(i, j) = s1(i, j) + gL(i, k)*g_C(k, j + (transport%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, transport%num_cr
do k = 1, transport%num_cr
s2(i + (KC - transport%num_cr), j) = s2(i + (KC - transport%num_cr), j) &
+ gR(i, k)*conjg(g_C(j, k &
+ (transport%num_cc - transport%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, transport%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) then
call set_error_dealloc(error, 'Error in deallocating ipiv in tran_lcr', comm)
return
end if
deallocate (s2, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating s2 in tran_lcr', comm)
return
end if
deallocate (s1, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating s1 in tran_lcr', comm)
return
end if
deallocate (c2, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating c2 in tran_lcr', comm)
return
end if
deallocate (c1, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating c1 in tran_lcr', comm)
return
end if
deallocate (gR, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating gR in tran_lcr', comm)
return
end if
deallocate (gL, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating gL in tran_lcr', comm)
return
end if
deallocate (sRr, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating sRr in tran_lcr', comm)
return
end if
deallocate (sLr, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating sLr in tran_lcr', comm)
return
end if
deallocate (g_C, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating g_C in tran_lcr', comm)
return
end if
deallocate (g_C_inv, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating g_C_inv in tran_lcr', comm)
return
end if
deallocate (g_surf_R, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating g_surf_R in tran_lcr', comm)
return
end if
deallocate (g_surf_L, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating g_surf_L in tran_lcr', comm)
return
end if
if (allocated(tottR)) then
deallocate (tottR, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating tottR in tran_lcr', comm)
return
end if
end if
if (allocated(totR)) then
deallocate (totR, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating totR in tran_lcr', comm)
return
end if
end if
deallocate (tottL, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating tottL in tran_lcr', comm)
return
end if
deallocate (totL, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating totL in tran_lcr', comm)
return
end if
deallocate (hCR_cmp, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating hCR_cmp in tran_lcr', comm)
return
end if
deallocate (hLC_cmp, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating hLC_cmp in tran_lcr', comm)
return
end if
deallocate (hCband, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating hCband in tran_lcr', comm)
return
end if
if (timing_level > 1) call io_stopwatch_stop('tran: lcr', timer)
return
end subroutine tran_lcr
