subroutine tran_lcr_2c2_build_ham(pl_warning)
!==============================================!
! 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_error, stdout, seedname, io_file_unit, io_date, io_stopwatch
use w90_parameters, only: tran_num_cell_ll, num_wann, tran_num_ll, kpt_cart, nfermi, fermi_energy_list, &
tran_write_ht, tran_num_rr, tran_num_lc, tran_num_cr, tran_num_cc, &
tran_num_bandc, timing_level, dist_cutoff_mode, dist_cutoff, &
dist_cutoff_hc
use w90_hamiltonian, only: wannier_centres_translated
implicit none
logical, intent(in) :: pl_warning
integer :: i, j, k, num_wann_cell_ll, file_unit, ierr, band_size
real(dp), allocatable, dimension(:, :) :: sub_block
real(dp) :: PL_length, dist, dist_vec(3)
character(len=9) :: cdate, ctime
if (timing_level > 1) call io_stopwatch('tran: lcr_2c2_build_ham', 1)
if (nfermi > 1) call io_error("Error in tran_lcr_2c2_build_ham: nfermi>1. " &
//"Set the fermi level using the input parameter 'fermi_evel'")
allocate (hL0(tran_num_ll, tran_num_ll), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating hL0 in tran_lcr_2c2_build_ham')
allocate (hL1(tran_num_ll, tran_num_ll), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating hL1 in tran_lcr_2c2_build_ham')
allocate (hR0(tran_num_ll, tran_num_ll), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating hR0 in tran_lcr_2c2_build_ham')
allocate (hR1(tran_num_ll, tran_num_ll), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating hR1 in tran_lcr_2c2_build_ham')
allocate (hLC(tran_num_ll, tran_num_ll), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating hLC in tran_lcr_2c2_build_ham')
allocate (hCR(tran_num_ll, tran_num_ll), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating hCR in tran_lcr_2c2_build_ham')
allocate (hC(num_wann - (2*tran_num_ll), num_wann - (2*tran_num_ll)), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating hC in tran_lcr_2c2_build_ham')
!
!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.
!
if ((size(kpt_cart, 2) .ne. 1) .and. (kpt_cart(1, 1) .eq. 0.0_dp) &
.and. (kpt_cart(2, 1) .eq. 0.0_dp) &
.and. (kpt_cart(3, 1) .eq. 0.0_dp)) then
call io_error('Calculation must be performed at gamma only')
endif
num_wann_cell_ll = tran_num_ll/tran_num_cell_ll
allocate (sub_block(num_wann_cell_ll, num_wann_cell_ll), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating sub_block in tran_lcr_2c2_build_ham')
!
!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, tran_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)
enddo
enddo
!
!Filling up hL0 sub_block by sub_block
!
do j = 1, tran_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)
endif
enddo
!
!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((tran_num_cell_ll - (i - j))*num_wann_cell_ll + 1:(tran_num_cell_ll - (i - 1 - j))*num_wann_cell_ll, &
(j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll) = sub_block
enddo
endif
!
! 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 - tran_num_ll + j), &
tran_sorted_idx((i - 1)*num_wann_cell_ll + k), 0)
enddo
enddo
!
! MS: Now fill subblocks of hL1
!
do j = 1, tran_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
enddo
endif
enddo
!
!Special case tran_num_cell_ll=1, the diagonal sub-block of hL1 is hL1, so cannot be left as zero
!
if (tran_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)
enddo
enddo
endif
!
!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, tran_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)
enddo
enddo
!
!Filling up hR0 sub_block by sub_block
!
do j = 1, tran_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)
endif
enddo
!
!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((tran_num_cell_ll - (i - j))*num_wann_cell_ll + 1:(tran_num_cell_ll - (i - 1 - j))*num_wann_cell_ll, &
(j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll) = sub_block
enddo
endif
!
! 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 - tran_num_ll + j), 0)
enddo
enddo
!
! MS: Now fill subblocks of hR1
!
do j = 1, tran_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
enddo
endif
enddo
!
!Special case tran_num_cell_ll=1, the diagonal sub-block of hR1 is hR1, so cannot be left as zero
!
if (tran_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)
enddo
enddo
endif
!
!Building hLC
!
hLC = 0.0_dp
do i = 1, tran_num_ll
do j = tran_num_ll + 1, 2*tran_num_ll
hLC(i, j - tran_num_ll) = hr_one_dim(tran_sorted_idx(i), tran_sorted_idx(j), 0)
enddo
enddo
!----!
! MS ! Rely on dist_cutoff doing the work here, as it cuts element-wise, not block wise (incorrect)
!----!
! if (tran_num_cell_ll .gt. 1) then
! do j=1,tran_num_cell_ll
! do k=1,tran_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 (dist_cutoff_hc .ne. dist_cutoff) then
dist_cutoff = dist_cutoff_hc
write (stdout, *) 'Applying dist_cutoff_hc to Hamiltonian for construction of hC'
deallocate (hr_one_dim, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating hr_one_dim in tran_lcr_2c2_sort')
call tran_reduce_hr()
call tran_cut_hr_one_dim()
endif
do i = tran_num_ll + 1, num_wann - tran_num_ll
do j = tran_num_ll + 1, num_wann - tran_num_ll
hC(i - tran_num_ll, j - tran_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 - tran_num_ll, j - tran_num_ll)) .lt. 10.0_dp*eps5) then
hC(i - tran_num_ll, j - tran_num_ll) = 0.0_dp
band_size = max(band_size, abs(i - j))
endif
enddo
enddo
!
!Building hCR
!
hCR = 0.0_dp
do i = num_wann - 2*tran_num_ll + 1, num_wann - tran_num_ll
do j = num_wann - tran_num_ll + 1, num_wann
hCR(i - (num_wann - 2*tran_num_ll), j - (num_wann - tran_num_ll)) = hr_one_dim(tran_sorted_idx(i), tran_sorted_idx(j), 0)
enddo
enddo
!----!
! MS ! Rely on dist_cutoff doing the work here, as it cuts element-wise, not block wise (incorrect)
!----!
! if (tran_num_cell_ll .gt. 1) then
! do j=1,tran_num_cell_ll
! do k=1,tran_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, tran_num_ll
hL0(i, i) = hL0(i, i) - fermi_energy_list(1)
hR0(i, i) = hR0(i, i) - fermi_energy_list(1)
enddo
do i = 1, num_wann - (2*tran_num_ll)
hC(i, i) = hC(i, i) - fermi_energy_list(1)
enddo
!
!Define tran_num_** parameters that are used later in tran_lcr
!
tran_num_rr = tran_num_ll
tran_num_lc = tran_num_ll
tran_num_cr = tran_num_ll
tran_num_cc = num_wann - (2*tran_num_ll)
!
! Set appropriate tran_num_bandc if has not been set (0.0_dp is default value)
!
if (tran_num_bandc .eq. 0.0_dp) then
tran_num_bandc = min(band_size + 1, (tran_num_cc + 1)/2 + 1)
endif
!
! MS: Find and print effective PL length
!
if (.not. pl_warning) then
PL_length = 0.0_dp
do i = 1, tran_num_ll
do j = 1, tran_num_ll
if (abs(hL1(i, j)) .gt. 0.0_dp) then
if (index(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 + tran_num_ll)))
else
dist_vec(:) = wannier_centres_translated(:, tran_sorted_idx(i)) &
- wannier_centres_translated(:, tran_sorted_idx(j + tran_num_ll))
dist = sqrt(dot_product(dist_vec, dist_vec))
endif
PL_length = max(PL_length, dist)
endif
if (abs(hR1(i, j)) .gt. 0.0_dp) then
if (index(dist_cutoff_mode, 'one_dim') .gt. 0) then
dist = abs(wannier_centres_translated(coord(1), tran_sorted_idx(num_wann - 2*tran_num_ll + i)) &
- wannier_centres_translated(coord(1), tran_sorted_idx(num_wann - tran_num_ll + j)))
else
dist_vec(:) = wannier_centres_translated(:, tran_sorted_idx(num_wann - 2*tran_num_ll + i)) &
- wannier_centres_translated(:, tran_sorted_idx(num_wann - tran_num_ll + j))
dist = sqrt(dot_product(dist_vec, dist_vec))
endif
PL_length = max(PL_length, dist)
endif
enddo
enddo
write (stdout, '(1x,a,f12.6,a)') 'Approximate effective principal layer length is: ', PL_length, ' Ang.'
endif
!
!Writing to file:
!
if (tran_write_ht) then
write (stdout, *) '------------------------------- Writing ht files ----------------------------'
!
file_unit = io_file_unit()
open (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)') tran_num_ll
write (file_unit, '(6F12.6)') ((hL0(j, i), j=1, tran_num_ll), i=1, tran_num_ll)
write (file_unit, '(I6)') tran_num_ll
write (file_unit, '(6F12.6)') ((hL1(j, i), j=1, tran_num_ll), i=1, tran_num_ll)
close (file_unit)
write (stdout, *) ' '//trim(seedname)//'_htL.dat written'
!
!hR
!
file_unit = io_file_unit()
open (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)') tran_num_rr
write (file_unit, '(6F12.6)') ((hR0(j, i), j=1, tran_num_rr), i=1, tran_num_rr)
write (file_unit, '(I6)') tran_num_rr
write (file_unit, '(6F12.6)') ((hR1(j, i), j=1, tran_num_rr), i=1, tran_num_rr)
close (file_unit)
write (stdout, *) ' '//trim(seedname)//'_htR.dat written'
!
!hLC
!
file_unit = io_file_unit()
open (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)') tran_num_ll, tran_num_lc
write (file_unit, '(6F12.6)') ((hLC(j, i), j=1, tran_num_lc), i=1, tran_num_lc)
close (file_unit)
write (stdout, *) ' '//trim(seedname)//'_htLC.dat written'
!
!hCR
!
file_unit = io_file_unit()
open (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)') tran_num_cr, tran_num_rr
write (file_unit, '(6F12.6)') ((hCR(j, i), j=1, tran_num_cr), i=1, tran_num_cr)
close (file_unit)
write (stdout, *) ' '//trim(seedname)//'_htCR.dat written'
!
!hC
!
file_unit = io_file_unit()
open (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)') tran_num_cc
write (file_unit, '(6F12.6)') ((hC(j, i), j=1, tran_num_cc), i=1, tran_num_cc)
close (file_unit)
write (stdout, *) ' '//trim(seedname)//'_htC.dat written'
write (stdout, *) '------------------------------------------------------------------------------'
end if
deallocate (sub_block, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating sub_block in tran_lcr_2c2_build_ham')
if (timing_level > 1) call io_stopwatch('tran: lcr_2c2_build_ham', 2)
return
end subroutine tran_lcr_2c2_build_ham