kmesh_shell_automatic Subroutine

private subroutine kmesh_shell_automatic(kmesh_input, print_output, bweight, dnn, kpt_cart, recip_lattice, lmn, multi, num_kpts, stdout, timer, error, comm)

Uses

  • proc~~kmesh_shell_automatic~~UsesGraph proc~kmesh_shell_automatic kmesh_shell_automatic module~w90_constants w90_constants proc~kmesh_shell_automatic->module~w90_constants module~w90_io w90_io proc~kmesh_shell_automatic->module~w90_io module~w90_types w90_types proc~kmesh_shell_automatic->module~w90_types module~w90_io->module~w90_constants module~w90_types->module~w90_constants

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

Arguments

Type IntentOptional Attributes Name
type(kmesh_input_type), intent(inout) :: kmesh_input
type(print_output_type), intent(in) :: print_output
real(kind=dp), intent(out) :: bweight(kmesh_input%max_shells_h)
real(kind=dp), intent(in) :: dnn(kmesh_input%search_shells)
real(kind=dp), intent(in) :: kpt_cart(:,:)
real(kind=dp), intent(in) :: recip_lattice(3,3)
integer, intent(in) :: lmn(:,:)
integer, intent(in) :: multi(kmesh_input%search_shells)
integer, intent(in) :: num_kpts
integer, intent(in) :: stdout
type(timer_list_type), intent(inout) :: timer
type(w90_error_type), intent(out), allocatable :: error
type(w90_comm_type), intent(in) :: comm

Calls

proc~~kmesh_shell_automatic~~CallsGraph proc~kmesh_shell_automatic kmesh_shell_automatic dgesvd dgesvd proc~kmesh_shell_automatic->dgesvd proc~io_stopwatch_start io_stopwatch_start proc~kmesh_shell_automatic->proc~io_stopwatch_start proc~io_stopwatch_stop io_stopwatch_stop proc~kmesh_shell_automatic->proc~io_stopwatch_stop proc~kmesh_check_condition kmesh_check_condition proc~kmesh_shell_automatic->proc~kmesh_check_condition proc~kmesh_get_amat kmesh_get_amat proc~kmesh_shell_automatic->proc~kmesh_get_amat proc~kmesh_get_bvectors kmesh_get_bvectors proc~kmesh_shell_automatic->proc~kmesh_get_bvectors proc~set_error_alloc set_error_alloc proc~kmesh_shell_automatic->proc~set_error_alloc proc~set_error_dealloc set_error_dealloc proc~kmesh_shell_automatic->proc~set_error_dealloc proc~set_error_fatal set_error_fatal proc~kmesh_shell_automatic->proc~set_error_fatal proc~kmesh_get_bvectors->proc~io_stopwatch_start proc~kmesh_get_bvectors->proc~io_stopwatch_stop proc~kmesh_get_bvectors->proc~set_error_fatal proc~comms_sync_error comms_sync_error proc~set_error_alloc->proc~comms_sync_error proc~set_base_error set_base_error proc~set_error_alloc->proc~set_base_error proc~set_error_dealloc->proc~comms_sync_error proc~set_error_dealloc->proc~set_base_error proc~set_error_fatal->proc~comms_sync_error proc~set_error_fatal->proc~set_base_error

Called by

proc~~kmesh_shell_automatic~~CalledByGraph proc~kmesh_shell_automatic kmesh_shell_automatic proc~kmesh_get kmesh_get proc~kmesh_get->proc~kmesh_shell_automatic proc~w90_create_kmesh w90_create_kmesh proc~w90_create_kmesh->proc~kmesh_get program~postw90 postw90 program~postw90->proc~kmesh_get proc~w90_get_gkpb~2 w90_get_gkpb proc~w90_get_gkpb~2->proc~w90_create_kmesh proc~w90_get_nnkp~2 w90_get_nnkp proc~w90_get_nnkp~2->proc~w90_create_kmesh proc~w90_get_nn~2 w90_get_nn proc~w90_get_nn~2->proc~w90_create_kmesh proc~write_kmesh write_kmesh proc~write_kmesh->proc~w90_create_kmesh proc~w90_get_gkpb w90_get_gkpb proc~w90_get_gkpb->proc~w90_get_gkpb~2 proc~w90_get_nn w90_get_nn proc~w90_get_nn->proc~w90_get_nn~2 proc~w90_get_nnkp w90_get_nnkp proc~w90_get_nnkp->proc~w90_get_nnkp~2 program~wannier wannier program~wannier->proc~w90_get_nn~2 program~wannier->proc~write_kmesh

Source Code

  subroutine kmesh_shell_automatic(kmesh_input, print_output, bweight, dnn, kpt_cart, &
                                   recip_lattice, lmn, multi, num_kpts, stdout, timer, error, &
                                   comm)
    !================================================
    !! 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_stopwatch_start, io_stopwatch_stop
    use w90_types, only: kmesh_input_type, print_output_type, timer_list_type

    implicit none

    ! arguments
    type(print_output_type), intent(in) :: print_output
    type(kmesh_input_type), intent(inout) :: kmesh_input
    type(timer_list_type), intent(inout) :: timer
    type(w90_error_type), allocatable, intent(out) :: error
    type(w90_comm_type), intent(in) :: comm

    integer, intent(in) :: num_kpts
    integer, intent(in) :: stdout
    integer, intent(in) :: lmn(:, :)
    integer, intent(in) :: multi(kmesh_input%search_shells)   ! the number of kpoints in the shell

    real(kind=dp), intent(in) :: recip_lattice(3, 3)
    real(kind=dp), intent(in) :: kpt_cart(:, :)
    real(kind=dp), intent(in) :: dnn(kmesh_input%search_shells) ! the bvectors
    real(kind=dp), intent(out) :: bweight(kmesh_input%max_shells_h)

    ! local variables
    integer :: loop_bn, loop_b, loop_s, info, cur_shell, ierr, loop, shell
    logical :: lpar
    real(kind=dp), allocatable :: amat(:, :), umat(:, :), vmat(:, :), smat(:, :), tmp0(:, :)
    real(kind=dp), allocatable :: bvector(:, :, :) ! the bvectors
    real(kind=dp), allocatable :: singv(:), tmp1(:), tmp2(:), tmp3(:)
    real(kind=dp) :: delta
    real(kind=dp) :: target(kmesh_input%max_shells_aux)
    real(kind=dp) :: work((kmesh_input%max_shells_aux)*10)

    ! variables for higher-order finite-difference
    integer, dimension(:, :), allocatable :: num_x, num_y, num_z
    logical :: bsat
    integer :: loop_order, num_of_eqs, higher_order_n_local

    if (kmesh_input%higher_order_nearest_shells) then
      higher_order_n_local = kmesh_input%higher_order_n
    else
      higher_order_n_local = 1 !find 1st-order b and weights first in this subroutine
    end if
    target = 0.0_dp; target(1) = 1.0_dp; target(3) = 1.0_dp; target(6) = 1.0_dp

    if (print_output%timing_level > 1) call io_stopwatch_start('kmesh: shell_automatic', timer)
    allocate (bvector(3, maxval(multi), kmesh_input%max_shells_h), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error allocating bvector in kmesh_shell_automatic', comm)
      return
    end if
    bvector = 0.0_dp; bweight = 0.0_dp

    if (print_output%iprint > 0) then
      write (stdout, '(1x,a)') '| The b-vectors are chosen automatically                                     |'
    end if

    ! note allocation of kmesh_input%shell list in subroutine w90_readwrite_read_kmesh_data()
    ! kmesh_input%num_shells = 0 in same place
    bsat = .false.

    do shell = 1, kmesh_input%search_shells
      cur_shell = kmesh_input%num_shells + 1

      ! get the b vectors for the new shell
      call kmesh_get_bvectors(kmesh_input, print_output, bvector(:, 1:multi(shell), cur_shell), &
                              kpt_cart, recip_lattice, dnn(shell), lmn, 1, multi(shell), num_kpts, &
                              timer, error, comm)
      if (allocated(error)) return

      if (print_output%iprint >= 3) 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)/print_output%lenconfac, '('//trim(print_output%length_unit)//'^-1)', '|'
        end do
      end if

      ! We check that the new shell is not parrallel to an existing shell (cosine=1)
      if (higher_order_n_local == 1) then
        lpar = .false.
        if (kmesh_input%num_shells > 0) then
          do loop_bn = 1, multi(shell)
            do loop_s = 1, kmesh_input%num_shells
              do loop_b = 1, multi(kmesh_input%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 (print_output%iprint >= 3) then
            write (stdout, '(1x,a)') '| This shell is linearly dependent on existing shells: Trying next shell     |'
          end if
          cycle
        end if
      end if

      kmesh_input%num_shells = kmesh_input%num_shells + 1
      kmesh_input%shell_list(kmesh_input%num_shells) = shell

      allocate (tmp0(kmesh_input%max_shells_aux, kmesh_input%max_shells_aux), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating tmp0 in kmesh_shell_automatic', comm)
        return
      end if
      allocate (tmp1(kmesh_input%max_shells_aux), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating tmp1 in kmesh_shell_automatic', comm)
        return
      end if
      allocate (tmp2(kmesh_input%num_shells), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating tmp2 in kmesh_shell_automatic', comm)
        return
      end if
      allocate (tmp3(kmesh_input%num_shells), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating tmp3 in kmesh_shell_automatic', comm)
        return
      end if
      allocate (amat(kmesh_input%max_shells_aux, kmesh_input%num_shells), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating amat in kmesh_shell_automatic', comm)
        return
      end if
      allocate (umat(kmesh_input%max_shells_aux, kmesh_input%max_shells_aux), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating umat in kmesh_shell_automatic', comm)
        return
      end if
      allocate (vmat(kmesh_input%num_shells, kmesh_input%num_shells), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating vmat in kmesh_shell_automatic', comm)
        return
      end if
      allocate (smat(kmesh_input%num_shells, kmesh_input%max_shells_aux), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating smat in kmesh_shell_automatic', comm)
        return
      end if
      allocate (singv(kmesh_input%num_shells), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating singv in kmesh_shell_automatic', comm)
        return
      end if
      amat(:, :) = 0.0_dp; umat(:, :) = 0.0_dp; vmat(:, :) = 0.0_dp; smat(:, :) = 0.0_dp; singv(:) = 0.0_dp

      num_of_eqs = (1 + higher_order_n_local)*(1 + 2*higher_order_n_local)
      allocate (num_x(higher_order_n_local, num_of_eqs), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating num_x in kmesh_shell_automatic', comm)
        return
      end if
      allocate (num_y(higher_order_n_local, num_of_eqs), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating num_y in kmesh_shell_automatic', comm)
        return
      end if
      allocate (num_z(higher_order_n_local, num_of_eqs), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating num_z in kmesh_shell_automatic', comm)
        return
      end if

      !find higher finite-diff weights
      ! make test suite(compare nnkp files)
      do loop_order = 1, higher_order_n_local
        call kmesh_get_amat(kmesh_input, amat, bvector, multi, loop_order, &
                            num_x(loop_order, :), num_y(loop_order, :), num_z(loop_order, :))
      end do

      info = 0
      call dgesvd('A', 'A', kmesh_input%max_shells_aux, kmesh_input%num_shells, amat, &
                  kmesh_input%max_shells_aux, singv, umat, &
                  kmesh_input%max_shells_aux, vmat, kmesh_input%num_shells, work, kmesh_input%max_shells_aux*10, info)
      if (info < 0) then
        if (print_output%iprint > 0) then
          write (stdout, '(1x,a,1x,I1,1x,a)') 'kmesh_shell_automatic: Argument', abs(info), &
            'of dgesvd is incorrect'
        end if
        call set_error_fatal(error, 'kmesh_shell_automatic: Problem with Singular Value Decomposition', comm)
        return
      else if (info > 0) then
        call set_error_fatal(error, 'kmesh_shell_automatic: Singular Value Decomposition did not converge', comm)
        return
      end if

      if (any(abs(singv) < eps5)) then
        if (kmesh_input%num_shells == 1) then
          call set_error_fatal(error, &
                               'kmesh_shell_automatic: Singular Value Decomposition has found a very small singular value', comm)
          return
        else
          if (print_output%iprint > 0) then
            write (stdout, '(1x,a)') '| SVD found small singular value, Rejecting this shell and trying the next   |'
          end if
          bsat = .false.
          kmesh_input%num_shells = kmesh_input%num_shells - 1
          goto 200
        end if
      end if

      smat = 0.0_dp
      do loop_s = 1, kmesh_input%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:kmesh_input%num_shells) = tmp3

      if (print_output%iprint >= 2) then
        do loop_s = 1, kmesh_input%num_shells
          write (stdout, '(1x,a,I2,a,f12.7,5x,a8,36x,a)') '| Shell: ', loop_s, &
            ' w_b ', bweight(loop_s)*print_output%lenconfac**2, '('//trim(print_output%length_unit)//'^2)', '|'
        end do
      end if

      !check if the conditions including (B1) for finite-difference are satisfied
      bsat = .true.
      do loop_order = 1, higher_order_n_local
        call kmesh_check_condition(kmesh_input, bsat, bvector, bweight, multi, loop_order, &
                                   num_x(loop_order, :), num_y(loop_order, :), num_z(loop_order, :))
      end do

      if (.not. bsat) then
        if (shell < kmesh_input%search_shells .and. print_output%iprint >= 3) then
          if (print_output%iprint > 0) write (stdout, '(1x,a,24x,a1)') '| B1 condition is not satisfied: Adding another shell', '|'

        elseif (shell == kmesh_input%search_shells) then

          if (print_output%iprint > 0) then
            write (stdout, *) ' '
            write (stdout, '(1x,a,i3,a)') 'Unable to satisfy the higher-order version of B1 with any of the first ' &
              , kmesh_input%search_shells, ' shells'
            write (stdout, '(1x,a)') 'Check that you have specified your unit cell to a high precision'
            write (stdout, '(1x,a)') 'Low precision might cause a loss of symmetry.'
            write (stdout, '(1x,a)') ' '
            write (stdout, '(1x,a)') 'If your cell is very long, or you have an irregular MP grid'
            write (stdout, '(1x,a)') 'Try increasing the parameter search_shells in the win file (default=30)'
            write (stdout, *) ' '
            call set_error_fatal(error, 'kmesh_shell_automatic: unable to satisfy the higher-order version of B1', comm)
            return
          end if

        end if
      end if

200   continue

      deallocate (tmp0, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating tmp0 in kmesh_shell_automatic', comm)
        return
      end if
      deallocate (tmp1, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating tmp1 in kmesh_shell_automatic', comm)
        return
      end if
      deallocate (tmp2, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating tmp2 in kmesh_shell_automatic', comm)
        return
      end if
      deallocate (tmp3, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating tmp3 in kmesh_shell_automatic', comm)
        return
      end if
      deallocate (amat, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating amat in kmesh_shell_automatic', comm)
        return
      end if
      deallocate (umat, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating umat in kmesh_shell_automatic', comm)
        return
      end if
      deallocate (vmat, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating vmat in kmesh_shell_automatic', comm)
        return
      end if
      deallocate (smat, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating smat in kmesh_shell_automatic', comm)
        return
      end if
      deallocate (singv, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating singv in kmesh_shell_automatic', comm)
        return
      end if

      deallocate (num_x, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating num_x in kmesh_shell_automatic', comm)
        return
      end if
      deallocate (num_y, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating num_y in kmesh_shell_automatic', comm)
        return
      end if
      deallocate (num_z, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating num_z in kmesh_shell_automatic', comm)
        return
      end if

      if (bsat) exit

    end do

    if (.not. bsat) then
      if (print_output%iprint > 0) then
        write (stdout, *) ' '
        write (stdout, '(1x,a,i3,a)') 'Unable to satisfy B1 with any of the first ', &
          kmesh_input%search_shells, ' shells'
        write (stdout, '(1x,a)') 'Your cell might be very long, or you may have an irregular MP grid'
        write (stdout, '(1x,a)') 'Try increasing the parameter search_shells in the win file (default=36)'
        write (stdout, *) ' '
      end if
      call set_error_fatal(error, 'kmesh_shell_automatic: Unable to satisfy B1 condition', comm)
      return
    end if

    if (print_output%timing_level > 1) call io_stopwatch_stop('kmesh: shell_automatic', timer)

    return

  end subroutine kmesh_shell_automatic