berry_main Subroutine

public subroutine berry_main(pw90_berry, dis_manifold, fermi_energy_list, kmesh_info, kpoint_dist, kpt_latt, pw90_band_deriv_degen, pw90_oper_read, pw90_spin, physics, ws_region, pw90_spin_hall, wannier_data, ws_distance, wigner_seitz, print_output, AA_R, BB_R, CC_R, HH_R, SH_R, SHR_R, SR_R, SS_R, SAA_R, SBB_R, u_matrix, v_matrix, eigval, real_lattice, scissors_shift, mp_grid, fermi_n, num_wann, num_kpts, num_bands, num_valence_bands, effective_model, have_disentangled, spin_decomp, seedname, stdout, timer, error, comm)

Uses

  • proc~~berry_main~~UsesGraph proc~berry_main berry_main module~w90_comms w90_comms proc~berry_main->module~w90_comms module~w90_constants w90_constants proc~berry_main->module~w90_constants module~w90_get_oper w90_get_oper proc~berry_main->module~w90_get_oper module~w90_io w90_io proc~berry_main->module~w90_io module~w90_postw90_types w90_postw90_types proc~berry_main->module~w90_postw90_types module~w90_tetrahedron w90_tetrahedron proc~berry_main->module~w90_tetrahedron module~w90_types w90_types proc~berry_main->module~w90_types module~w90_utility w90_utility proc~berry_main->module~w90_utility module~w90_comms->module~w90_constants module~w90_error_base w90_error_base module~w90_comms->module~w90_error_base module~w90_get_oper->module~w90_comms module~w90_get_oper->module~w90_constants module~w90_get_oper->module~w90_io module~w90_error w90_error module~w90_get_oper->module~w90_error module~w90_io->module~w90_constants module~w90_postw90_types->module~w90_comms module~w90_postw90_types->module~w90_constants module~w90_tetrahedron->module~w90_constants 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

Computes the following quantities: (i) Anomalous Hall conductivity (from Berry curvature) (ii) Complex optical conductivity (Kubo-Greenwood) & JDOS (iii) Orbital magnetization (iv) Nonlinear shift current (v) Spin Hall conductivity

Arguments

Type IntentOptional Attributes Name
type(pw90_berry_mod_type), intent(inout) :: 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(kpoint_dist_type), intent(in) :: kpoint_dist
real(kind=dp), intent(in) :: kpt_latt(:,:)
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(pw90_physical_constants_type), intent(in) :: physics
type(ws_region_type), intent(in) :: ws_region
type(pw90_spin_hall_type), intent(in) :: pw90_spin_hall
type(wannier_data_type), intent(in) :: wannier_data
type(ws_distance_type), intent(inout) :: ws_distance
type(wigner_seitz_type), intent(inout) :: wigner_seitz
type(print_output_type), intent(in) :: print_output
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) :: u_matrix(:,:,:)
complex(kind=dp), intent(in) :: v_matrix(:,:,:)
real(kind=dp), intent(in) :: eigval(:,:)
real(kind=dp), intent(in) :: real_lattice(3,3)
real(kind=dp), intent(in) :: scissors_shift
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: fermi_n
integer, intent(in) :: num_wann
integer, intent(in) :: num_kpts
integer, intent(in) :: num_bands
integer, intent(in) :: num_valence_bands
logical, intent(in) :: effective_model
logical, intent(in) :: have_disentangled
logical, intent(in) :: spin_decomp
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~~berry_main~~CallsGraph proc~berry_main berry_main interface~comms_reduce comms_reduce proc~berry_main->interface~comms_reduce proc~berry_get_imf_klist berry_get_imf_klist proc~berry_main->proc~berry_get_imf_klist proc~berry_get_imfgh_klist berry_get_imfgh_klist proc~berry_main->proc~berry_get_imfgh_klist proc~berry_get_kdotp berry_get_kdotp proc~berry_main->proc~berry_get_kdotp proc~berry_get_kubo_k berry_get_kubo_k proc~berry_main->proc~berry_get_kubo_k proc~berry_get_sc_klist berry_get_sc_klist proc~berry_main->proc~berry_get_sc_klist proc~berry_get_shc_klist berry_get_shc_klist proc~berry_main->proc~berry_get_shc_klist proc~berry_get_shc_tetrahedron berry_get_shc_tetrahedron proc~berry_main->proc~berry_get_shc_tetrahedron proc~berry_print_progress berry_print_progress proc~berry_main->proc~berry_print_progress proc~comms_array_split comms_array_split proc~berry_main->proc~comms_array_split proc~get_aa_r get_AA_R proc~berry_main->proc~get_aa_r proc~get_aa_r_effective get_AA_R_effective proc~berry_main->proc~get_aa_r_effective proc~get_bb_r get_BB_R proc~berry_main->proc~get_bb_r proc~get_cc_r get_CC_R proc~berry_main->proc~get_cc_r proc~get_hh_r get_HH_R proc~berry_main->proc~get_hh_r proc~get_saa_r get_SAA_R proc~berry_main->proc~get_saa_r proc~get_sbb_r get_SBB_R proc~berry_main->proc~get_sbb_r proc~get_sh_r get_SH_R proc~berry_main->proc~get_sh_r proc~get_shc_r get_SHC_R proc~berry_main->proc~get_shc_r proc~get_ss_r get_SS_R proc~berry_main->proc~get_ss_r proc~io_stopwatch_start io_stopwatch_start proc~berry_main->proc~io_stopwatch_start proc~io_stopwatch_stop io_stopwatch_stop proc~berry_main->proc~io_stopwatch_stop proc~mpirank mpirank proc~berry_main->proc~mpirank proc~mpisize mpisize proc~berry_main->proc~mpisize proc~set_error_alloc set_error_alloc proc~berry_main->proc~set_error_alloc proc~set_error_input set_error_input proc~berry_main->proc~set_error_input proc~tetrahedron_array_init tetrahedron_array_init proc~berry_main->proc~tetrahedron_array_init proc~tetrahedron_p_matrix_init tetrahedron_P_matrix_init proc~berry_main->proc~tetrahedron_p_matrix_init proc~tetrahedron_spinhall tetrahedron_spinhall proc~berry_main->proc~tetrahedron_spinhall 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~berry_get_imf_klist->proc~berry_get_imfgh_klist proc~pw90common_fourier_r_to_k pw90common_fourier_R_to_k 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~utility_rotate utility_rotate proc~berry_get_kdotp->proc~utility_rotate proc~wham_get_d_h_p_value wham_get_D_h_P_value proc~berry_get_kdotp->proc~wham_get_d_h_p_value proc~wham_get_eig_deleig wham_get_eig_deleig proc~berry_get_kdotp->proc~wham_get_eig_deleig proc~wham_get_eig_uu_hh_aa_sc wham_get_eig_UU_HH_AA_sc proc~berry_get_kdotp->proc~wham_get_eig_uu_hh_aa_sc interface~pw90common_kmesh_spacing pw90common_kmesh_spacing proc~berry_get_kubo_k->interface~pw90common_kmesh_spacing proc~pw90common_fourier_r_to_k_new pw90common_fourier_R_to_k_new proc~berry_get_kubo_k->proc~pw90common_fourier_r_to_k_new proc~berry_get_kubo_k->proc~pw90common_fourier_r_to_k_vec proc~pw90common_get_occ pw90common_get_occ proc~berry_get_kubo_k->proc~pw90common_get_occ proc~spin_get_nk spin_get_nk proc~berry_get_kubo_k->proc~spin_get_nk proc~utility_diagonalize utility_diagonalize proc~berry_get_kubo_k->proc~utility_diagonalize proc~utility_recip_lattice_base utility_recip_lattice_base proc~berry_get_kubo_k->proc~utility_recip_lattice_base proc~berry_get_kubo_k->proc~utility_rotate proc~utility_w0gauss utility_w0gauss proc~berry_get_kubo_k->proc~utility_w0gauss proc~wham_get_d_h wham_get_D_h proc~berry_get_kubo_k->proc~wham_get_d_h proc~berry_get_kubo_k->proc~wham_get_eig_deleig dger dger proc~berry_get_sc_klist->dger proc~berry_get_sc_klist->interface~pw90common_kmesh_spacing proc~pw90common_fourier_r_to_k_vec_dadb pw90common_fourier_R_to_k_vec_dadb proc~berry_get_sc_klist->proc~pw90common_fourier_r_to_k_vec_dadb proc~pw90common_fourier_r_to_k_vec_dadb_tb_conv pw90common_fourier_R_to_k_vec_dadb_TB_conv proc~berry_get_sc_klist->proc~pw90common_fourier_r_to_k_vec_dadb_tb_conv proc~berry_get_sc_klist->proc~pw90common_get_occ proc~berry_get_sc_klist->proc~utility_recip_lattice_base proc~berry_get_sc_klist->proc~utility_rotate proc~utility_w0gauss_vec utility_w0gauss_vec proc~berry_get_sc_klist->proc~utility_w0gauss_vec proc~utility_zdotu utility_zdotu proc~berry_get_sc_klist->proc~utility_zdotu proc~berry_get_sc_klist->proc~wham_get_d_h proc~berry_get_sc_klist->proc~wham_get_d_h_p_value proc~berry_get_sc_klist->proc~wham_get_eig_deleig proc~wham_get_eig_deleig_tb_conv wham_get_eig_deleig_TB_conv proc~berry_get_sc_klist->proc~wham_get_eig_deleig_tb_conv proc~berry_get_sc_klist->proc~wham_get_eig_uu_hh_aa_sc proc~wham_get_eig_uu_hh_aa_sc_tb_conv wham_get_eig_UU_HH_AA_sc_TB_conv proc~berry_get_sc_klist->proc~wham_get_eig_uu_hh_aa_sc_tb_conv proc~berry_get_shc_klist->interface~pw90common_kmesh_spacing 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~berry_get_shc_klist->proc~pw90common_get_occ proc~berry_get_shc_klist->proc~utility_recip_lattice_base proc~berry_get_shc_klist->proc~utility_rotate proc~berry_get_shc_klist->proc~wham_get_d_h proc~berry_get_shc_klist->proc~wham_get_eig_deleig proc~berry_get_shc_tetrahedron->proc~pw90common_fourier_r_to_k_new proc~berry_get_shc_tetrahedron->proc~pw90common_fourier_r_to_k_vec proc~berry_get_shc_tetrahedron->proc~utility_diagonalize proc~berry_get_shc_tetrahedron->proc~utility_rotate proc~io_wallclocktime io_wallclocktime proc~berry_print_progress->proc~io_wallclocktime proc~comms_array_split->proc~mpisize proc~get_aa_r->interface~comms_reduce proc~get_aa_r->proc~comms_array_split proc~get_aa_r->proc~io_stopwatch_start proc~get_aa_r->proc~io_stopwatch_stop proc~get_aa_r->proc~mpirank proc~get_aa_r->proc~mpisize interface~comms_bcast comms_bcast proc~get_aa_r->interface~comms_bcast 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~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->proc~io_stopwatch_start proc~get_aa_r_effective->proc~io_stopwatch_stop proc~get_aa_r_effective->proc~mpirank proc~get_aa_r_effective->interface~comms_bcast proc~get_aa_r_effective->proc~set_error_fatal proc~get_aa_r_effective->proc~set_error_file proc~get_bb_r->interface~comms_reduce proc~get_bb_r->proc~comms_array_split proc~get_bb_r->proc~io_stopwatch_start proc~get_bb_r->proc~io_stopwatch_stop proc~get_bb_r->proc~mpirank proc~get_bb_r->proc~mpisize proc~get_bb_r->interface~comms_bcast 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~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_reduce proc~get_cc_r->proc~comms_array_split proc~get_cc_r->proc~io_stopwatch_start proc~get_cc_r->proc~io_stopwatch_stop proc~get_cc_r->proc~mpirank proc~get_cc_r->proc~mpisize proc~get_cc_r->interface~comms_bcast 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~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->proc~io_stopwatch_start proc~get_hh_r->proc~io_stopwatch_stop proc~get_hh_r->proc~mpirank proc~get_hh_r->proc~set_error_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~operator_wigner_setup proc~get_hh_r->proc~set_error_fatal proc~get_hh_r->proc~set_error_file proc~get_saa_r->interface~comms_reduce proc~get_saa_r->proc~comms_array_split proc~get_saa_r->proc~io_stopwatch_start proc~get_saa_r->proc~io_stopwatch_stop proc~get_saa_r->proc~mpirank proc~get_saa_r->proc~mpisize proc~get_saa_r->interface~comms_bcast proc~get_saa_r->interface~comms_scatterv proc~get_saa_r->proc~fourier_loc_q_to_r proc~get_saa_r->proc~get_win_min proc~get_saa_r->proc~operator_wigner_setup proc~get_saa_r->proc~set_error_fatal proc~get_saa_r->proc~set_error_file proc~get_sbb_r->interface~comms_reduce proc~get_sbb_r->proc~comms_array_split proc~get_sbb_r->proc~io_stopwatch_start proc~get_sbb_r->proc~io_stopwatch_stop proc~get_sbb_r->proc~mpirank proc~get_sbb_r->proc~mpisize proc~get_sbb_r->interface~comms_bcast proc~get_sbb_r->interface~comms_scatterv proc~get_sbb_r->proc~fourier_loc_q_to_r proc~get_sbb_r->proc~get_win_min proc~get_sbb_r->proc~operator_wigner_setup proc~get_sbb_r->proc~set_error_fatal proc~get_sbb_r->proc~set_error_file proc~get_sh_r->proc~io_stopwatch_start proc~get_sh_r->proc~io_stopwatch_stop proc~get_sh_r->proc~mpirank proc~get_sh_r->proc~set_error_alloc proc~get_sh_r->interface~comms_bcast proc~get_sh_r->proc~fourier_q_to_r proc~get_sh_r->proc~get_gauge_overlap_matrix proc~get_sh_r->proc~operator_wigner_setup proc~set_error_dealloc set_error_dealloc proc~get_sh_r->proc~set_error_dealloc proc~get_sh_r->proc~set_error_fatal proc~get_sh_r->proc~set_error_file proc~get_shc_r->proc~io_stopwatch_start proc~get_shc_r->proc~io_stopwatch_stop proc~get_shc_r->proc~mpirank proc~get_shc_r->proc~set_error_alloc 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~operator_wigner_setup 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->proc~io_stopwatch_start proc~get_ss_r->proc~io_stopwatch_stop proc~get_ss_r->proc~mpirank proc~get_ss_r->proc~set_error_alloc 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~operator_wigner_setup 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~comms_sync_error comms_sync_error proc~set_error_alloc->proc~comms_sync_error proc~set_base_error set_base_error proc~set_error_alloc->proc~set_base_error proc~set_error_input->proc~comms_sync_error proc~set_error_input->proc~set_base_error proc~tetrahedron_fermidirac tetrahedron_fermidirac proc~tetrahedron_spinhall->proc~tetrahedron_fermidirac 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_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_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~get_gauge_overlap_matrix->proc~get_win_min proc~utility_zgemmm utility_zgemmm proc~get_gauge_overlap_matrix->proc~utility_zgemmm proc~utility_cart_to_frac utility_cart_to_frac proc~pw90common_fourier_r_to_k_vec_dadb_tb_conv->proc~utility_cart_to_frac proc~utility_inverse_mat utility_inverse_mat proc~pw90common_fourier_r_to_k_vec_dadb_tb_conv->proc~utility_inverse_mat 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~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~tetrahedron_integral tetrahedron_integral proc~tetrahedron_fermidirac->proc~tetrahedron_integral proc~tetrahedron_sort tetrahedron_sort proc~tetrahedron_fermidirac->proc~tetrahedron_sort 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~utility_w0gauss->proc~set_error_input proc~utility_w0gauss_vec->proc~set_error_alloc proc~utility_w0gauss_vec->proc~set_error_input zgemm zgemm proc~utility_zgemm_new->zgemm proc~wham_get_d_h->proc~utility_rotate proc~wham_get_d_h_p_value->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_deleig_tb_conv->proc~wham_get_deleig_a proc~wham_get_eig_uu_hh_aa_sc->proc~get_hh_r proc~wham_get_eig_uu_hh_aa_sc->proc~utility_diagonalize proc~pw90common_fourier_r_to_k_new_second_d pw90common_fourier_R_to_k_new_second_d proc~wham_get_eig_uu_hh_aa_sc->proc~pw90common_fourier_r_to_k_new_second_d proc~wham_get_eig_uu_hh_aa_sc_tb_conv->proc~get_aa_r proc~wham_get_eig_uu_hh_aa_sc_tb_conv->proc~get_aa_r_effective proc~wham_get_eig_uu_hh_aa_sc_tb_conv->proc~get_hh_r proc~wham_get_eig_uu_hh_aa_sc_tb_conv->proc~utility_diagonalize proc~pw90common_fourier_r_to_k_new_second_d_tb_conv pw90common_fourier_R_to_k_new_second_d_TB_conv proc~wham_get_eig_uu_hh_aa_sc_tb_conv->proc~pw90common_fourier_r_to_k_new_second_d_tb_conv proc~wham_get_eig_uu_hh_jjlist->proc~get_hh_r proc~wham_get_eig_uu_hh_jjlist->proc~pw90common_fourier_r_to_k_new proc~wham_get_eig_uu_hh_jjlist->proc~utility_diagonalize 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 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~pw90common_fourier_r_to_k_new_second_d_tb_conv->proc~utility_cart_to_frac proc~pw90common_fourier_r_to_k_new_second_d_tb_conv->proc~utility_inverse_mat proc~tetrahedron_integral->proc~tetrahedron_sort proc~tetrahedron_integral->proc~utility_inv3 proc~tetrahedron_jacobian tetrahedron_jacobian proc~tetrahedron_integral->proc~tetrahedron_jacobian proc~tetrahedron_log1p tetrahedron_log1p proc~tetrahedron_integral->proc~tetrahedron_log1p proc~utility_inverse_mat->proc~utility_inv3 proc~utility_rotate_diag->proc~utility_zgemm_new proc~utility_matmul_diag utility_matmul_diag proc~utility_rotate_diag->proc~utility_matmul_diag 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~~berry_main~~CalledByGraph proc~berry_main berry_main program~postw90 postw90 program~postw90->proc~berry_main

Source Code

  subroutine berry_main(pw90_berry, dis_manifold, fermi_energy_list, kmesh_info, kpoint_dist, &
                        kpt_latt, pw90_band_deriv_degen, pw90_oper_read, pw90_spin, physics, &
                        ws_region, pw90_spin_hall, wannier_data, ws_distance, wigner_seitz, &
                        print_output, AA_R, BB_R, CC_R, HH_R, SH_R, SHR_R, SR_R, SS_R, SAA_R, &
                        SBB_R, u_matrix, v_matrix, eigval, real_lattice, scissors_shift, mp_grid, &
                        fermi_n, num_wann, num_kpts, num_bands, num_valence_bands, &
                        effective_model, have_disentangled, spin_decomp, seedname, stdout, timer, &
                        error, comm)
    !================================================!
    !
    !! Computes the following quantities:
    !!   (i) Anomalous Hall conductivity (from Berry curvature)
    !!  (ii) Complex optical conductivity (Kubo-Greenwood) & JDOS
    !! (iii) Orbital magnetization
    !!  (iv) Nonlinear shift current
    !!   (v) Spin Hall conductivity
    !
    !================================================!

    use w90_comms, only: comms_reduce, w90_comm_type, mpirank, mpisize, comms_array_split
    use w90_constants, only: dp, cmplx_0, cmplx_i, pi, pw90_physical_constants_type
    use w90_utility, only: utility_recip_lattice_base
    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, &
                            get_SH_R, get_SAA_R, get_SBB_R
    use w90_io, only: io_stopwatch_start, io_stopwatch_stop
    use w90_types, only: 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_mod_type, &
                                 pw90_spin_hall_type, pw90_band_deriv_degen_type, pw90_oper_read_type, wigner_seitz_type, &
                                 kpoint_dist_type
    use w90_tetrahedron

    implicit none

    ! arguments
    type(pw90_berry_mod_type), intent(inout) :: pw90_berry
    type(dis_manifold_type), intent(in) :: dis_manifold
    type(kmesh_info_type), intent(in) :: kmesh_info
    type(kpoint_dist_type), intent(in) :: kpoint_dist
    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(pw90_physical_constants_type), intent(in) :: physics
    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(:, :, :, :) ! <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>
    !spin Hall using Ryoo's method
    complex(kind=dp), allocatable, intent(inout) :: SAA_R(:, :, :, :, :) ! <0n|sigma_x,y,z.(r-R)_alpha|Rm>
    !! $$\langle 0n | \sigma_{x,y,z}.(\hat{r}-R)_{\alpha}  | Rm \rangle$$
    complex(kind=dp), allocatable, intent(inout) :: SBB_R(:, :, :, :, :) ! <0n|sigma_x,y,z.H.(r-R)_alpha|Rm>
    !! $$\langle 0n | \sigma_{x,y,z}.H.(\hat{r}-R)_{\alpha}  | Rm \rangle$$
    complex(kind=dp), intent(in) :: u_matrix(:, :, :), v_matrix(:, :, :)

    real(kind=dp), intent(in) :: eigval(:, :)
    real(kind=dp), intent(in) :: real_lattice(3, 3)
    real(kind=dp), intent(in) :: scissors_shift
    real(kind=dp), 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_kpts, num_bands, num_valence_bands, fermi_n
    integer, intent(in) :: stdout

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

    ! local variables
    real(kind=dp), allocatable :: adkpt(:, :)

    ! AHC and orbital magnetization, calculated for a list of Fermi levels
    !
    ! First index labels J0,J1,J2 terms, second labels the Cartesian component
    !
    real(kind=dp) :: imf_k_list(3, 3, fermi_n), imf_list(3, 3, fermi_n), imf_list2(3, 3, fermi_n)
    real(kind=dp) :: img_k_list(3, 3, fermi_n), img_list(3, 3, fermi_n)
    real(kind=dp) :: imh_k_list(3, 3, fermi_n), imh_list(3, 3, fermi_n)
    real(kind=dp) :: ahc_list(3, 3, fermi_n)
    real(kind=dp) :: LCtil_list(3, 3, fermi_n), ICtil_list(3, 3, fermi_n), Morb_list(3, 3, fermi_n)
    real(kind=dp) :: imf_k_list_dummy(3, 3, fermi_n) ! adaptive refinement of AHC
    ! shift current
    real(kind=dp), allocatable :: sc_k_list(:, :, :)
    real(kind=dp), allocatable :: sc_list(:, :, :)
    ! kdotp
    complex(kind=dp), allocatable :: kdotp(:, :, :, :, :)
    ! Complex optical conductivity, dividided into Hermitean and
    ! anti-Hermitean parts
    !
    complex(kind=dp), allocatable :: kubo_H_k(:, :, :)
    complex(kind=dp), allocatable :: kubo_H(:, :, :)
    complex(kind=dp), allocatable :: kubo_AH_k(:, :, :)
    complex(kind=dp), allocatable :: kubo_AH(:, :, :)
    ! decomposition into up-up, down-down and spin-flip transitions
    complex(kind=dp), allocatable :: kubo_H_k_spn(:, :, :, :)
    complex(kind=dp), allocatable :: kubo_H_spn(:, :, :, :)
    complex(kind=dp), allocatable :: kubo_AH_k_spn(:, :, :, :)
    complex(kind=dp), allocatable :: kubo_AH_spn(:, :, :, :)

    ! Joint density of states
    !
    real(kind=dp), allocatable :: jdos_k(:)
    real(kind=dp), allocatable :: jdos(:)
    ! decomposition into up-up, down-down and spin-flip transitions
    real(kind=dp), allocatable :: jdos_k_spn(:, :)
    real(kind=dp), allocatable :: jdos_spn(:, :)

    ! Spin Hall conductivity
    real(kind=dp), allocatable :: shc_fermi(:), shc_k_fermi(:)
    complex(kind=dp), allocatable :: shc_freq(:), shc_k_freq(:)
    ! Tetrahedron method
    real(kind=dp), allocatable :: imjv(:, :, :, :, :)
    real(kind=dp), allocatable :: eig(:, :, :, :)
    real(kind=dp), allocatable :: imjv_tet(:, :, :)
    real(kind=dp), allocatable :: eig_tet(:, :)
    integer, allocatable :: counts(:)
    integer, allocatable :: displs(:)
    real(kind=dp)     :: kptc(3, 64), kptv(4, 3), &
                         Ftet(4), E1tet(4), E2tet(4), ttet(3, 3), omega, Ef
    complex(kind=dp)  :: shc_k_tet
    integer           :: itet, m, l, nfreq, tet_array(6, 20)
    real(kind=dp)     :: E1_opt(20), E2_opt(20), F_opt(20), P_matrix(4, 20)
    real(kind=dp), parameter :: mesh_shift = 0.5_dp

    ! for fermi energy scan, adaptive kmesh
    real(kind=dp), allocatable :: shc_k_fermi_dummy(:)

    real(kind=dp) :: cell_volume
    real(kind=dp) :: kweight, kweight_adpt, kpt(3), db1, db2, db3, fac, rdum, vdum(3)

    integer :: n, i, j, k, jk, ikpt, if, ierr, loop_x, loop_y, loop_z, kdotp_nbands
    integer :: loop_xyz, loop_adpt, adpt_counter_list(fermi_n), ifreq, file_unit
    integer :: my_node_id, num_nodes

    character(len=120) :: file_name

    logical :: eval_ahc, eval_morb, eval_kubo, not_scannable, eval_sc, eval_shc, eval_kdotp
    logical :: ladpt_kmesh
    logical :: ladpt(fermi_n)

    my_node_id = mpirank(comm)
    num_nodes = mpisize(comm)

    if (fermi_n == 0) then
      call set_error_input(error, 'Must specify one or more Fermi levels when berry=true', comm)
      return
    end if

    if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
      call io_stopwatch_start('berry: prelims', timer)

    cell_volume = real_lattice(1, 1)*(real_lattice(2, 2)*real_lattice(3, 3) - &
                                      real_lattice(3, 2)*real_lattice(2, 3)) + &
                  real_lattice(1, 2)*(real_lattice(2, 3)*real_lattice(3, 1) - &
                                      real_lattice(3, 3)*real_lattice(2, 1)) + &
                  real_lattice(1, 3)*(real_lattice(2, 1)*real_lattice(3, 2) - &
                                      real_lattice(3, 1)*real_lattice(2, 2))

    ! Mesh spacing in reduced coordinates
    !
    db1 = 1.0_dp/real(pw90_berry%kmesh%mesh(1), dp)
    db2 = 1.0_dp/real(pw90_berry%kmesh%mesh(2), dp)
    db3 = 1.0_dp/real(pw90_berry%kmesh%mesh(3), dp)

    eval_ahc = .false.
    eval_morb = .false.
    eval_kubo = .false.
    eval_sc = .false.
    eval_shc = .false.
    eval_kdotp = .false.

    if (index(pw90_berry%task, 'ahc') > 0) eval_ahc = .true.
    if (index(pw90_berry%task, 'morb') > 0) eval_morb = .true.
    if (index(pw90_berry%task, 'kubo') > 0) eval_kubo = .true.
    if (index(pw90_berry%task, 'sc') > 0) eval_sc = .true.
    if (index(pw90_berry%task, 'shc') > 0) eval_shc = .true.
    if (index(pw90_berry%task, 'kdotp') > 0) eval_kdotp = .true.

    ! Wannier matrix elements, allocations and initializations
    !
    if (eval_ahc) then
      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 (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
      imf_list = 0.0_dp
      adpt_counter_list = 0
    end if

    if (eval_morb) then
      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 (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_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

      imf_list2 = 0.0_dp
      img_list = 0.0_dp
      imh_list = 0.0_dp
    end if

    ! List here berry_tasks that assume nfermi=1
    !
    not_scannable = eval_kubo .or. (eval_shc .and. pw90_spin_hall%freq_scan)
    if (not_scannable .and. fermi_n .ne. 1) then
      call set_error_input(error, 'The berry_task(s, comm, comm) you chose require that you specify a single ' &
                           //'Fermi energy: scanning the Fermi energy is not implemented', comm)
      return
    end if

    if (eval_kubo) then
      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 (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
      allocate (kubo_H_k(3, 3, pw90_berry%kubo_nfreq))
      allocate (kubo_H(3, 3, pw90_berry%kubo_nfreq))
      allocate (kubo_AH_k(3, 3, pw90_berry%kubo_nfreq))
      allocate (kubo_AH(3, 3, pw90_berry%kubo_nfreq))
      allocate (jdos_k(pw90_berry%kubo_nfreq))
      allocate (jdos(pw90_berry%kubo_nfreq))
      kubo_H = cmplx_0
      kubo_AH = cmplx_0
      jdos = 0.0_dp
      if (spin_decomp) 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
        allocate (kubo_H_k_spn(3, 3, 3, pw90_berry%kubo_nfreq))
        allocate (kubo_H_spn(3, 3, 3, pw90_berry%kubo_nfreq))
        allocate (kubo_AH_k_spn(3, 3, 3, pw90_berry%kubo_nfreq))
        allocate (kubo_AH_spn(3, 3, 3, pw90_berry%kubo_nfreq))
        allocate (jdos_k_spn(3, pw90_berry%kubo_nfreq))
        allocate (jdos_spn(3, pw90_berry%kubo_nfreq))
        ! fixme, check these allocs for failure
        kubo_H_spn = cmplx_0
        kubo_AH_spn = cmplx_0
        jdos_spn = 0.0_dp
      end if
    end if

    if (eval_sc) then
      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 (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
      allocate (sc_k_list(3, 6, pw90_berry%kubo_nfreq))
      allocate (sc_list(3, 6, pw90_berry%kubo_nfreq))
      sc_k_list = 0.0_dp
      sc_list = 0.0_dp
    end if

    if (eval_shc) then

      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 (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

      if (index(pw90_spin_hall%method, 'qiao') > 0) then
        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
      else
        call get_SH_R(dis_manifold, kmesh_info, kpt_latt, print_output, pw90_oper_read, &
                      pw90_spin_hall, SH_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
        call get_SAA_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, SS_R, SAA_R, &
                       v_matrix, 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_SBB_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, SH_R, SBB_R, &
                       v_matrix, 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 (pw90_spin_hall%freq_scan) then
        allocate (shc_freq(pw90_berry%kubo_nfreq))
        allocate (shc_k_freq(pw90_berry%kubo_nfreq))
        shc_freq = 0.0_dp
        shc_k_freq = 0.0_dp
      else
        allocate (shc_fermi(fermi_n))
        allocate (shc_k_fermi(fermi_n))
        allocate (shc_k_fermi_dummy(fermi_n))
        shc_fermi = 0.0_dp
        shc_k_fermi = 0.0_dp
        !only used for fermiscan & adpt kmesh
        shc_k_fermi_dummy = 0.0_dp
        adpt_counter_list = 0
      end if

      if (pw90_berry%tetrahedron_method) then
        if (pw90_berry%tetrahedron_higher_correction) then
          allocate (imjv(num_wann, num_wann, 0:pw90_berry%kmesh%mesh(1) + 2, 0:pw90_berry%kmesh%mesh(2) + 2, 0:3))
          allocate (eig(num_wann, 0:pw90_berry%kmesh%mesh(1) + 2, 0:pw90_berry%kmesh%mesh(2) + 2, 0:3))
          allocate (imjv_tet(num_wann, num_wann, 64))
          allocate (eig_tet(num_wann, 64))
          allocate (counts(0:num_nodes - 1))
          allocate (displs(0:num_nodes - 1))
          call tetrahedron_P_matrix_init(P_matrix)
          call tetrahedron_array_init(tet_array)
          if (pw90_spin_hall%freq_scan) then
            nfreq = pw90_berry%kubo_nfreq
          else
            nfreq = fermi_n
          end if
        else !w/o correction: not implemented
          call set_error_input(error, 'Error: tetrahedron method without higher-order correction not implemented', comm)
          !  allocate (imjv(num_wann, num_wann, pw90_berry%kmesh%mesh(1) + 1, pw90_berry%kmesh%mesh(2) + 1, 2))
          !  allocate (eig(num_wann, pw90_berry%kmesh%mesh(1) + 1, pw90_berry%kmesh%mesh(2) + 1, 2))
          !  allocate (imjv_tet(num_wann, num_wann, 8))
          !  allocate (eig_tet(num_wann, 8))
          !  allocate (counts(0:num_nodes - 1))
          !  allocate (displs(0:num_nodes - 1))
          !  !tetrahedron_array_small
        end if
        call comms_array_split(pw90_berry%kmesh%mesh(3), counts, displs, comm)
      end if

    end if

    if (eval_kdotp) then
      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
      kdotp_nbands = size(pw90_berry%kdotp_bands)
      allocate (kdotp(kdotp_nbands, kdotp_nbands, 3, 3, 3))
      kdotp = cmplx_0
    end if

    if (print_output%iprint > 0) then

      write (stdout, '(/,/,1x,a)') &
        'Properties calculated in module  b e r r y'
      write (stdout, '(1x,a)') &
        '------------------------------------------'

      if (eval_ahc) write (stdout, '(/,3x,a)') &
        '* Anomalous Hall conductivity'

      if (eval_morb) write (stdout, '(/,3x,a)') '* Orbital magnetization'

      if (eval_kubo) then
        if (spin_decomp) then
          write (stdout, '(/,3x,a)') &
            '* Complex optical conductivity and its spin-decomposition'
          write (stdout, '(/,3x,a)') &
            '* Joint density of states and its spin-decomposition'
        else
          write (stdout, '(/,3x,a)') '* Complex optical conductivity'
          write (stdout, '(/,3x,a)') '* Joint density of states'
        end if
      end if

      if (eval_sc) write (stdout, '(/,3x,a)') &
        '* Shift current'

      if (eval_shc) then
        write (stdout, '(/,3x,a)') '* Spin Hall Conductivity'
        if (index(pw90_spin_hall%method, 'qiao') > 0) then
          write (stdout, '(/,3x,a)') '  Qiao''s SHC (Phys.Rev.B 98.214402)'
        else
          write (stdout, '(/,3x,a)') '  Ryoo''s SHC (Phys.Rev.B 99.235113)'
        end if
        if (pw90_spin_hall%freq_scan) then
          write (stdout, '(/,3x,a)') '  Frequency scan'
        else
          write (stdout, '(/,3x,a)') '  Fermi energy scan'
        end if
      end if

      if (eval_kdotp) write (stdout, '(/,3x,a)') '* k.p expansion coefficients'

      if (pw90_berry%transl_inv) then
        if (eval_morb) then
          call set_error_input(error, 'transl_inv=T disabled for morb', comm)
          return
        end if
        write (stdout, '(/,1x,a)') 'Using a translationally-invariant discretization for the'
        write (stdout, '(1x,a)') 'band-diagonal Wannier matrix elements of r, etc.'
      end if

      if (pw90_berry%tetrahedron_method) then
        if (pw90_berry%tetrahedron_higher_correction) then
          write (stdout, '(/,3x,a)') '  Tetrahedron method with higher-order correction(PRB 89, 094515)'
        else
          write (stdout, '(/,3x,a)') '  Tetrahedron method without correction(PRB 89, 094515)'
          call set_error_input(error, 'Not yet implemented', comm)
        end if
      end if

      if (print_output%timing_level > 1) then
        call io_stopwatch_stop('berry: prelims', timer)
        call io_stopwatch_start('berry: k-interpolation', timer)
      end if

      if (eval_kdotp) then
        ! JJ pw90_berry%kdotp_bands is only allocated on process 0
        ! this causes a segfault at 2895 (accessing nonexistent element zero)
        ! moving to process 0 (on_root) only
        call berry_get_kdotp(kdotp, dis_manifold, kpt_latt, print_output, pw90_berry, &
                             pw90_band_deriv_degen, wannier_data, ws_distance, wigner_seitz, &
                             ws_region, HH_R, u_matrix, v_matrix, eigval, real_lattice, &
                             scissors_shift, mp_grid, num_bands, num_kpts, num_wann, &
                             num_valence_bands, effective_model, have_disentangled, seedname, &
                             stdout, timer, error, comm)
        if (allocated(error)) return
      end if
    end if ! print_output%iprint > 0, aka "on_root"

    ! Set up adaptive refinement mesh
    !
    allocate (adkpt(3, pw90_berry%curv_adpt_kmesh**3), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating adkpt in berry', comm)
      return
    end if
    ikpt = 0
    !
    ! OLD VERSION (only works correctly for odd grids including original point)
    !
    ! do i=-(pw90_berry_curv_adpt_kmesh-1)/2,(pw90_berry_curv_adpt_kmesh-1)/2
    !    do j=-(pw90_berry_curv_adpt_kmesh-1)/2,(pw90_berry_curv_adpt_kmesh-1)/2
    !       do k=-(pw90_berry_curv_adpt_kmesh-1)/2,(pw90_berry_curv_adpt_kmesh-1)/2
    !          ikpt=ikpt+1
    !          adkpt(1,ikpt)=i*db1/pw90_berry_curv_adpt_kmesh
    !          adkpt(2,ikpt)=j*db2/pw90_berry_curv_adpt_kmesh
    !          adkpt(3,ikpt)=k*db3/pw90_berry_curv_adpt_kmesh
    !       end do
    !    end do
    ! end do
    !
    ! NEW VERSION (both even and odd grids)
    !
    do i = 0, pw90_berry%curv_adpt_kmesh - 1
      do j = 0, pw90_berry%curv_adpt_kmesh - 1
        do k = 0, pw90_berry%curv_adpt_kmesh - 1
          ikpt = ikpt + 1
          adkpt(1, ikpt) = db1*((i + 0.5_dp)/pw90_berry%curv_adpt_kmesh - 0.5_dp)
          adkpt(2, ikpt) = db2*((j + 0.5_dp)/pw90_berry%curv_adpt_kmesh - 0.5_dp)
          adkpt(3, ikpt) = db3*((k + 0.5_dp)/pw90_berry%curv_adpt_kmesh - 0.5_dp)
        end do
      end do
    end do

    ! Loop over interpolation k-points
    !
    if (pw90_berry%wanint_kpoint_file) then

      if (pw90_berry%tetrahedron_method) call set_error_input &
        (error, 'Tetrahedron method not implemented with wanint_kpoint_file', comm)
      if (allocated(error)) return
      ! NOTE: still need to specify pw90_pw90_berry%kmesh%mesh in the input file
      !
      !        - Must use the correct nominal value in order to
      !          correctly set up adaptive smearing in kubo

      if (print_output%iprint > 0) write (stdout, '(/,1x,a,i10,a)') &
        'Reading interpolation grid from file kpoint.dat: ', &
        sum(kpoint_dist%num_int_kpts_on_node), ' points'

      ! Loop over k-points on the irreducible wedge of the Brillouin
      ! zone, read from file 'kpoint.dat'
      !
      do loop_xyz = 1, kpoint_dist%num_int_kpts_on_node(my_node_id)
        kpt(:) = kpoint_dist%int_kpts(:, loop_xyz)
        kweight = kpoint_dist%weight(loop_xyz)
        kweight_adpt = kweight/pw90_berry%curv_adpt_kmesh**3
        !               .
        ! ***BEGIN COPY OF CODE BLOCK 1***
        !
        if (eval_ahc) 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

          ladpt = .false.
          do if = 1, fermi_n
            vdum(1) = sum(imf_k_list(:, 1, if))
            vdum(2) = sum(imf_k_list(:, 2, if))
            vdum(3) = sum(imf_k_list(:, 3, if))
            if (pw90_berry%curv_unit == 'bohr2') vdum = vdum/physics%bohr**2
            rdum = sqrt(dot_product(vdum, vdum))
            if (rdum > pw90_berry%curv_adpt_kmesh_thresh) then
              adpt_counter_list(if) = adpt_counter_list(if) + 1
              ladpt(if) = .true.
            else
              imf_list(:, :, if) = imf_list(:, :, if) + imf_k_list(:, :, if)*kweight
            end if
          end do
          if (any(ladpt)) then
            do loop_adpt = 1, pw90_berry%curv_adpt_kmesh**3
              ! Using imf_k_list here would corrupt values for other
              ! frequencies, hence dummy. Only if-th element is used
              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(:) + adkpt(:, loop_adpt), real_lattice, &
                                       imf_k_list_dummy, scissors_shift, mp_grid, num_bands, &
                                       num_kpts, num_wann, num_valence_bands, effective_model, &
                                       have_disentangled, seedname, stdout, timer, error, comm, &
                                       ladpt=ladpt)
              if (allocated(error)) return

              do if = 1, fermi_n
                if (ladpt(if)) then
                  imf_list(:, :, if) = imf_list(:, :, if) &
                                       + imf_k_list_dummy(:, :, if)*kweight_adpt
                end if
              end do
            end do
          end if
        end if

        if (eval_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

          imf_list2 = imf_list2 + imf_k_list*kweight
          img_list = img_list + img_k_list*kweight
          imh_list = imh_list + imh_k_List*kweight
        end if

        if (eval_kubo) then
          if (spin_decomp) then
            call berry_get_kubo_k(pw90_berry, dis_manifold, fermi_energy_list, kpt_latt, &
                                  pw90_band_deriv_degen, pw90_spin, ws_region, print_output, &
                                  wannier_data, ws_distance, wigner_seitz, AA_R, HH_R, kubo_AH_k, &
                                  kubo_H_k, SS_R, u_matrix, v_matrix, eigval, kpt, real_lattice, &
                                  jdos_k, scissors_shift, mp_grid, num_bands, num_kpts, num_wann, &
                                  num_valence_bands, effective_model, have_disentangled, &
                                  spin_decomp, seedname, stdout, timer, error, comm, &
                                  kubo_AH_k_spn, kubo_H_k_spn, jdos_k_spn)
            if (allocated(error)) return
          else
            call berry_get_kubo_k(pw90_berry, dis_manifold, fermi_energy_list, kpt_latt, &
                                  pw90_band_deriv_degen, pw90_spin, ws_region, print_output, &
                                  wannier_data, ws_distance, wigner_seitz, AA_R, HH_R, kubo_AH_k, &
                                  kubo_H_k, SS_R, u_matrix, v_matrix, eigval, kpt, real_lattice, &
                                  jdos_k, scissors_shift, mp_grid, num_bands, num_kpts, num_wann, &
                                  num_valence_bands, effective_model, have_disentangled, &
                                  spin_decomp, seedname, stdout, timer, error, comm)
            if (allocated(error)) return
          end if
          kubo_H = kubo_H + kubo_H_k*kweight
          kubo_AH = kubo_AH + kubo_AH_k*kweight
          jdos = jdos + jdos_k*kweight
          if (spin_decomp) then
            kubo_H_spn = kubo_H_spn + kubo_H_k_spn*kweight
            kubo_AH_spn = kubo_AH_spn + kubo_AH_k_spn*kweight
            jdos_spn = jdos_spn + jdos_k_spn*kweight
          end if
        end if

        if (eval_sc) then
          call berry_get_sc_klist(pw90_berry, dis_manifold, fermi_energy_list, kmesh_info, &
                                  kpt_latt, ws_region, print_output, pw90_band_deriv_degen, &
                                  wannier_data, ws_distance, wigner_seitz, AA_R, HH_R, u_matrix, &
                                  v_matrix, eigval, kpt, real_lattice, sc_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
          sc_list = sc_list + sc_k_list*kweight
        end if

        ! ***END COPY OF CODE BLOCK 1***

        if (eval_shc) then
          ! print calculation progress, from 0%, 10%, ... to 100%
          ! Note the 1st call to berry_get_shc_klist will be much longer
          ! than later calls due to the time spent on
          !   berry_get_shc_klist -> wham_get_eig_deleig ->
          !   pw90common_fourier_R_to_k -> ws_translate_dist
          if (print_output%iprint > 0) then
            call berry_print_progress(kpoint_dist%num_int_kpts_on_node(my_node_id), loop_xyz, &
                                      1, 1, stdout)
          end if
          if (.not. pw90_spin_hall%freq_scan) 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

            !check whether needs to tigger adpt kmesh or not.
            !Since the calculated shc_k at one Fermi energy can be reused
            !by all the Fermi energies, if we find out that at a specific
            !Fermi energy shc_k(if) > thresh, then we will update shc_k at
            !all the Fermi energies as well.
            !This also avoids repeated calculation if shc_k(if) > thresh
            !is satisfied at more than one Fermi energy.
            ladpt_kmesh = .false.
            !if adpt_kmesh==1, no need to calculate on the same kpt again.
            !This happens if adpt_kmesh==1 while adpt_kmesh_thresh is low.
            if (pw90_berry%curv_adpt_kmesh > 1) then
              do if = 1, fermi_n
                rdum = abs(shc_k_fermi(if))
                if (pw90_berry%curv_unit == 'bohr2') rdum = rdum/physics%bohr**2
                if (rdum > pw90_berry%curv_adpt_kmesh_thresh) then
                  adpt_counter_list(1) = adpt_counter_list(1) + 1
                  ladpt_kmesh = .true.
                  exit
                end if
              end do
            else
              ladpt_kmesh = .false.
            end if
            if (ladpt_kmesh) then
              do loop_adpt = 1, pw90_berry%curv_adpt_kmesh**3
                !Using shc_k here would corrupt values for other
                !kpt, hence dummy. Only if-th element is used.
                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(:) + adkpt(:, loop_adpt), &
                                         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_dummy)
                if (allocated(error)) return

                shc_fermi = shc_fermi + kweight_adpt*shc_k_fermi_dummy
              end do
            else
              shc_fermi = shc_fermi + kweight*shc_k_fermi
            end if
          else ! freq_scan, no adaptive kmesh
            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_freq=shc_k_freq)
            if (allocated(error)) return

            shc_freq = shc_freq + kweight*shc_k_freq
          end if
        end if

      end do !loop_xyz

    else if (.not. pw90_berry%tetrahedron_method) then! Do not read 'kpoint.dat'. Loop over a regular grid in the full BZ

      kweight = db1*db2*db3
      kweight_adpt = kweight/pw90_berry%curv_adpt_kmesh**3

      do loop_xyz = my_node_id, PRODUCT(pw90_berry%kmesh%mesh) - 1, num_nodes
        loop_x = loop_xyz/(pw90_berry%kmesh%mesh(2)*pw90_berry%kmesh%mesh(3))
        loop_y = (loop_xyz - loop_x*(pw90_berry%kmesh%mesh(2) &
                                     *pw90_berry%kmesh%mesh(3)))/pw90_berry%kmesh%mesh(3)
        loop_z = loop_xyz - loop_x*(pw90_berry%kmesh%mesh(2)*pw90_berry%kmesh%mesh(3)) &
                 - loop_y*pw90_berry%kmesh%mesh(3)
        kpt(1) = loop_x*db1
        kpt(2) = loop_y*db2
        kpt(3) = loop_z*db3

        ! ***BEGIN CODE BLOCK 1***
        if (eval_ahc) 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

          ladpt = .false.
          do if = 1, fermi_n
            vdum(1) = sum(imf_k_list(:, 1, if))
            vdum(2) = sum(imf_k_list(:, 2, if))
            vdum(3) = sum(imf_k_list(:, 3, if))
            if (pw90_berry%curv_unit == 'bohr2') vdum = vdum/physics%bohr**2
            rdum = sqrt(dot_product(vdum, vdum))
            if (rdum > pw90_berry%curv_adpt_kmesh_thresh) then
              adpt_counter_list(if) = adpt_counter_list(if) + 1
              ladpt(if) = .true.
            else
              imf_list(:, :, if) = imf_list(:, :, if) + imf_k_list(:, :, if)*kweight
            end if
          end do
          if (any(ladpt)) then
            do loop_adpt = 1, pw90_berry%curv_adpt_kmesh**3
              ! Using imf_k_list here would corrupt values for other
              ! frequencies, hence dummy. Only if-th element is used
              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(:) + adkpt(:, loop_adpt), real_lattice, &
                                       imf_k_list_dummy, scissors_shift, mp_grid, num_bands, &
                                       num_kpts, num_wann, num_valence_bands, effective_model, &
                                       have_disentangled, seedname, stdout, timer, error, comm, &
                                       ladpt=ladpt)
              if (allocated(error)) return

              do if = 1, fermi_n
                if (ladpt(if)) then
                  imf_list(:, :, if) = imf_list(:, :, if) &
                                       + imf_k_list_dummy(:, :, if)*kweight_adpt
                end if
              end do
            end do
          end if
        end if

        if (eval_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

          imf_list2 = imf_list2 + imf_k_list*kweight
          img_list = img_list + img_k_list*kweight
          imh_list = imh_list + imh_k_List*kweight
        end if

        if (eval_kubo) then
          if (spin_decomp) then
            call berry_get_kubo_k(pw90_berry, dis_manifold, fermi_energy_list, kpt_latt, &
                                  pw90_band_deriv_degen, pw90_spin, ws_region, print_output, &
                                  wannier_data, ws_distance, wigner_seitz, AA_R, HH_R, kubo_AH_k, &
                                  kubo_H_k, SS_R, u_matrix, v_matrix, eigval, kpt, real_lattice, &
                                  jdos_k, scissors_shift, mp_grid, num_bands, num_kpts, num_wann, &
                                  num_valence_bands, effective_model, have_disentangled, &
                                  spin_decomp, seedname, stdout, timer, error, comm, &
                                  kubo_AH_k_spn, kubo_H_k_spn, jdos_k_spn)
            if (allocated(error)) return

          else
            call berry_get_kubo_k(pw90_berry, dis_manifold, fermi_energy_list, kpt_latt, &
                                  pw90_band_deriv_degen, pw90_spin, ws_region, print_output, &
                                  wannier_data, ws_distance, wigner_seitz, AA_R, HH_R, kubo_AH_k, &
                                  kubo_H_k, SS_R, u_matrix, v_matrix, eigval, kpt, real_lattice, &
                                  jdos_k, scissors_shift, mp_grid, num_bands, num_kpts, num_wann, &
                                  num_valence_bands, effective_model, have_disentangled, &
                                  spin_decomp, seedname, stdout, timer, error, comm)
            if (allocated(error)) return

          end if
          kubo_H = kubo_H + kubo_H_k*kweight
          kubo_AH = kubo_AH + kubo_AH_k*kweight
          jdos = jdos + jdos_k*kweight
          if (spin_decomp) then
            kubo_H_spn = kubo_H_spn + kubo_H_k_spn*kweight
            kubo_AH_spn = kubo_AH_spn + kubo_AH_k_spn*kweight
            jdos_spn = jdos_spn + jdos_k_spn*kweight
          end if
        end if

        if (eval_sc) then
          call berry_get_sc_klist(pw90_berry, dis_manifold, fermi_energy_list, kmesh_info, &
                                  kpt_latt, ws_region, print_output, pw90_band_deriv_degen, &
                                  wannier_data, ws_distance, wigner_seitz, AA_R, HH_R, u_matrix, &
                                  v_matrix, eigval, kpt, real_lattice, sc_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

          sc_list = sc_list + sc_k_list*kweight
        end if

        ! ***END CODE BLOCK 1***

        if (eval_shc) then
          ! print calculation progress, from 0%, 10%, ... to 100%
          ! Note the 1st call to berry_get_shc_klist will be much longer
          ! than later calls due to the time spent on
          !   berry_get_shc_klist -> wham_get_eig_deleig ->
          !   pw90common_fourier_R_to_k -> ws_translate_dist
          if (print_output%iprint > 0) then
            call berry_print_progress(PRODUCT(pw90_berry%kmesh%mesh) - 1, loop_xyz, my_node_id, &
                                      num_nodes, stdout)
          end if
          if (.not. pw90_spin_hall%freq_scan) 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

            !check whether needs to tigger adpt kmesh or not.
            !Since the calculated shc_k at one Fermi energy can be reused
            !by all the Fermi energies, if we find out that at a specific
            !Fermi energy shc_k(if) > thresh, then we will update shc_k at
            !all the Fermi energies as well.
            !This also avoids repeated calculation if shc_k(if) > thresh
            !is satisfied at more than one Fermi energy.
            ladpt_kmesh = .false.
            !if adpt_kmesh==1, no need to calculate on the same kpt again.
            !This happens if adpt_kmesh==1 while adpt_kmesh_thresh is low.
            if (pw90_berry%curv_adpt_kmesh > 1) then
              do if = 1, fermi_n
                rdum = abs(shc_k_fermi(if))
                if (pw90_berry%curv_unit == 'bohr2') rdum = rdum/physics%bohr**2
                if (rdum > pw90_berry%curv_adpt_kmesh_thresh) then
                  adpt_counter_list(1) = adpt_counter_list(1) + 1
                  ladpt_kmesh = .true.
                  exit
                end if
              end do
            else
              ladpt_kmesh = .false.
            end if
            if (ladpt_kmesh) then
              do loop_adpt = 1, pw90_berry%curv_adpt_kmesh**3
                !Using shc_k here would corrupt values for other
                !kpt, hence dummy. Only if-th element is used.
                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(:) + adkpt(:, loop_adpt), &
                                         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_dummy)
                if (allocated(error)) return

                shc_fermi = shc_fermi + kweight_adpt*shc_k_fermi_dummy
              end do
            else
              shc_fermi = shc_fermi + kweight*shc_k_fermi
            end if
          else ! freq_scan, no adaptive kmesh
            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_freq=shc_k_freq)
            if (allocated(error)) return

            shc_freq = shc_freq + kweight*shc_k_freq
          end if
        end if

      end do !loop_xyz

    else !tetrahedron_method
      if (eval_shc) then
        ! to do: move it to tetrahedron.f90?
        do loop_z = displs(my_node_id), displs(my_node_id) + counts(my_node_id) - 1 ! current implementation limits max of nodes = pw90_berry%kmesh%mesh(3)
          !obtaining energy eigenvalues and matrix elements
          !parallelization by loop_z
          if (print_output%iprint > 0) &
            write (stdout, '(a,i0,a,i0)') "! obtaining energy eigenvalues and matrix elements... ", &
            loop_z - displs(my_node_id) + 1, "/", counts(my_node_id)
          do loop_x = -1, pw90_berry%kmesh%mesh(1) + 1
            do loop_y = -1, pw90_berry%kmesh%mesh(2) + 1
              kpt(1) = loop_x*db1 + mesh_shift*db1; kpt(2) = loop_y*db2 + mesh_shift*db2
              !boundary of BZ
              if (kpt(1) < 0) kpt(1) = kpt(1) + pw90_berry%kmesh%mesh(1)*db1
              if (kpt(2) < 0) kpt(2) = kpt(2) + pw90_berry%kmesh%mesh(2)*db2
              if (kpt(1) > pw90_berry%kmesh%mesh(1)*db1 - 0.001_dp*db1) kpt(1) = kpt(1) - pw90_berry%kmesh%mesh(1)*db1
              if (kpt(2) > pw90_berry%kmesh%mesh(2)*db2 - 0.001_dp*db2) kpt(2) = kpt(2) - pw90_berry%kmesh%mesh(2)*db2
              !optimized tetrahedron method(Kawamura, PRB 89 094515)
              !1st layer for points 6, 9, 10, 13 (in Fig. 3)
              !2nd layer for points 1, 2, 5, 14, 19, 20
              !3rd layer for points 3, 4, 7, 16, 17, 18
              !4th layer for points 8, 11, 12, 15
              do i = 0, 3 ! main tetrahedra are between 2nd(i=1) and 3rd(i=2) layers
                kpt(3) = (loop_z - 1 + i)*db3 + mesh_shift*db3
                if (kpt(3) < 0) kpt(3) = kpt(3) + pw90_berry%kmesh%mesh(3)*db3 !boundary of BZ
                if (kpt(3) > pw90_berry%kmesh%mesh(3)*db3 - 0.001_dp*db3) kpt(3) = kpt(3) - pw90_berry%kmesh%mesh(3)*db3 !boundary of BZ
                if (loop_z == displs(my_node_id) .or. i == 3) then
                  call berry_get_shc_tetrahedron(pw90_berry, ws_region, pw90_spin_hall, wannier_data, ws_distance, &
                                                 wigner_seitz, AA_R, HH_R, SH_R, SHR_R, SR_R, SS_R, SAA_R, SBB_R, &
                                                 kpt, imjv(:, :, loop_x + 1, loop_y + 1, i), eig(:, loop_x + 1, loop_y + 1, i), &
                                                 real_lattice, mp_grid, num_wann, seedname, stdout, error, comm)
                else
                  imjv(:, :, loop_x + 1, loop_y + 1, i) = imjv(:, :, loop_x + 1, loop_y + 1, i + 1)
                  eig(:, loop_x + 1, loop_y + 1, i) = eig(:, loop_x + 1, loop_y + 1, i + 1)
                end if
              end do
            end do
          end do
          !summation
          do loop_x = 0, pw90_berry%kmesh%mesh(1) - 1
            do loop_y = 0, pw90_berry%kmesh%mesh(2) - 1
              ! writing progress - summation is the main bottleneck
              loop_xyz = (loop_z - displs(my_node_id))*pw90_berry%kmesh%mesh(1)*pw90_berry%kmesh%mesh(2) &
                         + loop_x*pw90_berry%kmesh%mesh(2) + loop_y

              if (print_output%iprint > 0) then ! only print from root
                call berry_print_progress(loop_xyz, 0, counts(my_node_id)*pw90_berry%kmesh%mesh(1) &
                                          *pw90_berry%kmesh%mesh(2) - 1, 1, stdout)
              end if

              ! setting 8 vertices and surrounding points
              do i = 0, 3 !16*l+4*k+i+1 = 1,2,3,...,64, eight vertices of a mesh(22, 23, 26, 27, 38, 39, 42, 43) and their surrounding points
                do k = 0, 3
                  do l = 0, 3
                    kptc(1, 16*l + 4*k + i + 1) = (loop_x + i - 1)*db1 + mesh_shift*db1
                    kptc(2, 16*l + 4*k + i + 1) = (loop_y + k - 1)*db2 + mesh_shift*db2
                    kptc(3, 16*l + 4*k + i + 1) = (loop_z + l - 1)*db3 + mesh_shift*db3
                    imjv_tet(:, :, 16*l + 4*k + i + 1) = imjv(:, :, loop_x + i, loop_y + k, l)
                    eig_tet(:, 16*l + 4*k + i + 1) = eig(:, loop_x + i, loop_y + k, l)
                  end do
                end do
              end do
              do itet = 1, 6 ! 6 tetrahedra
                do i = 1, 4 ! four vertices
                  kptv(i, :) = kptc(:, tet_array(itet, i))
                end do
                do i = 1, 3 ! xyz
                  do k = 1, 3 ! three vectors forming a tetrahedron
                    ttet(i, k) = kptv(k + 1, i) - kptv(1, i)
                  end do
                end do

                do n = 1, num_wann
                  do m = 1, num_wann
                    if (n == m) cycle
                    do i = 1, 20 ! four vertices + 16 points for optimization
                      F_opt(i) = imjv_tet(n, m, tet_array(itet, i))
                      E1_opt(i) = eig_tet(n, tet_array(itet, i))
                      E2_opt(i) = eig_tet(m, tet_array(itet, i))
                    end do

                    E1tet = 0.0_dp
                    E2tet = 0.0_dp
                    Ftet = 0.0_dp
                    do i = 1, 4
                      do k = 1, 20 !Eq. (16) of Kawamura, PRB 89 094515
                        E1tet(i) = E1tet(i) + P_matrix(i, k)*E1_opt(k)
                        E2tet(i) = E2tet(i) + P_matrix(i, k)*E2_opt(k)
                        Ftet(i) = Ftet(i) + P_matrix(i, k)*F_opt(k)
                      end do
                    end do

                    ! do i = 1, 4
                    !   Ftet(i) = imjv_tet(n, m, tet_array(itet, i))
                    !   E1tet(i) = eig_tet(n, tet_array(itet, i))
                    !   E2tet(i) = eig_tet(m, tet_array(itet, i))
                    ! enddo
                    do ifreq = 1, nfreq !fermiscan or freqscan
                      if (.not. pw90_spin_hall%freq_scan) then
                        omega = real(pw90_berry%kubo_freq_list(1), dp)
                        Ef = fermi_energy_list(ifreq)
                      else
                        omega = real(pw90_berry%kubo_freq_list(ifreq), dp)
                        Ef = fermi_energy_list(1)
                      end if

                      if (omega == 0.0) then
                        shc_k_tet = &
                          tetrahedron_spinhall(Ftet, E1tet, E2tet, ttet, &
                                               0.0_dp, Ef, 3, pw90_berry%tetrahedron_cutoff, &
                                               pw90_berry%tetrahedron_avoid_degeneracy)
                      else
                        shc_k_tet = &
                          (tetrahedron_spinhall(Ftet, E1tet, E2tet, ttet, &
                                                -omega, Ef, 1, pw90_berry%tetrahedron_cutoff, &
                                                pw90_berry%tetrahedron_avoid_degeneracy) &
                           - tetrahedron_spinhall(Ftet, E1tet, E2tet, ttet, &
                                                  omega, Ef, 1, pw90_berry%tetrahedron_cutoff, &
                                                  pw90_berry%tetrahedron_avoid_degeneracy)) &
                          /(2.0_dp*omega) + (cmplx_i*pi* &
                                             (tetrahedron_spinhall(Ftet, E1tet, E2tet, ttet, &
                                                                   -omega, Ef, 2, pw90_berry%tetrahedron_cutoff, &
                                                                   pw90_berry%tetrahedron_avoid_degeneracy) &
                                              + tetrahedron_spinhall(Ftet, E1tet, E2tet, ttet, &
                                                                     omega, Ef, 2, pw90_berry%tetrahedron_cutoff, &
                                                                     pw90_berry%tetrahedron_avoid_degeneracy))) &
                          /(2.0_dp*omega)
                      end if

                      if (.not. pw90_spin_hall%freq_scan) then
                        shc_fermi(ifreq) = shc_fermi(ifreq) - real(shc_k_tet, dp)
                      else
                        shc_freq(ifreq) = shc_freq(ifreq) - shc_k_tet
                      end if
                    end do !ifreq
                  end do ! m
                end do ! n
              end do ! itet
            end do ! loop_y for summation
          end do ! loop_x for summation
        end do ! loop_z
      end if
    end if !wanint_kpoint_file

    ! Collect contributions from all nodes
    if (eval_ahc) then
      call comms_reduce(imf_list(1, 1, 1), 3*3*fermi_n, 'SUM', error, comm)
      if (allocated(error)) return
      call comms_reduce(adpt_counter_list(1), fermi_n, 'SUM', error, comm)
      if (allocated(error)) return
    end if

    if (eval_morb) then
      call comms_reduce(imf_list2(1, 1, 1), 3*3*fermi_n, 'SUM', error, comm)
      if (allocated(error)) return
      call comms_reduce(img_list(1, 1, 1), 3*3*fermi_n, 'SUM', error, comm)
      if (allocated(error)) return
      call comms_reduce(imh_list(1, 1, 1), 3*3*fermi_n, 'SUM', error, comm)
      if (allocated(error)) return
    end if

    if (eval_kubo) then
      call comms_reduce(kubo_H(1, 1, 1), 3*3*pw90_berry%kubo_nfreq, 'SUM', error, comm)
      if (allocated(error)) return
      call comms_reduce(kubo_AH(1, 1, 1), 3*3*pw90_berry%kubo_nfreq, 'SUM', error, comm)
      if (allocated(error)) return
      call comms_reduce(jdos(1), pw90_berry%kubo_nfreq, 'SUM', error, comm)
      if (allocated(error)) return
      if (spin_decomp) then
        call comms_reduce(kubo_H_spn(1, 1, 1, 1), 3*3*3*pw90_berry%kubo_nfreq, 'SUM', error, comm)
        if (allocated(error)) return
        call comms_reduce(kubo_AH_spn(1, 1, 1, 1), 3*3*3*pw90_berry%kubo_nfreq, 'SUM', error, comm)
        if (allocated(error)) return
        call comms_reduce(jdos_spn(1, 1), 3*pw90_berry%kubo_nfreq, 'SUM', error, comm)
        if (allocated(error)) return
      end if
    end if

    if (eval_sc) then
      call comms_reduce(sc_list(1, 1, 1), 3*6*pw90_berry%kubo_nfreq, 'SUM', error, comm)
      if (allocated(error)) return
    end if

    if (eval_shc) then
      if (pw90_spin_hall%freq_scan) then
        call comms_reduce(shc_freq(1), pw90_berry%kubo_nfreq, 'SUM', error, comm)
        if (allocated(error)) return
      else
        call comms_reduce(shc_fermi(1), fermi_n, 'SUM', error, comm)
        if (allocated(error)) return
        call comms_reduce(adpt_counter_list(1), fermi_n, 'SUM', error, comm)
        if (allocated(error)) return
      end if
    end if

    if (print_output%iprint > 0) then

      if (print_output%timing_level > 1) call io_stopwatch_stop('berry: k-interpolation', timer)
      write (stdout, '(1x,a)') ' '
      if (eval_ahc .and. pw90_berry%curv_adpt_kmesh .ne. 1) then
        if (.not. pw90_berry%wanint_kpoint_file) write (stdout, '(1x,a28,3(i0,1x))') &
          'Regular interpolation grid: ', pw90_berry%kmesh%mesh
        write (stdout, '(1x,a28,3(i0,1x))') 'Adaptive refinement grid: ', &
          pw90_berry%curv_adpt_kmesh, pw90_berry%curv_adpt_kmesh, pw90_berry%curv_adpt_kmesh
        if (pw90_berry%curv_unit == 'ang2') then
          write (stdout, '(1x,a28,a17,f6.2,a)') &
            'Refinement threshold: ', 'Berry curvature >', &
            pw90_berry%curv_adpt_kmesh_thresh, ' Ang^2'
        elseif (pw90_berry%curv_unit == 'bohr2') then
          write (stdout, '(1x,a28,a17,f6.2,a)') &
            'Refinement threshold: ', 'Berry curvature >', &
            pw90_berry%curv_adpt_kmesh_thresh, ' bohr^2'
        end if
        if (fermi_n == 1) then
          if (pw90_berry%wanint_kpoint_file) then
            write (stdout, '(1x,a30,i5,a,f5.2,a)') &
              ' Points triggering refinement: ', &
              adpt_counter_list(1), '(', &
              100*real(adpt_counter_list(1), dp) &
              /sum(kpoint_dist%num_int_kpts_on_node), '%)'
          else
            write (stdout, '(1x,a30,i5,a,f5.2,a)') &
              ' Points triggering refinement: ', &
              adpt_counter_list(1), '(', &
              100*real(adpt_counter_list(1), dp)/product(pw90_berry%kmesh%mesh), '%)'
          end if
        end if
      elseif (eval_shc) then
        if (pw90_berry%curv_adpt_kmesh .ne. 1) then
          if (.not. pw90_berry%wanint_kpoint_file) write (stdout, '(1x,a28,3(i0,1x))') &
            'Regular interpolation grid: ', pw90_berry%kmesh%mesh
          if (.not. pw90_spin_hall%freq_scan) then
            write (stdout, '(1x,a28,3(i0,1x))') &
              'Adaptive refinement grid: ', &
              pw90_berry%curv_adpt_kmesh, pw90_berry%curv_adpt_kmesh, pw90_berry%curv_adpt_kmesh
            if (pw90_berry%curv_unit == 'ang2') then
              write (stdout, '(1x,a28,f12.2,a)') &
                'Refinement threshold: ', &
                pw90_berry%curv_adpt_kmesh_thresh, ' Ang^2'
            elseif (pw90_berry%curv_unit == 'bohr2') then
              write (stdout, '(1x,a28,f12.2,a)') &
                'Refinement threshold: ', &
                pw90_berry%curv_adpt_kmesh_thresh, ' bohr^2'
            end if
            if (pw90_berry%wanint_kpoint_file) then
              write (stdout, '(1x,a30,i8,a,f6.2,a)') &
                ' Points triggering refinement: ', adpt_counter_list(1), '(', &
                100*real(adpt_counter_list(1), dp)/sum(kpoint_dist%num_int_kpts_on_node), '%)'
            else
              write (stdout, '(1x,a30,i8,a,f6.2,a)') &
                ' Points triggering refinement: ', adpt_counter_list(1), '(', &
                100*real(adpt_counter_list(1), dp)/product(pw90_berry%kmesh%mesh), '%)'
            end if
          end if
        else
          if (.not. pw90_berry%wanint_kpoint_file) write (stdout, &
                                                          '(1x,a20,3(i0,1x))') 'Interpolation grid: ', pw90_berry%kmesh%mesh(1:3)
        end if
        write (stdout, '(a)') ''
        if (pw90_berry%kubo_smearing%use_adaptive) then
          write (stdout, '(1x,a)') 'Using adaptive smearing'
          write (stdout, '(7x,a,f8.3)') 'adaptive smearing prefactor ', &
            pw90_berry%kubo_smearing%adaptive_prefactor
          write (stdout, '(7x,a,f8.3,a)') 'adaptive smearing max width ', &
            pw90_berry%kubo_smearing%adaptive_max_width, ' eV'
        else
          write (stdout, '(1x,a)') 'Using fixed smearing'
          write (stdout, '(7x,a,f8.3,a)') 'fixed smearing width ', &
            pw90_berry%kubo_smearing%fixed_width, ' eV'
        end if
        write (stdout, '(a)') ''
        if (abs(scissors_shift) > 1.0e-7_dp) then
          write (stdout, '(1X,A,I0,A,G18.10,A)') "Using scissors_shift to shift energy bands with index > ", &
            num_valence_bands, " by ", scissors_shift, " eV."
        end if
        if (pw90_spin_hall%bandshift) then
          write (stdout, '(1X,A,I0,A,G18.10,A)') "Using shc_bandshift to shift energy bands with index >= ", &
            pw90_spin_hall%bandshift_firstband, " by ", pw90_spin_hall%bandshift_energyshift, " eV."
        end if
      else
        if (.not. pw90_berry%wanint_kpoint_file) write (stdout, &
                                                        '(1x,a20,3(i0,1x))') 'Interpolation grid: ', pw90_berry%kmesh%mesh(1:3)
      end if

      if (eval_ahc) then
        !
        ! --------------------------------------------------------------------
        ! At this point imf contains
        !
        ! (1/N) sum_k Omega_{alpha beta}(k),
        !
        ! an approximation to
        !
        ! V_c.int dk/(2.pi)^3 Omega_{alpha beta}(k) dk
        !
        ! (V_c is the cell volume). We want
        !
        ! sigma_{alpha beta}=-(e^2/hbar) int dk/(2.pi)^3 Omega(k) dk
        !
        ! Hence need to multiply by -(e^2/hbar.V_c).
        ! To get a conductivity in units of S/cm,
        !
        ! (i)   Divide by V_c to obtain (1/N) sum_k omega(k)/V_c, with units
        !       of [L]^{-1} (Berry curvature Omega(k) has units of [L]^2)
        ! (ii)  [L] = Angstrom. Multiply by 10^8 to convert to (cm)^{-1}
        ! (iii) Multiply by -e^2/hbar in SI, with has units ofconductance,
        !       (Ohm)^{-1}, or Siemens (S), to get the final result in S/cm
        !
        !==================================================
        ! fac = -e^2/(hbar.V_c*10^-8)
        !==================================================
        !
        ! with 'V_c' in Angstroms^3, and 'e', 'hbar' in SI units
        ! --------------------------------------------------------------------
        !
        fac = -1.0e8_dp*physics%elem_charge_SI**2/(physics%hbar_SI*cell_volume)
        ahc_list(:, :, :) = imf_list(:, :, :)*fac
        if (fermi_n > 1) then
          write (stdout, '(/,1x,a)') &
            '---------------------------------'
          write (stdout, '(1x,a)') &
            'Output data files related to AHC:'
          write (stdout, '(1x,a)') &
            '---------------------------------'
          file_name = trim(seedname)//'-ahc-fermiscan.dat'
          write (stdout, '(/,3x,a)') '* '//file_name
          open (newunit=file_unit, FILE=file_name, STATUS='UNKNOWN', FORM='FORMATTED')
        end if
        do if = 1, fermi_n
          if (fermi_n > 1) write (file_unit, '(4(F12.6,1x))') &
            fermi_energy_list(if), sum(ahc_list(:, 1, if)), &
            sum(ahc_list(:, 2, if)), sum(ahc_list(:, 3, if))
          write (stdout, '(/,1x,a18,F10.4)') 'Fermi energy (ev):', &
            fermi_energy_list(if)
          if (fermi_n > 1) then
            if (pw90_berry%wanint_kpoint_file) then
              write (stdout, '(1x,a30,i5,a,f5.2,a)') &
                ' Points triggering refinement: ', &
                adpt_counter_list(if), '(', &
                100*real(adpt_counter_list(if), dp) &
                /sum(kpoint_dist%num_int_kpts_on_node), '%)'
            else
              write (stdout, '(1x,a30,i5,a,f5.2,a)') &
                ' Points triggering refinement: ', &
                adpt_counter_list(if), '(', &
                100*real(adpt_counter_list(if), dp) &
                /product(pw90_berry%kmesh%mesh), '%)'
            end if
          end if
          write (stdout, '(/,1x,a)') &
            'AHC (S/cm)       x          y          z'
          if (print_output%iprint > 1) then
            write (stdout, '(1x,a)') &
              '=========='
            write (stdout, '(1x,a9,2x,3(f10.4,1x))') 'J0 term :', &
              ahc_list(1, 1, if), ahc_list(1, 2, if), ahc_list(1, 3, if)
            write (stdout, '(1x,a9,2x,3(f10.4,1x))') 'J1 term :', &
              ahc_list(2, 1, if), ahc_list(2, 2, if), ahc_list(2, 3, if)
            write (stdout, '(1x,a9,2x,3(f10.4,1x))') 'J2 term :', &
              ahc_list(3, 1, if), ahc_list(3, 2, if), ahc_list(3, 3, if)
            write (stdout, '(1x,a)') &
              '-------------------------------------------'
            write (stdout, '(1x,a9,2x,3(f10.4,1x),/)') 'Total   :', &
              sum(ahc_list(:, 1, if)), sum(ahc_list(:, 2, if)), &
              sum(ahc_list(:, 3, if))
          else
            write (stdout, '(1x,a10,1x,3(f10.4,1x),/)') '==========', &
              sum(ahc_list(:, 1, if)), sum(ahc_list(:, 2, if)), &
              sum(ahc_list(:, 3, if))
          end if
        end do
        if (fermi_n > 1) close (file_unit)
      end if

      if (eval_morb) then
        !
        ! --------------------------------------------------------------------
        ! At this point X=img_ab(:)-fermi_energy*imf_ab(:) and
        !               Y=imh_ab(:)-fermi_energy*imf_ab(:)
        ! contain, eg,
        !
        ! (1/N) sum_k X(k), where X(k)=-2*Im[g(k)-E_F.f(k)]
        !
        ! This is an approximation to
        !
        ! V_c.int dk/(2.pi)^3 X(k) dk
        !
        ! (V_c is the cell volume). We want a magnetic moment per cell,
        ! in units of the Bohr magneton. The magnetization-like quantity is
        !
        ! \tilde{M}^LC=-(e/2.hbar) int dk/(2.pi)^3 X(k) dk
        !
        ! So we take X and
        !
        !  (i)  The summand is an energy in eV times a Berry curvature in
        !       Ang^2. To convert to a.u., divide by 27.2 and by 0.529^2
        !  (ii) Multiply by -(e/2.hbar)=-1/2 in atomic units
        ! (iii) At this point we have a magnetic moment (per cell) in atomic
        !       units. 1 Bohr magneton = 1/2 atomic unit, so need to multiply
        !       by 2 to convert it to Bohr magnetons
        ! --------------------------------------------------------------------
        !
        fac = -physics%eV_au/physics%bohr**2
        if (fermi_n > 1) then
          write (stdout, '(/,1x,a)') &
            '---------------------------------'
          write (stdout, '(1x,a)') &
            'Output data files related to the orbital magnetization:'
          write (stdout, '(1x,a)') &
            '---------------------------------'
          file_name = trim(seedname)//'-morb-fermiscan.dat'
          write (stdout, '(/,3x,a)') '* '//file_name
          open (newunit=file_unit, FILE=file_name, STATUS='UNKNOWN', FORM='FORMATTED')
        end if
        do if = 1, fermi_n
          LCtil_list(:, :, if) = (img_list(:, :, if) &
                                  - fermi_energy_list(if)*imf_list2(:, :, if))*fac
          ICtil_list(:, :, if) = (imh_list(:, :, if) &
                                  - fermi_energy_list(if)*imf_list2(:, :, if))*fac
          Morb_list(:, :, if) = LCtil_list(:, :, if) + ICtil_list(:, :, if)
          if (fermi_n > 1) write (file_unit, '(4(F12.6,1x))') &
            fermi_energy_list(if), sum(Morb_list(1:3, 1, if)), &
            sum(Morb_list(1:3, 2, if)), sum(Morb_list(1:3, 3, if))
          write (stdout, '(/,/,1x,a,F12.6)') 'Fermi energy (ev) =', &
            fermi_energy_list(if)
          write (stdout, '(/,/,1x,a)') &
            'M_orb (bohr magn/cell)        x          y          z'
          if (print_output%iprint > 1) then
            write (stdout, '(1x,a)') &
              '======================'
            write (stdout, '(1x,a22,2x,3(f16.10,1x))') 'Local circulation :', &
              sum(LCtil_list(1:3, 1, if)), sum(LCtil_list(1:3, 2, if)), &
              sum(LCtil_list(1:3, 3, if))
            write (stdout, '(1x,a22,2x,3(f16.10,1x))') &
              'Itinerant circulation:', &
              sum(ICtil_list(1:3, 1, if)), sum(ICtil_list(1:3, 2, if)), &
              sum(ICtil_list(1:3, 3, if))
            write (stdout, '(1x,a)') &
              '--------------------------------------------------------'
            write (stdout, '(1x,a22,2x,3(f16.10,1x),/)') 'Total   :', &
              sum(Morb_list(1:3, 1, if)), sum(Morb_list(1:3, 2, if)), &
              sum(Morb_list(1:3, 3, if))
          else
            write (stdout, '(1x,a22,2x,3(f16.10,1x),/)') &
              '======================', &
              sum(Morb_list(1:3, 1, if)), sum(Morb_list(1:3, 2, if)), &
              sum(Morb_list(1:3, 3, if))
          end if
        end do
        if (fermi_n > 1) close (file_unit)
      end if

      ! -----------------------------!
      ! Complex optical conductivity !
      ! -----------------------------!
      !
      if (eval_kubo) then
        !
        ! Convert to S/cm
        fac = 1.0e8_dp*physics%elem_charge_SI**2/(physics%hbar_SI*cell_volume)
        kubo_H = kubo_H*fac
        kubo_AH = kubo_AH*fac
        if (spin_decomp) then
          kubo_H_spn = kubo_H_spn*fac
          kubo_AH_spn = kubo_AH_spn*fac
        end if
        !
        write (stdout, '(/,1x,a)') &
          '----------------------------------------------------------'
        write (stdout, '(1x,a)') &
          'Output data files related to complex optical conductivity:'
        write (stdout, '(1x,a)') &
          '----------------------------------------------------------'
        !
        ! Symmetric: real (imaginary) part is Hermitean (anti-Hermitean)
        !
        do n = 1, 6
          i = alpha_S(n)
          j = beta_S(n)
          file_name = trim(seedname)//'-kubo_S_'// &
                      achar(119 + i)//achar(119 + j)//'.dat'
          file_name = trim(file_name)
          write (stdout, '(/,3x,a)') '* '//file_name
          open (newunit=file_unit, FILE=file_name, STATUS='UNKNOWN', FORM='FORMATTED')
          do ifreq = 1, pw90_berry%kubo_nfreq
            if (spin_decomp) then
              write (file_unit, '(9E16.8)') real(pw90_berry%kubo_freq_list(ifreq), dp), &
                real(0.5_dp*(kubo_H(i, j, ifreq) + kubo_H(j, i, ifreq)), dp), &
                aimag(0.5_dp*(kubo_AH(i, j, ifreq) + kubo_AH(j, i, ifreq))), &
                real(0.5_dp*(kubo_H_spn(i, j, 1, ifreq) &
                             + kubo_H_spn(j, i, 1, ifreq)), dp), &
                aimag(0.5_dp*(kubo_AH_spn(i, j, 1, ifreq) &
                              + kubo_AH_spn(j, i, 1, ifreq))), &
                real(0.5_dp*(kubo_H_spn(i, j, 2, ifreq) &
                             + kubo_H_spn(j, i, 2, ifreq)), dp), &
                aimag(0.5_dp*(kubo_AH_spn(i, j, 2, ifreq) &
                              + kubo_AH_spn(j, i, 2, ifreq))), &
                real(0.5_dp*(kubo_H_spn(i, j, 3, ifreq) &
                             + kubo_H_spn(j, i, 3, ifreq)), dp), &
                aimag(0.5_dp*(kubo_AH_spn(i, j, 3, ifreq) &
                              + kubo_AH_spn(j, i, 3, ifreq)))
            else
              write (file_unit, '(3E16.8)') real(pw90_berry%kubo_freq_list(ifreq), dp), &
                real(0.5_dp*(kubo_H(i, j, ifreq) + kubo_H(j, i, ifreq)), dp), &
                aimag(0.5_dp*(kubo_AH(i, j, ifreq) + kubo_AH(j, i, ifreq)))
            end if
          end do
          close (file_unit)
        end do
        !
        ! Antisymmetric: real (imaginary) part is anti-Hermitean (Hermitean)
        !
        do n = 1, 3
          i = alpha_A(n)
          j = beta_A(n)
          file_name = trim(seedname)//'-kubo_A_'// &
                      achar(119 + i)//achar(119 + j)//'.dat'
          file_name = trim(file_name)
          write (stdout, '(/,3x,a)') '* '//file_name
          open (newunit=file_unit, FILE=file_name, STATUS='UNKNOWN', FORM='FORMATTED')
          do ifreq = 1, pw90_berry%kubo_nfreq
            if (spin_decomp) then
              write (file_unit, '(9E16.8)') real(pw90_berry%kubo_freq_list(ifreq), dp), &
                real(0.5_dp*(kubo_AH(i, j, ifreq) - kubo_AH(j, i, ifreq)), dp), &
                aimag(0.5_dp*(kubo_H(i, j, ifreq) - kubo_H(j, i, ifreq))), &
                real(0.5_dp*(kubo_AH_spn(i, j, 1, ifreq) &
                             - kubo_AH_spn(j, i, 1, ifreq)), dp), &
                aimag(0.5_dp*(kubo_H_spn(i, j, 1, ifreq) &
                              - kubo_H_spn(j, i, 1, ifreq))), &
                real(0.5_dp*(kubo_AH_spn(i, j, 2, ifreq) &
                             - kubo_AH_spn(j, i, 2, ifreq)), dp), &
                aimag(0.5_dp*(kubo_H_spn(i, j, 2, ifreq) &
                              - kubo_H_spn(j, i, 2, ifreq))), &
                real(0.5_dp*(kubo_AH_spn(i, j, 3, ifreq) &
                             - kubo_AH_spn(j, i, 3, ifreq)), dp), &
                aimag(0.5_dp*(kubo_H_spn(i, j, 3, ifreq) &
                              - kubo_H_spn(j, i, 3, ifreq)))
            else
              write (file_unit, '(3E16.8)') real(pw90_berry%kubo_freq_list(ifreq), dp), &
                real(0.5_dp*(kubo_AH(i, j, ifreq) - kubo_AH(j, i, ifreq)), dp), &
                aimag(0.5_dp*(kubo_H(i, j, ifreq) - kubo_H(j, i, ifreq)))
            end if
          end do
          close (file_unit)
        end do
        !
        ! Joint density of states
        !
        file_name = trim(seedname)//'-jdos.dat'
        write (stdout, '(/,3x,a)') '* '//file_name
        open (newunit=file_unit, FILE=file_name, STATUS='UNKNOWN', FORM='FORMATTED')
        do ifreq = 1, pw90_berry%kubo_nfreq
          if (spin_decomp) then
            write (file_unit, '(5E16.8)') real(pw90_berry%kubo_freq_list(ifreq), dp), &
              jdos(ifreq), jdos_spn(:, ifreq)
          else
            write (file_unit, '(2E16.8)') real(pw90_berry%kubo_freq_list(ifreq), dp), &
              jdos(ifreq)
          end if
        end do
        close (file_unit)
      end if

      if (eval_sc) then
        ! -----------------------------!
        ! Nonlinear shift current
        ! -----------------------------!

        ! --------------------------------------------------------------------
        ! At this point sc_list contains
        !
        ! (1/N) sum_k (r_^{b}r^{c}_{a}+r_^{c}r^{b}_{a})(k) delta(w),
        !
        ! an approximation to
        !
        ! V_c.int dk/(2.pi)^3 (r_^{b}r^{c}_{a}+r_^{c}r^{b}_{a})(k) delta(w) dk
        !
        ! (V_c is the cell volume). We want
        !
        ! sigma_{abc}=( pi.e^3/(4.hbar^2) ) int dk/(2.pi)^3 Im[ (r_^{b}r^{c}_{a}+r_^{c}r^{b}_{a})(k) delta(w) ] dk
        !
        ! Note factor 1/4 instead of 1/2 as compared to SS PRB 61 5337 (2000) (Eq. 57),
        ! because we introduce 2 delta functions instead of 1.
        ! Hence we need to multiply by  pi.e^3/(4.hbar^2.V_c).
        ! To get the nonlinear response in units of A/V^2,
        !
        ! (i)   Divide by V_c to obtain (1/N) sum_k (r_^{b}r^{c}_{a}+r_^{c}r^{b}_{a})delta(w)/V_c, with units
        !       of [T] (integrand terms r_^{b}r^{c}_{a} delta(w) have units of [T].[L]^3)
        ! (ii)  Multiply by eV_seconds to convert the units of [T] from eV to seconds (coming from delta function)
        ! (iii) Multiply by ( pi.e^3/(4.hbar^2) ) in SI, which multiplied by [T] in seconds from (ii), gives final
        !       units of A/V^2
        !
        !==================================================
        ! fac = eV_seconds.( pi.e^3/(4.hbar^2.V_c) )
        !==================================================
        !
        ! with 'V_c' in Angstroms^3, and 'e', 'hbar' in SI units
        ! --------------------------------------------------------------------

        fac = physics%eV_seconds*pi*physics%elem_charge_SI**3/(4*physics%hbar_SI**(2)*cell_volume)
        write (stdout, '(/,1x,a)') &
          '----------------------------------------------------------'
        write (stdout, '(1x,a)') &
          'Output data files related to shift current:               '
        write (stdout, '(1x,a)') &
          '----------------------------------------------------------'

        do i = 1, 3
          do jk = 1, 6
            j = alpha_S(jk)
            k = beta_S(jk)
            file_name = trim(seedname)//'-sc_'// &
                        achar(119 + i)//achar(119 + j)//achar(119 + k)//'.dat'
            file_name = trim(file_name)
            write (stdout, '(/,3x,a)') '* '//file_name
            open (newunit=file_unit, FILE=file_name, STATUS='UNKNOWN', FORM='FORMATTED')
            do ifreq = 1, pw90_berry%kubo_nfreq
              write (file_unit, '(2E18.8E3)') real(pw90_berry%kubo_freq_list(ifreq), dp), &
                fac*sc_list(i, jk, ifreq)
            end do
            close (file_unit)
          end do
        end do

      end if

      ! -----------------------!
      ! Spin Hall conductivity !
      ! -----------------------!
      !
      if (eval_shc) then
        !
        ! Convert to the unit: (hbar/e) S/cm
        ! at this point, we need to
        ! (i)   multiply -e^2/hbar/(V*N_k) as in the QZYZ18 Eq.(5),
        !       note 1/N_k has already been applied by the kweight
        ! (ii)  convert charge current to spin current:
        !       divide the result by -e and multiply hbar/2 to
        !       recover the spin current, so the overall
        !       effect is -hbar/2/e
        ! (iii) multiply 1e8 to convert it to the unit S/cm
        ! So, the overall factor is
        !   fac = 1.0e8 * e^2 / hbar / V / 2.0
        ! and the final unit of spin Hall conductivity is (hbar/e)S/cm
        !
        fac = 1.0e8_dp*physics%elem_charge_SI**2/(physics%hbar_SI*cell_volume)/2.0_dp
        if (pw90_spin_hall%freq_scan) then
          shc_freq = shc_freq*fac
        else
          shc_fermi = shc_fermi*fac
        end if
        !
        write (stdout, '(/,1x,a)') &
          '----------------------------------------------------------'
        write (stdout, '(1x,a)') &
          'Output data files related to Spin Hall conductivity:'
        write (stdout, '(1x,a)') &
          '----------------------------------------------------------'
        !
        if (.not. pw90_spin_hall%freq_scan) then
          file_name = trim(seedname)//'-shc-fermiscan'//'.dat'
        else
          file_name = trim(seedname)//'-shc-freqscan'//'.dat'
        end if
        file_name = trim(file_name)
        write (stdout, '(/,3x,a)') '* '//file_name
        open (newunit=file_unit, FILE=file_name, STATUS='UNKNOWN', FORM='FORMATTED')
        if (.not. pw90_spin_hall%freq_scan) then
          write (file_unit, '(a,3x,a,3x,a)') &
            '#No.', 'Fermi energy(eV)', 'SHC((hbar/e)*S/cm)'
          do n = 1, fermi_n
            write (file_unit, '(I4,1x,F12.6,1x,E17.8)') &
              n, fermi_energy_list(n), shc_fermi(n)
          end do
        else
          write (file_unit, '(a,3x,a,3x,a,3x,a)') '#No.', 'Frequency(eV)', &
            'Re(sigma)((hbar/e)*S/cm)', 'Im(sigma)((hbar/e)*S/cm)'
          do n = 1, pw90_berry%kubo_nfreq
            write (file_unit, '(I4,1x,F12.6,1x,1x,2(E17.8,1x))') n, &
              real(pw90_berry%kubo_freq_list(n), dp), real(shc_freq(n), dp), aimag(shc_freq(n))
          end do
        end if
        close (file_unit)

      end if

      if (eval_kdotp) then
        ! -----------------------------!
        ! k.p expansion coefficients
        ! -----------------------------!

        write (stdout, '(/,1x,a)') &
          '----------------------------------------------------------'
        write (stdout, '(1x,a)') &
          'Output data files related to k.p:                         '
        write (stdout, '(1x,a)') &
          '----------------------------------------------------------'
        ! zeroth order in k
        file_name = trim(seedname)//'-kdotp_0.dat'
        file_name = trim(file_name)
        write (stdout, '(/,3x,a)') '* '//file_name
        open (newunit=file_unit, FILE=file_name, STATUS='UNKNOWN', FORM='FORMATTED')
        write (file_unit, '(2E18.8E3)') kdotp(:, :, 1, 1, 1)
        close (file_unit)

        ! first order in k
        file_name = trim(seedname)//'-kdotp_1.dat'
        write (stdout, '(/,3x,a)') '* '//file_name
        open (newunit=file_unit, FILE=file_name, STATUS='UNKNOWN', FORM='FORMATTED')
        do i = 1, 3
          write (file_unit, '(2E18.8E3)') kdotp(:, :, 2, i, 1)
        end do
        close (file_unit)

        ! second order in k
        file_name = trim(seedname)//'-kdotp_2.dat'
        write (stdout, '(/,3x,a)') '* '//file_name
        open (newunit=file_unit, FILE=file_name, STATUS='UNKNOWN', FORM='FORMATTED')
        do i = 1, 3
          do j = 1, 3
            write (file_unit, '(2E18.8E3)') kdotp(:, :, 3, i, j)
          end do
        end do
        close (file_unit)

      end if

    end if !print_output%iprint >0, aka "on_root"

  end subroutine berry_main