!-*- mode: F90 -*-!
!------------------------------------------------------------!
! Copyright (C) 2026 Wannier Developer Group !
! !
! This library is free software; you can redistribute it !
! and/or modify it under the terms of the GNU Lesser General !
! Public License as published by the Free Software !
! Foundation; either version 2.1 of the License, or (at your !
! option) any later version. !
! !
! This library is distributed in the hope that it will be !
! useful,but WITHOUT ANY WARRANTY; without even the implied !
! warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR !
! PURPOSE. See the GNU Lesser General Public License for !
! more details. !
! !
! You should have received a copy of the GNU Lesser General !
! Public License along with this library; if not, see !
! <https://www.gnu.org/licenses/>. !
! !
! The webpage of the Wannier90 code is !
! <https://www.wannier.org>. !
! !
! The Wannier90 code is hosted on GitHub !
! <https://github.com/wannier-developers/wannier90> !
!------------------------------------------------------------!
!-----------------------------------------------------------------------!
! Based on !
! < dosqc_1.0 > !
! Density Of States and Quantum Conductance - Version 1.0 !
! Marco Buongiorno Nardelli, January 2000. !
! !
! Reference: !
! - M. Buongiorno Nardelli, "Electronic transport in extended systems: !
! application to carbon nanotubes", Phys. Rev. B, vol. 60(11), 7828 !
! (1999) !
!-----------------------------------------------------------------------!
!=======================================================================!
! Definition of parameters used in w90_transport !
!=======================================================================!
! !
! transport_mode = 'bulk' or 'lcr' !
! tran_win_min = minimum E !
! tran_win_max = maximum E !
! tran_energy_step = delta E !
! tran_num_bb = # of WFs in a principal layer of a perfectly !
! periodic bulk system !
! tran_num_ll = # of WFs in a principal layer of a left lead !
! tran_num_rr = # of WFs in a principal layer of a right lead !
! tran_num_cc = # of WFs in a disordered conducter cell !
! tran_num_lc = # of WFs in a disordered conducter cell that !
! are used to calculate interaction with a !
! left lead !
! tran_num_cr = # of WFs in a disordered conducter cell that !
! are used to calculate interaction with a !
! right lead !
! tran_num_bandc = width of band-diagonal hC matrix !
! tran_read_ht = .true. => read H matrix from h*.dat files !
! tran_write_ht = .true. => write H matrix from h*.dat files !
! tran_use_same_lead = .true. => in L-C-R construction, left and !
! right lead are the same kind !
! tran_num_cell_ll = # of unit cells in a PL of left lead !
! tran_num_cell_rr = # of unit cells in a PL of right lead !
! (equal to tran_num_cell_ll for now) !
! tran_group_threshold = distance defining the grouping of WFs !
!=======================================================================!
module w90_transport_mod
!! Module to handle ballistic transport.
!! Based on
!! < dosqc_1.0 >
!! Density Of States and Quantum Conductance - Version 1.0
!! Marco Buongiorno Nardelli, January 2000.
!! Reference:
!! - M. Buongiorno Nardelli, "Electronic transport in extended systems:
!! application to carbon nanotubes", Phys. Rev. B, vol. 60(11), 7828
!! (1999)
use w90_constants, only: dp
use w90_error
use w90_comms
implicit none
private
complex(kind=dp), parameter :: eta = (0.0_dp, 0.0005_dp) !! small complex number
integer, parameter :: nterx = 50 !! nterx = # of maximum iteration to calculate transfer matrix
public :: tran_main
contains
!================================================!
subroutine tran_main(atom_data, dis_manifold, fermi_energy_list, ham_logical, kpt_latt, &
output_file, real_space_ham, transport, print_output, wannier_data, &
ws_region, w90_calculation, ham_k, ham_r, u_matrix, u_matrix_opt, eigval, &
real_lattice, wannier_centres_translated, irvec, mp_grid, ndegen, &
shift_vec, nrpts, num_bands, num_kpts, num_wann, rpt_origin, &
bands_plot_mode, have_disentangled, lsitesymmetry, seedname, stdout, &
timer, error, comm)
!================================================!
!
!! Main transport subroutine
!
!================================================!
use w90_io, only: io_stopwatch_start, io_stopwatch_stop
use w90_error, only: w90_error_type, set_error_dealloc
use w90_hamiltonian, only: hamiltonian_get_hr, hamiltonian_write_hr, hamiltonian_setup
use w90_types, only: wannier_data_type, print_output_type, ws_region_type, &
atom_data_type, dis_manifold_type, timer_list_type
use w90_wannier90_types, only: w90_calculation_type, transport_type, output_file_type, &
real_space_ham_type, ham_logical_type
implicit none
! arguments
type(transport_type), intent(inout) :: transport
type(real_space_ham_type), intent(inout) :: real_space_ham
type(ws_region_type), intent(inout) :: ws_region
type(print_output_type), intent(in) :: print_output
type(w90_calculation_type), intent(in) :: w90_calculation
type(output_file_type), intent(in) :: output_file
type(wannier_data_type), intent(in) :: wannier_data
type(atom_data_type), intent(in) :: atom_data
type(dis_manifold_type), intent(in) :: dis_manifold
real(kind=dp), intent(in) :: kpt_latt(:, :)
real(kind=dp), allocatable, intent(in) :: fermi_energy_list(:)
type(ham_logical_type), intent(inout) :: ham_logical
type(timer_list_type), intent(inout) :: timer
type(w90_error_type), allocatable, intent(out) :: error
type(w90_comm_type), intent(in) :: comm
integer, intent(inout) :: rpt_origin
integer, intent(inout) :: nrpts
integer, intent(inout), allocatable :: ndegen(:)
integer, intent(inout), allocatable :: shift_vec(:, :)
integer, intent(inout), allocatable :: irvec(:, :)
integer, intent(in) :: num_wann
integer, intent(in) :: num_bands
integer, intent(in) :: num_kpts
integer, intent(in) :: stdout
integer, intent(in) :: mp_grid(3)
real(kind=dp), intent(inout), allocatable :: wannier_centres_translated(:, :)
real(kind=dp), intent(in) :: real_lattice(3, 3)
real(kind=dp), intent(in) :: eigval(:, :)
complex(kind=dp), intent(in) :: u_matrix(:, :, :)
complex(kind=dp), intent(in) :: u_matrix_opt(:, :, :)
complex(kind=dp), allocatable, intent(inout) :: ham_k(:, :, :)
complex(kind=dp), intent(inout), allocatable :: ham_r(:, :, :)
character(len=*), intent(in) :: bands_plot_mode
character(len=50), intent(in) :: seedname
logical, intent(in) :: have_disentangled
logical, intent(in) :: lsitesymmetry !YN:
! local variables
integer :: one_dim_vec
!! cartesian axis to which real_lattice(:,one_dim_vec) is parallel
integer :: nrpts_one_dim
integer :: num_pl
!! number of unit cell in a principal layer
integer :: coord(3)
!! coord : coord(1) defines the conduction direction according to
!1=x,2=y,3=z
!! coord(2),coord(3) define the other directions during sorting routines
integer, allocatable :: tran_sorted_idx(:)
!! index of sorted WF centres to unsorted
integer :: num_G
integer :: irvec_max ! size of hr_one_dim's last dimension (:, :, -irvec_max:irvec_max)
integer :: ierr
real(kind=dp), allocatable :: hr_one_dim(:, :, :)
real(kind=dp), allocatable :: signatures(:, :)
real(kind=dp), allocatable :: hB0(:, :)
real(kind=dp), allocatable :: hB1(:, :)
real(kind=dp), allocatable :: hC(:, :)
real(kind=dp), allocatable :: hCR(:, :)
real(kind=dp), allocatable :: hL0(:, :)
real(kind=dp), allocatable :: hL1(:, :)
real(kind=dp), allocatable :: hLC(:, :)
real(kind=dp), allocatable :: hR0(:, :)
real(kind=dp), allocatable :: hR1(:, :)
logical :: pl_warning
! fixme printout guards as elsewhere please (even if not yet parallel)
if (print_output%timing_level > 0) call io_stopwatch_start('tran: main', timer)
write (stdout, '(/1x,a)') '*---------------------------------------------------------------------------*'
write (stdout, '(1x,a)') '| TRANSPORT |'
write (stdout, '(1x,a)') '*---------------------------------------------------------------------------*'
write (stdout, *)
if (index(transport%mode, 'bulk') > 0) then
write (stdout, '(/1x,a/)') 'Calculation of Quantum Conductance and DoS: bulk mode'
if (.not. transport%read_ht) then
call hamiltonian_setup(ham_logical, print_output, ws_region, w90_calculation, ham_k, &
ham_r, real_lattice, wannier_centres_translated, irvec, mp_grid, &
ndegen, num_kpts, num_wann, nrpts, rpt_origin, bands_plot_mode, &
stdout, timer, error, transport%mode, comm)
if (allocated(error)) return
call hamiltonian_get_hr(atom_data, dis_manifold, ham_logical, real_space_ham, &
print_output, ham_k, ham_r, u_matrix, u_matrix_opt, eigval, &
kpt_latt, real_lattice, wannier_data%centres, &
wannier_centres_translated, irvec, shift_vec, nrpts, num_bands, &
num_kpts, num_wann, have_disentangled, stdout, timer, error, &
lsitesymmetry, comm)
if (allocated(error)) return
if (output_file%write_hr) then
! is this redundant after the plot call?
call hamiltonian_write_hr(ham_r, irvec, ndegen, nrpts, num_wann, &
print_output%timing_level, seedname, timer, error, comm)
end if
if (allocated(error)) return
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)
call tran_get_ht(fermi_energy_list, transport, hB0, hB1, hr_one_dim, irvec_max, num_pl, &
num_wann, print_output%timing_level, seedname, timer, error, comm)
if (allocated(error)) return
if (output_file%write_xyz) then
call tran_write_xyz(atom_data, transport, wannier_centres_translated, tran_sorted_idx, &
num_wann, seedname, stdout)
end if
end if
call tran_bulk(transport, hB0, hB1, print_output%timing_level, stdout, seedname, timer, error, comm)
if (allocated(error)) return
end if
if (index(transport%mode, 'lcr') > 0) then
write (stdout, '(/1x,a/)') 'Calculation of Quantum Conductance and DoS: lead-conductor-lead mode'
if (.not. transport%read_ht) then
call hamiltonian_setup(ham_logical, print_output, ws_region, w90_calculation, ham_k, &
ham_r, real_lattice, wannier_centres_translated, irvec, mp_grid, &
ndegen, num_kpts, num_wann, nrpts, rpt_origin, bands_plot_mode, &
stdout, timer, error, transport%mode, comm)
if (allocated(error)) return
call hamiltonian_get_hr(atom_data, dis_manifold, ham_logical, real_space_ham, &
print_output, ham_k, ham_r, u_matrix, u_matrix_opt, eigval, &
kpt_latt, real_lattice, wannier_data%centres, &
wannier_centres_translated, irvec, shift_vec, nrpts, num_bands, &
num_kpts, num_wann, have_disentangled, stdout, timer, error, &
lsitesymmetry, comm)
if (allocated(error)) return
if (output_file%write_hr) then
call hamiltonian_write_hr(ham_r, irvec, ndegen, nrpts, num_wann, &
print_output%timing_level, seedname, timer, error, comm)
if (allocated(error)) 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)
write (stdout, *) '------------------------- 2c2 Calculation Type: ------------------------------'
write (stdout, *) ' '
call tran_find_integral_signatures(signatures, num_G, print_output, real_lattice, &
u_matrix_opt, u_matrix, num_bands, num_wann, &
have_disentangled, wannier_centres_translated, stdout, &
seedname, timer, error, comm)
if (allocated(error)) return
call tran_lcr_2c2_sort(signatures, num_G, pl_warning, transport, atom_data, wannier_data, &
real_space_ham, print_output, real_lattice, num_wann, mp_grid, &
ham_r, irvec, nrpts, wannier_centres_translated, one_dim_vec, &
nrpts_one_dim, num_pl, coord, tran_sorted_idx, hr_one_dim, &
irvec_max, output_file%write_xyz, stdout, seedname, timer, error, comm)
if (allocated(error)) return
if (output_file%write_xyz) call tran_write_xyz(atom_data, transport, &
wannier_centres_translated, &
tran_sorted_idx, num_wann, seedname, stdout)
call tran_parity_enforce(signatures, print_output, transport, num_wann, tran_sorted_idx, &
hr_one_dim, irvec_max, stdout, timer)
call 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)
if (allocated(error)) return
end if
call tran_lcr(transport, hC, hCR, hL0, hL1, hLC, hR0, hR1, print_output%timing_level, &
stdout, seedname, timer, error, comm)
if (allocated(error)) return
end if
if (print_output%timing_level > 0) call io_stopwatch_stop('tran: main', timer)
if (allocated(hR1)) then
deallocate (hR1, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating hR1 in tran_main', comm)
return
end if
end if
if (allocated(hR0)) then
deallocate (hR0, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating hR0 in tran_main', comm)
return
end if
end if
if (allocated(hL1)) then
deallocate (hL1, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating hL1 in tran_main', comm)
return
end if
end if
if (allocated(hB1)) then
deallocate (hB1, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating hB1 in tran_main', comm)
return
end if
end if
if (allocated(hB0)) then
deallocate (hB0, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating hB0 in tran_main', comm)
return
end if
end if
if (allocated(hr_one_dim)) then
deallocate (hr_one_dim, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating hr_one_dim in tran_main', comm)
return
end if
end if
end subroutine tran_main
!================================================!
subroutine 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, timing_level, &
stdout, timer, error, comm)
!================================================!
!
! reduce ham_r from 3-d to 1-d
!
!================================================!
use w90_constants, only: dp, eps8
use w90_io, only: io_stopwatch_start, io_stopwatch_stop
use w90_wannier90_types, only: real_space_ham_type
use w90_error, only: w90_error_type, set_error_alloc, set_error_fatal
use w90_types, only: timer_list_type
implicit none
! passed vars
type(real_space_ham_type), intent(in) :: real_space_ham
type(timer_list_type), intent(inout) :: timer
type(w90_error_type), allocatable, intent(out) :: error
type(w90_comm_type), intent(in) :: comm
integer, intent(in) :: irvec(:, :)
integer, intent(inout) :: irvec_max ! limits of hr_one_dim final dim
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: nrpts
integer, intent(in) :: num_wann
integer, intent(inout) :: nrpts_one_dim
integer, intent(inout) :: one_dim_vec
integer, intent(in) :: timing_level
integer, intent(in) :: stdout
real(kind=dp), allocatable, intent(inout) :: hr_one_dim(:, :, :)
real(kind=dp), intent(in) :: real_lattice(3, 3)
complex(kind=dp), intent(in) :: ham_r(:, :, :)
! local variables
integer :: ierr
integer :: irvec_tmp(3), two_dim_vec(2)
integer :: i, j
integer :: i1, i2, i3, n1, nrpts_tmp, loop_rpt
if (timing_level > 1) call io_stopwatch_start('tran: reduce_hr', timer)
! Find one_dim_vec which is parallel to one_dim_dir
! two_dim_vec - the other two lattice vectors
j = 0
do i = 1, 3
if (abs(abs(real_lattice(real_space_ham%one_dim_dir, i)) &
- sqrt(dot_product(real_lattice(:, i), real_lattice(:, i)))) .lt. eps8) then
one_dim_vec = i
j = j + 1
end if
end do
if (j .ne. 1) then
write (stdout, '(i3,a)') j, ' : 1-D LATTICE VECTOR NOT DEFINED'
call set_error_fatal(error, 'Error: 1-d lattice vector not defined in tran_reduce_hr', comm)
return
end if
j = 0
do i = 1, 3
if (i .ne. one_dim_vec) then
j = j + 1
two_dim_vec(j) = i
end if
end do
! starting H matrix should include all W-S supercell where
! the center of the cell spans the full space of the home cell
! adding one more buffer layer when mp_grid(one_dim_vec) is an odd number
!irvec_max = (mp_grid(one_dim_vec)+1)/2
irvec_tmp = maxval(irvec, DIM=2) + 1
irvec_max = irvec_tmp(one_dim_vec)
nrpts_one_dim = 2*irvec_max + 1
allocate (hr_one_dim(num_wann, num_wann, -irvec_max:irvec_max), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hr_one_dim in tran_reduce_hr', comm)
return
end if
hr_one_dim = 0.0_dp
! check imaginary part
write (stdout, '(1x,a,F12.6)') 'Maximum imaginary part of the real-space Hamiltonian: ', &
maxval(abs(aimag(ham_r)))
! select a subset of ham_r, where irvec is 0 along the two other lattice vectors
nrpts_tmp = 0
loop_n1: do n1 = -irvec_max, irvec_max
do loop_rpt = 1, nrpts
i1 = mod(n1 - irvec(one_dim_vec, loop_rpt), mp_grid(one_dim_vec))
i2 = irvec(two_dim_vec(1), loop_rpt)
i3 = irvec(two_dim_vec(2), loop_rpt)
if (i1 .eq. 0 .and. i2 .eq. 0 .and. i3 .eq. 0) then
nrpts_tmp = nrpts_tmp + 1
hr_one_dim(:, :, n1) = real(ham_r(:, :, loop_rpt), dp)
cycle loop_n1
end if
end do
end do loop_n1
if (nrpts_tmp .ne. nrpts_one_dim) then
write (stdout, '(a)') 'FAILED TO EXTRACT 1-D HAMILTONIAN'
call set_error_fatal(error, 'Error: cannot extract 1d hamiltonian in tran_reduce_hr', comm)
return
end if
if (timing_level > 1) call io_stopwatch_stop('tran: reduce_hr', timer)
return
end subroutine tran_reduce_hr
!================================================!
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
!================================================!
subroutine tran_get_ht(fermi_energy_list, transport, hB0, hB1, hr_one_dim, irvec_max, num_pl, &
num_wann, timing_level, seedname, timer, error, comm)
!================================================!
!
!! Construct h00 and h01
!
!================================================!
use w90_constants, only: dp
use w90_io, only: io_stopwatch_start, io_stopwatch_stop, io_date
use w90_wannier90_types, only: transport_type
use w90_error, only: w90_error_type, set_error_alloc, set_error_fatal
use w90_types, only: timer_list_type
implicit none
! arguments
real(kind=dp), allocatable, intent(in) :: fermi_energy_list(:)
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
integer, intent(in) :: num_pl
integer, intent(in) :: num_wann
integer, intent(in) :: irvec_max
integer, intent(in) :: timing_level
real(kind=dp), intent(in) :: hr_one_dim(:, :, -irvec_max:)
real(kind=dp), allocatable, intent(inout) :: hB0(:, :)
real(kind=dp), allocatable, intent(inout) :: hB1(:, :)
character(len=50), intent(in) :: seedname
! local variables
integer :: ierr, file_unit
integer :: i, j, n1, im, jm
integer :: fermi_n
character(len=9) :: cdate, ctime
if (timing_level > 1) call io_stopwatch_start('tran: get_ht', 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_get_ht: nfermi>1. " &
//"Set the fermi level using the input parameter 'fermi_evel'", comm)
return
end if
transport%num_bb = num_pl*num_wann
allocate (hB0(transport%num_bb, transport%num_bb), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hB0 in tran_get_ht', comm)
return
end if
allocate (hB1(transport%num_bb, transport%num_bb), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hB1 in tran_get_ht', comm)
return
end if
hB0 = 0.0_dp
hB1 = 0.0_dp
! h00
do j = 0, num_pl - 1
do i = 0, num_pl - 1
n1 = i - j
im = i*num_wann
jm = j*num_wann
hB0(jm + 1:jm + num_wann, im + 1:im + num_wann) = hr_one_dim(:, :, n1)
end do
end do
! h01
do j = 1, num_pl
do i = 0, j - 1
n1 = i - (j - 1) + num_pl
im = i*num_wann
jm = (j - 1)*num_wann
hB1(jm + 1:jm + num_wann, im + 1:im + num_wann) = hr_one_dim(:, :, n1)
end do
end do
! shift by fermi_energy
do i = 1, transport%num_bb
hB0(i, i) = hB0(i, i) - fermi_energy_list(1)
end do
if (transport%write_ht) then
open (newunit=file_unit, file=trim(seedname)//'_htB.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_bb
write (file_unit, '(6F12.6)') ((hB0(j, i), j=1, transport%num_bb), i=1, transport%num_bb)
write (file_unit, '(I6)') transport%num_bb
write (file_unit, '(6F12.6)') ((hB1(j, i), j=1, transport%num_bb), i=1, transport%num_bb)
close (file_unit)
end if
if (timing_level > 1) call io_stopwatch_stop('tran: get_ht', timer)
return
end subroutine tran_get_ht
!================================================!
subroutine tran_bulk(transport, hB0, hB1, timing_level, stdout, seedname, timer, error, comm)
!================================================!
use w90_constants, only: dp, cmplx_0, cmplx_1, cmplx_i, pi
use w90_io, only: io_stopwatch_start, io_stopwatch_stop, io_date
use w90_wannier90_types, only: transport_type
use w90_error, only: w90_error_type, set_error_alloc, set_error_dealloc
use w90_types, only: timer_list_type
implicit none
! arguments
integer, intent(in) :: stdout
integer, intent(in) :: timing_level
real(kind=dp), allocatable :: hB0(:, :)
real(kind=dp), allocatable :: hB1(:, :)
type(transport_type), intent(in) :: 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
! local variables
integer :: qc_unit, dos_unit
integer :: ierr
integer :: n_e, n, i
real(kind=dp) :: qc, dos
real(kind=dp) :: e_scan
complex(kind=dp) :: e_scan_cmp
complex(kind=dp), allocatable, dimension(:, :) :: tot, tott
complex(kind=dp), allocatable, dimension(:, :) :: g_B, gR, gL
complex(kind=dp), allocatable, dimension(:, :) :: sLr, sRr
complex(kind=dp), allocatable, dimension(:, :) :: s1, s2, c1
character(len=50) :: filename
character(len=9) :: cdate, ctime
if (timing_level > 1) call io_stopwatch_start('tran: bulk', timer)
allocate (tot(transport%num_bb, transport%num_bb), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating tot in tran_bulk', comm)
return
end if
allocate (tott(transport%num_bb, transport%num_bb), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating tott in tran_bulk', comm)
return
end if
allocate (g_B(transport%num_bb, transport%num_bb), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating g_B in tran_bulk', comm)
return
end if
allocate (gL(transport%num_bb, transport%num_bb), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating gL in tran_bulk', comm)
return
end if
allocate (gR(transport%num_bb, transport%num_bb), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating gR in tran_bulk', comm)
return
end if
allocate (sLr(transport%num_bb, transport%num_bb), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating sLr in tran_bulk', comm)
return
end if
allocate (sRr(transport%num_bb, transport%num_bb), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating sRr in tran_bulk', comm)
return
end if
allocate (s1(transport%num_bb, transport%num_bb), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating s1 in tran_bulk', comm)
return
end if
allocate (s2(transport%num_bb, transport%num_bb), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating s2 in tran_bulk', comm)
return
end if
allocate (c1(transport%num_bb, transport%num_bb), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating c1 in tran_bulk', comm)
return
end if
call io_date(cdate, ctime)
open (newunit=qc_unit, file=trim(seedname)//'_qc.dat', status='unknown', &
form='formatted', action='write')
write (qc_unit, *) '## written on '//cdate//' at '//ctime ! Date and time
open (newunit=dos_unit, file=trim(seedname)//'_dos.dat', status='unknown', &
form='formatted', action='write')
write (dos_unit, *) '## written on '//cdate//' at '//ctime ! Date and time
! set up the layer hamiltonians
if (transport%read_ht) then
allocate (hB0(transport%num_bb, transport%num_bb), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hB0 in tran_bulk', comm)
return
end if
allocate (hB1(transport%num_bb, transport%num_bb), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hB1 in tran_bulk', comm)
return
end if
filename = trim(seedname)//'_htB.dat'
call tran_read_htX(transport%num_bb, hB0, hB1, filename, stdout, error, comm)
if (allocated(error)) return
end if
! loop over the energies
n_e = floor((transport%win_max - transport%win_min)/transport%energy_step) + 1
write (stdout, '(/1x,a)', advance='no') 'Calculating quantum&
& conductance and density of states...'
do n = 1, n_e
e_scan = transport%win_min + real(n - 1, dp)*transport%energy_step
! if (mod(n,nint(0.1*n_e)).eq.0) write(stdout,'(a)',advance='no') '.'
! compute conductance according to Fisher and Lee
! retarded Green
e_scan_cmp = e_scan + eta
call tran_transfer(tot, tott, hB0, hB1, e_scan_cmp, transport%num_bb, stdout, error, comm)
if (allocated(error)) return
call tran_green(tot, tott, hB0, hB1, e_scan, g_B, 0, 1, transport%num_bb, stdout, error, comm)
if (allocated(error)) return
! compute S_Lr and S_Rr
c1(:, :) = cmplx(hB1(:, :), kind=dp)
! Self-energy (Sigma_L^r) : sLr = (hB1)^+ * tott
! Self-energy (Sigma_R^r) : sRr = (hB1) * tot
sLr = cmplx_0
sRr = cmplx_0
call ZGEMM('C', 'N', transport%num_bb, transport%num_bb, transport%num_bb, cmplx_1, c1, &
transport%num_bb, tott, transport%num_bb, cmplx_0, sLr, transport%num_bb)
call ZGEMM('N', 'N', transport%num_bb, transport%num_bb, transport%num_bb, cmplx_1, c1, &
transport%num_bb, tot, transport%num_bb, cmplx_0, sRr, transport%num_bb)
! Gamma_L = i(Sigma_L^r-Sigma_L^a)
gL = cmplx_i*(sLr - conjg(transpose(sLr)))
! Gamma_R = i(Sigma_R^r-Sigma_R^a)
gR = cmplx_i*(sRr - conjg(transpose(sRr)))
s1 = cmplx_0
s2 = cmplx_0
c1 = cmplx_0
! s1 = Gamma_L * g_B^r
call ZGEMM('N', 'N', transport%num_bb, transport%num_bb, transport%num_bb, cmplx_1, gL, &
transport%num_bb, g_B, transport%num_bb, cmplx_0, s1, transport%num_bb)
! s2 = Gamma_L * g_B^r * Gamma_R
call ZGEMM('N', 'N', transport%num_bb, transport%num_bb, transport%num_bb, cmplx_1, s1, &
transport%num_bb, gR, transport%num_bb, cmplx_0, s2, transport%num_bb)
! c1 = Gamma_L * g_B^r * Gamma_R * g_B^a
call ZGEMM('N', 'C', transport%num_bb, transport%num_bb, transport%num_bb, cmplx_1, s2, &
transport%num_bb, g_B, transport%num_bb, cmplx_0, c1, transport%num_bb)
qc = 0.0_dp
do i = 1, transport%num_bb
qc = qc + real(c1(i, i), dp)
end do
write (qc_unit, '(f15.9,f18.9)') e_scan, qc
dos = 0.0_dp
do i = 1, transport%num_bb
dos = dos - aimag(g_B(i, i))
end do
dos = dos/pi
write (dos_unit, '(f15.9,f18.9)') e_scan, dos
end do
write (stdout, '(a/)') ' done'
close (qc_unit)
close (dos_unit)
deallocate (c1, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating c1 in tran_bulk', comm)
return
end if
deallocate (s2, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating s2 in tran_bulk', comm)
return
end if
deallocate (s1, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating s1 in tran_bulk', comm)
return
end if
deallocate (sRr, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating sRr in tran_bulk', comm)
return
end if
deallocate (sLr, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating sLr in tran_bulk', comm)
return
end if
deallocate (gR, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating gR in tran_bulk', comm)
return
end if
deallocate (gL, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating gL in tran_bulk', comm)
return
end if
deallocate (g_B, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating g_B in tran_bulk', comm)
return
end if
deallocate (tott, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating tott in tran_bulk', comm)
return
end if
deallocate (tot, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating tot in tran_bulk', comm)
return
end if
if (timing_level > 1) call io_stopwatch_stop('tran: bulk', timer)
return
end subroutine tran_bulk
!================================================!
subroutine tran_lcr(transport, hC, hCR, hL0, hL1, hLC, hR0, hR1, timing_level, stdout, seedname, &
timer, error, comm)
!================================================!
use w90_constants, only: dp, cmplx_0, cmplx_1, cmplx_i, pi
use w90_io, only: io_stopwatch_start, io_stopwatch_stop, io_date
use w90_wannier90_types, only: transport_type
use w90_error, only: w90_error_type, set_error_alloc, set_error_dealloc, set_error_fatal
use w90_types, only: timer_list_type
implicit none
! arguments
integer, intent(in) :: stdout
integer, intent(in) :: timing_level
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(:, :)
type(transport_type), intent(in) :: 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
! local variables
integer :: qc_unit, dos_unit
integer :: ierr
integer :: KL, KU, KC
integer :: n_e, n, i, j, k, info
integer, allocatable :: ipiv(:)
real(kind=dp) :: qc, dos
real(kind=dp) :: e_scan
real(kind=dp), allocatable :: hCband(:, :)
complex(kind=dp) :: e_scan_cmp
complex(kind=dp), allocatable :: hLC_cmp(:, :), hCR_cmp(:, :)
complex(kind=dp), allocatable :: totL(:, :), tottL(:, :), totR(:, :), tottR(:, :)
complex(kind=dp), allocatable :: g_surf_L(:, :), g_surf_R(:, :)
complex(kind=dp), allocatable :: g_C(:, :), g_C_inv(:, :)
complex(kind=dp), allocatable :: gR(:, :), gL(:, :)
complex(kind=dp), allocatable :: sLr(:, :), sRr(:, :)
complex(kind=dp), allocatable :: s1(:, :), s2(:, :)
complex(kind=dp), allocatable :: c1(:, :), c2(:, :)
character(len=50) :: filename
character(len=9) :: cdate, ctime
if (timing_level > 1) call io_stopwatch_start('tran: lcr', timer)
call io_date(cdate, ctime)
open (newunit=qc_unit, file=trim(seedname)//'_qc.dat', status='unknown', &
form='formatted', action='write')
write (qc_unit, *) '## written on '//cdate//' at '//ctime ! Date and time
open (newunit=dos_unit, file=trim(seedname)//'_dos.dat', status='unknown', &
form='formatted', action='write')
write (dos_unit, *) '## written on '//cdate//' at '//ctime ! Date and time
KL = max(transport%num_lc, transport%num_cr, transport%num_bandc) - 1
KU = KL
KC = max(transport%num_lc, transport%num_cr)
allocate (hCband(2*KL + KU + 1, transport%num_cc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hCband in tran_lcr', comm)
return
end if
allocate (hLC_cmp(transport%num_ll, transport%num_lc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hLC_cmp in tran_lcr', comm)
return
end if
allocate (hCR_cmp(transport%num_cr, transport%num_rr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hCR_cmp in tran_lcr', comm)
return
end if
!If construct used only when reading matrices from file
if (transport%read_ht) then
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', 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', comm)
return
end if
allocate (hC(transport%num_cc, transport%num_cc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hC in tran_lcr', comm)
return
end if
allocate (hLC(transport%num_ll, transport%num_lc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hLC in tran_lcr', comm)
return
end if
allocate (hCR(transport%num_cr, transport%num_rr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hCR in tran_lcr', comm)
return
end if
filename = trim(seedname)//'_htL.dat'
call tran_read_htX(transport%num_ll, hL0, hL1, filename, stdout, error, comm)
if (allocated(error)) return
if (.not. transport%use_same_lead) then
allocate (hR0(transport%num_rr, transport%num_rr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hR0 in tran_lcr', comm)
return
end if
allocate (hR1(transport%num_rr, transport%num_rr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating hR1 in tran_lcr', comm)
return
end if
filename = trim(seedname)//'_htR.dat'
call tran_read_htX(transport%num_rr, hR0, hR1, filename, stdout, error, comm)
if (allocated(error)) return
end if
filename = trim(seedname)//'_htC.dat'
call tran_read_htC(transport%num_cc, hC, filename, stdout, error, comm)
if (allocated(error)) return
filename = trim(seedname)//'_htLC.dat'
call tran_read_htXY(transport%num_ll, transport%num_lc, hLC, filename, stdout, error, comm)
if (allocated(error)) return
filename = trim(seedname)//'_htCR.dat'
call tran_read_htXY(transport%num_cr, transport%num_rr, hCR, filename, stdout, error, comm)
if (allocated(error)) return
end if
! Banded matrix H_C : save memory !
do j = 1, transport%num_cc
do i = max(1, j - KU), min(transport%num_cc, j + KL)
hCband(KL + KU + 1 + i - j, j) = hC(i, j)
end do
end do
deallocate (hC, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating hC in tran_lcr', comm)
return
end if
! H_LC : to a complex matrix
hLC_cmp(:, :) = cmplx(hLC(:, :), kind=dp)
deallocate (hLC, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating hLC in tran_lcr', comm)
return
end if
! H_CR : to a complex matrix
hCR_cmp(:, :) = cmplx(hCR(:, :), kind=dp)
deallocate (hCR, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating hCR in tran_lcr', comm)
return
end if
allocate (totL(transport%num_ll, transport%num_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating totL in tran_lcr', comm)
return
end if
allocate (tottL(transport%num_ll, transport%num_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating tottL in tran_lcr', comm)
return
end if
if (.not. transport%use_same_lead) then
allocate (totR(transport%num_rr, transport%num_rr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating totR in tran_lcr', comm)
return
end if
allocate (tottR(transport%num_rr, transport%num_rr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating tottR in tran_lcr', comm)
return
end if
end if
allocate (g_surf_L(transport%num_ll, transport%num_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating g_surf_L in tran_lcr', comm)
return
end if
allocate (g_surf_R(transport%num_rr, transport%num_rr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating g_surf_R in tran_lcr', comm)
return
end if
allocate (g_C_inv(2*KL + KU + 1, transport%num_cc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating g_C_inv in tran_lcr', comm)
return
end if
allocate (g_C(transport%num_cc, transport%num_cc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating g_C in tran_lcr', comm)
return
end if
allocate (sLr(transport%num_lc, transport%num_lc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating sLr in tran_lcr', comm)
return
end if
allocate (sRr(transport%num_cr, transport%num_cr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating sRr in tran_lcr', comm)
return
end if
allocate (gL(transport%num_lc, transport%num_lc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating gL in tran_lcr', comm)
return
end if
allocate (gR(transport%num_cr, transport%num_cr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating gR in tran_lcr', comm)
return
end if
allocate (c1(transport%num_lc, transport%num_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating c1 in tran_lcr', comm)
return
end if
allocate (c2(transport%num_cr, transport%num_rr), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating c2 in tran_lcr', comm)
return
end if
allocate (s1(KC, KC), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating s1 in tran_lcr', comm)
return
end if
allocate (s2(KC, KC), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating s2 in tran_lcr', comm)
return
end if
allocate (ipiv(transport%num_cc), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating ipiv in tran_lcr', comm)
return
end if
! Loop over the energies
n_e = floor((transport%win_max - transport%win_min)/transport%energy_step) + 1
write (stdout, '(/1x,a)', advance='no') 'Calculating quantum conductance and &
&density of states...'
do n = 1, n_e
e_scan = transport%win_min + real(n - 1, dp)*transport%energy_step
! compute conductance according to Fisher and Lee
! compute self-energies following Datta
e_scan_cmp = e_scan + eta
! Surface green function for the left lead : g_surf_L
call tran_transfer(totL, tottL, hL0, hL1, e_scan_cmp, transport%num_ll, stdout, error, comm)
if (allocated(error)) return
call tran_green(totL, tottL, hL0, hL1, e_scan, g_surf_L, -1, 1, transport%num_ll, stdout, &
error, comm)
if (allocated(error)) return
! Self-energy (Sigma_L) : sLr = (hLC_cmp)^+ * g_surf_L * hLC_cmp
c1 = cmplx_0
sLr = cmplx_0
call ZGEMM('C', 'N', transport%num_lc, transport%num_ll, transport%num_ll, cmplx_1, &
hLC_cmp, transport%num_ll, g_surf_L, transport%num_ll, cmplx_0, c1, transport%num_lc)
call ZGEMM('N', 'N', transport%num_lc, transport%num_lc, transport%num_ll, cmplx_1, &
c1, transport%num_lc, hLC_cmp, transport%num_ll, cmplx_0, sLr, transport%num_lc)
! Surface green function for the right lead : g_surf_R
if (transport%use_same_lead) then
call tran_green(totL, tottL, hL0, hL1, e_scan, g_surf_R, 1, 1, transport%num_rr, stdout, &
error, comm)
if (allocated(error)) return
else
call tran_transfer(totR, tottR, hR0, hR1, e_scan_cmp, transport%num_rr, stdout, error, comm)
if (allocated(error)) return
call tran_green(totR, tottR, hR0, hR1, e_scan, g_surf_R, 1, 1, transport%num_rr, stdout, &
error, comm)
if (allocated(error)) return
end if
! Self-energy (Sigma_R) : sRr = hCR_cmp * g_surf_R * (hCR_cmp)^+
c2 = cmplx_0
sRr = cmplx_0
call ZGEMM('N', 'N', transport%num_cr, transport%num_rr, transport%num_rr, cmplx_1, &
hCR_cmp, transport%num_cr, g_surf_R, transport%num_rr, cmplx_0, c2, transport%num_cr)
call ZGEMM('N', 'C', transport%num_cr, transport%num_cr, transport%num_rr, cmplx_1, &
c2, transport%num_cr, hCR_cmp, transport%num_cr, cmplx_0, sRr, transport%num_cr)
! g_C^-1 = -H
g_C_inv(:, :) = cmplx(-hCband(:, :), kind=dp)
! g_C^-1 = -H - Sigma_L^r
do j = 1, transport%num_lc
do i = max(1, j - KU), min(transport%num_lc, j + KL)
g_C_inv(KL + KU + 1 + i - j, j) = g_C_inv(KL + KU + 1 + i - j, j) - sLr(i, j)
end do
end do
! g_C^-1 = -H - Sigma_L^r - Sigma_R^r
do j = (transport%num_cc - transport%num_cr) + 1, transport%num_cc
do i = max((transport%num_cc - transport%num_cr) + 1, j - (transport%num_cr - 1)), &
min(transport%num_cc, j + (transport%num_cr - 1))
g_C_inv(KL + KU + 1 + i - j, j) = &
g_C_inv(KL + KU + 1 + i - j, j) - &
sRr(i - (transport%num_cc - transport%num_cr), j - (transport%num_cc - transport%num_cr))
end do
end do
! g_C^-1 = eI - H - Sigma_L^r - Sigma_R^r
do i = 1, transport%num_cc
g_C_inv(KL + KU + 1, i) = e_scan + g_C_inv(KL + KU + 1, i)
end do
! invert g_C^-1 => g_C
g_C = cmplx_0
do i = 1, transport%num_cc
g_C(i, i) = cmplx_1
end do
call zgbsv(transport%num_cc, KL, KU, transport%num_cc, g_C_inv, 2*KL + KU + 1, ipiv, g_C, &
transport%num_cc, info)
if (info .ne. 0) then
write (stdout, *) 'ERROR: IN ZGBSV IN tran_lcr, INFO=', info
call set_error_fatal(error, 'tran_lcr: problem in ZGBSV', comm)
return
end if
! Gamma_L = i(Sigma_L^r-Sigma_L^a)
gL = cmplx_i*(sLr - conjg(transpose(sLr)))
! s1 = Gamma_L * g_C^r
s1 = cmplx_0
do j = 1, KC
do i = 1, transport%num_lc
do k = 1, transport%num_lc
s1(i, j) = s1(i, j) + gL(i, k)*g_C(k, j + (transport%num_cc - KC))
end do
end do
end do
! Gamma_R = i(Sigma_R^r-Sigma_R^a)
gR = cmplx_i*(sRr - conjg(transpose(sRr)))
! s2 = Gamma_R * g_C^a
s2 = cmplx_0
do j = 1, KC
do i = 1, transport%num_cr
do k = 1, transport%num_cr
s2(i + (KC - transport%num_cr), j) = s2(i + (KC - transport%num_cr), j) &
+ gR(i, k)*conjg(g_C(j, k &
+ (transport%num_cc - transport%num_cr)))
end do
end do
end do
qc = 0.0_dp
do i = 1, KC
do j = 1, KC
qc = qc + real(s1(i, j)*s2(j, i), dp)
end do
end do
write (qc_unit, '(f15.9,f18.9)') e_scan, qc
! compute density of states for the conductor layer
dos = 0.0_dp
do i = 1, transport%num_cc
dos = dos - aimag(g_C(i, i))
end do
dos = dos/pi
write (dos_unit, '(f15.9,f18.9)') e_scan, dos
end do
write (stdout, '(a)') ' done'
close (qc_unit)
close (dos_unit)
deallocate (ipiv, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating ipiv in tran_lcr', comm)
return
end if
deallocate (s2, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating s2 in tran_lcr', comm)
return
end if
deallocate (s1, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating s1 in tran_lcr', comm)
return
end if
deallocate (c2, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating c2 in tran_lcr', comm)
return
end if
deallocate (c1, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating c1 in tran_lcr', comm)
return
end if
deallocate (gR, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating gR in tran_lcr', comm)
return
end if
deallocate (gL, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating gL in tran_lcr', comm)
return
end if
deallocate (sRr, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating sRr in tran_lcr', comm)
return
end if
deallocate (sLr, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating sLr in tran_lcr', comm)
return
end if
deallocate (g_C, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating g_C in tran_lcr', comm)
return
end if
deallocate (g_C_inv, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating g_C_inv in tran_lcr', comm)
return
end if
deallocate (g_surf_R, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating g_surf_R in tran_lcr', comm)
return
end if
deallocate (g_surf_L, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating g_surf_L in tran_lcr', comm)
return
end if
if (allocated(tottR)) then
deallocate (tottR, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating tottR in tran_lcr', comm)
return
end if
end if
if (allocated(totR)) then
deallocate (totR, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating totR in tran_lcr', comm)
return
end if
end if
deallocate (tottL, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating tottL in tran_lcr', comm)
return
end if
deallocate (totL, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating totL in tran_lcr', comm)
return
end if
deallocate (hCR_cmp, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating hCR_cmp in tran_lcr', comm)
return
end if
deallocate (hLC_cmp, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating hLC_cmp in tran_lcr', comm)
return
end if
deallocate (hCband, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating hCband in tran_lcr', comm)
return
end if
if (timing_level > 1) call io_stopwatch_stop('tran: lcr', timer)
return
end subroutine tran_lcr
!================================================!
subroutine tran_transfer(tot, tott, h_00, h_01, e_scan_cmp, nxx, stdout, error, comm)
!================================================!
!
! iterative construction of the transfer matrix
! as Lopez-Sancho^2&Rubio, J.Phys.F:Met.Phys., v.14, 1205 (1984)
! and ibid. v.15, 851 (1985)
!
!================================================
use w90_constants, only: dp, cmplx_0, cmplx_1, eps7
use w90_error, only: w90_error_type, set_error_alloc, set_error_fatal, set_error_unconv, &
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
complex(kind=dp), intent(in) :: e_scan_cmp
complex(kind=dp), intent(out) :: tot(nxx, nxx)
complex(kind=dp), intent(out) :: tott(nxx, nxx)
real(kind=dp), intent(in) :: h_00(nxx, nxx)
real(kind=dp), intent(in) :: h_01(nxx, nxx)
integer :: ierr, info
integer :: i, j, n, nxx2
integer, allocatable :: ipiv(:)
real(kind=dp) :: conver, conver2
complex(kind=dp), allocatable, dimension(:, :) :: tsum, tsumt
complex(kind=dp), allocatable, dimension(:, :) :: t11, t12
complex(kind=dp), allocatable, dimension(:, :) :: s1, s2
complex(kind=dp), allocatable, dimension(:, :, :) :: tau, taut
allocate (ipiv(nxx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating ipiv in tran_transfer', comm)
return
end if
allocate (tsum(nxx, nxx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating tsum in tran_transfer', comm)
return
end if
allocate (tsumt(nxx, nxx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating tsumt in tran_transfer', comm)
return
end if
allocate (t11(nxx, nxx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating t11 in tran_transfer', comm)
return
end if
allocate (t12(nxx, nxx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating t12 in tran_transfer', 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_transfer', 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_transfer', comm)
return
end if
allocate (tau(nxx, nxx, 2), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating tau in tran_transfer', comm)
return
end if
allocate (taut(nxx, nxx, 2), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating taut in tran_transfer', comm)
return
end if
nxx2 = nxx*nxx
tot = cmplx_0
tott = cmplx_0
! construction of the transfer matrix
! t12 = e - h_00
t12(:, :) = cmplx(-h_00(:, :), kind=dp)
do i = 1, nxx
t12(i, i) = e_scan_cmp + t12(i, i)
end do
! compute (e - h_00)^-1 and store it in t11
t11 = cmplx_0
do i = 1, nxx
t11(i, i) = cmplx_1
end do
! inverse of t12 -> t11
call ZGESV(nxx, nxx, t12, nxx, ipiv, t11, nxx, info)
if (info .ne. 0) then
write (stdout, *) 'ERROR: IN ZGESV IN tran_transfer, INFO=', info
call set_error_fatal(error, 'tran_transfer: problem in ZGESV 1', comm)
return
end if
! compute intermediate t-matrices (defined as tau(nxx,nxx,niter)
! and taut(...)):
tau = cmplx_0
taut = cmplx_0
! t_0:
t12(:, :) = cmplx(h_01(:, :), kind=dp)
! tau = ( e - H_00 )^-1 * H_01^+
call ZGEMM('N', 'C', nxx, nxx, nxx, cmplx_1, t11, nxx, t12, nxx, cmplx_0, tau(1, 1, 1), nxx)
! taut = ( e - H_00 )^-1 * H_01
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, t11, nxx, t12, nxx, cmplx_0, taut(1, 1, 1), nxx)
! initialize T:
tot(:, :) = tau(:, :, 1)
tsum(:, :) = taut(:, :, 1)
! initialize T^bar:
tott(:, :) = taut(:, :, 1)
tsumt(:, :) = tau(:, :, 1)
! main loop:
do n = 1, nterx
t11 = cmplx_0
t12 = cmplx_0
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, tau(1, 1, 1), nxx, taut(1, 1, 1), nxx, cmplx_0, &
t11, nxx)
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, taut(1, 1, 1), nxx, tau(1, 1, 1), nxx, cmplx_0, &
t12, nxx)
s1(:, :) = -t11(:, :) - t12(:, :)
do i = 1, nxx
s1(i, i) = cmplx_1 + s1(i, i)
end do
s2 = cmplx_0
do i = 1, nxx
s2(i, i) = cmplx_1
end do
call ZGESV(nxx, nxx, s1, nxx, ipiv, s2, nxx, info)
if (info .ne. 0) then
write (stdout, *) 'ERROR: IN ZGESV IN tran_transfer, INFO=', info
call set_error_fatal(error, 'tran_transfer: problem in ZGESV 2', comm)
return
end if
t11 = cmplx_0
t12 = cmplx_0
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, tau(1, 1, 1), nxx, tau(1, 1, 1), nxx, cmplx_0, &
t11, nxx)
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, taut(1, 1, 1), nxx, taut(1, 1, 1), nxx, &
cmplx_0, t12, nxx)
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, s2, nxx, t11, nxx, cmplx_0, tau(1, 1, 2), nxx)
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, s2, nxx, t12, nxx, cmplx_0, taut(1, 1, 2), nxx)
! put the transfer matrices together
t11 = cmplx_0
s1 = cmplx_0
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, tsum, nxx, tau(1, 1, 2), nxx, cmplx_0, t11, nxx)
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, tsum, nxx, taut(1, 1, 2), nxx, cmplx_0, s1, nxx)
call ZCOPY(nxx2, t11, 1, s2, 1)
call ZAXPY(nxx2, cmplx_1, tot, 1, s2, 1)
tot(:, :) = s2(:, :)
tsum(:, :) = s1(:, :)
t11 = cmplx_0
s1 = cmplx_0
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, tsumt, nxx, taut(1, 1, 2), nxx, cmplx_0, t11, &
nxx)
call ZGEMM('N', 'N', nxx, nxx, nxx, cmplx_1, tsumt, nxx, tau(1, 1, 2), nxx, cmplx_0, s1, nxx)
call ZCOPY(nxx2, t11, 1, s2, 1)
call ZAXPY(nxx2, cmplx_1, tott, 1, s2, 1)
tott(:, :) = s2(:, :)
tsumt(:, :) = s1(:, :)
tau(:, :, 1) = tau(:, :, 2)
taut(:, :, 1) = taut(:, :, 2)
! convergence check on the t-matrices
conver = 0.0_dp
conver2 = 0.0_dp
do j = 1, nxx
do i = 1, nxx
conver = conver + sqrt(real(tau(i, j, 2), dp)**2 + aimag(tau(i, j, 2))**2)
conver2 = conver2 + sqrt(real(taut(i, j, 2), dp)**2 + aimag(taut(i, j, 2))**2)
end do
end do
if (conver .lt. eps7 .and. conver2 .lt. eps7) return
end do
if (conver .gt. eps7 .or. conver2 .gt. eps7) then
call set_error_unconv(error, 'Error in converging transfer matrix in tran_transfer', comm)
return
end if
deallocate (taut, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating taut in tran_transfer', comm)
return
end if
deallocate (tau, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating tau in tran_transfer', comm)
return
end if
deallocate (s2, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating s2 in tran_transfer', comm)
return
end if
deallocate (s1, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating s1 in tran_transfer', comm)
return
end if
deallocate (t12, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating t12 in tran_transfer', comm)
return
end if
deallocate (t11, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating t11 in tran_transfer', comm)
return
end if
deallocate (tsumt, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating tsumt in tran_transfer', comm)
return
end if
deallocate (tsum, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating tsum in tran_transfer', comm)
return
end if
deallocate (ipiv, stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in deallocating ipiv in tran_transfer', comm)
return
end if
return
end subroutine tran_transfer
!================================================!
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
!================================================!
subroutine tran_read_htX(nxx, h_00, h_01, h_file, stdout, error, comm)
!================================================!
use w90_constants, only: dp, maxlen
use w90_error, only: w90_error_type, set_error_file, set_error_file
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
real(kind=dp), intent(out) :: h_00(nxx, nxx), h_01(nxx, nxx)
character(len=50), intent(in) :: h_file
integer :: i, j, nw, file_unit
character(len=maxlen) :: dummy
open (newunit=file_unit, file=h_file, form='formatted', &
status='old', action='read', err=101)
write (stdout, '(/a)', advance='no') ' Reading H matrix from '//h_file//' : '
read (file_unit, '(a)', err=102, end=102) dummy
write (stdout, '(a)') trim(dummy)
read (file_unit, *, err=102, end=102) nw
if (nw .ne. nxx) then
call set_error_file(error, 'wrong matrix size in transport: read_htX', comm)
return
end if
read (file_unit, *) ((h_00(i, j), i=1, nxx), j=1, nxx)
read (file_unit, *, err=102, end=102) nw
if (nw .ne. nxx) then
call set_error_file(error, 'wrong matrix size in transport: read_htX', comm)
return
end if
read (file_unit, *, err=102, end=102) ((h_01(i, j), i=1, nxx), j=1, nxx)
close (unit=file_unit)
return
101 call set_error_file(error, 'Error: Problem opening input file '//h_file, comm)
return
102 call set_error_file(error, 'Error: Problem reading input file '//h_file, comm)
return
end subroutine tran_read_htX
!================================================!
subroutine tran_read_htC(nxx, h_00, h_file, stdout, error, comm)
!================================================!
use w90_constants, only: dp, maxlen
use w90_error, only: w90_error_type, set_error_file, set_error_file
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
real(kind=dp), intent(out) :: h_00(nxx, nxx)
character(len=50), intent(in) :: h_file
integer :: i, j, nw, file_unit, ierr
character(len=maxlen) :: dummy
open (newunit=file_unit, file=h_file, form='formatted', status='old', action='read', iostat=ierr)
if (ierr /= 0) then
call set_error_file(error, 'Error: Problem opening input file '//h_file, comm)
return
end if
write (stdout, '(/a)', advance='no') ' Reading H matrix from '//h_file//' : '
read (file_unit, '(a)', err=102, end=102) dummy
write (stdout, '(a)') trim(dummy)
read (file_unit, *, err=102, end=102) nw
if (nw .ne. nxx) then
call set_error_file(error, 'wrong matrix size in transport: read_htC', comm)
return
end if
read (file_unit, *, err=102, end=102) ((h_00(i, j), i=1, nxx), j=1, nxx)
close (unit=file_unit)
return
102 call set_error_file(error, 'Error: Problem reading input file '//h_file, comm)
end subroutine tran_read_htC
!================================================!
subroutine tran_read_htXY(nxx1, nxx2, h_01, h_file, stdout, error, comm)
!================================================!
use w90_constants, only: dp, maxlen
use w90_error, only: w90_error_type, set_error_file, set_error_file
implicit none
type(w90_error_type), allocatable, intent(out) :: error
type(w90_comm_type), intent(in) :: comm
integer, intent(in) :: nxx1, nxx2
integer, intent(in) :: stdout
real(kind=dp), intent(out) :: h_01(nxx1, nxx2)
character(len=50), intent(in) :: h_file
integer :: i, j, nw1, nw2, file_unit
character(len=maxlen) :: dummy
open (newunit=file_unit, file=h_file, form='formatted', &
status='old', action='read', err=101)
write (stdout, '(/a)', advance='no') ' Reading H matrix from '//h_file//' : '
read (file_unit, '(a)', err=102, end=102) dummy
write (stdout, '(a)') trim(dummy)
read (file_unit, *, err=102, end=102) nw1, nw2
if (nw1 .ne. nxx1 .or. nw2 .ne. nxx2) then
call set_error_file(error, 'wrong matrix size in transport: read_htXY', comm)
return
end if
read (file_unit, *, err=102, end=102) ((h_01(i, j), i=1, nxx1), j=1, nxx2)
close (unit=file_unit)
return
101 call set_error_file(error, 'Error: Problem opening input file '//h_file, comm)
return
102 call set_error_file(error, 'Error: Problem reading input file '//h_file, comm)
return
end subroutine tran_read_htXY
!================================================
subroutine tran_find_integral_signatures(signatures, num_G, print_output, real_lattice, &
u_matrix_opt, u_matrix, num_bands, num_wann, &
have_disentangled, wannier_centres_translated, stdout, &
seedname, timer, error, comm)
!================================================!
! Reads <seedname>.unkg file that contains the u_nk(G) and calculate
! Fourier components of each wannier function. Linear combinations of
! these provide integral of different spatial dependence.
! The set of these integrals provide a signature for distinguishing the
! type and 'parity' of each wannier function.
!================================================!
use w90_constants, only: dp, cmplx_0, twopi, cmplx_i
use w90_io, only: io_date, io_stopwatch_start, io_stopwatch_stop
use w90_types, only: print_output_type, timer_list_type
use w90_error, only: w90_error_type, set_error_alloc, set_error_file, set_error_file, &
set_error_dealloc
implicit none
type(print_output_type), intent(in) :: print_output
type(timer_list_type), intent(inout) :: timer
type(w90_error_type), allocatable, intent(out) :: error
type(w90_comm_type), intent(in) :: comm
integer, intent(in) :: num_bands
integer, intent(in) :: num_wann
integer, intent(in) :: stdout
integer, intent(out) :: num_G
complex(kind=dp), intent(in) :: u_matrix(:, :, :)
complex(kind=dp), intent(in) :: u_matrix_opt(:, :, :)
real(kind=dp), allocatable, intent(out) :: signatures(:, :)
real(kind=dp), intent(in) :: real_lattice(3, 3)
real(kind=dp), intent(in) :: wannier_centres_translated(:, :)
character(len=50), intent(in) :: seedname
logical, intent(in) :: have_disentangled
integer, allocatable :: g_abc(:, :)
integer :: i, ibnd, file_unit, ierr, p, p_max, n, m, ig, a, b, c, ig_idx(32)
real(kind=dp) :: i_unkg, r_unkg, wf_frac(3), det_rl, inv_t_rl(3, 3), mag_signature_sq
complex(kind=dp), allocatable :: unkg(:, :), tran_u_matrix(:, :)
complex(kind=dp) :: phase_factor, signature_basis(32)
logical :: have_file
if (print_output%timing_level > 1) call io_stopwatch_start('tran: find_sigs_unkg_int', timer)
inquire (file=trim(seedname)//'.unkg', exist=have_file)
if (.not. have_file) then
call set_error_file(error, 'tran_hr_parity_unkg: file '//trim(seedname)// &
'.unkg not found', comm)
return
end if
open (newunit=file_unit, file=trim(seedname)//'.unkg', form='formatted', action='read')
write (stdout, '(3a)') ' Reading '//trim(seedname)//'.unkg file'
read (file_unit, *) num_G
allocate (signatures(20, num_bands), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating signatures in tran_find_integral_signatures', comm)
return
end if
allocate (unkg(num_G, num_bands), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating unkg in tran_find_integral_signatures', comm)
return
end if
allocate (g_abc(num_G, 3), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating g_abc in tran_find_integral_signatures', comm)
return
end if
if (have_disentangled) then
allocate (tran_u_matrix(num_bands, num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating tran_u_matrix in tran_find_integral_signatures', comm)
return
end if
else
allocate (tran_u_matrix(num_wann, num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating tran_u_matrix in tran_find_integral_signatures', comm)
return
end if
end if
unkg = cmplx_0
do m = 1, num_bands
do i = 1, num_G
read (file_unit, *) ibnd, ig, a, b, c, r_unkg, i_unkg
if ((ig .ne. i) .OR. (ibnd .ne. m)) then
call set_error_file(error, 'tran_find_integral_signatures: Incorrect bands or g vectors', comm)
return
end if
unkg(i, m) = cmplx(r_unkg, i_unkg, kind=dp)
g_abc(i, :) = (/a, b, c/)
end do
end do
close (file_unit)
! Computing inverse of transpose of real_lattice
det_rl = real_lattice(1, 1)*(real_lattice(2, 2)*real_lattice(3, 3) &
- real_lattice(2, 3)*real_lattice(3, 2)) &
- real_lattice(2, 2)*(real_lattice(2, 1)*real_lattice(3, 3) &
- real_lattice(2, 3)*real_lattice(3, 1)) &
+ real_lattice(3, 3)*(real_lattice(2, 1)*real_lattice(3, 2) &
- real_lattice(2, 2)*real_lattice(3, 1))
inv_t_rl(1, 1) = (real_lattice(2, 2)*real_lattice(3, 3) - real_lattice(3, 2)*real_lattice(2, 3))
inv_t_rl(1, 2) = (real_lattice(2, 1)*real_lattice(3, 3) - real_lattice(3, 1)*real_lattice(2, 3))
inv_t_rl(1, 3) = (real_lattice(2, 1)*real_lattice(3, 2) - real_lattice(3, 1)*real_lattice(2, 2))
inv_t_rl(2, 1) = (real_lattice(1, 2)*real_lattice(3, 3) - real_lattice(3, 2)*real_lattice(1, 3))
inv_t_rl(2, 2) = (real_lattice(1, 1)*real_lattice(3, 3) - real_lattice(3, 1)*real_lattice(1, 3))
inv_t_rl(2, 3) = (real_lattice(1, 1)*real_lattice(3, 2) - real_lattice(3, 1)*real_lattice(1, 2))
inv_t_rl(3, 1) = (real_lattice(1, 2)*real_lattice(2, 3) - real_lattice(2, 2)*real_lattice(1, 3))
inv_t_rl(3, 2) = (real_lattice(1, 1)*real_lattice(2, 3) - real_lattice(2, 1)*real_lattice(1, 3))
inv_t_rl(3, 3) = (real_lattice(1, 1)*real_lattice(2, 2) - real_lattice(2, 1)*real_lattice(1, 2))
inv_t_rl = inv_t_rl/det_rl
!Loop over wannier functions to compute generalised U matrix
signatures = 0.0_dp
tran_u_matrix = cmplx_0
do n = 1, num_wann
!Disentanglement step
if (have_disentangled) then
do p = 1, num_bands
do m = 1, num_wann
tran_u_matrix(p, n) = tran_u_matrix(p, n) + u_matrix(m, n, 1)*u_matrix_opt(p, m, 1)
end do
end do
p_max = num_bands
else
tran_u_matrix = u_matrix(:, :, 1)
p_max = num_wann
end if
end do
if (print_output%iprint .ge. 5) write (stdout, *) 'Printing integral signatures for each &
&wannier function:'
! Loop over all wannier functions
do n = 1, num_wann
! Find fraction coordinate of wannier function in lattice vector basis
! wf_frac(:)=(transpose(real_lattice))^(-1) * wannier_centres_translated(:,n)
wf_frac = 0.0_dp
wf_frac = matmul(inv_t_rl, wannier_centres_translated(:, n))
! Loop over all g vectors, find a,b,c's required
do ig = 1, num_G
! 0th Order
if ((g_abc(ig, 1) .eq. 0) .and. (g_abc(ig, 2) .eq. 0) .and. (g_abc(ig, 3) .eq. 0)) ig_idx(1) = ig ! 1
! 1st Order
if ((g_abc(ig, 1) .eq. 1) .and. (g_abc(ig, 2) .eq. 0) .and. (g_abc(ig, 3) .eq. 0)) ig_idx(2) = ig ! x
if ((g_abc(ig, 1) .eq. 0) .and. (g_abc(ig, 2) .eq. 1) .and. (g_abc(ig, 3) .eq. 0)) ig_idx(3) = ig ! y
if ((g_abc(ig, 1) .eq. 0) .and. (g_abc(ig, 2) .eq. 0) .and. (g_abc(ig, 3) .eq. 1)) ig_idx(4) = ig ! z
! 2nd Order
if ((g_abc(ig, 1) .eq. 2) .and. (g_abc(ig, 2) .eq. 0) .and. (g_abc(ig, 3) .eq. 0)) ig_idx(5) = ig ! x^2
if ((g_abc(ig, 1) .eq. 1) .and. (g_abc(ig, 2) .eq. 1) .and. (g_abc(ig, 3) .eq. 0)) ig_idx(6) = ig ! xy
if ((g_abc(ig, 1) .eq. 1) .and. (g_abc(ig, 2) .eq. -1) .and. (g_abc(ig, 3) .eq. 0)) ig_idx(7) = ig ! xy
if ((g_abc(ig, 1) .eq. 1) .and. (g_abc(ig, 2) .eq. 0) .and. (g_abc(ig, 3) .eq. 1)) ig_idx(8) = ig ! xz
if ((g_abc(ig, 1) .eq. 1) .and. (g_abc(ig, 2) .eq. 0) .and. (g_abc(ig, 3) .eq. -1)) ig_idx(9) = ig ! xz
if ((g_abc(ig, 1) .eq. 0) .and. (g_abc(ig, 2) .eq. 2) .and. (g_abc(ig, 3) .eq. 0)) ig_idx(10) = ig ! y^2
if ((g_abc(ig, 1) .eq. 0) .and. (g_abc(ig, 2) .eq. 1) .and. (g_abc(ig, 3) .eq. 1)) ig_idx(11) = ig ! yz
if ((g_abc(ig, 1) .eq. 0) .and. (g_abc(ig, 2) .eq. 1) .and. (g_abc(ig, 3) .eq. -1)) ig_idx(12) = ig ! yz
if ((g_abc(ig, 1) .eq. 0) .and. (g_abc(ig, 2) .eq. 0) .and. (g_abc(ig, 3) .eq. 2)) ig_idx(13) = ig ! z^2
! 3rd Order
if ((g_abc(ig, 1) .eq. 3) .and. (g_abc(ig, 2) .eq. 0) .and. (g_abc(ig, 3) .eq. 0)) ig_idx(14) = ig ! x^3
if ((g_abc(ig, 1) .eq. 2) .and. (g_abc(ig, 2) .eq. 1) .and. (g_abc(ig, 3) .eq. 0)) ig_idx(15) = ig ! x^2y
if ((g_abc(ig, 1) .eq. 2) .and. (g_abc(ig, 2) .eq. -1) .and. (g_abc(ig, 3) .eq. 0)) ig_idx(16) = ig ! x^2y
if ((g_abc(ig, 1) .eq. 2) .and. (g_abc(ig, 2) .eq. 0) .and. (g_abc(ig, 3) .eq. 1)) ig_idx(17) = ig ! x^2z
if ((g_abc(ig, 1) .eq. 2) .and. (g_abc(ig, 2) .eq. 0) .and. (g_abc(ig, 3) .eq. -1)) ig_idx(18) = ig ! x^2z
if ((g_abc(ig, 1) .eq. 1) .and. (g_abc(ig, 2) .eq. 2) .and. (g_abc(ig, 3) .eq. 0)) ig_idx(19) = ig ! xy^2
if ((g_abc(ig, 1) .eq. 1) .and. (g_abc(ig, 2) .eq. -2) .and. (g_abc(ig, 3) .eq. 0)) ig_idx(20) = ig ! xy^2
if ((g_abc(ig, 1) .eq. 1) .and. (g_abc(ig, 2) .eq. 1) .and. (g_abc(ig, 3) .eq. 1)) ig_idx(21) = ig ! xyz
if ((g_abc(ig, 1) .eq. 1) .and. (g_abc(ig, 2) .eq. 1) .and. (g_abc(ig, 3) .eq. -1)) ig_idx(22) = ig ! xyz
if ((g_abc(ig, 1) .eq. 1) .and. (g_abc(ig, 2) .eq. -1) .and. (g_abc(ig, 3) .eq. 1)) ig_idx(23) = ig ! xyz
if ((g_abc(ig, 1) .eq. 1) .and. (g_abc(ig, 2) .eq. -1) .and. (g_abc(ig, 3) .eq. -1)) ig_idx(24) = ig ! xyz
if ((g_abc(ig, 1) .eq. 1) .and. (g_abc(ig, 2) .eq. 0) .and. (g_abc(ig, 3) .eq. 2)) ig_idx(25) = ig ! xz^2
if ((g_abc(ig, 1) .eq. 1) .and. (g_abc(ig, 2) .eq. 0) .and. (g_abc(ig, 3) .eq. -2)) ig_idx(26) = ig ! xz^2
if ((g_abc(ig, 1) .eq. 0) .and. (g_abc(ig, 2) .eq. 3) .and. (g_abc(ig, 3) .eq. 0)) ig_idx(27) = ig ! y^3
if ((g_abc(ig, 1) .eq. 0) .and. (g_abc(ig, 2) .eq. 2) .and. (g_abc(ig, 3) .eq. 1)) ig_idx(28) = ig ! y^2z
if ((g_abc(ig, 1) .eq. 0) .and. (g_abc(ig, 2) .eq. 2) .and. (g_abc(ig, 3) .eq. -1)) ig_idx(29) = ig ! y^2z
if ((g_abc(ig, 1) .eq. 0) .and. (g_abc(ig, 2) .eq. 1) .and. (g_abc(ig, 3) .eq. 2)) ig_idx(30) = ig ! yz^2
if ((g_abc(ig, 1) .eq. 0) .and. (g_abc(ig, 2) .eq. 1) .and. (g_abc(ig, 3) .eq. -2)) ig_idx(31) = ig ! yz^2
if ((g_abc(ig, 1) .eq. 0) .and. (g_abc(ig, 2) .eq. 0) .and. (g_abc(ig, 3) .eq. 3)) ig_idx(32) = ig ! z^3
end do
! Loop over the 32 required g-vectors
signature_basis = cmplx_0
do ig = 1, 32
phase_factor = cmplx_0
! Compute the phase factor exp(-i*G*x_c)
phase_factor = exp(-twopi*cmplx_i*(g_abc(ig_idx(ig), 1)*wf_frac(1) &
+ g_abc(ig_idx(ig), 2)*wf_frac(2) &
+ g_abc(ig_idx(ig), 3)*wf_frac(3)))
! Compute integrals that form the basis of the spatial integrals that form the signature
do p = 1, p_max
signature_basis(ig) = signature_basis(ig) + tran_u_matrix(p, n)*conjg(unkg(ig_idx(ig), p))
end do
signature_basis(ig) = signature_basis(ig)*phase_factor
end do
! Definitions of the signature integrals
! 0th Order
signatures(1, n) = real(signature_basis(1)) ! 1
! 1st Order
signatures(2, n) = aimag(signature_basis(2)) ! x
signatures(3, n) = aimag(signature_basis(3)) ! y
signatures(4, n) = aimag(signature_basis(4)) ! z
! 2nd Orde r
signatures(5, n) = real(signature_basis(1) - signature_basis(5))/2 ! x^2
signatures(6, n) = real(signature_basis(7) - signature_basis(6))/2 ! xy
signatures(7, n) = real(signature_basis(9) - signature_basis(8))/2 ! xz
signatures(8, n) = real(signature_basis(1) - signature_basis(10))/2 ! y^2
signatures(9, n) = real(signature_basis(12) - signature_basis(11))/2 ! yz
signatures(10, n) = real(signature_basis(1) - signature_basis(13))/2 ! z^2
! 3rd Order
signatures(11, n) = aimag(3*signature_basis(2) - signature_basis(14))/4 ! x^3
signatures(12, n) = aimag(2*signature_basis(3) + signature_basis(16) - signature_basis(15))/4 ! x^2y
signatures(13, n) = aimag(2*signature_basis(4) + signature_basis(18) - signature_basis(17))/4 ! x^2z
signatures(14, n) = aimag(2*signature_basis(2) - signature_basis(20) - signature_basis(19))/4 ! xy^2
signatures(15, n) = aimag(signature_basis(23) + signature_basis(22) - signature_basis(21) &
- signature_basis(24))/4 ! xyz
signatures(16, n) = aimag(2*signature_basis(2) - signature_basis(26) - signature_basis(25))/4 ! xz^2
signatures(17, n) = aimag(3*signature_basis(3) - signature_basis(27))/4 ! y^3
signatures(18, n) = aimag(2*signature_basis(4) + signature_basis(29) - signature_basis(28))/4 ! y^2z
signatures(19, n) = aimag(2*signature_basis(3) - signature_basis(31) - signature_basis(30))/4 ! yz^2
signatures(20, n) = aimag(3*signature_basis(4) - signature_basis(32))/4 ! z^3
if (print_output%iprint .ge. 5) then
write (stdout, *) ' '
write (stdout, *) ' Wannier function: ', n
do ig = 1, 20
write (stdout, *) ig - 1, signatures(ig, n)
end do
end if
!Normalise signature of each wannier function to a unit vector
mag_signature_sq = 0.0_dp
do ig = 1, 20
mag_signature_sq = mag_signature_sq + abs(signatures(ig, n))**2
end do
signatures(:, n) = signatures(:, n)/sqrt(mag_signature_sq)
end do ! Wannier Function loop
! Set num_G = 20 to ensure later subroutines work correctly
num_G = 20
deallocate (tran_u_matrix, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating tran_u_matrix in tran_find_integral_signatures', comm)
return
end if
deallocate (g_abc, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating g_abc in tran_find_integral_signatures', comm)
return
end if
deallocate (unkg, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating unkg in tran_find_integral_signatures', comm)
return
end if
if (print_output%timing_level > 1) call io_stopwatch_stop('tran: find_sigs_unkg_int', timer)
return
end subroutine tran_find_integral_signatures
!================================================!
subroutine tran_lcr_2c2_sort(signatures, num_G, pl_warning, transport, atom_data, wannier_data, &
real_space_ham, print_output, real_lattice, num_wann, mp_grid, &
ham_r, irvec, nrpts, wannier_centres_translated, one_dim_vec, &
nrpts_one_dim, num_pl, coord, tran_sorted_idx, hr_one_dim, &
irvec_max, write_xyz, stdout, seedname, timer, error, comm)
!================================================!
! This is the main subroutine controling the sorting
! for the 2c2 geometry. We first sort in the conduction
! direction, group, sort in 2nd direction, group and
! sort in 3rd direction. Rigourous checks are performed
! to ensure group and subgroup structure is consistent
! between principal layers (PLs), and unit cells. Once
! checks are passed we consider the possibility of
! multiple wannier functions are of similar centre, and
! sort those
!================================================!
use w90_constants, only: dp
use w90_io, only: io_stopwatch_start, io_stopwatch_stop
use w90_types, only: wannier_data_type, atom_data_type, 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_alloc, set_error_dealloc, set_error_fatal
implicit none
! arguments
integer, intent(in) :: irvec(:, :)
integer, intent(inout) :: irvec_max ! maybe modified
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: num_G
integer, intent(in) :: num_wann
integer, intent(inout), allocatable :: tran_sorted_idx(:)
integer, intent(inout) :: coord(3)
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
real(kind=dp), intent(in) :: signatures(:, :)
real(kind=dp), allocatable, intent(inout) :: hr_one_dim(:, :, :)
real(kind=dp), intent(inout) :: wannier_centres_translated(:, :)
real(kind=dp), intent(in) :: real_lattice(3, 3)
complex(kind=dp), intent(in) :: ham_r(:, :, :)
type(atom_data_type), intent(in) :: atom_data
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(wannier_data_type), intent(in) :: wannier_data
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) :: write_xyz
logical, intent(out) :: pl_warning
! local variables
real(kind=dp), dimension(2, num_wann) :: centres_non_sorted, centres_initial_sorted
real(kind=dp), dimension(2, transport%num_ll) :: PL1, PL2, PL3, PL4, PL
real(kind=dp), dimension(2, num_wann - (4*transport%num_ll)) :: central_region
real(kind=dp) :: reference_position, cell_length, distance, PL_max_val, PL_min_val
integer :: i, j, k, PL_selector, sort_iterator, sort_iterator2
integer :: ierr, temp_coord_2, temp_coord_3, n
integer :: num_wann_cell_ll, num_wf_group1, num_wf_last_group
integer, allocatable :: PL_groups(:), PL1_groups(:), PL2_groups(:), PL3_groups(:), &
PL4_groups(:), central_region_groups(:)
integer, allocatable :: PL_subgroup_info(:, :), PL1_subgroup_info(:, :), &
PL2_subgroup_info(:, :), PL3_subgroup_info(:, :), PL4_subgroup_info(:, :), &
central_subgroup_info(:, :), temp_subgroup(:, :)
character(30) :: fmt_1
allocate (tran_sorted_idx(num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating tran_sorted_idx in tran_lcr_2c2_sort', comm)
return
end if
num_wann_cell_ll = transport%num_ll/transport%num_cell_ll
if (print_output%timing_level > 1) call io_stopwatch_start('tran: lcr_2c2_sort', timer)
sort_iterator = 0
!Check translated centres have been found
if (size(wannier_centres_translated) .eq. 0) then
call set_error_fatal(error, 'Translated centres not known : required perform lcr transport, &
&try restart=plot', comm)
return
end if
!read one_dim_dir and creates an array (coord) that correspond to the
!conduction direction (coord(1)) and the two perpendicular directions
!(coord(2),coord(3)), such that a right-handed set is formed
if (real_space_ham%one_dim_dir .eq. 1) then
coord(1) = 1
coord(2) = 2
coord(3) = 3
elseif (real_space_ham%one_dim_dir .eq. 2) then
coord(1) = 2
coord(2) = 3
coord(3) = 1
elseif (real_space_ham%one_dim_dir .eq. 3) then
coord(1) = 3
coord(2) = 1
coord(3) = 2
end if
!Check
if (((real_lattice(coord(1), coord(2)) .ne. 0) .or. (real_lattice(coord(1), coord(3)) .ne. 0)) .or. &
((real_lattice(coord(2), coord(1)) .ne. 0) .or. (real_lattice(coord(3), coord(1)) .ne. 0))) then
call set_error_fatal(error, &
'Lattice vector in conduction direction must point along x,y or z &
& direction and be orthogonal to the remaining lattice vectors.', comm)
return
end if
!Check
if (num_wann .le. 4*transport%num_ll) then
call set_error_fatal(error, 'Principle layers are too big.', comm)
return
end if
100 continue
!Extract a 2d array of the wannier_indices and their coord(1) from wannier_centers_translated
do i = 1, num_wann
centres_non_sorted(1, i) = i
centres_non_sorted(2, i) = wannier_centres_translated(coord(1), i)
end do
write (stdout, '(/a)') ' Sorting WFs into principal layers'
!Initial sorting according to coord(1).
call sort(centres_non_sorted, centres_initial_sorted)
!Extract principal layers. WARNING: This extraction implies the structure of the supercell is
!2 principal layers of lead on the left and on the right of a central conductor.
PL1 = centres_initial_sorted(:, 1:transport%num_ll)
PL2 = centres_initial_sorted(:, transport%num_ll + 1:2*transport%num_ll)
PL3 = centres_initial_sorted(:, num_wann - (2*transport%num_ll - 1):num_wann - (transport%num_ll))
PL4 = centres_initial_sorted(:, num_wann - (transport%num_ll - 1):)
if (sort_iterator .eq. 1) then
temp_coord_2 = coord(2)
temp_coord_3 = coord(3)
coord(2) = temp_coord_3
coord(3) = temp_coord_2
end if
if (print_output%iprint .ge. 4) then
write (stdout, *) ' Group Breakdown of each principal layer'
end if
!Loop over principal layers
do i = 1, 4
!Creating a variable PL_selector which choose the appropriate PL
PL_selector = i
select case (PL_selector)
case (1)
PL = PL1
case (2)
PL = PL2
case (3)
PL = PL3
case (4)
PL = PL4
end select
!Grouping wannier functions with similar coord(1)
call group(PL, PL_groups, transport%group_threshold, error, comm)
if (allocated(error)) return
if (print_output%iprint .ge. 4) then
!Print group breakdown
write (fmt_1, '(i5)') size(PL_groups)
fmt_1 = adjustl(fmt_1)
fmt_1 = '(a3,i1,a1,i5,a2,'//trim(fmt_1)//'i4,a1)'
write (stdout, fmt_1) ' PL', i, ' ', size(PL_groups), ' (', (PL_groups(j), j=1, &
size(PL_groups)), ')'
end if
!Returns the sorted PL and informations on this PL
allocate (PL_subgroup_info(size(PL_groups), maxval(PL_groups)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating PL_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
call master_sort_and_group(PL, PL_groups, PL_subgroup_info, transport%group_threshold, &
print_output, wannier_centres_translated, coord, stdout, timer, &
error, comm)
if (allocated(error)) return
select case (PL_selector)
case (1)
allocate (PL1_groups(size(PL_groups)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating PL1_groups in tran_lcr_2c2_sort', comm)
return
end if
allocate (PL1_subgroup_info(size(PL_groups), maxval(PL_groups)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating PL1_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
PL1 = PL
PL1_groups = PL_groups
PL1_subgroup_info = PL_subgroup_info
deallocate (PL_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL1_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
case (2)
allocate (PL2_groups(size(PL_groups)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating PL2_groups in tran_lcr_2c2_sort', comm)
return
end if
allocate (PL2_subgroup_info(size(PL_groups), maxval(PL_groups)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating PL2_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
PL2 = PL
PL2_groups = PL_groups
PL2_subgroup_info = PL_subgroup_info
deallocate (PL_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL2_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
case (3)
allocate (PL3_groups(size(PL_groups)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating PL3_groups in tran_lcr_2c2_sort', comm)
return
end if
allocate (PL3_subgroup_info(size(PL_groups), maxval(PL_groups)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating PL3_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
PL3 = PL
PL3_groups = PL_groups
PL3_subgroup_info = PL_subgroup_info
deallocate (PL_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL3_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
case (4)
allocate (PL4_groups(size(PL_groups)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating PL4_groups in tran_lcr_2c2_sort', comm)
return
end if
allocate (PL4_subgroup_info(size(PL_groups), maxval(PL_groups)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating PL4_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
PL4 = PL
PL4_groups = PL_groups
PL4_subgroup_info = PL_subgroup_info
deallocate (PL_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL4_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
end select
deallocate (PL_groups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL_groups in tran_lcr_2c2_sort', comm)
return
end if
end do ! Principal layer loop
!Grouping and sorting of central conductor region
!Define central region
central_region = centres_initial_sorted(:, 2*transport%num_ll + 1:num_wann - (2*transport%num_ll))
!Group central region
call group(central_region, central_region_groups, transport%group_threshold, error, comm)
if (allocated(error)) return
!Print central region group breakdown
if (print_output%iprint .ge. 4) then
write (stdout, *) ' Group Breakdown of central region'
write (fmt_1, '(i5)') size(central_region_groups)
fmt_1 = adjustl(fmt_1)
fmt_1 = '(a5,i5,a2,'//trim(fmt_1)//'i4,a1)'
write (stdout, fmt_1) ' ', size(central_region_groups), ' (', &
(central_region_groups(j), j=1, size(central_region_groups)), ')'
end if
!Returns sorted central group region
allocate (central_subgroup_info(size(central_region_groups), maxval(central_region_groups)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating central_group_info in tran_lcr_2c2_sort', comm)
return
end if
call master_sort_and_group(central_region, central_region_groups, central_subgroup_info, &
transport%group_threshold, print_output, &
wannier_centres_translated, coord, stdout, timer, error, comm)
if (allocated(error)) return
deallocate (central_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating central_group_info in tran_lcr_2c2_sort', comm)
return
end if
write (stdout, *) ' '
!Build the sorted index array
tran_sorted_idx = nint(centres_initial_sorted(1, :))
tran_sorted_idx(1:transport%num_ll) = nint(PL1(1, :))
tran_sorted_idx(transport%num_ll + 1:2*transport%num_ll) = nint(PL2(1, :))
tran_sorted_idx(2*transport%num_ll + 1:num_wann - (2*transport%num_ll)) = nint(central_region(1, :))
tran_sorted_idx(num_wann - (2*transport%num_ll - 1):num_wann - (transport%num_ll)) = nint(PL3(1, :))
tran_sorted_idx(num_wann - (transport%num_ll - 1):) = nint(PL4(1, :))
sort_iterator = sort_iterator + 1
!Checks:
if ((size(PL1_groups) .ne. size(PL2_groups)) .or. &
(size(PL2_groups) .ne. size(PL3_groups)) .or. &
(size(PL3_groups) .ne. size(PL4_groups))) then
if (sort_iterator .ge. 2) then
if (write_xyz) call tran_write_xyz(atom_data, transport, wannier_centres_translated, &
tran_sorted_idx, num_wann, seedname, stdout)
call set_error_fatal(error, 'Sorting techniques exhausted:&
& Inconsistent number of groups among principal layers', comm)
return
end if
if (print_output%iprint > 0) write (stdout, *) 'Inconsistent number of groups among principal layers: restarting sorting...'
deallocate (PL1_groups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL1_groups in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL2_groups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL2_groups in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL3_groups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL3_groups in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL4_groups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL4_groups in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL1_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL1_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL2_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL2_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL3_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL3_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL4_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL4_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
goto 100
end if
do i = 1, size(PL1_groups)
if ((PL1_groups(i) .ne. PL2_groups(i)) .or. &
(PL2_groups(i) .ne. PL3_groups(i)) .or. &
(PL3_groups(i) .ne. PL4_groups(i))) then
if (sort_iterator .ge. 2) then
if (write_xyz) call tran_write_xyz(atom_data, transport, wannier_centres_translated, &
tran_sorted_idx, num_wann, seedname, stdout)
call set_error_fatal(error, &
'Sorting techniques exhausted: Inconsitent number of wannier function among &
& similar groups within principal layers', comm)
return
end if
write (stdout, *) 'Inconsitent number of wannier function among &
&similar groups within& principal layers: restarting sorting...'
deallocate (PL1_groups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL1_groups in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL2_groups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL2_groups in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL3_groups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL3_groups in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL4_groups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL4_groups in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL1_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL1_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL2_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL2_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL3_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL3_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL4_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL4_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
goto 100
end if
end do
!Now we check that the leftmost group and the rightmost group aren't
!supposed to be the same group
reference_position = wannier_centres_translated(coord(1), tran_sorted_idx(1))
cell_length = real_lattice(coord(1), coord(1))
sort_iterator2 = 1
do i = 1, transport%num_ll
distance = abs(abs(reference_position - wannier_centres_translated(coord(1), tran_sorted_idx(num_wann - i + 1))) &
- cell_length)
if (distance .lt. transport%group_threshold) then
wannier_centres_translated(coord(1), tran_sorted_idx(num_wann - i + 1)) = &
wannier_centres_translated(coord(1), tran_sorted_idx(num_wann - i + 1)) - cell_length
sort_iterator2 = sort_iterator2 + 1
end if
end do
if (sort_iterator2 .gt. 1) then
write (stdout, *) ' Grouping inconsistency found between first and last unit cells: '
write (stdout, *) ' suspect Wannier functions have been translated. '
write (stdout, *) ' '
write (stdout, *) ' Rebuilding Hamiltonian...'
write (stdout, *) ' '
deallocate (hr_one_dim, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating hr_one_dim in tran_lcr_2c2_sort', 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)
write (stdout, *) ' '
write (stdout, *) ' Restarting sorting...'
write (stdout, *) ' '
sort_iterator = sort_iterator - 1
deallocate (PL1_groups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL1_groups in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL2_groups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL2_groups in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL3_groups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL3_groups in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL4_groups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL4_groups in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL1_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL1_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL2_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL2_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL3_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL3_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL4_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL4_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
goto 100
end if
! if we reach this point, we don't have any left/right problems anymore. So we now
! check for inconsistencies in subgroups
allocate (temp_subgroup(size(PL1_subgroup_info, 1), size(PL1_subgroup_info, 2)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating tmp_subgroup in tran_lcr_2c2_sort', comm)
return
end if
do i = 2, 4
select case (i)
case (2)
temp_subgroup = PL1_subgroup_info - PL2_subgroup_info
case (3)
temp_subgroup = PL1_subgroup_info - PL3_subgroup_info
case (4)
temp_subgroup = PL1_subgroup_info - PL4_subgroup_info
end select
do j = 1, size(temp_subgroup, 1)
do k = 1, size(temp_subgroup, 2)
if (temp_subgroup(j, k) .ne. 0) then
if (sort_iterator .ge. 2) then
if (write_xyz) call tran_write_xyz(atom_data, transport, &
wannier_centres_translated, &
tran_sorted_idx, num_wann, seedname, &
stdout)
call set_error_fatal(error, &
'Sorting techniques exhausted: Inconsitent subgroup structures among principal layers', comm)
return
end if
if (print_output%iprint > 0) &
write (stdout, *) 'Inconsistent subgroup structure among principal layers: restarting sorting...'
deallocate (temp_subgroup, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating tmp_subgroup in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL1_groups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL1_groups in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL2_groups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL2_groups in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL3_groups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL3_groups in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL4_groups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL4_groups in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL1_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL1_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL2_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL2_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL3_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL3_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
deallocate (PL4_subgroup_info, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating PL4_subgroup_info in tran_lcr_2c2_sort', comm)
return
end if
goto 100
end if
end do
end do
end do
! At this point, every check has been cleared, and we need to use
! the parity signatures of the WFs for the possibility of equal centres
call check_and_sort_similar_centres(signatures, num_G, atom_data, transport, print_output, &
num_wann, wannier_centres_translated, coord, &
tran_sorted_idx, write_xyz, stdout, seedname, timer, error, comm)
if (allocated(error)) return
write (stdout, *) ' '
write (stdout, *) '------------------------- Sorted Wannier Centres -----------------------------'
write (stdout, *) 'Sorted index Unsorted index x y z Spread '
write (stdout, *) '==================================PL1========================================='
n = 0
do k = 1, 4
if (k .eq. 2) write (stdout, *) '==================================PL2========================================='
if (k .eq. 3) then
write (stdout, *) '===========================Central Region==================================='
do i = 1, num_wann - 4*transport%num_ll
n = n + 1
write (stdout, FMT='(2x,i6,10x,i6,6x,4F12.6)') n, tran_sorted_idx(n), &
wannier_centres_translated(1, tran_sorted_idx(n)), &
wannier_centres_translated(2, tran_sorted_idx(n)), &
wannier_centres_translated(3, tran_sorted_idx(n)), &
wannier_data%spreads(tran_sorted_idx(n))*print_output%lenconfac**2
end do
write (stdout, *) '==================================PL3========================================='
end if
if (k .eq. 4) write (stdout, *) '==================================PL4========================================='
do i = 1, transport%num_cell_ll
do j = 1, num_wann_cell_ll
n = n + 1
write (stdout, FMT='(2x,i6,10x,i6,6x,4F12.6)') n, tran_sorted_idx(n), &
wannier_centres_translated(1, tran_sorted_idx(n)), &
wannier_centres_translated(2, tran_sorted_idx(n)), &
wannier_centres_translated(3, tran_sorted_idx(n)), &
wannier_data%spreads(tran_sorted_idx(n))*print_output%lenconfac**2
end do
if (i .ne. transport%num_cell_ll) write (stdout, *) '---------------------&
&---------------------------------------------------------'
end do
end do
write (stdout, *) '=============================================================================='
write (stdout, *) ' '
! MS: Use sorting to assess whether dist_cutoff is suitable for correct PL cut
! by using limits of coord(1) values in 1st and last groups of PL1, & 1st group of PL2
pl_warning = .false.
num_wf_group1 = size(PL1_subgroup_info(1, :))
if (size(PL1_groups) .ge. 1) then
num_wf_last_group = size(PL1_subgroup_info(size(PL1_groups), :))
else
!Exception for 1 group in unit cell.
num_wf_last_group = num_wann_cell_ll
end if
PL_min_val = maxval(wannier_centres_translated(coord(1), tran_sorted_idx(transport%num_ll &
- num_wf_last_group + 1:transport%num_ll))) &
- minval(wannier_centres_translated(coord(1), tran_sorted_idx(1:num_wf_group1)))
PL_max_val = minval(wannier_centres_translated(coord(1), tran_sorted_idx(transport%num_ll + &
1:transport%num_ll + num_wf_group1))) &
- minval(wannier_centres_translated(coord(1), tran_sorted_idx(1:num_wf_group1)))
if ((real_space_ham%dist_cutoff .lt. PL_min_val) .or. (real_space_ham%dist_cutoff .gt. PL_max_val)) then
write (stdout, '(a)') ' WARNING: Expected dist_cutoff to be a PL length, I think this'
write (stdout, '(2(a,f10.6),a)') ' WARNING: is somewhere between ', PL_min_val, ' and ', PL_max_val, ' Ang'
pl_warning = .true.
end if
! End MS.
if (print_output%timing_level > 1) call io_stopwatch_stop('tran: lcr_2c2_sort', timer)
return
end subroutine tran_lcr_2c2_sort
!================================================!
subroutine master_sort_and_group(Array, Array_groups, subgroup_info, tran_group_threshold, &
print_output, wannier_centres_translated, coord, stdout, timer, &
error, comm)
!================================================!
! General sorting and grouping subroutine which takes Array,
! an ordered in conduction direction array of wannier function
! indexes and positions, and returns the ordered (and grouped)
! indexes and positions after considering the other two
! directions. Sub group info is also return for later checks.
!================================================!
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_error, only: w90_error_type, set_error_alloc, set_error_dealloc
implicit none
! arguments
type(print_output_type), intent(in) :: print_output
type(timer_list_type), intent(inout) :: timer
type(w90_error_type), allocatable, intent(out) :: error
type(w90_comm_type), intent(in) :: comm
integer, intent(in) :: stdout !, Array_size => size(Array_groups) so not needed
integer, intent(in) :: Array_groups(:), coord(3)
integer, intent(out), allocatable :: subgroup_info(:, :)
real(kind=dp), intent(inout) :: Array(:, :) !(2, Array_size)
real(kind=dp), intent(in) :: tran_group_threshold, wannier_centres_translated(:, :)
! local variables
integer :: i, j, k, Array_num_groups, increment, ierr, subgroup_increment, group_num_subgroups
integer, allocatable :: group_subgroups(:)
real(kind=dp), allocatable :: group_array(:, :), sorted_group_array(:, :)
real(kind=dp), allocatable :: subgroup_array(:, :), sorted_subgroup_array(:, :)
character(30) :: fmt_2
if (print_output%timing_level > 2) call io_stopwatch_start('tran: lcr_2c2_sort: master_sort', timer)
allocate (subgroup_info(size(Array_groups), maxval(Array_groups)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating subgroup_info in master_sort_and_group', comm)
return
end if
subgroup_info = 0
!Number of groups inside the principal layer
Array_num_groups = size(Array_groups)
!Convenient variable which will be amended later. Used to appropriately extract the group array from the Array
increment = 1
!Loop over groups inside Array
do j = 1, Array_num_groups
allocate (group_array(2, Array_groups(j)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating group_array in master_sort_and_group', comm)
return
end if
allocate (sorted_group_array(2, Array_groups(j)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating sorted_group_array in master_sort_and_group', comm)
return
end if
!Extract the group from the Array
group_array = Array(:, increment:increment + Array_groups(j) - 1)
!Updating group_array to contain coord(2)
do k = 1, Array_groups(j)
group_array(2, k) = wannier_centres_translated(coord(2), int(group_array(1, k)))
end do
call sort(group_array, sorted_group_array)
call group(sorted_group_array, group_subgroups, tran_group_threshold, error, comm)
if (allocated(error)) return
group_num_subgroups = size(group_subgroups)
if (print_output%iprint .ge. 4) then
!Printing subgroup breakdown
write (fmt_2, '(i5)') group_num_subgroups
fmt_2 = adjustl(fmt_2)
fmt_2 = '(a7,i3,a1,i5,a2,'//trim(fmt_2)//'i4,a1)'
write (stdout, fmt_2) ' Group ', j, ' ', group_num_subgroups, ' (', (group_subgroups(i), i=1, group_num_subgroups), ')'
end if
! filling up subgroup_info
do k = 1, group_num_subgroups
subgroup_info(j, k) = group_subgroups(k)
end do
!Convenient variable which will be amended later. Used to appropriately extract the subgroup array from the group_array
subgroup_increment = 1
!Loop over subgroups inside group
do k = 1, group_num_subgroups
allocate (subgroup_array(2, group_subgroups(k)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating subgroup_array in master_sort_and_group', comm)
return
end if
allocate (sorted_subgroup_array(2, group_subgroups(k)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating sorted_subgroup_array in master_sort_and_group', comm)
return
end if
!Extract the subgroup from the group
subgroup_array = sorted_group_array(:, subgroup_increment:subgroup_increment + group_subgroups(k) - 1)
!Updating subgroup_array to contain coord(3)
do i = 1, group_subgroups(k)
subgroup_array(2, i) = wannier_centres_translated(coord(3), int(subgroup_array(1, i)))
end do
call sort(subgroup_array, sorted_subgroup_array)
!Update sorted_group array with the sorted subgroup array
sorted_group_array(:, subgroup_increment:subgroup_increment + group_subgroups(k) - 1) = sorted_subgroup_array
!Update the subgroup_increment
subgroup_increment = subgroup_increment + group_subgroups(k)
deallocate (sorted_subgroup_array, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating sorted_subgroup_array in master_sort_and_group', comm)
return
end if
deallocate (subgroup_array, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating subgroup_array in master_sort_and_group', comm)
return
end if
end do
!Update Array with the sorted group array
Array(:, increment:increment + Array_groups(j) - 1) = sorted_group_array
!Update the group increment
increment = increment + Array_groups(j)
deallocate (group_array, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating group_array in master_sort_and_group', comm)
return
end if
deallocate (sorted_group_array, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating sorted_group_array in master_sort_and_group', comm)
return
end if
deallocate (group_subgroups, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating group_subgroups in master_sort_and_group', comm)
return
end if
end do
if (print_output%timing_level > 2) call io_stopwatch_stop('tran: lcr_2c2_sort: master_sort', timer)
return
end subroutine master_sort_and_group
!================================================!
subroutine sort(non_sorted, sorted)
!================================================!
use w90_constants, only: dp
implicit none
real(kind=dp), intent(inout) :: non_sorted(:, :)
real(kind=dp), intent(out) :: sorted(:, :)
integer, dimension(1) :: min_loc
integer :: num_col, i
num_col = size(non_sorted, 2)
do i = 1, num_col
!look for the location of the minimum value of the coordinates in non_sorted
min_loc = minloc(non_sorted(2, :))
!now the index in the first row of sorted is the index non_sorted(1,min_loc)
sorted(1, i) = non_sorted(1, min_loc(1))
!here is the corresponding coordinate
sorted(2, i) = non_sorted(2, min_loc(1))
!here one replaces the minimum coordinate with 10**10 such that this value
!will not be picked-up again by minloc
non_sorted(2, min_loc(1)) = 10.0**10
end do
return
end subroutine sort
!================================================!
subroutine group(array, array_groups, tran_group_threshold, error, comm)
!================================================!
use w90_constants, only: dp
use w90_error, only: w90_error_type, set_error_alloc, set_error_dealloc
implicit none
! arguments
type(w90_error_type), allocatable, intent(out) :: error
type(w90_comm_type), intent(in) :: comm
real(kind=dp), intent(in) :: tran_group_threshold
real(kind=dp), intent(in) :: array(:, :)
integer, intent(out), allocatable :: array_groups(:)
! local variables
integer, allocatable :: dummy_array(:)
integer :: array_idx, i, j, group_number, array_size, ierr
logical, allocatable :: logic(:)
array_size = size(array, 2)
allocate (dummy_array(array_size), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating dummy_array in group', comm)
return
end if
allocate (logic(array_size), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating logic in group', comm)
return
end if
!Initialise dummy array
dummy_array = 0
!Initialise logic to false
logic = .false.
!Define counter of number of groups
array_idx = 1
!Loop over columns of array (ie array_size)
do i = 1, array_size
!If an element of logic is true then it means the wannier function has already been grouped
if (logic(i) .eqv. .false.) then
!Create a group for the wannier function
logic(i) = .true.
!Initialise the number of wannier functions in this group to be 1
group_number = 1
!Loop over the rest of wannier functions in array
do j = min(i + 1, array_size), array_size
!Special termination cases
if ((j .eq. 1) .or. (i .eq. array_size)) then
dummy_array(array_idx) = group_number
exit
end if
if (j .eq. array_size .and. (abs(array(2, j) - array(2, i)) .le. tran_group_threshold)) then
group_number = group_number + 1
dummy_array(array_idx) = group_number
logic(j) = .true.
exit
end if
!Check distance between wannier function_i and wannier function_j
if (abs(array(2, j) - array(2, i)) .le. tran_group_threshold) then
!Increment number of wannier functions in group
group_number = group_number + 1
!Assigns wannier function to the group
logic(j) = .true.
else
!Group is finished and store number of wanniers in the group to dummy_array
dummy_array(array_idx) = group_number
!Increment number of groups
array_idx = array_idx + 1
exit
end if
end do
end if
end do
!Copy elements of dummy_array to array_groups
allocate (array_groups(array_idx), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating array_groups in group', comm)
return
end if
array_groups = dummy_array(:array_idx)
deallocate (dummy_array, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating dummy_array in group', comm)
return
end if
deallocate (logic, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating logic in group', comm)
return
end if
return
end subroutine group
!================================================
subroutine check_and_sort_similar_centres(signatures, num_G, atom_data, transport, print_output, &
num_wann, wannier_centres_translated, coord, &
tran_sorted_idx, write_xyz, stdout, seedname, timer, &
error, comm)
!================================================!
! Here, we consider the possiblity of wannier functions
! with similar centres, such as a set of d-orbitals
! centred on an atom. We use tran_group_threshold to
! identify them, if they exist, then use the signatures
! to dishinguish and sort then consistently from unit
! cell to unit cell.
!
! MS: For 2two-shot and beyond, some parameters,
! eg, first_group_element will need changing to consider
! the geometry of the new systems.
!================================================!
use w90_constants, only: dp
use w90_io, only: io_stopwatch_start, io_stopwatch_stop
use w90_types, only: atom_data_type, print_output_type, timer_list_type
use w90_wannier90_types, only: transport_type
use w90_error, only: w90_error_type, set_error_alloc, set_error_dealloc, set_error_fatal
implicit none
! arguments
type(atom_data_type), intent(in) :: atom_data
type(print_output_type), intent(in) :: print_output
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
integer, intent(in) :: coord(3)
integer, intent(in) :: num_G
integer, intent(in) :: num_wann
integer, intent(inout) :: tran_sorted_idx(:)
integer, intent(in) :: stdout
real(kind=dp), intent(in) :: wannier_centres_translated(:, :)
real(kind=dp), intent(in) :: signatures(:, :)
character(len=50), intent(in) :: seedname
logical, intent(in) :: write_xyz
! local variables
integer, allocatable :: idx_similar_wf(:), group_verifier(:), sorted_idx(:), centre_id(:)
integer, allocatable :: tmp_wf_verifier(:, :), wf_verifier(:, :), first_group_element(:, :)
integer, allocatable :: ref_similar_centres(:, :), unsorted_similar_centres(:, :)
integer, allocatable :: wf_similar_centres(:, :, :)
logical, allocatable :: has_similar_centres(:)
integer :: i, j, k, l, ierr, group_iterator, coord_iterator, num_wf_iterator, num_wann_cell_ll
integer :: iterator, max_position(1), p, num_wf_cell_iter
real(kind=dp), allocatable :: dot_p(:)
if (print_output%timing_level > 2) call io_stopwatch_start('tran: lcr_2c2_sort: similar_centres', timer)
num_wann_cell_ll = transport%num_ll/transport%num_cell_ll
allocate (wf_similar_centres(transport%num_cell_ll*4, num_wann_cell_ll, num_wann_cell_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating wf_similar_centre in check_and_sort_similar_centres', comm)
return
end if
allocate (idx_similar_wf(num_wann_cell_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating idx_similar_wf in check_and_sort_similar_centres', comm)
return
end if
allocate (has_similar_centres(num_wann_cell_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating has_similar_centres in check_and_sort_similar_centres', comm)
return
end if
allocate (tmp_wf_verifier(4*transport%num_cell_ll, num_wann_cell_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating tmp_wf_verifier in check_and_sort_similar_centres', comm)
return
end if
allocate (group_verifier(4*transport%num_cell_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating group_verifier in check_and_sort_similar_centres', comm)
return
end if
allocate (first_group_element(4*transport%num_cell_ll, num_wann_cell_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating first_group_element in check_and_sort_similar_centres', comm)
return
end if
allocate (centre_id(num_wann_cell_ll), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating centre_id in check_and_sort_similar_centres', comm)
return
end if
! First find WFs with similar centres: store in wf_similar_centres(cell#,group#,WF#)
group_verifier = 0
tmp_wf_verifier = 0
first_group_element = 0
centre_id = 0
! Loop over unit cells in PL1,PL2,PL3 and PL4
do i = 1, 4*transport%num_cell_ll
group_iterator = 0
has_similar_centres = .false.
! Loops over wannier functions in present unit cell
num_wf_cell_iter = 0
do j = 1, num_wann_cell_ll
num_wf_iterator = 0
! 2nd Loop over wannier functions in the present unit cell
do k = 1, num_wann_cell_ll
if ((.not. has_similar_centres(k)) .and. (j .ne. k)) then
coord_iterator = 0
! Loop over x,y,z to find similar centres
do l = 1, 3
if (i .le. 2*transport%num_cell_ll) then
if (abs(wannier_centres_translated(coord(l), tran_sorted_idx(j + (i - 1)*num_wann_cell_ll)) - &
wannier_centres_translated(coord(l), tran_sorted_idx(k + (i - 1)*num_wann_cell_ll))) &
.le. transport%group_threshold) then
coord_iterator = coord_iterator + 1
else
exit
end if
else
if (abs(wannier_centres_translated(coord(l), &
tran_sorted_idx(num_wann - 2*transport%num_ll &
+ j + (i - 2*transport%num_cell_ll - 1)*num_wann_cell_ll)) &
- wannier_centres_translated(coord(l), &
tran_sorted_idx(num_wann - 2*transport%num_ll &
+ k + (i - 2*transport%num_cell_ll - 1)*num_wann_cell_ll))) &
.le. transport%group_threshold) then
coord_iterator = coord_iterator + 1
else
exit
end if
end if
end do
if (coord_iterator .eq. 3) then
if (.not. has_similar_centres(j)) then
num_wf_iterator = num_wf_iterator + 1
if (i .le. 2*transport%num_cell_ll) then
idx_similar_wf(num_wf_iterator) = tran_sorted_idx(j + (i - 1)*num_wann_cell_ll)
else
idx_similar_wf(num_wf_iterator) = tran_sorted_idx(j + num_wann - 2*transport%num_ll + &
(i - 2*transport%num_cell_ll - 1)*num_wann_cell_ll)
end if
if (i .le. 2*transport%num_cell_ll) then
first_group_element(i, j) = j + (i - 1)*num_wann_cell_ll
else
first_group_element(i, j) = num_wann - 2*transport%num_ll + &
j + (i - 2*transport%num_cell_ll - 1)*num_wann_cell_ll
end if
num_wf_cell_iter = num_wf_cell_iter + 1
centre_id(num_wf_cell_iter) = j
end if
has_similar_centres(k) = .true.
has_similar_centres(j) = .true.
num_wf_iterator = num_wf_iterator + 1
if (i .le. 2*transport%num_cell_ll) then
idx_similar_wf(num_wf_iterator) = tran_sorted_idx(k + (i - 1)*num_wann_cell_ll)
else
idx_similar_wf(num_wf_iterator) = tran_sorted_idx(k + num_wann - 2*transport%num_ll + &
(i - 2*transport%num_cell_ll - 1)*num_wann_cell_ll)
end if
end if
end if
end do ! loop over k
if (num_wf_iterator .gt. 0) then
group_iterator = group_iterator + 1
wf_similar_centres(i, group_iterator, :) = idx_similar_wf(:)
!Save number of WFs in each group
tmp_wf_verifier(i, group_iterator) = num_wf_iterator
end if
end do
if ((count(has_similar_centres) .eq. 0) .and. (i .eq. 1)) then
write (stdout, '(a)') ' No wannier functions found with similar centres: sorting completed'
exit
elseif (i .eq. 1) then
write (stdout, *) ' Wannier functions found with similar centres: '
write (stdout, *) ' -> using signatures to complete sorting '
end if
!Save number of group of WFs in each unit cell and compare to previous unit cell
group_verifier(i) = group_iterator
if (print_output%iprint .ge. 4) write (stdout, '(a11,i4,a13,i4)') ' Unit cell:', i, ' Num groups:', group_verifier(i)
if (i .ne. 1) then
if (group_verifier(i) .ne. group_verifier(i - 1)) then
if (write_xyz) call tran_write_xyz(atom_data, transport, wannier_centres_translated, &
tran_sorted_idx, num_wann, seedname, stdout)
call set_error_fatal(error, 'Inconsistent number of groups of similar centred wannier functions between unit cells', comm)
return
elseif (i .eq. 4*transport%num_cell_ll) then
write (stdout, *) ' Consistent groups of similar centred wannier functions between '
write (stdout, *) ' unit cells found'
write (stdout, *) ' '
end if
end if
end do !Loop over all unit cells in PL1,PL2,PL3,PL4
! Perform check to ensure consistent number of WFs between equivalent groups in different unit cells
if (any(has_similar_centres)) then
allocate (wf_verifier(4*transport%num_cell_ll, group_verifier(1)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating wf_verifier in check_and_sort_similar_centres', comm)
return
end if
if (print_output%iprint .ge. 4) write (stdout, *) 'Unit cell Group number Num WFs'
wf_verifier = 0
wf_verifier = tmp_wf_verifier(:, 1:group_verifier(1))
do i = 1, 4*transport%num_cell_ll
do j = 1, group_verifier(1)
if (print_output%iprint .ge. 4) write (stdout, '(a3,i4,a9,i4,a7,i4)') ' ', i, ' ', &
j, ' ', wf_verifier(i, j)
if (i .ne. 1) then
if (wf_verifier(i, j) .ne. wf_verifier(i - 1, j)) then
call set_error_fatal(error, 'Inconsistent number of wannier &
&functions between equivalent groups of similar &
¢red wannier functions', comm)
return
end if
end if
end do
end do
write (stdout, *) ' Consistent number of wannier functions between equivalent groups of similar'
write (stdout, *) ' centred wannier functions'
write (stdout, *) ' '
write (stdout, *) ' Fixing order of similar centred wannier functions using parity signatures'
do i = 2, 4*transport%num_cell_ll
do j = 1, group_verifier(1)
! Make array of WF numbers which act as a reference to sort against
! and an array which need sorting
allocate (ref_similar_centres(group_verifier(1), wf_verifier(1, j)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating ref_similar_centres in check_and_sort_similar_centres', comm)
return
end if
allocate (unsorted_similar_centres(group_verifier(1), wf_verifier(1, j)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating unsorted_similar_centres in check_and_sort_similar_centres', comm)
return
end if
allocate (sorted_idx(wf_verifier(1, j)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating sorted_idx in check_and_sort_similar_centres', comm)
return
end if
allocate (dot_p(wf_verifier(1, j)), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating dot_p in check_and_sort_similar_centres', comm)
return
end if
do k = 1, wf_verifier(1, j)
ref_similar_centres(j, k) = wf_similar_centres(1, j, k)
unsorted_similar_centres(j, k) = wf_similar_centres(i, j, k)
end do
sorted_idx = 0
do k = 1, wf_verifier(1, j)
dot_p = 0.0_dp
! building the array of positive dot products of signatures between unsorted_similar_centres(j,k)
! and all the ref_similar_centres(j,:)
do l = 1, wf_verifier(1, j)
do p = 1, num_G
dot_p(l) = dot_p(l) + abs(signatures(p, unsorted_similar_centres(j, k)))* &
abs(signatures(p, ref_similar_centres(j, l)))
end do
end do
max_position = maxloc(dot_p)
sorted_idx(max_position(1)) = unsorted_similar_centres(j, k)
end do
! we have the properly ordered indexes for group j in unit cell i, now we need
! to overwrite the tran_sorted_idx array at the proper position
tran_sorted_idx(first_group_element(i, centre_id(j)):first_group_element(i, centre_id(j)) +&
&wf_verifier(i, j) - 1) = sorted_idx(:)
deallocate (dot_p, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating dot_p in check_and_sort_similar_centres', comm)
return
end if
deallocate (sorted_idx, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating sorted_idx in check_and_sort_similar_centres', comm)
return
end if
deallocate (unsorted_similar_centres, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating unsorted_similar_centres in check_and_sort_similar_centres', comm)
return
end if
deallocate (ref_similar_centres, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating ref_similar_centres in check_and_sort_similar_centres', comm)
return
end if
end do
end do
! checking that all the indices of WFs in the new tran_sorted_idx are distinct
! Remark: physically, no two WFs with similar centres can have the same type so we should expect
! this check to always pass unless the signatures/wf are very weird !!
do k = 1, num_wann
iterator = 0
do l = 1, num_wann
if (tran_sorted_idx(l) .eq. k) then
iterator = iterator + 1
end if
end do
if ((iterator .ge. 2) .or. (iterator .eq. 0)) then
call set_error_fatal(error, &
'A Wannier Function appears either zero times or twice after sorting, this may be due to a &
&poor wannierisation and/or disentanglement', comm)
return
end if
!write(stdout,*) ' WF : ',k,' appears ',iterator,' time(s)'
end do
deallocate (wf_verifier, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating wf_verifier in check_and_sort_similar_centres', comm)
return
end if
end if
deallocate (centre_id, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating centre_id in check_and_sort_similar_centres', comm)
return
end if
deallocate (first_group_element, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating first_group_element in check_and_sort_similar_centres', comm)
return
end if
deallocate (group_verifier, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating group_verifier in check_and_sort_similar_centres', comm)
return
end if
deallocate (tmp_wf_verifier, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating tmp_wf_verifier in check_and_sort_similar_centres', comm)
return
end if
deallocate (has_similar_centres, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating has_similar_centres in check_and_sort_similar_centres', comm)
return
end if
deallocate (idx_similar_wf, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating idx_similar_wf in check_and_sort_similar_centres', comm)
return
end if
deallocate (wf_similar_centres, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating wf_similar_centre in check_and_sort_similar_centres', comm)
return
end if
if (print_output%timing_level > 2) call io_stopwatch_stop('tran: lcr_2c2_sort: similar_centres', timer)
return
end subroutine check_and_sort_similar_centres
!================================================!
subroutine tran_write_xyz(atom_data, transport, wannier_centres_translated, tran_sorted_idx, &
num_wann, seedname, stdout)
!================================================!
!
! Write xyz file with Wannier centres
! and atomic positions
!
!================================================!
use w90_io, only: io_date
use w90_types, only: atom_data_type
use w90_wannier90_types, only: transport_type
implicit none
! arguments
type(transport_type), intent(inout) :: transport
type(atom_data_type), intent(in) :: atom_data
integer, intent(in) :: tran_sorted_idx(:)
integer, intent(in) :: num_wann
integer, intent(in) :: stdout
real(kind=dp), intent(in) :: wannier_centres_translated(:, :)
character(len=50), intent(in) :: seedname
! local variables
integer :: iw, ind, xyz_unit, nat, nsp
character(len=9) :: cdate, ctime
real(kind=dp) :: wc(3, num_wann)
if (index(transport%mode, 'bulk') > 0) wc = wannier_centres_translated
if (index(transport%mode, 'lcr') > 0) then
do iw = 1, num_wann
wc(:, iw) = wannier_centres_translated(:, tran_sorted_idx(iw))
end do
end if
open (newunit=xyz_unit, file=trim(seedname)//'_centres.xyz', form='formatted')
write (xyz_unit, '(i6)') num_wann + atom_data%num_atoms
call io_date(cdate, ctime)
write (xyz_unit, '(a84)') 'Wannier centres and atomic positions, written by Wannier90 on '//cdate//' at '//ctime
do iw = 1, num_wann
write (xyz_unit, '("X",6x,3(f14.8,3x))') (wc(ind, iw), ind=1, 3)
end do
do nsp = 1, atom_data%num_species
do nat = 1, atom_data%species_num(nsp)
write (xyz_unit, '(a2,5x,3(f14.8,3x))') atom_data%symbol(nsp), atom_data%pos_cart(:, nat, nsp)
end do
end do
write (stdout, *) ' Wannier centres written to file '//trim(seedname)//'_centres.xyz'
return
end subroutine tran_write_xyz
!================================================!
subroutine tran_parity_enforce(signatures, print_output, transport, num_wann, tran_sorted_idx, &
hr_one_dim, irvec_max, stdout, timer)
!================================================!
! Here, the signatures of the each wannier fucntion (stored in
! signatures) is used to determine its relavite parity
! with respect to the first unit cell. The parity pattern of
! first unit cell is then enforced.
!================================================!
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
implicit none
! arguments
integer, intent(in) :: irvec_max
integer, intent(in) :: num_wann
integer, intent(in) :: stdout
integer, intent(in) :: tran_sorted_idx(:)
real(kind=dp), intent(inout) :: signatures(:, :)
real(kind=dp), intent(inout) :: hr_one_dim(:, :, -irvec_max:)
type(print_output_type), intent(in) :: print_output
type(transport_type), intent(in) :: transport
type(timer_list_type), intent(inout) :: timer
! local variables
integer :: i, j, k, wf_idx, num_wann_cell_ll
real(kind=dp) :: signature_dot_p
if (print_output%timing_level > 1) call io_stopwatch_start('tran: parity_enforce', timer)
! NP: special "easy" fix of the parities by switching the sign
! of the Wannier Functions if the first element of the signature
! is found negative. Then updating the signature and the Hamiltonian
! matrix element for the corresponding line and column
if (transport%easy_fix) then
do i = 1, num_wann
if (real(signatures(1, i)) .lt. 0.0_dp) then
signatures(:, i) = -signatures(:, i)
do k = 1, num_wann
hr_one_dim(k, i, 0) = -hr_one_dim(k, i, 0)
hr_one_dim(i, k, 0) = -hr_one_dim(i, k, 0)
end do
end if
end do
end if
num_wann_cell_ll = transport%num_ll/transport%num_cell_ll
if (print_output%iprint .eq. 5) write (stdout, '(a101)') 'Unit cell Sorted WF index Unsort WF index &
&Unsorted WF Equiv Signature Dot Product'
! Loop over unit cell in principal layers
do i = 2, 4*transport%num_cell_ll
! Loop over wannier functions in unit cell
do j = 1, num_wann_cell_ll
if (i .le. 2*transport%num_cell_ll) then
wf_idx = j + (i - 1)*num_wann_cell_ll
else
wf_idx = num_wann - 2*transport%num_ll + j + (i - 1 - 2*transport%num_cell_ll)*num_wann_cell_ll
end if
signature_dot_p = dot_product(signatures(:, tran_sorted_idx(j)), signatures(:, tran_sorted_idx(wf_idx)))
if (print_output%iprint .eq. 5) then
write (stdout, '(2x,i4,3(13x,i5),12x,f20.17)') &
i, wf_idx, tran_sorted_idx(wf_idx), tran_sorted_idx(j), signature_dot_p
end if
if (abs(signature_dot_p) .le. 0.8_dp) then
write (stdout, '(a28,i4,a64,i4,a20)') ' WARNING: Wannier function (', tran_sorted_idx(wf_idx), &
') seems to has poor resemblance to equivalent wannier function (', tran_sorted_idx(j), ') in first unit cell'
if (print_output%iprint .lt. 5) write (stdout, *) 'Dot product of signatures: ', signature_dot_p
end if
if (signature_dot_p .lt. 0.0_dp) then
do k = 1, num_wann
hr_one_dim(k, tran_sorted_idx(wf_idx), 0) = -hr_one_dim(k, tran_sorted_idx(wf_idx), 0)
hr_one_dim(tran_sorted_idx(wf_idx), k, 0) = -hr_one_dim(tran_sorted_idx(wf_idx), k, 0)
end do
end if
end do
end do
if (print_output%timing_level > 1) call io_stopwatch_stop('tran: parity_enforce', timer)
return
end subroutine tran_parity_enforce
!================================================!
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
end module w90_transport_mod
