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