subroutine tran_lcr_2c2_build_ham(pl_warning, real_space_ham, fermi_energy_list, kpt_latt, &
num_wann, transport, print_output, real_lattice, mp_grid, &
ham_r, irvec, nrpts, wannier_centres_translated, one_dim_vec, &
nrpts_one_dim, num_pl, coord, tran_sorted_idx, hC, hCR, hL0, &
hL1, hLC, hR0, hR1, hr_one_dim, irvec_max, stdout, seedname, &
timer, error, comm)
!================================================!
! Builds hamiltonians blocks required for the
! Greens function caclulations of the quantum
! conductance according to the 2c2 geometry.
! Leads are also symmetrised, in that unit cell
! sub-blocks are copied to create truely ideal
! leads.
!================================================!
use w90_constants, only: dp, eps5
use w90_io, only: io_date, io_stopwatch_start, io_stopwatch_stop
use w90_types, only: print_output_type, timer_list_type
use w90_wannier90_types, only: transport_type, real_space_ham_type
use w90_error, only: w90_error_type, set_error_fatal, set_error_alloc, set_error_dealloc
implicit none
! arguments
integer, intent(in) :: coord(3)
integer, intent(in) :: irvec(:, :)
integer, intent(inout) :: irvec_max
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: num_wann
integer, intent(inout) :: nrpts
integer, intent(inout) :: nrpts_one_dim
integer, intent(inout) :: num_pl
integer, intent(inout) :: one_dim_vec
integer, intent(in) :: stdout
integer, intent(in) :: tran_sorted_idx(:)
real(kind=dp), allocatable, intent(inout) :: hr_one_dim(:, :, :) ! de/realloc'd below
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(:, :)
real(kind=dp), intent(in) :: real_lattice(3, 3)
real(kind=dp), intent(in) :: wannier_centres_translated(:, :)
complex(kind=dp), intent(in) :: ham_r(:, :, :)
real(kind=dp), allocatable, intent(in) :: fermi_energy_list(:)
real(kind=dp), intent(in) :: kpt_latt(:, :)
type(print_output_type), intent(in) :: print_output
type(real_space_ham_type), intent(inout) :: real_space_ham
type(transport_type), intent(inout) :: 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
logical, intent(in) :: pl_warning
! local variables
integer :: i, j, k, num_wann_cell_ll, file_unit, ierr, band_size
integer :: fermi_n
real(kind=dp), allocatable :: sub_block(:, :)
real(kind=dp) :: PL_length, dist, dist_vec(3)
character(len=9) :: cdate, ctime
if (print_output%timing_level > 1) call io_stopwatch_start('tran: lcr_2c2_build_ham', timer)
fermi_n = 0
if (allocated(fermi_energy_list)) fermi_n = size(fermi_energy_list)
if (fermi_n > 1) then
call set_error_fatal(error, "Error in tran_lcr_2c2_build_ham: nfermi>1. " &
//"Set the fermi level using the input parameter 'fermi_evel'", comm)
return
end if
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_2c2_build_ham', 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_2c2_build_ham', comm)
return
end if
allocate (hR0(transport%num_ll, transport%num_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hR0 in tran_lcr_2c2_build_ham', comm)
return
end if
allocate (hR1(transport%num_ll, transport%num_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hR1 in tran_lcr_2c2_build_ham', comm)
return
end if
allocate (hLC(transport%num_ll, transport%num_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hLC in tran_lcr_2c2_build_ham', comm)
return
end if
allocate (hCR(transport%num_ll, transport%num_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hCR in tran_lcr_2c2_build_ham', comm)
return
end if
allocate (hC(num_wann - (2*transport%num_ll), num_wann - (2*transport%num_ll)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hC in tran_lcr_2c2_build_ham', comm)
return
end if
!This checks that only the gamma point is used in wannierisation
!This is necessary since this calculation only makes sense if we
!have periodicity over the supercell.
!BGS, I think (0, 0, 0) in kpt_latt should work as well as in kpt_cart
if ((size(kpt_latt, 2) .ne. 1) .and. (kpt_latt(1, 1) .eq. 0.0_dp) &
.and. (kpt_latt(2, 1) .eq. 0.0_dp) .and. (kpt_latt(3, 1) .eq. 0.0_dp)) then
call set_error_fatal(error, 'Calculation must be performed at gamma only', comm)
return
end if
num_wann_cell_ll = transport%num_ll/transport%num_cell_ll
allocate (sub_block(num_wann_cell_ll, num_wann_cell_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating sub_block in tran_lcr_2c2_build_ham', comm)
return
end if
!Build hL0 & hL1
hL0 = 0.0_dp
hL1 = 0.0_dp
!Loop over the sub_blocks corresponding to distinct unit cells inside the principal layer
do i = 1, transport%num_cell_ll
!Each sub_block will be duplicated along the corresponding diagonal. This ensures the correct symmetry for the leads.
sub_block = 0.0_dp
!Extract matrix elements from hr_one_dim needed for hL0 (and lower triangular sub_blocks of hL1)
do j = 1, num_wann_cell_ll
do k = 1, num_wann_cell_ll
sub_block(j, k) = hr_one_dim(tran_sorted_idx(j), tran_sorted_idx((i - 1)*num_wann_cell_ll + k), 0)
end do
end do
!Filling up hL0 sub_block by sub_block
do j = 1, transport%num_cell_ll - i + 1
!Fill diagonal and upper diagonal sub_blocks
hL0((j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll, &
(j - 1)*num_wann_cell_ll + 1 + (i - 1)*num_wann_cell_ll:j*num_wann_cell_ll + (i - 1)*num_wann_cell_ll) = sub_block
!Fill lower diagonal sub_blocks
if (i .gt. 1) then
hL0((j - 1)*num_wann_cell_ll + 1 + (i - 1)*num_wann_cell_ll:j*num_wann_cell_ll + (i - 1)*num_wann_cell_ll, &
(j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll) = transpose(sub_block)
end if
end do
!Filling up non-diagonal hL1 sub_blocks (nothing need be done for i=1)
if (i .gt. 1) then
do j = 1, i - 1
hL1((transport%num_cell_ll - (i - j))*num_wann_cell_ll + 1:(transport%num_cell_ll - (i - 1 - j))*num_wann_cell_ll, &
(j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll) = sub_block
end do
end if
! MS: Get diagonal and upper triangular sublocks for hL1 - use periodic image of PL4
sub_block = 0.0_dp
if (i == 1) then !Do diagonal only
do j = 1, num_wann_cell_ll
do k = 1, num_wann_cell_ll
sub_block(j, k) = hr_one_dim( &
tran_sorted_idx(num_wann - transport%num_ll + j), &
tran_sorted_idx((i - 1)*num_wann_cell_ll + k), 0)
end do
end do
! MS: Now fill subblocks of hL1
do j = 1, transport%num_cell_ll - i + 1
hL1((j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll, &
(j - 1)*num_wann_cell_ll + 1 + (i - 1)*num_wann_cell_ll:j*num_wann_cell_ll + (i - 1)* &
num_wann_cell_ll) = sub_block
end do
end if
end do
!Special case tran_num_cell_ll=1, the diagonal sub-block of hL1 is hL1, so cannot be left as zero
if (transport%num_cell_ll .eq. 1) then
do j = num_wann - num_wann_cell_ll + 1, num_wann
do k = 1, num_wann_cell_ll
hL1(j - num_wann + num_wann_cell_ll, k) = hr_one_dim(tran_sorted_idx(j), tran_sorted_idx(k), 0)
end do
end do
end if
!Build hR0 & hR1
hR0 = 0.0_dp
hR1 = 0.0_dp
!Loop over the sub_blocks corresponding to distinct unit cells inside the principal layer
do i = 1, transport%num_cell_ll
!Each sub_block will be duplicated along the corresponding diagonal. This ensures the correct symmetry for the leads.
sub_block = 0.0_dp
!Extract matrix elements from hr_one_dim needed for hR0 (and lower triangular sub_blocks of hR1)
do j = 1, num_wann_cell_ll
do k = 1, num_wann_cell_ll
sub_block(j, k) = hr_one_dim(tran_sorted_idx(num_wann - i*(num_wann_cell_ll) + j), &
tran_sorted_idx(num_wann - num_wann_cell_ll + k), 0)
end do
end do
!Filling up hR0 sub_block by sub_block
do j = 1, transport%num_cell_ll - i + 1
!Fill diagonal and upper diagonal sub_blocks
hR0((j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll, &
(j - 1)*num_wann_cell_ll + 1 + (i - 1)*num_wann_cell_ll:j*num_wann_cell_ll + (i - 1)*num_wann_cell_ll) = sub_block
!Fill lower diagonal sub_blocks
if (i .gt. 1) then
hR0((j - 1)*num_wann_cell_ll + 1 + (i - 1)*num_wann_cell_ll:j*num_wann_cell_ll + (i - 1)*num_wann_cell_ll, &
(j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll) = transpose(sub_block)
end if
end do
!Filling up non-diagonal hR1 sub_blocks (nothing need be done for i=1)
if (i .gt. 1) then
do j = 1, i - 1
hR1((transport%num_cell_ll - (i - j))*num_wann_cell_ll + 1:(transport%num_cell_ll - (i - 1 - j))*num_wann_cell_ll, &
(j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll) = sub_block
end do
end if
! MS: Get diagonal and upper triangular sublocks for hR1 - use periodic image of PL1
sub_block = 0.0_dp
if (i == 1) then !Do diagonal only
do j = 1, num_wann_cell_ll
do k = 1, num_wann_cell_ll
sub_block(j, k) = hr_one_dim(tran_sorted_idx((i - 1)*num_wann_cell_ll + k), &
tran_sorted_idx(num_wann - transport%num_ll + j), 0)
end do
end do
! MS: Now fill subblocks of hR1
do j = 1, transport%num_cell_ll - i + 1
hR1((j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll, &
(j - 1)*num_wann_cell_ll + 1 + (i - 1)*num_wann_cell_ll:j*num_wann_cell_ll + (i - 1)*num_wann_cell_ll) = sub_block
end do
end if
end do
!Special case tran_num_cell_ll=1, the diagonal sub-block of hR1 is hR1, so cannot be left as zero
if (transport%num_cell_ll .eq. 1) then
do j = 1, num_wann_cell_ll
do k = num_wann - num_wann_cell_ll + 1, num_wann
hR1(k - num_wann + num_wann_cell_ll, j) = hr_one_dim(tran_sorted_idx(j), tran_sorted_idx(k), 0)
end do
end do
end if
!Building hLC
hLC = 0.0_dp
do i = 1, transport%num_ll
do j = transport%num_ll + 1, 2*transport%num_ll
hLC(i, j - transport%num_ll) = hr_one_dim(tran_sorted_idx(i), tran_sorted_idx(j), 0)
end do
end do
!----!
! MS ! Rely on dist_cutoff doing the work here, as it cuts element-wise, not block wise (incorrect)
!----!
! if (transport%num_cell_ll .gt. 1) then
! do j=1,transport%num_cell_ll
! do k=1,transport%num_cell_ll
! if (k .ge. j) then
! hLC((j-1)*num_wann_cell_ll+1:j*num_wann_cell_ll,(k-1)*num_wann_cell_ll+1:k*num_wann_cell_ll)=0.0_dp
! endif
! enddo
! enddo
! endif
!---!
!end!
!---!
!Building hC
hC = 0.0_dp
band_size = 0
if (real_space_ham%dist_cutoff_hc .ne. real_space_ham%dist_cutoff) then
real_space_ham%dist_cutoff = real_space_ham%dist_cutoff_hc
write (stdout, *) 'Applying dist_cutoff_hc to Hamiltonian for construction of hC'
deallocate (hr_one_dim, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating hr_one_dim in tran_lcr_2c2_build_ham', comm)
return
end if
call tran_reduce_hr(real_space_ham, ham_r, hr_one_dim, real_lattice, irvec, mp_grid, &
irvec_max, nrpts, nrpts_one_dim, num_wann, one_dim_vec, &
print_output%timing_level, stdout, timer, error, comm)
if (allocated(error)) return
call tran_cut_hr_one_dim(real_space_ham, transport, print_output, hr_one_dim, real_lattice, &
wannier_centres_translated, mp_grid, irvec_max, num_pl, num_wann, &
one_dim_vec, stdout, timer)
end if
do i = transport%num_ll + 1, num_wann - transport%num_ll
do j = transport%num_ll + 1, num_wann - transport%num_ll
hC(i - transport%num_ll, j - transport%num_ll) = hr_one_dim(tran_sorted_idx(i), tran_sorted_idx(j), 0)
! Impose a ham_cutoff of 1e-4 eV to reduce tran_num_bandc (and in turn hCband, and speed up transport)
if (abs(hC(i - transport%num_ll, j - transport%num_ll)) .lt. 10.0_dp*eps5) then
hC(i - transport%num_ll, j - transport%num_ll) = 0.0_dp
band_size = max(band_size, abs(i - j))
end if
end do
end do
!Building hCR
hCR = 0.0_dp
do i = num_wann - 2*transport%num_ll + 1, num_wann - transport%num_ll
do j = num_wann - transport%num_ll + 1, num_wann
hCR(i - (num_wann - 2*transport%num_ll), j - (num_wann - transport%num_ll)) = &
hr_one_dim(tran_sorted_idx(i), tran_sorted_idx(j), 0)
end do
end do
!----!
! MS ! Rely on dist_cutoff doing the work here, as it cuts element-wise, not block wise (incorrect)
!----!
! if (transport%num_cell_ll .gt. 1) then
! do j=1,transport%num_cell_ll
! do k=1,transport%num_cell_ll
! if (k .ge. j) then
! hCR((j-1)*num_wann_cell_ll+1:j*num_wann_cell_ll,(k-1)*num_wann_cell_ll+1:k*num_wann_cell_ll)=0.0_dp
! endif
! enddo
! enddo
! endif
!---!
!end!
!---!
!Subtract the Fermi energy from the diagonal elements of hC,hL0,hR0
do i = 1, transport%num_ll
hL0(i, i) = hL0(i, i) - fermi_energy_list(1)
hR0(i, i) = hR0(i, i) - fermi_energy_list(1)
end do
do i = 1, num_wann - (2*transport%num_ll)
hC(i, i) = hC(i, i) - fermi_energy_list(1)
end do
!Define tran_num_** parameters that are used later in tran_lcr
transport%num_rr = transport%num_ll
transport%num_lc = transport%num_ll
transport%num_cr = transport%num_ll
transport%num_cc = num_wann - (2*transport%num_ll)
! Set appropriate tran_num_bandc if has not been set (0.0_dp is default value)
if (transport%num_bandc .eq. 0.0_dp) then
transport%num_bandc = min(band_size + 1, (transport%num_cc + 1)/2 + 1)
end if
! MS: Find and print effective PL length
if (.not. pl_warning) then
PL_length = 0.0_dp
do i = 1, transport%num_ll
do j = 1, transport%num_ll
if (abs(hL1(i, j)) .gt. 0.0_dp) then
if (index(real_space_ham%dist_cutoff_mode, 'one_dim') .gt. 0) then
dist = abs(wannier_centres_translated(coord(1), tran_sorted_idx(i)) &
- wannier_centres_translated(coord(1), tran_sorted_idx(j + transport%num_ll)))
else
dist_vec(:) = wannier_centres_translated(:, tran_sorted_idx(i)) &
- wannier_centres_translated(:, tran_sorted_idx(j + transport%num_ll))
dist = sqrt(dot_product(dist_vec, dist_vec))
end if
PL_length = max(PL_length, dist)
end if
if (abs(hR1(i, j)) .gt. 0.0_dp) then
if (index(real_space_ham%dist_cutoff_mode, 'one_dim') .gt. 0) then
dist = abs(wannier_centres_translated(coord(1), tran_sorted_idx(num_wann - 2*transport%num_ll + i)) &
- wannier_centres_translated(coord(1), tran_sorted_idx(num_wann - transport%num_ll + j)))
else
dist_vec(:) = wannier_centres_translated(:, tran_sorted_idx(num_wann - 2*transport%num_ll + i)) &
- wannier_centres_translated(:, tran_sorted_idx(num_wann - transport%num_ll + j))
dist = sqrt(dot_product(dist_vec, dist_vec))
end if
PL_length = max(PL_length, dist)
end if
end do
end do
write (stdout, '(1x,a,f12.6,a)') 'Approximate effective principal layer length is: ', PL_length, ' Ang.'
end if
!Writing to file:
if (transport%write_ht) then
write (stdout, *) '------------------------------- Writing ht files ----------------------------'
open (newunit=file_unit, file=trim(seedname)//'_htL.dat', status='unknown', form='formatted', action='write')
call io_date(cdate, ctime)
write (file_unit, *) 'written on '//cdate//' at '//ctime ! Date and time
write (file_unit, '(I6)') transport%num_ll
write (file_unit, '(6F12.6)') ((hL0(j, i), j=1, transport%num_ll), i=1, transport%num_ll)
write (file_unit, '(I6)') transport%num_ll
write (file_unit, '(6F12.6)') ((hL1(j, i), j=1, transport%num_ll), i=1, transport%num_ll)
close (file_unit)
write (stdout, *) ' '//trim(seedname)//'_htL.dat written'
!hR
open (newunit=file_unit, file=trim(seedname)//'_htR.dat', status='unknown', form='formatted', action='write')
call io_date(cdate, ctime)
write (file_unit, *) 'written on '//cdate//' at '//ctime ! Date and time
write (file_unit, '(I6)') transport%num_rr
write (file_unit, '(6F12.6)') ((hR0(j, i), j=1, transport%num_rr), i=1, transport%num_rr)
write (file_unit, '(I6)') transport%num_rr
write (file_unit, '(6F12.6)') ((hR1(j, i), j=1, transport%num_rr), i=1, transport%num_rr)
close (file_unit)
write (stdout, *) ' '//trim(seedname)//'_htR.dat written'
!hLC
open (newunit=file_unit, file=trim(seedname)//'_htLC.dat', status='unknown', form='formatted', action='write')
call io_date(cdate, ctime)
write (file_unit, *) 'written on '//cdate//' at '//ctime ! Date and time
write (file_unit, '(2I6)') transport%num_ll, transport%num_lc
write (file_unit, '(6F12.6)') ((hLC(j, i), j=1, transport%num_lc), i=1, transport%num_lc)
close (file_unit)
write (stdout, *) ' '//trim(seedname)//'_htLC.dat written'
!hCR
open (newunit=file_unit, file=trim(seedname)//'_htCR.dat', status='unknown', form='formatted', action='write')
call io_date(cdate, ctime)
write (file_unit, *) 'written on '//cdate//' at '//ctime ! Date and time
write (file_unit, '(2I6)') transport%num_cr, transport%num_rr
write (file_unit, '(6F12.6)') ((hCR(j, i), j=1, transport%num_cr), i=1, transport%num_cr)
close (file_unit)
write (stdout, *) ' '//trim(seedname)//'_htCR.dat written'
!hC
open (newunit=file_unit, file=trim(seedname)//'_htC.dat', status='unknown', form='formatted', action='write')
call io_date(cdate, ctime)
write (file_unit, *) 'written on '//cdate//' at '//ctime ! Date and time
write (file_unit, '(I6)') transport%num_cc
write (file_unit, '(6F12.6)') ((hC(j, i), j=1, transport%num_cc), i=1, transport%num_cc)
close (file_unit)
write (stdout, *) ' '//trim(seedname)//'_htC.dat written'
write (stdout, *) '------------------------------------------------------------------------------'
end if
deallocate (sub_block, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating sub_block in tran_lcr_2c2_build_ham', comm)
return
end if
if (print_output%timing_level > 1) call io_stopwatch_stop('tran: lcr_2c2_build_ham', timer)
return
end subroutine tran_lcr_2c2_build_ham
