Find the correct set of shells to satisfy B1 The stratagy is: 1) Take the bvectors from the next shell 2) Reject them if they are parallel to exisiting b vectors 3) Test to see if we satisfy B1, if not add another shell and repeat
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | multi(search_shells) | |||
| real(kind=dp), | intent(in) | :: | dnn(search_shells) | |||
| real(kind=dp), | intent(out) | :: | bweight(max_shells) |
subroutine kmesh_shell_automatic(multi, dnn, bweight)
!==========================================================================
!
!! Find the correct set of shells to satisfy B1
!! The stratagy is:
!! 1) Take the bvectors from the next shell
!! 2) Reject them if they are parallel to exisiting b vectors
!! 3) Test to see if we satisfy B1, if not add another shell and repeat
!
!==========================================================================
use w90_constants, only: eps5, eps6
use w90_io, only: io_error, stdout, io_stopwatch
implicit none
integer, intent(in) :: multi(search_shells) ! the number of kpoints in the shell
real(kind=dp), intent(in) :: dnn(search_shells) ! the bvectors
real(kind=dp), intent(out) :: bweight(max_shells)
real(kind=dp), allocatable :: bvector(:, :, :) ! the bvectors
real(kind=dp), dimension(:), allocatable :: singv, tmp1, tmp2, tmp3
real(kind=dp), dimension(:, :), allocatable :: amat, umat, vmat, smat, tmp0
integer, parameter :: lwork = max_shells*10
real(kind=dp) :: work(lwork)
real(kind=dp), parameter :: target(6) = (/1.0_dp, 1.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/)
logical :: b1sat, lpar
integer :: loop_i, loop_j, loop_bn, loop_b, loop_s, info, cur_shell, ierr
real(kind=dp) :: delta
integer :: loop, shell
if (timing_level > 1) call io_stopwatch('kmesh: shell_automatic', 1)
allocate (bvector(3, maxval(multi), max_shells), stat=ierr)
if (ierr /= 0) call io_error('Error allocating bvector in kmesh_shell_automatic')
bvector = 0.0_dp; bweight = 0.0_dp
if (on_root) write (stdout, '(1x,a)') '| The b-vectors are chosen automatically |'
b1sat = .false.
do shell = 1, search_shells
cur_shell = num_shells + 1
! get the b vectors for the new shell
call kmesh_get_bvectors(multi(shell), 1, dnn(shell), bvector(:, 1:multi(shell), cur_shell))
if (iprint >= 3 .and. on_root) then
write (stdout, '(1x,a8,1x,I2,a14,1x,I2,49x,a)') '| Shell:', shell, ' Multiplicity:', multi(shell), '|'
do loop = 1, multi(shell)
write (stdout, '(1x,a10,I2,1x,a1,4x,3f12.6,5x,a9,9x,a)') '| b-vector ', loop, ':', &
bvector(:, loop, cur_shell)/lenconfac, '('//trim(length_unit)//'^-1)', '|'
end do
end if
! We check that the new shell is not parrallel to an existing shell (cosine=1)
lpar = .false.
if (num_shells > 0) then
do loop_bn = 1, multi(shell)
do loop_s = 1, num_shells
do loop_b = 1, multi(shell_list(loop_s))
delta = dot_product(bvector(:, loop_bn, cur_shell), bvector(:, loop_b, loop_s))/ &
sqrt(dot_product(bvector(:, loop_bn, cur_shell), bvector(:, loop_bn, cur_shell))* &
dot_product(bvector(:, loop_b, loop_s), bvector(:, loop_b, loop_s)))
if (abs(abs(delta) - 1.0_dp) < eps6) lpar = .true.
end do
end do
end do
end if
if (lpar) then
if (iprint >= 3 .and. on_root) then
write (stdout, '(1x,a)') '| This shell is linearly dependent on existing shells: Trying next shell |'
end if
cycle
end if
num_shells = num_shells + 1
shell_list(num_shells) = shell
allocate (tmp0(max_shells, max_shells), stat=ierr)
if (ierr /= 0) call io_error('Error allocating amat in kmesh_shell_automatic')
allocate (tmp1(max_shells), stat=ierr)
if (ierr /= 0) call io_error('Error allocating amat in kmesh_shell_automatic')
allocate (tmp2(num_shells), stat=ierr)
if (ierr /= 0) call io_error('Error allocating amat in kmesh_shell_automatic')
allocate (tmp3(num_shells), stat=ierr)
if (ierr /= 0) call io_error('Error allocating amat in kmesh_shell_automatic')
allocate (amat(max_shells, num_shells), stat=ierr)
if (ierr /= 0) call io_error('Error allocating amat in kmesh_shell_automatic')
allocate (umat(max_shells, max_shells), stat=ierr)
if (ierr /= 0) call io_error('Error allocating umat in kmesh_shell_automatic')
allocate (vmat(num_shells, num_shells), stat=ierr)
if (ierr /= 0) call io_error('Error allocating vmat in kmesh_shell_automatic')
allocate (smat(num_shells, max_shells), stat=ierr)
if (ierr /= 0) call io_error('Error allocating smat in kmesh_shell_automatic')
allocate (singv(num_shells), stat=ierr)
if (ierr /= 0) call io_error('Error allocating singv in kmesh_shell_automatic')
amat(:, :) = 0.0_dp; umat(:, :) = 0.0_dp; vmat(:, :) = 0.0_dp; smat(:, :) = 0.0_dp; singv(:) = 0.0_dp
amat = 0.0_dp
do loop_s = 1, num_shells
do loop_b = 1, multi(shell_list(loop_s))
amat(1, loop_s) = amat(1, loop_s) + bvector(1, loop_b, loop_s)*bvector(1, loop_b, loop_s)
amat(2, loop_s) = amat(2, loop_s) + bvector(2, loop_b, loop_s)*bvector(2, loop_b, loop_s)
amat(3, loop_s) = amat(3, loop_s) + bvector(3, loop_b, loop_s)*bvector(3, loop_b, loop_s)
amat(4, loop_s) = amat(4, loop_s) + bvector(1, loop_b, loop_s)*bvector(2, loop_b, loop_s)
amat(5, loop_s) = amat(5, loop_s) + bvector(2, loop_b, loop_s)*bvector(3, loop_b, loop_s)
amat(6, loop_s) = amat(6, loop_s) + bvector(3, loop_b, loop_s)*bvector(1, loop_b, loop_s)
end do
end do
info = 0
call dgesvd('A', 'A', max_shells, num_shells, amat, max_shells, singv, umat, max_shells, vmat, num_shells, work, lwork, info)
if (info < 0) then
if (on_root) write (stdout, '(1x,a,1x,I1,1x,a)') 'kmesh_shell_automatic: Argument', abs(info), 'of dgesvd is incorrect'
call io_error('kmesh_shell_automatic: Problem with Singular Value Decomposition')
else if (info > 0) then
call io_error('kmesh_shell_automatic: Singular Value Decomposition did not converge')
end if
if (any(abs(singv) < eps5)) then
if (num_shells == 1) then
call io_error('kmesh_shell_automatic: Singular Value Decomposition has found a very small singular value')
else
if (on_root) write (stdout, '(1x,a)') '| SVD found small singular value, Rejecting this shell and trying the next |'
b1sat = .false.
num_shells = num_shells - 1
goto 200
end if
end if
smat = 0.0_dp
do loop_s = 1, num_shells
smat(loop_s, loop_s) = 1.0_dp/singv(loop_s)
end do
! S. Ponce: The following below is correct but had to be unpacked because of PGI-15
! bweight(1:num_shells)=matmul(transpose(vmat),matmul(smat,matmul(transpose(umat),target)))
tmp0 = transpose(umat)
tmp1 = matmul(tmp0, target)
tmp2 = matmul(smat, tmp1)
tmp3 = matmul(transpose(vmat), tmp2)
bweight(1:num_shells) = tmp3
if (iprint >= 2 .and. on_root) then
do loop_s = 1, num_shells
write (stdout, '(1x,a,I2,a,f12.7,5x,a8,36x,a)') '| Shell: ', loop_s, &
' w_b ', bweight(loop_s)*lenconfac**2, '('//trim(length_unit)//'^2)', '|'
end do
end if
!check b1
b1sat = .true.
do loop_i = 1, 3
do loop_j = loop_i, 3
delta = 0.0_dp
do loop_s = 1, num_shells
do loop_b = 1, multi(shell_list(loop_s))
delta = delta + bweight(loop_s)*bvector(loop_i, loop_b, loop_s)*bvector(loop_j, loop_b, loop_s)
end do
end do
if (loop_i == loop_j) then
if (abs(delta - 1.0_dp) > kmesh_tol) b1sat = .false.
end if
if (loop_i /= loop_j) then
if (abs(delta) > kmesh_tol) b1sat = .false.
end if
end do
end do
if (.not. b1sat) then
if (shell < search_shells .and. iprint >= 3) then
if (on_root) write (stdout, '(1x,a,24x,a1)') '| B1 condition is not satisfied: Adding another shell', '|'
elseif (shell == search_shells) then
if (on_root) write (stdout, *) ' '
if (on_root) write (stdout, '(1x,a,i3,a)') 'Unable to satisfy B1 with any of the first ', search_shells, ' shells'
if (on_root) write (stdout, '(1x,a)') 'Check that you have specified your unit cell to a high precision'
if (on_root) write (stdout, '(1x,a)') 'Low precision might cause a loss of symmetry.'
if (on_root) write (stdout, '(1x,a)') ' '
if (on_root) write (stdout, '(1x,a)') 'If your cell is very long, or you have an irregular MP grid'
if (on_root) write (stdout, '(1x,a)') 'Try increasing the parameter search_shells in the win file (default=30)'
if (on_root) write (stdout, *) ' '
call io_error('kmesh_get_automatic')
end if
end if
200 continue
deallocate (tmp0, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating amat in kmesh_shell_automatic')
deallocate (tmp1, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating amat in kmesh_shell_automatic')
deallocate (tmp2, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating amat in kmesh_shell_automatic')
deallocate (tmp3, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating amat in kmesh_shell_automatic')
deallocate (amat, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating amat in kmesh_shell_automatic')
deallocate (umat, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating umat in kmesh_shell_automatic')
deallocate (vmat, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating vmat in kmesh_shell_automatic')
deallocate (smat, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating smat in kmesh_shell_automatic')
deallocate (singv, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating singv in kmesh_shell_automatic')
if (b1sat) exit
end do
if (.not. b1sat) then
if (on_root) write (stdout, *) ' '
if (on_root) write (stdout, '(1x,a,i3,a)') 'Unable to satisfy B1 with any of the first ', search_shells, ' shells'
if (on_root) write (stdout, '(1x,a)') 'Your cell might be very long, or you may have an irregular MP grid'
if (on_root) write (stdout, '(1x,a)') 'Try increasing the parameter search_shells in the win file (default=12)'
if (on_root) write (stdout, *) ' '
call io_error('kmesh_get_automatic')
end if
if (timing_level > 1) call io_stopwatch('kmesh: shell_automatic', 2)
return
end subroutine kmesh_shell_automatic