gyrotropic_main Subroutine

public subroutine gyrotropic_main(pw90_berry, dis_manifold, fermi_energy_list, pw90_gyrotropic, kmesh_info, kpt_latt, physics, pw90_oper_read, pw90_band_deriv_degen, ws_region, w90_system, print_output, wannier_data, wigner_seitz, ws_distance, AA_R, BB_R, CC_R, HH_R, SS_R, u_matrix, v_matrix, eigval, real_lattice, scissors_shift, mp_grid, num_bands, num_kpts, num_wann, effective_model, have_disentangled, seedname, stdout, timer, error, comm)

Uses

  • proc~~gyrotropic_main~~UsesGraph proc~gyrotropic_main gyrotropic_main module~w90_comms w90_comms proc~gyrotropic_main->module~w90_comms module~w90_constants w90_constants proc~gyrotropic_main->module~w90_constants module~w90_get_oper w90_get_oper proc~gyrotropic_main->module~w90_get_oper module~w90_io w90_io proc~gyrotropic_main->module~w90_io module~w90_postw90_types w90_postw90_types proc~gyrotropic_main->module~w90_postw90_types module~w90_types w90_types proc~gyrotropic_main->module~w90_types module~w90_utility w90_utility proc~gyrotropic_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_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) D tensor (ii) K tensor (iii) C tensor (iv) current-induced optical activity (v) natural optical activity

Arguments

Type IntentOptional Attributes Name
type(pw90_berry_mod_type), intent(in) :: pw90_berry
type(dis_manifold_type), intent(in) :: dis_manifold
real(kind=dp), intent(in), allocatable :: fermi_energy_list(:)
type(pw90_gyrotropic_type), intent(in) :: pw90_gyrotropic
type(kmesh_info_type), intent(in) :: kmesh_info
real(kind=dp), intent(in) :: kpt_latt(:,:)
type(pw90_physical_constants_type), intent(in) :: physics
type(pw90_oper_read_type), intent(in) :: pw90_oper_read
type(pw90_band_deriv_degen_type), intent(in) :: pw90_band_deriv_degen
type(ws_region_type), intent(in) :: ws_region
type(w90_system_type), intent(in) :: w90_system
type(print_output_type), intent(in) :: print_output
type(wannier_data_type), intent(in) :: wannier_data
type(wigner_seitz_type), intent(inout) :: wigner_seitz
type(ws_distance_type), intent(inout) :: ws_distance
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 :: SS_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) :: num_bands
integer, intent(in) :: num_kpts
integer, intent(in) :: num_wann
logical, intent(in) :: effective_model
logical, intent(in) :: have_disentangled
character(len=50), intent(in) :: seedname
integer, intent(in) :: stdout
type(timer_list_type), intent(inout) :: timer
type(w90_error_type), intent(out), allocatable :: error
type(w90_comm_type), intent(in) :: comm

Calls

proc~~gyrotropic_main~~CallsGraph proc~gyrotropic_main gyrotropic_main interface~comms_reduce comms_reduce proc~gyrotropic_main->interface~comms_reduce proc~get_aa_r get_AA_R proc~gyrotropic_main->proc~get_aa_r proc~get_aa_r_effective get_AA_R_effective proc~gyrotropic_main->proc~get_aa_r_effective proc~get_bb_r get_BB_R proc~gyrotropic_main->proc~get_bb_r proc~get_cc_r get_CC_R proc~gyrotropic_main->proc~get_cc_r proc~get_hh_r get_HH_R proc~gyrotropic_main->proc~get_hh_r proc~get_ss_r get_SS_R proc~gyrotropic_main->proc~get_ss_r proc~gyrotropic_get_k_list gyrotropic_get_k_list proc~gyrotropic_main->proc~gyrotropic_get_k_list proc~gyrotropic_outprint_tensor gyrotropic_outprint_tensor proc~gyrotropic_main->proc~gyrotropic_outprint_tensor proc~io_stopwatch_start io_stopwatch_start proc~gyrotropic_main->proc~io_stopwatch_start proc~io_stopwatch_stop io_stopwatch_stop proc~gyrotropic_main->proc~io_stopwatch_stop proc~mpirank mpirank proc~gyrotropic_main->proc~mpirank proc~mpisize mpisize proc~gyrotropic_main->proc~mpisize proc~set_error_input set_error_input proc~gyrotropic_main->proc~set_error_input proc~utility_det3 utility_det3 proc~gyrotropic_main->proc~utility_det3 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~get_aa_r->interface~comms_reduce 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~comms_array_split comms_array_split proc~get_aa_r->proc~comms_array_split 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~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~comms_array_split 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~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~comms_array_split 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_ss_r->proc~io_stopwatch_start proc~get_ss_r->proc~io_stopwatch_stop proc~get_ss_r->proc~mpirank 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~set_error_alloc set_error_alloc proc~get_ss_r->proc~set_error_alloc proc~set_error_dealloc set_error_dealloc 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~gyrotropic_get_k_list->proc~set_error_input proc~berry_get_imf_klist berry_get_imf_klist proc~gyrotropic_get_k_list->proc~berry_get_imf_klist proc~berry_get_imfgh_klist berry_get_imfgh_klist proc~gyrotropic_get_k_list->proc~berry_get_imfgh_klist proc~gyrotropic_get_curv_w_k gyrotropic_get_curv_w_k proc~gyrotropic_get_k_list->proc~gyrotropic_get_curv_w_k proc~gyrotropic_get_noa_k gyrotropic_get_NOA_k proc~gyrotropic_get_k_list->proc~gyrotropic_get_noa_k proc~pw90common_fourier_r_to_k_vec pw90common_fourier_R_to_k_vec proc~gyrotropic_get_k_list->proc~pw90common_fourier_r_to_k_vec proc~spin_get_s spin_get_S proc~gyrotropic_get_k_list->proc~spin_get_s proc~utility_rotate utility_rotate proc~gyrotropic_get_k_list->proc~utility_rotate proc~utility_w0gauss utility_w0gauss proc~gyrotropic_get_k_list->proc~utility_w0gauss proc~wham_get_d_h wham_get_D_h proc~gyrotropic_get_k_list->proc~wham_get_d_h proc~wham_get_eig_deleig wham_get_eig_deleig proc~gyrotropic_get_k_list->proc~wham_get_eig_deleig proc~gyrotropic_outprint_tensor_w gyrotropic_outprint_tensor_w proc~gyrotropic_outprint_tensor->proc~gyrotropic_outprint_tensor_w proc~comms_sync_error comms_sync_error proc~set_error_input->proc~comms_sync_error proc~set_base_error set_base_error proc~set_error_input->proc~set_base_error proc~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~berry_get_imf_klist->proc~berry_get_imfgh_klist proc~berry_get_imfgh_klist->proc~pw90common_fourier_r_to_k_vec proc~pw90common_fourier_r_to_k pw90common_fourier_R_to_k proc~berry_get_imfgh_klist->proc~pw90common_fourier_r_to_k 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~comms_array_split->proc~mpisize 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~gyrotropic_get_noa_k->proc~utility_rotate proc~gyrotropic_get_noa_bnl_orb gyrotropic_get_NOA_Bnl_orb proc~gyrotropic_get_noa_k->proc~gyrotropic_get_noa_bnl_orb proc~gyrotropic_get_noa_bnl_spin gyrotropic_get_NOA_Bnl_spin proc~gyrotropic_get_noa_k->proc~gyrotropic_get_noa_bnl_spin proc~pw90common_fourier_r_to_k_new pw90common_fourier_R_to_k_new proc~gyrotropic_get_noa_k->proc~pw90common_fourier_r_to_k_new proc~set_error_alloc->proc~comms_sync_error proc~set_error_alloc->proc~set_base_error proc~set_error_dealloc->proc~comms_sync_error proc~set_error_dealloc->proc~set_base_error proc~set_error_fatal->proc~comms_sync_error proc~set_error_fatal->proc~set_base_error proc~set_error_file->proc~comms_sync_error proc~set_error_file->proc~set_base_error proc~spin_get_s->proc~pw90common_fourier_r_to_k proc~utility_diagonalize utility_diagonalize proc~spin_get_s->proc~utility_diagonalize proc~utility_rotate_diag utility_rotate_diag proc~spin_get_s->proc~utility_rotate_diag proc~utility_w0gauss->proc~set_error_input proc~wham_get_d_h->proc~utility_rotate proc~wham_get_eig_deleig->proc~get_hh_r proc~wham_get_eig_deleig->proc~pw90common_fourier_r_to_k proc~wham_get_eig_deleig->proc~utility_diagonalize proc~wham_get_deleig_a wham_get_deleig_a proc~wham_get_eig_deleig->proc~wham_get_deleig_a proc~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~utility_diagonalize->proc~set_error_fatal zhpevx zhpevx proc~utility_diagonalize->zhpevx proc~utility_rotate_diag->proc~utility_zgemm_new proc~utility_matmul_diag utility_matmul_diag proc~utility_rotate_diag->proc~utility_matmul_diag zgemm zgemm proc~utility_zgemm_new->zgemm proc~utility_zgemmm->proc~utility_zgemm_new proc~wham_get_deleig_a->proc~utility_rotate proc~wham_get_deleig_a->proc~utility_diagonalize proc~wham_get_deleig_a->proc~utility_rotate_diag proc~wham_get_eig_uu_hh_jjlist->proc~get_hh_r proc~wham_get_eig_uu_hh_jjlist->proc~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~pw90common_get_occ pw90common_get_occ proc~wham_get_occ_mat_list->proc~pw90common_get_occ 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~~gyrotropic_main~~CalledByGraph proc~gyrotropic_main gyrotropic_main program~postw90 postw90 program~postw90->proc~gyrotropic_main

Source Code

  subroutine gyrotropic_main(pw90_berry, dis_manifold, fermi_energy_list, pw90_gyrotropic, &
                             kmesh_info, kpt_latt, physics, pw90_oper_read, pw90_band_deriv_degen, &
                             ws_region, w90_system, print_output, wannier_data, wigner_seitz, &
                             ws_distance, AA_R, BB_R, CC_R, HH_R, SS_R, u_matrix, v_matrix, &
                             eigval, real_lattice, scissors_shift, mp_grid, num_bands, num_kpts, &
                             num_wann, effective_model, have_disentangled, seedname, stdout, &
                             timer, error, comm)
    !================================================!
    !
    !! Computes the following quantities:
    !!   (i) D tensor
    !!  (ii) K tensor
    !! (iii) C tensor
    !!  (iv) current-induced optical activity
    !!   (v) natural optical activity
    !
    !================================================!

    use w90_comms, only: comms_reduce, w90_comm_type, mpirank, mpisize
    use w90_constants, only: dp, twopi, pw90_physical_constants_type
    use w90_get_oper, only: get_HH_R, get_AA_R_effective, get_AA_R, get_BB_R, get_CC_R, get_SS_R
    use w90_io, only: io_stopwatch_start, io_stopwatch_stop
    use w90_postw90_types, only: pw90_gyrotropic_type, pw90_berry_mod_type, pw90_oper_read_type, &
                                 pw90_band_deriv_degen_type, wigner_seitz_type
    use w90_types, only: dis_manifold_type, print_output_type, timer_list_type, &
                         kmesh_info_type, wannier_data_type, ws_region_type, w90_system_type, ws_distance_type
    use w90_utility, only: utility_det3

    implicit none

    ! arguments
    type(pw90_berry_mod_type), intent(in) :: pw90_berry
    type(dis_manifold_type), intent(in) :: dis_manifold
    type(pw90_gyrotropic_type), intent(in) :: pw90_gyrotropic
    type(kmesh_info_type), intent(in) :: kmesh_info
    type(pw90_band_deriv_degen_type), intent(in) :: pw90_band_deriv_degen
    type(pw90_oper_read_type), intent(in) :: pw90_oper_read
    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(w90_comm_type), intent(in) :: comm
    type(w90_system_type), intent(in) :: w90_system
    type(wannier_data_type), intent(in) :: wannier_data
    type(wigner_seitz_type), intent(inout) :: wigner_seitz
    type(ws_distance_type), intent(inout) :: ws_distance
    type(timer_list_type), intent(inout) :: timer
    type(w90_error_type), allocatable, intent(out) :: error

    complex(kind=dp), allocatable, intent(inout) :: AA_R(:, :, :, :)
    complex(kind=dp), allocatable, intent(inout) :: BB_R(:, :, :, :)
    complex(kind=dp), allocatable, intent(inout) :: CC_R(:, :, :, :, :)
    complex(kind=dp), allocatable, intent(inout) :: HH_R(:, :, :)
    complex(kind=dp), allocatable, intent(inout) :: SS_R(:, :, :, :)
    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_bands, num_kpts, num_wann
    integer, intent(in) :: stdout

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

    ! local variables
    real(kind=dp), allocatable    :: gyro_K_spn(:, :, :)
    real(kind=dp), allocatable    :: gyro_DOS(:)
    real(kind=dp), allocatable    :: gyro_K_orb(:, :, :)
    real(kind=dp), allocatable    :: gyro_C(:, :, :)
    real(kind=dp), allocatable    :: gyro_D(:, :, :)
    real(kind=dp), allocatable    :: gyro_Dw(:, :, :, :)
    real(kind=dp), allocatable    :: gyro_NOA_spn(:, :, :, :)
    real(kind=dp), allocatable    :: gyro_NOA_orb(:, :, :, :)

    character(len=30) :: f_out_name_tmp
    character(len=30) :: units_tmp
    character(len=120) :: comment_tmp

    real(kind=dp) :: cell_volume
    real(kind=dp)     :: kweight, kpt(3), &
                         db1, db2, db3, fac
    integer           :: loop_x, loop_y, loop_z, loop_xyz
    integer           :: fermi_n
    logical           :: eval_K, eval_C, eval_D, eval_Dw, eval_NOA, eval_spn, eval_DOS

    integer :: my_node_id, num_nodes
    my_node_id = mpirank(comm)
    num_nodes = mpisize(comm)

    if (.not. allocated(fermi_energy_list)) then
      call set_error_input(error, 'Must specify one or more Fermi levels when gyrotropic=true', comm)
      return
    end if

    if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
      call io_stopwatch_start('gyrotropic: 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_gyrotropic%kmesh%mesh(1), dp)
    db2 = 1.0_dp/real(pw90_gyrotropic%kmesh%mesh(2), dp)
    db3 = 1.0_dp/real(pw90_gyrotropic%kmesh%mesh(3), dp)

    eval_K = .false.
    eval_C = .false.
    eval_D = .false.
    eval_Dw = .false.
    eval_spn = .false.
    eval_NOA = .false.
    eval_DOS = .false.

    if (index(pw90_gyrotropic%task, '-k') > 0) eval_K = .true.
    if (index(pw90_gyrotropic%task, '-c') > 0) eval_C = .true.
    if (index(pw90_gyrotropic%task, '-d0') > 0) eval_D = .true.
    if (index(pw90_gyrotropic%task, '-dw') > 0) eval_Dw = .true.
    if (index(pw90_gyrotropic%task, '-spin') > 0) eval_spn = .true.
    if (index(pw90_gyrotropic%task, '-noa') > 0) eval_NOA = .true.
    if (index(pw90_gyrotropic%task, '-dos') > 0) eval_DOS = .true.
    if (index(pw90_gyrotropic%task, 'all') > 0) then
      eval_K = .true.
      eval_C = .true.
      eval_D = .true.
      eval_Dw = .true.
      if (w90_system%spinors) eval_spn = .true.
      eval_NOA = .true.
      eval_DOS = .true.
    end if

    if (.not. (eval_K .or. eval_noa)) eval_spn = .false.

    if ((.not. w90_system%spinors) .and. eval_spn) then
      call set_error_input(error, "spin contribution requested for gyrotropic, but the wavefunctions are not spinors", comm)
      return
    end if

    ! Wannier matrix elements, allocations and initializations

    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, &
                  w90_system%num_valence_bands, effective_model, have_disentangled, seedname, &
                  ws_distance, ws_region, stdout, timer, error, comm)
    if (allocated(error)) return

    if (eval_D .or. eval_Dw .or. eval_K .or. eval_NOA) then

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

    end if

    if (eval_spn) then

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

    end if

    ! not allocated was tested at start of routine
    fermi_n = size(fermi_energy_list)
    if (eval_K) then
      call get_BB_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, HH_R, BB_R, v_matrix, &
                    eigval, scissors_shift, wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, &
                    num_wann, have_disentangled, seedname, stdout, timer, error, comm)
      if (allocated(error)) return

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

      allocate (gyro_K_orb(3, 3, fermi_n))
      gyro_K_orb = 0.0_dp
      if (eval_spn) then
        allocate (gyro_K_spn(3, 3, fermi_n))
        gyro_K_spn = 0.0_dp
      end if
    end if

    if (eval_D) then
      allocate (gyro_D(3, 3, fermi_n))
      gyro_D = 0.0_dp
    end if

    if (eval_DOS) then
      allocate (gyro_DOS(fermi_n))
      gyro_DOS = 0.0_dp
    end if

    if (eval_C) then
      allocate (gyro_C(3, 3, fermi_n))
      gyro_C = 0.0_dp
    end if

    if (eval_Dw) then
      allocate (gyro_Dw(3, 3, fermi_n, pw90_gyrotropic%nfreq))
      gyro_Dw = 0.0_dp
    end if

    if (eval_NOA) then
      allocate (gyro_NOA_orb(3, 3, fermi_n, pw90_gyrotropic%nfreq))
      gyro_NOA_orb = 0.0_dp
      if (eval_spn) then
        allocate (gyro_NOA_spn(3, 3, fermi_n, pw90_gyrotropic%nfreq))
        gyro_NOA_spn = 0.0_dp
      end if
    end if

    if (print_output%iprint > 0) then
      flush (stdout)
      write (stdout, '(/,/,1x,a)') 'Properties calculated in module  g y r o t r o p i c'
      write (stdout, '(1x,a)') '------------------------------------------'

      if (eval_D) write (stdout, '(/,3x,a)') '* D-tensor  --- Eq.2 of TAS17 '

      if (eval_dos) write (stdout, '(/,3x,a)') '* density of states '

      if (eval_K) then
        write (stdout, '(/,3x,a)') '* K-tensor  --- Eq.3 of TAS17 '
        if (eval_spn) then
          write (stdout, '(3x,a)') '    * including spin component '
        else
          write (stdout, '(3x,a)') '    * excluding spin component '
        end if
      end if

      if (eval_Dw) write (stdout, '(/,3x,a)') '* Dw-tensor  --- Eq.12 of TAS17 '

      if (eval_C) write (stdout, '(/,3x,a)') '* C-tensor  --- Eq.B6 of TAS17 '

      if (eval_NOA) then
        write (stdout, '(/,3x,a)') '* gamma-tensor of NOA --- Eq.C12 of TAS17 '
        if (eval_spn) then
          write (stdout, '(3x,a)') '    * including spin component '
        else
          write (stdout, '(3x,a)') '    * excluding spin component '
        end if
      end if

      if (pw90_berry%transl_inv) then
        if (eval_K) then
          call set_error_input(error, 'transl_inv=T disabled for K-tensor', 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 (print_output%timing_level > 1) then
        call io_stopwatch_stop('gyrotropic: prelims', timer)
        call io_stopwatch_start('gyrotropic: k-interpolation', timer)
      end if

      write (stdout, '(1x,a20,3(i0,1x))') 'Interpolation grid: ', pw90_gyrotropic%kmesh%mesh(1:3)

      flush (stdout)

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

    ! Do not read 'kpoint.dat'. Loop over a regular grid in the full BZ

    kweight = db1*db2*db3*utility_det3(pw90_gyrotropic%box)

    do loop_xyz = my_node_id, PRODUCT(pw90_gyrotropic%kmesh%mesh) - 1, num_nodes
      loop_x = loop_xyz/(pw90_gyrotropic%kmesh%mesh(2)*pw90_gyrotropic%kmesh%mesh(3))
      loop_y = (loop_xyz - loop_x*(pw90_gyrotropic%kmesh%mesh(2) &
                                   *pw90_gyrotropic%kmesh%mesh(3)))/pw90_gyrotropic%kmesh%mesh(3)
      loop_z = loop_xyz - loop_x*(pw90_gyrotropic%kmesh%mesh(2)*pw90_gyrotropic%kmesh%mesh(3)) &
               - loop_y*pw90_gyrotropic%kmesh%mesh(3)
      kpt(1) = loop_x*db1
      kpt(2) = loop_y*db2
      kpt(3) = loop_z*db3
      kpt(:) = pw90_gyrotropic%box_corner(:) + matmul(kpt, pw90_gyrotropic%box)

      call gyrotropic_get_k_list(ws_region, w90_system%num_valence_bands, have_disentangled, kpt, &
                                 kweight, gyro_K_spn, gyro_K_orb, gyro_D, gyro_Dw, gyro_C, &
                                 gyro_DOS, gyro_NOA_orb, gyro_NOA_spn, eval_K, eval_D, eval_Dw, &
                                 eval_NOA, eval_spn, eval_C, eval_dos, num_wann, print_output, &
                                 fermi_energy_list, wannier_data, eigval, real_lattice, mp_grid, &
                                 num_bands, num_kpts, u_matrix, v_matrix, dis_manifold, kpt_latt, &
                                 pw90_gyrotropic, scissors_shift, effective_model, &
                                 pw90_band_deriv_degen, ws_distance, wigner_seitz, stdout, &
                                 seedname, timer, error, comm, HH_R, AA_R, BB_R, CC_R, SS_R)
      if (allocated(error)) return

    end do !loop_xyz

    ! Collect contributions from all nodes
    if (eval_K) then
      call comms_reduce(gyro_K_orb(1, 1, 1), 3*3*fermi_n, 'SUM', error, comm)
      if (allocated(error)) return
      if (eval_spn) then
        call comms_reduce(gyro_K_spn(1, 1, 1), 3*3*fermi_n, 'SUM', error, comm)
        if (allocated(error)) return
      end if
    end if

    if (eval_D) then
      call comms_reduce(gyro_D(1, 1, 1), 3*3*fermi_n, 'SUM', error, comm)
      if (allocated(error)) return
    end if

    if (eval_C) then
      call comms_reduce(gyro_C(1, 1, 1), 3*3*fermi_n, 'SUM', error, comm)
      if (allocated(error)) return
    end if

    if (eval_Dw) then
      call comms_reduce(gyro_Dw(1, 1, 1, 1), 3*3*fermi_n*pw90_gyrotropic%nfreq, 'SUM', error, comm)
      if (allocated(error)) return
    end if

    if (eval_dos) then
      call comms_reduce(gyro_DOS(1), fermi_n, 'SUM', error, comm)
      if (allocated(error)) return
    end if

    if (eval_NOA) then
      call comms_reduce(gyro_NOA_orb(1, 1, 1, 1), 3*3*fermi_n*pw90_gyrotropic%nfreq, &
                        'SUM', error, comm)
      if (allocated(error)) return
      if (eval_spn) then
        call comms_reduce(gyro_NOA_spn(1, 1, 1, 1), 3*3*fermi_n*pw90_gyrotropic%nfreq, &
                          '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('gyrotropic: k-interpolation', timer)
      write (stdout, '(1x,a)') ' '
      write (stdout, *) 'Calculation finished, writing results'
      flush (stdout)

      if (eval_K) then
        if (eval_spn) then
          ! At this point gme_spn_list contains
          ! (1/N) sum_k delta(E_kn-E_f).(d E_{kn}/d k_i).sigma_{kn,j}
          ! (units of length) in Angstroms.
          !        ====================================
          ! To get K in units of Ampere do the following:
          !        ====================================
          !   * Divide by V_c in Ang^3 to get a quantity with units of [L]^{-2}
          !   * Multiply by 10^20 to convert to SI
          !   * Multiply by -g_s.e.hbar/(4m_e) \simeq e.hbar/(2.m_e) in SI units
          !==================================================
          ! fac = 10^20*e*hbar/(2.m_e.V_c)
          !==================================================
          fac = -1.0e20_dp*physics%elem_charge_SI*physics%hbar_SI/(2.*physics%elec_mass_SI &
                                                                   *cell_volume)
          gyro_K_spn(:, :, :) = gyro_K_spn(:, :, :)*fac
          f_out_name_tmp = 'K_spin'
          units_tmp = "Ampere"
          comment_tmp = "spin part of the K tensor -- Eq. 3 of TAS17"
          call gyrotropic_outprint_tensor(stdout, seedname, pw90_gyrotropic, fermi_energy_list, f_out_name_tmp, &
                                          arrEf=gyro_K_spn, units=units_tmp, comment=comment_tmp)
        end if  ! eval_K && eval_spin

        ! At this point gme_orb_list contains
        ! (1/N)sum_{k,n} delta(E_kn-E_f).(d E_{kn}/d k_i)
        !                           .Im[<del_k u_kn| x (H_k-E_kn)|del_k u_kn>]
        ! (units of energy times length^3) in eV.Ang^3.
        !        ====================================
        ! To get K  in units of Ampere do the following:
        !        ====================================
        !   * Divide by V_c in Ang^3 to get a quantity with units of eV
        !   * Multiply by 'e' in SI to convert to SI (Joules)
        !   * Multiply by e/(2.hbar) to get K in Ampere
        !==================================================
        ! fac = e^2/(2.hbar.V_c)
        !==================================================
        fac = physics%elem_charge_SI**2/(2.*physics%hbar_SI*cell_volume)
        gyro_K_orb(:, :, :) = gyro_K_orb(:, :, :)*fac

        f_out_name_tmp = 'K_orb'
        units_tmp = "Ampere"
        comment_tmp = "orbital part of the K tensor -- Eq. 3 of TAS17"
        call gyrotropic_outprint_tensor(stdout, seedname, pw90_gyrotropic, fermi_energy_list, f_out_name_tmp, &
                                        arrEf=gyro_K_orb, units=units_tmp, comment=comment_tmp)
      end if ! eval_K

      if (eval_D) then
        fac = 1./cell_volume
        gyro_D(:, :, :) = gyro_D(:, :, :)*fac

        f_out_name_tmp = 'D'
        units_tmp = "dimensionless"
        comment_tmp = "the D tensor -- Eq. 2 of TAS17"
        call gyrotropic_outprint_tensor(stdout, seedname, pw90_gyrotropic, fermi_energy_list, f_out_name_tmp, &
                                        arrEf=gyro_D, units=units_tmp, comment=comment_tmp)
      end if

      if (eval_Dw) then
        fac = 1./cell_volume
        gyro_Dw(:, :, :, :) = gyro_Dw(:, :, :, :)*fac

        f_out_name_tmp = 'tildeD'
        units_tmp = "dimensionless"
        comment_tmp = "the tildeD tensor -- Eq. 12 of TAS17"
        call gyrotropic_outprint_tensor(stdout, seedname, pw90_gyrotropic, fermi_energy_list, f_out_name_tmp, &
                                        arrEfW=gyro_Dw, units=units_tmp, comment=comment_tmp)
      end if

      if (eval_C) then
        ! At this point gyro_C contains
        ! (1/N)sum_{k,n} delta(E_kn-E_f).(d E_{kn}/d k_i).(d E_{kn}/d k_j)
        ! (units of energy*length^2) in eV*Ang^2
        !
        ! To get it in Cab = e/h * (1/N*V_cell)sum_{k,n} delta(E_kn-E_f).(d E_{kn}/d k_i).(d E_{kn}/d k_j)
        ! in units Ampere/cm
        !
        ! divide by V_c in Ang^3 to get  eV/Ang
        ! multiply by 10^8*e in SI to get J/cm
        ! multiply by e/h in SI
        !
        fac = 1.0e+8_dp*physics%elem_charge_SI**2/(twopi*physics%hbar_SI*cell_volume)
        gyro_C(:, :, :) = gyro_C(:, :, :)*fac

        f_out_name_tmp = 'C'
        units_tmp = "Ampere/cm"
        comment_tmp = "the C tensor -- Eq. B6 of TAS17"
        call gyrotropic_outprint_tensor(stdout, seedname, pw90_gyrotropic, fermi_energy_list, f_out_name_tmp, &
                                        arrEf=gyro_C, units=units_tmp, comment=comment_tmp)
      end if

      if (eval_noa) then
        ! at this point gyro_NOA_orb  is in eV^-1.Ang^3   !
        !  We want the result in angstrems  !
        !   * Divide by V_c in Ang^3 to make it eV^{-1}
        !   * Divide by e in SI to get J^{-1}
        !   * multiply by e^2/eps_0 to get meters
        !   *multiply dy 1e10 to get Ang
        fac = 1e+10_dp*physics%elem_charge_SI/(cell_volume*physics%eps0_SI)
        gyro_NOA_orb = gyro_NOA_orb*fac
        f_out_name_tmp = 'NOA_orb'
        units_tmp = "Ang"
        comment_tmp = "the tensor $gamma_{abc}^{orb}$ (Eq. C12,C14 of TAS17)"
        call gyrotropic_outprint_tensor(stdout, seedname, pw90_gyrotropic, fermi_energy_list, f_out_name_tmp, &
                                        arrEfW=gyro_NOA_orb, units=units_tmp, comment=comment_tmp, &
                                        symmetrize=.false.)

        if (eval_spn) then
          ! at this point gyro_NOA_spn  is in eV^-2.Ang   !
          !  We want the result in angstrems  !
          !   * Divide by V_c in Ang^3 to make it (eV.Ang)^{-2}
          !   * multiply by 1e20/e^2 in SI to get (J.m)^{-2}
          !   * multiply by e^2/eps_0 to get (J.m)^{-1}
          !   *multiply dy hbar^2/m_e to get m
          !   *multiply by 1e10 to get Ang
          fac = 1e+30_dp*physics%hbar_SI**2/(cell_volume*physics%eps0_SI*physics%elec_mass_SI)
          gyro_NOA_spn = gyro_NOA_spn*fac
          f_out_name_tmp = 'NOA_spin'
          units_tmp = "Ang"
          comment_tmp = "the tensor $gamma_{abc}^{spin}$ (Eq. C12,C15 of TAS17)"
          call gyrotropic_outprint_tensor(stdout, seedname, pw90_gyrotropic, fermi_energy_list, f_out_name_tmp, &
                                          arrEfW=gyro_NOA_spn, units=units_tmp, &
                                          comment=comment_tmp, symmetrize=.false.)
        end if
      end if  !eval_NOA

      if (eval_DOS) then
        ! At this point gyro_C contains
        ! (1/N)sum_{k,n} delta(E_kn-E_f)
        ! in units of eV^{-1}
        ! divide by V_c in Ang^3 to get units 1./(eV*Ang^3)
        gyro_DOS(:) = gyro_DOS(:)/cell_volume
        f_out_name_tmp = 'DOS'
        units_tmp = "eV^{-1}.Ang^{-3}"
        comment_tmp = "density of states"
        call gyrotropic_outprint_tensor(stdout, seedname, pw90_gyrotropic, fermi_energy_list, f_out_name_tmp, &
                                        arrEf1d=gyro_DOS, units=units_tmp, comment=comment_tmp)
      end if

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

  end subroutine gyrotropic_main