subroutine tran_green(tot, tott, h_00, h_01, e_scan, g, igreen, invert, nxx)
!==================================================================!
! construct green's functions
!
! igreen = -1 left surface
! igreen = 1 right surface
! igreen = 0 bulk
! invert = 0 computes g^-1
! invert = 1 computes g^-1 and g
!==================================================================!
use w90_constants, only: dp, cmplx_0, cmplx_1
use w90_io, only: stdout, io_error
implicit none
integer, intent(in) :: nxx
integer, intent(in) :: igreen
integer, intent(in) :: invert
real(kind=dp), intent(in) :: e_scan
real(kind=dp), intent(in) :: h_00(nxx, nxx), h_01(nxx, nxx)
complex(kind=dp), intent(in) :: tot(nxx, nxx), tott(nxx, nxx)
complex(kind=dp), intent(out) :: g(nxx, nxx)
integer :: ierr, info
integer :: i
integer, allocatable :: ipiv(:)
complex(kind=dp), allocatable, dimension(:, :) :: g_inv, eh_00
complex(kind=dp), allocatable, dimension(:, :) :: s1, s2, c1
allocate (ipiv(nxx), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating ipiv in tran_green')
allocate (g_inv(nxx, nxx))
if (ierr /= 0) call io_error('Error in allocating g_inv in tran_green')
allocate (eh_00(nxx, nxx))
if (ierr /= 0) call io_error('Error in allocating eh_00 in tran_green')
allocate (c1(nxx, nxx))
if (ierr /= 0) call io_error('Error in allocating c1 in tran_green')
allocate (s1(nxx, nxx))
if (ierr /= 0) call io_error('Error in allocating s1 in tran_green')
allocate (s2(nxx, nxx))
if (ierr /= 0) call io_error('Error in allocating s2 in tran_green')
c1(:, :) = cmplx(h_01(:, :), kind=dp)
select case (igreen)
case (1)
! construct the surface green's function g00
s1 = cmplx_0
! s1 = H_01 * T
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, c1, nxx, tot, nxx, cmplx_0, s1, nxx)
! eh_00 = -H_00 - H_01*T
eh_00(:, :) = cmplx(-h_00(:, :), kind=dp) - s1(:, :)
! eh_00 = e_scan -H_00 - H_01*T
do i = 1, nxx
eh_00(i, i) = cmplx(e_scan, kind=dp) + eh_00(i, i)
end do
g_inv(:, :) = eh_00(:, :)
! identity
g = cmplx_0
do i = 1, nxx
g(i, i) = cmplx_1
end do
if (invert .eq. 1) then
call ZGESV(nxx, nxx, eh_00, nxx, ipiv, g, nxx, info)
if (info .ne. 0) then
write (stdout, *) 'ERROR: IN ZGESV IN tran_green, INFO=', info
call io_error('tran_green: problem in ZGESV 1')
end if
end if
case (-1)
! construct the dual surface green's function gbar00
s1 = cmplx_0
! s1 = H_01^+ * T^bar
call ZGEMM('C', 'N', nxx, nxx, nxx, cmplx_1, c1, nxx, tott, nxx, cmplx_0, s1, nxx)
! s1 = -H_00 - H_01^+ * T^bar
eh_00(:, :) = cmplx(-h_00(:, :), kind=dp) - s1(:, :)
! s1 = e_scan - H_00 - H_01^+ * T^bar
do i = 1, nxx
eh_00(i, i) = cmplx(e_scan, kind=dp) + eh_00(i, i)
end do
g_inv(:, :) = eh_00(:, :)
! identity
g = cmplx_0
do i = 1, nxx
g(i, i) = cmplx_1
end do
if (invert .eq. 1) then
call ZGESV(nxx, nxx, eh_00, nxx, ipiv, g, nxx, info)
if (info .ne. 0) then
write (stdout, *) 'ERROR: IN ZGESV IN tran_green, INFO=', info
call io_error('tran_green: problem in ZGESV 2')
end if
end if
case (0)
! construct the bulk green's function gnn or (if surface=.true.) the
! sub-surface green's function
s1 = cmplx_0
s2 = cmplx_0
! s1 = H_01 * T
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, c1, nxx, tot, nxx, cmplx_0, s1, nxx)
! s2 = H_01^+ * T^bar
call ZGEMM('C', 'N', nxx, nxx, nxx, cmplx_1, c1, nxx, tott, nxx, cmplx_0, s2, nxx)
eh_00(:, :) = cmplx(-h_00(:, :), kind=dp) - s1(:, :) - s2(:, :)
do i = 1, nxx
eh_00(i, i) = cmplx(e_scan, kind=dp) + eh_00(i, i)
end do
g_inv(:, :) = eh_00(:, :)
! identity
g = cmplx_0
do i = 1, nxx
g(i, i) = cmplx_1
end do
if (invert .eq. 1) then
call ZGESV(nxx, nxx, eh_00, nxx, ipiv, g, nxx, info)
if (info .ne. 0) then
write (stdout, *) 'ERROR: IN ZGESV IN tran_green, INFO=', info
call io_error('tran_green: problem in ZGESV 3')
end if
end if
end select
deallocate (s2)
if (ierr /= 0) call io_error('Error in deallocating s2 in tran_green')
deallocate (s1)
if (ierr /= 0) call io_error('Error in deallocating s1 in tran_green')
deallocate (c1)
if (ierr /= 0) call io_error('Error in deallocating c1 in tran_green')
deallocate (eh_00)
if (ierr /= 0) call io_error('Error in deallocating eh_00 in tran_green')
deallocate (g_inv)
if (ierr /= 0) call io_error('Error in deallocating g_inv in tran_green')
deallocate (ipiv)
if (ierr /= 0) call io_error('Error in deallocating ipiv in tran_green')
return
end subroutine tran_green