geninterp_main Subroutine

public subroutine geninterp_main(dis_manifold, pw90_geninterp, kpt_latt, pw90_band_deriv_degen, ws_region, print_output, wannier_data, ws_distance, wigner_seitz, HH_R, v_matrix, u_matrix, eigval, real_lattice, scissors_shift, mp_grid, num_bands, num_kpts, num_wann, num_valence_bands, effective_model, have_disentangled, seedname, stdout, timer, error, comm)

Uses

  • proc~~geninterp_main~~UsesGraph proc~geninterp_main geninterp_main module~w90_comms w90_comms proc~geninterp_main->module~w90_comms module~w90_constants w90_constants proc~geninterp_main->module~w90_constants module~w90_get_oper w90_get_oper proc~geninterp_main->module~w90_get_oper module~w90_io w90_io proc~geninterp_main->module~w90_io module~w90_postw90_common w90_postw90_common proc~geninterp_main->module~w90_postw90_common module~w90_postw90_types w90_postw90_types proc~geninterp_main->module~w90_postw90_types module~w90_types w90_types proc~geninterp_main->module~w90_types module~w90_utility w90_utility proc~geninterp_main->module~w90_utility module~w90_wan_ham w90_wan_ham proc~geninterp_main->module~w90_wan_ham module~w90_comms->module~w90_constants module~w90_error_base w90_error_base module~w90_comms->module~w90_error_base module~w90_get_oper->module~w90_comms module~w90_get_oper->module~w90_constants module~w90_get_oper->module~w90_io module~w90_error w90_error module~w90_get_oper->module~w90_error module~w90_io->module~w90_constants module~w90_postw90_common->module~w90_constants module~w90_postw90_common->module~w90_error module~w90_postw90_types->module~w90_comms module~w90_postw90_types->module~w90_constants module~w90_types->module~w90_constants module~w90_utility->module~w90_comms module~w90_utility->module~w90_constants module~w90_wan_ham->module~w90_constants module~w90_wan_ham->module~w90_error module~w90_error->module~w90_comms module~w90_error->module~w90_error_base

This routine prints the band energies (and possibly the band derivatives)

This routine is parallel, even if the scaling is very bad since at the moment everything must be written by the root node (we need that the output is sorted). But at least if works independently of the number of processors. I think that a way to write in parallel to the output would help a lot, so that we don't have to send all eigenvalues to the root node.

Arguments

Type IntentOptional Attributes Name
type(dis_manifold_type), intent(in) :: dis_manifold
type(pw90_geninterp_mod_type), intent(in) :: pw90_geninterp
real(kind=dp), intent(in) :: kpt_latt(:,:)
type(pw90_band_deriv_degen_type), intent(in) :: pw90_band_deriv_degen
type(ws_region_type), intent(in) :: ws_region
type(print_output_type), intent(in) :: print_output
type(wannier_data_type), intent(in) :: wannier_data
type(ws_distance_type), intent(inout) :: ws_distance
type(wigner_seitz_type), intent(inout) :: wigner_seitz
complex(kind=dp), intent(inout), allocatable :: HH_R(:,:,:)
complex(kind=dp), intent(in) :: v_matrix(:,:,:)
complex(kind=dp), intent(in) :: u_matrix(:,:,:)
real(kind=dp), intent(in) :: eigval(:,:)
real(kind=dp), intent(in) :: real_lattice(3,3)
real(kind=dp), intent(in) :: scissors_shift
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: num_bands
integer, intent(in) :: num_kpts
integer, intent(in) :: num_wann
integer, intent(in) :: num_valence_bands
logical, intent(in) :: effective_model
logical, intent(in) :: have_disentangled
character(len=50), intent(in) :: seedname
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~~geninterp_main~~CallsGraph proc~geninterp_main geninterp_main interface~comms_bcast comms_bcast proc~geninterp_main->interface~comms_bcast interface~comms_gatherv comms_gatherv proc~geninterp_main->interface~comms_gatherv interface~comms_scatterv comms_scatterv proc~geninterp_main->interface~comms_scatterv proc~comms_array_split comms_array_split proc~geninterp_main->proc~comms_array_split proc~get_hh_r get_HH_R proc~geninterp_main->proc~get_hh_r proc~internal_write_header internal_write_header proc~geninterp_main->proc~internal_write_header proc~io_stopwatch_start io_stopwatch_start proc~geninterp_main->proc~io_stopwatch_start proc~io_stopwatch_stop io_stopwatch_stop proc~geninterp_main->proc~io_stopwatch_stop proc~mpirank mpirank proc~geninterp_main->proc~mpirank proc~mpisize mpisize proc~geninterp_main->proc~mpisize proc~pw90common_fourier_r_to_k pw90common_fourier_R_to_k proc~geninterp_main->proc~pw90common_fourier_r_to_k proc~set_error_alloc set_error_alloc proc~geninterp_main->proc~set_error_alloc proc~set_error_file set_error_file proc~geninterp_main->proc~set_error_file proc~set_error_input set_error_input proc~geninterp_main->proc~set_error_input proc~utility_diagonalize utility_diagonalize proc~geninterp_main->proc~utility_diagonalize proc~utility_recip_lattice_base utility_recip_lattice_base proc~geninterp_main->proc~utility_recip_lattice_base proc~wham_get_eig_deleig wham_get_eig_deleig proc~geninterp_main->proc~wham_get_eig_deleig proc~comms_bcast_char comms_bcast_char interface~comms_bcast->proc~comms_bcast_char proc~comms_bcast_cmplx comms_bcast_cmplx interface~comms_bcast->proc~comms_bcast_cmplx proc~comms_bcast_int comms_bcast_int interface~comms_bcast->proc~comms_bcast_int proc~comms_bcast_logical comms_bcast_logical interface~comms_bcast->proc~comms_bcast_logical proc~comms_bcast_real comms_bcast_real interface~comms_bcast->proc~comms_bcast_real proc~comms_gatherv_cmplx_1 comms_gatherv_cmplx_1 interface~comms_gatherv->proc~comms_gatherv_cmplx_1 proc~comms_gatherv_cmplx_2 comms_gatherv_cmplx_2 interface~comms_gatherv->proc~comms_gatherv_cmplx_2 proc~comms_gatherv_cmplx_3 comms_gatherv_cmplx_3 interface~comms_gatherv->proc~comms_gatherv_cmplx_3 proc~comms_gatherv_cmplx_3_4 comms_gatherv_cmplx_3_4 interface~comms_gatherv->proc~comms_gatherv_cmplx_3_4 proc~comms_gatherv_cmplx_4 comms_gatherv_cmplx_4 interface~comms_gatherv->proc~comms_gatherv_cmplx_4 proc~comms_gatherv_logical comms_gatherv_logical interface~comms_gatherv->proc~comms_gatherv_logical proc~comms_gatherv_real_1 comms_gatherv_real_1 interface~comms_gatherv->proc~comms_gatherv_real_1 proc~comms_gatherv_real_2 comms_gatherv_real_2 interface~comms_gatherv->proc~comms_gatherv_real_2 proc~comms_gatherv_real_2_3 comms_gatherv_real_2_3 interface~comms_gatherv->proc~comms_gatherv_real_2_3 proc~comms_gatherv_real_3 comms_gatherv_real_3 interface~comms_gatherv->proc~comms_gatherv_real_3 proc~comms_scatterv_cmplx_3 comms_scatterv_cmplx_3 interface~comms_scatterv->proc~comms_scatterv_cmplx_3 proc~comms_scatterv_cmplx_4 comms_scatterv_cmplx_4 interface~comms_scatterv->proc~comms_scatterv_cmplx_4 proc~comms_scatterv_int_1 comms_scatterv_int_1 interface~comms_scatterv->proc~comms_scatterv_int_1 proc~comms_scatterv_int_2 comms_scatterv_int_2 interface~comms_scatterv->proc~comms_scatterv_int_2 proc~comms_scatterv_int_3 comms_scatterv_int_3 interface~comms_scatterv->proc~comms_scatterv_int_3 proc~comms_scatterv_real_1 comms_scatterv_real_1 interface~comms_scatterv->proc~comms_scatterv_real_1 proc~comms_scatterv_real_2 comms_scatterv_real_2 interface~comms_scatterv->proc~comms_scatterv_real_2 proc~comms_scatterv_real_3 comms_scatterv_real_3 interface~comms_scatterv->proc~comms_scatterv_real_3 proc~comms_array_split->proc~mpisize proc~get_hh_r->interface~comms_bcast proc~get_hh_r->proc~io_stopwatch_start proc~get_hh_r->proc~io_stopwatch_stop proc~get_hh_r->proc~mpirank proc~get_hh_r->proc~set_error_file proc~get_hh_r->proc~set_error_input proc~fourier_q_to_r fourier_q_to_R proc~get_hh_r->proc~fourier_q_to_r proc~get_win_min get_win_min proc~get_hh_r->proc~get_win_min proc~operator_wigner_setup operator_wigner_setup proc~get_hh_r->proc~operator_wigner_setup proc~set_error_fatal set_error_fatal proc~get_hh_r->proc~set_error_fatal proc~io_date io_date proc~internal_write_header->proc~io_date 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_file->proc~comms_sync_error proc~set_error_file->proc~set_base_error proc~set_error_input->proc~comms_sync_error proc~set_error_input->proc~set_base_error proc~utility_diagonalize->proc~set_error_fatal zhpevx zhpevx proc~utility_diagonalize->zhpevx proc~utility_inv3 utility_inv3 proc~utility_recip_lattice_base->proc~utility_inv3 proc~wham_get_eig_deleig->proc~get_hh_r proc~wham_get_eig_deleig->proc~pw90common_fourier_r_to_k proc~wham_get_eig_deleig->proc~utility_diagonalize proc~wham_get_deleig_a wham_get_deleig_a proc~wham_get_eig_deleig->proc~wham_get_deleig_a proc~comms_bcast_char->proc~comms_sync_error proc~comms_no_sync_bcast_char comms_no_sync_bcast_char proc~comms_bcast_char->proc~comms_no_sync_bcast_char proc~comms_bcast_cmplx->proc~comms_sync_error proc~comms_no_sync_bcast_cmplx comms_no_sync_bcast_cmplx proc~comms_bcast_cmplx->proc~comms_no_sync_bcast_cmplx proc~comms_bcast_int->proc~comms_sync_error proc~comms_no_sync_bcast_int comms_no_sync_bcast_int proc~comms_bcast_int->proc~comms_no_sync_bcast_int proc~comms_bcast_logical->proc~comms_sync_error proc~comms_no_sync_bcast_logical comms_no_sync_bcast_logical proc~comms_bcast_logical->proc~comms_no_sync_bcast_logical proc~comms_bcast_real->proc~comms_sync_error proc~comms_no_sync_bcast_real comms_no_sync_bcast_real proc~comms_bcast_real->proc~comms_no_sync_bcast_real proc~comms_gatherv_cmplx_1->proc~comms_sync_error proc~comms_no_sync_gatherv_cmplx_1 comms_no_sync_gatherv_cmplx_1 proc~comms_gatherv_cmplx_1->proc~comms_no_sync_gatherv_cmplx_1 proc~comms_gatherv_cmplx_2->proc~comms_sync_error proc~comms_no_sync_gatherv_cmplx_2 comms_no_sync_gatherv_cmplx_2 proc~comms_gatherv_cmplx_2->proc~comms_no_sync_gatherv_cmplx_2 proc~comms_gatherv_cmplx_3->proc~comms_sync_error proc~comms_no_sync_gatherv_cmplx_3 comms_no_sync_gatherv_cmplx_3 proc~comms_gatherv_cmplx_3->proc~comms_no_sync_gatherv_cmplx_3 proc~comms_gatherv_cmplx_3_4->proc~comms_sync_error proc~comms_no_sync_gatherv_cmplx_3_4 comms_no_sync_gatherv_cmplx_3_4 proc~comms_gatherv_cmplx_3_4->proc~comms_no_sync_gatherv_cmplx_3_4 proc~comms_gatherv_cmplx_4->proc~comms_sync_error proc~comms_no_sync_gatherv_cmplx_4 comms_no_sync_gatherv_cmplx_4 proc~comms_gatherv_cmplx_4->proc~comms_no_sync_gatherv_cmplx_4 proc~comms_gatherv_logical->proc~comms_sync_error proc~comms_no_sync_gatherv_logical comms_no_sync_gatherv_logical proc~comms_gatherv_logical->proc~comms_no_sync_gatherv_logical proc~comms_gatherv_real_1->proc~comms_sync_error proc~comms_no_sync_gatherv_real_1 comms_no_sync_gatherv_real_1 proc~comms_gatherv_real_1->proc~comms_no_sync_gatherv_real_1 proc~comms_gatherv_real_2->proc~comms_sync_error proc~comms_no_sync_gatherv_real_2 comms_no_sync_gatherv_real_2 proc~comms_gatherv_real_2->proc~comms_no_sync_gatherv_real_2 proc~comms_gatherv_real_2_3->proc~comms_sync_error proc~comms_no_sync_gatherv_real_2_3 comms_no_sync_gatherv_real_2_3 proc~comms_gatherv_real_2_3->proc~comms_no_sync_gatherv_real_2_3 proc~comms_gatherv_real_3->proc~comms_sync_error proc~comms_no_sync_gatherv_real_3 comms_no_sync_gatherv_real_3 proc~comms_gatherv_real_3->proc~comms_no_sync_gatherv_real_3 proc~comms_scatterv_cmplx_3->proc~comms_sync_error proc~comms_no_sync_scatterv_cmplx_3 comms_no_sync_scatterv_cmplx_3 proc~comms_scatterv_cmplx_3->proc~comms_no_sync_scatterv_cmplx_3 proc~comms_scatterv_cmplx_4->proc~comms_sync_error proc~comms_no_sync_scatterv_cmplx_4 comms_no_sync_scatterv_cmplx_4 proc~comms_scatterv_cmplx_4->proc~comms_no_sync_scatterv_cmplx_4 proc~comms_scatterv_int_1->proc~comms_sync_error proc~comms_no_sync_scatterv_int_1 comms_no_sync_scatterv_int_1 proc~comms_scatterv_int_1->proc~comms_no_sync_scatterv_int_1 proc~comms_scatterv_int_2->proc~comms_sync_error proc~comms_no_sync_scatterv_int_2 comms_no_sync_scatterv_int_2 proc~comms_scatterv_int_2->proc~comms_no_sync_scatterv_int_2 proc~comms_scatterv_int_3->proc~comms_sync_error proc~comms_no_sync_scatterv_int_3 comms_no_sync_scatterv_int_3 proc~comms_scatterv_int_3->proc~comms_no_sync_scatterv_int_3 proc~comms_scatterv_real_1->proc~comms_sync_error proc~comms_no_sync_scatterv_real_1 comms_no_sync_scatterv_real_1 proc~comms_scatterv_real_1->proc~comms_no_sync_scatterv_real_1 proc~comms_scatterv_real_2->proc~comms_sync_error proc~comms_no_sync_scatterv_real_2 comms_no_sync_scatterv_real_2 proc~comms_scatterv_real_2->proc~comms_no_sync_scatterv_real_2 proc~comms_scatterv_real_3->proc~comms_sync_error proc~comms_no_sync_scatterv_real_3 comms_no_sync_scatterv_real_3 proc~comms_scatterv_real_3->proc~comms_no_sync_scatterv_real_3 proc~set_error_fatal->proc~comms_sync_error proc~set_error_fatal->proc~set_base_error proc~wham_get_deleig_a->proc~utility_diagonalize proc~utility_rotate utility_rotate proc~wham_get_deleig_a->proc~utility_rotate proc~utility_rotate_diag utility_rotate_diag proc~wham_get_deleig_a->proc~utility_rotate_diag zcopy zcopy proc~comms_no_sync_gatherv_cmplx_2->zcopy proc~comms_no_sync_gatherv_cmplx_3_4->zcopy dcopy dcopy proc~comms_no_sync_gatherv_real_2_3->dcopy proc~utility_matmul_diag utility_matmul_diag proc~utility_rotate_diag->proc~utility_matmul_diag proc~utility_zgemm_new utility_zgemm_new proc~utility_rotate_diag->proc~utility_zgemm_new zgemm zgemm proc~utility_zgemm_new->zgemm

Called by

proc~~geninterp_main~~CalledByGraph proc~geninterp_main geninterp_main program~postw90 postw90 program~postw90->proc~geninterp_main

Source Code

  subroutine geninterp_main(dis_manifold, pw90_geninterp, kpt_latt, pw90_band_deriv_degen, &
                            ws_region, print_output, wannier_data, ws_distance, wigner_seitz, &
                            HH_R, v_matrix, u_matrix, eigval, real_lattice, scissors_shift, &
                            mp_grid, num_bands, num_kpts, num_wann, num_valence_bands, &
                            effective_model, have_disentangled, seedname, stdout, timer, error, &
                            comm)
    !==================================================
    !! This routine prints the band energies (and possibly the band derivatives)
    !!
    !! This routine is parallel, even if ***the scaling is very bad*** since at the moment
    !! everything must be written by the root node (we need that the output is sorted).
    !! But at least if works independently of the number of processors.
    !! I think that a way to write in parallel to the output would help a lot,
    !! so that we don't have to send all eigenvalues to the root node.
    !==================================================

    use w90_constants, only: dp, pi
    use w90_postw90_types, only: pw90_geninterp_mod_type, &
                                 pw90_band_deriv_degen_type, wigner_seitz_type
    use w90_types, only: dis_manifold_type, print_output_type, &
                         wannier_data_type, ws_region_type, ws_distance_type, timer_list_type
    use w90_io, only: io_stopwatch_start, io_stopwatch_stop
    use w90_postw90_common, only: pw90common_fourier_R_to_k
    use w90_utility, only: utility_diagonalize, utility_recip_lattice_base
    use w90_wan_ham, only: wham_get_eig_deleig
    use w90_get_oper, only: get_HH_R
    use w90_comms, only: mpirank, mpisize, comms_bcast, comms_array_split, comms_scatterv, &
                         comms_gatherv, w90_comm_type

    ! arguments
    type(dis_manifold_type), intent(in)          :: dis_manifold
    type(pw90_geninterp_mod_type), intent(in)    :: pw90_geninterp
    type(pw90_band_deriv_degen_type), intent(in) :: pw90_band_deriv_degen
    type(ws_region_type), intent(in)             :: ws_region
    type(print_output_type), intent(in)          :: print_output
    type(wannier_data_type), intent(in)          :: wannier_data
    type(ws_distance_type), intent(inout)        :: ws_distance
    type(wigner_seitz_type), intent(inout)       :: wigner_seitz
    type(timer_list_type), intent(inout)         :: timer
    type(w90_comm_type), intent(in)               :: comm
    type(w90_error_type), allocatable, intent(out) :: error

    complex(kind=dp), allocatable, intent(inout) :: HH_R(:, :, :)
    complex(kind=dp), intent(in) :: v_matrix(:, :, :), u_matrix(:, :, :)

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

    integer, intent(in) :: mp_grid(3)
    integer, intent(in) :: num_bands, num_kpts, num_wann, num_valence_bands, stdout

    character(len=50), intent(in) :: seedname
    logical, intent(in) :: have_disentangled
    logical, intent(in) :: effective_model

    ! local variables
    real(kind=dp) :: recip_lattice(3, 3), volume
    integer :: kpt_unit, outdat_unit, ierr, i, j, enidx
    integer :: nkinterp ! number of kpoints for which we perform the interpolation
    integer, allocatable :: counts(:), displs(:)
    integer :: my_node_id, num_nodes
    integer, allocatable :: kpointidx(:), localkpointidx(:)
    real(kind=dp), allocatable :: kpoints(:, :), localkpoints(:, :)
    real(kind=dp)  :: kpt(3), frac(3)
    real(kind=dp), allocatable :: localdeleig(:, :, :)
    real(kind=dp), allocatable :: globaldeleig(:, :, :)
    real(kind=dp), allocatable :: localeig(:, :)
    real(kind=dp), allocatable :: globaleig(:, :)
    complex(kind=dp), allocatable :: HH(:, :)
    complex(kind=dp), allocatable :: UU(:, :)
    complex(kind=dp), allocatable :: delHH(:, :, :)
    character(len=500) :: commentline
    character(len=50) :: cdum
    character(len=200) :: outdat_filename
    logical :: absoluteCoords
    logical :: on_root = .false.

    my_node_id = mpirank(comm)
    num_nodes = mpisize(comm)
    if (my_node_id == 0) on_root = .true.
    allocate (counts(0:num_nodes - 1))
    allocate (displs(0:num_nodes - 1))

    if (print_output%iprint > 0 .and. (print_output%timing_level > 0)) &
      call io_stopwatch_start('geninterp_main', timer)

    if (on_root) then
      write (stdout, *)
      write (stdout, '(1x,a)') '*---------------------------------------------------------------------------*'
      write (stdout, '(1x,a)') '|                      Generic Band Interpolation routines                  |'
      write (stdout, '(1x,a)') '*---------------------------------------------------------------------------*'

      open (newunit=kpt_unit, file=trim(seedname)//'_geninterp.kpt', form='formatted', status='old', &
            err=105)

      ! First line: comment (e.g. creation date, author, ...)
      read (kpt_unit, '(A500)', err=106, end=106) commentline
      read (kpt_unit, *, err=106, end=106) cdum

      if (index(cdum, 'crystal') > 0) then
        absoluteCoords = .false.
      elseif (index(cdum, 'frac') > 0) then
        absoluteCoords = .false.
      elseif (index(cdum, 'cart') > 0) then
        absoluteCoords = .true.
      elseif (index(cdum, 'abs') > 0) then
        absoluteCoords = .true.
      else
        call set_error_input(error, 'Error on second line of file '//trim(seedname) &
                             //'_geninterp.kpt: unable to recognize keyword', comm)
        return
      end if

      ! Third line: number of following kpoints
      read (kpt_unit, *, err=106, end=106) nkinterp
    end if

    call comms_bcast(nkinterp, 1, error, comm)
    if (allocated(error)) return

    allocate (HH(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating HH in calcTDF', comm)
      return
    end if
    allocate (UU(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating UU in calcTDF', comm)
      return
    end if
    if (pw90_geninterp%alsofirstder) then
      allocate (delHH(num_wann, num_wann, 3), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating delHH in calcTDF', comm)
        return
      end if
    end if

    ! I call once the routine to calculate the Hamiltonian in real-space <0n|H|Rm>
    call get_HH_R(dis_manifold, kpt_latt, print_output, wigner_seitz, HH_R, u_matrix, v_matrix, &
                  eigval, real_lattice, scissors_shift, num_bands, num_kpts, num_wann, &
                  num_valence_bands, effective_model, have_disentangled, seedname, ws_distance, ws_region, &
                  stdout, timer, error, comm)
    if (allocated(error)) return

    if (on_root) then
      allocate (kpointidx(nkinterp), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating kpointidx in geinterp_main.', comm)
        return
      end if
      allocate (kpoints(3, nkinterp), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating kpoints in geinterp_main.', comm)
        return
      end if
      if (pw90_geninterp%single_file) then
        allocate (globaleig(num_wann, nkinterp), stat=ierr)
        if (ierr /= 0) then
          call set_error_alloc(error, 'Error allocating globaleig in geinterp_main.', comm)
          return
        end if
        allocate (globaldeleig(num_wann, 3, nkinterp), stat=ierr)
        if (ierr /= 0) then
          call set_error_alloc(error, 'Error allocating globaldeleig in geinterp_main.', comm)
          return
        end if
      end if
    else
      ! On the other nodes, I still allocate them with size 1 to avoid
      ! that some compilers still try to access the memory
      allocate (kpointidx(1), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating kpointidx in geinterp_main.', comm)
        return
      end if
      allocate (kpoints(1, 1), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating kpoints in geinterp_main.', comm)
        return
      end if
      if (pw90_geninterp%single_file) then
        allocate (globaleig(num_wann, 1), stat=ierr)
        if (ierr /= 0) then
          call set_error_alloc(error, 'Error allocating globaleig in geinterp_main.', comm)
          return
        end if
        allocate (globaldeleig(num_wann, 3, 1), stat=ierr)
        if (ierr /= 0) then
          call set_error_alloc(error, 'Error allocating globaldeleig in geinterp_main.', comm)
          return
        end if
      end if
    end if

    ! I precalculate how to split on different nodes
    call comms_array_split(nkinterp, counts, displs, comm)

    allocate (localkpoints(3, max(1, counts(my_node_id))), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error allocating localkpoints in geinterp_main.', comm)
      return
    end if
    allocate (localeig(num_wann, max(1, counts(my_node_id))), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error allocating localeig in geinterp_main.', comm)
      return
    end if
    allocate (localdeleig(num_wann, 3, max(1, counts(my_node_id))), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error allocating localdeleig in geinterp_main.', comm)
      return
    end if

    ! On root, I read numpoints_thischunk points
    if (on_root) then
      ! Lines with integer identifier and three coordinates
      ! (in crystallographic coordinates relative to the reciprocal lattice vectors)
      do i = 1, nkinterp
        read (kpt_unit, *, err=106, end=106) kpointidx(i), kpt
        ! Internally, I need the relative (fractional) coordinates in units of the reciprocal-lattice vectors
        if (absoluteCoords .eqv. .false.) then
          kpoints(:, i) = kpt
        else
          kpoints(:, i) = 0._dp
          ! I use the real_lattice (transposed) and a factor of 2pi instead of inverting again recip_lattice
          do j = 1, 3
            kpoints(j, i) = real_lattice(j, 1)*kpt(1) + real_lattice(j, 2)*kpt(2) + &
                            real_lattice(j, 3)*kpt(3)
          end do
          kpoints(:, i) = kpoints(:, i)/(2._dp*pi)
        end if
      end do
      close (kpt_unit)
    end if

    ! Now, I distribute the kpoints; 3* because I send kx, ky, kz
    call comms_scatterv(localkpoints, 3*counts(my_node_id), kpoints, 3*counts, 3*displs, error, comm)
    if (allocated(error)) return
    if (.not. pw90_geninterp%single_file) then
      ! Allocate at least one entry, even if we don't use it
      allocate (localkpointidx(max(1, counts(my_node_id))), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating localkpointidx in geinterp_main.', comm)
        return
      end if
      call comms_scatterv(localkpointidx, counts(my_node_id), kpointidx, counts, displs, error, comm)
      if (allocated(error)) return
    end if

    ! I open the output file(s)
    if (pw90_geninterp%single_file) then
      if (on_root) then
        outdat_filename = trim(seedname)//'_geninterp.dat'
        open (newunit=outdat_unit, file=trim(outdat_filename), form='formatted', err=107)

        call internal_write_header(outdat_unit, commentline, pw90_geninterp)
      end if
    else
      if (num_nodes > 99999) then
        write (outdat_filename, '(a,a,I0,a)') trim(seedname), '_geninterp_', my_node_id, '.dat'
      else
        write (outdat_filename, '(a,a,I5.5,a)') trim(seedname), '_geninterp_', my_node_id, '.dat'
      end if
      open (newunit=outdat_unit, file=trim(outdat_filename), form='formatted', err=107)

      call comms_bcast(commentline, len(commentline), error, comm)
      if (allocated(error)) return

      call internal_write_header(outdat_unit, commentline, pw90_geninterp)
    end if

    ! And now, each node calculates its own k points
    do i = 1, counts(my_node_id)
      kpt = localkpoints(:, i)
      ! Here I get the band energies and the velocities (if required)
      if (pw90_geninterp%alsofirstder) then
        call wham_get_eig_deleig(dis_manifold, kpt_latt, pw90_band_deriv_degen, ws_region, &
                                 print_output, wannier_data, ws_distance, wigner_seitz, delHH, HH, &
                                 HH_R, u_matrix, UU, v_matrix, localdeleig(:, :, i), &
                                 localeig(:, i), eigval, kpt, real_lattice, &
                                 scissors_shift, mp_grid, num_bands, num_kpts, num_wann, &
                                 num_valence_bands, effective_model, have_disentangled, &
                                 seedname, stdout, timer, error, comm)
        if (allocated(error)) return

      else
        call pw90common_fourier_R_to_k(ws_region, wannier_data, ws_distance, wigner_seitz, HH, &
                                       HH_R, kpt, real_lattice, mp_grid, 0, num_wann, error, comm)
        if (allocated(error)) return

        call utility_diagonalize(HH, num_wann, localeig(:, i), UU, error, comm)
        if (allocated(error)) return

      end if
    end do

    call utility_recip_lattice_base(real_lattice, recip_lattice, volume)
    if (pw90_geninterp%single_file) then
      ! Now, I get the results from the different nodes
      call comms_gatherv(localeig, num_wann*counts(my_node_id), globaleig, &
                         num_wann*counts, num_wann*displs, error, comm)
      if (allocated(error)) return

      if (pw90_geninterp%alsofirstder) then
        call comms_gatherv(localdeleig, 3*num_wann*counts(my_node_id), globaldeleig, &
                           3*num_wann*counts, 3*num_wann*displs, error, comm)
        if (allocated(error)) return
      end if

      ! Now the printing, only on root node
      if (on_root) then
        do i = 1, nkinterp
          kpt = kpoints(:, i)
          ! First calculate the absolute coordinates for printing
          frac = 0._dp
          do j = 1, 3
            frac(j) = recip_lattice(1, j)*kpt(1) + recip_lattice(2, j)*kpt(2) &
                      + recip_lattice(3, j)*kpt(3)
          end do

          ! I print each line
          if (pw90_geninterp%alsofirstder) then
            do enidx = 1, num_wann
              write (outdat_unit, '(I10,7G18.10)') kpointidx(i), frac, &
                globaleig(enidx, i), globaldeleig(enidx, :, i)
            end do
          else
            do enidx = 1, num_wann
              write (outdat_unit, '(I10,4G18.10)') kpointidx(i), frac, globaleig(enidx, i)
            end do
          end if
        end do
        close (outdat_unit)
      end if
    else
      ! Each node simply writes to its own file
      do i = 1, counts(my_node_id)
        kpt = localkpoints(:, i)
        ! First calculate the absolute coordinates for printing
        frac = 0._dp
        do j = 1, 3
          frac(j) = recip_lattice(1, j)*kpt(1) + recip_lattice(2, j)*kpt(2) &
                    + recip_lattice(3, j)*kpt(3)
        end do

        ! I print each line
        if (pw90_geninterp%alsofirstder) then
          do enidx = 1, num_wann
            write (outdat_unit, '(I10,7G18.10)') localkpointidx(i), frac, &
              localeig(enidx, i), localdeleig(enidx, :, i)
          end do
        else
          do enidx = 1, num_wann
            write (outdat_unit, '(I10,4G18.10)') localkpointidx(i), frac, localeig(enidx, i)
          end do
        end if
      end do
      close (outdat_unit)
    end if

    ! All k points processed: Final processing/deallocations

    if (on_root) then
      write (stdout, '(1x,a)') '|                                All done.                                  |'
      write (stdout, '(1x,a)') '*---------------------------------------------------------------------------*'

    end if

    if (allocated(kpointidx)) deallocate (kpointidx)
    if (allocated(kpoints)) deallocate (kpoints)
    if (allocated(localkpoints)) deallocate (localkpoints)
    if (allocated(HH)) deallocate (HH)
    if (allocated(UU)) deallocate (UU)
    if (allocated(delHH)) deallocate (delHH)
    if (allocated(localeig)) deallocate (localeig)
    if (allocated(localdeleig)) deallocate (localdeleig)
    if (allocated(globaleig)) deallocate (globaleig)
    if (allocated(globaldeleig)) deallocate (globaldeleig)

    if (on_root .and. (print_output%timing_level > 0)) call io_stopwatch_stop('geninterp_main', timer)

    return

105 call set_error_file(error, 'Error: Problem opening k-point file '//trim(seedname)//'_geninterp.kpt', comm)
    return
106 call set_error_file(error, 'Error: Problem reading k-point file '//trim(seedname)//'_geninterp.kpt', comm)
    return
107 call set_error_file(error, 'Error: Problem opening output file '//trim(outdat_filename), comm)
    return

  end subroutine geninterp_main