Line data Source code
1 : !!****m* ABINIT/m_rta
2 : !! NAME
3 : !! m_rta
4 : !!
5 : !! FUNCTION
6 : !! This module provides objects and procedures to compute transport properties by solving the
7 : !! linearized Boltzmann transport equation (BTE) in the relaxation-time approximation (RTA).
8 : !! or within the Iterative Bolztmann Equation (IBTE).
9 : !! The RTA has two different flavors: Self-energy Relaxation Time Approximation (SERTA)
10 : !! in which the back-scattering term is completely ignored and the Momentum-Relaxation Time Approximation (MRTA)
11 : !! in which backscattering is partly accounter for by multiplying the e-ph self-energy integrand function
12 : !! by the efficiency factor alpha that depends on the incoming/outgoing electron group velocity.
13 : !! The implementation assumes e-ph scattering although additional scattering mechanisms (e.g. ionized impurities)
14 : !! can be easily included once an appropriate model is added to the ab-initio e-ph scattering rates.
15 : !!
16 : !! COPYRIGHT
17 : !! Copyright (C) 2008-2026 ABINIT group (HM, MG)
18 : !! This file is distributed under the terms of the
19 : !! GNU General Public License, see ~abinit/COPYING
20 : !! or http://www.gnu.org/copyleft/gpl.txt .
21 : !!
22 : !! NOTES
23 : !! Dimensional analysis for conductivity (sigma), mobility (mu).
24 : !!
25 : !! [sigma] = Siemens/m with S = Ampere/Volt = Ohm^-1
26 : !! [mu] = S L^2 Q
27 : !!
28 : !! SOURCE
29 :
30 : #if defined HAVE_CONFIG_H
31 : #include "config.h"
32 : #endif
33 :
34 : #include "abi_common.h"
35 :
36 : module m_rta
37 :
38 : use defs_basis
39 : use m_abicore
40 : use m_xmpi
41 : use m_errors
42 : use m_copy
43 : use m_nctk
44 : use m_wfk
45 : use m_ephtk
46 : use m_sigmaph
47 : use m_dtset
48 : use m_dtfil
49 : use m_krank
50 : use netcdf
51 :
52 : use defs_datatypes, only : pseudopotential_type
53 : use m_ebands, only : ebands_t, gaps_t, edos_t, klinterp_t, klinterp_new
54 : use m_io_tools, only : flush_unit, open_file
55 : use m_time, only : cwtime, cwtime_report
56 : use m_crystal, only : crystal_t
57 : use m_numeric_tools, only : bisect, simpson_int, safe_div, arth
58 : use m_fstrings, only : strcat, sjoin, itoa, ltoa, stoa, ftoa, yesno
59 : use m_kpts, only : kpts_timrev_from_kptopt, kpts_map
60 : use m_occ, only : occ_fd, occ_dfde
61 : use m_pawtab, only : pawtab_type
62 : use m_ddk, only : ddkstore_t
63 :
64 : implicit none
65 :
66 : private
67 : !!****
68 :
69 : public :: rta_driver ! Compute transport properties within the RTA (SERTA and MRTA)
70 : public :: ibte_driver ! Compute transport properties within the IBTE.
71 : !!****
72 :
73 : !----------------------------------------------------------------------
74 :
75 : !!****t* m_rta/rta_t
76 : !! NAME
77 : !! rta_t
78 : !!
79 : !! FUNCTION
80 : !! Container for transport quantities in the RTA
81 : !!
82 : !! SOURCE
83 :
84 : type,public :: rta_t
85 :
86 : integer :: nsppol
87 : ! Number of independent spin polarizations.
88 :
89 : integer :: nspinor
90 : ! Number of spinor components.
91 :
92 : integer :: nkcalc
93 : ! Number of computed k-points i.e. k-points inside the sigma_erange energy window.
94 :
95 : integer :: ntemp
96 : ! Number of temperatures.
97 :
98 : integer :: nw
99 : ! Number of energies (chemical potentials) at which transport quantities are computed
100 : ! Same number of energies used in DOS.
101 :
102 : integer :: bmin, bmax, bsize
103 : ! Only bands between bmin and bmax are considered in the integrals
104 : ! as we don't compute linewidths for all states.
105 : ! bmin = minval(%bstart_ks); bmax = maxval(%bstop_ks)
106 : ! bisze = bmax - bmin + 1
107 :
108 : integer :: nrta
109 : ! Number of relaxation-time approximations used (1 for SERTA, 2 for MRTA)
110 :
111 : real(dp) :: eph_extrael
112 : ! Extra electrons per unit cell used to compute SERTA lifetimes in sigmaph.
113 :
114 : real(dp) :: eph_fermie
115 : ! Fermi level specified in the input file when computing the SIGEPH file.
116 :
117 : real(dp) :: transport_extrael
118 : ! Extra electrons per unit cell specified in the input file when computing the SIGEPH file.
119 :
120 : real(dp) :: transport_fermie
121 : ! Fermi level specified in the input file when computing the SIGEPH file.
122 :
123 : logical :: assume_gap
124 : ! True if we are dealing with a semiconductor.
125 : ! This parameter is initialized from the value of sigma_erange provided by the user.
126 :
127 : integer,allocatable :: bstart_ks(:,:)
128 : ! bstart_ks(nkcalc, nsppol)
129 : ! Initial KS band index included in self-energy matrix elements for each k-point in kcalc.
130 : ! Depends on spin because all degenerate states should be included when symmetries are used.
131 :
132 : integer,allocatable :: bstop_ks(:,:)
133 : ! bstop_ks(nkcalc, nsppol)
134 :
135 : integer,allocatable :: nbcalc_ks(:,:)
136 : ! nbcalc_ks(nkcalc, nsppol)
137 : ! Number of bands included in self-energy matrix elements for each k-point in kcalc.
138 : ! Depends on spin because all degenerate states should be included when symmetries are used.
139 :
140 : integer,allocatable :: kcalc2ibz(:,:)
141 : !kcalc2ibz(nkcalc, 6))
142 : ! Mapping ikcalc --> IBZ as reported by listkk.
143 :
144 : integer,allocatable :: kcalc2ebands(:,:)
145 : ! (6, nkcalc)
146 : ! Mapping ikcalc --> ebands IBZ
147 : ! Note that this array is not necessarily equal to kcalc2ibz computed in sigmaph
148 : ! because we may have used sigma_ngkpt to downsample the initial nkpt mesh.
149 : ! This array is computed in get_ebands and is equal to kcalc2ibz if sigma_nkpt == ngkpt
150 :
151 : real(dp),allocatable :: kTmesh(:)
152 : ! (%ntemp)
153 : ! List of k * T temperatures at which to compute the transport
154 :
155 : real(dp),allocatable :: eph_mu_e(:)
156 : ! (%ntemp)
157 : ! Chemical potential at this carrier concentration and temperature from sigeph (lifetime)
158 :
159 : real(dp),allocatable :: transport_mu_e(:)
160 : ! (%ntemp)
161 : ! Chemical potential at this carrier concentration and temperature
162 :
163 : real(dp),allocatable :: eminmax_spin(:,:)
164 : ! (2, %nsppol))
165 : ! min/Max energy of the original ebands object
166 :
167 : real(dp),allocatable :: linewidths(:,:,:,:,:)
168 : ! (ntemp, bmin:bmax, nkpt, nsppol, nrta)
169 : ! Linewidth in the IBZ computed in the SERTA/MRTA.
170 : ! Non-zero only for the kcalc k-points.
171 :
172 : real(dp),allocatable :: vbks(:,:,:,:)
173 : ! (3, bmin:bmax, nkpt, nsppol))
174 : ! band velocity in Cartesian coordinates in the IBZ
175 : ! Non-zero only for the kcalc k-points.
176 :
177 : type(gaps_t) :: gaps
178 : ! gaps of original ebands object. Only if assume_gap
179 :
180 : type(ebands_t) :: ebands
181 : ! bandstructure object used to compute the transport properties
182 : ! Allocate using only the relevant bands for transport
183 : ! including valence states to allow to compute different doping
184 :
185 : type(edos_t) :: edos
186 : ! electronic density of states
187 : ! edos%mesh is the mesh used for vv_dos, vvtau_dos and tau_dos
188 : ! (%nw)
189 :
190 : real(dp),allocatable :: tau_dos(:,:,:,:)
191 : ! tau(e) (isotropic average for tau_nk for SERTA and MRTA.
192 : ! (nw, ntemp, nsppol, nrta)
193 :
194 : real(dp),allocatable :: vv_dos(:,:,:,:)
195 : ! (v x v) DOS
196 : ! (nw, 3, 3, nsppol)
197 :
198 : real(dp),allocatable :: vvtau_dos(:,:,:,:,:,:)
199 : ! (v x v * tau) DOS
200 : ! (nw, 3, 3, ntemp, nsppol, nrta)
201 :
202 : real(dp),allocatable :: n_ehst(:,:,:)
203 : ! (2, %nsppol, %ntemp)
204 : ! Number of electrons (e) and holes (h) per unit cell
205 : ! The first dimension is for electrons/holes.
206 : ! If nsppol == 2, the second dimension is the number of e/h for spin else the total number of e/h summed over spins.
207 :
208 : real(dp),allocatable :: l0(:,:,:,:,:,:), l1(:,:,:,:,:,:), l2(:,:,:,:,:,:)
209 : ! (3, 3, nw, nsppol, ntemp, nrta)
210 : ! Onsager coefficients in Cartesian coordinates
211 :
212 : real(dp),allocatable :: l11_mu(:,:,:,:,:), l12_mu(:,:,:,:,:), l22_mu(:,:,:,:,:)
213 : ! (3, 3, nsppol, ntemp, nrta)
214 : ! Onsager coefficients in Cartesian coordinates, calculated at the correct mu (at the exact temperature)
215 :
216 : real(dp),allocatable :: sigma(:,:,:,:,:,:)
217 : real(dp),allocatable :: seebeck(:,:,:,:,:,:)
218 : real(dp),allocatable :: kappa(:,:,:,:,:,:)
219 : real(dp),allocatable :: pi(:,:,:,:,:,:)
220 : real(dp),allocatable :: zte(:,:,:,:,:,:)
221 : ! (3, 3, nw, nsppol, ntemp, nrta)
222 : ! Transport coefficients in Cartesian coordinates
223 :
224 : real(dp),allocatable :: mobility(:,:,:,:,:,:,:)
225 : ! Mobility
226 : ! (3, 3, nw, %ntemp, 2, %nsppol, nrta)
227 : ! 5-th index is for e-h
228 :
229 : real(dp),allocatable :: conductivity(:,:,:,:,:)
230 : ! Conductivity at the Fermi level
231 : ! (3, 3, %ntemp, %nsppol, nrta)
232 :
233 : real(dp),allocatable :: resistivity(:,:,:,:)
234 : ! (3, 3, %ntemp, nrta)
235 :
236 : real(dp),allocatable :: n(:,:,:)
237 : ! (nw, ntemp, 2) carrier density for e/h (n/cm^3)
238 :
239 : real(dp),allocatable :: mobility_mu(:,:,:,:,:,:)
240 : ! (3, 3, 2, nsppol, ntemp, nrta)
241 : ! mobility for electrons and holes (third dimension) at transport_mu_e(ntemp)
242 : ! Third dimension is for electron/hole
243 :
244 : real(dp),allocatable :: conductivity_mu(:,:,:,:,:,:)
245 : ! (3, 3, 2, nsppol, ntemp, nrta)
246 : ! Conductivity in Siemens * cm-1
247 : ! computed by summing over k-points rather that by performing an energy integration).
248 :
249 : real(dp),allocatable :: resistivity_mu(:,:,:,:)
250 : ! (3, 3, ntemp, nrta)
251 : ! Resistivity obtained by inverting conductivity after summing of e- and holes and spins
252 : ! computed by summing over k-points rather that by performing an energy integration).
253 :
254 :
255 :
256 : contains
257 :
258 : procedure :: compute_rta
259 : procedure :: compute_rta_mobility
260 : procedure :: print_rta_txt_files
261 : procedure :: write_tensor
262 : procedure :: free => rta_free
263 : procedure :: rta_ncwrite
264 :
265 : end type rta_t
266 : !!***
267 :
268 : !----------------------------------------------------------------------
269 :
270 : contains !=====================================================
271 : !!***
272 :
273 : !----------------------------------------------------------------------
274 :
275 : !!****f* m_rta/rta_driver
276 : !! NAME
277 : !! rta_driver
278 : !!
279 : !! FUNCTION
280 : !! Driver to compute transport properties within the RTA.
281 : !!
282 : !! INPUTS
283 : !! dtfil<datafiles_type>=variables related to files.
284 : !! ngfftc(18)=Coarse FFT mesh.
285 : !! dtset<dataset_type>=All input variables for this dataset.
286 : !! ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...)
287 : !! cryst<crystal_t>=Crystalline structure
288 : !! pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
289 : !! psps<pseudopotential_type>=Variables related to pseudopotentials.
290 : !! comm=MPI communicator.
291 : !!
292 : !! SOURCE
293 :
294 12 : subroutine rta_driver(dtfil, ngfftc, dtset, ebands, cryst, pawtab, psps, comm)
295 :
296 : !Arguments ------------------------------------
297 : !scalars
298 : integer, intent(in) :: comm
299 : type(datafiles_type),intent(in) :: dtfil
300 : type(dataset_type),intent(in) :: dtset
301 : type(crystal_t),intent(in) :: cryst
302 : type(ebands_t),intent(in) :: ebands
303 : type(pseudopotential_type),intent(in) :: psps
304 : !arrays
305 : integer,intent(in) :: ngfftc(18)
306 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
307 :
308 : !Local variables ------------------------------
309 : type(rta_t) :: rta
310 : !arrays
311 : integer :: units(2)
312 :
313 : ! *************************************************************************
314 :
315 36 : units = [std_out, ab_out]
316 12 : call wrtout(units, ch10//' Entering transport RTA computation driver.')
317 12 : call wrtout(units, sjoin("- Reading carrier lifetimes from:", dtfil%filsigephin), newlines=1, do_flush=.True.)
318 :
319 : ! Initialize RTA object
320 12 : rta = rta_new(dtset, dtfil, ngfftc, cryst, ebands, pawtab, psps, comm)
321 :
322 : ! Compute RTA transport quantities
323 12 : call rta%compute_rta(cryst, dtset, dtfil, comm)
324 :
325 12 : call rta%free()
326 :
327 12 : end subroutine rta_driver
328 : !!***
329 :
330 : !----------------------------------------------------------------------
331 :
332 : !!****f* m_rta/rta_new
333 : !! NAME
334 : !! rta_new
335 : !!
336 : !! FUNCTION
337 : !! Build object to compute RTA transport quantities.
338 : !!
339 : !! INPUTS
340 : !! dtset<dataset_type>=All input variables for this dataset.
341 : !! dtfil<datafiles_type>=variables related to files.
342 : !! cryst<crystal_t>=Crystalline structure
343 : !! ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...)
344 : !! comm=MPI communicator.
345 : !!
346 : !! SOURCE
347 :
348 13 : type(rta_t) function rta_new(dtset, dtfil, ngfftc, cryst, ebands, pawtab, psps, comm) result (new)
349 :
350 : !Arguments -------------------------------------
351 : integer, intent(in) :: comm
352 : type(dataset_type),intent(in) :: dtset
353 : type(datafiles_type),intent(in) :: dtfil
354 : type(crystal_t),intent(in) :: cryst
355 : type(ebands_t),intent(in) :: ebands
356 : type(pseudopotential_type),intent(in) :: psps
357 : !arrays
358 : integer,intent(in) :: ngfftc(18)
359 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
360 :
361 : !Local variables ------------------------------
362 : integer,parameter :: sppoldbl1 = 1, master = 0
363 : integer :: ierr, spin, nprocs, my_rank, ik_ibz, ib, irta, itemp, ndat, nsppol, idat, mband, ikpt
364 : real(dp) :: cpu, wall, gflops
365 : character(len=500) :: msg
366 : character(len=fnlen) :: wfk_fname_dense
367 13 : type(ebands_t) :: tmp_ebands, ebands_dense
368 13 : type(klinterp_t) :: klinterp
369 13 : type(ddkstore_t) :: ds
370 13 : type(sigmaph_t) :: sigmaph
371 13 : type(krank_t) :: krank
372 : !arrays
373 : integer :: kptrlatt(3,3), units(2), sigma_ngkpt(3)
374 13 : integer,allocatable :: indkk(:,:)
375 : real(dp) :: extrael_fermie(2), sigma_erange(2)
376 13 : real(dp),allocatable :: values_bksd(:,:,:,:), vals_bsd(:,:,:), tmp_array4(:,:,:,:), tmp_array5(:,:,:,:,:)
377 : !************************************************************************
378 :
379 13 : call cwtime(cpu, wall, gflops, "start")
380 39 : units = [std_out, ab_out]
381 :
382 13 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
383 :
384 : ! Use sigma_erange to understand if we are dealing with a metal or a semiconductor.
385 15 : new%assume_gap = (.not. all(dtset%sigma_erange < zero) .or. dtset%gw_qprange /= 0)
386 :
387 : ! Read data from SIGEPH file.
388 : sigmaph = sigmaph_read(dtfil%filsigephin, dtset, xmpi_comm_self, msg, ierr, keep_open=.true., &
389 13 : extrael_fermie=extrael_fermie, sigma_ngkpt=sigma_ngkpt, sigma_erange=sigma_erange)
390 13 : ABI_CHECK(ierr == 0, msg)
391 :
392 : !if (any(sigma_erange /= zero)) then
393 : ! ABI_CHECK(all(dtset%sigma_erange /= zero), "sigma_erange is required in input with a value compatible with SIGEPH.nc")
394 : ! ! Make sure that the two values are consistent
395 : ! ! Cannot switch to metallic case if SigmaPH file was produced assuming gapped-system
396 : ! if (.not. (all(dtset%sigma_erange < zero) .eqv. all(sigma_erange < zero))) then
397 : ! ABI_ERROR("The values of sigma_erange from input and SIGEPH are not compatible")
398 : ! end if
399 : !end if
400 :
401 : ! How many RTA approximations have we computed in sigmaph? (SERTA, MRTA ...?)
402 13 : new%nrta = 2; if (sigmaph%mrta == 0) new%nrta = 1
403 :
404 : ! Copy important arrays from sigmaph file.
405 : ! Allocate temperature arrays (use same values as the ones used in the SIGEPH calculation).
406 13 : new%ntemp = sigmaph%ntemp
407 13 : call alloc_copy(sigmaph%kTmesh, new%kTmesh)
408 :
409 13 : new%nkcalc = sigmaph%nkcalc
410 13 : call alloc_copy(sigmaph%bstart_ks, new%bstart_ks)
411 13 : call alloc_copy(sigmaph%bstop_ks, new%bstop_ks)
412 13 : call alloc_copy(sigmaph%nbcalc_ks, new%nbcalc_ks)
413 13 : call alloc_copy(sigmaph%kcalc2ibz, new%kcalc2ibz)
414 :
415 141 : new%bmin = minval(sigmaph%bstart_ks); new%bmax = maxval(sigmaph%bstop_ks)
416 : !new%bmin = 1; new%bmax = ebands%mband ! This for debugging purposes, results should not change
417 13 : new%bsize = new%bmax - new%bmin + 1
418 :
419 13 : new%nsppol = ebands%nsppol; new%nspinor = ebands%nspinor
420 13 : nsppol = new%nsppol
421 :
422 39 : ABI_MALLOC(new%eminmax_spin, (2, nsppol))
423 65 : new%eminmax_spin = ebands%get_minmax("eig")
424 :
425 13 : if (new%assume_gap) then
426 : ! Get gaps
427 12 : new%gaps = ebands%get_gaps(ierr)
428 12 : if (ierr /= 0) then
429 0 : do spin=1, nsppol
430 0 : ABI_WARNING(trim(new%gaps%errmsg_spin(spin)))
431 0 : new%gaps%vb_max(spin) = ebands%fermie - 1 * eV_Ha
432 0 : new%gaps%cb_min(spin) = ebands%fermie + 1 * eV_Ha
433 : end do
434 : !ABI_ERROR("ebands_get_gaps returned non-zero exit status. See above warning messages...")
435 0 : ABI_WARNING("ebands_get_gaps returned non-zero exit status. See above warning messages...")
436 : end if
437 12 : if (my_rank == master) call new%gaps%print(units)
438 : end if
439 :
440 : ! =================================================
441 : ! Read lifetimes and new%ebands from SIGEPH.nc file
442 : ! =================================================
443 : ! After this point we have:
444 : !
445 : ! vbks(3, bmin:bmax, nkpt, nsppol)
446 : ! linewidths(self%ntemp, bmin:bmax, nkpt, nsppol, 2)
447 : !
448 46 : if (any(dtset%sigma_ngkpt /= 0)) then
449 : ! If integrals are computed with the sigma_ngkpt k-mesh, we need to downsample ebands.
450 : !call wrtout(units, sjoin(" SIGMAPH file used sigma_ngkpt:", ltoa(sigma_ngkpt)))
451 2 : call wrtout(units, sjoin(" Computing integrals with downsampled sigma_ngkpt:", ltoa(dtset%sigma_ngkpt)))
452 2 : kptrlatt = 0
453 2 : kptrlatt(1,1) = dtset%sigma_ngkpt(1); kptrlatt(2,2) = dtset%sigma_ngkpt(2); kptrlatt(3,3) = dtset%sigma_ngkpt(3)
454 :
455 2 : tmp_ebands = ebands%downsample(cryst, kptrlatt, dtset%sigma_nshiftk, dtset%sigma_shiftk)
456 : new%ebands = sigmaph%get_ebands(cryst, tmp_ebands, [new%bmin, new%bmax], &
457 6 : new%kcalc2ebands, new%linewidths, new%vbks, xmpi_comm_self)
458 2 : call tmp_ebands%free()
459 : else
460 : !call wrtout(units, sjoin(" Computing integrals with SIGEPH k-mesh:", ebands_kmesh2str(ebands))
461 : new%ebands = sigmaph%get_ebands(cryst, ebands, [new%bmin, new%bmax], &
462 33 : new%kcalc2ebands, new%linewidths, new%vbks, xmpi_comm_self)
463 143 : kptrlatt = new%ebands%kptrlatt
464 : end if
465 :
466 : !print *, "linewidth_serta", maxval(abs(new%linewidths(:,:,:,:,1)))
467 : !print *, "linewidth_mrta", maxval(abs(new%linewidths(:,:,:,:,2)))
468 : !print *, "max velocities", maxval(abs(new%vbks))
469 :
470 : if ( &
471 13 : dtset%useria == 888 .and. &
472 : (dtset%getwfkfine /= 0 .or. dtset%irdwfkfine /= 0 .or. dtset%getwfkfine_filepath /= ABI_NOFILE)) then
473 :
474 : ! In principle only getwfkfine_filepath is used here
475 0 : wfk_fname_dense = trim(dtfil%fnameabi_wfkfine)
476 0 : ABI_CHECK(nctk_try_fort_or_ncfile(wfk_fname_dense, msg) == 0, msg)
477 :
478 0 : call wrtout(units, " EPH double grid interpolation: will read energies from: "//trim(wfk_fname_dense), newlines=1)
479 0 : mband = new%ebands%mband
480 :
481 : !ebands_dense = wfk_read_ebands(wfk_fname_dense, comm)
482 :
483 0 : tmp_ebands = wfk_read_ebands(wfk_fname_dense, comm)
484 0 : ebands_dense = tmp_ebands%chop(1, mband)
485 0 : call tmp_ebands%free()
486 0 : if (my_rank == master) then
487 0 : write(std_out, *)" Using kptrlatt: ", ebands_dense%kptrlatt
488 0 : write(std_out, *)" shiftk: ", ebands_dense%shiftk
489 : end if
490 0 : ABI_CHECK_IEQ(mband, ebands_dense%mband, "Inconsistent number of bands for the fine and dense grid:")
491 :
492 : ! Compute v_{nk} on the dense grid in Cartesian coordinates.
493 : ! vdiago(3, bmin:bmax, nkpt, nsppol)
494 : ! NB: We select bands in [bmin:bmax] but all k-points in the IBZ are computed!
495 0 : ds%only_diago = .True.; ds%bmin = new%bmin; ds%bmax = new%bmax; ds%mode = "cart"
496 0 : call ds%compute_ddk(wfk_fname_dense, "", dtset, psps, pawtab, ngfftc, comm)
497 :
498 : ! Transfer data to new%vbks
499 0 : ABI_MOVE_ALLOC(ds%vdiago, new%vbks)
500 0 : call ds%free()
501 : !print *, "vbks:", new%vbks
502 :
503 : ! Linear interpolation in k-space of the linewidths from input SIGEPH to the dense IBZ provided by fine WFK file.
504 : ! First of all transfer linewidths to values_bksd to prepare call to klinterp_new.
505 0 : ndat = new%ntemp * new%nrta
506 0 : ABI_MALLOC(values_bksd, (new%bmin:new%bmax, new%ebands%nkpt, nsppol, ndat))
507 :
508 0 : do irta=1,new%nrta
509 0 : do spin=1,nsppol
510 0 : do ik_ibz=1,new%ebands%nkpt
511 0 : do ib=new%bmin,new%bmax
512 0 : do itemp=1,new%ntemp
513 0 : idat = itemp + new%ntemp * (irta - 1)
514 0 : values_bksd(ib, ik_ibz, spin, idat) = new%linewidths(itemp, ib, ik_ibz, spin, irta) ! - t(e)
515 : end do
516 : end do
517 : end do
518 : end do
519 : end do
520 :
521 : ! Build linear interpolator for linewidths (use only bsize bands)
522 : klinterp = klinterp_new(cryst, new%ebands%kptrlatt, new%ebands%nshiftk, new%ebands%shiftk, new%ebands%kptopt, &
523 0 : new%ebands%kptns, new%bsize, new%ebands%nkpt, nsppol, ndat, values_bksd, comm)
524 0 : ABI_FREE(values_bksd)
525 :
526 : ! HERE we re-malloc new%ebands and %linewidths on the fine k-mesh.
527 : ! The call must be executed here, once klinterp has been built.
528 : ! After this point we can use new%ebands to allocate stuff.
529 0 : call ebands_dense%move_alloc(new%ebands)
530 :
531 : ! Unlinke the ebands stored in SIGEPH, the eigens read from WFK_FINE have not been
532 : ! shifted with the scissors operator or updated according to extrael_fermie so do it now.
533 0 : call ephtk_update_ebands(dtset, new%ebands, "GS energies read from WFK_FINE")
534 :
535 : ! And now interpolate linewidths on the fine k-mesh
536 : ! Note: k-points that close to the edge of the pocket may get zero linewidths
537 : ! One may fix the problem by using lw(e).
538 0 : ABI_REMALLOC(new%linewidths, (new%ntemp, new%bmin:new%bmax, new%ebands%nkpt, nsppol, new%nrta))
539 0 : ABI_MALLOC(vals_bsd, (new%bmin:new%bmax, nsppol, ndat))
540 :
541 0 : ierr = 0
542 0 : do ik_ibz=1,new%ebands%nkpt
543 :
544 0 : call klinterp%eval_bsd(new%ebands%kptns(:, ik_ibz), vals_bsd)
545 : !vals_bsd = vals_bsd + lw(e)
546 :
547 0 : if (any(vals_bsd < zero)) then
548 0 : ierr = ierr + 1
549 0 : where (vals_bsd < zero) vals_bsd = zero
550 : end if
551 :
552 : ! Transfer data.
553 0 : do spin=1,nsppol
554 0 : do irta=1,new%nrta
555 0 : do itemp=1,new%ntemp
556 0 : idat = itemp + new%ntemp * (irta - 1)
557 0 : do ib=new%bmin,new%bmax
558 0 : new%linewidths(itemp, ib, ik_ibz, spin, irta) = vals_bsd(ib, spin, idat)
559 : end do
560 : end do
561 : end do
562 : end do
563 :
564 : end do ! ik_ibz
565 :
566 0 : if (ierr /= 0) then
567 : ! This should never happen for linear interpolation.
568 0 : ABI_WARNING(sjoin("Linear interpolation produced:", itoa(ierr), " k-points with negative linewidths"))
569 : end if
570 :
571 0 : ABI_FREE(vals_bsd)
572 0 : call klinterp%free()
573 : end if
574 :
575 : ! FIXME: I think transport_ngkpt is buggy, wrong ne(T), weird zeros if MRTA ...
576 : ! Do we really need this option? Can't we replace it with sigma_ngkpt and eph_task 7?
577 :
578 49 : if (any(dtset%transport_ngkpt /= 0)) then
579 : ! Perform further downsampling (useful for debugging purposes)
580 1 : call wrtout(units, " Downsampling the k-mesh before computing transport:")
581 1 : call wrtout(units, sjoin(" Using transport_ngkpt: ", ltoa(dtset%transport_ngkpt)))
582 1 : kptrlatt = 0
583 1 : kptrlatt(1, 1) = dtset%transport_ngkpt(1)
584 1 : kptrlatt(2, 2) = dtset%transport_ngkpt(2)
585 1 : kptrlatt(3, 3) = dtset%transport_ngkpt(3)
586 1 : tmp_ebands = new%ebands%downsample(cryst, kptrlatt, 1, [zero, zero, zero])
587 :
588 : ! Map the points of the downsampled bands to dense ebands
589 3 : ABI_MALLOC(indkk, (6, tmp_ebands%nkpt))
590 :
591 1 : call krank%from_kptrlatt(new%ebands%nkpt, new%ebands%kptns, new%ebands%kptrlatt, compute_invrank=.False.)
592 :
593 1 : if (kpts_map("symrec", ebands%kptopt, cryst, krank, tmp_ebands%nkpt, tmp_ebands%kptns, indkk) /= 0) then
594 : write(msg, '(3a)' ) &
595 0 : "Error while downsampling ebands in the transport driver",ch10, &
596 0 : "The k-point could not be generated from a symmetrical one."
597 0 : ABI_ERROR(msg)
598 : end if
599 :
600 1 : call krank%free()
601 :
602 : ! Downsampling linewidths and velocities.
603 1 : ABI_MOVE_ALLOC(new%linewidths, tmp_array5)
604 7 : ABI_REMALLOC(new%linewidths, (new%ntemp, new%bmin:new%bmax, tmp_ebands%nkpt, nsppol, new%nrta))
605 73 : do ikpt=1,tmp_ebands%nkpt
606 1369 : new%linewidths(:,:,ikpt,:,:) = tmp_array5(:,:,indkk(1, ikpt),:,:)
607 : end do
608 1 : ABI_FREE(tmp_array5)
609 :
610 1 : ABI_MOVE_ALLOC(new%vbks, tmp_array4)
611 5 : ABI_REMALLOC(new%vbks, (3, new%bmin:new%bmax, tmp_ebands%nkpt, nsppol))
612 73 : do ikpt=1,tmp_ebands%nkpt
613 433 : new%vbks(:,:,ikpt,:) = tmp_array4(:,:,indkk(1, ikpt),:)
614 : end do
615 1 : ABI_FREE(tmp_array4)
616 :
617 : !print *, "after downsampling linewidths"
618 : !print *, "linewidth_serta", maxval(abs(new%linewidths(:,:,:,:,1)))
619 : !print *, "linewidth_mrta", maxval(abs(new%linewidths(:,:,:,:,2)))
620 :
621 1 : ABI_FREE(indkk)
622 1 : call tmp_ebands%move_alloc(new%ebands)
623 : end if
624 :
625 : ! Same doping case as in sigmaph file.
626 39 : ABI_MALLOC(new%eph_mu_e, (new%ntemp))
627 26 : ABI_MALLOC(new%transport_mu_e, (new%ntemp))
628 :
629 13 : new%eph_extrael = extrael_fermie(1)
630 13 : new%eph_fermie = extrael_fermie(2)
631 13 : new%transport_fermie = dtset%eph_fermie
632 13 : new%transport_extrael = dtset%eph_extrael
633 88 : new%eph_mu_e = sigmaph%mu_e
634 88 : new%transport_mu_e = sigmaph%mu_e
635 :
636 13 : if (new%transport_fermie /= zero) new%transport_mu_e = new%transport_fermie
637 :
638 13 : if (new%transport_fermie == zero .and. new%transport_extrael /= new%eph_extrael) then
639 :
640 : if (new%transport_extrael /= new%eph_extrael) then
641 : write(msg,'(2(a,e18.8),3a)') &
642 0 : ' extrael from SIGEPH: ',new%transport_extrael, ' and input file: ',new%eph_extrael, "differ", ch10, &
643 0 : ' Will recompute the chemical potential'
644 0 : call wrtout(std_out, msg)
645 : end if
646 :
647 : ! Compute Fermi level for different T values.
648 0 : call ebands%get_muT_with_fd(new%ntemp, new%kTmesh, dtset%spinmagntarget, dtset%prtvol, new%transport_mu_e, comm)
649 : end if
650 :
651 : ! TODO: Implement possible change of sigma_erange, useful for convergence studies
652 : ! 1) Run sigmaph with relatively large sigma_erange.
653 : ! 2) Decrease energy window in the transport part to analyze the behaviour of transport tensors.
654 :
655 : ! sigmaph is not needed anymore. Free it.
656 13 : sigmaph%ncid = nctk_noid
657 13 : call sigmaph%free()
658 :
659 13 : call cwtime_report(" rta_new", cpu, wall, gflops)
660 :
661 26 : end function rta_new
662 : !!***
663 :
664 : !----------------------------------------------------------------------
665 :
666 : !!****f* m_rta/compute_rta
667 : !! NAME
668 : !! compute_rta
669 : !!
670 : !! FUNCTION
671 : !!
672 : !! INPUTS
673 : !! cryst<crystal_t>=Crystalline structure
674 : !! dtset<dataset_type>=All input variables for this dataset.
675 : !! dtfil<datafiles_type>=variables related to files.
676 : !! comm=MPI communicator.
677 : !!
678 : !! SOURCE
679 :
680 13 : subroutine compute_rta(self, cryst, dtset, dtfil, comm)
681 :
682 : !Arguments ------------------------------------
683 : integer,intent(in) :: comm
684 : class(rta_t),intent(inout) :: self
685 : type(dataset_type),intent(in) :: dtset
686 : type(datafiles_type),intent(in) :: dtfil
687 : type(crystal_t),intent(in) :: cryst
688 :
689 : !Local variables ------------------------------
690 : integer,parameter :: nvecs0 = 0, master = 0
691 : integer :: nsppol, nkibz, ib, ik_ibz, iw, spin, ii, jj, itemp, irta, itens_, iscal, cnt
692 : integer :: ntens, edos_intmeth, ifermi, iel, nvals, my_rank
693 : integer :: ncid
694 : !character(len=500) :: msg
695 : character(len=fnlen) :: path
696 : real(dp) :: emin, emax, edos_broad, edos_step, max_occ, kT, Tkelv, linewidth, fact0, cpu, wall, gflops
697 : !arrays
698 : integer :: units(2)
699 : real(dp) :: vr(3), dummy_vecs(1,1,1,1,1), work_33(3,3), S_33(3,3), mat33(3,3)
700 13 : real(dp),allocatable :: vv_tens(:,:,:,:,:,:,:), out_valsdos(:,:,:,:), dummy_dosvecs(:,:,:,:,:)
701 13 : real(dp),allocatable :: out_tensdos(:,:,:,:,:,:), tau_vals(:,:,:,:,:), l0inv_33nw(:,:,:)
702 : !************************************************************************
703 :
704 13 : call cwtime(cpu, wall, gflops, "start")
705 13 : my_rank = xmpi_comm_rank(comm)
706 39 : units = [std_out, ab_out]
707 :
708 : ! Basic dimensions
709 13 : nsppol = self%ebands%nsppol; nkibz = self%ebands%nkpt
710 :
711 : ! Allocate v x v tensors with and without the lifetimes. Eq 8 of [[cite:Madsen2018]]
712 : ! The total number of tensorial entries is ntens and accounts for nrta
713 : ! Remember that we haven't computed all the k-points in the IBZ hence we can have zero linewidths
714 : ! or very small values when the states are at the band edge so we use safe_dif to avoid SIGFPE.
715 : ! Also, note how we store only the states in the energy window.
716 :
717 13 : nvals = self%ntemp * self%nrta
718 56295 : ABI_CALLOC(tau_vals, (self%ntemp, self%nrta, self%bmin:self%bmax, nkibz, nsppol))
719 :
720 13 : ntens = (1 + self%ntemp) * self%nrta
721 651318 : ABI_CALLOC(vv_tens, (3, 3, 1 + self%ntemp, self%nrta, self%bmin:self%bmax, nkibz, nsppol))
722 :
723 : cnt = 0
724 26 : do spin=1,nsppol
725 3582 : do ik_ibz=1,nkibz
726 : !cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! MPI parallelism.
727 7230 : do ib=self%bmin,self%bmax
728 :
729 14644 : vr(:) = self%vbks(:, ib, ik_ibz, spin)
730 : ! Store outer product (v_bks x v_bks) in vv_tens. This part does not depend on T and irta.
731 14644 : do ii=1,3
732 47593 : do jj=1,3
733 109830 : vv_tens(ii, jj, 1, 1:self%nrta, ib, ik_ibz, spin) = vr(ii) * vr(jj)
734 : end do
735 : end do
736 :
737 : ! Multiply by the lifetime (SERTA and MRTA)
738 14539 : do irta=1,self%nrta
739 52635 : do itemp=1,self%ntemp
740 41652 : linewidth = self%linewidths(itemp, ib, ik_ibz, spin, irta)
741 : call safe_div(vv_tens(:,:, 1, irta, ib, ik_ibz, spin), two * linewidth, zero, &
742 541476 : vv_tens(:,:, 1 + itemp, irta, ib, ik_ibz, spin))
743 48974 : call safe_div(one, two * linewidth, zero, tau_vals(itemp, irta, ib, ik_ibz, spin))
744 : end do
745 : end do
746 :
747 : end do
748 : end do
749 : end do
750 :
751 : !call xmpi_sum(vv_tens, comm, ierr)
752 : !call xmpi_sum(tau_vals, comm, ierr)
753 13 : call cwtime_report(" compute_rta_loop1", cpu, wall, gflops)
754 :
755 : ! Compute DOS and VV_DOS and VV_TAU_DOS
756 : ! Define integration method and mesh step.
757 13 : edos_intmeth = 2; if (dtset%prtdos /= 0) edos_intmeth = dtset%prtdos
758 13 : edos_step = dtset%dosdeltae
759 13 : if (edos_step == zero) edos_step = 0.001
760 : !if (edos_step == zero) edos_step = ten / Ha_meV
761 13 : edos_broad = dtset%tsmear
762 :
763 : ! Set default energy range for DOS
764 : ! If sigma_erange is set, get emin and emax from this variable
765 : ! MG: TODO This value should be read from SIGEPH
766 : ! Recheck metals
767 13 : if (self%assume_gap) then
768 12 : emin = huge(one); emax = -huge(one)
769 24 : do spin=1,self%ebands%nsppol
770 12 : if (dtset%sigma_erange(1) >= zero) emin = min(emin, self%gaps%vb_max(spin) + tol2 * eV_Ha - dtset%sigma_erange(1))
771 24 : if (dtset%sigma_erange(2) >= zero) emax = max(emax, self%gaps%cb_min(spin) - tol2 * eV_Ha + dtset%sigma_erange(2))
772 : end do
773 12 : ABI_CHECK(emin /= huge(one), "Cannot initialize emin")
774 12 : ABI_CHECK(emax /= -huge(one), "Cannot initialize emax")
775 : else
776 3 : emin = minval(self%eminmax_spin(1, :)); emin = emin - tol1 * abs(emin)
777 3 : emax = maxval(self%eminmax_spin(2, :)); emax = emax + tol1 * abs(emax)
778 : end if
779 :
780 : ! Compute DOS, vv_dos and vvtau_DOS (v x v tau)
781 : !
782 : ! out_valsdos: (nw, 2, nvals, nsppol) array with DOS for scalar quantities if nvals > 0
783 : ! out_tensdos: (nw, 2, 3, 3, ntens, nsppol) array with DOS weighted by tensorial terms if ntens > 0
784 : !
785 : ! Vectors and tensors are in Cartesian coordinates.
786 : ! Note how we compute the DOS only between [emin, emax] to save time and memory
787 : ! this implies that IDOS and edos%ifermi are ill-defined
788 :
789 : self%edos = self%ebands%get_edos_matrix_elements(cryst, self%bsize, &
790 : nvals, tau_vals, nvecs0, dummy_vecs, ntens, vv_tens, &
791 : edos_intmeth, edos_step, edos_broad, &
792 : out_valsdos, dummy_dosvecs, out_tensdos, comm, &
793 65 : brange=[self%bmin, self%bmax], erange=[emin, emax])
794 :
795 13 : if (my_rank == master) then
796 39 : call self%edos%print([std_out, ab_out], header="Computation of DOS, VV_DOS and VVTAU_DOS")
797 : end if
798 :
799 13 : call cwtime_report(" compute_rta_edos", cpu, wall, gflops)
800 :
801 : ! Unpack data stored in out_tensdos with shape (nw, 2, 3, 3, ntens, nsppol)
802 13 : self%nw = self%edos%nw
803 78 : ABI_MALLOC(self%tau_dos, (self%nw, self%ntemp, nsppol, self%nrta))
804 : ! TODO: Exchange dims?
805 65 : ABI_MALLOC(self%vv_dos, (self%nw, 3, 3, nsppol))
806 91 : ABI_MALLOC(self%vvtau_dos, (self%nw, 3, 3, self%ntemp, nsppol, self%nrta))
807 :
808 39 : do irta=1,self%nrta
809 65 : do spin=1,nsppol
810 202 : do itemp=1,self%ntemp+1
811 :
812 150 : itens_ = itemp + (irta - 1) * (self%ntemp + 1)
813 176 : if (itemp == 1) then
814 100238 : self%vv_dos(:,:,:,spin) = out_tensdos(:, 1, :, :, itens_, spin)
815 : else
816 233092 : self%vvtau_dos(:,:,:, itemp-1, spin, irta) = out_tensdos(:, 1, :, :, itens_, spin)
817 : end if
818 :
819 : end do
820 : end do
821 : end do
822 :
823 : ! Transfer data for tau(e)
824 39 : do irta=1,self%nrta
825 65 : do spin=1,nsppol
826 176 : do itemp=1,self%ntemp
827 124 : iscal = itemp + (irta - 1) * self%ntemp
828 25870 : self%tau_dos(:, itemp, spin, irta) = out_valsdos(:, 1, iscal, spin)
829 : end do
830 : end do
831 : end do
832 :
833 : ! Free memory
834 13 : ABI_SFREE(out_tensdos)
835 13 : ABI_SFREE(tau_vals)
836 13 : ABI_SFREE(out_valsdos)
837 13 : ABI_SFREE(dummy_dosvecs)
838 13 : ABI_SFREE(vv_tens)
839 :
840 : ! Compute Onsager coefficients. Eq 9 of [[cite:Madsen2018]]
841 : ! See also Eqs 41, page 11 of https://arxiv.org/pdf/1402.6979.pdf
842 : !
843 : ! L^\alpha(\mu, T) = \int de \sigma(e, T) (e - mu)^\alpha (-df/de)
844 : !
845 : ! with \sigma(e, T) stored in vvtau_dos
846 :
847 78 : ABI_MALLOC(self%l0, (3, 3, self%nw, self%nsppol, self%ntemp, self%nrta))
848 78 : ABI_MALLOC(self%l1, (3, 3, self%nw, self%nsppol, self%ntemp, self%nrta))
849 78 : ABI_MALLOC(self%l2, (3, 3, self%nw, self%nsppol, self%ntemp, self%nrta))
850 :
851 13 : call onsager(0, self%l0)
852 13 : call onsager(1, self%l1)
853 13 : call onsager(2, self%l2)
854 :
855 13 : call cwtime_report(" compute_rta_onsanger", cpu, wall, gflops)
856 :
857 : ! Compute transport tensors, Eqs 12-15 of [[cite:Madsen2018]] and convert to SI units.
858 334712 : ABI_CALLOC(self%sigma, (3, 3, self%nw, self%nsppol, self%ntemp, self%nrta))
859 334712 : ABI_CALLOC(self%seebeck, (3, 3, self%nw, self%nsppol, self%ntemp, self%nrta))
860 334712 : ABI_CALLOC(self%kappa, (3, 3, self%nw, self%nsppol, self%ntemp, self%nrta))
861 334712 : ABI_CALLOC(self%pi, (3, 3, self%nw, self%nsppol, self%ntemp, self%nrta))
862 334712 : ABI_CALLOC(self%zte, (3, 3, self%nw, self%nsppol, self%ntemp, self%nrta))
863 :
864 : ! Sigma = L0
865 : !TODO : missing maxocc??
866 13 : fact0 = (siemens_SI / Bohr_meter / cryst%ucvol)
867 334660 : self%sigma = fact0 * self%l0
868 :
869 : ! Used to stored L0^-1
870 39 : ABI_MALLOC(l0inv_33nw, (3, 3, self%nw))
871 :
872 39 : do irta=1,self%nrta
873 65 : do spin=1,nsppol
874 176 : do itemp=1,self%ntemp
875 :
876 124 : TKelv = self%kTmesh(itemp) / kb_HaK; if (TKelv < one) Tkelv = one
877 :
878 : ! S = -1/T L0^-1 L1 = -1/T sigma L1
879 25844 : do iw=1,self%nw
880 25720 : call inv33(self%l0(:, :, iw, spin, itemp, irta), work_33)
881 334360 : l0inv_33nw(:,:,iw) = work_33
882 3729524 : self%seebeck(:,:,iw,spin,itemp,irta) = - (volt_SI / TKelv) * matmul(work_33, self%l1(:,:,iw,spin,itemp,irta))
883 : end do
884 :
885 : ! kappa = 1/T [L2 - L1 L0^-1 L1]
886 : ! HM: Check why do we need minus sign here to get consistent results with Boltztrap!
887 : ! MG: Likely because of a different definition of kappa.
888 25844 : do iw=1,self%nw
889 334360 : work_33 = self%l1(:, :, iw, spin, itemp, irta)
890 2957800 : work_33 = self%l2(:, :, iw, spin, itemp, irta) - matmul(work_33, matmul(l0inv_33nw(:, :, iw), work_33))
891 : !self%kappa(:,:,iw,spin, itemp,spin,irta) = - (volt_SI**2 * fact0 / TKelv) * work_33
892 334484 : self%kappa(:,:,iw,spin,itemp,irta) = + (volt_SI**2 * fact0 / TKelv) * work_33
893 : end do
894 :
895 : ! Peltier pi = -L1 L0^-1
896 25844 : do iw=1,self%nw
897 334360 : work_33 = self%l1(:, :, iw, spin, itemp, irta)
898 2340644 : self%pi(:,:,iw,spin,itemp,irta) = - volt_SI * matmul(work_33, l0inv_33nw(:, :, iw))
899 : end do
900 :
901 : ! ZT: S^T sigma S k^-1 T (tensor form with k=k_electronic only):
902 25870 : do iw=1,self%nw
903 334360 : S_33 = self%seebeck(:,:,iw,spin,itemp,irta)
904 3780840 : S_33 = matmul(matmul(transpose(S_33), self%sigma(:,:,iw,spin,itemp,irta)), S_33)
905 25720 : call inv33(self%kappa(:,:,iw,spin,itemp,irta), work_33)
906 2340644 : self%zte(:,:,iw,spin,itemp,irta) = matmul(S_33, work_33) * TKelv
907 : end do
908 :
909 : end do
910 : end do
911 : end do
912 :
913 13 : ABI_FREE(l0inv_33nw)
914 :
915 : !Here are computed the transport tensors at exact mu (exact temperature and not based on a grid)
916 :
917 : ! do irta=1,self%nrta
918 : ! do spin=1,nsppol
919 : ! do itemp=1,self%ntemp
920 : !
921 : ! TKelv = self%kTmesh(itemp) / kb_HaK; if (TKelv < one) Tkelv = one
922 : !
923 : ! ! S = -1/T L11^-1 L12 = -1/T sigma^-1 L12
924 : ! call inv33(self%l11_mu(:, :, spin, itemp, irta), work_33)
925 : ! ! l0inv_33nw(:,:) = work_33
926 : ! self%seebeck_mu(:,:,spin,itemp,irta) = - (volt_SI / TKelv) * matmul(work_33, self%l12_mu(:,:,spin,itemp,irta))
927 : !
928 : !
929 : ! ! kappa = 1/T [L22 - L21 L11^-1 L12] (in RTA L12=L21)
930 : ! call inv33(self%l11_mu(:, :, spin, itemp, irta), work_33)
931 : ! work_33 = matmul(work_33,self%l12_mu(:, :, spin, itemp, irta))
932 : ! work_33 = self%l22_mu(:, :, spin, itemp, irta) - matmul(self%l12_mu(:, :, spin, itemp, irta), work_33)
933 : ! self%kappa_mu(:,:,spin,itemp,irta) = + (volt_SI**2 * fact0 / TKelv) * work_33
934 : !
935 : !
936 : ! ! Peltier pi = -L12 L11^-1 (in RTA and only in RTA, if TR sym. is not broken)
937 : !
938 : ! call inv33(self%l11_mu(:, :, spin, itemp, irta), work_33)
939 : ! self%pi_mu(:,:,spin,itemp,irta) = - volt_SI * matmul(self%l12_mu(:, :, spin, itemp, irta), work_33)
940 : !
941 : ! end do ! itemp
942 : ! end do !spin
943 : ! end do !irta
944 :
945 :
946 :
947 : ! Compute the index of the Fermi level and handle possible out of range condition.
948 13 : ifermi = bisect(self%edos%mesh, self%ebands%fermie)
949 13 : if (ifermi == 0 .or. ifermi == self%nw) then
950 0 : ABI_ERROR("Bisection could not find the index of the Fermi level in edos%mesh!")
951 : end if
952 :
953 13 : max_occ = two / (self%nspinor * self%nsppol)
954 :
955 : ! Conductivity
956 65 : ABI_MALLOC(self%conductivity, (3, 3, self%ntemp, self%nsppol, self%nrta))
957 26 : do spin=1,self%nsppol
958 88 : do itemp=1,self%ntemp
959 199 : do irta=1,self%nrta
960 558 : do jj=1,3
961 1612 : do ii=1,3
962 1488 : self%conductivity(ii,jj,itemp,spin,irta) = self%sigma(ii, jj, ifermi, spin, itemp, irta) * 0.01 !m^-1 to cm^-1
963 : end do
964 : end do
965 : end do
966 : end do ! itemp
967 : end do ! spin
968 :
969 52 : ABI_MALLOC(self%resistivity, (3, 3, self%ntemp, self%nrta))
970 39 : do irta=1,self%nrta
971 163 : do itemp=1,self%ntemp
972 2728 : work_33 = sum(self%conductivity(:,:,itemp,:,irta), dim=3)
973 1612 : call inv33(work_33, mat33); mat33 = 1e+6_dp * mat33
974 1638 : self%resistivity(:, :, itemp, irta) = mat33
975 : end do
976 : end do
977 :
978 :
979 :
980 : ! Mobility
981 65 : ABI_MALLOC(self%n, (self%nw, self%ntemp, 2))
982 91 : ABI_MALLOC(self%mobility, (3, 3, self%nw, self%ntemp, 2, self%nsppol, self%nrta))
983 :
984 26 : do spin=1,self%nsppol
985 88 : do itemp=1,self%ntemp
986 : ! Compute carrier density
987 62 : kT = self%kTmesh(itemp)
988 :
989 : ! MG TODO: I think that here we should use mu_e instead of ifermi.
990 : ! Compute carrier density of electrons (ifermi:self%nw)
991 12922 : do iw=1,self%nw ! doping
992 : self%n(iw,itemp,1) = carriers(self%edos%mesh, self%edos%dos(:,spin) * max_occ, ifermi, self%nw, &
993 34127822 : kT, self%edos%mesh(iw)) / cryst%ucvol / Bohr_meter**3
994 : end do
995 :
996 : ! Compute carrier density of holes (1:ifermi)
997 12922 : do iw=1,self%nw ! doping
998 : self%n(iw,itemp,2) = carriers(self%edos%mesh, self%edos%dos(:,spin) * max_occ, 1, ifermi, &
999 34127822 : kT, self%edos%mesh(iw)) / cryst%ucvol / Bohr_meter**3
1000 : end do
1001 :
1002 12922 : self%n(:,itemp,2) = self%n(self%nw,itemp,2) - self%n(:,itemp,2)
1003 :
1004 : ! Compute mobility
1005 199 : do irta=1,self%nrta
1006 434 : do iel=1,2
1007 51812 : do iw=1,self%nw
1008 206008 : do jj=1,3
1009 668720 : do ii=1,3
1010 : call safe_div(self%sigma(ii, jj, iw, spin, itemp, irta) * 100**2, &
1011 : e_Cb * self%n(iw, itemp, iel), &
1012 617280 : zero, self%mobility(ii, jj, iw, itemp, iel, spin, irta))
1013 : end do
1014 : end do
1015 : end do
1016 : end do
1017 : end do
1018 : end do ! itemp
1019 : end do ! spin
1020 :
1021 : ! Compute RTA mobility
1022 13 : call self%compute_rta_mobility(cryst, comm)
1023 :
1024 13 : if (my_rank == master) then
1025 : ! Print RTA results to stdout and other external txt files (for the test suite)
1026 13 : call self%print_rta_txt_files(cryst, dtset, dtfil)
1027 :
1028 : ! Creates the netcdf file used to store the results of the calculation.
1029 13 : path = strcat(dtfil%filnam_ds(4), "_RTA.nc")
1030 13 : call wrtout(units, ch10//sjoin("- Writing RTA transport results to:", path))
1031 13 : NCF_CHECK(nctk_open_create(ncid, path , xmpi_comm_self))
1032 13 : call self%rta_ncwrite(cryst, dtset, ncid)
1033 13 : NCF_CHECK(nf90_close(ncid))
1034 : end if
1035 :
1036 65 : call cwtime_report(" compute_rta", cpu, wall, gflops)
1037 :
1038 : contains
1039 :
1040 25720 : real(dp) function carriers(wmesh, dos, istart, istop, kT, mu)
1041 :
1042 : !Arguments -------------------------------------------
1043 : real(dp),intent(in) :: kT, mu
1044 : real(dp),intent(in) :: wmesh(self%nw), dos(self%nw)
1045 : integer,intent(in) :: istart, istop
1046 :
1047 : !Local variables -------------------------------------
1048 : integer :: iw
1049 25720 : real(dp) :: kernel(self%nw), integral(self%nw)
1050 :
1051 68255520 : kernel = zero
1052 34153480 : do iw=istart,istop
1053 34153480 : kernel(iw) = dos(iw) * occ_fd(wmesh(iw), kT, mu)
1054 : end do
1055 25720 : call simpson_int(self%nw, edos_step, kernel, integral)
1056 25720 : carriers = integral(self%nw)
1057 :
1058 25720 : end function carriers
1059 :
1060 : ! Compute L^\alpha(\mu, T) = \int de \sigma(e, T) (e - mu)^\alpha (-df/de)
1061 39 : subroutine onsager(order, lorder)
1062 :
1063 : !Arguments -------------------------------------------
1064 : integer,intent(in) :: order
1065 : real(dp),intent(out) :: lorder(3, 3, self%nw, self%nsppol, self%ntemp, self%nrta)
1066 :
1067 : !Local variables -------------------------------------
1068 : integer :: spin, iw, imu, irta
1069 : real(dp) :: mu, ee, kT
1070 78 : real(dp) :: kernel(self%nw,3,3,self%nsppol), integral(self%nw)
1071 :
1072 : ! Get spin degeneracy
1073 39 : max_occ = two / (self%nspinor * self%nsppol)
1074 :
1075 117 : do irta=1,self%nrta
1076 489 : do itemp=1,self%ntemp
1077 372 : kT = self%kTmesh(itemp)
1078 : ! Loop over chemical potentials mu
1079 77610 : do imu=1,self%nw
1080 77160 : mu = self%edos%mesh(imu)
1081 :
1082 : ! Build integrand for given mu
1083 204766560 : do iw=1,self%nw
1084 204689400 : ee = self%edos%mesh(iw)
1085 204766560 : if (order > 0) then
1086 1910434400 : kernel(iw,:,:,:) = - max_occ * self%vvtau_dos(iw,:,:,itemp,:,irta) * (ee - mu)** order * occ_dfde(ee, kT, mu)
1087 : else
1088 955217200 : kernel(iw,:,:,:) = - max_occ * self%vvtau_dos(iw,:,:,itemp,:,irta) * occ_dfde(ee, kT, mu)
1089 : end if
1090 : end do
1091 :
1092 : ! Integrate with simpson_int
1093 154692 : do spin=1,self%nsppol
1094 385800 : do jj=1,3
1095 1003080 : do ii=1,3
1096 694440 : call simpson_int(self%nw, edos_step, kernel(:,ii,jj, spin), integral)
1097 925920 : lorder(ii, jj, imu, spin, itemp, irta) = integral(self%nw)
1098 : end do
1099 : end do
1100 : end do
1101 :
1102 : end do ! imu
1103 : end do ! itemp
1104 : end do ! irta
1105 :
1106 39 : end subroutine onsager
1107 :
1108 : end subroutine compute_rta
1109 : !!***
1110 :
1111 : !----------------------------------------------------------------------
1112 :
1113 : !!****f* m_rta/compute_rta_mobility
1114 : !! NAME
1115 : !! compute_rta_mobility
1116 : !!
1117 : !! FUNCTION
1118 : !!
1119 : !! INPUTS
1120 : !! cryst<crystal_t>=Crystalline structure
1121 : !! comm=MPI communicator.
1122 : !!
1123 : !! SOURCE
1124 :
1125 13 : subroutine compute_rta_mobility(self, cryst, comm)
1126 :
1127 : !Arguments ------------------------------------
1128 : class(rta_t),intent(inout) :: self
1129 : type(crystal_t),intent(in) :: cryst
1130 : integer,intent(in) :: comm
1131 :
1132 : !Local variables ------------------------------
1133 : integer :: nsppol, nkibz, ib, ik_ibz, spin, ii, jj, itemp, ieh, cnt, nprocs, irta, time_opt
1134 : real(dp) :: eig_nk, mu_e, linewidth, fact, fact0, max_occ, kT, wtk, cpu, wall, gflops
1135 : real(dp) :: vr(3), vv_tens(3,3), vv_tenslw(3,3), work_33(3,3), mat33(3,3) !, tmp_tens(3,3)
1136 : !************************************************************************
1137 :
1138 13 : call cwtime(cpu, wall, gflops, "start")
1139 :
1140 13 : nprocs = xmpi_comm_size(comm)
1141 13 : nkibz = self%ebands%nkpt; nsppol = self%ebands%nsppol
1142 :
1143 13 : time_opt = 0 ! This to preserve the previous behaviour in which TR was not used.
1144 : !time_opt = -1 ! This to preserve the previous behaviour in which TR was not used.
1145 :
1146 3563 : ABI_CALLOC(self%mobility_mu, (3, 3, 2, nsppol, self%ntemp, self%nrta))
1147 3550 : ABI_CALLOC(self%conductivity_mu, (3, 3, 2, nsppol, self%ntemp, self%nrta))
1148 1827 : ABI_CALLOC(self%l11_mu, (3, 3, nsppol, self%ntemp, self%nrta))
1149 1814 : ABI_CALLOC(self%l12_mu, (3, 3, nsppol, self%ntemp, self%nrta))
1150 1814 : ABI_CALLOC(self%l22_mu, (3, 3, nsppol, self%ntemp, self%nrta))
1151 1690 : ABI_CALLOC(self%resistivity_mu, (3, 3, self%ntemp, self%nrta))
1152 3511 : self%conductivity_mu=zero
1153 1775 : self%l11_mu=zero
1154 1775 : self%l12_mu=zero
1155 1775 : self%l22_mu=zero
1156 1651 : self%resistivity_mu=zero
1157 : ! ABI_CALLOC(self%seebeck_mu, (3, 3, 2, nsppol, self%ntemp, self%nrta))
1158 : ! ABI_CALLOC(self%kappa_mu, (3, 3, 2, nsppol, self%ntemp, self%nrta))
1159 : ! ABI_CALLOC(self%peltier_mu, (3, 3, 2, nsppol, self%ntemp, self%nrta))
1160 : ! Compute (e/h) carriers per unit cell at the different temperatures.
1161 300 : ABI_CALLOC(self%n_ehst, (2, self%nsppol, self%ntemp))
1162 13 : call self%ebands%get_carriers(self%ntemp, self%kTmesh, self%transport_mu_e, self%n_ehst)
1163 :
1164 : ! Compute conductivity_mu i.e. results in which lifetimes have been computed in a consistent way
1165 : ! with the same the Fermi level. In all the other cases, indeed, we assume that tau does not depend on ef.
1166 : !
1167 : ! sigma_RTA = -S e^2 / (N_k Omega) sum_\nk (v_\nk \otimes v_\nk) \tau_\nk (df^0/de_\nk)
1168 : !
1169 : ! with S the spin degeneracy factor.
1170 : !
1171 : ! TODO: Implement other tensors. Compare these results with the ones obtained with spectral sigma
1172 : ! In principle, they should be the same, in practice the integration of sigma requires enough resolution
1173 : ! around the band edge.
1174 : !print *, "in RTA max velocities", maxval(abs(self%vbks))
1175 13 : cnt = 0
1176 26 : do spin=1,nsppol
1177 3582 : do ik_ibz=1,nkibz
1178 : !cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! MPI parallelism.
1179 3556 : wtk = self%ebands%wtk(ik_ibz)
1180 :
1181 7230 : do ib=self%bmin,self%bmax
1182 3661 : eig_nk = self%ebands%eig(ib, ik_ibz, spin)
1183 :
1184 : ! Store outer product in vv_tens
1185 14644 : vr(:) = self%vbks(:, ib, ik_ibz, spin)
1186 : ! Don't remove this if: it makes the loop a bit faster and, most importantly,
1187 : ! it prevents intel from miscompiling the code.
1188 14472 : if (all(abs(vr) == zero)) cycle
1189 :
1190 232 : do ii=1,3
1191 754 : do jj=1,3
1192 696 : vv_tens(ii, jj) = vr(ii) * vr(jj)
1193 : end do
1194 : end do
1195 :
1196 : ! Symmetrize tensor.
1197 : !print *, "intens", vv_tens
1198 754 : vv_tens = cryst%symmetrize_cart_tens33(vv_tens, time_opt)
1199 : !print *, "out_tens", vv_tens
1200 :
1201 : ! Multiply by the lifetime (SERTA or MRTA)
1202 3730 : do irta=1,self%nrta
1203 4273 : do itemp=1,self%ntemp
1204 496 : kT = self%kTmesh(itemp)
1205 496 : mu_e = self%transport_mu_e(itemp)
1206 496 : ieh = 2; if (eig_nk >= mu_e) ieh = 1
1207 496 : linewidth = self%linewidths(itemp, ib, ik_ibz, spin, irta)
1208 : !print *, linewidth, wtk, occ_dfde(eig_nk, kT, mu_e), "tens", vv_tens
1209 6448 : call safe_div( - wtk * vv_tens * occ_dfde(eig_nk, kT, mu_e), two * linewidth, zero, vv_tenslw)
1210 : self%conductivity_mu(:, :, ieh, spin, itemp, irta) = self%conductivity_mu(:, :, ieh, spin, itemp, irta) &
1211 6448 : + vv_tenslw(:, :)
1212 15376 : self%l11_mu(:, :, spin, itemp, irta) = sum(self%conductivity_mu(:,:,:, spin, itemp, irta), dim=3) !sum over ieh (holes and e-)
1213 : !Here I implement the other onsager coefficients at the correct mu
1214 6448 : call safe_div( - wtk * vv_tens * occ_dfde(eig_nk, kT, mu_e)*(eig_nk-mu_e), two * linewidth, zero, vv_tenslw)
1215 6448 : self%l12_mu(:, :, spin, itemp, irta)=self%l12_mu(:, :, spin, itemp, irta)+vv_tenslw(:, :)
1216 6448 : call safe_div( - wtk * vv_tens * occ_dfde(eig_nk, kT, mu_e)*(eig_nk-mu_e)**2, two * linewidth, zero, vv_tenslw)
1217 6564 : self%l22_mu(:, :, spin, itemp, irta)=self%l22_mu(:, :, spin, itemp, irta)+vv_tenslw(:, :)
1218 : end do
1219 : end do
1220 : end do
1221 :
1222 : end do ! ik_ibz
1223 : end do ! spin
1224 :
1225 : !call xmpi_sum(self%conductivity_mu, comm, ierr)
1226 :
1227 : ! Get units conversion factor including spin degeneracy.
1228 13 : max_occ = two / (self%nspinor * self%nsppol)
1229 13 : fact0 = max_occ * (siemens_SI / Bohr_meter / cryst%ucvol) / 100
1230 3511 : self%conductivity_mu = fact0 * self%conductivity_mu ! siemens cm^-1
1231 1775 : self%l11_mu = max_occ * self%l11_mu
1232 1775 : self%l12_mu = max_occ * self%l12_mu
1233 1775 : self%l22_mu = max_occ * self%l22_mu
1234 :
1235 : !Same for resistivity_mu (at correct mu)
1236 :
1237 : ! ABI_MALLOC(self%resistivity_mu, (3, 3, self%ntemp, self%nrta))
1238 39 : do irta=1,self%nrta
1239 163 : do itemp=1,self%ntemp
1240 :
1241 2728 : work_33 = sum(self%l11_mu(:,:,:,itemp,irta), dim=3) !sum over spins
1242 1612 : work_33 = work_33*fact0
1243 1612 : call inv33(work_33, mat33); mat33 = 1e+6_dp * mat33
1244 1638 : self%resistivity_mu(:, :, itemp, irta) = mat33
1245 : end do
1246 : end do
1247 :
1248 :
1249 : !TODO: implement for mobility as well (verify)
1250 : ! Scale by the carrier concentration
1251 : fact = 100**3 / e_Cb
1252 39 : do irta=1,self%nrta
1253 65 : do spin=1,nsppol
1254 176 : do itemp=1,self%ntemp
1255 398 : do ieh=1,2 ! e/h
1256 : call safe_div(fact * self%conductivity_mu(:,:,ieh,spin,itemp, irta), &
1257 : self%n_ehst(ieh, spin, itemp) / cryst%ucvol / Bohr_meter**3, zero, &
1258 3348 : self%mobility_mu(:,:,ieh,spin,itemp,irta))
1259 : end do
1260 : end do
1261 : end do
1262 : end do
1263 :
1264 13 : call cwtime_report(" compute_rta_mobility", cpu, wall, gflops)
1265 :
1266 26 : end subroutine compute_rta_mobility
1267 : !!***
1268 :
1269 : !----------------------------------------------------------------------
1270 :
1271 : !!****f* m_rta/rta_ncwrite
1272 : !! NAME
1273 : !! rta_ncwrite
1274 : !!
1275 : !! FUNCTION
1276 : !!
1277 : !! INPUTS
1278 : !! cryst<crystal_t>=Crystalline structure
1279 : !! dtset<dataset_type>=All input variables for this dataset.
1280 : !! ncid=Netcdf file handle.
1281 : !!
1282 : !! SOURCE
1283 :
1284 13 : subroutine rta_ncwrite(self, cryst, dtset, ncid)
1285 :
1286 : !Arguments --------------------------------------
1287 : class(rta_t),intent(in) :: self
1288 : type(crystal_t),intent(in) :: cryst
1289 : type(dataset_type),intent(in) :: dtset
1290 : integer,intent(in) :: ncid
1291 :
1292 : !Local variables --------------------------------
1293 : integer :: ncerr, ii
1294 : real(dp) :: cpu, wall, gflops
1295 26 : real(dp) :: work(dtset%nsppol)
1296 : !************************************************************************
1297 :
1298 13 : call cwtime(cpu, wall, gflops, "start")
1299 :
1300 : ! Write to netcdf file
1301 13 : NCF_CHECK(cryst%ncwrite(ncid))
1302 13 : NCF_CHECK(self%ebands%ncwrite(ncid))
1303 13 : NCF_CHECK(self%edos%ncwrite(ncid))
1304 :
1305 : !nctk_copy from sigeph?
1306 : ! nctkarr_t("eph_ngqpt_fine", "int", "three"), &
1307 :
1308 : ncerr = nctk_def_dims(ncid, [ &
1309 : nctkdim_t("ntemp", self%ntemp), nctkdim_t("nrta", self%nrta), nctkdim_t("nsppol", self%nsppol), &
1310 : nctkdim_t("nkcalc", self%nkcalc), nctkdim_t("nkibz", self%ebands%nkpt) &
1311 78 : ], defmode=.True.)
1312 13 : NCF_CHECK(ncerr)
1313 :
1314 : ncerr = nctk_def_arrays(ncid, [ &
1315 : nctkarr_t('transport_ngkpt', "int", "three"), &
1316 : nctkarr_t('sigma_erange', "dp", "two"), &
1317 : nctkarr_t("kcalc2ibz", "int", "nkcalc, six"), &
1318 : nctkarr_t("kcalc2ebands", "int", "six, nkcalc"), &
1319 : nctkarr_t("kibz", "dp", "three, nkibz"), &
1320 : nctkarr_t('kTmesh', "dp", "ntemp"), &
1321 : nctkarr_t('transport_mu_e', "dp", "ntemp"), &
1322 : nctkarr_t('n_ehst', "dp", "two, nsppol, ntemp"), &
1323 : nctkarr_t('eph_mu_e', "dp", "ntemp"), &
1324 : nctkarr_t('vb_max', "dp", "nsppol"), &
1325 : nctkarr_t('cb_min', "dp", "nsppol"), &
1326 : nctkarr_t('vv_dos', "dp", "edos_nw, three, three, nsppol"), &
1327 : nctkarr_t('vvtau_dos', "dp", "edos_nw, three, three, ntemp, nsppol, nrta"), &
1328 : nctkarr_t('tau_dos', "dp", "edos_nw, ntemp, nsppol, nrta"), &
1329 : nctkarr_t('L0', "dp", "three, three, edos_nw, nsppol, ntemp, nrta"), &
1330 : nctkarr_t('L1', "dp", "three, three, edos_nw, nsppol, ntemp, nrta"), &
1331 : nctkarr_t('L2', "dp", "three, three, edos_nw, nsppol, ntemp, nrta"), &
1332 : nctkarr_t('sigma', "dp", "three, three, edos_nw, nsppol, ntemp, nrta"), &
1333 : nctkarr_t('kappa', "dp", "three, three, edos_nw, nsppol, ntemp, nrta"), &
1334 : nctkarr_t('zte', "dp", "three, three, edos_nw, nsppol, ntemp, nrta"), &
1335 : nctkarr_t('seebeck', "dp", "three, three, edos_nw, nsppol, ntemp, nrta"), &
1336 : nctkarr_t('pi', "dp", "three, three, edos_nw, nsppol, ntemp, nrta"), &
1337 : nctkarr_t('mobility', "dp", "three, three, edos_nw, ntemp, two, nsppol, nrta"), &
1338 : nctkarr_t('conductivity', "dp", "three, three, ntemp, nsppol, nrta"), &
1339 : nctkarr_t('resistivity', "dp", "three, three, ntemp, nrta"), &
1340 : nctkarr_t('N', "dp", "edos_nw, ntemp, two"), &
1341 : !nctkarr_t('conductivity_mu',"dp", "three, three, two, nsppol, ntemp, nrta")], &
1342 : nctkarr_t('mobility_mu', "dp", "three, three, two, nsppol, ntemp, nrta")], &
1343 364 : defmode=.True.)
1344 13 : NCF_CHECK(ncerr)
1345 :
1346 26 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "assume_gap"])
1347 13 : NCF_CHECK(ncerr)
1348 : ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: &
1349 65 : "eph_extrael", "eph_fermie", "transport_extrael", "transport_fermie"])
1350 13 : NCF_CHECK(ncerr)
1351 :
1352 : ! Write data.
1353 13 : ii = 0; if (self%assume_gap) ii = 1
1354 39 : ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: "assume_gap"], [ii], datamode=.True.)
1355 :
1356 : ncerr = nctk_write_dpscalars(ncid, [character(len=nctk_slen) :: &
1357 : "eph_extrael", "eph_fermie", "transport_extrael", "transport_fermie"], &
1358 117 : [self%eph_extrael, self%eph_fermie, self%transport_extrael, self%transport_fermie])
1359 13 : NCF_CHECK(ncerr)
1360 :
1361 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "transport_ngkpt"), dtset%transport_ngkpt))
1362 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "sigma_erange"), dtset%sigma_erange))
1363 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kcalc2ibz"), self%kcalc2ibz))
1364 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kcalc2ebands"), self%kcalc2ebands))
1365 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kibz"), self%ebands%kptns))
1366 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kTmesh"), self%kTmesh))
1367 13 : if (self%assume_gap) then
1368 12 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vb_max"), self%gaps%vb_max))
1369 12 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "cb_min"), self%gaps%cb_min))
1370 : else
1371 : ! Set vbm and cbm to fermie if metal.
1372 2 : work(:) = self%ebands%fermie
1373 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vb_max"), work))
1374 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "cb_min"), work))
1375 : end if
1376 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eph_mu_e"), self%eph_mu_e))
1377 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "transport_mu_e"), self%transport_mu_e))
1378 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "n_ehst"), self%n_ehst))
1379 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vv_dos"), self%vv_dos))
1380 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vvtau_dos"), self%vvtau_dos))
1381 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "tau_dos"), self%tau_dos))
1382 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "L0"), self%l0))
1383 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "L1"), self%l1))
1384 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "L2"), self%l2))
1385 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "sigma"), self%sigma))
1386 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kappa"), self%kappa))
1387 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "zte"), self%zte))
1388 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "seebeck"), self%seebeck))
1389 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "pi"), self%pi))
1390 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "N"), self%n))
1391 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "mobility"), self%mobility))
1392 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "conductivity"), self%conductivity))
1393 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "resistivity"), self%resistivity))
1394 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "mobility_mu"), self%mobility_mu))
1395 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "conductivity_mu"), self%conductivity_mu))
1396 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "resistivity_mu"), self%resistivity_mu))
1397 :
1398 13 : call cwtime_report(" rta_ncwrite", cpu, wall, gflops)
1399 :
1400 13 : end subroutine rta_ncwrite
1401 : !!***
1402 :
1403 : !----------------------------------------------------------------------
1404 :
1405 : !!****f* m_rta/print_rta_txt_files
1406 : !! NAME
1407 : !! print_rta_txt_files
1408 : !!
1409 : !! FUNCTION
1410 : !!
1411 : !! INPUTS
1412 : !! cryst<crystal_t>=Crystalline structure
1413 : !! dtset<dataset_type>=All input variables for this dataset.
1414 : !! dtfil<datafiles_type>=variables related to files.
1415 : !!
1416 : !! SOURCE
1417 :
1418 13 : subroutine print_rta_txt_files(self, cryst, dtset, dtfil)
1419 :
1420 : !Arguments --------------------------------------
1421 : class(rta_t),intent(in) :: self
1422 : type(crystal_t),intent(in) :: cryst
1423 : type(dataset_type),intent(in) :: dtset
1424 : type(datafiles_type),intent(in) :: dtfil
1425 :
1426 : !Local variables --------------------------------
1427 : integer :: itemp, spin, irta, ii, nsp
1428 : real(dp) :: TKelv
1429 : character(len=500) :: msg, pre, rta_type
1430 : integer :: units(2)
1431 : character(len=2) :: components(3)
1432 : real(dp) :: mat33(3,3), work33(3,3)
1433 : !************************************************************************
1434 :
1435 39 : units = [std_out, ab_out]
1436 13 : call wrtout(units, ch10//' Transport (RTA) calculation results, chemical potential adjusted with T:', newlines=1)
1437 52 : components = ["xx", "yy", "zz"]
1438 :
1439 39 : do irta=1,self%nrta
1440 26 : if (irta == 1) rta_type = "SERTA"
1441 26 : if (irta == 2) rta_type = "MRTA"
1442 :
1443 26 : if (self%assume_gap) then
1444 : ! SemiConductor
1445 96 : do ii=1,3
1446 72 : call wrtout(units, sjoin(" Cartesian component of", rta_type, "mobility tensor:", components(ii)))
1447 72 : write(msg, "(a16,a32,a32)") 'Temperature [K]', 'e/h density [cm^-3]', 'e/h mobility [cm^2/Vs]'
1448 72 : call wrtout(units, msg)
1449 :
1450 144 : do spin=1,self%nsppol
1451 72 : if (self%nsppol == 2) call wrtout(units, sjoin(" For spin:", stoa(spin)), newlines=1)
1452 :
1453 504 : do itemp=1,self%ntemp
1454 : write(msg,"(f16.2,2e16.2,2f16.2)") &
1455 360 : self%kTmesh(itemp) / kb_HaK, &
1456 360 : self%n_ehst(1, spin, itemp) / cryst%ucvol / Bohr_cm**3, &
1457 360 : self%n_ehst(2, spin, itemp) / cryst%ucvol / Bohr_cm**3, &
1458 360 : self%mobility_mu(ii, ii, 1, spin, itemp, irta), &
1459 720 : self%mobility_mu(ii, ii, 2, spin, itemp, irta)
1460 432 : call wrtout(units, msg)
1461 : end do ! itemp
1462 : end do ! spin
1463 96 : call wrtout(units, ch10)
1464 : end do ! ii
1465 :
1466 : else
1467 : ! Metals. Print conductivity (spin resolved) and resistivity (no spin resolved)
1468 : ! do ii=1,2
1469 : ! if (ii == 1) msg = sjoin(" Conductivity [Siemens cm^-1] using ", rta_type, "approximation")
1470 : ! if (ii == 2) msg = sjoin(" Resistivity [micro-Ohm cm] using ", rta_type, "approximation")
1471 : ! call wrtout(units, msg)
1472 : !
1473 : ! nsp = self%nsppol; if (ii == 2) nsp = 1
1474 : ! do spin=1,nsp
1475 : ! if (nsp == 2) call wrtout(units, sjoin(" For spin:", stoa(spin)), newlines=1)
1476 : ! write(msg, "(4a16)") 'Temperature (K)', 'xx', 'yy', 'zz'
1477 : ! call wrtout(units, msg)
1478 : ! do itemp=1,self%ntemp
1479 : ! if (ii == 1) then
1480 : ! mat33 = self%conductivity(:,:,itemp,spin,irta)
1481 : ! else
1482 : ! mat33 = self%resistivity(:,:,itemp,irta)
1483 : ! end if
1484 : ! write(msg,"(f16.2,3e16.6)") self%kTmesh(itemp) / kb_HaK, mat33(1,1), mat33(2,2), mat33(3,3)
1485 : ! call wrtout(units, msg)
1486 : ! end do !itemp
1487 : ! end do !spin
1488 : ! call wrtout(units, ch10)
1489 : ! end do
1490 :
1491 : !For coefficients calculated at exact mu (i.e for mu calculated at exact temp.):
1492 : !TODO: do the same for SM and thus mobility
1493 :
1494 6 : do ii=1,2
1495 4 : if (ii == 1) msg = sjoin(" Conductivity [Siemens cm^-1] using ", rta_type, "approximation")
1496 4 : if (ii == 2) msg = sjoin(" Resistivity [micro-Ohm cm] using ", rta_type, "approximation")
1497 4 : call wrtout(units, msg)
1498 :
1499 12 : do itemp=1, self%ntemp
1500 8 : nsp = self%nsppol; if (ii == 2) nsp = 1
1501 20 : do spin=1,nsp
1502 8 : if (nsp == 2) call wrtout(units, sjoin(" For spin:", stoa(spin)), newlines=1)
1503 8 : write(msg, "(4a16)") 'Temperature (K)', 'xx', 'yy', 'zz'
1504 8 : call wrtout(units, msg)
1505 :
1506 8 : if (ii == 1) then
1507 124 : mat33 = sum(self%conductivity_mu(:,:,:,spin, itemp, irta), dim=3)
1508 : else
1509 52 : mat33 = self%resistivity_mu(:,:,itemp,irta)
1510 : end if
1511 8 : write(msg,"(f16.2,3e16.6)") self%kTmesh(itemp) / kb_HaK, mat33(1,1), mat33(2,2), mat33(3,3)
1512 16 : call wrtout(units, msg)
1513 : end do !spin
1514 : end do !itemp
1515 6 : call wrtout(units, ch10)
1516 : end do
1517 : end if
1518 :
1519 26 : msg = sjoin(" Seebeck [Volts / Kelvin] using ", rta_type, "approximation")
1520 26 : call wrtout(units, msg)
1521 150 : do itemp=1, self%ntemp
1522 : !TODO : what is the following line? Resets 0 Kelvin to 1 Kelvin minimum T??? a bit arbitrary. Should handle T=0 more cleanly
1523 124 : TKelv = self%kTmesh(itemp) / kb_HaK; if (TKelv < one) Tkelv = one
1524 124 : nsp = self%nsppol;
1525 274 : do spin=1,nsp
1526 124 : if (nsp == 2) call wrtout(units, sjoin(" For spin:", stoa(spin)), newlines=1)
1527 124 : write(msg, "(4a16)") 'Temperature (K)', 'xx', 'yy', 'zz'
1528 124 : call wrtout(units, msg)
1529 124 : call inv33(self%l11_mu(:, :, spin, itemp, irta), work33)
1530 9796 : mat33 = - (volt_SI / TKelv) * matmul(work33, self%l12_mu(:,:,spin,itemp,irta))
1531 124 : write(msg,"(f16.2,3e16.6)") TKelv, mat33(1,1), mat33(2,2), mat33(3,3)
1532 248 : call wrtout(units, msg)
1533 : end do !spin
1534 : end do !itemp
1535 26 : call wrtout(units, ch10)
1536 :
1537 26 : msg = sjoin(" Kappa [W/m*K] using ", rta_type, "approximation")
1538 26 : call wrtout(units, msg)
1539 150 : do itemp=1, self%ntemp
1540 124 : TKelv = self%kTmesh(itemp) / kb_HaK; if (TKelv < one) Tkelv = one
1541 124 : nsp = self%nsppol;
1542 274 : do spin=1,nsp
1543 124 : if (nsp == 2) call wrtout(units, sjoin(" For spin:", stoa(spin)), newlines=1)
1544 124 : write(msg, "(4a16)") 'Temperature (K)', 'xx', 'yy', 'zz'
1545 124 : call wrtout(units, msg)
1546 124 : call inv33(self%l11_mu(:, :, spin, itemp, irta), work33)
1547 26412 : mat33 = + (volt_SI**2 * (siemens_SI / Bohr_meter / cryst%ucvol) / TKelv) * (self%l22_mu(:,:,spin,itemp,irta) - matmul(self%l12_mu(:,:,spin,itemp,irta),matmul(work33,self%l12_mu(:,:,spin,itemp,irta))))
1548 124 : write(msg,"(f16.2,3e16.6)") TKelv, mat33(1,1), mat33(2,2), mat33(3,3)
1549 248 : call wrtout(units, msg)
1550 : end do !spin
1551 : end do !itemp
1552 26 : call wrtout(units, ch10)
1553 :
1554 26 : msg = sjoin(" Peltier [Volts] using ", rta_type, "approximation")
1555 26 : call wrtout(units, msg)
1556 150 : do itemp=1, self%ntemp
1557 124 : TKelv = self%kTmesh(itemp) / kb_HaK; if (TKelv < one) Tkelv = one
1558 124 : nsp = self%nsppol;
1559 274 : do spin=1,nsp
1560 124 : if (nsp == 2) call wrtout(units, sjoin(" For spin:", stoa(spin)), newlines=1)
1561 124 : write(msg, "(4a16)") 'Temperature (K)', 'xx', 'yy', 'zz'
1562 124 : call wrtout(units, msg)
1563 124 : call inv33(self%l11_mu(:, :, spin, itemp, irta), work33)
1564 9920 : mat33 = - volt_SI * matmul(self%l12_mu(:,:,spin,itemp,irta), work33)
1565 124 : write(msg,"(f16.2,3e16.6)") TKelv, mat33(1,1), mat33(2,2), mat33(3,3)
1566 248 : call wrtout(units, msg)
1567 : end do !spin
1568 : end do !itemp
1569 39 : call wrtout(units, ch10)
1570 :
1571 :
1572 :
1573 : end do ! irta
1574 :
1575 39 : do irta=1,self%nrta
1576 13 : select case (irta)
1577 : case (1)
1578 13 : pre = "_SERTA"
1579 : case (2)
1580 13 : pre = "_MRTA"
1581 : case default
1582 26 : ABI_ERROR(sjoin("Don't know how to handle irta:", itoa(irta)))
1583 : end select
1584 26 : call self%write_tensor(dtset, irta, "sigma", self%sigma(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_SIGMA"))
1585 26 : call self%write_tensor(dtset, irta, "seebeck", self%seebeck(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_SBK"))
1586 26 : call self%write_tensor(dtset, irta, "kappa", self%kappa(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_KAPPA"))
1587 26 : call self%write_tensor(dtset, irta, "zte", self%zte(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_ZTE"))
1588 39 : call self%write_tensor(dtset, irta, "pi", self%pi(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_PI"))
1589 : end do
1590 :
1591 13 : end subroutine print_rta_txt_files
1592 : !!***
1593 :
1594 : !----------------------------------------------------------------------
1595 :
1596 : !!****f* m_rta/write_tensor
1597 : !! NAME
1598 : !!
1599 : !! FUNCTION
1600 : !!
1601 : !! INPUTS
1602 : !!
1603 : !! SOURCE
1604 :
1605 130 : subroutine write_tensor(self, dtset, irta, header, values, path)
1606 :
1607 : !Arguments --------------------------------------
1608 : class(rta_t),intent(in) :: self
1609 : type(dataset_type),intent(in) :: dtset
1610 : integer,intent(in) :: irta
1611 : character(len=*),intent(in) :: header
1612 : real(dp),intent(in) :: values(:,:,:,:,:)
1613 : character(len=*),intent(in) :: path
1614 :
1615 : !Local variables --------------------------------
1616 : integer :: itemp, iw, ount
1617 : character(len=500) :: msg, rta_type
1618 130 : real(dp),allocatable :: tmp_values(:,:,:,:,:)
1619 : !************************************************************************
1620 :
1621 130 : if (open_file(trim(path), msg, newunit=ount, form="formatted", action="write", status='unknown') /= 0) then
1622 0 : ABI_ERROR(msg)
1623 : end if
1624 :
1625 130 : if (irta == 1) rta_type = "RTA type: Self-energy relaxation time approximation (SERTA)"
1626 130 : if (irta == 2) rta_type = "RTA type: Momentum relaxation time approximation (MRTA)"
1627 :
1628 : ! write header
1629 130 : write(ount, "(2a)")"# ", trim(header)
1630 130 : write(ount, "(2a)")"# ", trim(rta_type)
1631 : ! TODO: Units ?
1632 130 : write(ount, "(a, 3(i0, 1x))")"#", dtset%transport_ngkpt
1633 130 : write(ount, "(a)")"#"
1634 :
1635 : ! This to improve portability of the unit tests.
1636 130 : call alloc_copy(values, tmp_values)
1637 1673170 : where (abs(values) > tol30)
1638 : tmp_values = values
1639 : else where
1640 : tmp_values = zero
1641 : end where
1642 :
1643 : ! (nw, 3, 3, nsppol, ntemp)
1644 130 : if (self%nsppol == 1) then
1645 750 : do itemp=1, self%ntemp
1646 620 : write(ount, "(/, a, 1x, f16.2)")"# T = ", self%kTmesh(itemp) / kb_HaK
1647 620 : write(ount, "(a)")"# Energy [Ha], (xx, yx, zx, xy, yy, zy, xz, yz, zz) Cartesian components of tensor."
1648 129350 : do iw=1,self%nw
1649 129220 : write(ount, "(10(es16.6))")self%edos%mesh(iw), tmp_values(:, :, iw, 1, itemp)
1650 : end do
1651 : end do
1652 130 : write(ount, "(a)")""
1653 : else
1654 0 : do itemp=1, self%ntemp
1655 0 : write(ount, "(/, a, 1x, f16.2)")"# T = ", self%kTmesh(itemp) / kb_HaK
1656 : write(ount, "(a)") &
1657 0 : "# Energy [Ha], (xx, yx, zx, xy, yy, zy, xz, yz, zz) Cartesian components of tensor for spin up followed by spin down."
1658 0 : do iw=1,self%nw
1659 0 : write(ount, "(19(es16.6))")self%edos%mesh(iw), tmp_values(:, :, iw, 1, itemp), tmp_values(:, :, iw, 2, itemp)
1660 : end do
1661 : end do
1662 0 : write(ount, "(a)")""
1663 : end if
1664 :
1665 130 : close(ount)
1666 :
1667 130 : ABI_FREE(tmp_values)
1668 :
1669 130 : end subroutine write_tensor
1670 : !!***
1671 :
1672 : !----------------------------------------------------------------------
1673 :
1674 : !!****f* m_rta/rta_free
1675 : !! NAME
1676 : !! rta_free
1677 : !!
1678 : !! FUNCTION
1679 : !! Free dynamic memory.
1680 : !!
1681 : !! INPUTS
1682 : !!
1683 : !! SOURCE
1684 :
1685 13 : subroutine rta_free(self)
1686 :
1687 : !Arguments --------------------------------------
1688 : class(rta_t),intent(inout) :: self
1689 :
1690 13 : ABI_SFREE(self%n)
1691 13 : ABI_SFREE(self%vv_dos)
1692 13 : ABI_SFREE(self%vvtau_dos)
1693 13 : ABI_SFREE(self%tau_dos)
1694 13 : ABI_SFREE(self%bstart_ks)
1695 13 : ABI_SFREE(self%bstop_ks)
1696 13 : ABI_SFREE(self%nbcalc_ks)
1697 13 : ABI_SFREE(self%kcalc2ibz)
1698 13 : ABI_SFREE(self%kcalc2ebands)
1699 13 : ABI_SFREE(self%kTmesh)
1700 13 : ABI_SFREE(self%eminmax_spin)
1701 13 : ABI_SFREE(self%eph_mu_e)
1702 13 : ABI_SFREE(self%transport_mu_e)
1703 13 : ABI_SFREE(self%vbks)
1704 13 : ABI_SFREE(self%linewidths)
1705 13 : ABI_SFREE(self%l0)
1706 13 : ABI_SFREE(self%l1)
1707 13 : ABI_SFREE(self%l2)
1708 13 : ABI_SFREE(self%sigma)
1709 13 : ABI_SFREE(self%mobility)
1710 13 : ABI_SFREE(self%conductivity)
1711 13 : ABI_SFREE(self%resistivity)
1712 13 : ABI_SFREE(self%seebeck)
1713 13 : ABI_SFREE(self%kappa)
1714 13 : ABI_SFREE(self%zte)
1715 13 : ABI_SFREE(self%pi)
1716 13 : ABI_SFREE(self%mobility_mu)
1717 13 : ABI_SFREE(self%conductivity_mu)
1718 13 : ABI_SFREE(self%resistivity_mu)
1719 13 : ABI_SFREE(self%l11_mu)
1720 13 : ABI_SFREE(self%l12_mu)
1721 13 : ABI_SFREE(self%l22_mu)
1722 13 : ABI_SFREE(self%n_ehst)
1723 :
1724 13 : call self%ebands%free()
1725 13 : call self%gaps%free()
1726 13 : call self%edos%free()
1727 :
1728 13 : end subroutine rta_free
1729 : !!***
1730 :
1731 : !----------------------------------------------------------------------
1732 :
1733 : !!****f* m_rta/ibte_driver
1734 : !! NAME
1735 : !! ibte_driver
1736 : !!
1737 : !! FUNCTION
1738 : !! Driver to compute transport properties within the IBTE.
1739 : !!
1740 : !! INPUTS
1741 : !! dtfil<datafiles_type>=variables related to files.
1742 : !! ngfftc(18)=Coarse FFT mesh
1743 : !! dtset<dataset_type>=All input variables for this dataset.
1744 : !! ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...)
1745 : !! cryst<crystal_t>=Crystalline structure
1746 : !! pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
1747 : !! psps<pseudopotential_type>=Variables related to pseudopotentials.
1748 : !! comm=MPI communicator.
1749 : !!
1750 : !! SOURCE
1751 :
1752 1 : subroutine ibte_driver(dtfil, ngfftc, dtset, ebands, cryst, pawtab, psps, comm)
1753 :
1754 : !Arguments ------------------------------------
1755 : !scalars
1756 : integer, intent(in) :: comm
1757 : type(datafiles_type),intent(in) :: dtfil
1758 : type(dataset_type),intent(in) :: dtset
1759 : type(crystal_t),intent(in) :: cryst
1760 : type(ebands_t),intent(in) :: ebands
1761 : type(pseudopotential_type),intent(in) :: psps
1762 : !arrays
1763 : integer,intent(in) :: ngfftc(18)
1764 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
1765 :
1766 : !Local variables ------------------------------
1767 : integer,parameter :: master = 0
1768 : integer :: iet, max_et
1769 : integer :: spin, ikcalc, nkcalc, nbsum, nbcalc, itemp, iter, ierr, bsize
1770 : integer :: nkibz, nsppol, band_k, ik_ibz, bmin, bmax, band_sum, ntemp, ii, jj, iq_sum, nsp
1771 : integer :: ikq_ibz, isym_kq, trev_kq, cnt, tag, nprocs, receiver, my_rank, isym, itime, isym_lgk
1772 : integer :: ncid, grp_ncid, ncerr
1773 : real(dp) :: kT, mu_e, e_nk, dfde_nk, tau_nk, lw_nk, max_adiff, cpu, wall, gflops, abs_tol, rtmp
1774 : logical :: send_data
1775 : character(len=500) :: msg
1776 : character(len=fnlen) :: path
1777 1 : type(rta_t) :: ibte
1778 : !arrays
1779 : integer :: units(2), dims(4)
1780 1 : logical,allocatable :: converged(:)
1781 : real(dp) :: vec3(3), sym_vec(3), mat33(3,3), f_kq(3), work33(3,3)
1782 : real(dp) :: inv_sig_p(3,3)
1783 3 : real(dp) :: fsum_eh(3,2,ebands%nsppol), max_adiff_spin(ebands%nsppol)
1784 2 : real(dp) :: onsager(3,3,3,ebands%nsppol)
1785 1 : real(dp),pointer :: sig_p(:,:,:,:), mob_p(:,:,:,:), sbk_p(:,:,:), kappa_p(:,:,:), pi_p(:,:,:)
1786 1 : real(dp),target,allocatable :: ibte_sigma(:,:,:,:,:), ibte_mob(:,:,:,:,:), ibte_rho(:,:,:)
1787 1 : real(dp),target,allocatable :: ibte_seebeck(:,:,:,:), ibte_kappa(:,:,:,:), ibte_pi(:,:,:,:)
1788 1 : real(dp),allocatable :: grp_srate(:,:,:,:), fkn_in(:,:,:,:), fkn_out(:,:,:,:), fkn_efield(:,:,:,:), fkn_serta(:,:,:,:), taukn_serta(:,:,:,:)
1789 : real(dp) :: fact_sbk
1790 : character(len=2) :: components(3)
1791 1 : real(dp), allocatable :: sig_gen(:,:,:,:), mob_gen(:,:,:,:), sig_l21(:,:,:,:), sig_l22(:,:,:,:), mob_21(:,:,:,:), mob_22(:,:,:,:)
1792 : type :: scatk_t
1793 :
1794 : integer :: rank = xmpi_undefined_rank
1795 :
1796 : integer :: nq_ibzk_eff
1797 : ! Number of effective q-points in the IBZ(k)
1798 :
1799 : integer :: lgk_nsym
1800 : ! Number of symmetry operations in the little group of k.
1801 :
1802 : integer,allocatable :: lgk_sym2glob(:,:)
1803 : ! lgk_sym2glob(2, lgk_nsym)
1804 : ! Mapping isym_lg --> [isym, itime]
1805 : ! where isym is the index of the operation in the global array **crystal%symrec**
1806 : ! and itim is 2 if time-reversal T must be included else 1. Depends on ikcalc
1807 :
1808 : integer,allocatable :: kq_symtab(:,:)
1809 : ! kq_symtab(6, nq_ibzk_eff)
1810 :
1811 : real(dp),allocatable :: vals(:,:,:,:)
1812 : ! (nq_ibzk_eff, nbsum, nbcalc, ntemp)
1813 : end type scatk_t
1814 :
1815 1 : type(scatk_t),target,allocatable :: sr(:,:)
1816 : type(scatk_t),pointer :: sr_p
1817 :
1818 : ! *************************************************************************
1819 :
1820 1 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1821 3 : units = [std_out, ab_out]
1822 :
1823 1 : call wrtout(units, ch10//' Entering IBTE driver.')
1824 : call wrtout(units, sjoin("- Reading SERTA lifetimes and e-ph scattering operator from:", &
1825 1 : dtfil%filsigephin), newlines=1, do_flush=.True.)
1826 :
1827 : ! Initialize IBTE object
1828 1 : ibte = rta_new(dtset, dtfil, ngfftc, cryst, ebands, pawtab, psps, comm)
1829 :
1830 : ! Compute RTA transport quantities
1831 1 : call ibte%compute_rta(cryst, dtset, dtfil, comm)
1832 :
1833 1 : nkcalc = ibte%nkcalc
1834 1 : nkibz = ibte%ebands%nkpt; nsppol = ibte%nsppol; ntemp = ibte%ntemp
1835 1 : bmin = ibte%bmin; bmax = ibte%bmax
1836 : !call wrtout(std_out, sjoin(" nkcalc", itoa(nkcalc), "bmin:", itoa(bmin), "bmax:", itoa(bmax)))
1837 :
1838 : ! Get units conversion factor including spin degeneracy.
1839 : !max_occ = two / (self%nspinor * self%nsppol)
1840 : !fact_sigma = max_occ * (siemens_SI / Bohr_meter) / 100. ! Siemens / cm
1841 : ! fact_mob = fact_sigma * 100.**3 / e_Cb * Bohr_meter**3 ! cm^2 / V / s
1842 1 : fact_sbk = -volt_SI * kb_HaK ! Volt / Kelvin
1843 : !call wrtout(std_out,"factors: ")
1844 : ! write (std_out,*) fact_sigma, fact_sbk
1845 :
1846 : !call ibte%read_scattering()
1847 : ! Loops and memory are distributed over k-points and collinear spins
1848 11 : ABI_MALLOC(sr, (nkcalc, nsppol))
1849 1 : cnt = 0
1850 2 : do spin=1,nsppol
1851 8 : do ikcalc=1,nkcalc
1852 6 : cnt = cnt + 1
1853 7 : sr(ikcalc, spin)%rank = mod(cnt, nprocs)
1854 : end do
1855 : end do
1856 :
1857 1 : call cwtime(cpu, wall, gflops, "start")
1858 : ! Master reads and sends data to the rank treating (ikcalc, spin).
1859 1 : if (my_rank == master) then
1860 1 : NCF_CHECK(nctk_open_read(ncid, dtfil%filsigephin, xmpi_comm_self))
1861 : end if
1862 :
1863 2 : do spin=1,nsppol
1864 8 : do ikcalc=1,nkcalc
1865 6 : sr_p => sr(ikcalc, spin)
1866 6 : receiver = sr_p%rank
1867 6 : send_data = master /= receiver
1868 6 : if (.not. any(my_rank == [master, receiver])) cycle
1869 : !call wrtout(std_out, sjoin(" Sending data from my_rank:", itoa(my_rank), " to:", itoa(receiver)))
1870 :
1871 6 : if (my_rank == master) then
1872 : ! Get ncid of group used to store scattering rate for this k-point.
1873 6 : ncerr = nf90_inq_ncid(ncid, strcat("srate_k", itoa(ikcalc), "_s", itoa(spin)), grp_ncid)
1874 6 : if (ncerr /= NF90_NOERR) then
1875 0 : ABI_ERROR("Cannot find collision terms in SIGEPH file. Rerun eph_task -4 step with ibte_prep 1.")
1876 : end if
1877 6 : NCF_CHECK(nctk_get_dim(grp_ncid, "nq_ibzk_eff", sr_p%nq_ibzk_eff))
1878 6 : NCF_CHECK(nctk_get_dim(grp_ncid, "nbsum", nbsum))
1879 6 : NCF_CHECK(nctk_get_dim(grp_ncid, "nbcalc", nbcalc))
1880 6 : NCF_CHECK(nctk_get_dim(grp_ncid, "lgk_nsym", sr_p%lgk_nsym))
1881 30 : dims = [sr_p%nq_ibzk_eff, nbsum, nbcalc, sr_p%lgk_nsym]
1882 : end if
1883 :
1884 6 : if (send_data) then
1885 0 : tag = size(dims)
1886 0 : if (my_rank == master) call xmpi_send(dims, receiver, tag, comm, ierr)
1887 0 : if (my_rank == receiver) then
1888 0 : call xmpi_recv(dims, master, tag, comm, ierr)
1889 0 : sr_p%nq_ibzk_eff = dims(1); nbsum = dims(2); nbcalc = dims(3); sr_p%lgk_nsym = dims(4)
1890 : end if
1891 : end if
1892 :
1893 : ! Note that the size along the (n, m) axis does not depend on the kcalc index.
1894 18434 : ABI_CALLOC(sr_p%vals, (sr_p%nq_ibzk_eff, bmin:bmax, bmin:bmax, ntemp))
1895 18 : ABI_MALLOC(sr_p%kq_symtab, (6, sr_p%nq_ibzk_eff))
1896 18 : ABI_MALLOC(sr_p%lgk_sym2glob, (2, sr_p%lgk_nsym))
1897 :
1898 6 : if (my_rank == master) then
1899 6 : NCF_CHECK(nf90_get_var(grp_ncid, nctk_idname(grp_ncid, "kq_symtab"), sr_p%kq_symtab))
1900 6 : NCF_CHECK(nf90_get_var(grp_ncid, nctk_idname(grp_ncid, "lgk_sym2glob"), sr_p%lgk_sym2glob))
1901 :
1902 : ! Note that on file, we have:
1903 : !
1904 : ! nctkarr_t("srate", "dp", "nq_ibzk_eff, nbsum, nbcalc, ntemp")
1905 : !
1906 : ! but in terms of n, m indices we have that the:
1907 : !
1908 : ! n index: bstart_ks bstop_ks
1909 : ! m index: bsum_start up to bsum_stop to account for phonon emission/absorption.
1910 : !
1911 : ! so we have to insert the values in bmin:bmax slice.
1912 : ! TODO: Recheck this part.
1913 36 : ABI_MALLOC(grp_srate, (sr_p%nq_ibzk_eff, nbsum, nbcalc, ntemp))
1914 6 : NCF_CHECK(nf90_get_var(grp_ncid, nctk_idname(grp_ncid, "srate"), grp_srate))
1915 6 : ii = ibte%bstart_ks(ikcalc, spin)
1916 6 : jj = ibte%bstop_ks(ikcalc, spin)
1917 18404 : sr_p%vals(:, bmin:bmax, ii:jj, :) = grp_srate(:, 1:bmax-bmin+1, 1:nbcalc, :)
1918 6 : ABI_SFREE(grp_srate)
1919 : end if
1920 :
1921 6 : if (send_data) then
1922 0 : tag = size(sr_p%vals)
1923 0 : if (my_rank == master) then
1924 0 : tag = tag + 1; call xmpi_send(sr_p%vals, receiver, tag, comm, ierr)
1925 0 : tag = tag + 1; call xmpi_send(sr_p%kq_symtab, receiver, tag, comm, ierr)
1926 0 : tag = tag + 1; call xmpi_send(sr_p%lgk_sym2glob, receiver, tag, comm, ierr)
1927 : end if
1928 0 : if (my_rank == receiver) then
1929 0 : tag = tag + 1; call xmpi_recv(sr_p%vals, master, tag, comm, ierr)
1930 0 : tag = tag + 1; call xmpi_recv(sr_p%kq_symtab, master, tag, comm, ierr)
1931 0 : tag = tag + 1; call xmpi_recv(sr_p%lgk_sym2glob, master, tag, comm, ierr)
1932 : end if
1933 : end if
1934 :
1935 1 : if (send_data .and. my_rank /= receiver) call free_sr_ks(ikcalc, spin)
1936 :
1937 : end do ! spin
1938 : end do ! ikcalc
1939 :
1940 1 : if (my_rank == master) then
1941 1 : NCF_CHECK(nf90_close(ncid))
1942 : end if
1943 :
1944 1 : call cwtime_report(" sigeph IO", cpu, wall, gflops)
1945 :
1946 : ! Solve the linearized BTE with B = 0.
1947 : !
1948 : ! F_\nk = e df/de_\nk v_\nk \tau^0 + \tau^0 \sum_{mq} Srate_{nk,mq} F_{m,k+q}
1949 : !
1950 : ! where F is a vector in Cartesian coordinates and tau^0 is the SERTA relaxation time.
1951 : !
1952 : ! Take advantage of the following symmetry properties:
1953 : !
1954 : ! 1. F_k = F_Sk.
1955 : ! 2. F_{-k} = -F_k if TR symmetry.
1956 : ! 3. The q-space integration is reduced to the IBZ(k) using the symmetries of the little group of k.
1957 :
1958 : !call ibte%solve_ibte(solver_type=1)
1959 :
1960 1 : bsize = bmax - bmin + 1
1961 587 : ABI_CALLOC(fkn_in, (3, nkibz, bmin:bmax, nsppol))
1962 586 : ABI_CALLOC(fkn_out, (3, nkibz, bmin:bmax, nsppol))
1963 586 : ABI_CALLOC(fkn_serta, (3, nkibz, bmin:bmax, nsppol))
1964 586 : ABI_CALLOC(taukn_serta, (3, nkibz, bmin:bmax, nsppol))
1965 4 : ABI_MALLOC(ibte_sigma, (3, 3, 2, nsppol, ntemp))
1966 4 : ABI_MALLOC(ibte_seebeck, (3, 3, nsppol, ntemp))
1967 3 : ABI_MALLOC(ibte_pi, (3, 3, nsppol, ntemp))
1968 3 : ABI_MALLOC(ibte_mob, (3, 3, 2, nsppol, ntemp))
1969 3 : ABI_MALLOC(converged, (ntemp))
1970 3 : ABI_MALLOC(sig_gen, (3, 3, 2, nsppol))
1971 2 : ABI_MALLOC(sig_l21, (3, 3, 2, nsppol))
1972 2 : ABI_MALLOC(sig_l22, (3, 3, 2, nsppol))
1973 2 : ABI_MALLOC(mob_gen, (3, 3, 2, nsppol))
1974 2 : ABI_MALLOC(mob_21, (3, 3, 2, nsppol))
1975 2 : ABI_MALLOC(mob_22, (3, 3, 2, nsppol))
1976 3 : ABI_MALLOC(ibte_kappa, (3, 3, nsppol, ntemp))
1977 1 : abs_tol = dtset%ibte_abs_tol
1978 :
1979 : ! If the fermi level is inside the gap, F_k is gonna be very small
1980 : ! hence once should use a much smaller tolerance to converge.
1981 : ! According to numerical tests, a reasonable value of abs_tol can be estimated from the free carrier density using:
1982 : !
1983 : ! 1e-20 * e_density (in cm**-3)
1984 : !
1985 : ! These are the numerical values used to derive the fit:
1986 :
1987 : ! max_adiff = np.array([9e-12, 5.6e-6, 2.7e-6, 8.2e-60, 2.9e-46, 6.9e-6, 2.9e-8, 9.5E+00, 7.3E-04, 9.2E-04, 8.4E-04])
1988 : ! e_density = np.array([0.97e8, 0.86e13, 0.26e14, 0.1e-40, 0.89e-26, 0.45e14, 0.28e12, 0.10E+19, 0.12E+16, 0.10E+15, 0.90E+14])
1989 :
1990 1 : if (abs_tol <= zero) then
1991 0 : rtmp = minval(ibte%n_ehst, mask=ibte%n_ehst > zero) * ibte%nsppol
1992 0 : abs_tol = 1e-20 * rtmp / cryst%ucvol / Bohr_cm**3
1993 0 : call wrtout(std_out, " Input ibte_abs_tol <= zero ==> computing abs tolerance from minimal carrier density over all T")
1994 0 : call wrtout(std_out, " using: abs_tol = 1e-20 * e_density (in cm**-3)")
1995 0 : call wrtout(std_out, sjoin(" abs_tol:", ftoa(abs_tol), " from carrier_density:", ftoa(rtmp / cryst%ucvol / Bohr_cm**3)))
1996 : end if
1997 :
1998 1 : if (my_rank == master) then
1999 1 : path = strcat(dtfil%filnam_ds(4), "_RTA.nc")
2000 1 : call wrtout(units, ch10//sjoin("- Writing IBTE transport results to:", path))
2001 1 : NCF_CHECK(nctk_open_modify(ncid, path , xmpi_comm_self))
2002 :
2003 : ncerr = nctk_def_dims(ncid, [ &
2004 : nctkdim_t("nkibz", nkibz), nctkdim_t("bsize", bsize), nctkdim_t("nkcalc", ibte%nkcalc) &
2005 4 : ], defmode=.True.)
2006 1 : NCF_CHECK(ncerr)
2007 :
2008 : ncerr = nctk_def_arrays(ncid, [ &
2009 : nctkarr_t('fkn_out_sigma', "dp", "three, nkibz, bsize, nsppol, ntemp"), &
2010 : nctkarr_t('ibte_sigma', "dp", "three, three, two, nsppol, ntemp"), &
2011 : nctkarr_t('ibte_mob', "dp", "three, three, two, nsppol, ntemp"), &
2012 : nctkarr_t("kcalc2ibz", "int", "nkcalc, six"), &
2013 : nctkarr_t("kcalc2ebands", "int", "six, nkcalc"), &
2014 : nctkarr_t("kibz", "dp", "three, nkibz"), &
2015 : nctkarr_t('ibte_rho', "dp", "three, three, ntemp") &
2016 8 : ], defmode=.True.)
2017 1 : NCF_CHECK(ncerr)
2018 :
2019 1 : NCF_CHECK(nctk_set_datamode(ncid))
2020 : end if
2021 :
2022 1 : cnt = 0
2023 3 : do itemp=1,ntemp
2024 : !do itemp=ntemp, 1, -1
2025 2 : cnt = cnt + 1
2026 2 : kT = max(ibte%kTmesh(itemp), one * kb_HaK)
2027 2 : mu_e = ibte%eph_mu_e(itemp)
2028 2 : sig_p => ibte_sigma(:,:,:,:,itemp)
2029 2 : mob_p => ibte_mob(:,:,:,:,itemp)
2030 2 : sbk_p => ibte_seebeck(:,:,:,itemp)
2031 2 : kappa_p => ibte_kappa(:,:,:,itemp)
2032 2 : pi_p => ibte_pi(:,:,:,itemp)
2033 :
2034 2 : call wrtout(std_out," Value of kT: ", pre_newlines=1, newlines=1)
2035 2 : write (std_out,*) kT
2036 :
2037 :
2038 : ! Precompute tau_serta and fkn_serta for this T: f^'_nk v_\nk * \tau^0
2039 4 : do spin=1,nsppol
2040 16 : do ikcalc=1,nkcalc
2041 : !ik_ibz = ibte%kcalc2ibz(ikcalc, 1)
2042 12 : ik_ibz = ibte%kcalc2ebands(1, ikcalc)
2043 26 : do band_k=ibte%bstart_ks(ikcalc, spin), ibte%bstop_ks(ikcalc, spin)
2044 12 : lw_nk = ibte%linewidths(itemp, band_k, ik_ibz, spin, 1) ! 1 --> SERTA linewidths.
2045 12 : call safe_div(one, two * lw_nk, zero, tau_nk)
2046 48 : taukn_serta(:, ik_ibz, band_k, spin) = tau_nk
2047 12 : e_nk = ebands%eig(band_k, ik_ibz, spin)
2048 12 : dfde_nk = occ_dfde(e_nk, kT, mu_e)
2049 60 : fkn_serta(:, ik_ibz, band_k, spin) = tau_nk * dfde_nk * ibte%vbks(:, band_k, ik_ibz, spin)
2050 : end do
2051 : end do
2052 : end do
2053 :
2054 : call wrtout(std_out, sjoin(" Begin IBTE loop for itemp:", itoa(itemp), ", KT:", ftoa(kT / kb_HaK), "[K]"), &
2055 2 : pre_newlines=1, newlines=1)
2056 :
2057 : ! iter = 0 --> Compute SERTA transport tensors just for initial reference.
2058 2 : call ibte_calc_tensors(ibte, cryst, itemp, kT, mu_e, fkn_serta, onsager, sig_p, mob_p, fsum_eh, comm, iet)
2059 :
2060 : !TODO fix dtset here ans use maxocc, CAREFUL, it has to be included in all L coeff. I think
2061 : !max_occ = two / (self%nspinor * self%nsppol)
2062 : ! fact_sigma = max_occ * (siemens_SI / Bohr_meter) / 100. ! Siemens / cm
2063 : ! fact_mob = fact_sigma * 100.**3 / e_Cb * Bohr_meter**3 ! cm^2 / V / s
2064 : !fact_sbk = volt_SI * kb_HaK ! Volt / Kelvin
2065 : ! sig_p = fact_sigma * sig_p ! siemens cm^-1
2066 : ! mob_p = fact_mob * mob_p !
2067 :
2068 : ! Print mobility for semiconductors, conductivity for metals.
2069 2 : if (ibte%assume_gap) then
2070 0 : do spin=1,nsppol
2071 0 : mat33 = sum(mob_p(:,:,:,spin), dim=3)
2072 : write(msg, "(i5,1x,es9.1, *(1x, f16.2))") &
2073 0 : 0, zero, mat33(1,1), mat33(2,2), mat33(3,3), sum(fsum_eh(:,:,spin))
2074 : end do
2075 : else
2076 4 : do spin=1,nsppol
2077 62 : mat33 = sum(sig_p(:,:,:,spin), dim=3)
2078 : write(msg, "(i5,1x,es9.1, *(1x, es16.6))") &
2079 20 : 0, zero, mat33(1,1), mat33(2,2), mat33(3,3), sum(fsum_eh(:,:,spin))
2080 : end do
2081 : end if
2082 2 : call wrtout(std_out, msg)
2083 :
2084 1168 : fkn_in = fkn_serta
2085 : ! TODO: B-field
2086 : ! Initialize fkn_in either from SERTA or from previous T.
2087 : !if (cnt == 1) fkn_in = fkn_serta
2088 : !if (cnt > 1 ) fkn_in = fkn_out
2089 1166 : fkn_out = zero
2090 :
2091 2 : max_et=2 !max_et =1 pour F^E et max_et=2 pour F^E et F^T, iet = 3 pour calculer kappa
2092 :
2093 : ! start here with loop over E and T perturbations: once F^Eps is done, init F^T and converge it as well
2094 6 : electherm_loop: do iet = 1, max_et
2095 :
2096 4 : if (iet == 2) then
2097 : ! initialize F^T from F^E
2098 4 : do spin=1,nsppol
2099 16 : do ikcalc=1,nkcalc
2100 12 : ik_ibz = ibte%kcalc2ebands(1, ikcalc)
2101 26 : do band_k=ibte%bstart_ks(ikcalc, spin), ibte%bstop_ks(ikcalc, spin)
2102 12 : e_nk = ebands%eig(band_k, ik_ibz, spin)
2103 : !Which one of the following is correct ?
2104 : ! fkn_in(:, ik_ibz, band_k, spin) = fkn_serta(:, ik_ibz, band_k, spin) * (e_nk-mu_e)/kT
2105 : !fkn_in(:, ik_ibz, band_k, spin) = fkn_efield(:, ik_ibz, band_k, spin) * (e_nk-mu_e)/kT
2106 48 : fkn_in(:, ik_ibz, band_k, spin) = fkn_efield(:, ik_ibz, band_k, spin) * (e_nk-mu_e)
2107 60 : fkn_serta(:, ik_ibz, band_k, spin) = fkn_serta(:, ik_ibz, band_k, spin) * (e_nk-mu_e)
2108 : end do !band_k
2109 : end do !ikcalc
2110 : end do !spin
2111 2 : ABI_CHECK(converged(itemp), "Warning: E field BTE not converged, I will not try to converge the T gradient")
2112 2 : if (.not. converged(itemp)) exit electherm_loop
2113 : end if ! iet == 2 ?????
2114 :
2115 :
2116 : ! Begin iterative solver.
2117 25 : iter_loop: do iter=1,dtset%ibte_niter
2118 :
2119 : ! call wrtout(std_out," check beginning loop, iter, iet:", pre_newlines=1, newlines=1)
2120 : ! write (std_out,*) iter, iet
2121 :
2122 : ! Loop over the nk index in F_nk.
2123 50 : do spin=1,nsppol
2124 200 : do ikcalc=1,nkcalc
2125 150 : sr_p => sr(ikcalc, spin)
2126 150 : if (sr_p%rank /= my_rank) cycle ! MPI parallelism
2127 150 : ik_ibz = ibte%kcalc2ebands(1, ikcalc)
2128 325 : do band_k=ibte%bstart_ks(ikcalc, spin), ibte%bstop_ks(ikcalc, spin)
2129 :
2130 : ! Summing over the q-points in the effective IBZ(k) and the m band index.
2131 : ! Results stored in vec3. Integration weights are already included.
2132 150 : vec3 = zero
2133 300 : do band_sum=ibte%bmin, ibte%bmax
2134 229825 : do iq_sum=1, sr_p%nq_ibzk_eff
2135 229525 : ikq_ibz = sr_p%kq_symtab(1, iq_sum); isym_kq = sr_p%kq_symtab(2, iq_sum)
2136 229525 : trev_kq = sr_p%kq_symtab(6, iq_sum) !; g0_kq = sr_p%kq_symtab(3:5, iq_sum)
2137 : ! Build F_{m,k+q} in the effective IBZ(k) from fkn_in using symmetries (need k+q --> IBZ map)
2138 : ! Use transpose(R) because we are using the tables for the wavefunctions
2139 : ! In this case listkk has been called with symrec and use_symrec=False
2140 : ! so q_bz = S^T q_ibz where S is the isym_kq symmetry
2141 : ! vkq = matmul(transpose(cryst%symrel_cart(:,:,isym_kq)), vkq)
2142 2983825 : mat33 = transpose(cryst%symrel_cart(:,:,isym_kq))
2143 2983825 : f_kq = matmul(mat33, fkn_in(:, ikq_ibz, band_sum, spin))
2144 229525 : if (trev_kq == 1) f_kq = -f_kq
2145 918250 : vec3 = vec3 + sr_p%vals(iq_sum, band_sum, band_k, itemp) * f_kq(:)
2146 : end do ! iq_sum
2147 : end do ! band_k
2148 :
2149 : ! Symmetrize intermediate results using the operations of the little group of k.
2150 150 : sym_vec = zero
2151 800 : do isym_lgk=1,sr_p%lgk_nsym
2152 650 : isym = sr_p%lgk_sym2glob(1, isym_lgk)
2153 650 : itime = sr_p%lgk_sym2glob(2, isym_lgk)
2154 8450 : mat33 = transpose(cryst%symrel_cart(:,:,isym))
2155 : !if(itime == 1) mat33 = -mat33 ! FIXME: here there's a different convention for TR used in m_lgroup
2156 4550 : if (itime == 2) mat33 = -mat33
2157 10550 : sym_vec = sym_vec + matmul(mat33, vec3)
2158 : end do
2159 600 : sym_vec = taukn_serta(:, ik_ibz, band_k, spin) * sym_vec / sr_p%lgk_nsym
2160 : ! if iet 2 serta and fkn_in have been multiplied by e-mu / T
2161 750 : fkn_out(:, ik_ibz, band_k, spin) = fkn_serta(:, ik_ibz, band_k, spin) + sym_vec
2162 : end do ! band_k
2163 : end do ! ikcalc
2164 : end do ! spin
2165 :
2166 :
2167 25 : call xmpi_sum(fkn_out, comm, ierr)
2168 :
2169 : ! Write fkn_out_sigma to disk.
2170 25 : if (my_rank == master) then
2171 150 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "fkn_out_sigma"), fkn_out, start=[1,1,1,1,itemp]))
2172 : end if
2173 :
2174 50 : do spin=1,nsppol
2175 14575 : max_adiff_spin(spin) = maxval(abs(fkn_out(:,:,:,spin) - fkn_in(:,:,:,spin)))
2176 : end do
2177 75 : max_adiff = maxval(max_adiff_spin)
2178 :
2179 :
2180 725 : mob_gen=sig_gen
2181 725 : mob_21=sig_gen
2182 725 : mob_22=sig_gen
2183 :
2184 : !Here I test if F^T=(eps-mu)F^E
2185 : ! if (iet==2) then
2186 : ! do spin=1,nsppol
2187 : ! do ikcalc=1,nkcalc
2188 : ! ik_ibz = ibte%kcalc2ebands(1, ikcalc)
2189 : ! do band_k=ibte%bstart_ks(ikcalc, spin), ibte%bstop_ks(ikcalc, spin)
2190 : ! e_nk = ebands%eig(band_k, ik_ibz, spin)
2191 : !Which one of the following is correct ?
2192 : ! fkn_in(:, ik_ibz, band_k, spin) = fkn_serta(:, ik_ibz, band_k, spin) * (e_nk-mu_e)/kT
2193 : !fkn_in(:, ik_ibz, band_k, spin) = fkn_efield(:, ik_ibz, band_k, spin) * (e_nk-mu_e)/kT
2194 : ! fkn_out(:, ik_ibz, band_k, spin) = fkn_efield(:, ik_ibz, band_k, spin) * (e_nk-mu_e)
2195 : !fkn_serta(:, ik_ibz, band_k, spin) = fkn_serta(:, ik_ibz, band_k, spin) * (e_nk-mu_e)
2196 : ! end do
2197 : ! end do
2198 : ! end do
2199 : ! end if
2200 :
2201 :
2202 : ! Compute transport tensors from fkn_out (= F_eps or F_T)
2203 25 : if(iet==1) call ibte_calc_tensors(ibte, cryst, itemp, kT, mu_e, fkn_out, onsager, sig_p, mob_p, fsum_eh, comm, iet)
2204 : ! call flush_unit(std_out)
2205 :
2206 : !For seebeck, use another array(sig_gen) because sig_p points towards ibte_sigma
2207 25 : if(iet==2) call ibte_calc_tensors(ibte, cryst, itemp, kT, mu_e, fkn_out, onsager, sig_gen, mob_gen, fsum_eh, comm, iet)
2208 :
2209 : !check to find the bug
2210 : ! call wrtout(std_out," check 3", pre_newlines=1, newlines=1)
2211 : ! call flush_unit(std_out)
2212 :
2213 : ! Print mobility for semiconductors or conductivity for metals.
2214 25 : if (iet==1) then! for F_eps and charge transport
2215 : !sig_p = fact_sigma * sig_gen ! siemens cm^-1
2216 : !mob_p = fact_mob * mob_gen ! cm2/V/s
2217 :
2218 : !call wrtout(std_out, "check 3bis: sig_gen, fact_sigma, sig_p")
2219 : ! write(std_out,*) , sig_gen, fact_sigma, sig_p
2220 : !check to find the bug
2221 : ! call wrtout(std_out," check 4", pre_newlines=1, newlines=1)
2222 :
2223 23 : if (ibte%assume_gap) then
2224 0 : do spin=1,nsppol
2225 0 : mat33 = sum(mob_p(:,:,:,spin), dim=3)
2226 : write(msg, "(i5,1x,es9.1,*(1x, f16.2))") &
2227 0 : iter, max_adiff_spin(spin), mat33(1,1), mat33(2,2), mat33(3,3), sum(fsum_eh(:,:,spin))
2228 : end do
2229 : else
2230 46 : do spin=1,nsppol
2231 713 : mat33 = sum(sig_p(:,:,:,spin), dim=3)
2232 : write(msg, "(i5,1x,es9.1,*(1x, es16.6))") &
2233 230 : iter, max_adiff_spin(spin), mat33(1,1), mat33(2,2), mat33(3,3), sum(fsum_eh(:,:,spin))
2234 : end do
2235 : end if
2236 23 : call wrtout(std_out, msg)
2237 :
2238 2 : else if (iet == 2) then! for Seebeck, Seebeck=1/T sig_p^-1 * sig_gen
2239 : ! output sig_p^-1 * sig_gen
2240 : !pay attention when implementing if sigma=zero, bug ! impossible to get S.
2241 : ! use matr3inv to inverse the matrix sigma but verify also that the det is not equal to zero !
2242 : ! maybe type the command use "..." to be able to use the functions, check in the code
2243 :
2244 : !call wrtout(std_out, "check : sig_gen, sig_p")
2245 : !write(std_out,*) sig_gen, sig_p
2246 :
2247 : !TODO: implement for inv_sig_p with two spins
2248 4 : do spin=1,nsppol
2249 62 : call inv33(sum(sig_p(:,:,:,spin),dim=3), inv_sig_p)
2250 : !We divide by 100 because sig_p is in Siemens cm^-1 and we want to retrieve meters
2251 : ! work33 = matmul (inv_sig_p, sum(sig_gen(:,:,:,spin),dim=4))
2252 164 : mat33 = matmul (inv_sig_p, sum(sig_gen(:,:,:,spin),dim=3)) / (kT)
2253 : ! mat33 = (volt_SI / (ibte%kTmesh(itemp) / kb_HaK)) * matmul (inv_sig_p, sum(sig_gen(:,:,:,spin),dim=3))
2254 : write(msg, "(i5,1x,es9.1,*(1x, es16.6))") &
2255 18 : iter, max_adiff_spin(spin), mat33(1,1), mat33(2,2), mat33(3,3), sum(fsum_eh(:,:,spin))
2256 28 : sbk_p (:,:,spin) = mat33
2257 : end do
2258 2 : call wrtout(std_out, msg)
2259 : end if ! iet == 1, 2
2260 :
2261 25 : if (ibte%assume_gap) then
2262 0 : write(msg, "(a5,1x,a9,*(1x, a16))")" ITER", "max_adiff", "mobility_e+h", "sum_k(df_k)"
2263 : else
2264 25 : write(msg, "(a5,1x,a9,*(1x, a16))")" ITER", "max_adiff", "conductivity", "sum_k(df_k)"
2265 : end if
2266 25 : call wrtout(std_out, msg)
2267 :
2268 :
2269 : ! Check for convergence by testing max_k |F_k^i - F_k^{i-1}|.
2270 25 : call wrtout (std_out, "Print check convergence: max_adiff, abs_tol")
2271 25 : write(std_out,*) max_adiff, abs_tol
2272 25 : converged(itemp) = max_adiff < abs_tol
2273 50 : if (converged(itemp)) then
2274 : call wrtout(std_out, sjoin(" IBTE solver converged after:", itoa(iter), &
2275 4 : "iterations within ibte_abs_tol:", ftoa(abs_tol)), pre_newlines=1)
2276 :
2277 4 : exit iter_loop
2278 : else
2279 : ! Linear mixing of fkn_in and fkn_out.
2280 12264 : fkn_in = (one - dtset%ibte_alpha_mix) * fkn_in + dtset%ibte_alpha_mix * fkn_out
2281 12243 : fkn_out = zero
2282 : end if
2283 : !check where is the problem with iter loop
2284 :
2285 : !call wrtout(std_out," check iet:", pre_newlines=1, newlines=1)
2286 : !write (std_out,*) iet
2287 :
2288 : end do iter_loop
2289 :
2290 4 : if (.not. converged(itemp)) then
2291 0 : msg = sjoin("Not converged after:", itoa(dtset%ibte_niter), "max iterations")
2292 0 : call wrtout(ab_out, msg, pre_newlines=1, newlines=1)
2293 0 : ABI_WARNING(msg)
2294 : end if
2295 :
2296 : !Change of unit here because sig_p is used in to calculate sbk in SI units
2297 1174 : if (iet == 1) fkn_efield = fkn_out
2298 6 : if (iet==2) then
2299 : !iet =3
2300 2 : call ibte_calc_tensors(ibte, cryst, itemp, kT, mu_e, fkn_efield, onsager, sig_l21, mob_21, fsum_eh, comm, 3)
2301 2 : call ibte_calc_tensors(ibte, cryst, itemp, kT, mu_e, fkn_out, onsager, sig_l22, mob_22, fsum_eh, comm, 3)
2302 4 : do spin=1, nsppol
2303 62 : call inv33(sum(sig_p(:,:,:,spin),dim=3), inv_sig_p)
2304 252 : kappa_p(:,:,spin)=(1/kT)*(sum(sig_l22(:,:,:,spin), dim=3)-matmul(sum(sig_l21(:,:,:,spin),dim=3),sbk_p(:,:,spin)*kT))
2305 166 : pi_p(:,:,spin)= matmul (sum(sig_l21(:,:,:,spin),dim=3), inv_sig_p) ! PI=L21/L11 in IBTE
2306 : end do
2307 : end if
2308 : end do electherm_loop
2309 : ! end E/T loop
2310 :
2311 : !Units conversion for sigma, at this step in order to keep atomic units for the calculation of S
2312 : !Truncated in purpose for nspinor=2 and nsppol=1 compared to how its done in the calc_tensor routine because dont have access to self in this routine
2313 57 : sig_p=sig_p*(siemens_SI / Bohr_meter / cryst%ucvol) / 100
2314 : !check if sig_p and sbk_p are empty
2315 : !call wrtout(std_out, "sig_p, sbk_p")
2316 : !write(std_out,*) sig_p, sbk_p
2317 :
2318 : end do ! itemp
2319 :
2320 : !TODO: do it for spins as well
2321 3 : ABI_MALLOC(ibte_rho, (3, 3, ntemp))
2322 3 : do itemp=1,ntemp
2323 62 : work33 = sum(ibte_sigma(:,:,:,1,itemp), dim=3)
2324 2 : if (ibte%nsppol == 2) work33 = work33 + sum(ibte_sigma(:,:,:,2,itemp), dim=3)
2325 2 : call inv33(work33, mat33)
2326 27 : ibte_rho(:, :, itemp) = 1e+6_dp * mat33
2327 : end do
2328 :
2329 : !ibte_sigma rho and seebeck seem empty, let's check
2330 :
2331 :
2332 : !write(*,*) "sigma: ", ibte_sigma
2333 : !write(*,*) "rho: ", ibte_rho
2334 : !write(*,*) "seebeck: ", ibte_seebeck
2335 :
2336 :
2337 1 : if (my_rank == master) then
2338 : ! Write final results to main output.
2339 4 : components = ["xx", "yy", "zz"]
2340 1 : if (ibte%assume_gap) then
2341 : ! SemiConductor
2342 0 : do ii=1,3
2343 0 : call wrtout(units, sjoin(" Cartesian component of IBTE mobility tensor:", components(ii)))
2344 0 : write(msg, "(a16,2(a32),a16)") 'Temperature [K]', 'e/h density [cm^-3]', 'e/h mobility [cm^2/Vs]', "Converged"
2345 0 : call wrtout(units, msg)
2346 :
2347 0 : do spin=1,ibte%nsppol
2348 0 : if (ibte%nsppol == 2) call wrtout(units, sjoin(" For spin:", stoa(spin)), newlines=1)
2349 :
2350 0 : do itemp=1,ibte%ntemp
2351 : write(msg,"(f16.2,2e16.2,2f16.2,a16)") &
2352 0 : ibte%kTmesh(itemp) / kb_HaK, &
2353 0 : ibte%n_ehst(1, spin, itemp) / cryst%ucvol / Bohr_cm **3, &
2354 0 : ibte%n_ehst(2, spin, itemp) / cryst%ucvol / Bohr_cm **3, &
2355 0 : ibte_mob(ii, ii, 1, spin, itemp), ibte_mob(ii, ii, 2, spin, itemp), &
2356 0 : yesno(converged(itemp))
2357 0 : call wrtout(units, msg)
2358 : end do ! itemp
2359 : end do ! spin
2360 : ! TODO:HERE look into adding seebeck output for semiconductors as well, check what happens with e and h components in tensors
2361 : ! routine
2362 0 : call wrtout(units, ch10)
2363 : end do ! ii
2364 :
2365 : else
2366 : ! Metals. Print conductivity (spin resolved) and resistivity (no spin resolved)
2367 3 : do ii=1,2
2368 2 : if (ii == 1) msg = " Conductivity [Siemens cm^-1] using IBTE"
2369 2 : if (ii == 2) msg = " Resistivity [micro-Ohm cm] using IBTE"
2370 2 : call wrtout(units, msg)
2371 :
2372 2 : nsp = ibte%nsppol; if (ii == 2) nsp = 1
2373 4 : do spin=1,nsp
2374 2 : if (ibte%nsppol == 2) call wrtout(units, sjoin(" For spin:", stoa(spin)), newlines=1)
2375 2 : write(msg, "(5a16)") 'Temperature (K)', 'xx', 'yy', 'zz', "Converged"
2376 2 : call wrtout(units, msg)
2377 8 : do itemp=1,ibte%ntemp
2378 4 : if (ii == 1) then
2379 62 : mat33 = sum(ibte_sigma(:,:,:,spin,itemp), dim=3)
2380 : else
2381 26 : mat33 = ibte_rho(:,:,itemp)
2382 : end if
2383 : write(msg,"(f16.2,3e16.6,a16)") &
2384 4 : ibte%kTmesh(itemp) / kb_HaK, mat33(1,1), mat33(2,2), mat33(3,3), yesno(converged(itemp))
2385 6 : call wrtout(units, msg)
2386 : end do ! itemp
2387 : end do ! spin
2388 3 : call wrtout(units, ch10)
2389 : end do ! ii
2390 :
2391 : ! HERE add output of Seebeck and kappa_el coefficients
2392 : ! TODO: add off diagonal Seebeck and sigma coefficients
2393 1 : msg = " Seebeck [Volts / Kelvin] using IBTE"
2394 1 : call wrtout(units, msg)
2395 2 : do spin=1,ibte%nsppol
2396 1 : if (ibte%nsppol == 2) call wrtout(units, sjoin(" For spin:", stoa(spin)), newlines=1)
2397 1 : write(msg, "(5a16)") 'Temperature (K)', 'xx', 'yy', 'zz', "Converged"
2398 1 : call wrtout(units, msg)
2399 4 : do itemp=1,ibte%ntemp
2400 26 : mat33 = fact_sbk * ibte_seebeck(:,:,spin,itemp)
2401 : write(msg,"(f16.2,3e16.6,a16)") &
2402 2 : ibte%kTmesh(itemp) / kb_HaK, mat33(1,1), mat33(2,2), mat33(3,3), yesno(converged(itemp))
2403 3 : call wrtout(units, msg)
2404 : end do ! itemp
2405 : end do ! spin
2406 1 : call wrtout(units, ch10)
2407 :
2408 :
2409 1 : msg = " Kappa [W/m*K] using IBTE"
2410 1 : call wrtout(units, msg)
2411 2 : do spin=1,ibte%nsppol
2412 1 : if (ibte%nsppol == 2) call wrtout(units, sjoin(" For spin:", stoa(spin)), newlines=1)
2413 1 : write(msg, "(5a16)") 'Temperature (K)', 'xx', 'yy', 'zz', "Converged"
2414 1 : call wrtout(units, msg)
2415 4 : do itemp=1,ibte%ntemp
2416 26 : mat33 = volt_SI**2 * kb_HaK * (siemens_SI / Bohr_meter / cryst%ucvol)* ibte_kappa(:,:,spin,itemp)
2417 : write(msg,"(f16.2,3e16.6,a16)") &
2418 2 : ibte%kTmesh(itemp) / kb_HaK, mat33(1,1), mat33(2,2), mat33(3,3), yesno(converged(itemp))
2419 3 : call wrtout(units, msg)
2420 : end do ! itemp
2421 : end do ! spin
2422 1 : call wrtout(units, ch10)
2423 :
2424 1 : msg = " Peltier [Volts] using IBTE"
2425 1 : call wrtout(units, msg)
2426 2 : do spin=1,ibte%nsppol
2427 1 : if (ibte%nsppol == 2) call wrtout(units, sjoin(" For spin:", stoa(spin)), newlines=1)
2428 1 : write(msg, "(5a16)") 'Temperature (K)', 'xx', 'yy', 'zz', "Converged"
2429 1 : call wrtout(units, msg)
2430 4 : do itemp=1,ibte%ntemp
2431 26 : mat33 = -volt_SI * ibte_pi(:,:,spin,itemp)
2432 : write(msg,"(f16.2,3e16.6,a16)") &
2433 2 : ibte%kTmesh(itemp) / kb_HaK, mat33(1,1), mat33(2,2), mat33(3,3), yesno(converged(itemp))
2434 3 : call wrtout(units, msg)
2435 : end do ! itemp
2436 : end do ! spin
2437 1 : call wrtout(units, ch10)
2438 :
2439 1 : msg = "Carrier density: "
2440 1 : call wrtout(units, msg)
2441 1 : write(msg, "(a16,a16,a32)") 'Temperature [K]', ' chem pot [eV] ', 'e/h density [cm^-3]'
2442 1 : call wrtout(units, msg)
2443 :
2444 2 : do spin=1,ibte%nsppol
2445 1 : if (ibte%nsppol == 2) call wrtout(units, sjoin(" For spin:", stoa(spin)), newlines=1)
2446 :
2447 4 : do itemp=1,ibte%ntemp
2448 : write(msg,"(f16.2,3e16.2)") &
2449 2 : ibte%kTmesh(itemp) / kb_HaK, &
2450 2 : ibte%eph_mu_e(itemp) * eV_Ha, &
2451 2 : ibte%n_ehst(1, spin, itemp) / cryst%ucvol / Bohr_cm **3, &
2452 4 : ibte%n_ehst(2, spin, itemp) / cryst%ucvol / Bohr_cm **3
2453 3 : call wrtout(units, msg)
2454 : end do ! itemp
2455 : end do ! spin
2456 1 : call wrtout(units, ch10)
2457 :
2458 :
2459 :
2460 :
2461 : end if
2462 :
2463 : !pre = "_IBTE"
2464 : !call ibte%write_tensor(dtset, irta, "sigma", ibte%sigma(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_SIGMA"))
2465 : !call ibte%write_tensor(dtset, irta, "seebeck", ibte%seebeck(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_SBK"))
2466 : !call ibte%write_tensor(dtset, irta, "kappa", ibte%kappa(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_KAPPA"))
2467 : !call ibte%write_tensor(dtset, irta, "zte", ibte%zte(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_ZTE"))
2468 : !call ibte%write_tensor(dtset, irta, "pi", ibte%pi(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_PI"))
2469 :
2470 : ! Print IBTE results to stdout and other external txt files (for the test suite)
2471 : !call ibte%print_rta_txt_files(cryst, dtset, dtfil)
2472 : ! Creates the netcdf file used to store the results of the calculation.
2473 1 : path = strcat(dtfil%filnam_ds(4), "_RTA.nc")
2474 1 : call wrtout(units, ch10//sjoin("- Writing IBTE transport results to:", path))
2475 : !NCF_CHECK(nctk_open_modify(ncid, path , xmpi_comm_self))
2476 :
2477 : ncerr = nctk_def_arrays(ncid, [ &
2478 : nctkarr_t('ibte_sigma', "dp", "three, three, two, nsppol, ntemp"), &
2479 : nctkarr_t('ibte_mob', "dp", "three, three, two, nsppol, ntemp"), &
2480 : nctkarr_t('ibte_rho', "dp", "three, three, ntemp"), &
2481 : nctkarr_t('ibte_kappa', "dp", "three, three, nsppol, ntemp"), &
2482 : nctkarr_t('ibte_seebeck', "dp", "three, three, nsppol, ntemp"), &
2483 : nctkarr_t('ibte_pi', "dp", "three, three, nsppol, ntemp") &
2484 7 : ], defmode=.True.)
2485 1 : NCF_CHECK(ncerr)
2486 :
2487 : ! Write data.
2488 1 : NCF_CHECK(nctk_set_datamode(ncid))
2489 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ibte_sigma"), ibte_sigma))
2490 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ibte_mob"), ibte_mob))
2491 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ibte_rho"), ibte_rho))
2492 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ibte_seebeck"), ibte_seebeck))
2493 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ibte_pi"), ibte_pi))
2494 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ibte_kappa"), ibte_kappa))
2495 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kcalc2ibz"), ibte%kcalc2ibz))
2496 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kcalc2ebands"), ibte%kcalc2ebands))
2497 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kibz"), ibte%ebands%kptns))
2498 1 : NCF_CHECK(nf90_close(ncid))
2499 : end if ! master
2500 :
2501 : ! Free memory
2502 1 : ABI_FREE(fkn_serta)
2503 1 : ABI_FREE(taukn_serta)
2504 1 : ABI_FREE(fkn_in)
2505 1 : ABI_FREE(fkn_out)
2506 1 : ABI_FREE(ibte_sigma)
2507 1 : ABI_FREE(ibte_seebeck)
2508 1 : ABI_FREE(ibte_pi)
2509 1 : ABI_FREE(ibte_mob)
2510 1 : ABI_FREE(ibte_rho)
2511 1 : ABI_FREE(ibte_kappa)
2512 1 : ABI_FREE(converged)
2513 1 : ABI_FREE(sig_gen)
2514 1 : ABI_FREE(sig_l21)
2515 1 : ABI_FREE(sig_l22)
2516 1 : ABI_FREE(mob_gen)
2517 1 : ABI_FREE(mob_21)
2518 1 : ABI_FREE(mob_22)
2519 :
2520 2 : do spin=1,nsppol
2521 8 : do ikcalc=1,nkcalc
2522 7 : call free_sr_ks(ikcalc, spin)
2523 : end do
2524 : end do
2525 7 : ABI_FREE(sr)
2526 :
2527 4 : call ibte%free()
2528 :
2529 : contains
2530 :
2531 6 : subroutine free_sr_ks(ikc, isp)
2532 : integer,intent(in) :: ikc, isp
2533 6 : ABI_SFREE(sr(ikc, isp)%vals)
2534 6 : ABI_SFREE(sr(ikc, isp)%lgk_sym2glob)
2535 6 : ABI_SFREE(sr(ikc, isp)%kq_symtab)
2536 6 : end subroutine free_sr_ks
2537 :
2538 : end subroutine ibte_driver
2539 : !!***
2540 :
2541 : !----------------------------------------------------------------------
2542 :
2543 : !!****f* m_rta/ibte_calc_tensors
2544 : !! NAME
2545 : !! ibte_calc_tensors
2546 : !!
2547 : !! FUNCTION
2548 : !! calculate transport tensors within iBTE method, from v x F expressions for current and sigma etc
2549 : !! this routine now stays in atomic units to accommodate calculation of sigma and seebeck
2550 : !!
2551 : !! INPUTS
2552 : !! cryst<crystal_t>=Crystalline structure
2553 : !! comm=MPI communicator.
2554 : !!
2555 : !! SOURCE
2556 :
2557 31 : subroutine ibte_calc_tensors(self, cryst, itemp, kT, mu_e, fk, onsager, sigma_eh, mob_eh, fsum_eh, comm, iet)
2558 :
2559 : !Arguments ------------------------------------
2560 : class(rta_t),intent(inout) :: self
2561 : type(crystal_t),intent(in) :: cryst
2562 : integer,intent(in) :: itemp
2563 : real(dp),intent(in) :: kT, mu_e
2564 : real(dp),intent(in) :: fk(3, self%ebands%nkpt, self%bmin:self%bmax, self%nsppol)
2565 : real(dp),intent(out) :: sigma_eh(3,3,2,self%nsppol), mob_eh(3,3,2,self%nsppol)
2566 : real(dp),intent(out) :: fsum_eh(3,2,self%nsppol), onsager(3,3,3,self%nsppol)
2567 : integer,intent(in) :: comm, iet !iet to know which tensor we calculate, 1 for L11, 2 for L12 and 3 for L21 and L22
2568 :
2569 : !Local variables ------------------------------
2570 : !scalars
2571 : integer :: nsppol, nkibz, ib, ik_ibz, spin, ii, jj, ieh, cnt, nprocs, ia, time_opt
2572 : real(dp) :: eig_nk, max_occ, wtk, emu_alpha, fact, fact0
2573 : !real(dp) :: fact_sigma, fact_mob
2574 : !arrays
2575 : real(dp) :: vr(3), vv_tens(3,3)
2576 : !************************************************************************
2577 :
2578 : ABI_UNUSED(kt)
2579 :
2580 31 : nprocs = xmpi_comm_size(comm)
2581 :
2582 : ! Copy important dimensions
2583 31 : nkibz = self%ebands%nkpt; nsppol = self%ebands%nsppol
2584 31 : time_opt = 0 ! This to preserve the previous behaviour in which TR was not used.
2585 :
2586 : ! sigma_IBTE = (-S e^ / omega sum_\nk) (v_\nk \otimes F_\nk)
2587 : ! with S the spin degeneracy factor.
2588 2387 : sigma_eh = zero; fsum_eh = zero; onsager = zero
2589 :
2590 : ! Compute mobility_mu i.e. results in which lifetimes have been computed in a consistent way
2591 : ! with the same the Fermi level. In all the other cases, indeed, we assume that tau does not depend on ef.
2592 : !
2593 :
2594 : !TODO: rewrite this because every L (onsager coeff) is given with a minus sign and it's just the case for L11 and L22 normally
2595 : ! Fortunately these minus signs compensate each other in the code or are suppressed by putting a -1 factor in front of sbk and PI
2596 : ! It is not because of the sign of e !
2597 : cnt = 0
2598 62 : do spin=1,nsppol
2599 4557 : do ik_ibz=1,nkibz
2600 : !cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! MPI parallelism.
2601 4495 : wtk = self%ebands%wtk(ik_ibz)
2602 :
2603 9021 : do ib=self%bmin,self%bmax
2604 4495 : eig_nk = self%ebands%eig(ib, ik_ibz, spin)
2605 :
2606 : ! Compute outer product in vv_tens and symmetrize tensor.
2607 17980 : vr(:) = self%vbks(:, ib, ik_ibz, spin)
2608 22475 : do ia=1,3
2609 13485 : if (ia == 1) then
2610 : emu_alpha = one
2611 : else
2612 8990 : emu_alpha = (eig_nk - mu_e) ** (ia - 1)
2613 : end if
2614 :
2615 53940 : do ii=1,3
2616 175305 : do jj=1,3
2617 161820 : vv_tens(ii, jj) = vr(ii) * fk(jj, ik_ibz, ib, spin)
2618 : end do
2619 : end do
2620 175305 : vv_tens = vv_tens * emu_alpha
2621 175305 : vv_tens = cryst%symmetrize_cart_tens33(vv_tens, time_opt)
2622 : ! TODO: write the expression below in another way, without the useless sigma_eh - and without the sign - (adapt the
2623 : ! calculations of coefficients in the main routine then)
2624 13485 : if (ia == 1 .and. iet /= 3) then
2625 3770 : ieh = 2; if (eig_nk >= mu_e) ieh = 1
2626 49010 : sigma_eh(:,:,ieh,spin) = sigma_eh(:,:,ieh,spin) - wtk * vv_tens
2627 15080 : fsum_eh(:,ieh,spin) = fsum_eh(:,ieh,spin) + wtk * cryst%symmetrize_cart_vec3(fk(:, ik_ibz, ib, spin), time_opt)
2628 :
2629 : !Here the tensorial product in order to get L21 and L22 to calculate kappa
2630 9715 : else if (ia==2 .and. iet==3) then
2631 725 : ieh = 2; if (eig_nk >= mu_e) ieh = 1
2632 9425 : sigma_eh(:,:,ieh,spin) = sigma_eh(:,:,ieh,spin) - wtk * vv_tens
2633 2900 : fsum_eh(:,ieh,spin) = fsum_eh(:,ieh,spin) + wtk * cryst%symmetrize_cart_vec3(fk(:, ik_ibz, ib, spin), time_opt)
2634 : end if
2635 :
2636 :
2637 :
2638 179800 : onsager(:,:,ia,spin) = onsager(:,:,ia,spin) - wtk * vv_tens
2639 : end do ! ia
2640 :
2641 : end do ! ib
2642 : end do ! ik_ibz
2643 : end do ! spin
2644 :
2645 31 : if (iet==1) then ! In order to output the mobility (only for L11, the conductivity, thus iet==1)
2646 23 : max_occ = two / (self%nspinor * self%nsppol)
2647 23 : fact0 = max_occ * (siemens_SI / Bohr_meter / cryst%ucvol) / 100
2648 23 : fact = 100**3 / e_Cb
2649 :
2650 644 : sigma_eh = fact0 * sigma_eh ! siemens cm^-1
2651 230 : fsum_eh = fsum_eh / cryst%ucvol
2652 :
2653 : ! Scale by the carrier concentration.
2654 46 : do spin=1,nsppol
2655 92 : do ieh=1,2
2656 : call safe_div(sigma_eh(:,:,ieh,spin) * fact, &
2657 621 : self%n_ehst(ieh, spin, itemp) / cryst%ucvol / Bohr_meter**3, zero, mob_eh(:,:,ieh,spin))
2658 : end do
2659 : end do
2660 :
2661 : !Here I rescale sigma_eh to output correctly the Onsager coeff L11
2662 644 : sigma_eh = sigma_eh / fact0
2663 : end if
2664 : !call xmpi_sum(sigma_eh, comm, ierr)
2665 : !call xmpi_sum(onsager, comm, ierr)
2666 :
2667 : ! max_occ = two / (self%nspinor * self%nsppol)
2668 : ! sigma_eh = max_occ * sigma_eh / cryst%ucvol
2669 : !fsum_eh = fsum_eh / cryst%ucvol
2670 :
2671 31 : max_occ = two / (self%nspinor * self%nsppol)
2672 : !Take into account spin degeneracy for all Onsager coefficient:
2673 868 : sigma_eh= max_occ * sigma_eh
2674 :
2675 :
2676 :
2677 : ! fact0 = max_occ * (siemens_SI / Bohr_meter / cryst%ucvol) / 100
2678 : ! fact = 100**3 / e_Cb
2679 : ! sigma_eh = fact0 * sigma_eh
2680 :
2681 :
2682 :
2683 : ! Scale by the carrier concentration.
2684 : ! do spin=1,nsppol
2685 : ! do ieh=1,2
2686 : ! call safe_div(sigma_eh(:,:,ieh,spin) * fact, &
2687 : ! self%n_ehst(ieh, spin, itemp) / cryst%ucvol / Bohr_meter**3, zero, mob_eh(:,:,ieh,spin))
2688 :
2689 : !check to find the bug
2690 : !call wrtout(std_out, " ieh, shape(sigma_eh) vol shape(mob_eh) ", pre_newlines=1, newlines=1)
2691 : !write (std_out,*) ieh, shape(sigma_eh), cryst%ucvol, shape(mob_eh)
2692 : !call wrtout(std_out," check safe_div, in order: sigma_eh, self%n_ehst, cryst%ucvol, zero", pre_newlines=1, newlines=1)
2693 : ! write (std_out,*) sigma_eh(:,:,ieh,spin), self%n_ehst(ieh, spin, itemp), cryst%ucvol, zero
2694 : ! call flush_unit(std_out)
2695 :
2696 : ! call safe_div(sigma_eh(:,:,ieh,spin), &
2697 : ! self%n_ehst(ieh, spin, itemp) / cryst%ucvol, zero, mob_eh(:,:,ieh,spin))
2698 : !mob_eh(:,:,ieh,spin) = sigma_eh(:,:,ieh,spin) / (self%n_ehst(ieh, spin, itemp) / cryst%ucvol)
2699 : !check to find the bug
2700 : ! call wrtout(std_out," check safe_div, mob_eh: ", pre_newlines=1, newlines=1)
2701 : ! write (std_out,*) mob_eh(:,:,ieh,spin)
2702 : ! call flush_unit(std_out)
2703 :
2704 : ! end do
2705 : ! end do
2706 : !In order to be able to compile, to have a value for the dummy argument mob_eh
2707 :
2708 : !call wrtout(std_out, "shape(sigma_eh), sigma_eh, shape(mob_eh), mob_eh ", pre_newlines=1, newlines=1)
2709 : !write (std_out,*) shape(sigma_eh),sigma_eh, shape(mob_eh), mob_eh
2710 :
2711 : !mob_eh=zero
2712 :
2713 : !mob_eh=sigma_eh
2714 :
2715 :
2716 : !max_occ = two / (self%nspinor * self%nsppol)
2717 : ! fact_sigma = max_occ * (siemens_SI / Bohr_meter) / 100. ! Siemens / cm
2718 : ! fact_mob = fact_sigma * 100.**3 / e_Cb * Bohr_meter**3 ! cm^2 / V / s
2719 : ! fact_sbk = volt_SI * kb_HaK ! Volt / Kelvin
2720 : ! sigma_eh = fact_sigma * sigma_eh ! siemens cm^-1
2721 : ! mob_eh = fact_mob * mob_eh !
2722 : !call wrtout(std_out, "factors: fact_sigma, fact_mob", pre_newlines=1, newlines=1)
2723 : ! write (std_out,*) fact_sigma, fact_mob
2724 :
2725 31 : end subroutine ibte_calc_tensors
2726 : !!***
2727 :
2728 : ! Invert 3x3 matrix, copied from matr3inv
2729 52066 : pure subroutine inv33(aa, ait)
2730 :
2731 : !Arguments ------------------------------------
2732 : !arrays
2733 : real(dp),intent(in) :: aa(3,3)
2734 : real(dp),intent(out) :: ait(3,3)
2735 :
2736 : !Local variables-------------------------------
2737 : !scalars
2738 : real(dp) :: dd,det,t1,t2,t3
2739 :
2740 : ! *************************************************************************
2741 :
2742 52066 : t1 = aa(2,2) * aa(3,3) - aa(3,2) * aa(2,3)
2743 52066 : t2 = aa(3,2) * aa(1,3) - aa(1,2) * aa(3,3)
2744 52066 : t3 = aa(1,2) * aa(2,3) - aa(2,2) * aa(1,3)
2745 52066 : det = aa(1,1) * t1 + aa(2,1) * t2 + aa(3,1) * t3
2746 :
2747 : ! Make sure matrix is not singular
2748 52066 : if (abs(det) > 100 * tiny(one)) then
2749 15779 : dd = one / det
2750 15779 : ait(1,1) = t1 * dd
2751 15779 : ait(2,1) = t2 * dd
2752 15779 : ait(3,1) = t3 * dd
2753 15779 : ait(1,2) = (aa(3,1)*aa(2,3)-aa(2,1)*aa(3,3)) * dd
2754 15779 : ait(2,2) = (aa(1,1)*aa(3,3)-aa(3,1)*aa(1,3)) * dd
2755 15779 : ait(3,2) = (aa(2,1)*aa(1,3)-aa(1,1)*aa(2,3)) * dd
2756 15779 : ait(1,3) = (aa(2,1)*aa(3,2)-aa(3,1)*aa(2,2)) * dd
2757 15779 : ait(2,3) = (aa(3,1)*aa(1,2)-aa(1,1)*aa(3,2)) * dd
2758 15779 : ait(3,3) = (aa(1,1)*aa(2,2)-aa(2,1)*aa(1,2)) * dd
2759 410254 : ait = transpose(ait)
2760 : else
2761 36287 : ait = zero
2762 : end if
2763 :
2764 52066 : end subroutine inv33
2765 :
2766 51566 : end module m_rta
2767 : !!***
|
