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:
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.
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 at least 1 entry
allocate (LocalElCond(6, max(1, counts(my_node_id))), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating LocalElCond in boltzwann_main')
allocate (LocalSigmaS(6, max(1, counts(my_node_id))), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating LocalSigmaS in boltzwann_main')
allocate (LocalSeebeck(9, max(1, counts(my_node_id))), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating LocalSeebeck in boltzwann_main')
allocate (LocalKappa(6, max(1, 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, 6*counts(my_node_id), ElCond, 6*counts, 6*displs)
call comms_gatherv(LocalSigmaS, 6*counts(my_node_id), SigmaS, 6*counts, 6*displs)
call comms_gatherv(LocalSeebeck, 9*counts(my_node_id), Seebeck, 9*counts, 9*displs)
call comms_gatherv(LocalKappa, 6*counts(my_node_id), Kappa, 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