hamiltonian_wigner_seitz Subroutine

public subroutine hamiltonian_wigner_seitz(ws_region, print_output, real_lattice, irvec, mp_grid, ndegen, nrpts, rpt_origin, stdout, timer, error, count_pts, comm)

Uses

  • proc~~hamiltonian_wigner_seitz~~UsesGraph proc~hamiltonian_wigner_seitz hamiltonian_wigner_seitz module~w90_constants w90_constants proc~hamiltonian_wigner_seitz->module~w90_constants module~w90_io w90_io proc~hamiltonian_wigner_seitz->module~w90_io module~w90_types w90_types proc~hamiltonian_wigner_seitz->module~w90_types module~w90_utility w90_utility proc~hamiltonian_wigner_seitz->module~w90_utility module~w90_io->module~w90_constants module~w90_types->module~w90_constants module~w90_utility->module~w90_constants module~w90_comms w90_comms module~w90_utility->module~w90_comms module~w90_comms->module~w90_constants module~w90_error_base w90_error_base module~w90_comms->module~w90_error_base

Calculates a grid of points that fall inside of (and eventually on the surface of) the Wigner-Seitz supercell centered on the origin of the B lattice with primitive translations nmonkh(1)a_1+nmonkh(2)a_2+nmonkh(3)*a_3

Arguments

Type IntentOptional Attributes Name
type(ws_region_type), intent(in) :: ws_region
type(print_output_type), intent(in) :: print_output
real(kind=dp), intent(in) :: real_lattice(3,3)
integer, intent(inout), allocatable :: irvec(:,:)
integer, intent(in) :: mp_grid(3)
integer, intent(inout), allocatable :: ndegen(:)
integer, intent(inout) :: nrpts
integer, intent(inout) :: rpt_origin
integer, intent(in) :: stdout
type(timer_list_type), intent(inout) :: timer
type(w90_error_type), intent(out), allocatable :: error
logical, intent(in) :: count_pts
type(w90_comm_type), intent(in) :: comm

Calls

proc~~hamiltonian_wigner_seitz~~CallsGraph proc~hamiltonian_wigner_seitz hamiltonian_wigner_seitz proc~io_stopwatch_start io_stopwatch_start proc~hamiltonian_wigner_seitz->proc~io_stopwatch_start proc~io_stopwatch_stop io_stopwatch_stop proc~hamiltonian_wigner_seitz->proc~io_stopwatch_stop proc~set_error_alloc set_error_alloc proc~hamiltonian_wigner_seitz->proc~set_error_alloc proc~set_error_dealloc set_error_dealloc proc~hamiltonian_wigner_seitz->proc~set_error_dealloc proc~set_error_fatal set_error_fatal proc~hamiltonian_wigner_seitz->proc~set_error_fatal proc~utility_metric utility_metric proc~hamiltonian_wigner_seitz->proc~utility_metric 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~~hamiltonian_wigner_seitz~~CalledByGraph proc~hamiltonian_wigner_seitz hamiltonian_wigner_seitz proc~hamiltonian_setup hamiltonian_setup proc~hamiltonian_setup->proc~hamiltonian_wigner_seitz proc~plot_main plot_main proc~plot_main->proc~hamiltonian_setup proc~tran_main tran_main proc~tran_main->proc~hamiltonian_setup proc~wann_main wann_main proc~wann_main->proc~hamiltonian_setup proc~w90_plot w90_plot proc~w90_plot->proc~plot_main proc~w90_transport w90_transport proc~w90_transport->proc~tran_main proc~w90_wannierise~2 w90_wannierise proc~w90_wannierise~2->proc~wann_main proc~w90_wannierise w90_wannierise proc~w90_wannierise->proc~w90_wannierise~2 program~wannier wannier program~wannier->proc~w90_plot program~wannier->proc~w90_transport program~wannier->proc~w90_wannierise~2

Source Code

  subroutine hamiltonian_wigner_seitz(ws_region, print_output, real_lattice, irvec, mp_grid, &
                                      ndegen, nrpts, rpt_origin, stdout, timer, error, count_pts, &
                                      comm)
    !================================================!
    !! Calculates a grid of points that fall inside of (and eventually on the
    !! surface of) the Wigner-Seitz supercell centered on the origin of the B
    !! lattice with primitive translations nmonkh(1)*a_1+nmonkh(2)*a_2+nmonkh(3)*a_3
    !================================================!

    use w90_constants, only: eps8
    use w90_io, only: io_stopwatch_start, io_stopwatch_stop
    use w90_utility, only: utility_metric
    use w90_types, only: print_output_type, ws_region_type, timer_list_type

    ! irvec(i,irpt)     The irpt-th Wigner-Seitz grid point has components
    !                   irvec(1:3,irpt) in the basis of the lattice vectors
    ! ndegen(irpt)      Weight of the irpt-th point is 1/ndegen(irpt)
    ! nrpts             number of Wigner-Seitz grid points

    implicit none

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

    integer, intent(inout)              :: nrpts
    integer, intent(inout), allocatable :: ndegen(:)
    integer, intent(inout), allocatable :: irvec(:, :)
    integer, intent(inout)              :: rpt_origin
    integer, intent(in)                 :: mp_grid(3)
    integer, intent(in)                 :: stdout

    real(kind=dp), intent(in)           :: real_lattice(3, 3)

    logical, intent(in)                 :: count_pts

    ! local variables
    integer       :: ndiff(3)
    integer       :: n1, n2, n3, i1, i2, i3, icnt, i, j, ierr, dist_dim
    real(kind=dp)              :: tot, dist_min
    real(kind=dp), allocatable :: dist(:)
    real(kind=dp)              :: real_metric(3, 3)

    if (print_output%timing_level > 1) &
      call io_stopwatch_start('hamiltonian: wigner_seitz', timer)

    call utility_metric(real_lattice, real_metric)
    dist_dim = 1
    do i = 1, 3
      dist_dim = dist_dim*((ws_region%ws_search_size(i) + 1)*2 + 1)
    end do
    allocate (dist(dist_dim), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating dist in hamiltonian_wigner_seitz', comm)
      return
    end if

    ! The Wannier functions live in a supercell of the real space unit cell
    ! this supercell is mp_grid unit cells long in each direction
    !
    ! We loop over grid points r on a unit cell that is (2*ws_search_size+1)**3 times
    ! larger than this primitive supercell.
    !
    ! One of these points is in the W-S cell if it is closer to R=0 than any of the
    ! other points, R (where R are the translation vectors of the supercell)

    ! In the end nrpts contains the total number of grid
    ! points that have been found in the Wigner-Seitz cell

    nrpts = 0
    ! Loop over the lattice vectors of the primitive cell
    ! that live in a supercell which is (2*ws_search_size+1)**2
    ! larger than the Born-von Karman supercell.
    ! We need to find which among these live in the Wigner-Seitz cell
    do n1 = -ws_region%ws_search_size(1)*mp_grid(1), ws_region%ws_search_size(1)*mp_grid(1)
      do n2 = -ws_region%ws_search_size(2)*mp_grid(2), ws_region%ws_search_size(2)*mp_grid(2)
        do n3 = -ws_region%ws_search_size(3)*mp_grid(3), ws_region%ws_search_size(3)*mp_grid(3)
          ! Loop over the lattice vectors R of the Born-von Karman supercell
          ! that contains all the points of the previous loop.
          ! There are (2*(ws_search_size+1)+1)**3 points R. R=0 corresponds to
          ! i1=i2=i3=0, or icnt=((2*(ws_search_size+1)+1)**3 + 1)/2
          icnt = 0
          do i1 = -ws_region%ws_search_size(1) - 1, ws_region%ws_search_size(1) + 1
            do i2 = -ws_region%ws_search_size(2) - 1, ws_region%ws_search_size(2) + 1
              do i3 = -ws_region%ws_search_size(3) - 1, ws_region%ws_search_size(3) + 1
                icnt = icnt + 1
                ! Calculate distance squared |r-R|^2
                ndiff(1) = n1 - i1*mp_grid(1)
                ndiff(2) = n2 - i2*mp_grid(2)
                ndiff(3) = n3 - i3*mp_grid(3)
                dist(icnt) = 0.0_dp
                do i = 1, 3
                  do j = 1, 3
                    dist(icnt) = dist(icnt) + real(ndiff(i), dp)*real_metric(i, j) &
                                 *real(ndiff(j), dp)
                  end do
                end do
              end do
            end do
          end do
          ! AAM: On first pass, we reference unallocated variables (ndegen,irvec)
          dist_min = minval(dist)
          if (abs(dist((dist_dim + 1)/2) - dist_min) .lt. ws_region%ws_distance_tol**2) then
            nrpts = nrpts + 1
            if (.not. count_pts) then
              ndegen(nrpts) = 0
              do i = 1, dist_dim
                if (abs(dist(i) - dist_min) .lt. ws_region%ws_distance_tol**2) &
                  ndegen(nrpts) = ndegen(nrpts) + 1
              end do
              irvec(1, nrpts) = n1
              irvec(2, nrpts) = n2
              irvec(3, nrpts) = n3
              !
              ! Record index of r=0
              if (n1 == 0 .and. n2 == 0 .and. n3 == 0) rpt_origin = nrpts
            end if
          end if

          !n3
        end do
        !n2
      end do
      !n1
    end do
    !
    deallocate (dist, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating dist hamiltonian_wigner_seitz', comm)
      return
    end if
    if (count_pts) then
      if (print_output%timing_level > 1) &
        call io_stopwatch_stop('hamiltonian: wigner_seitz', timer)
      return
    end if

    ! Check the "sum rule"
    tot = 0.0_dp
    do i = 1, nrpts
      tot = tot + 1.0_dp/real(ndegen(i), dp)
    end do

    if (print_output%iprint >= 3) then
      write (stdout, '(1x,i4,a,/)') nrpts, ' lattice points in Wigner-Seitz supercell:'
      do i = 1, nrpts
        write (stdout, '(4x,a,3(i3,1x),a,i2)') '  vector ', irvec(1, i), irvec(2, i), &
          irvec(3, i), '  degeneracy: ', ndegen(i)
      end do
      write (stdout, '(1x,a,f12.3)') ' tot = ', tot
      write (stdout, '(1x,a,i12)') ' mp_grid product = ', mp_grid(1)*mp_grid(2)*mp_grid(3)
    end if
    if (abs(tot - real(mp_grid(1)*mp_grid(2)*mp_grid(3), dp)) > eps8) then
      call set_error_fatal(error, 'ERROR in hamiltonian_wigner_seitz: error in finding Wigner-Seitz points', comm)
      return
    end if

    if (print_output%timing_level > 1) call io_stopwatch_stop('hamiltonian: wigner_seitz', timer)

    return

  end subroutine hamiltonian_wigner_seitz