boltzwann_main Subroutine

public subroutine boltzwann_main()

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

None

Calls

proc~~boltzwann_main~~CallsGraph proc~boltzwann_main boltzwann_main 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 io_stopwatch proc~boltzwann_main->proc~io_stopwatch proc~io_error io_error proc~boltzwann_main->proc~io_error proc~io_file_unit io_file_unit proc~boltzwann_main->proc~io_file_unit proc~minusfermiderivative MinusFermiDerivative proc~boltzwann_main->proc~minusfermiderivative interface~comms_reduce comms_reduce proc~boltzwann_main->interface~comms_reduce proc~utility_inv3 utility_inv3 proc~boltzwann_main->proc~utility_inv3 proc~utility_inv2 utility_inv2 proc~boltzwann_main->proc~utility_inv2 interface~comms_gatherv comms_gatherv proc~boltzwann_main->interface~comms_gatherv proc~calctdfanddos->proc~io_stopwatch proc~calctdfanddos->proc~io_error proc~calctdfanddos->proc~io_file_unit proc~calctdfanddos->interface~comms_reduce proc~param_get_smearing_type param_get_smearing_type proc~calctdfanddos->proc~param_get_smearing_type proc~tdf_kpt TDF_kpt proc~calctdfanddos->proc~tdf_kpt interface~comms_allreduce comms_allreduce proc~calctdfanddos->interface~comms_allreduce get_ss_r get_ss_r proc~calctdfanddos->get_ss_r get_hh_r get_hh_r proc~calctdfanddos->get_hh_r proc~wham_get_eig_deleig wham_get_eig_deleig proc~calctdfanddos->proc~wham_get_eig_deleig 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~io_stopwatch->proc~io_error proc~comms_reduce_cmplx comms_reduce_cmplx interface~comms_reduce->proc~comms_reduce_cmplx proc~comms_reduce_int comms_reduce_int interface~comms_reduce->proc~comms_reduce_int proc~comms_reduce_real comms_reduce_real interface~comms_reduce->proc~comms_reduce_real proc~comms_gatherv_real comms_gatherv_real interface~comms_gatherv->proc~comms_gatherv_real proc~spin_get_nk spin_get_nk proc~tdf_kpt->proc~spin_get_nk proc~utility_w0gauss utility_w0gauss proc~tdf_kpt->proc~utility_w0gauss proc~comms_allreduce_real comms_allreduce_real interface~comms_allreduce->proc~comms_allreduce_real proc~comms_allreduce_cmplx comms_allreduce_cmplx interface~comms_allreduce->proc~comms_allreduce_cmplx proc~wham_get_eig_deleig->get_hh_r pw90common_fourier_r_to_k pw90common_fourier_r_to_k proc~wham_get_eig_deleig->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~dos_get_k->proc~io_error proc~dos_get_k->proc~spin_get_nk proc~dos_get_k->proc~utility_w0gauss interface~pw90common_kmesh_spacing pw90common_kmesh_spacing proc~dos_get_levelspacing->interface~pw90common_kmesh_spacing ss_r ss_r proc~spin_get_nk->ss_r proc~spin_get_nk->pw90common_fourier_r_to_k proc~utility_rotate_diag utility_rotate_diag proc~spin_get_nk->proc~utility_rotate_diag proc~spin_get_nk->proc~utility_diagonalize proc~utility_w0gauss->proc~io_error proc~utility_matmul_diag utility_matmul_diag proc~utility_rotate_diag->proc~utility_matmul_diag proc~utility_diagonalize->proc~io_error zhpevx zhpevx proc~utility_diagonalize->zhpevx proc~wham_get_deleig_a->proc~utility_rotate_diag proc~wham_get_deleig_a->proc~utility_diagonalize proc~utility_rotate utility_rotate proc~wham_get_deleig_a->proc~utility_rotate 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 dcopy dcopy proc~comms_gatherv_real->dcopy
Help

Source Code


Source Code

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

    ! Needed to split an array on different nodes
    integer, dimension(0:num_nodes-1) :: counts
    integer, dimension(0:num_nodes-1) :: displs
    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
    integer :: NumberZeroDet

    if(on_root .and. (timing_level>0)) call io_stopwatch('boltzwann_main',1)

    if(on_root) 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 (on_root) then
       if (boltz_2d_dir_num /= 0) then
          write(stdout,'(1x,a)') '>                                                                           <'
          write(stdout,'(1x,a)') '> NOTE! Using the 2D version for the calculation of the Seebeck             <'
          if (boltz_2d_dir_num == 1) then
             write(stdout,'(1x,a)') '>       coefficient, where the non-periodic direction is x.                 <'
          elseif (boltz_2d_dir_num == 2) then
             write(stdout,'(1x,a)') '>       coefficient, where the non-periodic direction is y.                 <'
          elseif (boltz_2d_dir_num == 3) then
             write(stdout,'(1x,a)') '>       coefficient, where the non-periodic direction is z.                 <'
          else
             call io_error('Unrecognized value of boltz_2d_dir_num')
          end if
          write(stdout,'(1x,a)') '>                                                                           <'
          write(stdout,'(1x,a)') ''
       end if
    end if


    ! I precalculate the TempArray and the MuArray
    TempNumPoints = int(floor((boltz_temp_max-boltz_temp_min)/boltz_temp_step))+1
    allocate(TempArray(TempNumPoints),stat=ierr)
    if (ierr/=0) call io_error('Error in allocating TempArray in boltzwann_main')
    do i=1,TempNumPoints
       TempArray(i) = boltz_temp_min + real(i-1,dp)*boltz_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) call io_error('Error in allocating KTArray in boltzwann_main')
    ! (k_B in eV/kelvin is equal to k_B_SI / elem_charge_SI)
    KTArray = TempArray * k_B_SI / elem_charge_SI

    MuNumPoints = int(floor((boltz_mu_max-boltz_mu_min)/boltz_mu_step))+1
    allocate(MuArray(MuNumPoints),stat=ierr)
    if (ierr/=0) call io_error('Error in allocating MuArray in boltzwann_main')
    do i=1,MuNumPoints
       MuArray(i) = boltz_mu_min + real(i-1,dp)*boltz_mu_step
    end do

    ! 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 * boltz_TDF_smr_fixed_en_width,0.2_dp)
    TDFEnergyNumPoints = int(floor((dis_win_max-dis_win_min+2._dp*TDF_exceeding_energy)/boltz_tdf_energy_step))+1
    if (TDFEnergyNumPoints .eq. 1) TDFEnergyNumPoints = 2
    allocate(TDFEnergyArray(TDFEnergyNumPoints),stat=ierr)
    if (ierr/=0) call io_error('Error in allocating TDFEnergyArray in boltzwann_main')
    do i=1,TDFEnergyNumPoints
       TDFEnergyArray(i) = dis_win_min - TDF_exceeding_energy + real(i-1,dp)*boltz_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) call io_error('Error in allocating TDF in boltzwann_main')

    ! I call the subroutine that calculates the Transport Distribution Function
    call calcTDFandDOS(TDF,TDFEnergyArray)
    ! 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
       tdf_unit=io_file_unit()
       open(unit=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 (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. (timing_level>0)) call io_stopwatch('boltzwann_main: calc_props',1)

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

    ! I allocate the arrays for the spectra
    allocate(LocalElCond(6,counts(my_node_id)),stat=ierr)
    if (ierr/=0) call io_error('Error in allocating LocalElCond in boltzwann_main')
    allocate(LocalSigmaS(6,counts(my_node_id)),stat=ierr)
    if (ierr/=0) call io_error('Error in allocating LocalSigmaS in boltzwann_main')
    allocate(LocalSeebeck(9,counts(my_node_id)),stat=ierr)
    if (ierr/=0) call io_error('Error in allocating LocalSeebeck in boltzwann_main')
    allocate(LocalKappa(6,counts(my_node_id)),stat=ierr)
    if (ierr/=0) call io_error('Error in allocating LocalKappa in boltzwann_main')
    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) call io_error('Error in allocating FermiDerivArray in boltzwann_main')

    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)*boltz_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 (boltz_2d_dir_num /= 0) then
          ! Invert only the appropriate 2x2 submatrix
          if (boltz_2d_dir_num == 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 (boltz_2d_dir_num == 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 (boltz_2d_dir_num == 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 (boltz_2d_dir_num == 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 (boltz_2d_dir_num == 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 (boltz_2d_dir_num == 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)*boltz_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)*boltz_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')
    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 boltz_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 * elem_charge_SI**3 / (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 * elem_charge_SI**3 / (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 * elem_charge_SI**3 / (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) call io_error('Error in allocating ElCond in boltzwann_main')
       allocate(SigmaS(6,TempNumPoints,MuNumPoints),stat=ierr)
       if (ierr/=0) call io_error('Error in allocating SigmaS in boltzwann_main')
       allocate(Seebeck(9,TempNumPoints,MuNumPoints),stat=ierr)
       if (ierr/=0) call io_error('Error in allocating Seebeck in boltzwann_main')
       allocate(Kappa(6,TempNumPoints,MuNumPoints),stat=ierr)
       if (ierr/=0) call io_error('Error in allocating Kappa in boltzwann_main')
    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) call io_error('Error in allocating ElCond in boltzwann_main (2)')
       allocate(SigmaS(1,1,1),stat=ierr)
       if (ierr/=0) call io_error('Error in allocating SigmaS in boltzwann_main (2)')
       allocate(Seebeck(1,1,1),stat=ierr)
       if (ierr/=0) call io_error('Error in allocating Seebeck in boltzwann_main (2)')
       allocate(Kappa(1,1,1),stat=ierr)
       if (ierr/=0) call io_error('Error in allocating Kappa in boltzwann_main (2)')
    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(1,1),6*counts(my_node_id),ElCond(1,1,1),6*counts,6*displs)
    call comms_gatherv(LocalSigmaS(1,1),6*counts(my_node_id),SigmaS(1,1,1),6*counts,6*displs)
    call comms_gatherv(LocalSeebeck(1,1),9*counts(my_node_id),Seebeck(1,1,1),9*counts,9*displs)
    call comms_gatherv(LocalKappa(1,1),6*counts(my_node_id),Kappa(1,1,1),6*counts,6*displs)

    if(on_root .and. (timing_level>0)) call io_stopwatch('boltzwann_main: calc_props',2)

    ! Open files and print
    if (on_root) then
       elcond_unit = io_file_unit()
       open(unit=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 (iprint > 1) write(stdout,'(3X,A)') "Electrical conductivity written on the " // trim(seedname)//"_elcond.dat file."

       sigmas_unit = io_file_unit()
       open(unit=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 (iprint > 1) write(stdout,'(3X,A)') &
            "sigma*S (sigma=el. conductivity, S=Seebeck coeff.) written on the " // trim(seedname)//"_sigmas.dat file."

       seebeck_unit =io_file_unit()
       open(unit=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 (iprint > 1)  write(stdout,'(3X,A)') "Seebeck coefficient written on the " // trim(seedname)//"_seebeck.dat file."

       kappa_unit =io_file_unit()
       open(unit=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 (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) call io_error('Error in deallocating TempArray in boltzwann_main')
    deallocate(KTArray,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating KTArray in boltzwann_main')
    deallocate(MuArray,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating MuArray in boltzwann_main')
    deallocate(TDFEnergyArray,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating TDFEnergyArray in boltzwann_main')
    deallocate(TDF,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating TDF in boltzwann_main')
    deallocate(LocalElCond,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating LocalElCond in boltzwann_main')
    deallocate(LocalSigmaS,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating LocalSigmaS in boltzwann_main')
    deallocate(LocalSeebeck,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating LocalSeebeck in boltzwann_main')
    deallocate(LocalKappa,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating LocalKappa in boltzwann_main')

    deallocate(ElCond,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating ElCond in boltzwann_main')
    deallocate(SigmaS,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating SigmaS in boltzwann_main')
    deallocate(Seebeck,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating Seebeck in boltzwann_main')
    deallocate(Kappa,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating Kappa in boltzwann_main')

    deallocate(IntegrandArray,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating IntegrandArray in boltzwann_main')

    if(on_root .and. (timing_level>0)) call io_stopwatch('boltzwann_main',2)

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

  end subroutine boltzwann_main


berry_get_imf_klist berry_get_imfgh_klist berry_get_kubo_k berry_main boltzwann_main calcTDFandDOS check_and_sort_similar_centres clean_ws_translate comms_allreduce comms_allreduce_cmplx comms_allreduce_real comms_array_split comms_barrier comms_bcast comms_bcast_char comms_bcast_cmplx comms_bcast_int comms_bcast_logical comms_bcast_real comms_end comms_gatherv comms_gatherv_real comms_recv comms_recv_char comms_recv_cmplx comms_recv_int comms_recv_logical comms_recv_real comms_reduce comms_reduce_cmplx comms_reduce_int comms_reduce_real comms_scatterv comms_scatterv_int_1 comms_scatterv_int_2 comms_scatterv_int_3 comms_scatterv_real comms_send comms_send_char comms_send_cmplx comms_send_int comms_send_logical comms_send_real comms_setup conv_get_seedname conv_read_chkpt conv_read_chkpt_fmt conv_write_chkpt conv_write_chkpt_fmt dis_extract dis_extract_gamma dis_main dis_proj_froz dis_project dis_windows dos_get_k dos_get_levelspacing dos_main fourier_q_to_R gauss_freq geninterp_main get_AA_R get_BB_R get_CC_R get_FF_R get_HH_R get_module_kmesh get_smearing_index get_SS_R get_win_min group hamiltonian_dealloc hamiltonian_get_hr hamiltonian_setup hamiltonian_wigner_seitz hamiltonian_write_hr hamiltonian_write_rmn hamiltonian_write_tb internal_maxloc internal_set_kmesh internal_write_header io_date io_error io_file_unit io_get_seedname io_print_timings io_stopwatch io_time k_path k_slice kmesh_dealloc kmesh_get kmesh_get_bvectors kmesh_shell_automatic kmesh_shell_fixed kmesh_shell_from_file kmesh_spacing_mesh kmesh_spacing_singleinteger kmesh_supercell_sort kmesh_write master_sort_and_group MinusFermiDerivative my_ICOPY orthogonalize_u overlap_dealloc overlap_project overlap_project_gamma overlap_read overlap_rotate param_dealloc param_get_atoms param_get_block_length param_get_keyword param_get_keyword_block param_get_keyword_kpath param_get_keyword_vector param_get_projections param_get_range_vector param_get_smearing_type param_get_vector_length param_in_file param_lib_set_atoms param_memory_estimate param_postw90_write param_read param_read_chkpt param_uppercase param_write param_write_chkpt param_write_header plot_fermi_surface plot_interpolate_bands plot_main plot_u_matrices plot_wannier print_usage pw90common_fourier_R_to_k pw90common_fourier_R_to_k_new pw90common_fourier_R_to_k_vec pw90common_get_occ pw90common_kmesh_spacing pw90common_wanint_data_dist pw90common_wanint_get_kpoint_file pw90common_wanint_param_dist pw90common_wanint_setup qe_erf qe_erfc R_wz_sc R_wz_sc_equiv script_common script_fermi_lines sitesym_dealloc sitesym_dis_extract_symmetry sitesym_read sitesym_replace_d_matrix_band sitesym_slim_d_matrix_band sitesym_symmetrize_gradient sitesym_symmetrize_rotation sitesym_symmetrize_u_matrix sitesym_symmetrize_zmatrix sort spin_get_moment spin_get_moment_k spin_get_nk symmetrize_ukirr TDF_kpt tran_bulk tran_cut_hr_one_dim tran_dealloc tran_find_integral_signatures tran_get_ht tran_green tran_lcr tran_lcr_2c2_build_ham tran_lcr_2c2_sort tran_main tran_parity_enforce tran_read_htC tran_read_htX tran_read_htXY tran_reduce_hr tran_transfer tran_write_xyz utility_cart_to_frac utility_commutator_diag utility_compar utility_diagonalize utility_frac_to_cart utility_im_tr utility_inv2 utility_inv3 utility_lowercase utility_matmul_diag utility_metric utility_re_tr utility_recip_lattice utility_rotate utility_rotate_diag utility_string_to_coord utility_strip utility_translate_home utility_w0gauss utility_wgauss utility_zgemm wann_calc_projection wann_check_unitarity wann_domega wann_main wann_main_gamma wann_omega wann_omega_gamma wann_phases wann_spread_copy wann_svd_omega_i wann_write_r2mn wann_write_vdw_data wann_write_xyz wannier_run wannier_setup wham_get_D_h wham_get_D_h_a wham_get_deleig_a wham_get_eig_deleig wham_get_eig_UU_HH_JJlist wham_get_JJm_list wham_get_JJp_list wham_get_occ_mat_list wigner_seitz ws_translate_dist ws_write_vec