Main routine to calculate the b-vectors
subroutine kmesh_get()
!=====================================================
!
!! Main routine to calculate the b-vectors
!
!=====================================================
use w90_io, only: stdout, io_error, io_stopwatch
use w90_utility, only: utility_compar
implicit none
! Variables that are private
integer :: nlist, nkp, nkp2, l, m, n, ndnn, ndnnx, ndnntot
integer :: nnsh, nn, nnx, loop, i, j
integer :: ifound, counter, na, nap, loop_s, loop_b, shell, nbvec, bnum
integer :: ifpos, ifneg, ierr, multi(search_shells)
integer :: nnshell(num_kpts, search_shells)
integer, allocatable :: nnlist_tmp(:, :), nncell_tmp(:, :, :) ![ysl]
real(kind=dp) :: vkpp(3), vkpp2(3)
real(kind=dp) :: dist, dnn0, dnn1, bb1, bbn, ddelta
real(kind=dp), parameter :: eta = 99999999.0_dp ! eta = very large
real(kind=dp) :: bweight(max_shells)
real(kind=dp) :: dnn(search_shells)
real(kind=dp) :: wb_local(num_nnmax)
real(kind=dp) :: bk_local(3, num_nnmax, num_kpts), kpbvec(3)
real(kind=dp), allocatable :: bvec_tmp(:, :)
! Integer arrays that are public
if (timing_level > 0) call io_stopwatch('kmesh: get', 1)
if (on_root) write (stdout, '(/1x,a)') &
'*---------------------------------- K-MESH ----------------------------------*'
! Sort the cell neighbours so we loop in order of distance from the home shell
call kmesh_supercell_sort
! find the distance between k-point 1 and its nearest-neighbour shells
! if we have only one k-point, the n-neighbours are its periodic images
dnn0 = 0.0_dp
dnn1 = eta
ndnntot = 0
do nlist = 1, search_shells
do nkp = 1, num_kpts
do loop = 1, (2*nsupcell + 1)**3
l = lmn(1, loop); m = lmn(2, loop); n = lmn(3, loop)
!
vkpp = kpt_cart(:, nkp) + matmul(lmn(:, loop), recip_lattice)
dist = sqrt((kpt_cart(1, 1) - vkpp(1))**2 &
+ (kpt_cart(2, 1) - vkpp(2))**2 + (kpt_cart(3, 1) - vkpp(3))**2)
!
if ((dist .gt. kmesh_tol) .and. (dist .gt. dnn0 + kmesh_tol)) then
if (dist .lt. dnn1 - kmesh_tol) then
dnn1 = dist ! found a closer shell
counter = 0
end if
if (dist .gt. (dnn1 - kmesh_tol) .and. dist .lt. (dnn1 + kmesh_tol)) then
counter = counter + 1 ! count the multiplicity of the shell
end if
end if
enddo
enddo
if (dnn1 .lt. eta - kmesh_tol) ndnntot = ndnntot + 1
dnn(nlist) = dnn1
multi(nlist) = counter
dnn0 = dnn1
dnn1 = eta
enddo
if (on_root) then
write (stdout, '(1x,a)') '+----------------------------------------------------------------------------+'
write (stdout, '(1x,a)') '| Distance to Nearest-Neighbour Shells |'
write (stdout, '(1x,a)') '| ------------------------------------ |'
if (lenconfac .eq. 1.0_dp) then
write (stdout, '(1x,a)') '| Shell Distance (Ang^-1) Multiplicity |'
write (stdout, '(1x,a)') '| ----- ----------------- ------------ |'
else
write (stdout, '(1x,a)') '| Shell Distance (Bohr^-1) Multiplicity |'
write (stdout, '(1x,a)') '| ----- ------------------ ------------ |'
endif
do ndnn = 1, ndnntot
write (stdout, '(1x,a,11x,i3,17x,f10.6,19x,i4,12x,a)') '|', ndnn, dnn(ndnn)/lenconfac, multi(ndnn), '|'
enddo
write (stdout, '(1x,a)') '+----------------------------------------------------------------------------+'
endif
if (iprint >= 4) then
! Write out all the bvectors
if (on_root) then
write (stdout, '(1x,"|",76(" "),"|")')
write (stdout, '(1x,a)') '| Complete list of b-vectors and their lengths |'
write (stdout, '(1x,"|",76(" "),"|")')
write (stdout, '(1x,"+",76("-"),"+")')
endif
allocate (bvec_tmp(3, maxval(multi)), stat=ierr)
if (ierr /= 0) call io_error('Error allocating bvec_tmp in kmesh_get')
bvec_tmp = 0.0_dp
counter = 0
do shell = 1, search_shells
call kmesh_get_bvectors(multi(shell), 1, dnn(shell), bvec_tmp(:, 1:multi(shell)))
do loop = 1, multi(shell)
counter = counter + 1
if (on_root) write (stdout, '(a,I4,1x,a,2x,3f12.6,2x,a,2x,f12.6,a)') ' | b-vector ', counter, ': (', &
bvec_tmp(:, loop)/lenconfac, ')', dnn(shell)/lenconfac, ' |'
end do
end do
deallocate (bvec_tmp)
if (ierr /= 0) call io_error('Error deallocating bvec_tmp in kmesh_get')
if (on_root) write (stdout, '(1x,"|",76(" "),"|")')
if (on_root) write (stdout, '(1x,"+",76("-"),"+")')
end if
! Get the shell weights to satisfy the B1 condition
if (index(devel_flag, 'kmesh_degen') > 0) then
call kmesh_shell_from_file(multi, dnn, bweight)
else
if (num_shells == 0) then
call kmesh_shell_automatic(multi, dnn, bweight)
elseif (num_shells > 0) then
call kmesh_shell_fixed(multi, dnn, bweight)
end if
if (on_root) then
write (stdout, '(1x,a)', advance='no') '| The following shells are used: '
do ndnn = 1, num_shells
if (ndnn .eq. num_shells) then
write (stdout, '(i3,1x)', advance='no') shell_list(ndnn)
else
write (stdout, '(i3,",")', advance='no') shell_list(ndnn)
endif
enddo
do l = 1, 11 - num_shells
write (stdout, '(4x)', advance='no')
enddo
write (stdout, '("|")')
endif
end if
nntot = 0
do loop_s = 1, num_shells
nntot = nntot + multi(shell_list(loop_s))
end do
if (nntot > num_nnmax) then
if (on_root) then
write (stdout, '(a,i2,a)') ' **WARNING: kmesh has found >', num_nnmax, ' nearest neighbours**'
write (stdout, '(a)') ' '
write (stdout, '(a)') ' This is probably caused by an error in your unit cell specification'
write (stdout, '(a)') ' '
write (stdout, '(a)') ' If you think this is not the problem; please send your *.win file to the '
write (stdout, '(a)') ' wannier90 developers'
write (stdout, '(a)') ' '
write (stdout, '(a)') ' The problem may be caused by having accidentally degenerate shells of '
write (stdout, '(a)') ' kpoints. The solution is then to rerun wannier90 specifying the b-vectors '
write (stdout, '(a)') ' in each shell. Give devel_flag=kmesh_degen in the *.win file'
write (stdout, '(a)') ' and create a *.kshell file:'
write (stdout, '(a)') ' '
write (stdout, '(a)') ' $> cat hexagonal.kshell'
write (stdout, '(a)') ' $> 1 2'
write (stdout, '(a)') ' $> 5 6 7 8'
write (stdout, '(a)') ' '
write (stdout, '(a)') ' Where each line is a new shell (so num_shells in total)'
write (stdout, '(a)') ' The elements are the bvectors labelled according to the following '
write (stdout, '(a)') ' list (last column is distance)'
write (stdout, '(a)') ' '
endif
allocate (bvec_tmp(3, maxval(multi)), stat=ierr)
if (ierr /= 0) call io_error('Error allocating bvec_tmp in kmesh_get')
bvec_tmp = 0.0_dp
counter = 0
do shell = 1, search_shells
call kmesh_get_bvectors(multi(shell), 1, dnn(shell), bvec_tmp(:, 1:multi(shell)))
do loop = 1, multi(shell)
counter = counter + 1
if (on_root) write (stdout, '(a,I4,1x,a,2x,3f12.6,2x,a,2x,f12.6,a)') ' | b-vector ', counter, ': (', &
bvec_tmp(:, loop)/lenconfac, ')', dnn(shell)/lenconfac, ' |'
end do
end do
if (on_root) write (stdout, '(a)') ' '
deallocate (bvec_tmp)
if (ierr /= 0) call io_error('Error deallocating bvec_tmp in kmesh_get')
Call io_error('kmesh_get: something wrong, found too many nearest neighbours')
end if
allocate (nnlist(num_kpts, nntot), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating nnlist in kmesh_get')
allocate (neigh(num_kpts, nntot/2), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating neigh in kmesh_get')
allocate (nncell(3, num_kpts, nntot), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating nncell in kmesh_get')
allocate (wb(nntot), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating wb in kmesh_get')
allocate (bka(3, nntot/2), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating bka in kmesh_get')
allocate (bk(3, nntot, num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating bk in kmesh_get')
nnx = 0
do loop_s = 1, num_shells
do loop_b = 1, multi(shell_list(loop_s))
nnx = nnx + 1
wb_local(nnx) = bweight(loop_s)
end do
end do
! Now build up the list of nearest-neighbour shells for each k-point.
! nnlist(nkp,1...nnx) points to the nnx neighbours (ordered along increa
! shells) of the k-point nkp. nncell(i,nkp,nnth) tells us in which BZ is
! nnth nearest-neighbour of the k-point nkp. Construct the nnx b-vectors
! go from k-point nkp to each neighbour bk(1:3,nkp,1...nnx).
! Comment: Now we have bk(3,nntot,num_kps) 09/04/2006
if (on_root) then
write (stdout, '(1x,a)') '+----------------------------------------------------------------------------+'
write (stdout, '(1x,a)') '| Shell # Nearest-Neighbours |'
write (stdout, '(1x,a)') '| ----- -------------------- |'
endif
if (index(devel_flag, 'kmesh_degen') == 0) then
!
! Standard routine
!
nnshell = 0
do nkp = 1, num_kpts
nnx = 0
ok: do ndnnx = 1, num_shells
ndnn = shell_list(ndnnx)
do loop = 1, (2*nsupcell + 1)**3
l = lmn(1, loop); m = lmn(2, loop); n = lmn(3, loop)
vkpp2 = matmul(lmn(:, loop), recip_lattice)
do nkp2 = 1, num_kpts
vkpp = vkpp2 + kpt_cart(:, nkp2)
dist = sqrt((kpt_cart(1, nkp) - vkpp(1))**2 &
+ (kpt_cart(2, nkp) - vkpp(2))**2 + (kpt_cart(3, nkp) - vkpp(3))**2)
if ((dist .ge. dnn(ndnn)*(1 - kmesh_tol)) .and. (dist .le. dnn(ndnn)*(1 + kmesh_tol))) then
nnx = nnx + 1
nnshell(nkp, ndnn) = nnshell(nkp, ndnn) + 1
nnlist(nkp, nnx) = nkp2
nncell(1, nkp, nnx) = l
nncell(2, nkp, nnx) = m
nncell(3, nkp, nnx) = n
bk_local(:, nnx, nkp) = vkpp(:) - kpt_cart(:, nkp)
endif
!if we have the right number of neighbours we can exit
if (nnshell(nkp, ndnn) == multi(ndnn)) cycle ok
enddo
enddo
! check to see if too few neighbours here
end do ok
end do
else
!
! incase we set the bvectors explicitly
!
nnshell = 0
do nkp = 1, num_kpts
nnx = 0
ok2: do loop = 1, (2*nsupcell + 1)**3
l = lmn(1, loop); m = lmn(2, loop); n = lmn(3, loop)
vkpp2 = matmul(lmn(:, loop), recip_lattice)
do nkp2 = 1, num_kpts
vkpp = vkpp2 + kpt_cart(:, nkp2)
bnum = 0
do ndnnx = 1, num_shells
do nbvec = 1, multi(ndnnx)
bnum = bnum + 1
kpbvec = kpt_cart(:, nkp) + bvec_inp(:, nbvec, ndnnx)
dist = sqrt((kpbvec(1) - vkpp(1))**2 &
+ (kpbvec(2) - vkpp(2))**2 + (kpbvec(3) - vkpp(3))**2)
if (abs(dist) < kmesh_tol) then
nnx = nnx + 1
nnshell(nkp, ndnnx) = nnshell(nkp, ndnnx) + 1
nnlist(nkp, bnum) = nkp2
nncell(1, nkp, bnum) = l
nncell(2, nkp, bnum) = m
nncell(3, nkp, bnum) = n
bk_local(:, bnum, nkp) = bvec_inp(:, nbvec, ndnnx)
endif
enddo
end do
if (nnx == sum(multi)) exit ok2
end do
enddo ok2
! check to see if too few neighbours here
end do
end if
do ndnnx = 1, num_shells
ndnn = shell_list(ndnnx)
if (on_root) write (stdout, '(1x,a,24x,i3,13x,i3,33x,a)') '|', ndnn, nnshell(1, ndnn), '|'
end do
if (on_root) write (stdout, '(1x,"+",76("-"),"+")')
do nkp = 1, num_kpts
nnx = 0
do ndnnx = 1, num_shells
ndnn = shell_list(ndnnx)
do nnsh = 1, nnshell(nkp, ndnn)
bb1 = 0.0_dp
bbn = 0.0_dp
nnx = nnx + 1
do i = 1, 3
bb1 = bb1 + bk_local(i, nnx, 1)*bk_local(i, nnx, 1)
bbn = bbn + bk_local(i, nnx, nkp)*bk_local(i, nnx, nkp)
enddo
if (abs(sqrt(bb1) - sqrt(bbn)) .gt. kmesh_tol) then
if (on_root) write (stdout, '(1x,2f10.6)') bb1, bbn
call io_error('Non-symmetric k-point neighbours!')
endif
enddo
enddo
enddo
! now check that the completeness relation is satisfied for every kpoint
! We know it is true for kpt=1; but we check the rest to be safe.
! Eq. B1 in Appendix B PRB 56 12847 (1997)
if (.not. skip_B1_tests) then
do nkp = 1, num_kpts
do i = 1, 3
do j = 1, 3
ddelta = 0.0_dp
nnx = 0
do ndnnx = 1, num_shells
ndnn = shell_list(ndnnx)
do nnsh = 1, nnshell(1, ndnn)
nnx = nnx + 1
ddelta = ddelta + wb_local(nnx)*bk_local(i, nnx, nkp)*bk_local(j, nnx, nkp)
enddo
enddo
if ((i .eq. j) .and. (abs(ddelta - 1.0_dp) .gt. kmesh_tol)) then
if (on_root) write (stdout, '(1x,2i3,f12.8)') i, j, ddelta
call io_error('Eq. (B1) not satisfied in kmesh_get (1)')
endif
if ((i .ne. j) .and. (abs(ddelta) .gt. kmesh_tol)) then
if (on_root) write (stdout, '(1x,2i3,f12.8)') i, j, ddelta
call io_error('Eq. (B1) not satisfied in kmesh_get (2)')
endif
enddo
enddo
enddo
end if
if (on_root) write (stdout, '(1x,a)') '| Completeness relation is fully satisfied [Eq. (B1), PRB 56, 12847 (1997)] |'
if (on_root) write (stdout, '(1x,"+",76("-"),"+")')
!
wbtot = 0.0_dp
nnx = 0
do ndnnx = 1, num_shells
ndnn = shell_list(ndnnx)
do nnsh = 1, nnshell(1, ndnn)
nnx = nnx + 1
wbtot = wbtot + wb_local(nnx)
enddo
enddo
nnh = nntot/2
! make list of bka vectors from neighbours of first k-point
! delete any inverse vectors as you collect them
na = 0
do nn = 1, nntot
ifound = 0
if (na .ne. 0) then
do nap = 1, na
call utility_compar(bka(1, nap), bk_local(1, nn, 1), ifpos, ifneg)
if (ifneg .eq. 1) ifound = 1
enddo
endif
if (ifound .eq. 0) then
! found new vector to add to set
na = na + 1
bka(1, na) = bk_local(1, nn, 1)
bka(2, na) = bk_local(2, nn, 1)
bka(3, na) = bk_local(3, nn, 1)
endif
enddo
if (na .ne. nnh) call io_error('Did not find right number of bk directions')
if (on_root) then
if (lenconfac .eq. 1.0_dp) then
write (stdout, '(1x,a)') '| b_k Vectors (Ang^-1) and Weights (Ang^2) |'
write (stdout, '(1x,a)') '| ---------------------------------------- |'
else
write (stdout, '(1x,a)') '| b_k Vectors (Bohr^-1) and Weights (Bohr^2) |'
write (stdout, '(1x,a)') '| ------------------------------------------ |'
endif
write (stdout, '(1x,a)') '| No. b_k(x) b_k(y) b_k(z) w_b |'
write (stdout, '(1x,a)') '| --- -------------------------------- -------- |'
do i = 1, nntot
write (stdout, '(1x,"|",11x,i3,5x,3f12.6,3x,f10.6,8x,"|")') &
i, (bk_local(j, i, 1)/lenconfac, j=1, 3), wb_local(i)*lenconfac**2
enddo
write (stdout, '(1x,"+",76("-"),"+")')
if (lenconfac .eq. 1.0_dp) then
write (stdout, '(1x,a)') '| b_k Directions (Ang^-1) |'
write (stdout, '(1x,a)') '| ----------------------- |'
else
write (stdout, '(1x,a)') '| b_k Directions (Bohr^-1) |'
write (stdout, '(1x,a)') '| ------------------------ |'
endif
write (stdout, '(1x,a)') '| No. x y z |'
write (stdout, '(1x,a)') '| --- -------------------------------- |'
do i = 1, nnh
write (stdout, '(1x,"|",11x,i3,5x,3f12.6,21x,"|")') i, (bka(j, i)/lenconfac, j=1, 3)
enddo
write (stdout, '(1x,"+",76("-"),"+")')
write (stdout, *) ' '
endif
! find index array
do nkp = 1, num_kpts
do na = 1, nnh
! first, zero the index array so we can check it gets filled
neigh(nkp, na) = 0
! now search through list of neighbours of this k-point
do nn = 1, nntot
call utility_compar(bka(1, na), bk_local(1, nn, nkp), ifpos, ifneg)
if (ifpos .eq. 1) neigh(nkp, na) = nn
enddo
! check found
if (neigh(nkp, na) .eq. 0) then
if (on_root) write (stdout, *) ' nkp,na=', nkp, na
call io_error('kmesh_get: failed to find neighbours for this kpoint')
endif
enddo
enddo
!fill in the global arrays from the local ones
do loop = 1, nntot
wb(loop) = wb_local(loop)
end do
do loop_s = 1, num_kpts
do loop = 1, nntot
bk(:, loop, loop_s) = bk_local(:, loop, loop_s)
end do
end do
![ysl-b]
if (gamma_only) then
! use half of the b-vectors
if (num_kpts .ne. 1) call io_error('Error in kmesh_get: wrong choice of gamma_only option')
! reassign nnlist, nncell, wb, bk
allocate (nnlist_tmp(num_kpts, nntot), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating nnlist_tmp in kmesh_get')
allocate (nncell_tmp(3, num_kpts, nntot), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating nncell_tmp in kmesh_get')
nnlist_tmp(:, :) = nnlist(:, :)
nncell_tmp(:, :, :) = nncell(:, :, :)
deallocate (nnlist, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating nnlist in kmesh_get')
deallocate (nncell, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating nncell in kmesh_get')
deallocate (wb, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating wb in kmesh_get')
deallocate (bk, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating bk in kmesh_get')
nntot = nntot/2
allocate (nnlist(num_kpts, nntot), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating nnlist in kmesh_get')
allocate (nncell(3, num_kpts, nntot), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating nncell in kmesh_get')
allocate (wb(nntot), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating wb in kmesh_get')
allocate (bk(3, nntot, num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating bk in kmesh_get')
na = 0
do nn = 1, 2*nntot
ifound = 0
if (na .ne. 0) then
do nap = 1, na
call utility_compar(bk(1, nap, 1), bk_local(1, nn, 1), ifpos, ifneg)
if (ifneg .eq. 1) ifound = 1
enddo
endif
if (ifound .eq. 0) then
! found new vector to add to set
na = na + 1
bk(1, na, 1) = bk_local(1, nn, 1)
bk(2, na, 1) = bk_local(2, nn, 1)
bk(3, na, 1) = bk_local(3, nn, 1)
wb(na) = 2.0_dp*wb_local(nn)
nnlist(1, na) = nnlist_tmp(1, nn)
nncell(1, 1, na) = nncell_tmp(1, 1, nn)
nncell(2, 1, na) = nncell_tmp(2, 1, nn)
nncell(3, 1, na) = nncell_tmp(3, 1, nn)
neigh(1, na) = na
! check bk.eq.bka
call utility_compar(bk(1, na, 1), bka(1, na), ifpos, ifneg)
if (ifpos .ne. 1) call io_error('Error in kmesh_get: bk is not identical to bka in gamma_only option')
endif
enddo
if (na .ne. nnh) call io_error('Did not find right number of b-vectors in gamma_only option')
if (on_root) then
write (stdout, '(1x,"+",76("-"),"+")')
write (stdout, '(1x,a)') '| Gamma-point: number of the b-vectors is reduced by half |'
write (stdout, '(1x,"+",76("-"),"+")')
if (lenconfac .eq. 1.0_dp) then
write (stdout, '(1x,a)') '| b_k Vectors (Ang^-1) and Weights (Ang^2) |'
write (stdout, '(1x,a)') '| ---------------------------------------- |'
else
write (stdout, '(1x,a)') '| b_k Vectors (Bohr^-1) and Weights (Bohr^2) |'
write (stdout, '(1x,a)') '| ------------------------------------------ |'
endif
write (stdout, '(1x,a)') '| No. b_k(x) b_k(y) b_k(z) w_b |'
write (stdout, '(1x,a)') '| --- -------------------------------- -------- |'
do i = 1, nntot
write (stdout, '(1x,"|",11x,i3,5x,3f12.6,3x,f10.6,8x,"|")') &
i, (bk(j, i, 1)/lenconfac, j=1, 3), wb(i)*lenconfac**2
enddo
write (stdout, '(1x,"+",76("-"),"+")')
write (stdout, *) ' '
endif
deallocate (nnlist_tmp, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating nnlist_tmp in kmesh_get')
deallocate (nncell_tmp, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating nncell_tmp in kmesh_get')
endif
![ysl-e]
if (timing_level > 0) call io_stopwatch('kmesh: get', 2)
return
end subroutine kmesh_get