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
