Wannier representation of the Pauli matrices: <0n|sigma_a|Rm> (a=x,y,z)
subroutine get_SS_R
!================================================================
!
!! Wannier representation of the Pauli matrices: <0n|sigma_a|Rm>
!! (a=x,y,z)
!
!================================================================
use w90_constants, only: dp, pi, cmplx_0
use w90_parameters, only: num_wann, ndimwin, num_kpts, num_bands, &
timing_level, have_disentangled, spn_formatted
use w90_postw90_common, only: nrpts, v_matrix
use w90_io, only: io_error, io_stopwatch, stdout, seedname, &
io_file_unit
use w90_comms, only: on_root, comms_bcast
implicit none
complex(kind=dp), allocatable :: spn_o(:, :, :, :), SS_q(:, :, :, :), spn_temp(:, :)
real(kind=dp) :: s_real, s_img
integer, allocatable :: num_states(:)
integer :: i, j, ii, jj, m, n, spn_in, ik, is, &
winmin, nb_tmp, nkp_tmp, ierr, s, counter
character(len=60) :: header
if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_SS_R', 1)
if (.not. allocated(SS_R)) then
allocate (SS_R(num_wann, num_wann, nrpts, 3))
else
return ! been here before
end if
if (on_root) then
allocate (spn_o(num_bands, num_bands, num_kpts, 3))
allocate (SS_q(num_wann, num_wann, num_kpts, 3))
allocate (num_states(num_kpts))
do ik = 1, num_kpts
if (have_disentangled) then
num_states(ik) = ndimwin(ik)
else
num_states(ik) = num_wann
endif
enddo
! Read from .spn file the original spin matrices <psi_nk|sigma_i|psi_mk>
! (sigma_i = Pauli matrix) between ab initio eigenstates
!
spn_in = io_file_unit()
if (spn_formatted) then
open (unit=spn_in, file=trim(seedname)//'.spn', form='formatted', &
status='old', err=109)
write (stdout, '(/a)', advance='no') &
' Reading spin matrices from '//trim(seedname)//'.spn in get_SS_R : '
read (spn_in, *, err=110, end=110) header
write (stdout, '(a)') trim(header)
read (spn_in, *, err=110, end=110) nb_tmp, nkp_tmp
else
open (unit=spn_in, file=trim(seedname)//'.spn', form='unformatted', &
status='old', err=109)
write (stdout, '(/a)', advance='no') &
' Reading spin matrices from '//trim(seedname)//'.spn in get_SS_R : '
read (spn_in, err=110, end=110) header
write (stdout, '(a)') trim(header)
read (spn_in, err=110, end=110) nb_tmp, nkp_tmp
endif
if (nb_tmp .ne. num_bands) &
call io_error(trim(seedname)//'.spn has wrong number of bands')
if (nkp_tmp .ne. num_kpts) &
call io_error(trim(seedname)//'.spn has wrong number of k-points')
if (spn_formatted) then
do ik = 1, num_kpts
do m = 1, num_bands
do n = 1, m
read (spn_in, *, err=110, end=110) s_real, s_img
spn_o(n, m, ik, 1) = cmplx(s_real, s_img, dp)
read (spn_in, *, err=110, end=110) s_real, s_img
spn_o(n, m, ik, 2) = cmplx(s_real, s_img, dp)
read (spn_in, *, err=110, end=110) s_real, s_img
spn_o(n, m, ik, 3) = cmplx(s_real, s_img, dp)
! Read upper-triangular part, now build the rest
spn_o(m, n, ik, 1) = conjg(spn_o(n, m, ik, 1))
spn_o(m, n, ik, 2) = conjg(spn_o(n, m, ik, 2))
spn_o(m, n, ik, 3) = conjg(spn_o(n, m, ik, 3))
end do
end do
enddo
else
allocate (spn_temp(3, (num_bands*(num_bands + 1))/2), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating spm_temp in get_SS_R')
do ik = 1, num_kpts
read (spn_in) ((spn_temp(s, m), s=1, 3), m=1, (num_bands*(num_bands + 1))/2)
counter = 0
do m = 1, num_bands
do n = 1, m
counter = counter + 1
spn_o(n, m, ik, 1) = spn_temp(1, counter)
spn_o(m, n, ik, 1) = conjg(spn_temp(1, counter))
spn_o(n, m, ik, 2) = spn_temp(2, counter)
spn_o(m, n, ik, 2) = conjg(spn_temp(2, counter))
spn_o(n, m, ik, 3) = spn_temp(3, counter)
spn_o(m, n, ik, 3) = conjg(spn_temp(3, counter))
end do
end do
end do
deallocate (spn_temp, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating spm_temp in get_SS_R')
endif
close (spn_in)
! Transform to projected subspace, Wannier gauge
!
SS_q(:, :, :, :) = cmplx_0
do ik = 1, num_kpts
do is = 1, 3
call get_gauge_overlap_matrix( &
ik, num_states(ik), &
ik, num_states(ik), &
spn_o(:, :, ik, is), SS_q(:, :, ik, is))
enddo !is
enddo !ik
call fourier_q_to_R(SS_q(:, :, :, 1), SS_R(:, :, :, 1))
call fourier_q_to_R(SS_q(:, :, :, 2), SS_R(:, :, :, 2))
call fourier_q_to_R(SS_q(:, :, :, 3), SS_R(:, :, :, 3))
endif !on_root
call comms_bcast(SS_R(1, 1, 1, 1), num_wann*num_wann*nrpts*3)
if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_SS_R', 2)
return
109 call io_error &
('Error: Problem opening input file '//trim(seedname)//'.spn')
110 call io_error &
('Error: Problem reading input file '//trim(seedname)//'.spn')
end subroutine get_SS_R