pw90common_wanint_setup Subroutine

  subroutine pw90common_wanint_setup(num_wann, print_output, real_lattice, mp_grid, ws_region, &
                                     ws_distance, effective_model, wigner_seitz, wannier_data, &
                                     stdout, seedname, timer, error, comm)
    !================================================!
    !
    !! Setup data ready for interpolation
    !
    !================================================!

    use w90_constants, only: dp
    use w90_io, only: io_stopwatch_start, io_stopwatch_stop
    use w90_types, only: ws_region_type, ws_distance_type, wannier_data_type, &
                         print_output_type, timer_list_type
    use w90_comms, only: mpirank, w90_comm_type, comms_bcast
    use w90_postw90_types, only: wigner_seitz_type
    use w90_ws_distance, only: ws_translate_dist

    type(print_output_type), intent(in) :: print_output
    type(wigner_seitz_type), intent(inout) :: wigner_seitz
    type(ws_distance_type), intent(inout) :: ws_distance
    type(ws_region_type), intent(in) :: ws_region
    type(wannier_data_type), intent(in) :: wannier_data
    type(timer_list_type), intent(inout) :: timer
    type(w90_comm_type), intent(in) :: comm
    type(w90_error_type), allocatable, intent(out) :: error

    real(kind=dp), intent(in) :: real_lattice(3, 3)
    integer, intent(in) :: num_wann
    integer, intent(in) :: stdout
    integer, intent(in) :: mp_grid(3)
    logical, intent(in) :: effective_model
    character(len=50), intent(in)  :: seedname

    integer :: ierr, ir, file_unit, num_wann_loc
    logical :: on_root = .false.

    if (mpirank(comm) == 0) on_root = .true.

    ! Find nrpts, the number of points in the Wigner-Seitz cell
    if (effective_model) then
      if (on_root) then
        ! nrpts is read from file, together with num_wann
        open (newunit=file_unit, file=trim(seedname)//'_HH_R.dat', form='formatted', &
              status='old', err=101)
        read (file_unit, *) !header
        read (file_unit, *) num_wann_loc
        if (num_wann_loc /= num_wann) then
          call set_error_fatal(error, 'Inconsistent values of num_wann in '//trim(seedname) &
                               //'_HH_R.dat and '//trim(seedname)//'.win', comm)
          return
        end if
        read (file_unit, *) wigner_seitz%nrpts
        close (file_unit)
      end if
      call comms_bcast(wigner_seitz%nrpts, 1, error, comm)
      if (allocated(error)) return
    else
      call wignerseitz(print_output, real_lattice, mp_grid, ws_region, wigner_seitz, stdout, &
                       .true., timer, error, comm)
      if (allocated(error)) return
    end if

    ! Now can allocate several arrays
    allocate (wigner_seitz%irvec(3, wigner_seitz%nrpts), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating irvec in pw90common_wanint_setup', comm)
      return
    end if
    wigner_seitz%irvec = 0
    allocate (wigner_seitz%crvec(3, wigner_seitz%nrpts), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating crvec in pw90common_wanint_setup', comm)
      return
    end if
    wigner_seitz%crvec = 0.0_dp
    allocate (wigner_seitz%ndegen(wigner_seitz%nrpts), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating ndegen in pw90common_wanint_setup', comm)
      return
    end if
    wigner_seitz%ndegen = 0

    ! Also rpt_origin, so that when effective_model=.true it is not
    ! passed to get_HH_R without being initialized.
    wigner_seitz%rpt_origin = 0

    ! If effective_model, this is done in get_HH_R
    if (.not. effective_model) then
      ! Set up the lattice vectors on the Wigner-Seitz supercell
      ! where the Wannier functions live

      call wignerseitz(print_output, real_lattice, mp_grid, ws_region, wigner_seitz, stdout, &
                       .false., timer, error, comm)
      if (allocated(error)) return

      ! Convert from reduced to Cartesian coordinates

      do ir = 1, wigner_seitz%nrpts
        ! Note that 'real_lattice' stores the lattice vectors as *rows*
        wigner_seitz%crvec(:, ir) = matmul(transpose(real_lattice), wigner_seitz%irvec(:, ir))
      end do

      if (ws_region%use_ws_distance) then
        call ws_translate_dist(ws_distance, ws_region, num_wann, wannier_data%centres, real_lattice, &
                               mp_grid, wigner_seitz%nrpts, wigner_seitz%irvec, error, comm)
        if (allocated(error)) return
      end if

      call wigner_seitz_opt_setup(ws_distance, ws_region, wigner_seitz, num_wann, real_lattice, &
                                  wigner_seitz%nrpts, error, comm)

    end if

    return

101 call set_error_file(error, 'Error in pw90common_wanint_setup: problem opening file '// &
                        trim(seedname)//'_HH_R.dat', comm)
    return

  end subroutine pw90common_wanint_setup