k_path Subroutine

public 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)

Uses

  • proc~~k_path~~UsesGraph proc~k_path k_path module~w90_berry w90_berry proc~k_path->module~w90_berry module~w90_comms w90_comms proc~k_path->module~w90_comms module~w90_constants w90_constants proc~k_path->module~w90_constants module~w90_get_oper w90_get_oper proc~k_path->module~w90_get_oper module~w90_io w90_io proc~k_path->module~w90_io module~w90_postw90_common w90_postw90_common proc~k_path->module~w90_postw90_common module~w90_postw90_types w90_postw90_types proc~k_path->module~w90_postw90_types module~w90_spin w90_spin proc~k_path->module~w90_spin module~w90_types w90_types proc~k_path->module~w90_types module~w90_utility w90_utility proc~k_path->module~w90_utility module~w90_berry->module~w90_constants module~w90_error w90_error module~w90_berry->module~w90_error 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_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_spin->module~w90_constants module~w90_spin->module~w90_error module~w90_types->module~w90_constants module~w90_utility->module~w90_comms module~w90_utility->module~w90_constants module~w90_error->module~w90_comms module~w90_error->module~w90_error_base

Main routine

Arguments

Type IntentOptional Attributes Name
type(pw90_berry_mod_type), intent(in) :: pw90_berry
type(dis_manifold_type), intent(in) :: dis_manifold
real(kind=dp), intent(in), allocatable :: fermi_energy_list(:)
type(kmesh_info_type), intent(in) :: kmesh_info
type(pw90_kpath_mod_type), intent(in) :: pw90_kpath
real(kind=dp), intent(in) :: kpt_latt(:,:)
type(pw90_oper_read_type), intent(in) :: pw90_oper_read
type(pw90_band_deriv_degen_type), intent(in) :: pw90_band_deriv_degen
type(pw90_spin_mod_type), intent(in) :: pw90_spin
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(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 :: AA_R(:,:,:,:)
complex(kind=dp), intent(inout), allocatable :: BB_R(:,:,:,:)
complex(kind=dp), intent(inout), allocatable :: CC_R(:,:,:,:,:)
complex(kind=dp), intent(inout), allocatable :: HH_R(:,:,:)
complex(kind=dp), intent(inout), allocatable :: SH_R(:,:,:,:)
complex(kind=dp), intent(inout), allocatable :: SHR_R(:,:,:,:,:)
complex(kind=dp), intent(inout), allocatable :: SR_R(:,:,:,:,:)
complex(kind=dp), intent(inout), allocatable :: SS_R(:,:,:,:)
complex(kind=dp), intent(inout), allocatable :: SAA_R(:,:,:,:,:)
complex(kind=dp), intent(inout), allocatable :: SBB_R(:,:,:,:,:)
complex(kind=dp), intent(in) :: v_matrix(:,:,:)
complex(kind=dp), intent(in) :: 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
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: fermi_n
integer, intent(in) :: num_wann
integer, intent(in) :: num_bands
integer, intent(in) :: num_kpts
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~~k_path~~CallsGraph proc~k_path k_path 20 20 proc~k_path->20 interface~comms_bcast comms_bcast proc~k_path->interface~comms_bcast interface~comms_gatherv comms_gatherv proc~k_path->interface~comms_gatherv interface~comms_scatterv comms_scatterv proc~k_path->interface~comms_scatterv proc~berry_get_imf_klist berry_get_imf_klist proc~k_path->proc~berry_get_imf_klist proc~berry_get_imfgh_klist berry_get_imfgh_klist proc~k_path->proc~berry_get_imfgh_klist proc~berry_get_shc_klist berry_get_shc_klist proc~k_path->proc~berry_get_shc_klist proc~comms_array_split comms_array_split proc~k_path->proc~comms_array_split proc~get_aa_r get_AA_R proc~k_path->proc~get_aa_r proc~get_aa_r_effective get_AA_R_effective proc~k_path->proc~get_aa_r_effective proc~get_bb_r get_BB_R proc~k_path->proc~get_bb_r proc~get_cc_r get_CC_R proc~k_path->proc~get_cc_r proc~get_hh_r get_HH_R proc~k_path->proc~get_hh_r proc~get_shc_r get_SHC_R proc~k_path->proc~get_shc_r proc~get_ss_r get_SS_R proc~k_path->proc~get_ss_r proc~k_path_get_points k_path_get_points proc~k_path->proc~k_path_get_points proc~k_path_print_info k_path_print_info proc~k_path->proc~k_path_print_info proc~mpirank mpirank proc~k_path->proc~mpirank proc~mpisize mpisize proc~k_path->proc~mpisize proc~pw90common_fourier_r_to_k pw90common_fourier_R_to_k proc~k_path->proc~pw90common_fourier_r_to_k proc~set_error_input set_error_input proc~k_path->proc~set_error_input proc~spin_get_nk spin_get_nk proc~k_path->proc~spin_get_nk proc~utility_diagonalize utility_diagonalize proc~k_path->proc~utility_diagonalize proc~utility_recip_lattice_base utility_recip_lattice_base proc~k_path->proc~utility_recip_lattice_base 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~berry_get_imf_klist->proc~berry_get_imfgh_klist proc~berry_get_imfgh_klist->proc~pw90common_fourier_r_to_k proc~pw90common_fourier_r_to_k_vec pw90common_fourier_R_to_k_vec proc~berry_get_imfgh_klist->proc~pw90common_fourier_r_to_k_vec proc~utility_im_tr_prod utility_im_tr_prod proc~berry_get_imfgh_klist->proc~utility_im_tr_prod proc~utility_re_tr_prod utility_re_tr_prod proc~berry_get_imfgh_klist->proc~utility_re_tr_prod proc~utility_zgemm_new utility_zgemm_new proc~berry_get_imfgh_klist->proc~utility_zgemm_new proc~wham_get_eig_uu_hh_jjlist wham_get_eig_UU_HH_JJlist proc~berry_get_imfgh_klist->proc~wham_get_eig_uu_hh_jjlist proc~wham_get_occ_mat_list wham_get_occ_mat_list proc~berry_get_imfgh_klist->proc~wham_get_occ_mat_list proc~berry_get_shc_klist->proc~utility_recip_lattice_base interface~pw90common_kmesh_spacing pw90common_kmesh_spacing proc~berry_get_shc_klist->interface~pw90common_kmesh_spacing proc~pw90common_fourier_r_to_k_new pw90common_fourier_R_to_k_new proc~berry_get_shc_klist->proc~pw90common_fourier_r_to_k_new proc~berry_get_shc_klist->proc~pw90common_fourier_r_to_k_vec proc~pw90common_get_occ pw90common_get_occ proc~berry_get_shc_klist->proc~pw90common_get_occ proc~utility_rotate utility_rotate proc~berry_get_shc_klist->proc~utility_rotate proc~wham_get_d_h wham_get_D_h proc~berry_get_shc_klist->proc~wham_get_d_h proc~wham_get_eig_deleig wham_get_eig_deleig proc~berry_get_shc_klist->proc~wham_get_eig_deleig proc~comms_array_split->proc~mpisize proc~get_aa_r->interface~comms_bcast proc~get_aa_r->interface~comms_scatterv proc~get_aa_r->proc~comms_array_split proc~get_aa_r->proc~mpirank proc~get_aa_r->proc~mpisize interface~comms_reduce comms_reduce proc~get_aa_r->interface~comms_reduce proc~fourier_loc_q_to_r fourier_loc_q_to_R proc~get_aa_r->proc~fourier_loc_q_to_r proc~get_gauge_overlap_matrix get_gauge_overlap_matrix proc~get_aa_r->proc~get_gauge_overlap_matrix proc~io_stopwatch_start io_stopwatch_start proc~get_aa_r->proc~io_stopwatch_start proc~io_stopwatch_stop io_stopwatch_stop proc~get_aa_r->proc~io_stopwatch_stop proc~operator_wigner_setup operator_wigner_setup proc~get_aa_r->proc~operator_wigner_setup proc~set_error_fatal set_error_fatal proc~get_aa_r->proc~set_error_fatal proc~set_error_file set_error_file proc~get_aa_r->proc~set_error_file proc~get_aa_r_effective->interface~comms_bcast proc~get_aa_r_effective->proc~mpirank proc~get_aa_r_effective->proc~io_stopwatch_start proc~get_aa_r_effective->proc~io_stopwatch_stop proc~get_aa_r_effective->proc~set_error_fatal proc~get_aa_r_effective->proc~set_error_file proc~get_bb_r->interface~comms_bcast proc~get_bb_r->interface~comms_scatterv proc~get_bb_r->proc~comms_array_split proc~get_bb_r->proc~mpirank proc~get_bb_r->proc~mpisize proc~get_bb_r->interface~comms_reduce proc~get_bb_r->proc~fourier_loc_q_to_r proc~get_bb_r->proc~get_gauge_overlap_matrix proc~get_win_min get_win_min proc~get_bb_r->proc~get_win_min proc~get_bb_r->proc~io_stopwatch_start proc~get_bb_r->proc~io_stopwatch_stop proc~get_bb_r->proc~operator_wigner_setup proc~get_bb_r->proc~set_error_fatal proc~get_bb_r->proc~set_error_file proc~get_cc_r->interface~comms_bcast proc~get_cc_r->interface~comms_scatterv proc~get_cc_r->proc~comms_array_split proc~get_cc_r->proc~mpirank proc~get_cc_r->proc~mpisize proc~get_cc_r->interface~comms_reduce proc~get_cc_r->proc~fourier_loc_q_to_r proc~get_cc_r->proc~get_gauge_overlap_matrix proc~get_cc_r->proc~get_win_min proc~get_cc_r->proc~io_stopwatch_start proc~get_cc_r->proc~io_stopwatch_stop proc~get_cc_r->proc~operator_wigner_setup proc~get_cc_r->proc~set_error_fatal proc~get_cc_r->proc~set_error_file proc~utility_compar utility_compar proc~get_cc_r->proc~utility_compar proc~get_hh_r->interface~comms_bcast proc~get_hh_r->proc~mpirank 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_hh_r->proc~get_win_min proc~get_hh_r->proc~io_stopwatch_start proc~get_hh_r->proc~io_stopwatch_stop proc~get_hh_r->proc~operator_wigner_setup proc~get_hh_r->proc~set_error_fatal proc~get_hh_r->proc~set_error_file proc~get_shc_r->interface~comms_bcast proc~get_shc_r->proc~mpirank proc~get_shc_r->proc~fourier_q_to_r proc~get_shc_r->proc~get_gauge_overlap_matrix proc~get_shc_r->proc~io_stopwatch_start proc~get_shc_r->proc~io_stopwatch_stop proc~get_shc_r->proc~operator_wigner_setup proc~set_error_alloc set_error_alloc proc~get_shc_r->proc~set_error_alloc proc~set_error_dealloc set_error_dealloc proc~get_shc_r->proc~set_error_dealloc proc~get_shc_r->proc~set_error_fatal proc~get_shc_r->proc~set_error_file proc~get_ss_r->interface~comms_bcast proc~get_ss_r->proc~mpirank proc~get_ss_r->proc~fourier_q_to_r proc~get_ss_r->proc~get_gauge_overlap_matrix proc~get_ss_r->proc~io_stopwatch_start proc~get_ss_r->proc~io_stopwatch_stop proc~get_ss_r->proc~operator_wigner_setup proc~get_ss_r->proc~set_error_alloc proc~get_ss_r->proc~set_error_dealloc proc~get_ss_r->proc~set_error_fatal proc~get_ss_r->proc~set_error_file proc~utility_metric utility_metric proc~k_path_get_points->proc~utility_metric proc~k_path_print_info->proc~mpirank proc~k_path_print_info->proc~set_error_input proc~comms_sync_error comms_sync_error proc~set_error_input->proc~comms_sync_error proc~set_base_error set_base_error proc~set_error_input->proc~set_base_error proc~spin_get_nk->proc~pw90common_fourier_r_to_k proc~spin_get_nk->proc~utility_diagonalize proc~utility_rotate_diag utility_rotate_diag proc~spin_get_nk->proc~utility_rotate_diag 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~comms_reduce_cmplx comms_reduce_cmplx interface~comms_reduce->proc~comms_reduce_cmplx proc~comms_reduce_int comms_reduce_int interface~comms_reduce->proc~comms_reduce_int proc~comms_reduce_real comms_reduce_real interface~comms_reduce->proc~comms_reduce_real proc~kmesh_spacing_mesh kmesh_spacing_mesh interface~pw90common_kmesh_spacing->proc~kmesh_spacing_mesh proc~kmesh_spacing_singleinteger kmesh_spacing_singleinteger interface~pw90common_kmesh_spacing->proc~kmesh_spacing_singleinteger 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~get_gauge_overlap_matrix->proc~get_win_min proc~utility_zgemmm utility_zgemmm proc~get_gauge_overlap_matrix->proc~utility_zgemmm proc~set_error_alloc->proc~comms_sync_error proc~set_error_alloc->proc~set_base_error proc~set_error_dealloc->proc~comms_sync_error proc~set_error_dealloc->proc~set_base_error proc~set_error_fatal->proc~comms_sync_error proc~set_error_fatal->proc~set_base_error proc~set_error_file->proc~comms_sync_error proc~set_error_file->proc~set_base_error proc~utility_rotate_diag->proc~utility_zgemm_new proc~utility_matmul_diag utility_matmul_diag proc~utility_rotate_diag->proc~utility_matmul_diag zgemm zgemm proc~utility_zgemm_new->zgemm proc~wham_get_d_h->proc~utility_rotate 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~wham_get_eig_uu_hh_jjlist->proc~get_hh_r proc~wham_get_eig_uu_hh_jjlist->proc~utility_diagonalize proc~wham_get_eig_uu_hh_jjlist->proc~pw90common_fourier_r_to_k_new proc~wham_get_jjp_jjm_list wham_get_JJp_JJm_list proc~wham_get_eig_uu_hh_jjlist->proc~wham_get_jjp_jjm_list proc~wham_get_occ_mat_list->proc~set_error_input proc~wham_get_occ_mat_list->proc~pw90common_get_occ 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~comms_reduce_cmplx->proc~comms_sync_error proc~comms_no_sync_reduce_cmplx comms_no_sync_reduce_cmplx proc~comms_reduce_cmplx->proc~comms_no_sync_reduce_cmplx proc~comms_reduce_int->proc~comms_sync_error proc~comms_no_sync_reduce_int comms_no_sync_reduce_int proc~comms_reduce_int->proc~comms_no_sync_reduce_int proc~comms_reduce_real->proc~comms_sync_error proc~comms_no_sync_reduce_real comms_no_sync_reduce_real proc~comms_reduce_real->proc~comms_no_sync_reduce_real proc~utility_zgemmm->proc~utility_zgemm_new proc~wham_get_deleig_a->proc~utility_diagonalize proc~wham_get_deleig_a->proc~utility_rotate proc~wham_get_deleig_a->proc~utility_rotate_diag proc~utility_rotate_new utility_rotate_new proc~wham_get_jjp_jjm_list->proc~utility_rotate_new proc~utility_rotate_new->proc~utility_zgemm_new

Called by

proc~~k_path~~CalledByGraph proc~k_path k_path program~postw90 postw90 program~postw90->proc~k_path

Source Code

  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