boltzwann_main Subroutine

public subroutine boltzwann_main(pw90_boltzwann, dis_manifold, pw90_dos, kpt_latt, pw90_band_deriv_degen, postw90_oper, pw90_spin, physics, ws_region, w90_system, wannier_data, ws_distance, wigner_seitz, print_output, HH_R, SS_R, v_matrix, u_matrix, eigval, real_lattice, scissors_shift, mp_grid, num_wann, num_bands, num_kpts, effective_model, have_disentangled, spin_decomp, seedname, stdout, timer, error, comm)

Uses

  • proc~~boltzwann_main~~UsesGraph proc~boltzwann_main boltzwann_main module~w90_comms w90_comms proc~boltzwann_main->module~w90_comms module~w90_constants w90_constants proc~boltzwann_main->module~w90_constants module~w90_io w90_io proc~boltzwann_main->module~w90_io module~w90_postw90_types w90_postw90_types proc~boltzwann_main->module~w90_postw90_types module~w90_types w90_types proc~boltzwann_main->module~w90_types module~w90_comms->module~w90_constants module~w90_error_base w90_error_base module~w90_comms->module~w90_error_base 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

This is the main routine of the BoltzWann module. It calculates the transport coefficients using the Boltzmann transport equation.

It produces six files that contain:

  1. the Transport Distribution function (TDF) in units of 1/hbar^2 * eV*fs/angstrom
  2. the electrical conductivity in SI units (1/Ohm/m)
  3. the tensor sigma*S (sigma=el.cond., S=seebeck) in SI units (Ampere/meter/K)
  4. the Seebeck coefficient in SI units (V/K)
  5. the thermal conductivity in SI units (W/meter/K)
  6. if requested, the density of states

Files from 2 to 4 are output on a grid of (mu,T) points, where mu is the chemical potential in eV and T is the temperature in Kelvin. The grid is defined in the input.

Arguments

Type IntentOptional Attributes Name
type(pw90_boltzwann_type), intent(in) :: pw90_boltzwann
type(dis_manifold_type), intent(in) :: dis_manifold
type(pw90_dos_mod_type), intent(in) :: pw90_dos
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) :: postw90_oper
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(w90_system_type), intent(in) :: w90_system
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 :: HH_R(:,:,:)
complex(kind=dp), intent(inout), allocatable :: SS_R(:,:,:,:)
complex(kind=dp), intent(in) :: v_matrix(:,:,:)
complex(kind=dp), intent(in) :: u_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_wann
integer, intent(in) :: num_bands
integer, intent(in) :: num_kpts
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~~boltzwann_main~~CallsGraph proc~boltzwann_main boltzwann_main interface~comms_gatherv comms_gatherv proc~boltzwann_main->interface~comms_gatherv interface~comms_reduce comms_reduce proc~boltzwann_main->interface~comms_reduce proc~calctdfanddos calcTDFandDOS proc~boltzwann_main->proc~calctdfanddos proc~comms_array_split comms_array_split proc~boltzwann_main->proc~comms_array_split proc~io_stopwatch_start io_stopwatch_start proc~boltzwann_main->proc~io_stopwatch_start proc~io_stopwatch_stop io_stopwatch_stop proc~boltzwann_main->proc~io_stopwatch_stop proc~minusfermiderivative MinusFermiDerivative proc~boltzwann_main->proc~minusfermiderivative proc~mpirank mpirank proc~boltzwann_main->proc~mpirank proc~mpisize mpisize proc~boltzwann_main->proc~mpisize proc~set_error_alloc set_error_alloc proc~boltzwann_main->proc~set_error_alloc proc~set_error_dealloc set_error_dealloc proc~boltzwann_main->proc~set_error_dealloc proc~set_error_input set_error_input proc~boltzwann_main->proc~set_error_input proc~utility_inv2 utility_inv2 proc~boltzwann_main->proc~utility_inv2 proc~utility_inv3 utility_inv3 proc~boltzwann_main->proc~utility_inv3 proc~comms_gatherv_cmplx_1 comms_gatherv_cmplx_1 interface~comms_gatherv->proc~comms_gatherv_cmplx_1 proc~comms_gatherv_cmplx_2 comms_gatherv_cmplx_2 interface~comms_gatherv->proc~comms_gatherv_cmplx_2 proc~comms_gatherv_cmplx_3 comms_gatherv_cmplx_3 interface~comms_gatherv->proc~comms_gatherv_cmplx_3 proc~comms_gatherv_cmplx_3_4 comms_gatherv_cmplx_3_4 interface~comms_gatherv->proc~comms_gatherv_cmplx_3_4 proc~comms_gatherv_cmplx_4 comms_gatherv_cmplx_4 interface~comms_gatherv->proc~comms_gatherv_cmplx_4 proc~comms_gatherv_logical comms_gatherv_logical interface~comms_gatherv->proc~comms_gatherv_logical proc~comms_gatherv_real_1 comms_gatherv_real_1 interface~comms_gatherv->proc~comms_gatherv_real_1 proc~comms_gatherv_real_2 comms_gatherv_real_2 interface~comms_gatherv->proc~comms_gatherv_real_2 proc~comms_gatherv_real_2_3 comms_gatherv_real_2_3 interface~comms_gatherv->proc~comms_gatherv_real_2_3 proc~comms_gatherv_real_3 comms_gatherv_real_3 interface~comms_gatherv->proc~comms_gatherv_real_3 proc~comms_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~calctdfanddos->interface~comms_reduce proc~calctdfanddos->proc~io_stopwatch_start proc~calctdfanddos->proc~io_stopwatch_stop proc~calctdfanddos->proc~mpirank proc~calctdfanddos->proc~mpisize proc~calctdfanddos->proc~set_error_alloc proc~calctdfanddos->proc~set_error_dealloc proc~calctdfanddos->proc~set_error_input interface~comms_allreduce comms_allreduce proc~calctdfanddos->interface~comms_allreduce proc~dos_get_k dos_get_k proc~calctdfanddos->proc~dos_get_k proc~dos_get_levelspacing dos_get_levelspacing proc~calctdfanddos->proc~dos_get_levelspacing proc~get_hh_r get_HH_R proc~calctdfanddos->proc~get_hh_r proc~get_ss_r get_SS_R proc~calctdfanddos->proc~get_ss_r proc~tdf_kpt TDF_kpt proc~calctdfanddos->proc~tdf_kpt proc~utility_recip_lattice_base utility_recip_lattice_base proc~calctdfanddos->proc~utility_recip_lattice_base proc~w90_readwrite_get_smearing_type w90_readwrite_get_smearing_type proc~calctdfanddos->proc~w90_readwrite_get_smearing_type proc~wham_get_eig_deleig wham_get_eig_deleig proc~calctdfanddos->proc~wham_get_eig_deleig proc~comms_array_split->proc~mpisize 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_dealloc->proc~comms_sync_error proc~set_error_dealloc->proc~set_base_error proc~set_error_input->proc~comms_sync_error proc~set_error_input->proc~set_base_error proc~comms_allreduce_cmplx comms_allreduce_cmplx interface~comms_allreduce->proc~comms_allreduce_cmplx proc~comms_allreduce_real comms_allreduce_real interface~comms_allreduce->proc~comms_allreduce_real proc~comms_gatherv_cmplx_1->proc~comms_sync_error proc~comms_no_sync_gatherv_cmplx_1 comms_no_sync_gatherv_cmplx_1 proc~comms_gatherv_cmplx_1->proc~comms_no_sync_gatherv_cmplx_1 proc~comms_gatherv_cmplx_2->proc~comms_sync_error proc~comms_no_sync_gatherv_cmplx_2 comms_no_sync_gatherv_cmplx_2 proc~comms_gatherv_cmplx_2->proc~comms_no_sync_gatherv_cmplx_2 proc~comms_gatherv_cmplx_3->proc~comms_sync_error proc~comms_no_sync_gatherv_cmplx_3 comms_no_sync_gatherv_cmplx_3 proc~comms_gatherv_cmplx_3->proc~comms_no_sync_gatherv_cmplx_3 proc~comms_gatherv_cmplx_3_4->proc~comms_sync_error proc~comms_no_sync_gatherv_cmplx_3_4 comms_no_sync_gatherv_cmplx_3_4 proc~comms_gatherv_cmplx_3_4->proc~comms_no_sync_gatherv_cmplx_3_4 proc~comms_gatherv_cmplx_4->proc~comms_sync_error proc~comms_no_sync_gatherv_cmplx_4 comms_no_sync_gatherv_cmplx_4 proc~comms_gatherv_cmplx_4->proc~comms_no_sync_gatherv_cmplx_4 proc~comms_gatherv_logical->proc~comms_sync_error proc~comms_no_sync_gatherv_logical comms_no_sync_gatherv_logical proc~comms_gatherv_logical->proc~comms_no_sync_gatherv_logical proc~comms_gatherv_real_1->proc~comms_sync_error proc~comms_no_sync_gatherv_real_1 comms_no_sync_gatherv_real_1 proc~comms_gatherv_real_1->proc~comms_no_sync_gatherv_real_1 proc~comms_gatherv_real_2->proc~comms_sync_error proc~comms_no_sync_gatherv_real_2 comms_no_sync_gatherv_real_2 proc~comms_gatherv_real_2->proc~comms_no_sync_gatherv_real_2 proc~comms_gatherv_real_2_3->proc~comms_sync_error proc~comms_no_sync_gatherv_real_2_3 comms_no_sync_gatherv_real_2_3 proc~comms_gatherv_real_2_3->proc~comms_no_sync_gatherv_real_2_3 proc~comms_gatherv_real_3->proc~comms_sync_error proc~comms_no_sync_gatherv_real_3 comms_no_sync_gatherv_real_3 proc~comms_gatherv_real_3->proc~comms_no_sync_gatherv_real_3 proc~comms_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~dos_get_k->proc~set_error_input proc~spin_get_nk spin_get_nk proc~dos_get_k->proc~spin_get_nk proc~utility_w0gauss utility_w0gauss proc~dos_get_k->proc~utility_w0gauss interface~pw90common_kmesh_spacing pw90common_kmesh_spacing proc~dos_get_levelspacing->interface~pw90common_kmesh_spacing 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 interface~comms_bcast comms_bcast 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_win_min get_win_min proc~get_hh_r->proc~get_win_min proc~operator_wigner_setup operator_wigner_setup proc~get_hh_r->proc~operator_wigner_setup proc~set_error_fatal set_error_fatal proc~get_hh_r->proc~set_error_fatal proc~set_error_file set_error_file 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->proc~set_error_alloc proc~get_ss_r->proc~set_error_dealloc proc~get_ss_r->interface~comms_bcast proc~get_ss_r->proc~fourier_q_to_r proc~get_gauge_overlap_matrix get_gauge_overlap_matrix proc~get_ss_r->proc~get_gauge_overlap_matrix proc~get_ss_r->proc~operator_wigner_setup proc~get_ss_r->proc~set_error_fatal proc~get_ss_r->proc~set_error_file proc~tdf_kpt->proc~spin_get_nk proc~tdf_kpt->proc~utility_w0gauss proc~utility_recip_lattice_base->proc~utility_inv3 proc~wham_get_eig_deleig->proc~get_hh_r proc~pw90common_fourier_r_to_k pw90common_fourier_R_to_k proc~wham_get_eig_deleig->proc~pw90common_fourier_r_to_k proc~utility_diagonalize utility_diagonalize proc~wham_get_eig_deleig->proc~utility_diagonalize proc~wham_get_deleig_a wham_get_deleig_a proc~wham_get_eig_deleig->proc~wham_get_deleig_a proc~comms_bcast_char comms_bcast_char interface~comms_bcast->proc~comms_bcast_char proc~comms_bcast_cmplx comms_bcast_cmplx interface~comms_bcast->proc~comms_bcast_cmplx proc~comms_bcast_int comms_bcast_int interface~comms_bcast->proc~comms_bcast_int proc~comms_bcast_logical comms_bcast_logical interface~comms_bcast->proc~comms_bcast_logical proc~comms_bcast_real comms_bcast_real interface~comms_bcast->proc~comms_bcast_real proc~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_allreduce_cmplx->proc~comms_sync_error proc~comms_no_sync_allreduce_cmplx comms_no_sync_allreduce_cmplx proc~comms_allreduce_cmplx->proc~comms_no_sync_allreduce_cmplx proc~comms_allreduce_real->proc~comms_sync_error proc~comms_no_sync_allreduce_real comms_no_sync_allreduce_real proc~comms_allreduce_real->proc~comms_no_sync_allreduce_real zcopy zcopy proc~comms_no_sync_gatherv_cmplx_2->zcopy proc~comms_no_sync_gatherv_cmplx_3_4->zcopy dcopy dcopy proc~comms_no_sync_gatherv_real_2_3->dcopy proc~get_gauge_overlap_matrix->proc~get_win_min proc~utility_zgemmm utility_zgemmm proc~get_gauge_overlap_matrix->proc~utility_zgemmm 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~utility_diagonalize->proc~set_error_fatal zhpevx zhpevx proc~utility_diagonalize->zhpevx proc~utility_w0gauss->proc~set_error_input proc~wham_get_deleig_a->proc~utility_diagonalize proc~utility_rotate utility_rotate proc~wham_get_deleig_a->proc~utility_rotate proc~wham_get_deleig_a->proc~utility_rotate_diag 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~utility_matmul_diag utility_matmul_diag proc~utility_rotate_diag->proc~utility_matmul_diag proc~utility_zgemm_new utility_zgemm_new proc~utility_rotate_diag->proc~utility_zgemm_new proc~utility_zgemmm->proc~utility_zgemm_new zgemm zgemm proc~utility_zgemm_new->zgemm

Called by

proc~~boltzwann_main~~CalledByGraph proc~boltzwann_main boltzwann_main program~postw90 postw90 program~postw90->proc~boltzwann_main

Source Code

  subroutine boltzwann_main(pw90_boltzwann, dis_manifold, pw90_dos, kpt_latt, &
                            pw90_band_deriv_degen, postw90_oper, pw90_spin, physics, ws_region, &
                            w90_system, wannier_data, ws_distance, wigner_seitz, print_output, &
                            HH_R, SS_R, v_matrix, u_matrix, eigval, real_lattice, scissors_shift, &
                            mp_grid, num_wann, num_bands, num_kpts, effective_model, &
                            have_disentangled, spin_decomp, seedname, stdout, timer, error, comm)
    !================================================!
    !! This is the main routine of the BoltzWann module.
    !! It calculates the transport coefficients using the Boltzmann transport equation.
    !!
    !! It produces six files that contain:
    !!
    !!  1. the Transport Distribution function (TDF) in units of 1/hbar^2 * eV*fs/angstrom
    !!  2. the electrical conductivity in SI units (1/Ohm/m)
    !!  3. the tensor sigma*S (sigma=el.cond., S=seebeck) in SI units (Ampere/meter/K)
    !!  4. the Seebeck coefficient in SI units (V/K)
    !!  5. the thermal conductivity in SI units (W/meter/K)
    !!  6. if requested, the density of states
    !!
    !! Files from 2 to 4 are output on a grid of (mu,T) points, where mu is the chemical potential in eV and
    !! T is the temperature in Kelvin. The grid is defined in the input.
    !================================================!

    use w90_constants, only: dp
    use w90_io, only: io_stopwatch_start, io_stopwatch_stop
    use w90_comms, only: comms_bcast, w90_comm_type, mpirank
    use w90_types, only: dis_manifold_type, print_output_type, wannier_data_type, &
                         ws_region_type, w90_system_type, ws_distance_type, timer_list_type
    use w90_postw90_types, only: pw90_boltzwann_type, pw90_spin_mod_type, &
                                 pw90_band_deriv_degen_type, pw90_dos_mod_type, pw90_oper_read_type, wigner_seitz_type

    implicit none

    ! arguments
    type(pw90_boltzwann_type), intent(in) :: pw90_boltzwann
    type(dis_manifold_type), intent(in) :: dis_manifold
    type(pw90_dos_mod_type), intent(in) :: pw90_dos
    type(pw90_band_deriv_degen_type), intent(in) :: pw90_band_deriv_degen
    type(pw90_oper_read_type), intent(in) :: postw90_oper
    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(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) :: HH_R(:, :, :) !  <0n|r|Rm>
    complex(kind=dp), allocatable, intent(inout) :: SS_R(:, :, :, :) ! <0n|sigma_x,y,z|Rm>
    complex(kind=dp), intent(in) :: v_matrix(:, :, :), u_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), intent(in) :: kpt_latt(:, :)

    integer, intent(in) :: mp_grid(3)
    integer, intent(in) :: num_wann, num_bands, num_kpts
    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
    integer :: TempNumPoints, MuNumPoints, TDFEnergyNumPoints
    integer :: i, j, ierr, EnIdx, TempIdx, MuIdx
    real(kind=dp), allocatable :: TempArray(:), MuArray(:), KTArray(:)
    real(kind=dp), allocatable :: TDF(:, :, :) ! (coordinate,Energy)
    real(kind=dp), allocatable :: TDFEnergyArray(:)
    real(kind=dp), allocatable :: IntegrandArray(:, :) ! (coordinate, Energy) at a given T and mu
    real(kind=dp) :: SigmaS_FP(3, 3), ThisElCond(3, 3), ElCondInverse(3, 3), ThisSeebeck(3, 3)
    real(kind=dp) :: ThisElCond2d(2, 2), ElCondInverse2d(2, 2)
    !real(kind=dp), dimension(6) :: ElCondTimesSeebeck
    real(kind=dp), allocatable :: ElCond(:, :, :) ! (coordinate,Temp, mu)
    real(kind=dp), allocatable :: SigmaS(:, :, :) ! (coordinate,Temp, mu)
    real(kind=dp), allocatable :: Seebeck(:, :, :) ! (coordinate,Temp, mu)
    real(kind=dp), allocatable :: Kappa(:, :, :) ! (coordinate,Temp, mu)
    real(kind=dp), allocatable :: LocalElCond(:, :) ! (coordinate,Temp+mu combined index)
    real(kind=dp), allocatable :: LocalSigmaS(:, :) ! (coordinate,Temp+mu combined index)
    real(kind=dp), allocatable :: LocalSeebeck(:, :) ! (coordinate,Temp+mu combined index)
    real(kind=dp), allocatable :: LocalKappa(:, :) ! (coordinate,Temp+mu combined index)
    real(kind=dp) :: Determinant
    integer :: tdf_unit, elcond_unit, sigmas_unit, seebeck_unit, kappa_unit, ndim

    integer :: LocalIdx, GlobalIdx
    ! I also add 3 times the smearing on each side of the TDF energy array to take into account also possible smearing effects
    real(kind=dp), parameter :: TDF_exceeding_energy_times_smr = 3._dp
    real(kind=dp) :: TDF_exceeding_energy
    real(kind=dp) :: cell_volume
    integer :: NumberZeroDet

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

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

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

    if (print_output%iprint > 0 .and. print_output%timing_level > 0) call io_stopwatch_start('boltzwann_main', timer)

    if (print_output%iprint > 0) then
      write (stdout, *)
      write (stdout, '(1x,a)') '*---------------------------------------------------------------------------*'
      write (stdout, '(1x,a)') '|                   Boltzmann Transport (BoltzWann module)                  |'
      write (stdout, '(1x,a)') '*---------------------------------------------------------------------------*'
      write (stdout, '(1x,a)') '| '//pub_string_1//'|'
      write (stdout, '(1x,a)') '| '//pub_string_2//'|'
      write (stdout, '(1x,a)') '| '//pub_string_3//'|'
      write (stdout, '(1x,a)') '| '//pub_string_4//'|'
      write (stdout, '(1x,a)') '*---------------------------------------------------------------------------*'
      write (stdout, *)
    end if

    if (print_output%iprint > 0) then
      if (pw90_boltzwann%dir_num_2d /= 0) then
        write (stdout, '(1x,a)') '>                                                                           <'
        write (stdout, '(1x,a)') '> NOTE! Using the 2D version for the calculation of the Seebeck             <'
        if (pw90_boltzwann%dir_num_2d == 1) then
          write (stdout, '(1x,a)') '>       coefficient, where the non-periodic direction is x.                 <'
        elseif (pw90_boltzwann%dir_num_2d == 2) then
          write (stdout, '(1x,a)') '>       coefficient, where the non-periodic direction is y.                 <'
        elseif (pw90_boltzwann%dir_num_2d == 3) then
          write (stdout, '(1x,a)') '>       coefficient, where the non-periodic direction is z.                 <'
        end if
        write (stdout, '(1x,a)') '>                                                                           <'
        write (stdout, '(1x,a)') ''
      end if
    end if

    ! separate error condition from info printout above (which may be avoided entirely in lib mode?)
    if (pw90_boltzwann%dir_num_2d < 0 .or. pw90_boltzwann%dir_num_2d > 3) then
      call set_error_input(error, 'Unrecognized value of pw90_boltzwann_2d_dir_num', comm)
      return
    end if

    ! I precalculate the TempArray and the MuArray
    TempNumPoints = int(floor((pw90_boltzwann%temp_max - pw90_boltzwann%temp_min)/pw90_boltzwann%temp_step)) + 1
    allocate (TempArray(TempNumPoints), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating TempArray in boltzwann_main', comm)
      return
    end if
    do i = 1, TempNumPoints
      TempArray(i) = pw90_boltzwann%temp_min + real(i - 1, dp)*pw90_boltzwann%temp_step
    end do

    ! This array contains the same temperatures of the TempArray, but multiplied by k_boltzmann, in units of eV
    allocate (KTArray(TempNumPoints), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating KTArray in boltzwann_main', comm)
      return
    end if
    ! (k_B in eV/kelvin is equal to k_B_SI / elem_charge_SI)
    KTArray = TempArray*physics%k_B_SI/physics%elem_charge_SI

    MuNumPoints = int(floor((pw90_boltzwann%mu_max - pw90_boltzwann%mu_min)/pw90_boltzwann%mu_step)) + 1
    allocate (MuArray(MuNumPoints), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating MuArray in boltzwann_main', comm)
      return
    end if
    do i = 1, MuNumPoints
      MuArray(i) = pw90_boltzwann%mu_min + real(i - 1, dp)*pw90_boltzwann%mu_step
    end do

    if (pw90_boltzwann%tdf_smearing%use_adaptive) then
      call set_error_input(error, 'Adaptive smearing not allowed in Boltzwann TDF', comm)
      return
    end if
    ! I precalculate the TDFEnergyArray
    ! I assume that dis_win_min and dis_win_max are set to sensible values, related to the max and min energy
    ! This is true if the .eig file is present. I can assume its presence since we need it to interpolate the
    ! bands.
    ! I also add 3 times the smearing on each side of the TDF energy array to take into account also possible smearing effects,
    ! or at least 0.2 eV
    TDF_exceeding_energy = max(TDF_exceeding_energy_times_smr*pw90_boltzwann%tdf_smearing%fixed_width, 0.2_dp)
    TDFEnergyNumPoints = int(floor((dis_manifold%win_max - dis_manifold%win_min &
                                    + 2._dp*TDF_exceeding_energy)/pw90_boltzwann%tdf_energy_step)) + 1
    if (TDFEnergyNumPoints .eq. 1) TDFEnergyNumPoints = 2
    allocate (TDFEnergyArray(TDFEnergyNumPoints), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating TDFEnergyArray in boltzwann_main', comm)
      return
    end if
    do i = 1, TDFEnergyNumPoints
      TDFEnergyArray(i) = dis_manifold%win_min - TDF_exceeding_energy &
                          + real(i - 1, dp)*pw90_boltzwann%tdf_energy_step
    end do

    if (spin_decomp) then
      ndim = 3
    else
      ndim = 1
    end if

    ! I allocate the array for the TDF
    allocate (TDF(6, TDFEnergyNumPoints, ndim), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating TDF in boltzwann_main', comm)
      return
    end if

    ! I call the subroutine that calculates the Transport Distribution Function
    call calcTDFandDOS(pw90_boltzwann, dis_manifold, pw90_dos, kpt_latt, postw90_oper, &
                       pw90_band_deriv_degen, pw90_spin, ws_region, print_output, wannier_data, &
                       ws_distance, wigner_seitz, HH_R, SS_R, u_matrix, v_matrix, eigval, &
                       real_lattice, TDF, TDFEnergyArray, cell_volume, scissors_shift, mp_grid, &
                       num_bands, num_kpts, num_wann, w90_system%num_valence_bands, &
                       w90_system%num_elec_per_state, effective_model, have_disentangled, &
                       spin_decomp, seedname, stdout, timer, error, comm)
    if (allocated(error)) return

    ! The TDF array contains now the TDF, or more precisely
    ! hbar^2 * TDF in units of eV * fs / angstrom

    ! I print on file the TDF
    if (on_root) then
      open (newunit=tdf_unit, file=trim(seedname)//'_tdf.dat')
      write (tdf_unit, '(A)') "# Written by the BoltzWann module of the Wannier90 code."
      write (tdf_unit, '(A)') "# Transport distribution function (in units of 1/hbar^2 * eV * fs / angstrom)"// &
        " vs energy in eV"
      write (tdf_unit, '(A)') "# Content of the columns:"
      write (tdf_unit, '(A)') "# Energy TDF_xx TDF_xy TDF_yy TDF_xz TDF_yz TDF_zz"
      write (tdf_unit, '(A)') '#   (if spin decomposition is required, 12 further columns are provided, with the 6'
      write (tdf_unit, '(A)') '#    components of the TDF for the spin up, followed by those for the spin down)'
      do i = 1, size(TDFEnergyArray)
        if (ndim .eq. 1) then
          write (tdf_unit, 101) TDFEnergyArray(i), TDF(:, i, 1)
        else
          write (tdf_unit, 102) TDFEnergyArray(i), TDF(:, i, :)
        end if
      end do
      close (tdf_unit)
      if (print_output%iprint > 1) &
        write (stdout, '(3X,A)') "Transport distribution function written on the "//trim(seedname)//"_tdf.dat file."
    end if

    ! *********************************************************************************
    ! I got the TDF and I printed it. Now I use it to calculate the transport properties.

    if (on_root .and. (print_output%timing_level > 0)) call io_stopwatch_start('boltzwann_main: calc_props', timer)

    ! I obtain the counts and displs arrays, which tell how I should partition a big array
    ! on the different nodes.
    call comms_array_split(TempNumPoints*MuNumPoints, counts, displs, comm)

    ! I allocate the arrays for the spectra
    ! Allocate at least 1 entry
    allocate (LocalElCond(6, max(1, counts(my_node_id))), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating LocalElCond in boltzwann_main', comm)
      return
    end if
    allocate (LocalSigmaS(6, max(1, counts(my_node_id))), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating LocalSigmaS in boltzwann_main', comm)
      return
    end if
    allocate (LocalSeebeck(9, max(1, counts(my_node_id))), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating LocalSeebeck in boltzwann_main', comm)
      return
    end if
    allocate (LocalKappa(6, max(1, counts(my_node_id))), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating LocalKappa in boltzwann_main', comm)
      return
    end if
    LocalSigmaS = 0._dp
    LocalElCond = 0._dp
    LocalSeebeck = 0._dp
    LocalKappa = 0._dp

    ! I allocate the array that I will use to store the functions to be integrated
    allocate (IntegrandArray(6, TDFEnergyNumPoints), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating FermiDerivArray in boltzwann_main', comm)
      return
    end if

    NumberZeroDet = 0
    ! Now, I calculate the various spectra for all mu and T values
    do LocalIdx = 1, counts(my_node_id)
      ! GlobalIdx is an index from 0 to TempNumPoints*MuNumPoints-1
      GlobalIdx = displs(my_node_id) + (LocalIdx - 1)
      ! MuIdx goes from 1 to MuNumPoints
      MuIdx = GlobalIdx/TempNumPoints + 1
      ! TempIdx goes from 1 to TempNumPoints
      TempIdx = GlobalIdx - TempNumPoints*(MuIdx - 1) + 1

      ! For the calculation of the properties, I only use the total TDF (i.e., unresolved for spin)
      IntegrandArray = TDF(:, :, 1)
      do EnIdx = 1, TDFEnergyNumPoints
        IntegrandArray(:, EnIdx) = IntegrandArray(:, EnIdx)* &
                                   MinusFermiDerivative(E=TDFEnergyArray(EnIdx), mu=MuArray(MuIdx), KT=KTArray(TempIdx))
      end do
      ! Now, IntegrandArray contains (-dn/dE) * TDF_ij(E), where n is the Fermi distribution function
      ! Its integral is ElCond_ij/e^2
      LocalElCond(:, LocalIdx) = sum(IntegrandArray, DIM=2)*pw90_boltzwann%tdf_energy_step
      ! ElCond contains now (hbar^2/e^2) * sigma in eV*fs/angstrom, where sigma is the conductivity tensor
      ! (note that MinusFermiDerivative is in units of 1/eV, so that when I perform the integration
      ! ElCond has the same units of TDF)

      ! I store in ThisElCond the conductivity tensor in standard format
      ! Store both the upper and lower triangle
      do j = 1, 3
        do i = 1, j
          ThisElCond(i, j) = LocalElCond(i + ((j - 1)*j)/2, LocalIdx)
          ThisElCond(j, i) = LocalElCond(i + ((j - 1)*j)/2, LocalIdx)
        end do
      end do

      ! I calculate the inverse matrix of the conductivity
      if (pw90_boltzwann%dir_num_2d /= 0) then
        ! Invert only the appropriate 2x2 submatrix
        if (pw90_boltzwann%dir_num_2d == 1) then
          ThisElCond2d(1, 1) = ThisElCond(2, 2)
          ThisElCond2d(1, 2) = ThisElCond(2, 3)
          ThisElCond2d(2, 1) = ThisElCond(3, 2)
          ThisElCond2d(2, 2) = ThisElCond(3, 3)
        elseif (pw90_boltzwann%dir_num_2d == 2) then
          ThisElCond2d(1, 1) = ThisElCond(1, 1)
          ThisElCond2d(1, 2) = ThisElCond(1, 3)
          ThisElCond2d(2, 1) = ThisElCond(3, 1)
          ThisElCond2d(2, 2) = ThisElCond(3, 3)
        elseif (pw90_boltzwann%dir_num_2d == 3) then
          ThisElCond2d(1, 1) = ThisElCond(1, 1)
          ThisElCond2d(1, 2) = ThisElCond(1, 2)
          ThisElCond2d(2, 1) = ThisElCond(2, 1)
          ThisElCond2d(2, 2) = ThisElCond(2, 2)
          ! I do not do the else case because this was already checked at the beginning
          ! of this routine
        end if
        call utility_inv2(ThisElCond2d, ElCondInverse2d, Determinant)
        ElCondInverse = 0._dp ! Other elements must be set to zero
        if (pw90_boltzwann%dir_num_2d == 1) then
          ElCondInverse(2, 2) = ElCondInverse2d(1, 1)
          ElCondInverse(2, 3) = ElCondInverse2d(1, 2)
          ElCondInverse(3, 2) = ElCondInverse2d(2, 1)
          ElCondInverse(3, 3) = ElCondInverse2d(2, 2)
        elseif (pw90_boltzwann%dir_num_2d == 2) then
          ElCondInverse(1, 1) = ElCondInverse2d(1, 1)
          ElCondInverse(1, 3) = ElCondInverse2d(1, 2)
          ElCondInverse(3, 1) = ElCondInverse2d(2, 1)
          ElCondInverse(3, 3) = ElCondInverse2d(2, 2)
        elseif (pw90_boltzwann%dir_num_2d == 3) then
          ElCondInverse(1, 1) = ElCondInverse2d(1, 1)
          ElCondInverse(1, 2) = ElCondInverse2d(1, 2)
          ElCondInverse(2, 1) = ElCondInverse2d(2, 1)
          ElCondInverse(2, 2) = ElCondInverse2d(2, 2)
          ! I do not do the else case because this was already checked at the beginning
          ! of this routine
        end if
      else
        call utility_inv3(ThisElCond, ElCondInverse, Determinant)
      end if
      if (Determinant .eq. 0._dp) then
        NumberZeroDet = NumberZeroDet + 1
        ! If the determinant is zero (i.e., zero conductivity along a given direction)
        ! I set the Inverse to zero, so that the Seebeck is zero. I will also issue
        ! a warning.
        ElCondInverse = 0._dp
      else
        ! The routine returns the adjoint of ThisElCond; in order to get
        ! the inverse, I have to calculate ElCondInverse / Determinant
        ElCondInverse = ElCondInverse/Determinant
      end if

      ! Now, I multiply IntegrandArray by (E-mu): then, IntegrandArray contains
      ! (-dn/dE) * TDF_ij(E) * (E-mu) and its integral is (ElCond*Seebeck)_ij * T / e, where
      ! T is the temperature
      do EnIdx = 1, TDFEnergyNumPoints
        IntegrandArray(:, EnIdx) = IntegrandArray(:, EnIdx)*(TDFEnergyArray(EnIdx) - MuArray(MuIdx))
      end do

      ! I store in SigmaS_FP the product of the two tensors in full-packed format
      LocalSigmaS(:, LocalIdx) = sum(IntegrandArray, DIM=2)*pw90_boltzwann%tdf_energy_step/TempArray(TempIdx)
      do j = 1, 3
        do i = 1, j
          ! Both upper and lower diagonal
          SigmaS_FP(i, j) = LocalSigmaS(i + ((j - 1)*j)/2, LocalIdx)
          SigmaS_FP(j, i) = LocalSigmaS(i + ((j - 1)*j)/2, LocalIdx)
        end do
      end do
      ! Now, LocalSigmaS (and SigmaS_FP) contain
      ! [ElCond*Seebeck] * hbar^2/e in units of eV^2*fs/angstrom/kelvin

      ! I calculate ElCond^(-1) . (LocalSigmaS) = Seebeck in fully-packed format and then
      ! store it in the LocalSeebeck array (not in packed format because, unless S and sigma
      ! commute, there is no a-priori reason for which S should be symmetric - even if
      ! probably one can find physical reasons).

      ! I invert the sign because the electron charge is < 0
      ThisSeebeck = -matmul(ElCondInverse, SigmaS_FP)
      ! Reshuffle in 1D; order: xx, xy, xz, yx, yy, yz, zx, zy, zz
      ! Note that this is different
      do j = 1, 3
        do i = 1, 3
          LocalSeebeck((i - 1)*3 + j, LocalIdx) = ThisSeebeck(i, j)
        end do
      end do

      ! Now, Seebeck contains the Seebeck coefficient in volt / Kelvin. In fact:
      ! - ElCond contains (hbar^2/e^2) * sigma in eV*fs/angstrom, where sigma is the value of the
      !   conductivity, i.e. it is the conductivity in units of (e^2/hbar^2) * eV * fs / angstrom
      ! - ElCondInverse is in thus in units of (hbar^2/e^2) / eV / fs * angstrom
      ! - ElCondTimesSeebeck is in units of  e / hbar^2 * eV^2 * fs / angstrom / Kelvin
      ! therefore ThisSeebeck, which has the units of ElCondInverse * ElCondTimesSeebeck, i.e.
      ! [(hbar^2/e^2) / eV / fs * angstrom] * [ e / hbar^2 * eV^2 * fs / angstrom / kelvin] =
      ! eV/e/kelvin = volt/kelvin

      ! Now, I multiply IntegrandArray by (E-mu): then, IntegrandArray contains
      ! (-dn/dE) * TDF_ij(E) * (E-mu)^2 and its integral is (Kappa)_ij * T
      do EnIdx = 1, TDFEnergyNumPoints
        IntegrandArray(:, EnIdx) = IntegrandArray(:, EnIdx)*(TDFEnergyArray(EnIdx) - MuArray(MuIdx))
      end do
      LocalKappa(:, LocalIdx) = sum(IntegrandArray, DIM=2)*pw90_boltzwann%tdf_energy_step/TempArray(TempIdx)
      ! Kappa contains now the thermal conductivity in units of
      ! 1/hbar^2 * eV^3*fs/angstrom/kelvin
    end do
    ! I check if there were (mu,T) pairs for which we got sigma = 0
    call comms_reduce(NumberZeroDet, 1, 'SUM', error, comm)
    if (allocated(error)) return
    if (on_root) then
      if ((NumberZeroDet .gt. 0)) then
        write (stdout, '(1X,A,I0,A)') "> WARNING! There are ", NumberZeroDet, " (mu,T) pairs for which the electrical"
        write (stdout, '(1X,A)') ">          conductivity has zero determinant."
        write (stdout, '(1X,A)') ">          Seebeck coefficient set to zero for those pairs."
        write (stdout, '(1X,A)') ">          Check if this is physical or not."
        write (stdout, '(1X,A)') ">          (If you are dealing with a 2D system, set the pw90_boltzwann_2d_dir flag.)"
        write (stdout, '(1X,A)') ""
      end if
    end if

    ! Now, I multiply by the correct factors to obtain the tensors in SI units

    ! **** Electrical conductity ****
    ! Now, ElCond is in units of (e^2/hbar^2) * eV * fs / angstrom
    ! I want the conductivity in units of 1/Ohm/meter (i.e., SI units). The conversion factor is then
    ! e^2/hbar^2 * eV*fs/angstrom * Ohm * meter; we use the constants in constants.F90 by noting that:
    ! Ohm = V/A = V*s/C, where A=ampere, C=coulomb
    ! Then we have: (e/C)^3/hbar^2 * V*s^2 * V * C^2 * (meter / angstrom)  * (fs / s)
    ! Now: e/C = elem_charge_SI; CV=Joule, CV/s=Watt, hbar/Watt = hbar_SI;
    ! moreover meter / angstrom = 1e10, fs / s = 1e-15 so that we finally get the
    ! CONVERSION FACTOR: elem_charge_SI**3 / (hbar_SI**2) * 1.e-5_dp
    LocalElCond = LocalElCond*physics%elem_charge_SI**3/(physics%hbar_SI**2)*1.e-5_dp
    ! THIS IS NOW THE ELECTRICAL CONDUCTIVITY IN SI UNITS, i.e. in 1/Ohm/meter

    ! *** Sigma * S ****
    ! Again, as above or below for Kappa, the conversion factor is
    ! * elem_charge_SI**3 / (hbar_SI**2) * 1.e-5_dp
    ! and brings the result to Ampere/m/K
    LocalSigmaS = LocalSigmaS*physics%elem_charge_SI**3/(physics%hbar_SI**2)*1.e-5_dp

    ! **** Seebeck coefficient ****
    ! THE SEEECK COEFFICIENTS IS ALREADY IN volt/kelvin, so nothing has to be done

    ! **** K coefficient (approx. the Thermal conductivity) ****
    ! Now, Kappa is in units of 1/hbar^2 * eV^3*fs/angstrom/K
    ! I want it to be in units of W/m/K (i.e., SI units). Then conversion factor is then
    ! 1/hbar^2 * eV^3 * fs/angstrom/K / W * m * K =
    ! 1/hbar^2 * eV^3 * fs / W * (m/angstrom) =
    ! 1/hbar^2 * C^3 * V^3 * s / W * [e/C]^3 * (m/angstrom) * (fs / s) =
    ! 1/hbar^2 * J^2 * [e/C]^3 * (m/angstrom) * (fs / s) =
    ! elem_charge_SI**3 / (hbar_SI**2) * 1.e-5_dp, i.e. the same conversion factor as above
    LocalKappa = LocalKappa*physics%elem_charge_SI**3/(physics%hbar_SI**2)*1.e-5_dp
    ! THIS IS NOW THE THERMAL CONDUCTIVITY IN SI UNITS, i.e. in W/meter/K

    ! Now I send the different pieces to the local node
    if (on_root) then
      allocate (ElCond(6, TempNumPoints, MuNumPoints), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating ElCond in boltzwann_main', comm)
        return
      end if
      allocate (SigmaS(6, TempNumPoints, MuNumPoints), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating SigmaS in boltzwann_main', comm)
        return
      end if
      allocate (Seebeck(9, TempNumPoints, MuNumPoints), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating Seebeck in boltzwann_main', comm)
        return
      end if
      allocate (Kappa(6, TempNumPoints, MuNumPoints), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating Kappa in boltzwann_main', comm)
        return
      end if
    else
      ! In principle, this should not be needed, because we use ElCond,
      ! Seebeck and Kappa only on the root node. However, since all
      ! processors call comms_gatherv a few lines below, and one argument
      ! is ElCond(1,1,1), some compilers complain.
      allocate (ElCond(1, 1, 1), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating ElCond in boltzwann_main (2)', comm)
        return
      end if
      allocate (SigmaS(1, 1, 1), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating SigmaS in boltzwann_main (2)', comm)
        return
      end if
      allocate (Seebeck(1, 1, 1), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating Seebeck in boltzwann_main (2)', comm)
        return
      end if
      allocate (Kappa(1, 1, 1), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating Kappa in boltzwann_main (2)', comm)
        return
      end if
    end if

    ! The 6* factors are due to the fact that for each (T,mu) pair we have 6 components (xx,xy,yy,xz,yz,zz)
    ! NOTE THAT INSTEAD SEEBECK IS A FULL MATRIX AND HAS 9 COMPONENTS!
    call comms_gatherv(LocalElCond, 6*counts(my_node_id), ElCond, 6*counts, 6*displs, error, comm)
    if (allocated(error)) return
    call comms_gatherv(LocalSigmaS, 6*counts(my_node_id), SigmaS, 6*counts, 6*displs, error, comm)
    if (allocated(error)) return
    call comms_gatherv(LocalSeebeck, 9*counts(my_node_id), Seebeck, 9*counts, 9*displs, error, comm)
    if (allocated(error)) return
    call comms_gatherv(LocalKappa, 6*counts(my_node_id), Kappa, 6*counts, 6*displs, error, comm)
    if (allocated(error)) return

    if (on_root .and. (print_output%timing_level > 0)) call io_stopwatch_stop('boltzwann_main: calc_props', timer)

    ! Open files and print
    if (on_root) then
      open (newunit=elcond_unit, file=trim(seedname)//'_elcond.dat')
      write (elcond_unit, '(A)') "# Written by the BoltzWann module of the Wannier90 code."
      write (elcond_unit, '(A)') "# [Electrical conductivity in SI units, i.e. in 1/Ohm/m]"
      write (elcond_unit, '(A)') "# Mu(eV) Temp(K) ElCond_xx ElCond_xy ElCond_yy ElCond_xz ElCond_yz ElCond_zz"
      do MuIdx = 1, MuNumPoints
        do TempIdx = 1, TempNumPoints
          write (elcond_unit, 103) MuArray(MuIdx), TempArray(TempIdx), ElCond(:, TempIdx, MuIdx)
        end do
      end do
      close (elcond_unit)
      if (print_output%iprint > 1) &
        write (stdout, '(3X,A)') "Electrical conductivity written on the "//trim(seedname)//"_elcond.dat file."

      open (newunit=sigmas_unit, file=trim(seedname)//'_sigmas.dat')
      write (sigmas_unit, '(A)') "# Written by the BoltzWann module of the Wannier90 code."
      write (sigmas_unit, '(A)') "# [(Electrical conductivity * Seebeck coefficient) in SI units, i.e. in Ampere/m/K]"
      write (sigmas_unit, '(A)') "# Mu(eV) Temp(K) (Sigma*S)_xx (Sigma*S)_xy (Sigma*S)_yy (Sigma*S)_xz (Sigma*S)_yz (Sigma*S)_zz"
      do MuIdx = 1, MuNumPoints
        do TempIdx = 1, TempNumPoints
          write (sigmas_unit, 103) MuArray(MuIdx), TempArray(TempIdx), SigmaS(:, TempIdx, MuIdx)
        end do
      end do
      close (sigmas_unit)
      if (print_output%iprint > 1) write (stdout, '(3X,A)') &
        "sigma*S (sigma=el. conductivity, S=Seebeck coeff.) written on the "//trim(seedname)//"_sigmas.dat file."

      open (newunit=seebeck_unit, file=trim(seedname)//'_seebeck.dat')
      write (seebeck_unit, '(A)') "# Written by the BoltzWann module of the Wannier90 code."
      write (seebeck_unit, '(A)') "# [Seebeck coefficient in SI units, i.e. in V/K]"
      write (seebeck_unit, '(A)') &
        "# Mu(eV) Temp(K) Seebeck_xx Seebeck_xy Seebeck_xz Seebeck_yx Seebeck_yy Seebeck_yz Seebeck_zx Seebeck_zy Seebeck_zz"
      do MuIdx = 1, MuNumPoints
        do TempIdx = 1, TempNumPoints
          write (seebeck_unit, 104) MuArray(MuIdx), TempArray(TempIdx), Seebeck(:, TempIdx, MuIdx)
        end do
      end do
      close (seebeck_unit)
      if (print_output%iprint > 1) &
        write (stdout, '(3X,A)') "Seebeck coefficient written on the "//trim(seedname)//"_seebeck.dat file."

      open (newunit=kappa_unit, file=trim(seedname)//'_kappa.dat')
      write (kappa_unit, '(A)') "# Written by the BoltzWann module of the Wannier90 code."
      write (kappa_unit, '(A)') "# [K coefficient in SI units, i.e. in W/m/K]"
      write (kappa_unit, '(A)') "# [the K coefficient is defined in the documentation, and is an ingredient of"
      write (kappa_unit, '(A)') "#  the thermal conductivity. See the docs for further information.]"
      write (kappa_unit, '(A)') "# Mu(eV) Temp(K) Kappa_xx Kappa_xy Kappa_yy Kappa_xz Kappa_yz Kappa_zz"
      do MuIdx = 1, MuNumPoints
        do TempIdx = 1, TempNumPoints
          write (kappa_unit, 103) MuArray(MuIdx), TempArray(TempIdx), Kappa(:, TempIdx, MuIdx)
        end do
      end do
      close (kappa_unit)
      if (print_output%iprint > 1) &
        write (stdout, '(3X,A)') "K coefficient written on the "//trim(seedname)//"_kappa.dat file."
    end if

    if (on_root) then
      write (stdout, '(3X,A)') "Transport properties calculated."
      write (stdout, *)
      write (stdout, '(1x,a)') '*---------------------------------------------------------------------------*'
      write (stdout, '(1x,a)') '|                        End of the BoltzWann module                        |'
      write (stdout, '(1x,a)') '*---------------------------------------------------------------------------*'
      write (stdout, *)
    end if

    ! Before ending, I deallocate memory
    deallocate (TempArray, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating TempArray in boltzwann_main', comm)
      return
    end if
    deallocate (KTArray, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating KTArray in boltzwann_main', comm)
      return
    end if
    deallocate (MuArray, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating MuArray in boltzwann_main', comm)
      return
    end if
    deallocate (TDFEnergyArray, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating TDFEnergyArray in boltzwann_main', comm)
      return
    end if
    deallocate (TDF, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating TDF in boltzwann_main', comm)
      return
    end if
    deallocate (LocalElCond, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating LocalElCond in boltzwann_main', comm)
      return
    end if
    deallocate (LocalSigmaS, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating LocalSigmaS in boltzwann_main', comm)
      return
    end if
    deallocate (LocalSeebeck, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating LocalSeebeck in boltzwann_main', comm)
      return
    end if
    deallocate (LocalKappa, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating LocalKappa in boltzwann_main', comm)
      return
    end if

    deallocate (ElCond, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating ElCond in boltzwann_main', comm)
      return
    end if
    deallocate (SigmaS, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating SigmaS in boltzwann_main', comm)
      return
    end if
    deallocate (Seebeck, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating Seebeck in boltzwann_main', comm)
      return
    end if
    deallocate (Kappa, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating Kappa in boltzwann_main', comm)
      return
    end if

    deallocate (IntegrandArray, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating IntegrandArray in boltzwann_main', comm)
      return
    end if

    if (on_root .and. (print_output%timing_level > 0)) call io_stopwatch_stop('boltzwann_main', timer)

101 FORMAT(7G18.10)
102 FORMAT(19G18.10)
103 FORMAT(8G18.10)
104 FORMAT(11G18.10)

  end subroutine boltzwann_main