subroutine tran_cut_hr_one_dim()
!==================================================================!
!
use w90_constants, only: dp
use w90_io, only: io_stopwatch, stdout
use w90_parameters, only: num_wann, mp_grid, timing_level, real_lattice, &
hr_cutoff, dist_cutoff, dist_cutoff_mode, &
one_dim_dir, length_unit, transport_mode, &
tran_num_cell_ll, tran_num_ll, dist_cutoff_hc
use w90_hamiltonian, only: wannier_centres_translated
implicit none
!
integer :: irvec_max
integer :: i, j, n1
real(kind=dp) :: hr_max
real(kind=dp) :: dist
real(kind=dp) :: dist_vec(3)
real(kind=dp) :: dist_ij_vec(3)
real(kind=dp) :: shift_vec(3, -nrpts_one_dim/2:nrpts_one_dim/2)
real(kind=dp) :: hr_tmp(num_wann, num_wann)
!
if (timing_level > 1) call io_stopwatch('tran: cut_hr_one_dim', 1)
!
irvec_max = nrpts_one_dim/2
! maximum possible dist_cutoff
dist = real(mp_grid(one_dim_vec), dp)*abs(real_lattice(one_dim_dir, one_dim_vec))/2.0_dp
if (dist_cutoff .gt. dist) then
write (stdout, '(1x,a,1x,F10.5,1x,a)') 'dist_cutoff', dist_cutoff, trim(length_unit), 'is too large'
dist_cutoff = dist
! aam_2012-04-13
dist_cutoff_hc = dist
write (stdout, '(4x,a,1x,F10.5,1x,a)') 'reset to', dist_cutoff, trim(length_unit)
end if
do n1 = -irvec_max, irvec_max
shift_vec(:, n1) = real(n1, dp)*(real_lattice(:, one_dim_vec))
! write(stdout,'(a,3f10.6)') 'shift_vec', shift_vec(:,n1)
end do
! apply dist_cutoff first
if (index(dist_cutoff_mode, 'one_dim') > 0) then
do i = 1, num_wann
do j = 1, num_wann
dist_ij_vec(one_dim_dir) = wannier_centres_translated(one_dim_dir, i) - wannier_centres_translated(one_dim_dir, j)
do n1 = -irvec_max, irvec_max
dist_vec(one_dim_dir) = dist_ij_vec(one_dim_dir) + shift_vec(one_dim_dir, n1)
!
!MS: Add special case for lcr: We must not cut the elements that are within
! dist_cutoff under PBC's (single kpt assumed) in order to build
! hamiltonians correctly in tran_2c2_build_hams
!
if ((index(transport_mode, 'lcr') > 0) .and. &
!~ (tran_num_cell_ll .eq. 1) .and. &
(abs(dist_vec(one_dim_dir)) .gt. dist_cutoff)) then
! Move to right
dist_vec(one_dim_dir) = dist_ij_vec(one_dim_dir) + real_lattice(one_dim_dir, one_dim_vec)
! Move to left
if (abs(dist_vec(one_dim_dir)) .gt. dist_cutoff) &
dist_vec(one_dim_dir) = dist_ij_vec(one_dim_dir) - real_lattice(one_dim_dir, one_dim_vec)
endif
!
!end MS
!
dist = abs(dist_vec(one_dim_dir))
if (dist .gt. dist_cutoff) hr_one_dim(j, i, n1) = 0.0_dp
end do
end do
end do
else
do i = 1, num_wann
do j = 1, num_wann
dist_ij_vec(:) = wannier_centres_translated(:, i) - wannier_centres_translated(:, j)
do n1 = -irvec_max, irvec_max
dist_vec(:) = dist_ij_vec(:) + shift_vec(:, n1)
dist = sqrt(dot_product(dist_vec, dist_vec))
!
! MS: Special case (as above) equivalent for alternate definition of cut off
!
if ((index(transport_mode, 'lcr') > 0) .and. &
!~ (tran_num_cell_ll .eq. 1) .and. &
(dist .gt. dist_cutoff)) then
! Move to right
dist_vec(:) = dist_ij_vec(:) + real_lattice(:, one_dim_vec)
dist = sqrt(dot_product(dist_vec, dist_vec))
! Move to left
if (dist .gt. dist_cutoff) then
dist_vec(:) = dist_ij_vec(:) - real_lattice(:, one_dim_vec)
dist = sqrt(dot_product(dist_vec, dist_vec))
endif
endif
!
! End MS
!
if (dist .gt. dist_cutoff) hr_one_dim(j, i, n1) = 0.0_dp
end do
end do
end do
end if
! output maximum to check a decay of H as a function of lattice vector R
write (stdout, '(/1x,a78)') repeat('-', 78)
write (stdout, '(1x,4x,a)') &
'Maximum real part of the real-space Hamiltonian at each lattice point'
write (stdout, '(1x,8x,a62)') repeat('-', 62)
write (stdout, '(1x,11x,a,11x,a)') 'Lattice point R', 'Max |H_ij(R)|'
! calculate number of units inside a principal layer
num_pl = 0
do n1 = -irvec_max, irvec_max
hr_tmp(:, :) = abs(hr_one_dim(:, :, n1))
hr_max = maxval(hr_tmp)
if (hr_max .gt. hr_cutoff) then
if (abs(n1) .gt. num_pl) num_pl = abs(n1)
else
hr_one_dim(:, :, n1) = 0.0_dp
end if
write (stdout, '(1x,9x,5x,I5,5x,12x,F12.6)') n1, hr_max
end do
write (stdout, '(1x,8x,a62)') repeat('-', 62)
if (index(transport_mode, 'lcr') > 0) then
write (stdout, '(/1x,a,I6)') 'Number of unit cells inside the principal layer:', tran_num_cell_ll
write (stdout, '(1x,a,I6)') 'Number of Wannier Functions inside the principal layer:', tran_num_ll
elseif (index(transport_mode, 'bulk') > 0) then
write (stdout, '(/1x,a,I6)') 'Number of unit cells inside the principal layer:', num_pl
write (stdout, '(1x,a,I6)') 'Number of Wannier Functions inside the principal layer:', num_pl*num_wann
endif
! apply hr_cutoff to each element inside the principal layer
do n1 = -num_pl, num_pl
do i = 1, num_wann
do j = 1, num_wann
if (abs(hr_one_dim(j, i, n1)) .lt. hr_cutoff) hr_one_dim(j, i, n1) = 0.0_dp
end do
end do
end do
if (timing_level > 1) call io_stopwatch('tran: cut_hr_one_dim', 2)
return
end subroutine tran_cut_hr_one_dim