subroutine 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)
!================================================!
use w90_constants, only: dp
use w90_io, only: 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
implicit none
! arguments
type(real_space_ham_type), intent(inout) :: real_space_ham
type(print_output_type), intent(in) :: print_output
type(transport_type), intent(inout) :: transport
type(timer_list_type), intent(inout) :: timer
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: irvec_max
integer, intent(in) :: num_wann
integer, intent(in) :: one_dim_vec
integer, intent(inout) :: num_pl
integer, intent(in) :: stdout
real(kind=dp), intent(inout) :: hr_one_dim(:, :, -irvec_max:)
real(kind=dp), intent(in) :: real_lattice(3, 3)
real(kind=dp), intent(in) :: wannier_centres_translated(:, :)
! local variables
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, -irvec_max:irvec_max)
real(kind=dp) :: hr_tmp(num_wann, num_wann)
if (print_output%timing_level > 1) call io_stopwatch_start('tran: cut_hr_one_dim', timer)
!irvec_max = nrpts_one_dim/2 ! now passed as arg
! maximum possible dist_cutoff
dist = real(mp_grid(one_dim_vec), dp)*abs(real_lattice(real_space_ham%one_dim_dir, one_dim_vec)) &
/2.0_dp
if (real_space_ham%dist_cutoff .gt. dist) then
write (stdout, '(1x,a,1x,F10.5,1x,a)') 'dist_cutoff', real_space_ham%dist_cutoff, &
trim(print_output%length_unit), 'is too large'
real_space_ham%dist_cutoff = dist
! aam_2012-04-13
real_space_ham%dist_cutoff_hc = dist
write (stdout, '(4x,a,1x,F10.5,1x,a)') 'reset to', real_space_ham%dist_cutoff, &
trim(print_output%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(real_space_ham%dist_cutoff_mode, 'one_dim') > 0) then
do i = 1, num_wann
do j = 1, num_wann
dist_ij_vec(real_space_ham%one_dim_dir) = wannier_centres_translated(real_space_ham%one_dim_dir, i) &
- wannier_centres_translated(real_space_ham%one_dim_dir, j)
do n1 = -irvec_max, irvec_max
dist_vec(real_space_ham%one_dim_dir) = dist_ij_vec(real_space_ham%one_dim_dir) &
+ shift_vec(real_space_ham%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. &
!~ (transport%num_cell_ll .eq. 1) .and. &
(abs(dist_vec(real_space_ham%one_dim_dir)) .gt. real_space_ham%dist_cutoff)) then
! Move to right
dist_vec(real_space_ham%one_dim_dir) = dist_ij_vec(real_space_ham%one_dim_dir) &
+ real_lattice(real_space_ham%one_dim_dir, one_dim_vec)
! Move to left
if (abs(dist_vec(real_space_ham%one_dim_dir)) .gt. real_space_ham%dist_cutoff) &
dist_vec(real_space_ham%one_dim_dir) = dist_ij_vec(real_space_ham%one_dim_dir) &
- real_lattice(real_space_ham%one_dim_dir, one_dim_vec)
end if
!end MS
dist = abs(dist_vec(real_space_ham%one_dim_dir))
if (dist .gt. real_space_ham%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. &
!~ (transport%num_cell_ll .eq. 1) .and. &
(dist .gt. real_space_ham%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. real_space_ham%dist_cutoff) then
dist_vec(:) = dist_ij_vec(:) - real_lattice(:, one_dim_vec)
dist = sqrt(dot_product(dist_vec, dist_vec))
end if
end if
! End MS
if (dist .gt. real_space_ham%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. real_space_ham%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:', &
transport%num_cell_ll
write (stdout, '(1x,a,I6)') 'Number of Wannier Functions inside the principal layer:', &
transport%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
end if
! 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. real_space_ham%hr_cutoff) hr_one_dim(j, i, n1) = 0.0_dp
end do
end do
end do
if (print_output%timing_level > 1) call io_stopwatch_stop('tran: cut_hr_one_dim', timer)
return
end subroutine tran_cut_hr_one_dim
