!-*- mode: F90 -*-!
!------------------------------------------------------------!
! Copyright (C) 2026 Wannier Developer Group                 !
!                                                            !
! This library is free software; you can redistribute it     !
! and/or modify it under the terms of the GNU Lesser General !
! Public License as published by the Free Software           !
! Foundation; either version 2.1 of the License, or (at your !
! option) any later version.                                 !
!                                                            !
! This library is distributed in the hope that it will be    !
! useful,but WITHOUT ANY WARRANTY; without even the implied  !
! warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR    !
! PURPOSE.  See the GNU Lesser General Public License for    !
! more details.                                              !
!                                                            !
! You should have received a copy of the GNU Lesser General  !
! Public License along with this library; if not, see        !
! <https://www.gnu.org/licenses/>.                           !
!                                                            !
! The webpage of the Wannier90 code is                       !
! <https://www.wannier.org>.                                 !
!                                                            !
! The Wannier90 code is hosted on GitHub                     !
! <https://github.com/wannier-developers/wannier90>          !
!------------------------------------------------------------!
!                                                            !
!  w90_kpath: evaluate properties along a path in k-space    !
!                                                            !
!------------------------------------------------------------!

module w90_kpath

  !! Calculates quantities along a specified k-path:
  !!
  !!  - Energy bands (eventually colored by the spin)
  !!  - Energy bands (colored by Berry curvature-like
  !!                  term of spin Hall conductivity)
  !!
  !!  - (Berry curvature)x(-1) summed over occupied bands
  !!  - Berry curvature-like term of spin Hall conductivity
  !!    summed over occupied bands
  !!
  !!  - Integrand of orbital magnetization Morb=LCtil+ICtil

  use w90_constants, only: dp
  use w90_error, only: w90_error_type, set_error_alloc, set_error_dealloc, set_error_fatal, &
                       set_error_input, set_error_fatal, set_error_file

  implicit none

  private

  public :: k_path

contains

  !================================================!
  !                   PUBLIC PROCEDURES
  !================================================!

  subroutine k_path(pw90_berry, dis_manifold, fermi_energy_list, kmesh_info, pw90_kpath, kpt_latt, &
                    pw90_oper_read, pw90_band_deriv_degen, pw90_spin, ws_region, kpoint_path, &
                    pw90_spin_hall, print_output, wannier_data, ws_distance, wigner_seitz, AA_R, &
                    BB_R, CC_R, HH_R, SH_R, SHR_R, SR_R, SS_R, SAA_R, SBB_R, v_matrix, u_matrix, &
                    bohr, eigval, real_lattice, scissors_shift, mp_grid, fermi_n, num_wann, &
                    num_bands, num_kpts, num_valence_bands, effective_model, have_disentangled, &
                    seedname, stdout, timer, error, comm)
    !================================================!
    !
    !! Main routine
    !
    !================================================!

    use w90_comms, only: w90_comm_type, mpirank, mpisize, comms_array_split, comms_scatterv, &
                         comms_gatherv, comms_bcast
    use w90_constants, only: dp, eps8
    use w90_get_oper, only: get_HH_R, get_AA_R_effective, get_AA_R, get_BB_R, get_CC_R, get_SS_R, get_SHC_R
    use w90_io, only: io_time
    use w90_postw90_common, only: pw90common_fourier_R_to_k
    use w90_types, only: kpoint_path_type, print_output_type, wannier_data_type, &
                         dis_manifold_type, kmesh_info_type, ws_region_type, ws_distance_type, timer_list_type
    use w90_postw90_types, only: pw90_berry_mod_type, pw90_spin_hall_type, pw90_kpath_mod_type, &
                                 pw90_spin_mod_type, pw90_band_deriv_degen_type, pw90_oper_read_type, wigner_seitz_type
    use w90_berry, only: berry_get_imf_klist, berry_get_imfgh_klist, berry_get_shc_klist
    use w90_spin, only: spin_get_nk
    use w90_utility, only: utility_diagonalize, utility_recip_lattice_base

    implicit none

    ! arguments
    type(pw90_berry_mod_type), intent(in) :: pw90_berry
    type(dis_manifold_type), intent(in) :: dis_manifold
    type(kmesh_info_type), intent(in) :: kmesh_info
    type(pw90_kpath_mod_type), intent(in) :: pw90_kpath
    type(pw90_band_deriv_degen_type), intent(in) :: pw90_band_deriv_degen
    type(pw90_oper_read_type), intent(in) :: pw90_oper_read
    type(pw90_spin_mod_type), intent(in) :: pw90_spin
    type(print_output_type), intent(in) :: print_output
    type(ws_region_type), intent(in) :: ws_region
    type(kpoint_path_type), intent(in) :: kpoint_path
    type(pw90_spin_hall_type), intent(in) :: pw90_spin_hall
    type(w90_comm_type), intent(in) :: comm
    type(wannier_data_type), intent(in) :: wannier_data
    type(wigner_seitz_type), intent(inout) :: wigner_seitz
    type(ws_distance_type), intent(inout) :: ws_distance
    type(timer_list_type), intent(inout) :: timer
    type(w90_error_type), allocatable, intent(out) :: error

    complex(kind=dp), allocatable, intent(inout) :: AA_R(:, :, :, :) ! <0n|r|Rm>
    complex(kind=dp), allocatable, intent(inout) :: BB_R(:, :, :, :) ! <0|H(r-R)|R>
    complex(kind=dp), allocatable, intent(inout) :: CC_R(:, :, :, :, :) ! <0|r_alpha.H(r-R)_beta|R>
    complex(kind=dp), allocatable, intent(inout) :: HH_R(:, :, :) !  <0n|r|Rm>
    complex(kind=dp), allocatable, intent(inout) :: SH_R(:, :, :, :) ! <0n|sigma_x,y,z.H|Rm>
    complex(kind=dp), allocatable, intent(inout) :: SHR_R(:, :, :, :, :) ! <0n|sigma_x,y,z.H.(r-R)_alpha|Rm>
    complex(kind=dp), allocatable, intent(inout) :: SR_R(:, :, :, :, :) ! <0n|sigma_x,y,z.(r-R)_alpha|Rm>
    complex(kind=dp), allocatable, intent(inout) :: SS_R(:, :, :, :) ! <0n|sigma_x,y,z|Rm>
    complex(kind=dp), allocatable, intent(inout) :: SAA_R(:, :, :, :, :) !<0n|sigma_x,y,z.(r-R)_alpha|Rm>
    complex(kind=dp), allocatable, intent(inout) :: SBB_R(:, :, :, :, :) !<0n|sigma_x,y,z.H.(r-R)_alpha|Rm>
    complex(kind=dp), intent(in) :: v_matrix(:, :, :), u_matrix(:, :, :)

    real(kind=dp), intent(in) :: bohr
    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), allocatable, intent(in) :: fermi_energy_list(:)
    real(kind=dp), intent(in) :: kpt_latt(:, :)

    integer, intent(in) :: mp_grid(3)
    integer, intent(in) :: num_wann, num_bands, num_kpts, num_valence_bands, fermi_n
    integer, intent(in) :: 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           :: i, j, n, num_paths, num_spts, loop_kpt, &
                         total_pts, loop_i, dataunit, gnuunit, pyunit, &
                         my_num_pts
    real(kind=dp)     :: ymin, ymax, kpt(3), spn_k(num_wann), &
                         imf_k_list(3, 3, fermi_n), img_k_list(3, 3, fermi_n), &
                         imh_k_list(3, 3, fermi_n), Morb_k(3, 3), &
                         range, zmin, zmax
    real(kind=dp)     :: shc_k_band(num_wann), shc_k_fermi(fermi_n)
    real(kind=dp), allocatable :: kpath_len(:)
    logical           :: plot_bands, plot_curv, plot_morb, plot_shc
    character(len=120) :: file_name

    complex(kind=dp), allocatable :: HH(:, :)
    complex(kind=dp), allocatable :: UU(:, :)
    real(kind=dp), allocatable    :: xval(:), eig(:, :), my_eig(:, :), &
                                     curv(:, :), my_curv(:, :), &
                                     morb(:, :), my_morb(:, :), &
                                     color(:, :), my_color(:, :), &
                                     plot_kpoint(:, :), my_plot_kpoint(:, :), &
                                     shc(:), my_shc(:)
    character(len=3), allocatable  :: glabel(:)

    integer, allocatable :: counts(:), displs(:)
    logical :: on_root = .false.
    integer :: my_node_id, num_nodes

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

    ! Everything is done on the root node (not worthwhile parallelizing)
    ! However, we still have to read and distribute the data if we
    ! are in parallel. So calls to get_oper are done on all nodes at the moment

    plot_bands = index(pw90_kpath%task, 'bands') > 0
    plot_curv = index(pw90_kpath%task, 'curv') > 0
    plot_morb = index(pw90_kpath%task, 'morb') > 0
    plot_shc = index(pw90_kpath%task, 'shc') > 0

    if (on_root) then
      if (plot_shc .or. (plot_bands .and. pw90_kpath%bands_colour == 'shc')) then
        ! not allowed to use adpt smr, since adpt smr needs berry_kmesh,
        ! see line 1837 of berry.F90
        if (pw90_berry%kubo_smearing%use_adaptive) then
          call set_error_input(error, 'Error: Must use fixed smearing when plotting spin Hall conductivity', comm)
          return
        end if
      end if
      if (plot_shc) then
        if (fermi_n == 0) then
          call set_error_input(error, 'Error: must specify Fermi energy', comm)
          return
        else if (fermi_n /= 1) then
          call set_error_input(error, 'Error: kpath plot only accept one Fermi energy, ' &
                               //'use fermi_energy instead of fermi_energy_min', comm)
          return
        end if
      end if
    end if

    call k_path_print_info(plot_bands, plot_curv, plot_morb, plot_shc, fermi_energy_list, &
                           pw90_kpath, pw90_berry%curv_unit, stdout, error, comm)
    if (allocated(error)) return

    ! Set up the needed Wannier matrix elements

    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 (plot_curv .or. plot_morb) then
      if (effective_model) then
        call get_AA_R_effective(print_output, AA_R, HH_R, wigner_seitz%nrpts, num_wann, seedname, &
                                stdout, timer, error, comm)
      else
        call get_AA_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, wannier_data, AA_R, &
                      v_matrix, eigval, wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, &
                      num_wann, have_disentangled, seedname, stdout, timer, error, comm)
      end if
      if (allocated(error)) return

    end if
    if (plot_morb) then

      call get_BB_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, HH_R, BB_R, v_matrix, &
                    eigval, scissors_shift, wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, &
                    num_wann, have_disentangled, seedname, stdout, timer, error, comm)
      if (allocated(error)) return

      call get_CC_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, pw90_oper_read, &
                    HH_R, BB_R, CC_R, v_matrix, eigval, scissors_shift, wigner_seitz, ws_distance, &
                    ws_region, num_bands, num_kpts, num_wann, have_disentangled, seedname, stdout, &
                    timer, error, comm)
      if (allocated(error)) return

    end if

    if (plot_shc .or. (plot_bands .and. pw90_kpath%bands_colour == 'shc')) then

      if (effective_model) then
        call get_AA_R_effective(print_output, AA_R, HH_R, wigner_seitz%nrpts, num_wann, seedname, &
                                stdout, timer, error, comm)
      else
        call get_AA_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, wannier_data, AA_R, &
                      v_matrix, eigval, wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, &
                      num_wann, have_disentangled, seedname, stdout, timer, error, comm)
      end if
      if (allocated(error)) return

      call get_SS_R(dis_manifold, kpt_latt, print_output, pw90_oper_read, SS_R, v_matrix, eigval, &
                    wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, num_wann, &
                    have_disentangled, seedname, stdout, timer, error, comm)
      if (allocated(error)) return

      call get_SHC_R(dis_manifold, kmesh_info, kpt_latt, print_output, pw90_oper_read, &
                     pw90_spin_hall, SH_R, SHR_R, SR_R, v_matrix, eigval, scissors_shift, &
                     wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, num_wann, &
                     num_valence_bands, have_disentangled, seedname, stdout, timer, error, comm)
      if (allocated(error)) return

    end if

    if (plot_bands .and. pw90_kpath%bands_colour == 'spin') then
      call get_SS_R(dis_manifold, kpt_latt, print_output, pw90_oper_read, SS_R, v_matrix, eigval, &
                    wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, num_wann, &
                    have_disentangled, seedname, stdout, timer, error, comm)
      if (allocated(error)) return

    end if

    num_paths = 0
    if (on_root) then
      call utility_recip_lattice_base(real_lattice, recip_lattice, volume)
      ! Determine the number of k-points (total_pts) as well as
      ! their reciprocal-lattice coordinates long the path (plot_kpoint)
      ! and their associated horizontal coordinate for the plot (xval)
      call k_path_get_points(num_paths, kpath_len, total_pts, xval, plot_kpoint, kpoint_path, &
                             recip_lattice, pw90_kpath)
      ! (paths)
      num_spts = num_paths + 1 ! number of path endpoints (special pts)
    else
      ! Dummy allocation for making scatterv work
      allocate (plot_kpoint(1, 1))
    end if

    ! Broadcast number of k-points on the path
    call comms_bcast(total_pts, 1, error, comm)
    if (allocated(error)) return

    ! Partition set of k-points into junks
!   call comms_array_split(total_pts, counts, displs)
    call comms_array_split(total_pts, counts, displs, comm)
    !kpt_lo = displs(my_node_id)+1
    !kpt_hi = displs(my_node_id)+counts(my_node_id)
    my_num_pts = counts(my_node_id)

    ! Distribute coordinates
    allocate (my_plot_kpoint(3, my_num_pts))
    call comms_scatterv(my_plot_kpoint, 3*my_num_pts, &
                        plot_kpoint, 3*counts, 3*displs, error, comm)
    if (allocated(error)) return

    ! Value of the vertical coordinate in the actual plots: energy bands
    !
    if (plot_bands) then
      allocate (HH(num_wann, num_wann))
      allocate (UU(num_wann, num_wann))
      allocate (my_eig(num_wann, my_num_pts))
      if (pw90_kpath%bands_colour /= 'none') allocate (my_color(num_wann, my_num_pts))
    end if

    ! Value of the vertical coordinate in the actual plots
    !
    if (plot_curv) allocate (my_curv(my_num_pts, 3))
    if (plot_morb) allocate (my_morb(my_num_pts, 3))
    if (plot_shc) allocate (my_shc(my_num_pts))

    ! Loop over local junk of k-points on the path and evaluate the requested quantities
    !
    do loop_kpt = 1, my_num_pts
      kpt(:) = my_plot_kpoint(:, loop_kpt)

      if (plot_bands) then
        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, my_eig(:, loop_kpt), UU, error, comm)
        if (allocated(error)) return

        !
        ! Color-code energy bands with the spin projection along the
        ! chosen spin quantization axis
        !
        if (pw90_kpath%bands_colour == 'spin') then
          call spin_get_nk(ws_region, pw90_spin, wannier_data, ws_distance, wigner_seitz, HH_R, &
                           SS_R, kpt, real_lattice, spn_k, mp_grid, num_wann, error, comm)
          if (allocated(error)) return

          my_color(:, loop_kpt) = spn_k(:)
          !
          ! The following is needed to prevent bands from disappearing
          ! when the magnitude of the Wannier interpolated spn_k (very
          ! slightly) exceeds 1.0 (e.g. in bcc Fe along N--G--H)
          !
          do n = 1, num_wann
            if (my_color(n, loop_kpt) > 1.0_dp - eps8) then
              my_color(n, loop_kpt) = 1.0_dp - eps8
            else if (my_color(n, loop_kpt) < -1.0_dp + eps8) then
              my_color(n, loop_kpt) = -1.0_dp + eps8
            end if
          end do
        else if (pw90_kpath%bands_colour == 'shc') then
          call berry_get_shc_klist(pw90_berry, dis_manifold, fermi_energy_list, kpt_latt, &
                                   pw90_band_deriv_degen, ws_region, pw90_spin_hall, print_output, &
                                   wannier_data, ws_distance, wigner_seitz, AA_R, HH_R, SH_R, &
                                   SHR_R, SR_R, SS_R, SAA_R, SBB_R, u_matrix, v_matrix, eigval, &
                                   kpt, real_lattice, scissors_shift, mp_grid, fermi_n, num_bands, &
                                   num_kpts, num_wann, num_valence_bands, effective_model, &
                                   have_disentangled, seedname, stdout, timer, error, comm, &
                                   shc_k_band=shc_k_band)
          if (allocated(error)) return

          my_color(:, loop_kpt) = shc_k_band
        end if
      end if

      if (plot_morb) then
        call berry_get_imfgh_klist(dis_manifold, fermi_energy_list, kpt_latt, ws_region, &
                                   print_output, wannier_data, ws_distance, wigner_seitz, AA_R, &
                                   BB_R, CC_R, HH_R, u_matrix, v_matrix, eigval, kpt, &
                                   real_lattice, scissors_shift, mp_grid, fermi_n, num_bands, &
                                   num_kpts, num_wann, num_valence_bands, effective_model, &
                                   have_disentangled, seedname, stdout, timer, error, comm, &
                                   imf_k_list, img_k_list, imh_k_list)
        if (allocated(error)) return

        Morb_k = img_k_list(:, :, 1) + imh_k_list(:, :, 1) &
                 - 2.0_dp*fermi_energy_list(1)*imf_k_list(:, :, 1)
        Morb_k = -Morb_k/2.0_dp ! differs by -1/2 from Eq.97 LVTS12
        my_morb(loop_kpt, 1) = sum(Morb_k(:, 1))
        my_morb(loop_kpt, 2) = sum(Morb_k(:, 2))
        my_morb(loop_kpt, 3) = sum(Morb_k(:, 3))
      end if

      if (plot_curv) then
        if (.not. plot_morb) then
          call berry_get_imf_klist(dis_manifold, fermi_energy_list, kpt_latt, ws_region, &
                                   print_output, wannier_data, ws_distance, wigner_seitz, AA_R, &
                                   BB_R, CC_R, HH_R, u_matrix, v_matrix, eigval, kpt, &
                                   real_lattice, imf_k_list, 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

        end if
        my_curv(loop_kpt, 1) = sum(imf_k_list(:, 1, 1))
        my_curv(loop_kpt, 2) = sum(imf_k_list(:, 2, 1))
        my_curv(loop_kpt, 3) = sum(imf_k_list(:, 3, 1))
      end if

      if (plot_shc) then
        call berry_get_shc_klist(pw90_berry, dis_manifold, fermi_energy_list, kpt_latt, &
                                 pw90_band_deriv_degen, ws_region, pw90_spin_hall, print_output, &
                                 wannier_data, ws_distance, wigner_seitz, AA_R, HH_R, SH_R, SHR_R, &
                                 SR_R, SS_R, SAA_R, SBB_R, u_matrix, v_matrix, eigval, kpt, &
                                 real_lattice, scissors_shift, mp_grid, fermi_n, num_bands, &
                                 num_kpts, num_wann, num_valence_bands, effective_model, &
                                 have_disentangled, seedname, stdout, timer, error, comm, &
                                 shc_k_fermi=shc_k_fermi)
        if (allocated(error)) return
        my_shc(loop_kpt) = shc_k_fermi(1)
      end if
    end do !loop_kpt

    ! Send results to root process
    if (plot_bands) then
      allocate (eig(num_wann, total_pts))
      call comms_gatherv(my_eig, num_wann*my_num_pts, &
                         eig, num_wann*counts, num_wann*displs, error, comm)
      if (allocated(error)) return
      if (pw90_kpath%bands_colour /= 'none') then
        allocate (color(num_wann, total_pts))
        call comms_gatherv(my_color, num_wann*my_num_pts, &
                           color, num_wann*counts, num_wann*displs, error, comm)
        if (allocated(error)) return
      end if
    end if

    if (plot_curv) then
      allocate (curv(total_pts, 3))
      do i = 1, 3
        call comms_gatherv(my_curv(:, i), my_num_pts, &
                           curv(:, i), counts, displs, error, comm)
        if (allocated(error)) return
      end do
    end if

    if (plot_morb) then
      allocate (morb(total_pts, 3))
      do i = 1, 3
        call comms_gatherv(my_morb(:, i), my_num_pts, &
                           morb(:, i), counts, displs, error, comm)
        if (allocated(error)) return
      end do
    end if

    if (plot_shc) then
      allocate (shc(total_pts))
      call comms_gatherv(my_shc, my_num_pts, shc, counts, displs, error, comm)
      if (allocated(error)) return
    end if

    if (on_root) then
      num_spts = num_paths + 1
      allocate (glabel(num_spts))

      if (plot_bands) then
        !
        ! Write out the kpoints in the path in a format that can be inserted
        ! directly in the pwscf input file (the '1.0_dp' in the second column is
        ! a k-point weight, expected by pwscf)
        !
        open (newunit=dataunit, file=trim(seedname)//'-path.kpt', form='formatted')
        write (dataunit, *) total_pts
        do loop_kpt = 1, total_pts
          write (dataunit, '(3f12.6,3x,f4.1)') &
            (plot_kpoint(loop_i, loop_kpt), loop_i=1, 3), 1.0_dp
        end do
        close (dataunit)
      end if
      if (plot_curv .and. pw90_berry%curv_unit == 'bohr2') curv = curv/bohr**2

      if (plot_bands .and. pw90_kpath%bands_colour == 'shc') then
        if (pw90_berry%curv_unit == 'bohr2') color = color/bohr**2
      end if
      if (plot_shc) then
        if (pw90_berry%curv_unit == 'bohr2') shc = shc/bohr**2
      end if

      ! Axis labels
      !
      glabel(1) = ' '//kpoint_path%labels(1)//' '
      do i = 2, num_paths
        if (kpoint_path%labels(2*(i - 1)) /= kpoint_path%labels(2*(i - 1) + 1)) then
          glabel(i) = kpoint_path%labels(2*(i - 1))//'/'//kpoint_path%labels(2*(i - 1) + 1)
        else
          glabel(i) = ' '//kpoint_path%labels(2*(i - 1))//' '
        end if
      end do
      glabel(num_spts) = ' '//kpoint_path%labels(num_paths*2)//' '

      ! Now write the plotting files

      write (stdout, '(/,1x,a)') 'Output files:'

      if (plot_bands) then
        file_name = trim(seedname)//'-path.kpt'
        write (stdout, '(/,3x,a)') file_name
        !
        ! Data file
        !
        file_name = trim(seedname)//'-bands.dat'
        write (stdout, '(/,3x,a)') file_name
        open (newunit=dataunit, file=file_name, form='formatted')
        do i = 1, num_wann
          do loop_kpt = 1, total_pts
            if (pw90_kpath%bands_colour == 'none') then
              write (dataunit, '(2E16.8)') xval(loop_kpt), eig(i, loop_kpt)
            else
              write (dataunit, '(3E16.8)') xval(loop_kpt), &
                eig(i, loop_kpt), color(i, loop_kpt)
            end if
          end do
          write (dataunit, *) ' '
        end do
        close (dataunit)
      end if

      if (plot_bands .and. .not. plot_curv .and. .not. plot_morb &
          .and. .not. plot_shc) then
        !
        ! Gnuplot script
        !
        ymin = minval(eig) - 1.0_dp
        ymax = maxval(eig) + 1.0_dp
        file_name = trim(seedname)//'-bands.gnu'
        write (stdout, '(/,3x,a)') file_name
        open (newunit=gnuunit, file=file_name, form='formatted')
        do i = 1, num_paths - 1
          write (gnuunit, 705) sum(kpath_len(1:i)), ymin, &
            sum(kpath_len(1:i)), ymax
        end do
        if (pw90_kpath%bands_colour == 'none') then
          write (gnuunit, 701) xval(total_pts), ymin, ymax
          write (gnuunit, 702, advance="no") glabel(1), 0.0_dp, &
            (glabel(i + 1), sum(kpath_len(1:i)), i=1, num_paths - 1)
          write (gnuunit, 703) glabel(1 + num_paths), sum(kpath_len(:))
          write (gnuunit, *) 'plot ', '"'//trim(seedname)//'-bands.dat', '"'
        else if (pw90_kpath%bands_colour == 'spin') then
          !
          ! Only works with gnuplot v4.2 and higher
          !
          write (gnuunit, 706) xval(total_pts), ymin, ymax
          write (gnuunit, 702, advance="no") glabel(1), 0.0_dp, &
            (glabel(i + 1), sum(kpath_len(1:i)), i=1, num_paths - 1)
          write (gnuunit, 703) glabel(1 + num_paths), sum(kpath_len(:))
          write (gnuunit, *) &
            'set palette defined (-1 "blue", 0 "green", 1 "red")'
          write (gnuunit, *) 'set pm3d map'
          write (gnuunit, *) 'set zrange [-1:1]'
          write (gnuunit, *) 'splot ', '"'//trim(seedname)//'-bands.dat', &
            '" with dots palette'
        else if (pw90_kpath%bands_colour == 'shc') then
          !
          ! Only works with gnuplot v4.2 and higher
          !
          write (gnuunit, 706) xval(total_pts), ymin, ymax
          write (gnuunit, 702, advance="no") glabel(1), 0.0_dp, &
            (glabel(i + 1), sum(kpath_len(1:i)), i=1, num_paths - 1)
          write (gnuunit, 703) glabel(1 + num_paths), sum(kpath_len(:))
          write (gnuunit, *) &
            'sgnlog10(x) = abs(x) > 10.0 ? log10(abs(x))*sgn(x) : x/10.0'
          ! merge: Fortran ternary operator: similar to ? : in C
          zmin = minval(color)
          zmin = merge(sign(log10(abs(zmin)), zmin), zmin/10.0_dp, abs(zmin) > 10.0_dp)
          zmax = maxval(color)
          zmax = merge(sign(log10(abs(zmax)), zmax), zmax/10.0_dp, abs(zmax) > 10.0_dp)
          write (gnuunit, *) &
            'set palette defined (', zmin, ' "blue", 0 "green", ', zmax, ' "red")'
          write (gnuunit, *) 'set pm3d map'
          write (gnuunit, *) 'set zrange [', zmin, ':', zmax, ']'
          write (gnuunit, *) 'splot ', '"'//trim(seedname)//'-bands.dat', &
            '" u 1:2:(sgnlog10($3)) with dots palette'
        end if
        close (gnuunit)
        !
        ! python script
        !
        file_name = trim(seedname)//'-bands.py'
        write (stdout, '(/,3x,a)') file_name
        open (newunit=pyunit, file=file_name, form='formatted')
        write (pyunit, '(a)') 'import pylab as pl'
        write (pyunit, '(a)') 'import numpy as np'
        write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
          "-bands.dat')"
        write (pyunit, '(a)') "x=data[:,0]"
        write (pyunit, '(a)') "y=data[:,1]"
        if (pw90_kpath%bands_colour == 'spin' &
            .or. pw90_kpath%bands_colour == 'shc') write (pyunit, '(a)') "z=data[:,2]"
        if (pw90_kpath%bands_colour == 'shc') write (pyunit, '(a)') &
          "z=np.array([np.log10(abs(elem))*np.sign(elem) " &
          //"if abs(elem)>10 else elem/10.0 for elem in z])"
        write (pyunit, '(a)') "tick_labels=[]"
        write (pyunit, '(a)') "tick_locs=[]"
        do j = 1, num_spts
          if (trim(glabel(j)) == ' G') then
            write (pyunit, '(a)') "tick_labels.append('$\Gamma$')"
          else
            write (pyunit, '(a)') "tick_labels.append('"//trim(glabel(j)) &
              //"'.strip())"
          end if
          if (j == 1) then
            write (pyunit, '(a,F12.6,a)') "tick_locs.append(0)"
          else
            write (pyunit, '(a,F12.6,a)') "tick_locs.append(", &
              sum(kpath_len(1:j - 1)), ")"
          end if
        end do
        if (pw90_kpath%bands_colour == 'none') then
          write (pyunit, '(a)') "pl.scatter(x,y,color='k',marker='+',s=0.1)"
        else if (pw90_kpath%bands_colour == 'spin' .or. &
                 pw90_kpath%bands_colour == 'shc') then
          write (pyunit, '(a)') &
            "pl.scatter(x,y,c=z,marker='+',s=1,cmap=pl.cm.jet)"
        end if
        write (pyunit, '(a)') "pl.xlim([0,max(x)])"
        write (pyunit, '(a)') "pl.ylim([min(y)-0.025*(max(y)-min(y))," &
          //"max(y)+0.025*(max(y)-min(y))])"
        write (pyunit, '(a)') "pl.xticks(tick_locs,tick_labels)"
        write (pyunit, '(a)') "for n in range(1,len(tick_locs)):"
        write (pyunit, '(a)') "   pl.plot([tick_locs[n],tick_locs[n]]," &
          //"[pl.ylim()[0],pl.ylim()[1]],color='gray'," &
          //"linestyle='-',linewidth=0.5)"
        write (pyunit, '(a)') "pl.ylabel('Energy [eV]')"
        if (pw90_kpath%bands_colour == 'spin' .or. pw90_kpath%bands_colour == 'shc') then
          write (pyunit, '(a)') &
            "pl.axes().set_aspect(aspect=0.65*max(x)/(max(y)-min(y)))"
          write (pyunit, '(a)') "pl.colorbar(shrink=0.7)"
        end if
        write (pyunit, '(a)') "outfile = '"//trim(seedname)//"-bands.pdf'"
        write (pyunit, '(a)') "pl.savefig(outfile,bbox_inches='tight')"
        write (pyunit, '(a)') "pl.show()"

      end if ! plot_bands .and. .not.plot_curv .and. .not.plot_morb .and. .not. plot_shc

      if (plot_curv) then
        ! It is conventional to plot the negative curvature
        curv = -curv
        file_name = trim(seedname)//'-curv.dat'
        write (stdout, '(/,3x,a)') file_name
        open (newunit=dataunit, file=file_name, form='formatted')
        do loop_kpt = 1, total_pts
          write (dataunit, '(4E16.8)') xval(loop_kpt), &
            curv(loop_kpt, :)
        end do
        write (dataunit, *) ' '
        close (dataunit)
      end if

      if (plot_curv .and. .not. plot_bands) then

        do i = 1, 3
          !
          ! gnuplot script
          !
          file_name = trim(seedname)//'-curv_'//achar(119 + i)//'.gnu'
          write (stdout, '(/,3x,a)') file_name
          open (newunit=gnuunit, file=file_name, form='formatted')
          ymin = minval(curv(:, i))
          ymax = maxval(curv(:, i))
          range = ymax - ymin
          ymin = ymin - 0.02_dp*range
          ymax = ymax + 0.02_dp*range
          write (gnuunit, 707) xval(total_pts), ymin, ymax
          do j = 1, num_paths - 1
            write (gnuunit, 705) sum(kpath_len(1:j)), ymin, &
              sum(kpath_len(1:j)), ymax
          end do
          write (gnuunit, 702, advance="no") glabel(1), 0.0_dp, &
            (glabel(j + 1), sum(kpath_len(1:j)), j=1, num_paths - 1)
          write (gnuunit, 703) glabel(1 + num_paths), sum(kpath_len(:))
          write (gnuunit, *) &
            'plot ', '"'//trim(seedname)//'-curv.dat', '" u 1:'//achar(49 + i)
          close (gnuunit)
          !
          ! python script
          !
          file_name = trim(seedname)//'-curv_'//achar(119 + i)//'.py'
          write (stdout, '(/,3x,a)') file_name
          open (newunit=pyunit, file=file_name, form='formatted')
          write (pyunit, '(a)') 'import pylab as pl'
          write (pyunit, '(a)') 'import numpy as np'
          write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
            "-curv.dat')"
          write (pyunit, '(a)') "x=data[:,0]"
          write (pyunit, '(a)') "y=data[:,"//achar(48 + i)//"]"
          write (pyunit, '(a)') "tick_labels=[]"
          write (pyunit, '(a)') "tick_locs=[]"
          do j = 1, num_spts
            if (trim(glabel(j)) == ' G') then
              write (pyunit, '(a)') "tick_labels.append('$\Gamma$')"
            else
              write (pyunit, '(a)') "tick_labels.append('" &
                //trim(glabel(j))//"'.strip())"
            end if
            if (j == 1) then
              write (pyunit, '(a,F12.6,a)') "tick_locs.append(0)"
            else
              write (pyunit, '(a,F12.6,a)') "tick_locs.append(", &
                sum(kpath_len(1:j - 1)), ")"
            end if
          end do
          write (pyunit, '(a)') "pl.plot(x,y,color='k')"
          write (pyunit, '(a)') "pl.xlim([0,max(x)])"
          write (pyunit, '(a)') "pl.ylim([min(y)-0.025*(max(y)-min(y))," &
            //"max(y)+0.025*(max(y)-min(y))])"
          write (pyunit, '(a)') "pl.xticks(tick_locs,tick_labels)"
          write (pyunit, '(a)') "for n in range(1,len(tick_locs)):"
          write (pyunit, '(a)') "   pl.plot([tick_locs[n],tick_locs[n]]," &
            //"[pl.ylim()[0],pl.ylim()[1]],color='gray'," &
            //"linestyle='-',linewidth=0.5)"
          if (pw90_berry%curv_unit == 'ang2') then
            write (pyunit, '(a)') "pl.ylabel('$-\Omega_"//achar(119 + i) &
              //"(\mathbf{k})$  [ $\AA^2$ ]')"
          else if (pw90_berry%curv_unit == 'bohr2') then
            write (pyunit, '(a)') "pl.ylabel('$-\Omega_"//achar(119 + i) &
              //"(\mathbf{k})$  [ bohr$^2$ ]')"
          end if
          write (pyunit, '(a)') "outfile = '"//trim(seedname)// &
            "-curv_"//achar(119 + i)//".pdf'"
          write (pyunit, '(a)') "pl.savefig(outfile,bbox_inches='tight')"
          write (pyunit, '(a)') "pl.show()"
        end do

      end if ! plot_curv .and. .not.plot_bands

      if (plot_morb) then
        file_name = trim(seedname)//'-morb.dat'
        write (stdout, '(/,3x,a)') file_name
        open (newunit=dataunit, file=file_name, form='formatted')
        do loop_kpt = 1, total_pts
          write (dataunit, '(4E16.8)') xval(loop_kpt), morb(loop_kpt, :)
        end do
        write (dataunit, *) ' '
        close (dataunit)
      end if

      if (plot_morb .and. .not. plot_bands) then
        do i = 1, 3
          !
          ! gnuplot script
          !
          file_name = trim(seedname)//'-morb_'//achar(119 + i)//'.gnu'
          write (stdout, '(/,3x,a)') file_name
          open (newunit=gnuunit, file=file_name, form='formatted')
          ymin = minval(morb(:, i))
          ymax = maxval(morb(:, i))
          range = ymax - ymin
          ymin = ymin - 0.02_dp*range
          ymax = ymax + 0.02_dp*range
          write (gnuunit, 707) xval(total_pts), ymin, ymax
          do j = 1, num_paths - 1
            write (gnuunit, 705) sum(kpath_len(1:j)), ymin, &
              sum(kpath_len(1:j)), ymax
          end do
          write (gnuunit, 702, advance="no") glabel(1), 0.0_dp, &
            (glabel(j + 1), sum(kpath_len(1:j)), j=1, num_paths - 1)
          write (gnuunit, 703) glabel(1 + num_paths), sum(kpath_len(:))
          write (gnuunit, *) &
            'plot ', '"'//trim(seedname)//'-morb.dat', '" u 1:'//achar(49 + i)
          close (gnuunit)
          !
          ! python script
          !
          file_name = trim(seedname)//'-morb_'//achar(119 + i)//'.py'
          write (stdout, '(/,3x,a)') file_name
          open (newunit=pyunit, file=file_name, form='formatted')
          write (pyunit, '(a)') 'import pylab as pl'
          write (pyunit, '(a)') 'import numpy as np'
          write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
            "-morb.dat')"
          write (pyunit, '(a)') "x=data[:,0]"
          write (pyunit, '(a)') "y=data[:,"//achar(48 + i)//"]"
          write (pyunit, '(a)') "tick_labels=[]"
          write (pyunit, '(a)') "tick_locs=[]"
          do j = 1, num_spts
            if (trim(glabel(j)) == ' G') then
              write (pyunit, '(a)') "tick_labels.append('$\Gamma$')"
            else
              write (pyunit, '(a)') &
                "tick_labels.append('"//trim(glabel(j))//"'.strip())"
            end if
            if (j == 1) then
              write (pyunit, '(a,F12.6,a)') "tick_locs.append(0)"
            else
              write (pyunit, '(a,F12.6,a)') "tick_locs.append(", &
                sum(kpath_len(1:j - 1)), ")"
            end if
          end do
          write (pyunit, '(a)') "pl.plot(x,y,color='k')"
          write (pyunit, '(a)') "pl.xlim([0,max(x)])"
          write (pyunit, '(a)') "pl.ylim([min(y)-0.025*(max(y)-min(y))," &
            //"max(y)+0.025*(max(y)-min(y))])"
          write (pyunit, '(a)') "pl.xticks(tick_locs,tick_labels)"
          write (pyunit, '(a)') "for n in range(1,len(tick_locs)):"
          write (pyunit, '(a)') "   pl.plot([tick_locs[n],tick_locs[n]]," &
            //"[pl.ylim()[0],pl.ylim()[1]],color='gray'," &
            //"linestyle='-',linewidth=0.5)"
          write (pyunit, '(a)') "pl.ylabel(r'$M^{\rm{orb}}_z(\mathbf{k})$" &
            //"  [ Ry$\cdot\AA^2$ ]')"
          write (pyunit, '(a)') "outfile = '"//trim(seedname)// &
            "-morb_"//achar(119 + i)//".pdf'"
          write (pyunit, '(a)') "pl.savefig(outfile,bbox_inches='tight')"
          write (pyunit, '(a)') "pl.show()"
        end do

      end if ! plot_morb .and. .not.plot_bands

      if (plot_shc) then
        file_name = trim(seedname)//'-shc.dat'
        write (stdout, '(/,3x,a)') file_name
        open (newunit=dataunit, file=file_name, form='formatted')
        do loop_kpt = 1, total_pts
          write (dataunit, '(2E16.8)') xval(loop_kpt), &
            shc(loop_kpt)
        end do
        write (dataunit, *) ' '
        close (dataunit)
      end if

      if (plot_shc .and. .not. plot_bands) then
        !
        ! gnuplot script
        !
        file_name = trim(seedname)//'-shc'//'.gnu'
        write (stdout, '(/,3x,a)') file_name
        open (newunit=gnuunit, file=file_name, form='formatted')
        ymin = minval(shc(:))
        ymax = maxval(shc(:))
        range = ymax - ymin
        ymin = ymin - 0.02_dp*range
        ymax = ymax + 0.02_dp*range
        write (gnuunit, 707) xval(total_pts), ymin, ymax
        do j = 1, num_paths - 1
          write (gnuunit, 705) sum(kpath_len(1:j)), ymin, &
            sum(kpath_len(1:j)), ymax
        end do
        write (gnuunit, 702, advance="no") glabel(1), 0.0_dp, &
          (glabel(j + 1), sum(kpath_len(1:j)), j=1, num_paths - 1)
        write (gnuunit, 703) glabel(1 + num_paths), sum(kpath_len(:))
        write (gnuunit, *) &
          'plot ', '"'//trim(seedname)//'-shc.dat', '" u 1:2'
        close (gnuunit)
        !
        ! python script
        !
        file_name = trim(seedname)//'-shc'//'.py'
        write (stdout, '(/,3x,a)') file_name
        open (newunit=pyunit, file=file_name, form='formatted')
        write (pyunit, '(a)') "# uncomment these two lines if you are " &
          //"running in non-GUI environment"
        write (pyunit, '(a)') "#import matplotlib"
        write (pyunit, '(a)') "#matplotlib.use('Agg')"
        write (pyunit, '(a)') 'import pylab as pl'
        write (pyunit, '(a)') 'import numpy as np'
        write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
          "-shc.dat')"
        write (pyunit, '(a)') "x=data[:,0]"
        write (pyunit, '(a)') "y=data[:,1]"
        write (pyunit, '(a)') "tick_labels=[]"
        write (pyunit, '(a)') "tick_locs=[]"
        do j = 1, num_spts
          if (trim(glabel(j)) == ' G') then
            write (pyunit, '(a)') "tick_labels.append('$\Gamma$')"
          else
            write (pyunit, '(a)') "tick_labels.append('" &
              //trim(glabel(j))//"'.strip())"
          end if
          if (j == 1) then
            write (pyunit, '(a,F12.6,a)') "tick_locs.append(0)"
          else
            write (pyunit, '(a,F12.6,a)') "tick_locs.append(", &
              sum(kpath_len(1:j - 1)), ")"
          end if
        end do
        write (pyunit, '(a)') "pl.plot(x,y,color='k')"
        write (pyunit, '(a)') "pl.xlim([0,max(x)])"
        write (pyunit, '(a)') "pl.ylim([min(y)-0.025*(max(y)-min(y))," &
          //"max(y)+0.025*(max(y)-min(y))])"
        write (pyunit, '(a)') "pl.xticks(tick_locs,tick_labels)"
        write (pyunit, '(a)') "for n in range(1,len(tick_locs)):"
        write (pyunit, '(a)') "   pl.plot([tick_locs[n],tick_locs[n]]," &
          //"[pl.ylim()[0],pl.ylim()[1]],color='gray'," &
          //"linestyle='-',linewidth=0.5)"
        if (pw90_berry%curv_unit == 'ang2') then
          write (pyunit, '(a)') "pl.ylabel('$\Omega_{" &
            //achar(119 + pw90_spin_hall%alpha)//achar(119 + pw90_spin_hall%beta) &
            //"}^{spin"//achar(119 + pw90_spin_hall%gamma) &
            //"}(\mathbf{k})$  [ $\AA^2$ ]')"
        else if (pw90_berry%curv_unit == 'bohr2') then
          write (pyunit, '(a)') "pl.ylabel('$\Omega_{" &
            //achar(119 + pw90_spin_hall%alpha)//achar(119 + pw90_spin_hall%beta) &
            //"}^{spin"//achar(119 + pw90_spin_hall%gamma) &
            //"}(\mathbf{k})$  [ bohr$^2$ ]')"
        end if
        write (pyunit, '(a)') "outfile = '"//trim(seedname)// &
          "-shc"//".pdf'"
        write (pyunit, '(a)') "pl.savefig(outfile,bbox_inches='tight')"
        write (pyunit, '(a)') "pl.show()"

      end if ! plot_shc .and. .not.plot_bands

      if (plot_bands .and. plot_shc) then
        !
        ! python script
        !
        file_name = trim(seedname)//'-bands+shc'//'.py'
        write (stdout, '(/,3x,a)') file_name
        open (newunit=pyunit, file=file_name, form='formatted')
        write (pyunit, '(a)') "# uncomment these two lines if you are " &
          //"running in non-GUI environment"
        write (pyunit, '(a)') "#import matplotlib"
        write (pyunit, '(a)') "#matplotlib.use('Agg')"
        write (pyunit, '(a)') 'import pylab as pl'
        write (pyunit, '(a)') 'import numpy as np'
        write (pyunit, '(a)') 'from matplotlib.gridspec import GridSpec'
        write (pyunit, '(a)') "tick_labels=[]"
        write (pyunit, '(a)') "tick_locs=[]"
        do j = 1, num_spts
          if (trim(glabel(j)) == ' G') then
            write (pyunit, '(a)') "tick_labels.append('$\Gamma$')"
          else
            write (pyunit, '(a)') "tick_labels.append('"//trim(glabel(j)) &
              //"'.strip())"
          end if
          if (j == 1) then
            write (pyunit, '(a,F12.6,a)') "tick_locs.append(0)"
          else
            write (pyunit, '(a,F12.6,a)') "tick_locs.append(", &
              sum(kpath_len(1:j - 1)), ")"
          end if
        end do
        write (pyunit, '(a)') "fig = pl.figure()"
        write (pyunit, '(a)') "gs = GridSpec(2,52,hspace=0.00,wspace=1)"
        !
        ! upper panel (energy bands)
        !
        write (pyunit, '(a)') "axes1 = pl.subplot(gs[0, :-2])"
        write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
          "-bands.dat')"
        write (pyunit, '(a)') "x=data[:,0]"
        write (pyunit, '(a,F12.6)') "y=data[:,1]-", fermi_energy_list(1)
        if (pw90_kpath%bands_colour == 'spin' .or. pw90_kpath%bands_colour == 'shc') &
          write (pyunit, '(a)') "z=data[:,2]"
        if (pw90_kpath%bands_colour == 'shc') &
          write (pyunit, '(a)') "z=np.array([np.log10(abs(elem))*np.sign(elem) " &
          //"if abs(elem)>10 else elem/10.0 for elem in z])"
        if (pw90_kpath%bands_colour == 'none') then
          write (pyunit, '(a)') "pl.scatter(x,y,color='k',marker='+',s=0.1)"
        else if (pw90_kpath%bands_colour == 'spin' &
                 .or. pw90_kpath%bands_colour == 'shc') then
          write (pyunit, '(a)') &
            "pl.scatter(x,y,c=z,marker='+',s=1,cmap=pl.cm.jet)"
        end if
        write (pyunit, '(a)') "pl.xlim([0,max(x)])"
        write (pyunit, '(a)') "pl.ylim([-0.65,0.65]) # Adjust this range as needed"
        write (pyunit, '(a)') "pl.plot([tick_locs[0],tick_locs[-1]],[0,0]," &
          //"color='black',linestyle='--',linewidth=0.5)"
        write (pyunit, '(a)') "pl.xticks(tick_locs,tick_labels)"
        write (pyunit, '(a)') "for n in range(1,len(tick_locs)):"
        write (pyunit, '(a)') "   pl.plot([tick_locs[n],tick_locs[n]]," &
          //"[pl.ylim()[0],pl.ylim()[1]],color='gray'," &
          //"linestyle='-',linewidth=0.5)"
        write (pyunit, '(a)') "pl.ylabel('Energy$-$E$_F$ [eV]')"
        write (pyunit, '(a)') "pl.tick_params(axis='x'," &
          //"which='both',bottom='off',top='off',labelbottom='off')"
        !
        ! color bar for upper panel
        !
        write (pyunit, '(a)') "# Now adding the colorbar"
        write (pyunit, '(a)') "cbaxes = pl.subplot(gs[:, -2:])"
        write (pyunit, '(a)') "cb = pl.colorbar(cax = cbaxes," &
          //"orientation='vertical') "
        write (pyunit, '(a)') "#cblim = int(min(abs(max(z)),abs(min(z))))"
        write (pyunit, '(a)') "#cb.set_ticks([-cblim,0,cblim])"
        write (pyunit, '(a)') "#cblim = [min(z),0,max(z)]"
        write (pyunit, '(a)') "#cb.set_ticks([min(z),0,max(z)])"
        write (pyunit, '(a)') "#cb.set_ticklabels(['-','0','+'])"
        !
        ! lower panel (SHC)
        !
        write (pyunit, '(a)') "axes2 = pl.subplot(gs[1, :-2])"
        write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
          "-shc.dat')"
        write (pyunit, '(a)') "x=data[:,0]"
        write (pyunit, '(a)') "y=data[:,1]"
        write (pyunit, '(a)') &
                "y=np.array([np.log10(abs(elem))*np.sign(elem) &
                        &if abs(elem)>10 else elem/10.0 for elem in y])"
        write (pyunit, '(a)') "pl.plot(x,y,color='k')"
        write (pyunit, '(a)') "pl.xlim([0,max(x)])"
        write (pyunit, '(a)') "pl.ylim([min(y)-0.025*(max(y)-min(y))," &
          //"max(y)+0.025*(max(y)-min(y))])"
        write (pyunit, '(a)') "pl.plot([0,max(x)],[0,0],color='black'," &
          //"linestyle='--',linewidth=0.5)"
        write (pyunit, '(a)') "pl.xticks(tick_locs,tick_labels)"
        write (pyunit, '(a)') "for n in range(1,len(tick_locs)):"
        write (pyunit, '(a)') "   pl.plot([tick_locs[n],tick_locs[n]]," &
          //"[pl.ylim()[0],pl.ylim()[1]],color='gray'," &
          //"linestyle='-',linewidth=0.5)"
        if (pw90_berry%curv_unit == 'ang2') then
          write (pyunit, '(a)') "pl.ylabel('$log_{10}|\Omega_{" &
            //achar(119 + pw90_spin_hall%alpha)//achar(119 + pw90_spin_hall%beta) &
            //"}^{spin"//achar(119 + pw90_spin_hall%gamma) &
            //"}(\mathbf{k})|$  [ $\AA^2$ ]')"
        else if (pw90_berry%curv_unit == 'bohr2') then
          write (pyunit, '(a)') "pl.ylabel('$\Omega_{" &
            //achar(119 + pw90_spin_hall%alpha)//achar(119 + pw90_spin_hall%beta) &
            //"}^{spin"//achar(119 + pw90_spin_hall%gamma) &
            //"}(\mathbf{k})$  [ bohr$^2$ ]')"
        end if
        write (pyunit, '(a)') "outfile = '"//trim(seedname)// &
          "-bands+shc"//".pdf'"
        write (pyunit, '(a)') "pl.savefig(outfile,bbox_inches='tight')"
        write (pyunit, '(a)') "pl.show()"

      end if ! plot_bands .and. plot_shc

      if (plot_bands .and. (plot_curv .or. plot_morb)) then
        !
        ! python script
        !
        do i = 1, 3
          if (plot_curv) then
            file_name = trim(seedname)//'-bands+curv_'//achar(119 + i)//'.py'
          else if (plot_morb) then
            file_name = trim(seedname)//'-bands+morb_'//achar(119 + i)//'.py'
          end if
          write (stdout, '(/,3x,a)') file_name
          open (newunit=pyunit, file=file_name, form='formatted')
          write (pyunit, '(a)') 'import pylab as pl'
          write (pyunit, '(a)') 'import numpy as np'
          write (pyunit, '(a)') 'from matplotlib.gridspec import GridSpec'
          write (pyunit, '(a)') "tick_labels=[]"
          write (pyunit, '(a)') "tick_locs=[]"
          do j = 1, num_spts
            if (trim(glabel(j)) == ' G') then
              write (pyunit, '(a)') "tick_labels.append('$\Gamma$')"
            else
              write (pyunit, '(a)') "tick_labels.append('"//trim(glabel(j)) &
                //"'.strip())"
            end if
            if (j == 1) then
              write (pyunit, '(a,F12.6,a)') "tick_locs.append(0)"
            else
              write (pyunit, '(a,F12.6,a)') "tick_locs.append(", &
                sum(kpath_len(1:j - 1)), ")"
            end if
          end do
          write (pyunit, '(a)') "fig = pl.figure()"
          write (pyunit, '(a)') "gs = GridSpec(2, 1,hspace=0.00)"
          !
          ! upper panel (energy bands)
          !
          write (pyunit, '(a)') "axes1 = pl.subplot(gs[0, 0:])"
          write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
            "-bands.dat')"
          write (pyunit, '(a)') "x=data[:,0]"
          write (pyunit, '(a,F12.6)') "y=data[:,1]-", fermi_energy_list(1)
          if (pw90_kpath%bands_colour == 'spin') write (pyunit, '(a)') "z=data[:,2]"
          if (pw90_kpath%bands_colour == 'shc') then
            write (pyunit, '(a)') "z=data[:,2]"
            write (pyunit, '(a)') "z=np.array([np.log10(abs(elem))*np.sign(elem) " &
              //"if abs(elem)>10 else elem/10.0 for elem in z])"
          end if
          if (pw90_kpath%bands_colour == 'none') then
            write (pyunit, '(a)') "pl.scatter(x,y,color='k',marker='+',s=0.1)"
          else if (pw90_kpath%bands_colour == 'spin' .or. pw90_kpath%bands_colour == 'shc') then
            write (pyunit, '(a)') &
              "pl.scatter(x,y,c=z,marker='+',s=1,cmap=pl.cm.jet)"
          end if
          write (pyunit, '(a)') "pl.xlim([0,max(x)])"
          write (pyunit, '(a)') "pl.ylim([-0.65,0.65]) # Adjust this range as needed"
          write (pyunit, '(a)') "pl.plot([tick_locs[0],tick_locs[-1]],[0,0]," &
            //"color='black',linestyle='--',linewidth=0.5)"
          write (pyunit, '(a)') "pl.xticks(tick_locs,tick_labels)"
          write (pyunit, '(a)') "for n in range(1,len(tick_locs)):"
          write (pyunit, '(a)') "   pl.plot([tick_locs[n],tick_locs[n]]," &
            //"[pl.ylim()[0],pl.ylim()[1]],color='gray'," &
            //"linestyle='-',linewidth=0.5)"
          write (pyunit, '(a)') "pl.ylabel('Energy$-$E$_F$ [eV]')"
          write (pyunit, '(a)') "pl.tick_params(axis='x'," &
            //"which='both',bottom='off',top='off',labelbottom='off')"
          !
          ! lower panel (curvature or orbital magnetization)
          !
          write (pyunit, '(a)') "axes2 = pl.subplot(gs[1, 0:])"
          if (plot_curv) then
            write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
              "-curv.dat')"
          else if (plot_morb) then
            write (pyunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
              "-morb.dat')"
          end if
          write (pyunit, '(a)') "x=data[:,0]"
          write (pyunit, '(a)') "y=data[:,"//achar(48 + i)//"]"
          write (pyunit, '(a)') "pl.plot(x,y,color='k')"
          write (pyunit, '(a)') "pl.xlim([0,max(x)])"
          write (pyunit, '(a)') "pl.ylim([min(y)-0.025*(max(y)-min(y))," &
            //"max(y)+0.025*(max(y)-min(y))])"
          write (pyunit, '(a)') "pl.xticks(tick_locs,tick_labels)"
          write (pyunit, '(a)') "for n in range(1,len(tick_locs)):"
          write (pyunit, '(a)') "   pl.plot([tick_locs[n],tick_locs[n]]," &
            //"[pl.ylim()[0],pl.ylim()[1]],color='gray'," &
            //"linestyle='-',linewidth=0.5)"
          if (plot_curv) then
            if (pw90_berry%curv_unit == 'ang2') then
              write (pyunit, '(a)') "pl.ylabel('$-\Omega_"//achar(119 + i) &
                //"(\mathbf{k})$  [ $\AA^2$ ]')"
            else if (pw90_berry%curv_unit == 'bohr2') then
              write (pyunit, '(a)') "pl.ylabel('$-\Omega_"//achar(119 + i) &
                //"(\mathbf{k})$  [ bohr$^2$ ]')"
            end if
            write (pyunit, '(a)') "outfile = '"//trim(seedname)// &
              "-bands+curv_"//achar(119 + i)//".pdf'"
          else if (plot_morb) then
            write (pyunit, '(a)') "pl.ylabel(r'$M^{\rm{orb}}_z(\mathbf{k})$" &
              //"  [ Ry$\cdot\AA^2$ ]')"
            write (pyunit, '(a)') "outfile = '"//trim(seedname)// &
              "-morb_"//achar(119 + i)//".pdf'"
          end if
          write (pyunit, '(a)') "pl.savefig(outfile,bbox_inches='tight')"
          write (pyunit, '(a)') "pl.show()"
        end do

      end if ! plot_bands .and. plot_curv

    end if ! on_root

701 format('set style data dots', /, 'unset key', /, &
           'set xrange [0:', F8.5, ']', /, 'set yrange [', F16.8, ' :', F16.8, ']')
702 format('set xtics (', :20('"', A3, '" ', F8.5, ','))
703 format(A3, '" ', F8.5, ')')
!704 format('set palette defined (', F8.5, ' "red", 0 "green", ', F8.5, ' "blue")') !not used
705 format('set arrow from ', F16.8, ',', F16.8, ' to ', F16.8, ',', F16.8, ' nohead')
706 format('unset key', /, &
           'set xrange [0:', F9.5, ']', /, 'set yrange [', F16.8, ' :', F16.8, ']')
707 format('set style data lines', /, 'set nokey', /, &
           'set xrange [0:', F8.5, ']', /, 'set yrange [', F16.8, ' :', F16.8, ']')

  end subroutine k_path

  !================================================!
  !                   PRIVATE PROCEDURES                      !
  !================================================!
  subroutine k_path_print_info(plot_bands, plot_curv, plot_morb, plot_shc, fermi_energy_list, &
                               pw90_kpath, berry_curv_unit, stdout, error, comm)
    !================================================!

    use w90_postw90_types, only: pw90_kpath_mod_type
    use w90_comms, only: w90_comm_type, mpirank

    ! arguments
    real(kind=dp), allocatable, intent(in) :: fermi_energy_list(:)
    type(pw90_kpath_mod_type), intent(in) :: pw90_kpath
    type(w90_comm_type), intent(in) :: comm
    type(w90_error_type), allocatable, intent(out) :: error
    integer, intent(in) :: stdout
    character(len=*), intent(in) :: berry_curv_unit
    logical, intent(in) :: plot_bands, plot_curv, plot_morb, plot_shc

    ! local variables
    integer :: fermi_n
    logical :: on_root = .false.
    if (mpirank(comm) == 0) on_root = .true.

    if (on_root) then
      write (stdout, '(/,/,1x,a)') &
        'Properties calculated in module  k p a t h'
      write (stdout, '(1x,a)') &
        '------------------------------------------'

      if (plot_bands) then
        select case (pw90_kpath%bands_colour)
        case ("none")
          write (stdout, '(/,3x,a)') '* Energy bands in eV'
        case ("spin")
          write (stdout, '(/,3x,a)') '* Energy bands in eV, coloured by spin'
        case ("shc")
          write (stdout, '(/,3x,a)') '* Energy bands in eV, coloured by SHC'
        end select
      end if
      fermi_n = 0
      if (allocated(fermi_energy_list)) fermi_n = size(fermi_energy_list)
      if (plot_curv) then
        if (berry_curv_unit == 'ang2') then
          write (stdout, '(/,3x,a)') '* Negative Berry curvature in Ang^2'
        else if (berry_curv_unit == 'bohr2') then
          write (stdout, '(/,3x,a)') '* Negative Berry curvature in Bohr^2'
        end if
        if (fermi_n /= 1) then
          call set_error_input(error, 'Must specify one Fermi level when kpath_task=curv', comm)
          return
        end if
      end if
      if (plot_morb) then
        write (stdout, '(/,3x,a)') &
          '* Orbital magnetization k-space integrand in eV.Ang^2'
        if (fermi_n /= 1) then
          call set_error_input(error, 'Must specify one Fermi level when kpath_task=morb', comm)
          return
        end if
      end if
      if (plot_shc) then
        if (berry_curv_unit == 'ang2') then
          write (stdout, '(/,3x,a)') '* Berry curvature-like term for' &
            //' spin Hall conductivity in Ang^2'
        else if (berry_curv_unit == 'bohr2') then
          write (stdout, '(/,3x,a)') '* Berry curvature-like term for' &
            //' spin Hall conductivity in Bohr^2'
        end if
        if (fermi_n /= 1) then
          call set_error_input(error, 'Must specify one Fermi level when kpath_task=shc', comm)
          return
        end if
      end if
    end if ! on_root

  end subroutine

  !================================================!
  subroutine k_path_get_points(num_paths, kpath_len, total_pts, xval, plot_kpoint, kpoint_path, &
                               recip_lattice, pw90_kpath)
    !================================================!
    ! Determine the number of k-points (total_pts) as well as           !
    ! their reciprocal-lattice coordinates long the path (plot_kpoint)  !
    ! and their associated horizontal coordinate for the plot (xval)    !
    !================================================!

    use w90_postw90_types, only: pw90_kpath_mod_type
    use w90_types, only: kpoint_path_type
    use w90_utility, only: utility_metric

    ! arguments
    type(pw90_kpath_mod_type), intent(in) :: pw90_kpath
    type(kpoint_path_type), intent(in) :: kpoint_path
    integer, intent(out) :: num_paths, total_pts
    real(kind=dp), allocatable, intent(out)   :: kpath_len(:), xval(:)
    real(kind=dp), allocatable, intent(out) :: plot_kpoint(:, :)
    real(kind=dp), intent(in) :: recip_lattice(3, 3)

    ! local variables
    integer :: counter, loop_path, loop_i
    integer, allocatable :: kpath_pts(:)
    real(kind=dp) :: vec(3)
    real(kind=dp) :: recip_metric(3, 3)

    call utility_metric(recip_lattice, recip_metric)
    ! Work out how many points there are in the total path, and the
    ! positions of the special points
    !
    if (allocated(kpoint_path%labels)) then
      !num_paths = kpoint_path%bands_num_spec_points/2 ! number of straight line segments
      num_paths = size(kpoint_path%labels)/2
    else
      num_paths = 0
    end if
    allocate (kpath_pts(num_paths))
    allocate (kpath_len(num_paths))
    do loop_path = 1, num_paths
      vec = kpoint_path%points(:, 2*loop_path) &
            - kpoint_path%points(:, 2*loop_path - 1)
      kpath_len(loop_path) = &
        sqrt(dot_product(vec, (matmul(recip_metric, vec))))
      !
      ! kpath_pts(loop_path) is the number of points in path number
      ! loop_path (all segments have the same density of points)
      !
      if (loop_path == 1) then
        kpath_pts(loop_path) = pw90_kpath%num_points
      else
        kpath_pts(loop_path) = nint(real(pw90_kpath%num_points, dp) &
                                    *kpath_len(loop_path)/kpath_len(1))
      end if
    end do
    total_pts = sum(kpath_pts) + 1

    ! Reciprocal-lattice coordinates of the k-points along the path
    !
    allocate (plot_kpoint(3, total_pts))

    ! Value of the horizontal coordinate in the actual plots (units of
    ! distance in k-space)
    !
    allocate (xval(total_pts))

    ! Find the position of each kpoint along the path
    !
    counter = 0
    do loop_path = 1, num_paths
      do loop_i = 1, kpath_pts(loop_path)
        counter = counter + 1
        if (counter == 1) then
          xval(counter) = 0.0_dp
        else
          xval(counter) = xval(counter - 1) &
                          + kpath_len(loop_path)/real(kpath_pts(loop_path), dp)
        end if
        plot_kpoint(:, counter) = kpoint_path%points(:, 2*loop_path - 1) &
                                  + (kpoint_path%points(:, 2*loop_path) &
                                     - kpoint_path%points(:, 2*loop_path - 1) &
                                     ) &
                                  *(real(loop_i - 1, dp)/real(kpath_pts(loop_path), dp))
      end do
    end do
    !
    ! Last point
    !
    xval(total_pts) = sum(kpath_len)
    plot_kpoint(:, total_pts) = kpoint_path%points(:, num_paths*2)

  end subroutine

end module w90_kpath