Branch data 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 : : !!***
|
