k_slice Subroutine

public subroutine k_slice(pw90_berry, dis_manifold, fermi_energy_list, kmesh_info, kpt_latt, pw90_kslice, pw90_oper_read, pw90_band_deriv_degen, pw90_spin, ws_region, 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_bands, num_kpts, num_wann, num_valence_bands, effective_model, have_disentangled, seedname, stdout, timer, error, comm)

Uses

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

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
real(kind=dp), intent(in) :: kpt_latt(:,:)
type(pw90_kslice_mod_type), intent(in) :: pw90_kslice
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(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_bands
integer, intent(in) :: num_kpts
integer, intent(in) :: num_wann
integer, intent(in) :: num_valence_bands
logical, intent(in) :: effective_model
logical, intent(in) :: have_disentangled
character(len=50), intent(in) :: seedname
integer, intent(in) :: stdout
type(timer_list_type), intent(inout) :: timer
type(w90_error_type), intent(out), allocatable :: error
type(w90_comm_type), intent(in) :: comm

Calls

proc~~k_slice~~CallsGraph proc~k_slice k_slice interface~comms_gatherv comms_gatherv proc~k_slice->interface~comms_gatherv proc~berry_get_imf_klist berry_get_imf_klist proc~k_slice->proc~berry_get_imf_klist proc~berry_get_imfgh_klist berry_get_imfgh_klist proc~k_slice->proc~berry_get_imfgh_klist proc~berry_get_shc_klist berry_get_shc_klist proc~k_slice->proc~berry_get_shc_klist proc~comms_array_split comms_array_split proc~k_slice->proc~comms_array_split proc~get_aa_r get_AA_R proc~k_slice->proc~get_aa_r proc~get_aa_r_effective get_AA_R_effective proc~k_slice->proc~get_aa_r_effective proc~get_bb_r get_BB_R proc~k_slice->proc~get_bb_r proc~get_cc_r get_CC_R proc~k_slice->proc~get_cc_r proc~get_hh_r get_HH_R proc~k_slice->proc~get_hh_r proc~get_shc_r get_SHC_R proc~k_slice->proc~get_shc_r proc~get_ss_r get_SS_R proc~k_slice->proc~get_ss_r proc~kslice_print_info kslice_print_info proc~k_slice->proc~kslice_print_info proc~mpirank mpirank proc~k_slice->proc~mpirank proc~mpisize mpisize proc~k_slice->proc~mpisize proc~pw90common_fourier_r_to_k pw90common_fourier_R_to_k proc~k_slice->proc~pw90common_fourier_r_to_k proc~script_common script_common proc~k_slice->proc~script_common proc~script_fermi_lines script_fermi_lines proc~k_slice->proc~script_fermi_lines proc~set_error_fatal set_error_fatal proc~k_slice->proc~set_error_fatal proc~set_error_input set_error_input proc~k_slice->proc~set_error_input proc~spin_get_nk spin_get_nk proc~k_slice->proc~spin_get_nk proc~utility_diagonalize utility_diagonalize proc~k_slice->proc~utility_diagonalize proc~utility_recip_lattice utility_recip_lattice proc~k_slice->proc~utility_recip_lattice proc~utility_recip_lattice_base utility_recip_lattice_base proc~k_slice->proc~utility_recip_lattice_base proc~wham_get_eig_deleig wham_get_eig_deleig proc~k_slice->proc~wham_get_eig_deleig proc~write_coords_file write_coords_file proc~k_slice->proc~write_coords_file proc~write_data_file write_data_file proc~k_slice->proc~write_data_file 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~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 proc~berry_get_shc_klist->proc~wham_get_eig_deleig 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~comms_array_split->proc~mpisize proc~get_aa_r->proc~comms_array_split proc~get_aa_r->proc~mpirank proc~get_aa_r->proc~mpisize proc~get_aa_r->proc~set_error_fatal interface~comms_bcast comms_bcast proc~get_aa_r->interface~comms_bcast interface~comms_reduce comms_reduce proc~get_aa_r->interface~comms_reduce interface~comms_scatterv comms_scatterv proc~get_aa_r->interface~comms_scatterv 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_file set_error_file proc~get_aa_r->proc~set_error_file proc~get_aa_r_effective->proc~mpirank proc~get_aa_r_effective->proc~set_error_fatal proc~get_aa_r_effective->interface~comms_bcast 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_file proc~get_bb_r->proc~comms_array_split proc~get_bb_r->proc~mpirank proc~get_bb_r->proc~mpisize proc~get_bb_r->proc~set_error_fatal proc~get_bb_r->interface~comms_bcast proc~get_bb_r->interface~comms_reduce proc~get_bb_r->interface~comms_scatterv 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_file proc~get_cc_r->proc~comms_array_split proc~get_cc_r->proc~mpirank proc~get_cc_r->proc~mpisize proc~get_cc_r->proc~set_error_fatal proc~get_cc_r->interface~comms_bcast proc~get_cc_r->interface~comms_reduce proc~get_cc_r->interface~comms_scatterv 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_file proc~utility_compar utility_compar proc~get_cc_r->proc~utility_compar proc~get_hh_r->proc~mpirank proc~get_hh_r->proc~set_error_fatal proc~get_hh_r->proc~set_error_input proc~get_hh_r->interface~comms_bcast 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_file proc~get_shc_r->proc~mpirank proc~get_shc_r->proc~set_error_fatal proc~get_shc_r->interface~comms_bcast 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_file proc~get_ss_r->proc~mpirank proc~get_ss_r->proc~set_error_fatal proc~get_ss_r->interface~comms_bcast 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_file proc~kslice_print_info->proc~set_error_input proc~comms_sync_error comms_sync_error proc~set_error_fatal->proc~comms_sync_error proc~set_base_error set_base_error proc~set_error_fatal->proc~set_base_error proc~set_error_input->proc~comms_sync_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_recip_lattice->proc~set_error_fatal proc~utility_recip_lattice->proc~utility_recip_lattice_base proc~utility_inv3 utility_inv3 proc~utility_recip_lattice_base->proc~utility_inv3 proc~wham_get_eig_deleig->proc~get_hh_r proc~wham_get_eig_deleig->proc~pw90common_fourier_r_to_k proc~wham_get_eig_deleig->proc~utility_diagonalize proc~wham_get_deleig_a wham_get_deleig_a proc~wham_get_eig_deleig->proc~wham_get_deleig_a proc~comms_bcast_char 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_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~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~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_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~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_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_deleig_a->proc~utility_diagonalize proc~wham_get_deleig_a->proc~utility_rotate proc~wham_get_deleig_a->proc~utility_rotate_diag 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 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 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~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~utility_zgemmm->proc~utility_zgemm_new 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_slice~~CalledByGraph proc~k_slice k_slice program~postw90 postw90 program~postw90->proc~k_slice

Source Code

  subroutine k_slice(pw90_berry, dis_manifold, fermi_energy_list, kmesh_info, kpt_latt, &
                     pw90_kslice, pw90_oper_read, pw90_band_deriv_degen, pw90_spin, ws_region, &
                     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_bands, &
                     num_kpts, num_wann, num_valence_bands, effective_model, have_disentangled, &
                     seedname, stdout, timer, error, comm)
    !================================================!
    !
    !! Main routine
    !
    !================================================!

    use w90_postw90_types, only: pw90_kslice_mod_type, pw90_berry_mod_type, pw90_spin_mod_type, &
                                 pw90_band_deriv_degen_type, pw90_oper_read_type, pw90_spin_hall_type, wigner_seitz_type
    use w90_berry, only: berry_get_imf_klist, berry_get_imfgh_klist, berry_get_shc_klist
    use w90_comms, only: comms_bcast, w90_comm_type, mpirank, mpisize, comms_gatherv, comms_array_split
    use w90_constants, only: dp, twopi, 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_types, only: dis_manifold_type, kmesh_info_type, print_output_type, &
                         wannier_data_type, ws_region_type, ws_distance_type, timer_list_type
    use w90_postw90_common, only: pw90common_fourier_R_to_k
    use w90_spin, only: spin_get_nk
    use w90_utility, only: utility_diagonalize, utility_recip_lattice, utility_recip_lattice_base
    use w90_wan_ham, only: wham_get_eig_deleig

    implicit none

    ! arguments
    type(pw90_berry_mod_type), intent(in) :: pw90_berry
    type(dis_manifold_type), intent(in) :: dis_manifold
    real(kind=dp), allocatable, intent(in) :: fermi_energy_list(:)
    type(kmesh_info_type), intent(in) :: kmesh_info
    real(kind=dp), intent(in) :: kpt_latt(:, :)
    type(pw90_kslice_mod_type), intent(in) :: pw90_kslice
    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(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(:, :, :, :)
    complex(kind=dp), allocatable, intent(inout) :: BB_R(:, :, :, :)
    complex(kind=dp), allocatable, intent(inout) :: CC_R(:, :, :, :, :)
    complex(kind=dp), allocatable, intent(inout) :: HH_R(:, :, :)
    complex(kind=dp), allocatable, intent(inout) :: SH_R(:, :, :, :)
    complex(kind=dp), allocatable, intent(inout) :: SHR_R(:, :, :, :, :)
    complex(kind=dp), allocatable, intent(inout) :: SR_R(:, :, :, :, :)
    complex(kind=dp), allocatable, intent(inout) :: SS_R(:, :, :, :)
    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

    integer, intent(in) :: mp_grid(3)
    integer, intent(in) :: num_bands, num_kpts, num_wann, 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           :: iloc, itot, i1, i2, n, n1, n2, n3, i, nkpts, my_nkpts
    integer           :: scriptunit, dataunit, loop_kpt
    real(kind=dp)     :: avec_2d(3, 3), bvec(3, 3), yvec(3), zvec(3), &
                         b1mod, b2mod, ymod, cosb1b2, kcorner_cart(3), &
                         areab1b2, cosyb2, kpt(3), kpt_x, kpt_y, k1, k2, &
                         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), curv(3), morb(3), &
                         spn_k(num_wann), del_eig(num_wann, 3), Delta_k, Delta_E, &
                         zhat(3), vdum(3), rdum, shc_k_fermi(fermi_n)
    logical           :: plot_fermi_lines, plot_curv, plot_morb, &
                         fermi_lines_color, heatmap, plot_shc
    character(len=120) :: filename, square

    complex(kind=dp), allocatable :: HH(:, :)
    complex(kind=dp), allocatable :: delHH(:, :, :)
    complex(kind=dp), allocatable :: UU(:, :)
    real(kind=dp), allocatable :: eig(:)

    ! Output data buffers
    real(kind=dp), allocatable    :: coords(:, :), my_coords(:, :), &
                                     spndata(:, :), my_spndata(:, :), &
                                     bandsdata(:, :), my_bandsdata(:, :), &
                                     zdata(:, :), my_zdata(:, :)
    logical, allocatable          :: spnmask(:, :), my_spnmask(:, :)

    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.

    plot_fermi_lines = index(pw90_kslice%task, 'fermi_lines') > 0
    plot_curv = index(pw90_kslice%task, 'curv') > 0
    plot_morb = index(pw90_kslice%task, 'morb') > 0
    plot_shc = index(pw90_kslice%task, 'shc') > 0
    fermi_lines_color = pw90_kslice%fermi_lines_colour /= 'none'
    heatmap = plot_curv .or. plot_morb .or. plot_shc
    if (plot_fermi_lines .and. fermi_lines_color .and. heatmap) then
      call set_error_input(error, 'Error: spin-colored Fermi lines not allowed in ' &
                           //'curv/morb/shc heatmap plots', comm)
      return
    end if
    if (plot_shc) then
      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
      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

    if (on_root) then
      call kslice_print_info(plot_fermi_lines, fermi_lines_color, plot_curv, plot_morb, plot_shc, &
                             stdout, pw90_berry, fermi_energy_list, error, comm)
      if (allocated(error)) return
    end if

    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) 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 (fermi_lines_color) 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
    call utility_recip_lattice_base(real_lattice, recip_lattice, volume)
    ! Set Cartesian components of the vectors (b1,b2) spanning the slice

    bvec(1, :) = matmul(pw90_kslice%b1(:), recip_lattice(:, :))
    bvec(2, :) = matmul(pw90_kslice%b2(:), recip_lattice(:, :))
    ! z_vec (orthogonal to the slice)
    zvec(1) = bvec(1, 2)*bvec(2, 3) - bvec(1, 3)*bvec(2, 2)
    zvec(2) = bvec(1, 3)*bvec(2, 1) - bvec(1, 1)*bvec(2, 3)
    zvec(3) = bvec(1, 1)*bvec(2, 2) - bvec(1, 2)*bvec(2, 1)
    ! y_vec (orthogonal to b1=x_vec)
    yvec(1) = zvec(2)*bvec(1, 3) - zvec(3)*bvec(1, 2)
    yvec(2) = zvec(3)*bvec(1, 1) - zvec(1)*bvec(1, 3)
    yvec(3) = zvec(1)*bvec(1, 2) - zvec(2)*bvec(1, 1)
    ! Area (modulus b1 x b2 = z_vec)
    areab1b2 = sqrt(zvec(1)**2 + zvec(2)**2 + zvec(3)**2)
    if (areab1b2 < eps8) then
      call set_error_fatal(error, 'Error in kslice: Vectors pw90_kslice_b1 and pw90_kslice_b2 ' &
                           //'not linearly independent', comm)
      return
    end if
    ! This is the unit vector zvec/|zvec| which completes the triad
    ! in the 2D case
    bvec(3, :) = zvec(:)/areab1b2
    ! Now that we have bvec(3,:), we can compute the dual vectors
    ! avec_2d as in the 3D case
    call utility_recip_lattice(bvec, avec_2d, rdum, error, comm)
    if (allocated(error)) return

    ! Moduli b1,b2,y_vec
    b1mod = sqrt(bvec(1, 1)**2 + bvec(1, 2)**2 + bvec(1, 3)**2)
    b2mod = sqrt(bvec(2, 1)**2 + bvec(2, 2)**2 + bvec(2, 3)**2)
    ymod = sqrt(yvec(1)**2 + yvec(2)**2 + yvec(3)**2)
    ! Cosine of the angle between y_vec and b2
    cosyb2 = yvec(1)*bvec(2, 1) + yvec(2)*bvec(2, 2) + yvec(3)*bvec(2, 3)
    cosyb2 = cosyb2/(ymod*b2mod)
    ! Cosine of the angle between b1=x_vec and b2
    cosb1b2 = bvec(1, 1)*bvec(2, 1) + bvec(1, 2)*bvec(2, 2) + bvec(1, 3)*bvec(2, 3)
    cosb1b2 = cosb1b2/(b1mod*b2mod)
    if (abs(cosb1b2) < eps8 .and. abs(b1mod - b2mod) < eps8) then
      square = 'True'
    else
      square = 'False'
    end if

    nkpts = (pw90_kslice%kmesh2d(1) + 1)*(pw90_kslice%kmesh2d(2) + 1)

    ! Partition set of k-points into junks
    call comms_array_split(nkpts, counts, displs, comm)
    my_nkpts = counts(my_node_id)

    allocate (my_coords(2, my_nkpts))
    if (heatmap) allocate (my_zdata(3, my_nkpts))

    if (plot_fermi_lines) then
      allocate (HH(num_wann, num_wann))
      allocate (UU(num_wann, num_wann))
      allocate (eig(num_wann))
      if (fermi_lines_color) then
        allocate (delHH(num_wann, num_wann, 3))
        allocate (my_spndata(num_wann, my_nkpts))
        allocate (my_spnmask(num_wann, my_nkpts))
        my_spnmask = .false.
      else
        allocate (my_bandsdata(num_wann, my_nkpts))
      end if
    end if

    ! Loop over local portion of uniform mesh of k-points covering the slice,
    ! including all four borders

    do iloc = 1, my_nkpts
      itot = iloc - 1 + displs(my_node_id)
      i2 = itot/(pw90_kslice%kmesh2d(1) + 1) ! slow
      i1 = itot - i2*(pw90_kslice%kmesh2d(1) + 1) !fast
      ! k1 and k2 are the coefficients of the k-point in the basis
      ! (pw90_kslice_b1,pw90_kslice_b2)
      k1 = i1/real(pw90_kslice%kmesh2d(1), dp)
      k2 = i2/real(pw90_kslice%kmesh2d(2), dp)
      kpt = pw90_kslice%corner + k1*pw90_kslice%b1 + k2*pw90_kslice%b2
      ! Add to (k1,k2) the projection of pw90_kslice_corner on the
      ! (pw90_kslice_b1,pw90_kslice_b2) plane, expressed as a linear
      ! combination of pw90_kslice_b1 and pw90_kslice_b2
      kcorner_cart(:) = matmul(pw90_kslice%corner(:), recip_lattice(:, :))
      k1 = k1 + dot_product(kcorner_cart, avec_2d(1, :))/twopi
      k2 = k2 + dot_product(kcorner_cart, avec_2d(2, :))/twopi
      ! Convert to (kpt_x,kpt_y), the 2D Cartesian coordinates
      ! with x along x_vec=b1 and y along y_vec
      kpt_x = k1*b1mod + k2*b2mod*cosb1b2
      kpt_y = k2*b2mod*cosyb2

      my_coords(:, iloc) = [kpt_x, kpt_y]

      if (plot_fermi_lines) then
        if (fermi_lines_color) 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

          do n = 1, num_wann
            if (spn_k(n) > 1.0_dp - eps8) then
              spn_k(n) = 1.0_dp - eps8
            elseif (spn_k(n) < -1.0_dp + eps8) then
              spn_k(n) = -1.0_dp + eps8
            end if
          end do

          call wham_get_eig_deleig(dis_manifold, kpt_latt, pw90_band_deriv_degen, ws_region, &
                                   print_output, wannier_data, ws_distance, wigner_seitz, delHH, &
                                   HH, HH_R, u_matrix, UU, v_matrix, del_eig, eig, eigval, kpt, &
                                   real_lattice, scissors_shift, mp_grid, num_bands, num_kpts, &
                                   num_wann, num_valence_bands, effective_model, &
                                   have_disentangled, seedname, stdout, timer, error, comm)
          if (allocated(error)) return

          Delta_k = max(b1mod/pw90_kslice%kmesh2d(1), b2mod/pw90_kslice%kmesh2d(2))
        else
          call pw90common_fourier_R_to_k(ws_region, wannier_data, ws_distance, wigner_seitz, HH, &
                                         HH_R, kpt, real_lattice, mp_grid, 0, num_wann, error, comm)
          if (allocated(error)) return

          call utility_diagonalize(HH, num_wann, eig, UU, error, comm)
          if (allocated(error)) return

        end if

        if (allocated(my_bandsdata)) then
          my_bandsdata(:, iloc) = eig(:)
        else
          my_spndata(:, iloc) = spn_k(:)
          do n = 1, num_wann
            ! vdum = dE/dk projected on the k-slice
            zhat = zvec/sqrt(dot_product(zvec, zvec))
            vdum(:) = del_eig(n, :) - dot_product(del_eig(n, :), zhat)*zhat(:)
            Delta_E = sqrt(dot_product(vdum, vdum))*Delta_k
!                   Delta_E=Delta_E*sqrt(2.0_dp) ! optimize this factor
            my_spnmask(n, iloc) = abs(eig(n) - fermi_energy_list(1)) < Delta_E
          end do
        end if
      end if

      if (plot_curv) 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

        curv(1) = sum(imf_k_list(:, 1, 1))
        curv(2) = sum(imf_k_list(:, 2, 1))
        curv(3) = sum(imf_k_list(:, 3, 1))
        if (pw90_berry%curv_unit == 'bohr2') curv = curv/bohr**2
        ! Print _minus_ the Berry curvature
        my_zdata(:, iloc) = -curv(:)
      else 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
        morb(1) = sum(Morb_k(:, 1))
        morb(2) = sum(Morb_k(:, 2))
        morb(3) = sum(Morb_k(:, 3))
        my_zdata(:, iloc) = morb(:)
      else 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_zdata(1, iloc) = shc_k_fermi(1)
      end if

    end do !iloc

    ! Send results to root process
    if (on_root) then
      allocate (coords(2, nkpts))
    else
      allocate (coords(1, 1))
    end if
    call comms_gatherv(my_coords, 2*my_nkpts, coords, 2*counts, 2*displs, error, comm)
    if (allocated(error)) return

    if (allocated(my_spndata)) then
      if (on_root) then
        allocate (spndata(num_wann, nkpts))
      else
        allocate (spndata(1, 1))
      end if
      call comms_gatherv(my_spndata, num_wann*my_nkpts, &
                         spndata, num_wann*counts, num_wann*displs, error, comm)
      if (allocated(error)) return
    end if

    if (allocated(my_spnmask)) then
      if (on_root) then
        allocate (spnmask(num_wann, nkpts))
      else
        allocate (spnmask(1, 1))
      end if
      call comms_gatherv(my_spnmask(1, 1), num_wann*my_nkpts, &
                         spnmask(1, 1), num_wann*counts, num_wann*displs, error, comm)
      if (allocated(error)) return
    end if

    if (allocated(my_bandsdata)) then
      if (on_root) then
        allocate (bandsdata(num_wann, nkpts))
      else
        allocate (bandsdata(1, 1))
      end if
      call comms_gatherv(my_bandsdata, num_wann*my_nkpts, &
                         bandsdata, num_wann*counts, num_wann*displs, error, comm)
      if (allocated(error)) return
    end if

    ! This holds either -curv or morb
    if (allocated(my_zdata)) then
      if (on_root) then
        allocate (zdata(3, nkpts))
      else
        allocate (zdata(1, 1))
      end if
      call comms_gatherv(my_zdata, 3*my_nkpts, &
                         zdata, 3*counts, 3*displs, error, comm)
      if (allocated(error)) return
    end if

    ! Write output files
    if (on_root) then
      ! set kpt_x and kpt_y to last evaluated point
      kpt_x = coords(1, nkpts)
      kpt_y = coords(2, nkpts)

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

      if (.not. fermi_lines_color) then
        filename = trim(seedname)//'-kslice-coord.dat'
        call write_data_file(stdout, filename, '(2E16.8)', coords)
      end if

      if (allocated(bandsdata)) then
        ! For python
        filename = trim(seedname)//'-kslice-bands.dat'
        call write_data_file(stdout, filename, '(E16.8)', &
                             reshape(bandsdata, [1, nkpts*num_wann]))

        ! For gnuplot, using 'grid data' format
        if (.not. heatmap) then
          do n = 1, num_wann
            n1 = n/100
            n2 = (n - n1*100)/10
            n3 = n - n1*100 - n2*10
            filename = trim(seedname)//'-bnd_' &
                       //achar(48 + n1)//achar(48 + n2)//achar(48 + n3)//'.dat'

            call write_coords_file(stdout, filename, '(3E16.8)', coords, &
                                   reshape(bandsdata(n, :), [1, 1, nkpts]), &
                                   blocklen=pw90_kslice%kmesh2d(1) + 1)
          end do
        end if
      end if

      if (allocated(spndata)) then
        filename = trim(seedname)//'-kslice-fermi-spn.dat'
        call write_coords_file(stdout, filename, '(3E16.8)', coords, &
                               reshape(spndata, [1, num_wann, nkpts]), &
                               spnmask)
      end if

      if (allocated(my_zdata)) then
        if (plot_curv .or. plot_morb .or. plot_shc) then
          if (plot_morb) then ! ugly. But to keep the logic the same as other places
            filename = trim(seedname)//'-kslice-morb.dat'
          elseif (plot_curv) then
            filename = trim(seedname)//'-kslice-curv.dat'
          elseif (plot_shc) then
            filename = trim(seedname)//'-kslice-shc.dat'
          end if
          write (stdout, '(/,3x,a)') filename
          open (newunit=dataunit, file=filename, form='formatted')
          if (plot_shc) then
            if (pw90_berry%curv_unit == 'bohr2') zdata = zdata/bohr**2
            do loop_kpt = 1, nkpts
              write (dataunit, '(1E16.8)') zdata(1, loop_kpt)
            end do
          else
            do loop_kpt = 1, nkpts
              write (dataunit, '(4E16.8)') zdata(:, loop_kpt)
            end do
          end if
          write (dataunit, *) ' '
          close (dataunit)
        end if
      end if

      if (plot_fermi_lines .and. .not. fermi_lines_color .and. .not. heatmap) then
        !
        ! gnuplot script for black Fermi lines
        !
        filename = trim(seedname)//'-kslice-fermi_lines.gnu'
        write (stdout, '(/,3x,a)') filename
        open (newunit=scriptunit, file=filename, form='formatted')
        write (scriptunit, '(a)') "unset surface"
        write (scriptunit, '(a)') "set contour"
        write (scriptunit, '(a)') "set view map"
        write (scriptunit, '(a,f9.5)') "set cntrparam levels discrete ", &
          fermi_energy_list(1)
        write (scriptunit, '(a)') "set cntrparam bspline"
        do n = 1, num_wann
          n1 = n/100
          n2 = (n - n1*100)/10
          n3 = n - n1*100 - n2*10
          write (scriptunit, '(a)') "set table 'bnd_" &
            //achar(48 + n1)//achar(48 + n2)//achar(48 + n3)//".dat'"
          write (scriptunit, '(a)') "splot '"//trim(seedname)//"-bnd_" &
            //achar(48 + n1)//achar(48 + n2)//achar(48 + n3)//".dat'"
          write (scriptunit, '(a)') "unset table"
        end do
        write (scriptunit, '(a)') &
          "#Uncomment next two lines to create postscript"
        write (scriptunit, '(a)') "#set term post eps enh"
        write (scriptunit, '(a)') &
          "#set output '"//trim(seedname)//"-kslice-fermi_lines.eps'"
        write (scriptunit, '(a)') "set size ratio -1"
        write (scriptunit, '(a)') "unset tics"
        write (scriptunit, '(a)') "unset key"
        write (scriptunit, '(a)') &
          "#For postscript try changing lw 1 --> lw 2 in the next line"
        write (scriptunit, '(a)') "set style line 1 lt 1 lw 1"
        if (num_wann == 1) then
          write (scriptunit, '(a)') &
            "plot 'bnd_001.dat' using 1:2 w lines ls 1"
        else
          write (scriptunit, '(a)') &
            "plot 'bnd_001.dat' using 1:2 w lines ls 1,"//achar(92)
        end if
        do n = 2, num_wann - 1
          n1 = n/100
          n2 = (n - n1*100)/10
          n3 = n - n1*100 - n2*10
          write (scriptunit, '(a)') "     'bnd_" &
            //achar(48 + n1)//achar(48 + n2)//achar(48 + n3) &
            //".dat' using 1:2 w lines ls 1,"//achar(92)
        end do
        n = num_wann
        n1 = n/100
        n2 = (n - n1*100)/10
        n3 = n - n1*100 - n2*10
        write (scriptunit, '(a)') "     'bnd_" &
          //achar(48 + n1)//achar(48 + n2)//achar(48 + n3) &
          //".dat' using 1:2 w lines ls 1"
        close (scriptunit)
        !
        ! Python script for black Fermi lines
        !
        filename = trim(seedname)//'-kslice-fermi_lines.py'
        write (stdout, '(/,3x,a)') filename
        open (newunit=scriptunit, file=filename, form='formatted')
        call script_common(scriptunit, areab1b2, square, seedname)
        call script_fermi_lines(scriptunit, seedname, fermi_energy_list)
        write (scriptunit, '(a)') " "
        write (scriptunit, '(a)') "# Remove the axes"
        write (scriptunit, '(a)') "ax = pl.gca()"
        write (scriptunit, '(a)') "ax.xaxis.set_visible(False)"
        write (scriptunit, '(a)') "ax.yaxis.set_visible(False)"
        write (scriptunit, '(a)') " "
        write (scriptunit, '(a)') "pl.axes().set_aspect('equal')"
        write (scriptunit, '(a)') " "
        write (scriptunit, '(a)') "outfile = '"//trim(seedname)// &
          "-fermi_lines.pdf'"
        write (scriptunit, '(a)') " "
        write (scriptunit, '(a)') " "
        write (scriptunit, '(a)') "pl.savefig(outfile,bbox_inches='tight')"
        write (scriptunit, '(a)') "pl.show()"
        close (scriptunit)
      end if !plot_fermi_lines .and. .not.fermi_lines_color .and. .not.heatmap

      if (plot_fermi_lines .and. fermi_lines_color .and. .not. heatmap) then
        !
        ! gnuplot script for spin-colored Fermi lines
        !
        filename = trim(seedname)//'-kslice-fermi_lines.gnu'
        write (stdout, '(/,3x,a)') filename
        open (newunit=scriptunit, file=filename, form='formatted')
        write (scriptunit, '(a)') "unset key"
        write (scriptunit, '(a)') "unset tics"
        write (scriptunit, '(a)') "set cbtics"
        write (scriptunit, '(a)') &
          "set palette defined (-1 'blue', 0 'green', 1 'red')"
        write (scriptunit, '(a)') "set pm3d map"
        write (scriptunit, '(a)') "set zrange [-1:1]"
        write (scriptunit, '(a)') "set size ratio -1"
        write (scriptunit, '(a)') &
          "#Uncomment next two lines to create postscript"
        write (scriptunit, '(a)') "#set term post eps enh"
        write (scriptunit, '(a)') "#set output '" &
          //trim(seedname)//"-kslice-fermi_lines.eps'"
        write (scriptunit, '(a)') "splot '" &
          //trim(seedname)//"-kslice-fermi-spn.dat' with dots palette"
        !
        ! python script for spin-colored Fermi lines
        !
        filename = trim(seedname)//'-kslice-fermi_lines.py'
        write (stdout, '(/,3x,a)') filename
        open (newunit=scriptunit, file=filename, form='formatted')
        write (scriptunit, '(a)') "import pylab as pl"
        write (scriptunit, '(a)') "import numpy as np"
        write (scriptunit, '(a)') "data = np.loadtxt('"//trim(seedname)// &
          "-kslice-fermi-spn.dat')"
        write (scriptunit, '(a)') "x=data[:,0]"
        write (scriptunit, '(a)') "y=data[:,1]"
        write (scriptunit, '(a)') "z=data[:,2]"
        write (scriptunit, '(a)') &
          "pl.scatter(x,y,c=z,marker='+',s=2,cmap=pl.cm.jet)"
        write (scriptunit, '(a,F12.6,a)') &
          "pl.plot([0,", kpt_x, "],[0,0],color='black',linestyle='-'," &
          //"linewidth=0.5)"
        write (scriptunit, '(a,F12.6,a,F12.6,a,F12.6,a)') &
          "pl.plot([", kpt_x, ",", kpt_x, "],[0,", kpt_y, "],color='black'," &
          //"linestyle='-',linewidth=0.5)"
        write (scriptunit, '(a,F12.6,a,F12.6,a,F12.6,a)') &
          "pl.plot([0,", kpt_x, "],[", kpt_y, ",", kpt_y, &
          "],color='black',linestyle='-',linewidth=0.5)"
        write (scriptunit, '(a,F12.6,a)') "pl.plot([0,0],[0,", kpt_y, &
          "],color='black',linestyle='-',linewidth=0.5)"
        write (scriptunit, '(a,F12.6,a)') "pl.xlim([0,", kpt_x, "])"
        write (scriptunit, '(a,F12.6,a)') "pl.ylim([0,", kpt_y, "])"
        write (scriptunit, '(a)') "cbar=pl.colorbar()"
        write (scriptunit, '(a)') "ax = pl.gca()"
        write (scriptunit, '(a)') "ax.xaxis.set_visible(False)"
        write (scriptunit, '(a)') "ax.yaxis.set_visible(False)"
        write (scriptunit, '(a)') "pl.savefig('"//trim(seedname)// &
          "-kslice-fermi_lines.pdf',bbox_inches='tight')"
        write (scriptunit, '(a)') "pl.show()"
        close (scriptunit)
      end if ! plot_fermi_lines .and. fermi_lines_color .and. .not.heatmap

      if (heatmap .and. (.not. plot_shc)) then
        !
        ! python script for curvature/Morb/SHC heatmaps [+ black Fermi lines]
        !
        do i = 1, 3

          if (plot_curv .and. .not. plot_fermi_lines) then
            filename = trim(seedname)//'-kslice-curv_'//achar(119 + i)//'.py'
            write (stdout, '(/,3x,a)') filename
            open (newunit=scriptunit, file=filename, form='formatted')
          elseif (plot_curv .and. plot_fermi_lines) then
            filename = trim(seedname)//'-kslice-curv_'//achar(119 + i)// &
                       '+fermi_lines.py'
            write (stdout, '(/,3x,a)') filename
            open (newunit=scriptunit, file=filename, form='formatted')
          elseif (plot_morb .and. .not. plot_fermi_lines) then
            filename = trim(seedname)//'-kslice-morb_'//achar(119 + i)//'.py'
            write (stdout, '(/,3x,a)') filename
            open (newunit=scriptunit, file=filename, form='formatted')
          elseif (plot_morb .and. plot_fermi_lines) then
            filename = trim(seedname)//'-kslice-morb_'//achar(119 + i)// &
                       '+fermi_lines.py'
            write (stdout, '(/,3x,a)') filename
            open (newunit=scriptunit, file=filename, form='formatted')
          end if
          call script_common(scriptunit, areab1b2, square, seedname)
          if (plot_fermi_lines) call script_fermi_lines(scriptunit, seedname, fermi_energy_list)

          if (plot_curv) then
            write (scriptunit, '(a)') " "
            write (scriptunit, '(a)') "outfile = '"//trim(seedname)// &
              "-kslice-curv_"//achar(119 + i)//".pdf'"
            write (scriptunit, '(a)') " "
            write (scriptunit, '(a)') &
              "val = np.loadtxt('"//trim(seedname)// &
              "-kslice-curv.dat', usecols=("//achar(47 + i)//",))"
            write (scriptunit, '(a)') " "
            write (scriptunit, '(a)') &
                 "val_log=np.array([np.log10(abs(elem))*np.sign(elem) &
                 &if abs(elem)>10 else elem/10.0 for elem in val])"
            write (scriptunit, '(a)') " "
            write (scriptunit, '(a)') "if square: "
            write (scriptunit, '(a)') "  Z=val_log.reshape(dimy,dimx)"
            write (scriptunit, '(a)') "  mn=int(np.floor(Z.min()))"
            write (scriptunit, '(a)') "  mx=int(np.ceil(Z.max()))"
            write (scriptunit, '(a)') "  ticks=range(mn,mx+1)"
            write (scriptunit, '(a)') "  pl.contourf(x_coord,y_coord,Z," &
              //"ticks,origin='lower')"
            write (scriptunit, '(a)') "  #pl.imshow(Z,origin='lower'," &
              //"extent=(min(x_coord),max(x_coord),min(y_coord)," &
              //"max(y_coord)))"
            write (scriptunit, '(a)') "else: "
            write (scriptunit, '(a)') "  valint = interpolate.griddata((points_x," &
              //"points_y), val_log, (grid_x,grid_y), method='nearest') # or 'cubic' or 'nearest'"
            write (scriptunit, '(a)') "  mn=int(np.floor(valint.min()))"
            write (scriptunit, '(a)') "  mx=int(np.ceil(valint.max()))"
            write (scriptunit, '(a)') "  ticks=range(mn,mx+1)"
            write (scriptunit, '(a)') "  pl.contourf(xint,yint,valint,ticks)"
            write (scriptunit, '(a)') "  #pl.imshow(valint,origin='lower'," &
              //"extent=(min(xint),max(xint),min(yint),max(yint)))"
            write (scriptunit, '(a)') " "
            write (scriptunit, '(a)') "ticklabels=[]"
            write (scriptunit, '(a)') "for n in ticks:"
            write (scriptunit, '(a)') " if n<0: "
            write (scriptunit, '(a)') &
              "  ticklabels.append('-$10^{%d}$' % abs(n))"
            write (scriptunit, '(a)') " elif n==0:"
            write (scriptunit, '(a)') "  ticklabels.append(' $%d$' %  n)"
            write (scriptunit, '(a)') " else:"
            write (scriptunit, '(a)') "  ticklabels.append(' $10^{%d}$' % n)"
            write (scriptunit, '(a)') " "
            write (scriptunit, '(a)') "cbar=pl.colorbar()"
            write (scriptunit, '(a)') "cbar.set_ticks(ticks)"
            write (scriptunit, '(a)') "cbar.set_ticklabels(ticklabels)"

          elseif (plot_morb) then

            write (scriptunit, '(a)') " "
            write (scriptunit, '(a)') "outfile = '"//trim(seedname)// &
              "-kslice-morb_"//achar(119 + i)//".pdf'"
            write (scriptunit, '(a)') " "
            write (scriptunit, '(a)') &
              "val = np.loadtxt('"//trim(seedname)// &
              "-kslice-morb.dat', usecols=("//achar(47 + i)//",))"
            write (scriptunit, '(a)') " "
            write (scriptunit, '(a)') "if square: "
            write (scriptunit, '(a)') "  Z=val.reshape(dimy,dimx)"
            write (scriptunit, '(a)') "  pl.imshow(Z,origin='lower'," &
              //"extent=(min(x_coord),max(x_coord),min(y_coord)," &
              //"max(y_coord)))"
            write (scriptunit, '(a)') "else: "
            write (scriptunit, '(a)') "  valint = interpolate.griddata((points_x," &
              //"points_y), val_log, (grid_x,grid_y), method='nearest') # or 'cubic' or 'nearest'"
            write (scriptunit, '(a)') "  pl.imshow(valint,origin='lower'," &
              //"extent=(min(xint),max(xint),min(yint),max(yint)))"
            write (scriptunit, '(a)') "cbar=pl.colorbar()"

          end if

          write (scriptunit, '(a)') " "
          write (scriptunit, '(a)') "ax = pl.gca()"
          write (scriptunit, '(a)') "ax.xaxis.set_visible(False)"
          write (scriptunit, '(a)') "ax.yaxis.set_visible(False)"
          write (scriptunit, '(a)') " "
          write (scriptunit, '(a)') "pl.savefig(outfile,bbox_inches='tight')"
          write (scriptunit, '(a)') "pl.show()"

          close (scriptunit)

        end do !i

      end if !heatmap

      if (heatmap .and. plot_shc) then
        if (.not. plot_fermi_lines) then
          filename = trim(seedname)//'-kslice-shc'//'.py'
          write (stdout, '(/,3x,a)') filename
          open (newunit=scriptunit, file=filename, form='formatted')
        elseif (plot_fermi_lines) then
          filename = trim(seedname)//'-kslice-shc'//'+fermi_lines.py'
          write (stdout, '(/,3x,a)') filename
          open (newunit=scriptunit, file=filename, form='formatted')
        end if
        write (scriptunit, '(a)') "# uncomment these two lines if you are " &
          //"running in non-GUI environment"
        write (scriptunit, '(a)') "#import matplotlib"
        write (scriptunit, '(a)') "#matplotlib.use('Agg')"
        write (scriptunit, '(a)') "import matplotlib.pyplot as plt"
        call script_common(scriptunit, areab1b2, square, seedname)
        if (plot_fermi_lines) call script_fermi_lines(scriptunit, seedname, fermi_energy_list)

        write (scriptunit, '(a)') " "
        write (scriptunit, '(a)') "def shiftedColorMap(cmap, start=0, " &
          //"midpoint=0.5, stop=1.0, name='shiftedcmap'):"
        write (scriptunit, '(a)') "  '''"
        write (scriptunit, '(a)') '  Function to offset the "center" ' &
          //'of a colormap. Useful for'
        write (scriptunit, '(a)') '  data with a negative min and ' &
          //'positive max and you want the'
        write (scriptunit, '(a)') "  middle of the colormap's dynamic " &
          //"range to be at zero."
        write (scriptunit, '(a)') '  '
        write (scriptunit, '(a)') '  Input'
        write (scriptunit, '(a)') '  -----'
        write (scriptunit, '(a)') '  cmap : The matplotlib colormap to ' &
          //'be altered'
        write (scriptunit, '(a)') "  start : Offset from lowest point in " &
          //"the colormap's range."
        write (scriptunit, '(a)') '    Defaults to 0.0 (no lower offset). ' &
          //'Should be between'
        write (scriptunit, '(a)') '    0.0 and `midpoint`.'
        write (scriptunit, '(a)') '  midpoint : The new center of the ' &
          //'colormap. Defaults to '
        write (scriptunit, '(a)') '    0.5 (no shift). Should be between ' &
          //'0.0 and 1.0. In'
        write (scriptunit, '(a)') '    general, this should be  1 - ' &
          //'vmax / (vmax + abs(vmin))'
        write (scriptunit, '(a)') '    For example if your data range from ' &
          //'-15.0 to +5.0 and'
        write (scriptunit, '(a)') '    you want the center of the colormap ' &
          //'at 0.0, `midpoint`'
        write (scriptunit, '(a)') '    should be set to  1 - 5/(5 + 15)) ' &
          //'or 0.75'
        write (scriptunit, '(a)') "  stop : Offset from highest point in " &
          //"the colormap's range."
        write (scriptunit, '(a)') '    Defaults to 1.0 (no upper offset). ' &
          //'Should be between'
        write (scriptunit, '(a)') '    `midpoint` and 1.0.'
        write (scriptunit, '(a)') "  '''"
        write (scriptunit, '(a)') "  cdict = {'red': [],'green': []," &
          //"'blue': [],'alpha': []}"
        write (scriptunit, '(a)') '  # regular index to compute the colors'
        write (scriptunit, '(a)') '  reg_index = np.linspace(start, stop, 257)'
        write (scriptunit, '(a)') '  # shifted index to match the data'
        write (scriptunit, '(a)') '  shift_index = np.hstack(['
        write (scriptunit, '(a)') '    np.linspace(0.0, midpoint, 128, ' &
          //'endpoint=False),'
        write (scriptunit, '(a)') '    np.linspace(midpoint, 1.0, 129, ' &
          //'endpoint=True)'
        write (scriptunit, '(a)') '  ])'
        write (scriptunit, '(a)') '  for ri, si in zip(reg_index, shift_index):'
        write (scriptunit, '(a)') '    r, g, b, a = cmap(ri)'
        write (scriptunit, '(a)') "    cdict['red'].append((si, r, r))"
        write (scriptunit, '(a)') "    cdict['green'].append((si, g, g))"
        write (scriptunit, '(a)') "    cdict['blue'].append((si, b, b))"
        write (scriptunit, '(a)') "    cdict['alpha'].append((si, a, a))"
        write (scriptunit, '(a)') '  newcmap = matplotlib.colors' &
          //'.LinearSegmentedColormap(name, cdict)'
        write (scriptunit, '(a)') '  plt.register_cmap(cmap=newcmap)'
        write (scriptunit, '(a)') '  return newcmap'
        write (scriptunit, '(a)') " "
        write (scriptunit, '(a)') "outfile = '"//trim(seedname)//"-kslice-shc.pdf'"
        write (scriptunit, '(a)') " "
        write (scriptunit, '(a)') "val = np.loadtxt('"//trim(seedname) &
          //"-kslice-shc.dat', usecols=(0,))"
        write (scriptunit, '(a)') " "
        write (scriptunit, '(a)') "val_log=np.array([np.log10(abs(elem))*np.sign(elem)" &
          //"if abs(elem)>10 else elem/10.0 for elem in val])"
        write (scriptunit, '(a)') "#val_log = val"
        write (scriptunit, '(a)') "valmax=max(val_log)"
        write (scriptunit, '(a)') "valmin=min(val_log)"
        write (scriptunit, '(a)') "#cmnew=shiftedColorMap(matplotlib.cm.bwr," &
          //"0,1-valmax/(valmax+abs(valmin)),1)"
        write (scriptunit, '(a)') " "
        write (scriptunit, '(a)') "if square: "
        write (scriptunit, '(a)') "  Z=val_log.reshape(dimy,dimx)"
        write (scriptunit, '(a)') "  mn=int(np.floor(Z.min()))"
        write (scriptunit, '(a)') "  mx=int(np.ceil(Z.max()))"
        write (scriptunit, '(a)') "  ticks=range(mn,mx+1)"
        write (scriptunit, '(a)') "  #pl.contourf(x_coord,y_coord,Z," &
          //"ticks,origin='lower')"
        write (scriptunit, '(a)') "  pl.imshow(Z,origin='lower'," &
          //"extent=(min(x_coord),max(x_coord),min(y_coord)," &
          //"max(y_coord)))#,cmap=cmnew)"
        write (scriptunit, '(a)') "else: "
        write (scriptunit, '(a)') "  grid_x, grid_y = np.meshgrid(xint,yint)"
        write (scriptunit, '(a)') "  valint = interpolate.griddata((points_x," &
          //"points_y), val_log, (grid_x,grid_y), method='nearest')"
        write (scriptunit, '(a)') "  mn=int(np.floor(valint.min()))"
        write (scriptunit, '(a)') "  mx=int(np.ceil(valint.max()))"
        write (scriptunit, '(a)') "  ticks=range(mn,mx+1)"
        write (scriptunit, '(a)') "  #pl.contourf(xint,yint,valint,ticks)"
        write (scriptunit, '(a)') "  pl.imshow(valint,origin='lower'," &
          //"extent=(min(xint),max(xint),min(yint),max(yint)))#,cmap=cmnew)"
        write (scriptunit, '(a)') " "
        write (scriptunit, '(a)') "ticklabels=[]"
        write (scriptunit, '(a)') "for n in ticks:"
        write (scriptunit, '(a)') " if n<0: "
        write (scriptunit, '(a)') "  ticklabels.append('-$10^{%d}$' % abs(n))"
        write (scriptunit, '(a)') " elif n==0:"
        write (scriptunit, '(a)') "  ticklabels.append(' $%d$' %  n)"
        write (scriptunit, '(a)') " else:"
        write (scriptunit, '(a)') "  ticklabels.append(' $10^{%d}$' % n)"
        write (scriptunit, '(a)') " "
        write (scriptunit, '(a)') "cbar=pl.colorbar()"
        write (scriptunit, '(a)') "#cbar.set_ticks(ticks)"
        write (scriptunit, '(a)') "#cbar.set_ticklabels(ticklabels)"
        write (scriptunit, '(a)') " "
        write (scriptunit, '(a)') "ax = pl.gca()"
        write (scriptunit, '(a)') "ax.xaxis.set_visible(False)"
        write (scriptunit, '(a)') "ax.yaxis.set_visible(False)"
        write (scriptunit, '(a)') " "
        write (scriptunit, '(a)') "pl.savefig(outfile,bbox_inches='tight')"
        write (scriptunit, '(a)') "pl.show()"
      end if

      write (stdout, *) ' '

    end if ! on_root

  end subroutine k_slice