Find the B1 weights for a set of b-vectors given in a file. This routine is only activated via a devel_flag and is not intended for regular use.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(inout) | :: | multi(search_shells) | |||
| real(kind=dp), | intent(in) | :: | dnn(search_shells) | |||
| real(kind=dp), | intent(out) | :: | bweight(max_shells) |
subroutine kmesh_shell_from_file(multi, dnn, bweight)
!=================================================================
!
!! Find the B1 weights for a set of b-vectors given in a file.
!! This routine is only activated via a devel_flag and is not
!! intended for regular use.
!
!=================================================================
use w90_constants, only: eps7
use w90_io, only: io_error, stdout, io_stopwatch, io_file_unit, seedname, maxlen
implicit none
integer, intent(inout) :: 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(:, :)
real(kind=dp), dimension(:), allocatable :: singv
real(kind=dp), dimension(:, :), allocatable :: amat, umat, vmat, smat
integer, parameter :: lwork = max_shells*10
real(kind=dp) :: work(lwork)
integer :: bvec_list(num_nnmax, max_shells)
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
integer :: ierr, loop_i, loop_j, loop_b, loop_s, info
real(kind=dp) :: delta
integer :: loop, shell, pos, kshell_in, counter, length, i, loop2, num_lines, tot_num_lines
character(len=maxlen) :: dummy, dummy2
if (timing_level > 1) call io_stopwatch('kmesh: shell_fixed', 1)
allocate (bvector(3, sum(multi)), stat=ierr)
if (ierr /= 0) call io_error('Error allocating bvector in kmesh_shell_fixed')
bvector = 0.0_dp; bweight = 0.0_dp
if (on_root) write (stdout, '(1x,a)') '| The b-vectors are defined in the kshell file |'
counter = 1
do shell = 1, search_shells
! get the b vectors
call kmesh_get_bvectors(multi(shell), 1, dnn(shell), &
bvector(:, counter:counter + multi(shell) - 1))
counter = counter + multi(shell)
end do
kshell_in = io_file_unit()
open (unit=kshell_in, file=trim(seedname)//'.kshell', &
form='formatted', status='old', action='read', err=101)
num_lines = 0; tot_num_lines = 0
do
read (kshell_in, '(a)', iostat=ierr, err=200, end=210) dummy
dummy = adjustl(dummy)
tot_num_lines = tot_num_lines + 1
if (.not. dummy(1:1) == '!' .and. .not. dummy(1:1) == '#') then
if (len(trim(dummy)) > 0) num_lines = num_lines + 1
endif
end do
200 call io_error('Error: Problem (1) reading input file '//trim(seedname)//'.kshell')
210 continue
rewind (kshell_in)
num_shells = num_lines
multi(:) = 0
bvec_list = 1
counter = 0
do loop = 1, tot_num_lines
read (kshell_in, '(a)', err=103, end=103) dummy2
dummy2 = adjustl(dummy2)
if (dummy2(1:1) == '!' .or. dummy2(1:1) == '#' .or. (len(trim(dummy2)) == 0)) cycle
counter = counter + 1
shell_list(counter) = counter
dummy = dummy2
length = 1
dummy = adjustl(dummy)
do
pos = index(dummy, ' ')
dummy = dummy(pos + 1:)
dummy = adjustl(dummy)
if (len_trim(dummy) > 0) then
length = length + 1
else
exit
endif
end do
multi(counter) = length
read (dummy2, *, err=230, end=230) (bvec_list(i, loop), i=1, length)
end do
bvec_inp = 0.0_dp
do loop = 1, num_shells
do loop2 = 1, multi(loop)
bvec_inp(:, loop2, loop) = bvector(:, bvec_list(loop2, loop))
end do
end do
if (iprint >= 3 .and. on_root) then
do shell = 1, num_shells
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, ':', &
bvec_inp(:, loop, shell)/lenconfac, '('//trim(length_unit)//'^-1)', '|'
end do
end do
end if
allocate (amat(max_shells, num_shells), stat=ierr)
if (ierr /= 0) call io_error('Error allocating amat in kmesh_shell_from_file')
allocate (umat(max_shells, max_shells), stat=ierr)
if (ierr /= 0) call io_error('Error allocating umat in kmesh_shell_from_file')
allocate (vmat(num_shells, num_shells), stat=ierr)
if (ierr /= 0) call io_error('Error allocating vmat in kmesh_shell_from_file')
allocate (smat(num_shells, max_shells), stat=ierr)
if (ierr /= 0) call io_error('Error allocating smat in kmesh_shell_from_file')
allocate (singv(num_shells), stat=ierr)
if (ierr /= 0) call io_error('Error allocating singv in kmesh_shell_from_file')
amat = 0.0_dp; umat = 0.0_dp; vmat = 0.0_dp; smat = 0.0_dp; singv = 0.0_dp
do loop_s = 1, num_shells
do loop_b = 1, multi(loop_s)
amat(1, loop_s) = amat(1, loop_s) + bvec_inp(1, loop_b, loop_s)*bvec_inp(1, loop_b, loop_s)
amat(2, loop_s) = amat(2, loop_s) + bvec_inp(2, loop_b, loop_s)*bvec_inp(2, loop_b, loop_s)
amat(3, loop_s) = amat(3, loop_s) + bvec_inp(3, loop_b, loop_s)*bvec_inp(3, loop_b, loop_s)
amat(4, loop_s) = amat(4, loop_s) + bvec_inp(1, loop_b, loop_s)*bvec_inp(2, loop_b, loop_s)
amat(5, loop_s) = amat(5, loop_s) + bvec_inp(2, loop_b, loop_s)*bvec_inp(3, loop_b, loop_s)
amat(6, loop_s) = amat(6, loop_s) + bvec_inp(3, loop_b, loop_s)*bvec_inp(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_fixed: Argument', abs(info), 'of dgesvd is incorrect'
call io_error('kmesh_shell_fixed: Problem with Singular Value Decomposition')
else if (info > 0) then
call io_error('kmesh_shell_fixed: Singular Value Decomposition did not converge')
end if
if (any(abs(singv) < eps7)) &
call io_error('kmesh_shell_fixed: Singular Value Decomposition has found a very small singular value')
smat = 0.0_dp
do loop_s = 1, num_shells
smat(loop_s, loop_s) = 1/singv(loop_s)
end do
bweight(1:num_shells) = matmul(transpose(vmat), matmul(smat, matmul(transpose(umat), target)))
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.
if (.not. skip_B1_tests) then
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(loop_s)
delta = delta + bweight(loop_s)*bvec_inp(loop_i, loop_b, loop_s)*bvec_inp(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
end if
if (.not. b1sat) call io_error('kmesh_shell_fixed: B1 condition not satisfied')
if (timing_level > 1) call io_stopwatch('kmesh: shell_fixed', 2)
return
101 call io_error('Error: Problem (2) reading input file '//trim(seedname)//'.kshell')
103 call io_error('Error: Problem (3) reading input file '//trim(seedname)//'.kshell')
230 call io_error('Error: Problem reading in param_get_keyword_vector')
end subroutine kmesh_shell_from_file