subroutine tran_lcr_2c2_sort(signatures, num_G, pl_warning)
!=======================================================!
! 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_error, stdout, io_stopwatch
use w90_parameters, only: one_dim_dir, tran_num_ll, num_wann, tran_num_cell_ll, &
real_lattice, tran_group_threshold, iprint, timing_level, lenconfac, &
wannier_spreads, write_xyz, dist_cutoff
use w90_hamiltonian, only: wannier_centres_translated
implicit none
integer, intent(in) :: num_G
real(dp), intent(in), dimension(:, :) :: signatures
logical, intent(out) :: pl_warning
real(dp), dimension(2, num_wann) :: centres_non_sorted, centres_initial_sorted
real(dp), dimension(2, tran_num_ll) :: PL1, PL2, PL3, PL4, PL
real(dp), dimension(2, num_wann - (4*tran_num_ll)) :: central_region
real(dp) :: reference_position, &
cell_length, distance, PL_max_val, PL_min_val
!~ integer :: l,max_i,iterator !aam: unused variables
integer :: i, j, k, PL_selector, &
sort_iterator, sort_iterator2, ierr, temp_coord_2, temp_coord_3, n, &
num_wann_cell_ll, num_wf_group1, num_wf_last_group
integer, allocatable, dimension(:) :: PL_groups, &
PL1_groups, PL2_groups, PL3_groups, PL4_groups, central_region_groups
integer, allocatable, dimension(:, :) :: 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) call io_error('Error in allocating tran_sorted_idx in tran_lcr_2c2_sort')
num_wann_cell_ll = tran_num_ll/tran_num_cell_ll
if (timing_level > 1) call io_stopwatch('tran: lcr_2c2_sort', 1)
sort_iterator = 0
!
!Check translated centres have been found
!
if (size(wannier_centres_translated) .eq. 0) then
call io_error('Translated centres not known : required perform lcr transport, try restart=plot')
endif
!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 (one_dim_dir .eq. 1) then
coord(1) = 1
coord(2) = 2
coord(3) = 3
elseif (one_dim_dir .eq. 2) then
coord(1) = 2
coord(2) = 3
coord(3) = 1
elseif (one_dim_dir .eq. 3) then
coord(1) = 3
coord(2) = 1
coord(3) = 2
endif
!
!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 io_error( &
'Lattice vector in conduction direction must point along x,y or z &
& direction and be orthogonal to the remaining lattice vectors.')
endif
!
!Check
!
if (num_wann .le. 4*tran_num_ll) then
call io_error('Principle layers are too big.')
endif
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)
enddo
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:tran_num_ll)
PL2 = centres_initial_sorted(:, tran_num_ll + 1:2*tran_num_ll)
PL3 = centres_initial_sorted(:, num_wann - (2*tran_num_ll - 1):num_wann - (tran_num_ll))
PL4 = centres_initial_sorted(:, num_wann - (tran_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
endif
!
if (iprint .ge. 4) then
write (stdout, *) ' Group Breakdown of each principal layer'
endif
!
!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
endselect
!
!Grouping wannier functions with similar coord(1)
!
call group(PL, PL_groups)
if (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)), ')'
endif
!
!Returns the sorted PL and informations on this PL
!
allocate (PL_subgroup_info(size(PL_groups), maxval(PL_groups)), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating PL_subgroup_info in tran_lcr_2c2_sort')
call master_sort_and_group(PL, PL_groups, tran_num_ll, PL_subgroup_info)
select case (PL_selector)
case (1)
allocate (PL1_groups(size(PL_groups)), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating PL1_groups in tran_lcr_2c2_sort')
allocate (PL1_subgroup_info(size(PL_groups), maxval(PL_groups)), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating PL1_subgroup_info in tran_lcr_2c2_sort')
!
PL1 = PL
PL1_groups = PL_groups
PL1_subgroup_info = PL_subgroup_info
!
deallocate (PL_subgroup_info, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL1_subgroup_info in tran_lcr_2c2_sort')
case (2)
allocate (PL2_groups(size(PL_groups)), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating PL2_groups in tran_lcr_2c2_sort')
allocate (PL2_subgroup_info(size(PL_groups), maxval(PL_groups)), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating PL2_subgroup_info in tran_lcr_2c2_sort')
!
PL2 = PL
PL2_groups = PL_groups
PL2_subgroup_info = PL_subgroup_info
!
deallocate (PL_subgroup_info, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL2_subgroup_info in tran_lcr_2c2_sort')
case (3)
allocate (PL3_groups(size(PL_groups)), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating PL3_groups in tran_lcr_2c2_sort')
allocate (PL3_subgroup_info(size(PL_groups), maxval(PL_groups)), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating PL3_subgroup_info in tran_lcr_2c2_sort')
!
PL3 = PL
PL3_groups = PL_groups
PL3_subgroup_info = PL_subgroup_info
!
deallocate (PL_subgroup_info, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL3_subgroup_info in tran_lcr_2c2_sort')
case (4)
allocate (PL4_groups(size(PL_groups)), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating PL4_groups in tran_lcr_2c2_sort')
allocate (PL4_subgroup_info(size(PL_groups), maxval(PL_groups)), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating PL4_subgroup_info in tran_lcr_2c2_sort')
!
PL4 = PL
PL4_groups = PL_groups
PL4_subgroup_info = PL_subgroup_info
!
deallocate (PL_subgroup_info, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL4_subgroup_info in tran_lcr_2c2_sort')
endselect
deallocate (PL_groups, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL_groups in tran_lcr_2c2_sort')
enddo ! Principal layer loop
!
!Grouping and sorting of central conductor region
!
!Define central region
!
central_region = centres_initial_sorted(:, 2*tran_num_ll + 1:num_wann - (2*tran_num_ll))
!
!Group central region
!
call group(central_region, central_region_groups)
!
!Print central region group breakdown
!
if (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)), ')'
endif
!
!Returns sorted central group region
!
allocate (central_subgroup_info(size(central_region_groups), maxval(central_region_groups)), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating central_group_info in tran_lcr_2c2_sort')
call master_sort_and_group(central_region, central_region_groups, num_wann - (4*tran_num_ll), central_subgroup_info)
deallocate (central_subgroup_info, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating central_group_info in tran_lcr_2c2_sort')
write (stdout, *) ' '
!
!Build the sorted index array
!
tran_sorted_idx = nint(centres_initial_sorted(1, :))
tran_sorted_idx(1:tran_num_ll) = nint(PL1(1, :))
tran_sorted_idx(tran_num_ll + 1:2*tran_num_ll) = nint(PL2(1, :))
tran_sorted_idx(2*tran_num_ll + 1:num_wann - (2*tran_num_ll)) = nint(central_region(1, :))
tran_sorted_idx(num_wann - (2*tran_num_ll - 1):num_wann - (tran_num_ll)) = nint(PL3(1, :))
tran_sorted_idx(num_wann - (tran_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()
call io_error('Sorting techniques exhausted:&
& Inconsistent number of groups among principal layers')
endif
write (stdout, *) 'Inconsistent number of groups among principal layers: restarting sorting...'
deallocate (PL1_groups, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL1_groups in tran_lcr_2c2_sort')
deallocate (PL2_groups, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL2_groups in tran_lcr_2c2_sort')
deallocate (PL3_groups, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL3_groups in tran_lcr_2c2_sort')
deallocate (PL4_groups, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL4_groups in tran_lcr_2c2_sort')
deallocate (PL1_subgroup_info)
if (ierr /= 0) call io_error('Error deallocating PL1_subgroup_info in tran_lcr_2c2_sort')
deallocate (PL2_subgroup_info, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL2_subgroup_info in tran_lcr_2c2_sort')
deallocate (PL3_subgroup_info, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL3_subgroup_info in tran_lcr_2c2_sort')
deallocate (PL4_subgroup_info, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL4_subgroup_info in tran_lcr_2c2_sort')
goto 100
endif
!
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()
call io_error &
('Sorting techniques exhausted: Inconsitent number of wannier function among &
& similar groups within principal layers')
endif
write (stdout, *) 'Inconsitent number of wannier function among &
&similar groups within& principal layers: restarting sorting...'
deallocate (PL1_groups, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL1_groups in tran_lcr_2c2_sort')
deallocate (PL2_groups, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL2_groups in tran_lcr_2c2_sort')
deallocate (PL3_groups, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL3_groups in tran_lcr_2c2_sort')
deallocate (PL4_groups, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL4_groups in tran_lcr_2c2_sort')
deallocate (PL1_subgroup_info)
if (ierr /= 0) call io_error('Error deallocating PL1_subgroup_info in tran_lcr_2c2_sort')
deallocate (PL2_subgroup_info, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL2_subgroup_info in tran_lcr_2c2_sort')
deallocate (PL3_subgroup_info, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL3_subgroup_info in tran_lcr_2c2_sort')
deallocate (PL4_subgroup_info, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL4_subgroup_info in tran_lcr_2c2_sort')
goto 100
endif
enddo
!
!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, tran_num_ll
distance = abs(abs(reference_position - wannier_centres_translated(coord(1), tran_sorted_idx(num_wann - i + 1))) &
- cell_length)
if (distance .lt. tran_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
endif
enddo
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) call io_error('Error deallocating hr_one_dim in tran_lcr_2c2_sort')
call tran_reduce_hr()
call tran_cut_hr_one_dim()
write (stdout, *) ' '
write (stdout, *) ' Restarting sorting...'
write (stdout, *) ' '
sort_iterator = sort_iterator - 1
deallocate (PL1_groups, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL1_groups in tran_lcr_2c2_sort')
deallocate (PL2_groups, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL2_groups in tran_lcr_2c2_sort')
deallocate (PL3_groups, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL3_groups in tran_lcr_2c2_sort')
deallocate (PL4_groups, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL4_groups in tran_lcr_2c2_sort')
deallocate (PL1_subgroup_info)
if (ierr /= 0) call io_error('Error deallocating PL1_subgroup_info in tran_lcr_2c2_sort')
deallocate (PL2_subgroup_info, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL2_subgroup_info in tran_lcr_2c2_sort')
deallocate (PL3_subgroup_info, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL3_subgroup_info in tran_lcr_2c2_sort')
deallocate (PL4_subgroup_info, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL4_subgroup_info in tran_lcr_2c2_sort')
goto 100
endif
!
! 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) call io_error('Error in allocating tmp_subgroup in tran_lcr_2c2_sort')
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()
call io_error &
('Sorting techniques exhausted: Inconsitent subgroup structures among principal layers')
endif
write (stdout, *) 'Inconsitent subgroup structure among principal layers: restarting sorting...'
deallocate (temp_subgroup, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating tmp_subgroup in tran_lcr_2c2_sort')
deallocate (PL1_groups, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL1_groups in tran_lcr_2c2_sort')
deallocate (PL2_groups, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL2_groups in tran_lcr_2c2_sort')
deallocate (PL3_groups, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL3_groups in tran_lcr_2c2_sort')
deallocate (PL4_groups, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL4_groups in tran_lcr_2c2_sort')
deallocate (PL1_subgroup_info)
if (ierr /= 0) call io_error('Error deallocating PL1_subgroup_info in tran_lcr_2c2_sort')
deallocate (PL2_subgroup_info, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL2_subgroup_info in tran_lcr_2c2_sort')
deallocate (PL3_subgroup_info, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL3_subgroup_info in tran_lcr_2c2_sort')
deallocate (PL4_subgroup_info, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating PL4_subgroup_info in tran_lcr_2c2_sort')
goto 100
endif
enddo
enddo
enddo
!
! 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)
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*tran_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_spreads(tran_sorted_idx(n))*lenconfac**2
enddo
write (stdout, *) '==================================PL3========================================='
endif
if (k .eq. 4) write (stdout, *) '==================================PL4========================================='
do i = 1, tran_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_spreads(tran_sorted_idx(n))*lenconfac**2
enddo
if (i .ne. tran_num_cell_ll) write (stdout, *) '---------------------&
&---------------------------------------------------------'
enddo
enddo
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
endif
PL_min_val = maxval(wannier_centres_translated(coord(1), tran_sorted_idx(tran_num_ll - num_wf_last_group + 1:tran_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(tran_num_ll + 1:tran_num_ll + num_wf_group1))) &
- minval(wannier_centres_translated(coord(1), tran_sorted_idx(1:num_wf_group1)))
if ((dist_cutoff .lt. PL_min_val) .or. (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.
endif
!
! End MS.
!
if (timing_level > 1) call io_stopwatch('tran: lcr_2c2_sort', 2)
return
end subroutine tran_lcr_2c2_sort