kmesh_shell_fixed Subroutine

  subroutine kmesh_shell_fixed(multi, dnn, bweight)
    !================================================================
    !
    !!  Find the B1 weights for a set of shells specified by the user
    !
    !================================================================

    use w90_constants, only: eps7
    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(:, :, :)

    real(kind=dp) :: singv(num_shells)
    real(kind=dp) :: amat(max_shells, num_shells)
    real(kind=dp) :: umat(max_shells, max_shells)
    real(kind=dp) :: vmat(num_shells, num_shells)
    real(kind=dp) :: smat(num_shells, max_shells)
    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
    integer :: ierr, loop_i, loop_j, loop_b, loop_s, info
    real(kind=dp) :: delta

    integer :: loop, shell

    if (timing_level > 1) call io_stopwatch('kmesh: shell_fixed', 1)

    allocate (bvector(3, maxval(multi), num_shells), stat=ierr)
    if (ierr /= 0) call io_error('Error allocating bvector in kmesh_shell_fixed')
    bvector = 0.0_dp; bweight = 0.0_dp
    amat = 0.0_dp; umat = 0.0_dp; vmat = 0.0_dp; smat = 0.0_dp; singv = 0.0_dp

    if (on_root) write (stdout, '(1x,a)') '| The b-vectors are set in the win file                                      |'

    do shell = 1, num_shells
      ! get the b vectors for this shell
      call kmesh_get_bvectors(multi(shell_list(shell)), 1, dnn(shell_list(shell)), &
                              bvector(:, 1:multi(shell_list(shell)), shell))
    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_list(shell)), '|'
        do loop = 1, multi(shell_list(shell))
          write (stdout, '(1x,a10,I2,1x,a1,4x,3f12.6,5x,a9,9x,a)') '| b-vector ', loop, ':', &
            bvector(:, loop, shell)/lenconfac, '('//trim(length_unit)//'^-1)', '|'
        end do
      end do
    end if

    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_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,49x,a)') '| Shell: ',loop_s,' w_b ', bweight(loop_s),'|'
        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(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
    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

  end subroutine kmesh_shell_fixed