subroutine tran_green(tot, tott, h_00, h_01, e_scan, g, igreen, invert, nxx, stdout, error, comm)
!================================================!
! 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_error, only: w90_error_type, set_error_alloc, set_error_fatal, set_error_dealloc
implicit none
type(w90_error_type), allocatable, intent(out) :: error
type(w90_comm_type), intent(in) :: comm
integer, intent(in) :: nxx
integer, intent(in) :: stdout
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) then
call set_error_alloc(error, 'Error in allocating ipiv in tran_green', comm)
return
end if
allocate (g_inv(nxx, nxx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating g_inv in tran_green', comm)
return
end if
allocate (eh_00(nxx, nxx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating eh_00 in tran_green', comm)
return
end if
allocate (c1(nxx, nxx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating c1 in tran_green', comm)
return
end if
allocate (s1(nxx, nxx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating s1 in tran_green', comm)
return
end if
allocate (s2(nxx, nxx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating s2 in tran_green', comm)
return
end if
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 set_error_fatal(error, 'tran_green: problem in ZGESV 1', comm)
return
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 set_error_fatal(error, 'tran_green: problem in ZGESV 2', comm)
return
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 set_error_fatal(error, 'tran_green: problem in ZGESV 3', comm)
return
end if
end if
end select
deallocate (s2, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating s2 in tran_green', comm)
return
end if
deallocate (s1, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating s1 in tran_green', comm)
return
end if
deallocate (c1, stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in deallocating c1 in tran_green', comm)
return
end if
deallocate (eh_00, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating eh_00 in tran_green', comm)
return
end if
deallocate (g_inv, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating g_inv in tran_green', comm)
return
end if
deallocate (ipiv, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating ipiv in tran_green', comm)
return
end if
return
end subroutine tran_green
