subroutine tran_find_integral_signatures(signatures, num_G)
!=========================================================================!
! 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_error, stdout, seedname, io_file_unit, io_date, &
io_stopwatch
use w90_parameters, only: num_wann, have_disentangled, num_bands, u_matrix, u_matrix_opt, &
real_lattice, iprint, timing_level
use w90_hamiltonian, only: wannier_centres_translated
implicit none
integer, intent(out) :: num_G
real(kind=dp), allocatable, dimension(:, :), intent(out) :: signatures
complex(kind=dp), allocatable :: unkg(:, :), tran_u_matrix(:, :)
complex(kind=dp) :: phase_factor, signature_basis(32)
real(kind=dp) :: i_unkg, r_unkg, wf_frac(3), det_rl, inv_t_rl(3, 3), &
mag_signature_sq
!~ character(len=11) :: unkg_file
logical :: have_file
integer, allocatable, dimension(:, :) :: g_abc
integer :: i, ibnd, file_unit, ierr, p, p_max, n, m, ig, a, b, c, ig_idx(32)
if (timing_level > 1) call io_stopwatch('tran: find_sigs_unkg_int', 1)
!
file_unit = io_file_unit()
inquire (file=trim(seedname)//'.unkg', exist=have_file)
if (.not. have_file) call io_error('tran_hr_parity_unkg: file '//trim(seedname)//'.unkg not found')
open (file_unit, file=trim(seedname)//'.unkg', form='formatted', action='read')
!
!Read unkg file
!
write (stdout, '(3a)') ' Reading '//trim(seedname)//'.unkg file'
read (file_unit, *) num_G
allocate (signatures(20, num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating signatures in tran_find_sigs_unkg_int')
allocate (unkg(num_G, num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating unkg in tran_find_sigs_unkg_int')
allocate (g_abc(num_G, 3), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating g_abc in tran_find_sigs_unkg_int')
if (have_disentangled) then
allocate (tran_u_matrix(num_bands, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating tran_u_matrix in tran_find_sigs_unkg_int')
else
allocate (tran_u_matrix(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating tran_u_matrix in tran_find_sigs_unkg_int')
endif
!
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 io_error('tran_find_sigs_unkg_int: Incorrect bands or g vectors')
endif
unkg(i, m) = cmplx(r_unkg, i_unkg, kind=dp)
g_abc(i, :) = (/a, b, c/)
enddo
enddo
!
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)
enddo
enddo
p_max = num_bands
else
tran_u_matrix = u_matrix(:, :, 1)
p_max = num_wann
endif
enddo
if (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
enddo
!
! 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))
enddo
signature_basis(ig) = signature_basis(ig)*phase_factor
enddo
!
! 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 (iprint .ge. 5) then
write (stdout, *) ' '
write (stdout, *) ' Wannier function: ', n
do ig = 1, 20
write (stdout, *) ig - 1, signatures(ig, n)
enddo
endif
!
!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
enddo
signatures(:, n) = signatures(:, n)/sqrt(mag_signature_sq)
!
enddo ! 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) call io_error('Error deallocating tran_u_matrix in tran_find_signatures')
deallocate (g_abc, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating g_abc in tran_find_signatures')
deallocate (unkg, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating unkg in tran_find_signatures')
if (timing_level > 1) call io_stopwatch('tran: find_sigs_unkg_int', 2)
return
end subroutine tran_find_integral_signatures