Line data Source code
1 : !!****m* ABINIT/m_a2ftr
2 : !! NAME
3 : !! m_a2ftr
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2004-2026 ABINIT group (JPC, MJV, BXU)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_a2ftr
23 :
24 : use defs_basis
25 : use defs_elphon
26 : use m_errors
27 : use m_abicore
28 : use m_xmpi
29 : use m_splines
30 : use m_ebands
31 :
32 : use m_io_tools, only : open_file
33 : use m_numeric_tools, only : simpson_int
34 : use m_hide_lapack, only : matrginv
35 : use m_geometry, only : phdispl_cart2red
36 : use m_crystal, only : crystal_t
37 : use m_ifc, only : ifc_type
38 : use m_dynmat, only : ftgam_init, ftgam
39 : use m_epweights, only : d2c_wtq, ep_ph_weights, ep_el_weights, ep_ph_weights
40 :
41 : implicit none
42 :
43 : private
44 : !!***
45 :
46 : public :: mka2f_tr
47 : public :: mka2f_tr_lova
48 : public :: get_tau_k
49 : !!***
50 :
51 : contains
52 : !!***
53 :
54 : !!****f* ABINIT/mka2f_tr
55 : !!
56 : !! NAME
57 : !! mka2f_tr
58 : !!
59 : !! FUNCTION
60 : !! calculates the FS averaged Transport alpha^2F_tr function
61 : !! calculates and outputs the associated electrical conductivity, relaxation time, and Seebeck coefficient
62 : !! and thermal conductivities
63 : !! for the first task : copied from mka2F
64 : !!
65 : !! INPUTS
66 : !! crystal<crystal_t>=data type gathering info on the crystalline structure.
67 : !! Ifc<ifc_type>=Object containing the interatomic force constants.
68 : !! elph_ds
69 : !! elph_ds%gkk2 = gkk2 matrix elements on full FS grid for each phonon mode
70 : !! elph_ds%nbranch = number of phonon branches = 3*natom
71 : !! elph_ds%nFSband = number of bands included in the FS integration
72 : !! elph_ds%k_phon%nkpt = number of kpts included in the FS integration
73 : !! elph_ds%k_phon%wtk = integration weights on the FS
74 : !! delph_ds%n0 = DOS at the Fermi level calculated from the k_phon integration weights
75 : !! elph_ds%k_phon%kpt = coordinates of all FS kpoints
76 : !! mustar = coulomb pseudopotential parameter eventually for 2 spin channels
77 : !! natom = number of atoms
78 : !! ntemper = number of temperature points to calculate, from tempermin to tempermin+ntemper*temperinc
79 : !! tempermin = minimum temperature at which resistivity etc are calculated (in K)
80 : !! temperinc = interval for temperature grid on which resistivity etc are calculated (in K)
81 : !! elph_tr_ds%dos_n0 = DOS at the Fermi level calculated from the k_phon integration
82 : !! weights, but has a temperature dependence
83 : !! elph_tr_ds%dos_n = DOS at varied energy levels around Fermi level
84 : !! elph_tr_ds%veloc_sq0 = Fermi velocity square with T dependence
85 : !!
86 : !! OUTPUT
87 : !! elph_ds
88 : !!
89 : !! NOTES
90 : !! copied from ftiaf9.f
91 : !!
92 : !! SOURCE
93 :
94 0 : subroutine mka2f_tr(crystal,ifc,elph_ds,ntemper,tempermin,temperinc,pair2red,elph_tr_ds)
95 :
96 : !Arguments ------------------------------------
97 : !scalars
98 : integer,intent(in) :: ntemper
99 : real(dp),intent(in) :: tempermin,temperinc
100 : type(ifc_type),intent(in) :: ifc
101 : type(crystal_t),intent(in) :: crystal
102 : type(elph_tr_type),intent(inout) :: elph_tr_ds
103 : type(elph_type),intent(inout) :: elph_ds
104 : !arrays
105 : integer,intent(in) :: pair2red(elph_ds%nenergy,elph_ds%nenergy)
106 :
107 : !Local variables -------------------------
108 : !x =w/(2kbT)
109 : !scalars
110 : integer :: iFSqpt,ibranch,iomega,isppol,jbranch,nerr
111 : integer :: unit_a2f_tr,natom,ii,jj
112 : integer :: idir, iatom, k1, kdir
113 : integer :: unit_lor,unit_rho,unit_tau,unit_sbk, unit_therm
114 : integer :: itemp, tmp_nenergy
115 : integer :: itrtensor, icomp, jcomp!, kcomp
116 : integer :: ie, ie_1, ie2, ie_2, ie1, ie_tmp, ssp, s1(4), s2(4)
117 : integer :: ie2_left, ie2_right
118 : integer :: ik_this_proc, ierr,nrpt
119 : logical,parameter :: debug=.False.
120 : real(dp) :: Temp,chgu,chwu,diagerr,ucvol
121 : real(dp) :: a2fprefactor, gtemp
122 : real(dp) :: lambda_tr,lor0,lorentz,maxerr,omega
123 : real(dp) :: rho,tau,wtherm,xtr
124 : real(dp) :: lambda_tr_trace
125 : real(dp) :: domega, omega_min, omega_max
126 : real(dp) :: gaussval, gaussprefactor, gaussfactor, gaussmaxarg, xx
127 : real(dp) :: qnorm2, tmp_fermie
128 : real(dp) :: e1, e2, diff, xe
129 : real(dp) :: occ_omega, occ_e1, occ_e2
130 : real(dp) :: nv1, nv2, sigma1, sigma2
131 : real(dp) :: dos_n_e2, veloc_sq_icomp, veloc_sq_jcomp
132 : real(dp) :: tointegq00_1, tointegq00_2, tointegq01_1, tointegq01_2,tointegq11_1, tointegq11_2
133 : real(dp) :: j00, j01, j11
134 : real(dp) :: tointegq00,tointegq01,tointegq11
135 : real(dp) :: pref_s, pref_w, tmp_veloc_sq0, tmp_veloc_sq1, tmp_veloc_sq2
136 : character(len=500) :: message
137 : character(len=fnlen) :: fname
138 : !arrays
139 : complex(dp),parameter :: c0=dcmplx(0.d0,0.d0),c1=dcmplx(1.d0,0.d0)
140 : real(dp) :: gprimd(3,3)
141 0 : real(dp) :: eigval(elph_ds%nbranch)
142 0 : real(dp) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
143 0 : real(dp) :: gam_now(2,elph_ds%nbranch*elph_ds%nbranch)
144 0 : real(dp) :: tmpa2f(elph_ds%na2f)
145 0 : real(dp) :: tmpgam1(2,elph_ds%nbranch,elph_ds%nbranch)
146 0 : real(dp) :: tmpgam2(2,elph_ds%nbranch,elph_ds%nbranch)
147 : real(dp) :: q11_inv(3,3)
148 : !real(dp) :: fullq(6,6)
149 0 : real(dp),allocatable :: phfrq(:,:)
150 0 : real(dp),allocatable :: tmp_a2f_1d_tr(:,:,:,:,:)
151 0 : real(dp),allocatable :: displ(:,:,:,:)
152 0 : real(dp),allocatable :: pheigvec(:,:)
153 0 : real(dp),allocatable :: tmp_wtq(:,:,:)
154 0 : real(dp),allocatable :: integrho(:), tointegrho(:)
155 0 : real(dp),allocatable :: integrand_q00(:),integrand_q01(:),integrand_q11(:)
156 0 : real(dp),allocatable :: q00(:,:,:,:), q01(:,:,:,:),q11(:,:,:,:)
157 0 : real(dp),allocatable :: seebeck(:,:,:,:)!, rho_nm(:,:,:,:)
158 0 : real(dp),allocatable :: rho_T(:)
159 0 : real(dp),allocatable :: coskr(:,:), sinkr(:,:), coskr_tmp(:), sinkr_tmp(:)
160 0 : real(dp),allocatable :: gam_rpt(:,:,:)
161 :
162 : ! *********************************************************************
163 : !calculate a2f_tr for frequencies between 0 and omega_max
164 :
165 :
166 0 : write(std_out,*) 'mka2f_tr : enter '
167 :
168 0 : ucvol = crystal%ucvol
169 0 : natom = crystal%natom
170 0 : gprimd = crystal%gprimd
171 :
172 : ! number of real-space points for FT interpolation
173 0 : nrpt = ifc%nrpt
174 : !
175 : !MG: the step should be calculated locally using nomega and the extrema of the spectrum.
176 : !One should not rely on previous calls for the setup of elph_ds%domega
177 : !I will remove elph_ds%domega since mka2f.F90 will become a method of gamma_t
178 0 : domega =elph_ds%domega
179 0 : omega_min = elph_ds%omega_min
180 0 : omega_max = elph_ds%omega_max
181 :
182 0 : if (elph_ds%ep_lova .eq. 1) then
183 : tmp_nenergy = 1
184 0 : else if (elph_ds%ep_lova .eq. 0) then
185 0 : tmp_nenergy = elph_ds%nenergy
186 : end if
187 :
188 : !! defaults for number of temperature steps and max T (all in Kelvin...)
189 : !ntemper=1000
190 : !tempermin=zero
191 : !temperinc=one
192 :
193 0 : ABI_MALLOC(rho_T,(ntemper))
194 :
195 :
196 0 : gaussprefactor = sqrt(piinv) / elph_ds%a2fsmear
197 0 : gaussfactor = one / elph_ds%a2fsmear
198 0 : gaussmaxarg = sqrt(-log(1.d-90))
199 : !lor0=(pi*kb_HaK)**2/3.
200 0 : lor0=pi**2/3.0_dp
201 :
202 : !maximum value of frequency (a grid has to be chosen for the representation of alpha^2 F)
203 : !WARNING! supposes this value has been set in mkelph_linwid.
204 :
205 : !ENDMG
206 :
207 0 : maxerr=0.
208 0 : nerr=0
209 :
210 0 : ABI_MALLOC(tmp_wtq,(elph_ds%nbranch, elph_ds%k_fine%nkpt, elph_ds%na2f+1))
211 0 : ABI_MALLOC(elph_ds%k_fine%wtq,(elph_ds%nbranch, elph_ds%k_fine%nkpt, elph_ds%na2f))
212 0 : ABI_MALLOC(elph_ds%k_phon%wtq,(elph_ds%nbranch, elph_ds%k_phon%nkpt, elph_ds%na2f))
213 :
214 0 : ABI_MALLOC(phfrq,(elph_ds%nbranch, elph_ds%k_fine%nkpt))
215 0 : ABI_MALLOC(displ,(2, elph_ds%nbranch, elph_ds%nbranch, elph_ds%k_fine%nkpt))
216 0 : ABI_MALLOC(pheigvec,(2*elph_ds%nbranch*elph_ds%nbranch, elph_ds%k_fine%nkpt))
217 :
218 0 : do iFSqpt=1,elph_ds%k_fine%nkpt
219 0 : call ifc%fourq(crystal,elph_ds%k_fine%kpt(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt),out_eigvec=pheigvec(:,iFSqpt))
220 : end do
221 0 : omega_min = omega_min - domega
222 :
223 : !bxu, obtain wtq for the q_fine, then condense to q_phon
224 : call ep_ph_weights(phfrq,elph_ds%a2fsmear,omega_min,omega_max,elph_ds%na2f+1,gprimd,elph_ds%kptrlatt_fine, &
225 0 : & elph_ds%nbranch,elph_ds%telphint,elph_ds%k_fine,tmp_wtq)
226 : !call ep_ph_weights(phfrq,elph_ds%a2fsmear,omega_min,omega_max,elph_ds%na2f+1,gprimd,elph_ds%kptrlatt_fine, &
227 : !& elph_ds%nbranch,1,elph_ds%k_fine,tmp_wtq)
228 0 : omega_min = omega_min + domega
229 :
230 0 : do iomega = 1, elph_ds%na2f
231 0 : elph_ds%k_fine%wtq(:,:,iomega) = tmp_wtq(:,:,iomega+1)
232 : end do
233 0 : ABI_FREE(tmp_wtq)
234 :
235 0 : if (elph_ds%use_k_fine == 1) then
236 0 : call d2c_wtq(elph_ds)
237 : end if
238 :
239 0 : ABI_FREE(phfrq)
240 0 : ABI_FREE(displ)
241 0 : ABI_FREE(pheigvec)
242 :
243 : !reduce the dimension from fine to phon for phfrq and pheigvec
244 0 : ABI_MALLOC(phfrq,(elph_ds%nbranch, elph_ds%k_phon%nkpt))
245 0 : ABI_MALLOC(displ,(2, elph_ds%nbranch, elph_ds%nbranch, elph_ds%k_phon%nkpt))
246 0 : ABI_MALLOC(pheigvec,(2*elph_ds%nbranch*elph_ds%nbranch, elph_ds%k_phon%nkpt))
247 :
248 0 : do iFSqpt=1,elph_ds%k_phon%nkpt
249 0 : call ifc%fourq(crystal,elph_ds%k_phon%kpt(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt),out_eigvec=pheigvec(:,iFSqpt))
250 : end do
251 :
252 0 : ABI_MALLOC(elph_tr_ds%a2f_1d_tr,(elph_ds%na2f,9,elph_ds%nsppol,4,tmp_nenergy**2,ntemper))
253 0 : ABI_MALLOC(tmp_a2f_1d_tr,(elph_ds%na2f,9,elph_ds%nsppol,4,tmp_nenergy**2))
254 :
255 : ! prepare phase factors
256 0 : ABI_MALLOC(coskr, (elph_ds%k_phon%nkpt, nrpt))
257 0 : ABI_MALLOC(sinkr, (elph_ds%k_phon%nkpt, nrpt))
258 0 : call ftgam_init(ifc%gprim, elph_ds%k_phon%nkpt, nrpt, elph_ds%k_phon%kpt, ifc%rpt, coskr, sinkr)
259 :
260 0 : elph_tr_ds%a2f_1d_tr = zero
261 0 : tmp_a2f_1d_tr = zero
262 :
263 0 : ABI_MALLOC(gam_rpt,(2,3*natom*3*natom,nrpt))
264 0 : ABI_MALLOC(coskr_tmp,(nrpt))
265 0 : ABI_MALLOC(sinkr_tmp,(nrpt))
266 :
267 0 : do ie = 1, elph_ds%n_pair
268 0 : do ssp = 1,4
269 0 : do isppol = 1, elph_ds%nsppol
270 :
271 : ! loop over qpoint in full kpt grid (presumably dense)
272 0 : do ik_this_proc =1,elph_ds%k_phon%my_nkpt
273 0 : iFSqpt = elph_ds%k_phon%my_ikpt(ik_this_proc)
274 :
275 0 : qnorm2 = sum(elph_ds%k_phon%kpt(:,iFSqpt)**2)
276 : ! if (flag_to_exclude_soft_modes = .false.) qnorm2 = zero
277 0 : do itrtensor=1,9
278 :
279 0 : if (elph_ds%ep_int_gkk == 1) then
280 0 : gam_now(:,:) = elph_tr_ds%gamma_qpt_tr(:,itrtensor,:,isppol,iFSqpt)
281 : else
282 : ! Do FT from real-space gamma grid to 1 qpt.
283 0 : gam_rpt(:,:,:)=elph_tr_ds%gamma_rpt_tr(:,itrtensor,:,isppol,:,ssp,ie)
284 0 : coskr_tmp(:)=coskr(iFSqpt,:)
285 0 : sinkr_tmp(:)=sinkr(iFSqpt,:)
286 0 : call ftgam(ifc%wghatm,gam_now,gam_rpt,natom,1,nrpt,0,coskr_tmp, sinkr_tmp)
287 : end if
288 :
289 : ! Diagonalize gamma matrix at this qpoint (complex matrix).
290 :
291 : ! if ep_scalprod==0 we have to dot in the displacement vectors here
292 0 : if (elph_ds%ep_scalprod==0) then
293 :
294 0 : displ_red(:,:,:) = zero
295 0 : do jbranch=1,elph_ds%nbranch
296 0 : do iatom=1,natom
297 0 : do idir=1,3
298 0 : ibranch=idir+3*(iatom-1)
299 0 : do kdir=1,3
300 0 : k1 = kdir+3*(iatom-1)
301 : displ_red(1,ibranch,jbranch) = displ_red(1,ibranch,jbranch) + &
302 0 : & gprimd(kdir,idir)*displ(1,k1,jbranch,iFSqpt)
303 : displ_red(2,ibranch,jbranch) = displ_red(2,ibranch,jbranch) + &
304 0 : & gprimd(kdir,idir)*displ(2,k1,jbranch,iFSqpt)
305 : end do
306 : end do
307 : end do
308 : end do
309 :
310 0 : tmpgam2 = reshape (gam_now, (/2,elph_ds%nbranch,elph_ds%nbranch/))
311 0 : call gam_mult_displ(elph_ds%nbranch, displ_red, tmpgam2, tmpgam1)
312 0 : do jbranch=1,elph_ds%nbranch
313 0 : eigval(jbranch) = tmpgam1(1, jbranch, jbranch)
314 : end do
315 :
316 0 : else if (elph_ds%ep_scalprod == 1) then
317 :
318 :
319 : ! NOTE: in these calls gam_now and pheigvec do not have the right rank, but blas usually does not care
320 :
321 : call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, c1, gam_now, 3*natom,&
322 0 : & pheigvec(:,iFSqpt), 3*natom, c0, tmpgam1, 3*natom)
323 : call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, c1, pheigvec(:,iFSqpt), 3*natom,&
324 0 : & tmpgam1, 3*natom, c0, tmpgam2, 3*natom)
325 0 : diagerr = zero
326 0 : do ibranch=1,elph_ds%nbranch
327 0 : eigval(ibranch) = tmpgam2(1,ibranch,ibranch)
328 0 : do jbranch=1,ibranch-1
329 0 : diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch))
330 : end do
331 0 : do jbranch=ibranch+1,elph_ds%nbranch
332 0 : diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch))
333 : end do
334 : end do
335 0 : if (diagerr > tol12) then
336 0 : nerr=nerr+1
337 0 : maxerr=max(diagerr, maxerr)
338 : end if
339 : end if ! end ep_scalprod if
340 :
341 : ! Add all contributions from the phonon modes at this qpoint to a2f and the phonon dos.
342 0 : do ibranch=1,elph_ds%nbranch
343 : ! if (abs(phfrq(ibranch,iFSqpt)) < tol10) then
344 0 : if ( abs(phfrq(ibranch,iFSqpt)) < tol7 .or. &
345 : & (phfrq(ibranch,iFSqpt) < tol4 .and. qnorm2 > 0.03 )) then
346 : ! note: this should depend on the velocity of sound, to accept acoustic modes!
347 : a2fprefactor = zero
348 : else
349 : ! a2fprefactor = eigval (ibranch)/(two_pi*abs(phfrq(ibranch,iFSqpt))*elph_ds%n0(isppol))
350 : ! Use the dos_n0 at the lowest input temperature, assuming to be low
351 0 : a2fprefactor = eigval (ibranch)/(two_pi*abs(phfrq(ibranch,iFSqpt)))
352 : end if
353 :
354 0 : omega = omega_min
355 0 : tmpa2f(:) = zero
356 0 : do iomega=1,elph_ds%na2f
357 0 : xx = (omega-phfrq(ibranch,iFSqpt))*gaussfactor
358 0 : omega = omega+domega
359 0 : if (abs(xx) > gaussmaxarg) cycle
360 :
361 0 : gaussval = gaussprefactor*exp(-xx*xx)
362 0 : gtemp = gaussval*a2fprefactor
363 :
364 0 : if (dabs(gtemp) < 1.0d-50) gtemp = zero
365 0 : tmpa2f(iomega) = tmpa2f(iomega) + gtemp
366 : end do
367 :
368 : ! tmpa2f(:) = zero
369 : ! do iomega=1,elph_ds%na2f
370 : ! gtemp = a2fprefactor*elph_ds%k_phon%wtq(ibranch,iFSqpt,iomega)
371 : ! if (dabs(gtemp) < 1.0d-50) gtemp = zero
372 : ! tmpa2f(iomega) = tmpa2f(iomega) + gtemp
373 : ! end do
374 :
375 0 : tmp_a2f_1d_tr (:,itrtensor,isppol,ssp,ie) = tmp_a2f_1d_tr (:,itrtensor,isppol,ssp,ie) + tmpa2f(:)
376 :
377 : end do ! end ibranch
378 : end do ! end itrtensor
379 : end do ! end iFSqpt - loop done in parallel
380 : end do ! end isppol
381 : end do ! ss'
382 : end do ! n_pair
383 :
384 : ! MG: FIXME: Why xmpi_world? besides only one CPU should perform IO (see below)
385 : ! Likely this routine is never executed in parallel
386 0 : call xmpi_sum (tmp_a2f_1d_tr, xmpi_world, ierr)
387 :
388 0 : ABI_FREE(gam_rpt)
389 0 : ABI_FREE(coskr_tmp)
390 0 : ABI_FREE(sinkr_tmp)
391 0 : ABI_FREE(coskr)
392 0 : ABI_FREE(sinkr)
393 :
394 0 : do itemp=1,ntemper ! runs over termperature in K
395 0 : do isppol=1,elph_ds%nsppol
396 0 : do jj=1,tmp_nenergy**2
397 0 : do ii=1,4
398 0 : elph_tr_ds%a2f_1d_tr(:,:,isppol,ii,jj,itemp) = tmp_a2f_1d_tr(:,:,isppol,ii,jj)/elph_tr_ds%dos_n0(itemp,isppol)
399 : end do
400 : end do
401 : end do
402 : end do
403 :
404 0 : ABI_FREE(tmp_a2f_1d_tr)
405 :
406 : !second 1 / elph_ds%k_phon%nkpt factor for the integration weights
407 0 : elph_tr_ds%a2f_1d_tr = elph_tr_ds%a2f_1d_tr / elph_ds%k_phon%nkpt
408 :
409 0 : if (elph_ds%ep_scalprod == 1) then
410 0 : write(std_out,*) 'mka2f_tr: errors in diagonalization of gamma_tr with phon eigenvectors: ', nerr,maxerr
411 : end if
412 :
413 : !output the elph_tr_ds%a2f_1d_tr
414 0 : fname = trim(elph_ds%elph_base_name) // '_A2F_TR'
415 0 : if (open_file(fname,message,newunit=unit_a2f_tr,status='unknown') /= 0) then
416 0 : ABI_ERROR(message)
417 : end if
418 :
419 0 : write (unit_a2f_tr,'(a)') '#'
420 0 : write (unit_a2f_tr,'(a)') '# ABINIT package : a2f_tr file'
421 0 : write (unit_a2f_tr,'(a)') '#'
422 0 : write (unit_a2f_tr,'(a)') '# a2f_tr function integrated over the FS. omega in a.u.'
423 0 : write (unit_a2f_tr,'(a,I10)') '# number of kpoints integrated over : ', elph_ds%k_phon%nkpt
424 0 : write (unit_a2f_tr,'(a,I10)') '# number of energy points : ',elph_ds%na2f
425 0 : write (unit_a2f_tr,'(a,E16.6,a,E16.6,a)') '# between omega_min = ', omega_min,' Ha and omega_max = ', omega_max, ' Ha'
426 0 : write (unit_a2f_tr,'(a,E16.6)') '# and the smearing width for gaussians is ', elph_ds%a2fsmear
427 0 : write (unit_a2f_tr,'(a)') '#'
428 :
429 : !done with header
430 0 : do isppol=1,elph_ds%nsppol
431 0 : write (unit_a2f_tr,'(a,E16.6)') '# The DOS at Fermi level is ', elph_tr_ds%dos_n0(1,isppol)
432 0 : omega = omega_min
433 0 : do iomega=1,elph_ds%na2f
434 : ! bxu, at which eps and eps' should I save it
435 : ! better to save them all, but could be too many
436 0 : write (unit_a2f_tr, '(10D16.6)') omega, elph_tr_ds%a2f_1d_tr(iomega,:,isppol,1,INT(elph_ds%n_pair/2)+1,1)
437 0 : omega=omega+domega
438 : end do
439 0 : write (unit_a2f_tr,*)
440 : end do !isppol
441 :
442 0 : close (unit=unit_a2f_tr)
443 :
444 : !calculation of transport properties
445 0 : ABI_MALLOC(integrho,(elph_ds%na2f))
446 0 : ABI_MALLOC(tointegrho,(elph_ds%na2f))
447 0 : ABI_MALLOC(integrand_q00,(elph_ds%na2f))
448 0 : ABI_MALLOC(integrand_q01,(elph_ds%na2f))
449 0 : ABI_MALLOC(integrand_q11,(elph_ds%na2f))
450 0 : ABI_MALLOC(q00,(ntemper,3,3,elph_ds%nsppol))
451 0 : ABI_MALLOC(q01,(ntemper,3,3,elph_ds%nsppol))
452 0 : ABI_MALLOC(q11,(ntemper,3,3,elph_ds%nsppol))
453 0 : ABI_MALLOC(seebeck,(elph_ds%nsppol,ntemper,3,3))
454 : !ABI_MALLOC(rho_nm,(elph_ds%nsppol,ntemper,3,3))
455 :
456 0 : fname = trim(elph_ds%elph_base_name) // '_RHO'
457 0 : if (open_file(fname,message,newunit=unit_rho,status='unknown') /= 0) then
458 0 : ABI_ERROR(message)
459 : end if
460 : !print header to resistivity file
461 0 : write (unit_rho,*) '# Resistivity as a function of temperature.'
462 0 : write (unit_rho,*) '# the formalism is isotropic, so non-cubic crystals may be wrong'
463 0 : write (unit_rho,*) '# '
464 0 : write (unit_rho,*) '# Columns are: '
465 0 : write (unit_rho,*) '# temperature[K] rho[au] rho [SI] rho/temp [au]'
466 0 : write (unit_rho,*) '# '
467 :
468 0 : fname = trim(elph_ds%elph_base_name) // '_TAU'
469 0 : if (open_file(fname,message,newunit=unit_tau,status='unknown') /= 0) then
470 0 : ABI_ERROR(message)
471 : end if
472 : !print header to relaxation time file
473 0 : write (unit_tau,*) '# Relaxation time as a function of temperature.'
474 0 : write (unit_tau,*) '# the formalism is isotropic, so non-cubic crystals may be wrong'
475 0 : write (unit_tau,*) '# '
476 0 : write (unit_tau,*) '# Columns are: '
477 0 : write (unit_tau,*) '# temperature[K] tau[au] tau [SI] '
478 0 : write (unit_tau,*) '# '
479 :
480 0 : fname = trim(elph_ds%elph_base_name) // '_SBK'
481 0 : if (open_file(fname,message,newunit=unit_sbk,status='unknown') /= 0) then
482 0 : ABI_ERROR(message)
483 : end if
484 : !print header to relaxation time file
485 0 : write (unit_sbk,*) '# Seebeck Coefficint as a function of temperature.'
486 0 : write (unit_sbk,*) '# the formalism is isotropic, so non-cubic crystals may be wrong'
487 0 : write (unit_sbk,*) '# '
488 0 : write (unit_sbk,*) '# Columns are: '
489 0 : write (unit_sbk,*) '# temperature[K] S [au] S [SI] '
490 0 : write (unit_sbk,*) '# '
491 :
492 0 : fname = trim(elph_ds%elph_base_name) // '_WTH'
493 0 : if (open_file(fname,message,newunit=unit_therm,status='unknown') /= 0) then
494 0 : ABI_ERROR(message)
495 : end if
496 :
497 : !print header to thermal conductivity file
498 0 : write (unit_therm,'(a)') '# Thermal conductivity/resistivity as a function of temperature.'
499 0 : write (unit_therm,'(a)') '# the formalism is isotropic, so non-cubic crystals may be wrong'
500 0 : write (unit_therm,'(a)') '# '
501 0 : write (unit_therm,'(a)') '# Columns are: '
502 0 : write (unit_therm,'(a)') '# temperature[K] thermal rho[au] thermal cond [au] thermal rho [SI] thermal cond [SI]'
503 0 : write (unit_therm,'(a)') '# '
504 :
505 0 : fname = trim(elph_ds%elph_base_name) // '_LOR'
506 0 : if (open_file(fname,message,newunit=unit_lor,status='unknown') /= 0) then
507 0 : ABI_ERROR(message)
508 : end if
509 :
510 : !print header to lorentz file
511 0 : write (unit_lor,*) '# Lorentz number as a function of temperature.'
512 0 : write (unit_lor,*) '# the formalism is isotropic, so non-cubic crystals may be wrong'
513 0 : write (unit_lor,*) '# '
514 0 : write (unit_lor,*) '# Columns are: '
515 0 : write (unit_lor,*) '# temperature[K] Lorentz number[au] Lorentz quantum = (pi*kb_HaK)**2/3'
516 0 : write (unit_lor,*) '# '
517 :
518 0 : do isppol=1,elph_ds%nsppol
519 0 : lambda_tr_trace = zero
520 0 : do itrtensor=1,9
521 0 : omega = omega_min
522 0 : tointegrho = zero
523 0 : do iomega=1,elph_ds%na2f
524 0 : if(omega<=0) then
525 0 : omega=omega+domega
526 0 : cycle
527 : end if
528 : ! bxu, agian, which eps and eps' to use?
529 : ! sometimes Ef is in the gap
530 0 : tointegrho(iomega)=two*elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1)/omega
531 0 : omega=omega+domega
532 : end do
533 :
534 0 : integrho = zero
535 0 : call simpson_int(elph_ds%na2f,domega,tointegrho,integrho)
536 0 : lambda_tr=integrho(elph_ds%na2f)
537 : write (message, '(a,2i3,a,es16.6)' )&
538 0 : & ' mka2f_tr: TRANSPORT lambda for isppol itrtensor', isppol, itrtensor, ' = ', lambda_tr
539 0 : call wrtout(std_out,message,'COLL')
540 0 : if (itrtensor == 1 .or. itrtensor == 5 .or. itrtensor == 9) lambda_tr_trace = lambda_tr_trace + lambda_tr
541 : end do !end itrtensor do
542 :
543 0 : lambda_tr_trace = lambda_tr_trace / three
544 : write (message, '(a,i3,a,es16.6)' )&
545 0 : & ' mka2f_tr: 1/3 trace of TRANSPORT lambda for isppol ', isppol, ' = ', lambda_tr_trace
546 0 : call wrtout(std_out,message,'COLL')
547 0 : call wrtout(ab_out,message,'COLL')
548 : end do !end isppol do
549 :
550 : !constant to change units of rho from au to SI
551 0 : chgu=2.173969*1.0d-7 ! REPLACE WITH DATA FROM DEFS_BASIS
552 0 : chwu=9.270955772*1.0d-5 ! au to mK/W ! REPLACE WITH DATA FROM DEFS_BASIS
553 :
554 : !change the fermi level to zero, as required for q01 to vanish.
555 0 : tmp_fermie = elph_ds%fermie
556 : !Get Q00, Q01, Q11, and derive rho, tau
557 0 : q00 = zero
558 0 : q01 = zero
559 0 : q11 = zero
560 : ! prepare s1 and s2 arrays
561 0 : s1 = (/1, 1, -1, -1/)
562 0 : s2 = (/1, -1, 1, -1/)
563 :
564 0 : do isppol=1,elph_ds%nsppol
565 0 : do icomp=1, 3
566 0 : do jcomp=1, 3
567 0 : itrtensor=(icomp-1)*3+jcomp
568 :
569 0 : write(unit_rho,*) '# Rho for isppol, itrten = ', isppol, itrtensor
570 0 : write(unit_tau,*) '# Tau for isppol, itrten = ', isppol, itrtensor
571 :
572 0 : do itemp=1,ntemper ! runs over termperature in K
573 0 : Temp=tempermin+temperinc*dble(itemp)
574 0 : tmp_veloc_sq0 = sqrt(elph_tr_ds%veloc_sq0(itemp,icomp,isppol)*elph_tr_ds%veloc_sq0(itemp,jcomp,isppol))
575 :
576 0 : integrand_q00 = zero
577 0 : integrand_q01 = zero
578 0 : integrand_q11 = zero
579 :
580 0 : omega = omega_min
581 0 : do iomega=1,elph_ds%na2f
582 0 : if(omega .le. 0) then
583 0 : omega=omega+domega
584 0 : cycle
585 : end if
586 0 : xtr=omega/(kb_HaK*Temp)
587 0 : occ_omega=1.0_dp/(exp(xtr)-1.0_dp)
588 :
589 0 : tmp_veloc_sq1 = zero
590 0 : tmp_veloc_sq2 = zero
591 0 : do ie1=1,elph_ds%nenergy
592 0 : e1 = elph_tr_ds%en_all(isppol,ie1)
593 :
594 : !! BXU, the tolerance here needs to be used with caution
595 : !! which depends on the dimensions of the system
596 : !! E.g. in 2D system, DOS can be much smaller
597 0 : if (elph_tr_ds%dos_n(ie1,isppol)/natom .lt. 0.1d0) cycle ! energy in the gap forbidden
598 :
599 0 : xtr=(e1-tmp_fermie)/(kb_HaK*Temp)
600 0 : occ_e1=1.0_dp/(exp(xtr)+1.0_dp)
601 :
602 0 : e2 = e1 + omega
603 0 : xtr=(e2-tmp_fermie)/(kb_HaK*Temp)
604 0 : occ_e2=1.0_dp/(exp(xtr)+1.0_dp)
605 : ! Do we need to change the fermie to the one with T dependence?
606 : ! find which ie2 give the closest energy
607 0 : if (e2 .gt. elph_tr_ds%en_all(isppol,elph_ds%nenergy)) then
608 0 : ie2 = 0
609 : cycle
610 : else
611 0 : ie_tmp = 1
612 0 : diff = dabs(e2-elph_tr_ds%en_all(isppol,1))
613 0 : do ie2 = 2, elph_ds%nenergy
614 0 : if (dabs(e2-elph_tr_ds%en_all(isppol,ie2)) .lt. diff) then
615 0 : diff = dabs(e2-elph_tr_ds%en_all(isppol,ie2))
616 0 : ie_tmp = ie2
617 : end if
618 : end do
619 0 : ie2 = ie_tmp
620 :
621 0 : if (e2 < elph_tr_ds%en_all(isppol,ie2)) then
622 0 : ie2_right = ie2
623 0 : ie2_left = ie2-1
624 : else
625 0 : ie2_right = ie2+1
626 0 : ie2_left = ie2
627 : end if
628 :
629 : end if
630 :
631 0 : if (elph_tr_ds%dos_n(ie2,isppol)/natom .lt. 0.1d0) cycle
632 :
633 0 : tointegq00 = zero
634 0 : tointegq01 = zero
635 0 : tointegq11 = zero
636 :
637 : ! BXU linear interpolation of volec_sq and dos_n
638 : xe=(e2-elph_tr_ds%en_all(isppol,ie2_left))/ &
639 0 : & (elph_tr_ds%en_all(isppol,ie2_right)-elph_tr_ds%en_all(isppol,ie2_left))
640 : veloc_sq_icomp = elph_tr_ds%veloc_sq(icomp,isppol,ie2_left)*(1.0d0-xe) + &
641 0 : & elph_tr_ds%veloc_sq(icomp,isppol,ie2_right)*xe
642 : veloc_sq_jcomp = elph_tr_ds%veloc_sq(jcomp,isppol,ie2_left)*(1.0d0-xe) + &
643 0 : & elph_tr_ds%veloc_sq(jcomp,isppol,ie2_right)*xe
644 : dos_n_e2 = elph_tr_ds%dos_n(ie2_left,isppol)*(1.0d0-xe) + &
645 0 : & elph_tr_ds%dos_n(ie2_right,isppol)*xe
646 :
647 0 : tmp_veloc_sq1 = sqrt(elph_tr_ds%veloc_sq(icomp,isppol,ie1)*elph_tr_ds%veloc_sq(jcomp,isppol,ie1))
648 : ! tmp_veloc_sq2 = sqrt(elph_tr_ds%veloc_sq(icomp,isppol,ie2)*elph_tr_ds%veloc_sq(jcomp,isppol,ie2))
649 0 : tmp_veloc_sq2 = sqrt(veloc_sq_icomp*veloc_sq_jcomp)
650 :
651 : ! ie_1 = (ie1-1)*elph_ds%nenergy + ie2
652 : ! ie_2 = (ie2-1)*elph_ds%nenergy + ie1
653 0 : ie_1 = pair2red(ie1,ie2)
654 0 : ie_2 = pair2red(ie2,ie1)
655 0 : if (ie_1 .eq. 0 .or. ie_2 .eq. 0) then
656 0 : ABI_BUG('CHECK pair2red!')
657 : end if
658 :
659 0 : do ssp=1,4 ! (s,s'=+/-1, condense the indices)
660 :
661 0 : nv1 = 1.0_dp/(elph_tr_ds%dos_n(ie1,isppol)*sqrt(tmp_veloc_sq1))
662 0 : sigma1 = sqrt(3.0_dp)*(e1-tmp_fermie)/(pi*Temp*kb_HaK)
663 : !DEBUG
664 0 : if (elph_ds%ep_lova .eq. 1) then
665 0 : nv1 = 1.0_dp/(elph_tr_ds%dos_n0(itemp,isppol)*sqrt(tmp_veloc_sq0))
666 0 : sigma1 = sqrt(3.0_dp)*(e1-tmp_fermie)/(pi*Temp*kb_HaK)
667 : end if
668 : !ENDDEBUG
669 :
670 0 : tointegq00_1 = zero
671 0 : tointegq01_1 = zero
672 0 : tointegq11_1 = zero
673 :
674 : !DEBUG
675 0 : if (elph_ds%ep_lova .eq. 1) then
676 0 : nv2 = 1.0_dp/(elph_tr_ds%dos_n0(itemp,isppol)*sqrt(tmp_veloc_sq0))
677 0 : sigma2 = sqrt(3.0_dp)*(e2-tmp_fermie)/(pi*Temp*kb_HaK)
678 0 : j00 = (nv1 + s1(ssp)*nv2)*(nv1 + s2(ssp)*nv2)/4.0_dp
679 0 : j01 = (nv1 + s1(ssp)*nv2)*(nv1*sigma1 + s2(ssp)*nv2*sigma2)/4.0_dp
680 0 : j11 = (nv1*sigma1 + s1(ssp)*nv2*sigma2)*(nv1*sigma1 + s2(ssp)*nv2*sigma2)/4.0_dp
681 : tointegq00_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* &
682 0 : & occ_e1*(1.0_dp-occ_e2)*j00*occ_omega
683 : tointegq01_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* &
684 0 : & occ_e1*(1.0_dp-occ_e2)*j01*occ_omega
685 : tointegq11_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* &
686 0 : & occ_e1*(1.0_dp-occ_e2)*j11*occ_omega
687 : !END DEBUG
688 0 : else if (elph_ds%ep_lova .eq. 0) then
689 0 : nv2 = 1.0_dp/(dos_n_e2*sqrt(tmp_veloc_sq2))
690 0 : sigma2 = sqrt(3.0_dp)*(e2-tmp_fermie)/(pi*Temp*kb_HaK)
691 0 : j00 = (nv1 + s1(ssp)*nv2)*(nv1 + s2(ssp)*nv2)/4.0_dp
692 0 : j01 = (nv1 + s1(ssp)*nv2)*(nv1*sigma1 + s2(ssp)*nv2*sigma2)/4.0_dp
693 0 : j11 = (nv1*sigma1 + s1(ssp)*nv2*sigma2)*(nv1*sigma1 + s2(ssp)*nv2*sigma2)/4.0_dp
694 : ! bxu TEST
695 : if (debug) then
696 : if (ssp .eq. 1 .and. itrtensor .eq. 1) then
697 : write(21,"(3i5,4E20.12)") iomega, ie1, ie2, &
698 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp), j01, &
699 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01, &
700 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01*occ_e1*(1.0_dp-occ_e2)*occ_omega
701 : end if
702 : if (ssp .eq. 2 .and. itrtensor .eq. 1) then
703 : write(22,"(3i5,4E20.12)") iomega, ie1, ie2, &
704 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp), j01, &
705 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01, &
706 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01*occ_e1*(1.0_dp-occ_e2)*occ_omega
707 : end if
708 : if (ssp .eq. 3 .and. itrtensor .eq. 1) then
709 : write(23,"(3i5,4E20.12)") iomega, ie1, ie2, &
710 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp), j01, &
711 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01, &
712 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01*occ_e1*(1.0_dp-occ_e2)*occ_omega
713 : end if
714 : if (ssp .eq. 4 .and. itrtensor .eq. 1) then
715 : write(24,"(3i5,4E20.12)") iomega, ie1, ie2, &
716 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp), j01, &
717 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01, &
718 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01*occ_e1*(1.0_dp-occ_e2)*occ_omega
719 : end if
720 : end if
721 : tointegq00_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)* &
722 0 : & occ_e1*(1.0_dp-occ_e2)*j00*occ_omega
723 : tointegq01_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)* &
724 0 : & occ_e1*(1.0_dp-occ_e2)*j01*occ_omega
725 : tointegq11_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)* &
726 0 : & occ_e1*(1.0_dp-occ_e2)*j11*occ_omega
727 : end if
728 :
729 0 : tointegq00_2 = zero
730 0 : tointegq01_2 = zero
731 0 : tointegq11_2 = zero
732 :
733 : !DEBUG
734 0 : if (elph_ds%ep_lova .eq. 1) then
735 0 : nv2 = 1.0_dp/(elph_tr_ds%dos_n0(itemp,isppol)*sqrt(tmp_veloc_sq0))
736 0 : sigma2 = sqrt(3.0_dp)*(e2-tmp_fermie)/(pi*Temp*kb_HaK)
737 0 : j00 = (nv2 + s1(ssp)*nv1)*(nv2 + s2(ssp)*nv1)/4.0_dp
738 0 : j01 = (nv2 + s1(ssp)*nv1)*(nv2*sigma2 + s2(ssp)*nv1*sigma1)/4.0_dp
739 0 : j11 = (nv2*sigma2 + s1(ssp)*nv1*sigma1)*(nv2*sigma2 + s2(ssp)*nv1*sigma1)/4.0_dp
740 : tointegq00_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* &
741 0 : & occ_e1*(1.0_dp-occ_e2)*j00*(occ_omega+1)
742 : tointegq01_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* &
743 0 : & occ_e1*(1.0_dp-occ_e2)*j01*(occ_omega+1)
744 : tointegq11_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* &
745 0 : & occ_e1*(1.0_dp-occ_e2)*j11*(occ_omega+1)
746 : !END DEBUG
747 0 : else if (elph_ds%ep_lova .eq. 0) then
748 0 : nv2 = 1.0_dp/(dos_n_e2*sqrt(tmp_veloc_sq2))
749 0 : sigma2 = sqrt(3.0_dp)*(e2-tmp_fermie)/(pi*Temp*kb_HaK)
750 0 : j00 = (nv2 + s1(ssp)*nv1)*(nv2 + s2(ssp)*nv1)/4.0_dp
751 0 : j01 = (nv2 + s1(ssp)*nv1)*(nv2*sigma2 + s2(ssp)*nv1*sigma1)/4.0_dp
752 0 : j11 = (nv2*sigma2 + s1(ssp)*nv1*sigma1)*(nv2*sigma2 + s2(ssp)*nv1*sigma1)/4.0_dp
753 : !DEBUG bxu TEST
754 : if (debug) then
755 : if (ssp .eq. 1 .and. itrtensor .eq. 1) then
756 : write(31,"(3i5,4E20.12)") iomega, ie2, ie1, &
757 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp), j01, &
758 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01, &
759 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01*occ_e2*(1.0_dp-occ_e1)*(occ_omega+1)
760 : end if
761 : if (ssp .eq. 2 .and. itrtensor .eq. 1) then
762 : write(32,"(3i5,4E20.12)") iomega, ie2, ie1, &
763 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp), j01, &
764 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01, &
765 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01*occ_e2*(1.0_dp-occ_e1)*(occ_omega+1)
766 : end if
767 : if (ssp .eq. 3 .and. itrtensor .eq. 1) then
768 : write(33,"(3i5,4E20.12)") iomega, ie2, ie1, &
769 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp), j01, &
770 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01, &
771 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01*occ_e2*(1.0_dp-occ_e1)*(occ_omega+1)
772 : end if
773 : if (ssp .eq. 4 .and. itrtensor .eq. 1) then
774 : write(34,"(3i5,4E20.12)") iomega, ie2, ie1, &
775 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp), j01, &
776 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01, &
777 : & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01*occ_e2*(1.0_dp-occ_e1)*(occ_omega+1)
778 : end if
779 : end if
780 : !ENDDEBUG
781 : tointegq00_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)* &
782 0 : & occ_e2*(1.0_dp-occ_e1)*j00*(occ_omega+1)
783 : tointegq01_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)* &
784 0 : & occ_e2*(1.0_dp-occ_e1)*j01*(occ_omega+1)
785 : tointegq11_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)* &
786 0 : & occ_e2*(1.0_dp-occ_e1)*j11*(occ_omega+1)
787 : end if ! elph_ds%ep_lova
788 :
789 0 : tointegq00 = tointegq00 + tointegq00_1 + tointegq00_2
790 0 : tointegq01 = tointegq01 + tointegq01_1 + tointegq01_2
791 0 : tointegq11 = tointegq11 + tointegq11_1 + tointegq11_2
792 :
793 : end do ! ss' = 4
794 0 : integrand_q00(iomega) = integrand_q00(iomega) + elph_tr_ds%de_all(isppol,ie1)*tointegq00
795 0 : integrand_q01(iomega) = integrand_q01(iomega) + elph_tr_ds%de_all(isppol,ie1)*tointegq01
796 0 : integrand_q11(iomega) = integrand_q11(iomega) + elph_tr_ds%de_all(isppol,ie1)*tointegq11
797 : end do ! ie1 ~ 20
798 0 : omega=omega+domega
799 : q00(itemp,icomp,jcomp,isppol) = q00(itemp,icomp,jcomp,isppol) +&
800 0 : & domega*integrand_q00(iomega)
801 : q01(itemp,icomp,jcomp,isppol) = q01(itemp,icomp,jcomp,isppol) +&
802 0 : & domega*integrand_q01(iomega)
803 : q11(itemp,icomp,jcomp,isppol) = q11(itemp,icomp,jcomp,isppol) +&
804 0 : & domega*integrand_q11(iomega)
805 : end do ! omega ~ 400
806 :
807 : q00(itemp,icomp,jcomp,isppol)=q00(itemp,icomp,jcomp,isppol)* &
808 0 : & ucvol*two_pi*elph_tr_ds%dos_n0(itemp,isppol)/(kb_HaK*Temp)
809 : q01(itemp,icomp,jcomp,isppol)=q01(itemp,icomp,jcomp,isppol)* &
810 0 : & ucvol*two_pi*elph_tr_ds%dos_n0(itemp,isppol)/(kb_HaK*Temp)
811 : q11(itemp,icomp,jcomp,isppol)=q11(itemp,icomp,jcomp,isppol)* &
812 0 : & ucvol*two_pi*elph_tr_ds%dos_n0(itemp,isppol)/(kb_HaK*Temp)
813 :
814 0 : rho = 0.5_dp*q00(itemp,icomp,jcomp,isppol)
815 : ! is tau energy dependent?
816 0 : tau = 2.0d0*ucvol/(q00(itemp,icomp,jcomp,isppol)*elph_tr_ds%dos_n0(itemp,isppol)*tmp_veloc_sq0)
817 0 : write(unit_rho,'(4D20.10)')temp,rho,rho*chgu,rho/temp
818 0 : write(unit_tau,'(3D20.10)')temp,tau,tau*Time_Sec
819 0 : rho_T(itemp)=rho
820 : end do ! temperature = 1?
821 0 : write(unit_rho,*)
822 0 : write(unit_tau,*)
823 :
824 : end do ! jcomp = 3
825 : end do ! icomp = 3
826 : end do ! isppol = 2
827 :
828 : !-----------------------------
829 :
830 0 : seebeck = zero
831 : !rho_nm = zero
832 :
833 : !do isppol=1,elph_ds%nsppol
834 : !do itemp=1,ntemper
835 : !q11_inv(:,:)=q11(itemp,:,:,isppol)
836 : !call matrginv(q11_inv,3,3)
837 : !do icomp=1,3
838 : !do jcomp=1,3
839 : !do kcomp=1,3
840 : !seebeck(isppol,itemp,icomp,jcomp) = seebeck(isppol,itemp,icomp,jcomp) + &
841 : !& q01(itemp,icomp,kcomp,isppol)*q11_inv(kcomp,jcomp)
842 : !end do
843 : !end do
844 : !end do
845 : !end do
846 : !end do
847 :
848 0 : do isppol=1,elph_ds%nsppol
849 0 : do itemp=1,ntemper
850 0 : q11_inv(:,:)=q11(itemp,:,:,isppol)
851 0 : call matrginv(q11_inv,3,3)
852 : call DGEMM('N','N',3,3,3,one,q01(itemp,:,:,isppol),3,q11_inv,&
853 0 : & 3,zero,seebeck(isppol,itemp,:,:),3)
854 : ! call DGEMM('N','N',3,3,3,one,seebeck(isppol,itemp,:,:),3,q01(itemp,:,:,isppol),&
855 : ! & 3,zero,rho_nm(isppol,itemp,:,:),3)
856 : end do
857 : end do
858 : pref_s = pi/sqrt(3.0_dp)
859 0 : seebeck=pref_s*seebeck
860 :
861 : !fullq = zero
862 : !do icomp=1,3
863 : !do jcomp=1,3
864 : !fullq(icomp,jcomp) = q00(1,icomp,jcomp,1)
865 : !end do
866 : !end do
867 : !do icomp=1,3
868 : !do jcomp=4,6
869 : !fullq(icomp,jcomp) = q01(1,icomp,jcomp-3,1)
870 : !end do
871 : !end do
872 : !do icomp=4,6
873 : !do jcomp=1,3
874 : !fullq(icomp,jcomp) = q01(1,icomp-3,jcomp,1)
875 : !end do
876 : !end do
877 : !do icomp=4,6
878 : !do jcomp=4,6
879 : !fullq(icomp,jcomp) = q11(1,icomp-3,jcomp-3,1)
880 : !end do
881 : !end do
882 : !write(*,*)' fullq'
883 : !write(*,"(6E20.12)") (fullq(1,jcomp),jcomp=1,6)
884 : !write(*,"(6E20.12)") (fullq(2,jcomp),jcomp=1,6)
885 : !write(*,"(6E20.12)") (fullq(3,jcomp),jcomp=1,6)
886 : !write(*,"(6E20.12)") (fullq(4,jcomp),jcomp=1,6)
887 : !write(*,"(6E20.12)") (fullq(5,jcomp),jcomp=1,6)
888 : !write(*,"(6E20.12)") (fullq(6,jcomp),jcomp=1,6)
889 0 : write(message,'(a)') 'q00:'
890 0 : call wrtout(std_out,message,'COLL')
891 0 : write(message,'(3E20.12)') (q00(1,1,jcomp,1),jcomp=1,3)
892 0 : call wrtout(std_out,message,'COLL')
893 0 : write(message,'(3E20.12)') (q00(1,2,jcomp,1),jcomp=1,3)
894 0 : call wrtout(std_out,message,'COLL')
895 0 : write(message,'(3E20.12)') (q00(1,3,jcomp,1),jcomp=1,3)
896 0 : call wrtout(std_out,message,'COLL')
897 0 : write(message,'(a)') 'q01:'
898 0 : call wrtout(std_out,message,'COLL')
899 0 : write(message,'(3E20.12)') (q01(1,1,jcomp,1),jcomp=1,3)
900 0 : call wrtout(std_out,message,'COLL')
901 0 : write(message,'(3E20.12)') (q01(1,2,jcomp,1),jcomp=1,3)
902 0 : call wrtout(std_out,message,'COLL')
903 0 : write(message,'(3E20.12)') (q01(1,3,jcomp,1),jcomp=1,3)
904 0 : call wrtout(std_out,message,'COLL')
905 0 : write(message,'(a)') 'q11, q11_inv:'
906 0 : call wrtout(std_out,message,'COLL')
907 0 : write(message,'(6E20.12)') (q11(1,1,jcomp,1),jcomp=1,3),(q11_inv(1,jcomp),jcomp=1,3)
908 0 : call wrtout(std_out,message,'COLL')
909 0 : write(message,'(6E20.12)') (q11(1,2,jcomp,1),jcomp=1,3),(q11_inv(2,jcomp),jcomp=1,3)
910 0 : call wrtout(std_out,message,'COLL')
911 0 : write(message,'(6E20.12)') (q11(1,3,jcomp,1),jcomp=1,3),(q11_inv(3,jcomp),jcomp=1,3)
912 0 : call wrtout(std_out,message,'COLL')
913 : !q11_inv = zero
914 : !do icomp = 1, 3
915 : !q11_inv(icomp,icomp) = 2.0_dp
916 : !end do
917 :
918 : !call matrginv(fullq,6,6)
919 :
920 : !do isppol=1,elph_ds%nsppol
921 : !do itemp=1,ntemper
922 : !rho_nm(isppol,itemp,:,:) = q00(itemp,:,:,isppol) - rho_nm(isppol,itemp,:,:)
923 : !end do
924 : !end do
925 : !rho_nm = 0.5_dp*rho_nm
926 :
927 : !Output of Seebeck coefficient
928 0 : do isppol=1,elph_ds%nsppol
929 0 : do icomp=1,3
930 0 : do jcomp=1,3
931 0 : itrtensor=(icomp-1)*3+jcomp
932 0 : write(unit_sbk,*) '# Seebeck for isppol, itrten = ', isppol, itrtensor
933 : ! write(88,*) '# Rho for isppol, itrten = ', isppol, itrtensor
934 : ! write(89,*) '# Rho for isppol, itrten = ', isppol, itrtensor
935 0 : do itemp=1,ntemper
936 0 : Temp=tempermin+temperinc*dble(itemp)
937 0 : write(unit_sbk,'(3D20.10)')temp, seebeck(isppol,itemp,icomp,jcomp), seebeck(isppol,itemp,icomp,jcomp)*kb_HaK
938 : ! write(88,'(3D20.10)')temp, rho_nm(isppol,itemp,icomp,jcomp), rho_nm(isppol,itemp,icomp,jcomp)*chgu
939 : ! write(89,'(3D20.10)')temp, 0.5_dp/fullq(1,1), 0.5_dp*chgu/fullq(1,1)
940 : end do ! temperature
941 0 : write(unit_sbk,*)
942 : ! write(88,*)
943 : ! write(89,*)
944 : end do ! jcomp
945 : end do ! icomp
946 : end do ! isppol
947 :
948 : !Get thermal resistivity, based on eqn. (52) in Allen's PRB 17, 3725 (1978) [[cite:Allen1978]]
949 : !WARNING: before 6.13.1 the thermal resistivity and Lorentz number were not in
950 : !atomic units, BUT the SI units are good.
951 0 : pref_w = 3.0_dp/(2.0_dp*pi**2.0d0)
952 0 : do isppol=1,elph_ds%nsppol
953 0 : do icomp=1, 3
954 0 : do jcomp=1, 3
955 0 : itrtensor=(icomp-1)*3+jcomp
956 :
957 0 : write(unit_therm,*) '# Thermal resistivity for isppol, itrten= ', isppol
958 0 : write(unit_lor,*) '# Lorentz coefficient for isppol, itrten= ', isppol
959 :
960 0 : do itemp=1,ntemper
961 :
962 0 : Temp=tempermin + temperinc*dble(itemp)
963 :
964 0 : wtherm = pref_w*q11(itemp,icomp,jcomp,isppol)/(kb_HaK*Temp)
965 :
966 : ! write(unit_therm,'(5D20.10)')temp,wtherm,1./wtherm,wtherm/3.4057d9,1./(wtherm) *3.4057d9
967 0 : write(unit_therm,'(5D20.10)')temp,wtherm,1.0_dp/wtherm,wtherm*chwu,1.0_dp/(wtherm*chwu)
968 :
969 0 : lorentz=rho_T(itemp)/(wtherm*kb_HaK*Temp)
970 0 : write(unit_lor,*)temp,lorentz,lor0
971 :
972 : end do
973 0 : write(unit_therm,*)
974 0 : write(unit_lor,*)
975 : end do ! jcomp
976 : end do ! icomp
977 : end do ! isppol
978 :
979 :
980 0 : ABI_FREE(phfrq)
981 0 : ABI_FREE(displ)
982 0 : ABI_FREE(pheigvec)
983 0 : ABI_FREE(integrand_q00)
984 0 : ABI_FREE(integrand_q01)
985 0 : ABI_FREE(integrand_q11)
986 0 : ABI_FREE(q00)
987 0 : ABI_FREE(q01)
988 0 : ABI_FREE(q11)
989 0 : ABI_FREE(seebeck)
990 0 : ABI_FREE(rho_T)
991 0 : ABI_FREE(integrho)
992 0 : ABI_FREE(tointegrho)
993 :
994 0 : close (unit=unit_lor)
995 0 : close (unit=unit_rho)
996 0 : close (unit=unit_tau)
997 0 : close (unit=unit_sbk)
998 0 : close (unit=unit_therm)
999 :
1000 0 : ABI_FREE(elph_ds%k_fine%wtq)
1001 0 : ABI_FREE(elph_ds%k_phon%wtq)
1002 :
1003 0 : ABI_FREE(elph_tr_ds%a2f_1d_tr)
1004 :
1005 0 : ABI_FREE(elph_tr_ds%gamma_qpt_tr)
1006 0 : ABI_FREE(elph_tr_ds%gamma_rpt_tr)
1007 0 : write(std_out,*) ' mka2f_tr : end '
1008 :
1009 0 : end subroutine mka2f_tr
1010 : !!***
1011 :
1012 :
1013 : !!****f* ABINIT/mka2f_tr_lova
1014 : !!
1015 : !! NAME
1016 : !! mka2f_tr_lova
1017 : !!
1018 : !! FUNCTION
1019 : !! calculates the FS averaged Transport alpha^2F_tr alpha^2F_trout alpha^2F_trin functions
1020 : !! calculates and outputs the associated electrical and thermal conductivities
1021 : !! for the first task: copied from mka2F
1022 : !!
1023 : !! INPUTS
1024 : !! crystal<crystal_t>=data type gathering info on the crystalline structure.
1025 : !! Ifc<ifc_type>=Object containing the interatomic force constants.
1026 : !! elph_ds
1027 : !! elph_ds%gkk2 = gkk2 matrix elements on full FS grid for each phonon mode
1028 : !! elph_ds%nbranch = number of phonon branches = 3*natom
1029 : !! elph_ds%nFSband = number of bands included in the FS integration
1030 : !! elph_ds%k_fine%nkpt = number of kpts included in the FS integration
1031 : !! elph_ds%k_fine%wtk = integration weights on the FS
1032 : !! delph_ds%n0 = DOS at the Fermi level calculated from the k_fine integration weights
1033 : !! elph_ds%k_fine%kpt = coordinates of all FS kpoints
1034 : !! mustar = coulomb pseudopotential parameter
1035 : !! eventually for 2 spin channels
1036 : !! ntemper = number of temperature points to calculate, from tempermin to tempermin+ntemper*temperinc
1037 : !! tempermin = minimum temperature at which resistivity etc are calculated (in K)
1038 : !! temperinc = interval for temperature grid on which resistivity etc are calculated (in K)
1039 : !!
1040 : !! OUTPUT
1041 : !! elph_ds
1042 : !!
1043 : !! NOTES
1044 : !! copied from ftiaf9.f
1045 : !!
1046 : !! SOURCE
1047 :
1048 5 : subroutine mka2f_tr_lova(crystal,ifc,elph_ds,ntemper,tempermin,temperinc,elph_tr_ds)
1049 :
1050 : !Arguments ------------------------------------
1051 : !scalars
1052 : integer,intent(in) :: ntemper
1053 : real(dp),intent(in) :: tempermin,temperinc
1054 : type(crystal_t),intent(in) :: crystal
1055 : type(ifc_type),intent(in) :: ifc
1056 : type(elph_tr_type),intent(inout) :: elph_tr_ds
1057 : type(elph_type),intent(inout) :: elph_ds
1058 :
1059 : !Local variables -------------------------
1060 : !x =w/(2kbT)
1061 : !scalars
1062 : integer :: iFSqpt,ibranch,iomega,isppol,jbranch,nerr
1063 : integer :: unit_a2f_tr, unit_a2f_trout, unit_a2f_trin, natom
1064 : integer :: idir, iatom, k1, kdir,unit_lor,unit_rho,unit_tau,unit_therm
1065 : integer :: itemp,nrpt,itrtensor, icomp, jcomp
1066 : real(dp) :: Temp,chgu,femto,diagerr,firh,firhT,gaussfactor,domega
1067 : real(dp) :: firh_tau,firhT_tau ! added by BX to get Tau
1068 : real(dp) :: a2fprefactor_in, temp_in
1069 : real(dp) :: a2fprefactor_out, temp_out
1070 : real(dp) :: gaussprefactor,gaussval,lambda_tr,lor0,lorentz,maxerr,maxx,omega
1071 : real(dp) :: rho,tau,tolexp,wtherm,xtr,xx
1072 : real(dp) :: lambda_tr_trace,omega_min, omega_max,qnorm2,spinfact
1073 : character(len=500) :: message
1074 : character(len=fnlen) :: fname
1075 : !arrays
1076 : complex(dp),parameter :: c0=dcmplx(0.d0,0.d0),c1=dcmplx(1.d0,0.d0)
1077 : real(dp) :: gprimd(3,3)
1078 10 : real(dp) :: eigval_in(elph_ds%nbranch)
1079 10 : real(dp) :: eigval_out(elph_ds%nbranch)
1080 10 : real(dp) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
1081 10 : real(dp) :: gam_now_in (2,elph_ds%nbranch*elph_ds%nbranch)
1082 10 : real(dp) :: gam_now_out(2,elph_ds%nbranch*elph_ds%nbranch)
1083 10 : real(dp) :: tmpa2f_in (elph_ds%na2f)
1084 10 : real(dp) :: tmpa2f_out(elph_ds%na2f)
1085 10 : real(dp) :: tmpgam1(2,elph_ds%nbranch,elph_ds%nbranch)
1086 10 : real(dp) :: tmpgam2(2,elph_ds%nbranch,elph_ds%nbranch)
1087 5 : real(dp),allocatable :: phfrq(:,:)
1088 5 : real(dp),allocatable :: displ(:,:,:,:)
1089 5 : real(dp),allocatable :: pheigvec(:,:)
1090 5 : real(dp),allocatable :: integrho(:),integtau(:),tointegrho(:),tointega2f(:),tointegtau(:)
1091 5 : real(dp),allocatable :: rho_T(:),tau_T(:)
1092 5 : real(dp),allocatable :: coskr(:,:)
1093 5 : real(dp),allocatable :: sinkr(:,:)
1094 : !real(dp),allocatable :: gam_rpt(:,:,:)
1095 : ! *********************************************************************
1096 :
1097 : !calculate a2f_tr for frequencies between 0 and omega_max
1098 5 : write(std_out,*) 'mka2f_tr_lova : enter '
1099 : !
1100 : !MG: the step should be calculated locally using nomega and the extrema of the spectrum.
1101 : !One should not rely on previous calls for the setup of elph_ds%domega
1102 : !I will remove elph_ds%domega since mka2f.F90 will become a method of gamma_t
1103 5 : domega =elph_ds%domega
1104 :
1105 : ! Number of points for FFT interpolation
1106 5 : nrpt = ifc%nrpt
1107 5 : natom = crystal%natom
1108 65 : gprimd = crystal%gprimd
1109 :
1110 20 : ABI_MALLOC(elph_tr_ds%a2f_1d_tr,(elph_ds%na2f,9,elph_ds%nsppol,1,1,1))
1111 15 : ABI_MALLOC(elph_tr_ds%a2f_1d_trin,(elph_ds%na2f,9,elph_ds%nsppol))
1112 15 : ABI_MALLOC(elph_tr_ds%a2f_1d_trout,(elph_ds%na2f,9,elph_ds%nsppol))
1113 :
1114 : !! defaults for number of temperature steps and max T (all in Kelvin...)
1115 : !ntemper=1000
1116 : !tempermin=zero
1117 : !temperinc=one
1118 15 : ABI_MALLOC(rho_T,(ntemper))
1119 10 : ABI_MALLOC(tau_T,(ntemper))
1120 :
1121 :
1122 : !tolerance on gaussian being = 0
1123 5 : tolexp = 1.d-100
1124 5 : maxx = sqrt(-log(tolexp))
1125 5 : lor0=(pi*kb_HaK)**2/3.
1126 :
1127 : !maximum value of frequency (a grid has to be chosen for the representation of alpha^2 F)
1128 : !WARNING! supposes this value has been set in mkelph_linwid.
1129 :
1130 5 : gaussprefactor = sqrt(piinv) / elph_ds%a2fsmear
1131 5 : gaussfactor = one / elph_ds%a2fsmear
1132 :
1133 : !spinfact should be 1 for a normal non sppol calculation without spinorbit
1134 : !for spinors it should also be 1 as bands are twice as numerous but n0 has been divided by 2
1135 : !for sppol 2 it should be 0.5 as we have 2 spin channels to sum
1136 5 : spinfact = one / elph_ds%nsppol !/ elph_ds%nspinor
1137 :
1138 : !ENDMG
1139 :
1140 18070 : elph_tr_ds%a2f_1d_tr = zero
1141 18055 : elph_tr_ds%a2f_1d_trin = zero
1142 18055 : elph_tr_ds%a2f_1d_trout = zero
1143 :
1144 5 : maxerr=0.
1145 5 : nerr=0
1146 :
1147 20 : ABI_MALLOC(phfrq,(elph_ds%nbranch, elph_ds%k_fine%nkpt))
1148 25 : ABI_MALLOC(displ,(2, elph_ds%nbranch, elph_ds%nbranch, elph_ds%k_fine%nkpt))
1149 20 : ABI_MALLOC(pheigvec,(2*elph_ds%nbranch*elph_ds%nbranch, elph_ds%k_fine%nkpt))
1150 :
1151 429 : do iFSqpt=1,elph_ds%k_fine%nkpt
1152 429 : call ifc%fourq(crystal,elph_ds%k_fine%kpt(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt),out_eigvec=pheigvec(:,iFSqpt))
1153 : end do
1154 :
1155 2133 : omega_min = minval(phfrq)
1156 2133 : omega_max = maxval(phfrq)
1157 :
1158 20 : ABI_MALLOC(coskr, (elph_ds%k_fine%nkpt,nrpt))
1159 15 : ABI_MALLOC(sinkr, (elph_ds%k_fine%nkpt,nrpt))
1160 5 : call ftgam_init(Ifc%gprim, elph_ds%k_fine%nkpt, nrpt, elph_ds%k_fine%kpt, Ifc%rpt, coskr, sinkr)
1161 :
1162 10 : do isppol=1,elph_ds%nsppol
1163 :
1164 : ! loop over qpoint in full kpt grid (presumably dense)
1165 434 : do iFSqpt=1,elph_ds%k_fine%nkpt
1166 1696 : qnorm2 = sum(elph_ds%k_fine%kpt(:,iFSqpt)**2)
1167 : ! if (flag_to_exclude_soft_modes = .false.) qnorm2 = zero
1168 4245 : do itrtensor=1,9
1169 : ! Do FT from real-space gamma grid to 1 qpt.
1170 :
1171 3816 : if (elph_ds%ep_int_gkk == 1) then
1172 0 : gam_now_in(:,:) = elph_tr_ds%gamma_qpt_trin(:,itrtensor,:,isppol,iFSqpt)
1173 0 : gam_now_out(:,:) = elph_tr_ds%gamma_qpt_trout(:,itrtensor,:,isppol,iFSqpt)
1174 : else
1175 : call ftgam(Ifc%wghatm,gam_now_in, elph_tr_ds%gamma_rpt_trin(:,itrtensor,:,isppol,:),natom,1,nrpt,0,&
1176 9908856 : & coskr(iFSqpt,:), sinkr(iFSqpt,:))
1177 : call ftgam(Ifc%wghatm,gam_now_out,elph_tr_ds%gamma_rpt_trout(:,itrtensor,:,isppol,:),natom,1,nrpt,0,&
1178 9908856 : & coskr(iFSqpt,:), sinkr(iFSqpt,:))
1179 : end if
1180 :
1181 : ! Diagonalize gamma matrix at this qpoint (complex matrix).
1182 :
1183 : ! if ep_scalprod==0 we have to dot in the displacement vectors here
1184 3816 : if (elph_ds%ep_scalprod==0) then
1185 :
1186 60264 : displ_red(:,:,:) = zero
1187 7776 : do jbranch=1,elph_ds%nbranch
1188 13608 : do iatom=1,natom
1189 29160 : do idir=1,3
1190 17496 : ibranch=idir+3*(iatom-1)
1191 75816 : do kdir=1,3
1192 52488 : k1 = kdir+3*(iatom-1)
1193 : displ_red(1,ibranch,jbranch) = displ_red(1,ibranch,jbranch) + &
1194 52488 : & gprimd(kdir,idir)*displ(1,k1,jbranch,iFSqpt)
1195 : displ_red(2,ibranch,jbranch) = displ_red(2,ibranch,jbranch) + &
1196 69984 : & gprimd(kdir,idir)*displ(2,k1,jbranch,iFSqpt)
1197 : end do
1198 : end do
1199 : end do
1200 : end do
1201 :
1202 7776 : tmpgam2 = reshape (gam_now_in, (/2,elph_ds%nbranch,elph_ds%nbranch/))
1203 1944 : call gam_mult_displ(elph_ds%nbranch, displ_red, tmpgam2, tmpgam1)
1204 7776 : do jbranch=1,elph_ds%nbranch
1205 7776 : eigval_in(jbranch) = tmpgam1(1, jbranch, jbranch)
1206 : end do
1207 :
1208 7776 : tmpgam2 = reshape (gam_now_out, (/2,elph_ds%nbranch,elph_ds%nbranch/))
1209 1944 : call gam_mult_displ(elph_ds%nbranch, displ_red, tmpgam2, tmpgam1)
1210 7776 : do jbranch=1,elph_ds%nbranch
1211 7776 : eigval_out(jbranch) = tmpgam1(1, jbranch, jbranch)
1212 : end do
1213 :
1214 1872 : else if (elph_ds%ep_scalprod == 1) then
1215 :
1216 : !
1217 : ! NOTE: in these calls gam_now and pheigvec do not have the right rank, but blas usually does not care
1218 : call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, c1, gam_now_in, 3*natom,&
1219 1872 : & pheigvec(:,iFSqpt), 3*natom, c0, tmpgam1, 3*natom)
1220 : call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, c1, pheigvec(:,iFSqpt), 3*natom,&
1221 1872 : & tmpgam1, 3*natom, c0, tmpgam2, 3*natom)
1222 1872 : diagerr = zero
1223 :
1224 11376 : do ibranch=1,elph_ds%nbranch
1225 9504 : eigval_in(ibranch) = tmpgam2(1,ibranch,ibranch)
1226 30672 : do jbranch=1,ibranch-1
1227 30672 : diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch))
1228 : end do
1229 32544 : do jbranch=ibranch+1,elph_ds%nbranch
1230 30672 : diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch))
1231 : end do
1232 : end do
1233 1872 : if (diagerr > tol12) then
1234 1036 : nerr=nerr+1
1235 1036 : maxerr=max(diagerr, maxerr)
1236 : end if
1237 :
1238 : call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, c1, gam_now_out, 3*natom,&
1239 1872 : & pheigvec(:,iFSqpt), 3*natom, c0, tmpgam1, 3*natom)
1240 : call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, c1, pheigvec(:,iFSqpt), 3*natom,&
1241 1872 : & tmpgam1, 3*natom, c0, tmpgam2, 3*natom)
1242 1872 : diagerr = zero
1243 :
1244 11376 : do ibranch=1,elph_ds%nbranch
1245 9504 : eigval_out(ibranch) = tmpgam2(1,ibranch,ibranch)
1246 30672 : do jbranch=1,ibranch-1
1247 30672 : diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch))
1248 : end do
1249 32544 : do jbranch=ibranch+1,elph_ds%nbranch
1250 30672 : diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch))
1251 : end do
1252 : end do
1253 1872 : if (diagerr > tol12) then
1254 1036 : nerr=nerr+1
1255 1036 : maxerr=max(diagerr, maxerr)
1256 : end if
1257 : end if
1258 : ! end ep_scalprod if
1259 :
1260 : ! Add all contributions from the phonon modes at this qpoint to
1261 : ! a2f and the phonon dos.
1262 19576 : do ibranch=1,elph_ds%nbranch
1263 : ! if (abs(phfrq(ibranch,iFSqpt)) < tol10) then
1264 15336 : if ( abs(phfrq(ibranch,iFSqpt)) < tol7 .or. &
1265 : & (phfrq(ibranch,iFSqpt) < tol4 .and. qnorm2 > 0.03 )) then !
1266 : ! note: this should depend on the velocity of sound, to accept acoustic
1267 : ! modes!
1268 : a2fprefactor_in = zero
1269 : a2fprefactor_out= zero
1270 : else
1271 14067 : a2fprefactor_in = eigval_in (ibranch)/(two_pi*abs(phfrq(ibranch,iFSqpt))*elph_ds%n0(isppol))
1272 14067 : a2fprefactor_out = eigval_out(ibranch)/(two_pi*abs(phfrq(ibranch,iFSqpt))*elph_ds%n0(isppol))
1273 : end if
1274 :
1275 15336 : omega = omega_min
1276 6149736 : tmpa2f_in (:) = zero
1277 6149736 : tmpa2f_out(:) = zero
1278 6149736 : do iomega=1,elph_ds%na2f
1279 6134400 : xx = (omega-phfrq(ibranch,iFSqpt))*gaussfactor
1280 6134400 : gaussval = gaussprefactor*exp(-xx*xx)
1281 :
1282 6134400 : temp_in = gaussval*a2fprefactor_in
1283 6134400 : temp_out = gaussval*a2fprefactor_out
1284 :
1285 6134400 : if (dabs(temp_in) < 1.0d-50) temp_in = zero
1286 6134400 : if (dabs(temp_out) < 1.0d-50) temp_out = zero
1287 6134400 : tmpa2f_in (iomega) = tmpa2f_in (iomega) + temp_in
1288 6134400 : tmpa2f_out(iomega) = tmpa2f_out(iomega) + temp_out
1289 6149736 : omega = omega+domega
1290 : end do
1291 :
1292 6149736 : elph_tr_ds%a2f_1d_trin (:,itrtensor,isppol) = elph_tr_ds%a2f_1d_trin (:,itrtensor,isppol) + tmpa2f_in(:)
1293 6153552 : elph_tr_ds%a2f_1d_trout(:,itrtensor,isppol) = elph_tr_ds%a2f_1d_trout(:,itrtensor,isppol) + tmpa2f_out(:)
1294 :
1295 : end do ! end ibranch do
1296 : end do ! end itrtensor do
1297 : end do ! end iFSqpt do
1298 : end do ! end isppol
1299 :
1300 5 : ABI_FREE(coskr)
1301 5 : ABI_FREE(sinkr)
1302 :
1303 : !second 1 / elph_ds%k_fine%nkpt factor for the integration weights
1304 18055 : elph_tr_ds%a2f_1d_trin = elph_tr_ds%a2f_1d_trin / elph_ds%k_fine%nkpt
1305 18055 : elph_tr_ds%a2f_1d_trout = elph_tr_ds%a2f_1d_trout / elph_ds%k_fine%nkpt
1306 :
1307 5 : if (elph_ds%ep_scalprod == 1) then
1308 4 : write(std_out,*) 'mka2f_tr_lova: errors in diagonalization of gamma_tr with phon eigenvectors: ', nerr,maxerr
1309 : end if
1310 :
1311 18055 : elph_tr_ds%a2f_1d_tr(:,:,:,1,1,1) = elph_tr_ds%a2f_1d_trout(:,:,:) - elph_tr_ds%a2f_1d_trin(:,:,:)
1312 :
1313 : !output the elph_tr_ds%a2f_1d_tr
1314 5 : fname = trim(elph_ds%elph_base_name) // '_A2F_TR'
1315 5 : if (open_file (fname,message,newunit=unit_a2f_tr,status='unknown') /= 0) then
1316 0 : ABI_ERROR(message)
1317 : end if
1318 :
1319 5 : fname = trim(elph_ds%elph_base_name) // '_A2F_TRIN'
1320 5 : if (open_file(fname,message,newunit=unit_a2f_trin,status='unknown') /= 0) then
1321 0 : ABI_ERROR(message)
1322 : end if
1323 :
1324 5 : fname = trim(elph_ds%elph_base_name) // '_A2F_TROUT'
1325 5 : if (open_file (fname,message,newunit=unit_a2f_trout,status='unknown') /=0) then
1326 0 : ABI_ERROR(message)
1327 : end if
1328 :
1329 5 : write (unit_a2f_tr,'(a)') '#'
1330 5 : write (unit_a2f_tr,'(a)') '# ABINIT package : a2f_tr file'
1331 5 : write (unit_a2f_tr,'(a)') '#'
1332 5 : write (unit_a2f_tr,'(a)') '# a2f_tr function integrated over the FS. omega in a.u.'
1333 5 : write (unit_a2f_tr,'(a,I10)') '# number of kpoints integrated over : ', elph_ds%k_fine%nkpt
1334 5 : write (unit_a2f_tr,'(a,I10)') '# number of energy points : ',elph_ds%na2f
1335 5 : write (unit_a2f_tr,'(a,E16.6,a,E16.6,a)') '# between omega_min = ', omega_min,' Ha and omega_max = ', omega_max, ' Ha'
1336 5 : write (unit_a2f_tr,'(a,E16.6)') '# and the smearing width for gaussians is ', elph_ds%a2fsmear
1337 5 : write (unit_a2f_tr,'(a)') '#'
1338 :
1339 5 : write (unit_a2f_trin,'(a)') '#'
1340 5 : write (unit_a2f_trin,'(a)') '# ABINIT package : a2f_trin file'
1341 5 : write (unit_a2f_trin,'(a)') '#'
1342 5 : write (unit_a2f_trin,'(a)') '# a2f_trin function integrated over the FS. omega in a.u.'
1343 5 : write (unit_a2f_trin,'(a,I10)') '# number of kpoints integrated over : ', elph_ds%k_fine%nkpt
1344 5 : write (unit_a2f_trin,'(a,I10)') '# number of energy points : ',elph_ds%na2f
1345 5 : write (unit_a2f_trin,'(a,E16.6,a,E16.6,a)') '# between omega_min = ', omega_min,' Ha and omega_max = ', omega_max, ' Ha'
1346 5 : write (unit_a2f_trin,'(a,E16.6)') '# and the smearing width for gaussians is ', elph_ds%a2fsmear
1347 5 : write (unit_a2f_trin,'(a)') '#'
1348 :
1349 5 : write (unit_a2f_trout,'(a)') '#'
1350 5 : write (unit_a2f_trout,'(a)') '# ABINIT package : a2f_trout file'
1351 5 : write (unit_a2f_trout,'(a)') '#'
1352 5 : write (unit_a2f_trout,'(a)') '# a2f_trout function integrated over the FS. omega in a.u.'
1353 5 : write (unit_a2f_trout,'(a,I10)') '# number of kpoints integrated over : ', elph_ds%k_fine%nkpt
1354 5 : write (unit_a2f_trout,'(a,I10)') '# number of energy points : ',elph_ds%na2f
1355 5 : write (unit_a2f_trout,'(a,E16.6,a,E16.6,a)') '# between omega_min = ', omega_min,' Ha and omega_max = ', omega_max, ' Ha'
1356 5 : write (unit_a2f_trout,'(a,E16.6)') '# and the smearing width for gaussians is ', elph_ds%a2fsmear
1357 5 : write (unit_a2f_trout,'(a)') '#'
1358 :
1359 : !done with header
1360 10 : do isppol=1,elph_ds%nsppol
1361 5 : write (unit_a2f_tr,'(a,E16.6)') '# The DOS at Fermi level is ', elph_ds%n0(isppol)
1362 5 : write (unit_a2f_trin,'(a,E16.6)') '# The DOS at Fermi level is ', elph_ds%n0(isppol)
1363 5 : write (unit_a2f_trout,'(a,E16.6)') '# The DOS at Fermi level is ', elph_ds%n0(isppol)
1364 : ! omega = zero
1365 5 : omega = omega_min
1366 2005 : do iomega=1,elph_ds%na2f
1367 2000 : write (unit_a2f_tr, '(10D16.6)') omega, elph_tr_ds%a2f_1d_tr (iomega,:,isppol,1,1,1)
1368 2000 : write (unit_a2f_trin, '(10D16.6)') omega, elph_tr_ds%a2f_1d_trin (iomega,:,isppol)
1369 2000 : write (unit_a2f_trout,'(10D16.6)') omega, elph_tr_ds%a2f_1d_trout(iomega,:,isppol)
1370 2005 : omega=omega+domega
1371 : end do
1372 5 : write (unit_a2f_tr,*)
1373 5 : write (unit_a2f_trin,*)
1374 10 : write (unit_a2f_trout,*)
1375 : end do !isppol
1376 :
1377 5 : close (unit=unit_a2f_tr)
1378 5 : close (unit=unit_a2f_trin)
1379 5 : close (unit=unit_a2f_trout)
1380 :
1381 : !calculation of transport properties
1382 15 : ABI_MALLOC(integrho,(elph_ds%na2f))
1383 10 : ABI_MALLOC(tointegrho,(elph_ds%na2f))
1384 10 : ABI_MALLOC(tointega2f,(elph_ds%na2f))
1385 10 : ABI_MALLOC(integtau,(elph_ds%na2f))
1386 10 : ABI_MALLOC(tointegtau,(elph_ds%na2f))
1387 :
1388 5 : fname = trim(elph_ds%elph_base_name) // '_RHO'
1389 5 : if (open_file(fname,message,newunit=unit_rho,status='unknown') /= 0) then
1390 0 : ABI_ERROR(message)
1391 : end if
1392 :
1393 : !print header to resistivity file
1394 5 : write (unit_rho,*) '# Resistivity as a function of temperature.'
1395 5 : write (unit_rho,*) '# the formalism is isotropic, so non-cubic crystals may be wrong'
1396 5 : write (unit_rho,*) '# '
1397 5 : write (unit_rho,*) '# Columns are: '
1398 5 : write (unit_rho,*) '# temperature[K] rho[au] rho [SI] rho/temp [au]'
1399 5 : write (unit_rho,*) '# '
1400 :
1401 5 : fname = trim(elph_ds%elph_base_name) // '_TAU'
1402 5 : if (open_file(fname,message,newunit=unit_tau,status='unknown') /= 0) then
1403 0 : ABI_ERROR(message)
1404 : end if
1405 :
1406 : !print header to relaxation time file
1407 5 : write (unit_tau,*) '# Relaxation time as a function of temperature.'
1408 5 : write (unit_tau,*) '# the formalism is isotropic, so non-cubic crystals may be wrong'
1409 5 : write (unit_tau,*) '# '
1410 5 : write (unit_tau,*) '# Columns are: '
1411 5 : write (unit_tau,*) '# temperature[K] tau[au] tau [femtosecond] '
1412 5 : write (unit_tau,*) '# '
1413 :
1414 5 : fname = trim(elph_ds%elph_base_name) // '_WTH'
1415 5 : if (open_file(fname,message,newunit=unit_therm,status='unknown') /= 0) then
1416 0 : ABI_ERROR(message)
1417 : end if
1418 :
1419 : !print header to thermal conductivity file
1420 5 : write (unit_therm,'(a)') '# Thermal conductivity/resistivity as a function of temperature.'
1421 5 : write (unit_therm,'(a)') '# the formalism is isotropic, so non-cubic crystals may be wrong'
1422 5 : write (unit_therm,'(a)') '# '
1423 5 : write (unit_therm,'(a)') '# Columns are: '
1424 5 : write (unit_therm,'(a)') '# temperature[K] thermal rho[au] thermal cond [au] thermal rho [SI] thermal cond [SI]'
1425 5 : write (unit_therm,'(a)') '# '
1426 :
1427 5 : fname = trim(elph_ds%elph_base_name) // '_LOR'
1428 5 : if (open_file(fname,message,newunit=unit_lor,status='unknown') /= 0) then
1429 0 : ABI_ERROR(message)
1430 : end if
1431 :
1432 : !print header to lorentz file
1433 5 : write (unit_lor,*) '# Lorentz number as a function of temperature.'
1434 5 : write (unit_lor,*) '# the formalism is isotropic, so non-cubic crystals may be wrong'
1435 5 : write (unit_lor,*) '# '
1436 5 : write (unit_lor,*) '# Columns are: '
1437 5 : write (unit_lor,*) '# temperature[K] Lorentz number[au] Lorentz quantum = (pi*kb_HaK)**2/3'
1438 5 : write (unit_lor,*) '# '
1439 :
1440 10 : do isppol=1,elph_ds%nsppol
1441 5 : lambda_tr_trace = zero
1442 50 : do itrtensor=1,9
1443 45 : omega = omega_min
1444 18045 : tointega2f = zero
1445 18045 : do iomega=1,elph_ds%na2f
1446 18000 : if(omega<=0) then
1447 2682 : omega=omega+domega
1448 2682 : cycle
1449 : end if
1450 15318 : tointega2f(iomega)=elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1)/omega
1451 15363 : omega=omega+domega
1452 : end do
1453 :
1454 18045 : integrho = zero
1455 45 : call simpson_int(elph_ds%na2f,domega,tointega2f,integrho)
1456 45 : lambda_tr = two * spinfact * integrho(elph_ds%na2f)
1457 : write (message, '(a,2i3,a,es16.6)' )&
1458 45 : & ' mka2f_tr_lova : TRANSPORT lambda for isppol itrtensor', isppol, itrtensor, ' = ', lambda_tr
1459 45 : call wrtout(std_out,message,'COLL')
1460 50 : if (itrtensor == 1 .or. itrtensor == 5 .or. itrtensor == 9) lambda_tr_trace = lambda_tr_trace + lambda_tr
1461 : end do !end itrtensor do
1462 :
1463 5 : lambda_tr_trace = lambda_tr_trace / three
1464 : write (message, '(a,i3,a,es16.6)' )&
1465 5 : & ' mka2f_tr_lova: 1/3 trace of TRANSPORT lambda for isppol ', isppol, ' = ', lambda_tr_trace
1466 5 : call wrtout(std_out,message,'COLL')
1467 10 : call wrtout(ab_out,message,'COLL')
1468 : end do !end isppol do
1469 :
1470 : !constant to change units of rho from au to SI
1471 5 : chgu=2.173969d-7
1472 5 : femto=1.0d-15
1473 :
1474 10 : do isppol=1,elph_ds%nsppol
1475 25 : do icomp=1, 3
1476 65 : do jcomp=1, 3
1477 45 : itrtensor=(icomp-1)*3+jcomp
1478 :
1479 : ! prefactor for resistivity integral
1480 : ! firh=6.d0*pi*crystal%ucvol*kb_HaK/(elph_ds%n0(isppol)*elph_tr_ds%FSelecveloc_sq(isppol))
1481 : ! FIXME: check factor of 2 which is different from Savrasov paper. 6 below for thermal conductivity is correct.
1482 : firh=2.d0*pi*crystal%ucvol*kb_HaK/elph_ds%n0(isppol)/&
1483 45 : & sqrt(elph_tr_ds%FSelecveloc_sq(icomp,isppol)*elph_tr_ds%FSelecveloc_sq(jcomp,isppol))
1484 :
1485 : ! Add by BX to get Tau_elph
1486 45 : firh_tau = 2.0d0*pi*kb_HaK
1487 : ! End Adding
1488 :
1489 45 : write(unit_rho,*) '# Rho for isppol, itrten = ', isppol, itrtensor
1490 45 : write(unit_tau,*) '# Tau for isppol, itrten = ', isppol, itrtensor
1491 :
1492 : ! jmb
1493 18045 : tointegtau(:)=0.
1494 18045 : tointegrho(:)=0.
1495 495 : do itemp=1,ntemper ! runs over termperature in K
1496 450 : Temp=tempermin+temperinc*dble(itemp)
1497 450 : firhT=firh*Temp
1498 450 : firhT_tau=firh_tau*Temp
1499 450 : omega = omega_min
1500 180450 : do iomega=1,elph_ds%na2f
1501 180000 : if(omega<=0) then
1502 26820 : omega=omega+domega
1503 26820 : cycle
1504 : end if
1505 153180 : xtr=omega/(2*kb_HaK*Temp)
1506 153180 : if(xtr < log(huge(zero)*tol16)/2)then
1507 : tointegrho(iomega)=spinfact*firhT*omega*elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1) &
1508 153180 : & /(((2*Temp*kb_HaK)**2)*((exp(xtr)-exp(-xtr))/2)**2)
1509 : ! Add by BX to get Tau
1510 : tointegtau(iomega)=spinfact*firhT_tau*omega*elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1) &
1511 153180 : & /(((2*Temp*kb_HaK)**2)*((exp(xtr)-exp(-xtr))/2)**2)
1512 : else
1513 0 : tointegrho(iomega)=zero
1514 0 : tointegtau(iomega)=zero
1515 : end if
1516 153630 : omega=omega+domega
1517 : end do
1518 :
1519 450 : call simpson_int(elph_ds%na2f,domega,tointegrho,integrho)
1520 450 : call simpson_int(elph_ds%na2f,domega,tointegtau,integtau)
1521 450 : rho=integrho(elph_ds%na2f)
1522 450 : tau=1.0d99
1523 450 : if(dabs(integtau(elph_ds%na2f)) < tol7) then
1524 301 : write(message,'(a)') ' Cannot get a physical relaxation time '
1525 301 : ABI_WARNING(message)
1526 : else
1527 149 : tau=1.0d0/integtau(elph_ds%na2f)
1528 : end if
1529 : ! if(elph_ds%na2f < 350.0) then
1530 : ! tau=1.0d0/integtau(elph_ds%na2f)
1531 : ! end if
1532 450 : write(unit_rho,'(4D20.10)')temp,rho,rho*chgu,rho/temp
1533 450 : write(unit_tau,'(3D20.10)')temp,tau,tau*Time_Sec/femto
1534 450 : rho_T(itemp)=rho
1535 495 : tau_T(itemp)=tau
1536 : end do ! temperature
1537 45 : write(unit_rho,*)
1538 60 : write(unit_tau,*)
1539 :
1540 : end do ! jcomp
1541 : end do ! icomp
1542 : end do ! isppol
1543 :
1544 : !-----------------------------
1545 :
1546 :
1547 10 : do isppol=1,elph_ds%nsppol
1548 25 : do icomp=1, 3
1549 65 : do jcomp=1, 3
1550 45 : itrtensor=(icomp-1)*3+jcomp
1551 : ! prefactor for integral of thermal conductivity
1552 : ! firh=(18.*crystal%ucvol)/(pi*kb_HaK*elph_ds%n0(isppol)*elph_tr_ds%FSelecveloc_sq(isppol))
1553 : firh=(6.d0*crystal%ucvol)/(pi*kb_HaK*elph_ds%n0(isppol))/ &
1554 45 : & sqrt(elph_tr_ds%FSelecveloc_sq(icomp,isppol)*elph_tr_ds%FSelecveloc_sq(jcomp,isppol))
1555 :
1556 :
1557 45 : write(unit_therm,*) '# Thermal resistivity for isppol, itrten= ', isppol
1558 45 : write(unit_lor,*) '# Lorentz coefficient for isppol, itrten= ', isppol
1559 :
1560 18045 : tointegrho(:)=0.
1561 495 : do itemp=1,ntemper
1562 :
1563 450 : Temp=tempermin + temperinc*dble(itemp)
1564 450 : omega = omega_min
1565 180450 : do iomega=1,elph_ds%na2f
1566 180000 : if(omega<=0) then
1567 26820 : omega=omega+domega
1568 26820 : cycle
1569 : end if
1570 153180 : xtr=omega/(2*kb_HaK*Temp)
1571 153180 : if(xtr < log(huge(zero)*tol16)/2)then
1572 : tointegrho(iomega) = spinfact*xtr**2/omega*&
1573 : & ( elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1)+&
1574 : & 4*xtr**2*elph_tr_ds%a2f_1d_trout(iomega,itrtensor,isppol)/pi**2+ &
1575 : & 2*xtr**2*elph_tr_ds%a2f_1d_trin(iomega,itrtensor,isppol)/pi**2) &
1576 153180 : & /(((exp(xtr)-exp(-xtr))/2)**2)
1577 : else
1578 0 : tointegrho(iomega) = zero
1579 : end if
1580 153630 : omega=omega+domega
1581 : end do
1582 :
1583 450 : call simpson_int(elph_ds%na2f,domega,tointegrho,integrho)
1584 450 : wtherm=integrho(elph_ds%na2f)*firh
1585 :
1586 495 : if(abs(wtherm) > tol12)then
1587 422 : write(unit_therm,'(5D20.10)')temp,wtherm,1./wtherm,wtherm/3.4057d9,1./(wtherm) *3.4057d9
1588 :
1589 422 : lorentz=rho_T(itemp)/(wtherm*temp)
1590 422 : write(unit_lor,*)temp,lorentz,lor0
1591 : else
1592 28 : write(unit_therm,'(5D20.10)')temp,zero,huge(one),zero,huge(one)
1593 28 : write(unit_lor,*)temp,huge(one),lor0
1594 : end if
1595 :
1596 : end do
1597 45 : write(unit_therm,*)
1598 60 : write(unit_lor,*)
1599 : end do ! jcomp
1600 : end do ! icomp
1601 : end do !end isppol do
1602 :
1603 :
1604 5 : ABI_FREE(phfrq)
1605 5 : ABI_FREE(displ)
1606 5 : ABI_FREE(pheigvec)
1607 5 : ABI_FREE(rho_T)
1608 5 : ABI_FREE(tau_T)
1609 :
1610 5 : close (unit=unit_lor)
1611 5 : close (unit=unit_rho)
1612 5 : close (unit=unit_tau)
1613 5 : close (unit=unit_therm)
1614 :
1615 5 : ABI_FREE(integrho)
1616 5 : ABI_FREE(integtau)
1617 5 : ABI_FREE(tointega2f)
1618 5 : ABI_FREE(tointegrho)
1619 5 : ABI_FREE(tointegtau)
1620 5 : ABI_FREE(elph_tr_ds%a2f_1d_tr)
1621 5 : ABI_FREE(elph_tr_ds%a2f_1d_trin)
1622 5 : ABI_FREE(elph_tr_ds%a2f_1d_trout)
1623 :
1624 5 : ABI_FREE(elph_tr_ds%gamma_qpt_trin)
1625 5 : ABI_FREE(elph_tr_ds%gamma_qpt_trout)
1626 5 : ABI_FREE(elph_tr_ds%gamma_rpt_trin)
1627 5 : ABI_FREE(elph_tr_ds%gamma_rpt_trout)
1628 :
1629 : !DEBUG
1630 5 : write(std_out,*) ' mka2f_tr_lova : end '
1631 : !ENDDEBUG
1632 :
1633 5 : end subroutine mka2f_tr_lova
1634 : !!***
1635 :
1636 : !!****f* ABINIT/get_tau_k
1637 : !! NAME
1638 : !! get_tau_k
1639 : !!
1640 : !! FUNCTION
1641 : !! Calculate the k-dependent relaxation time due to EPC. Impelementation based
1642 : !! on derivation from Grmvall's book or
1643 : !! OD Restrepo's paper (PRB 94 212103 (2009) [[cite:Restrepo2009]])
1644 : !!
1645 : !! INPUTS
1646 : !! Cryst<crystal_t>=Info on the unit cell and on its symmetries.
1647 : !! Ifc<ifc_type>=Object containing the interatomic force constants.
1648 : !! elph_ds = elphon datastructure with data and dimensions
1649 : !! eigenGS = Ground State eigenvalues
1650 : !! max_occ = maximal occupancy for a band
1651 : !!
1652 : !! OUTPUT
1653 : !! tau_k(nsppol,nkptirr,nband)=mode relaxation time due to electron phonono coupling
1654 : !! rate_e(nene)= scattering rate due to electron phonono coupling vs. energy
1655 : !!
1656 : !! SOURCE
1657 :
1658 0 : subroutine get_tau_k(Cryst,ifc,Bst,elph_ds,elph_tr_ds,eigenGS,max_occ)
1659 :
1660 : !Arguments ------------------------------------
1661 : type(crystal_t),intent(in) :: Cryst
1662 : type(ifc_type),intent(in) :: ifc
1663 : type(ebands_t),intent(inout) :: Bst
1664 : type(elph_type),intent(inout) :: elph_ds
1665 : type(elph_tr_type), intent(inout) :: elph_tr_ds
1666 : real(dp),intent(in) :: max_occ
1667 : real(dp),intent(in) :: eigenGS(elph_ds%nband,elph_ds%k_phon%nkpt,elph_ds%nsppol)
1668 :
1669 : !Local variables-------------------------------
1670 : !scalars
1671 : character(len=500) :: message
1672 : character(len=fnlen) :: fname
1673 : integer :: ntemper,nsppol,nbranch,nband,natom
1674 : integer :: nkpt,nqpt,nkptirr,nqptirr,new_nkptirr
1675 : integer :: isppol,iFSkpt,iFSqpt,iqpt,iqpt_fullbz,imqpt_fullbz,ikpt_kpq,ikpt_kmq
1676 : integer :: iband,jband,jpband,jbeff,ibranch,jbranch,itemp
1677 : integer :: irec,ierr,nrpt,ik_this_proc
1678 : integer :: unit_tau,unit_invtau
1679 : integer :: nene,nene_all,iene,iene_fine,unit_taue,unit_mfp
1680 : integer :: icomp,jcomp,itensor
1681 : integer :: ikpt_irr,iomega,unit_cond,unit_therm,unit_sbk
1682 : integer :: nskip,nspline
1683 : real(dp) :: occ_omega,occ_e
1684 : real(dp) :: xx,Temp,therm_factor
1685 : real(dp) :: factor,dfermide
1686 : real(dp) :: e_k,rate_e,mfp_e
1687 : real(dp) :: ene,enemin,enemax,deltaene
1688 : real(dp) :: omega,omega_min,omega_max,domega
1689 : real(dp) :: diagerr
1690 : real(dp) :: chu_cond,chu_cth,femto
1691 0 : real(dp) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
1692 0 : real(dp) :: eigval(elph_ds%nbranch),eigval2(elph_ds%nbranch)
1693 0 : real(dp) :: imeigval(elph_ds%nbranch)
1694 : real(dp) :: tmp_wtkpq, tmp_wtkmq, tol_wtk
1695 : real(dp) :: yp1,ypn
1696 : !arrays
1697 0 : integer,allocatable :: FSfullpktofull(:,:),mqtofull(:)
1698 0 : integer,allocatable :: kpttokpt(:,:,:)
1699 : real(dp) :: cond_inv(3,3)
1700 0 : real(dp),allocatable :: fermie(:)
1701 0 : real(dp),allocatable :: tmp_eigenGS(:,:,:)
1702 0 : real(dp),allocatable :: tmp_gkk_qpt(:,:,:),tmp_gkk_rpt(:,:,:),tmp_gkk_kpt(:,:)
1703 0 : real(dp),allocatable :: tmp_gkk_kpt2(:,:,:), gkk_kpt(:,:,:)
1704 0 : real(dp),allocatable :: tau_k(:,:,:,:),inv_tau_k(:,:,:,:),tmp_tau_k(:,:,:,:)
1705 0 : real(dp),allocatable :: phfrq(:,:),pheigvec(:,:)
1706 0 : real(dp),allocatable :: displ(:,:,:,:)
1707 0 : real(dp),allocatable :: a2f_2d(:),a2f_2d2(:)
1708 0 : real(dp),allocatable :: tmp_wtk(:,:,:,:),tmp2_wtk(:),tmp_wtk1(:),tmp_wtk2(:)
1709 0 : real(dp),allocatable :: ene_pt(:),ene_ptfine(:),ff2(:)
1710 0 : real(dp),allocatable :: wtq(:,:,:),tmp_wtq(:,:,:),tmp2_wtq(:,:)
1711 0 : real(dp),allocatable :: dos_e(:,:)
1712 0 : real(dp),allocatable :: coskr1(:,:),sinkr1(:,:)
1713 0 : real(dp),allocatable :: coskr2(:,:),sinkr2(:,:)
1714 0 : real(dp),allocatable :: cond_e(:,:,:,:),cond(:,:,:,:),sbk(:,:,:,:),seebeck(:,:,:,:),cth(:,:,:,:)
1715 :
1716 : ! *************************************************************************
1717 :
1718 0 : write(std_out,*) 'get_tau_k : enter '
1719 :
1720 0 : nrpt = ifc%nrpt
1721 0 : natom = cryst%natom
1722 :
1723 0 : nsppol = elph_ds%nsppol
1724 0 : nbranch = elph_ds%nbranch
1725 0 : nband = elph_ds%ngkkband
1726 0 : nkpt = elph_ds%k_phon%nkpt
1727 0 : nqpt = elph_ds%nqpt_full
1728 0 : nkptirr = elph_ds%k_phon%nkptirr
1729 0 : new_nkptirr = elph_ds%k_phon%new_nkptirr
1730 0 : nqptirr = elph_ds%nqptirred
1731 0 : ntemper = elph_ds%ntemper
1732 0 : nene = 2*elph_ds%na2f-1 ! only need e_k +- omega_max range, take deltaene=delta_oemga
1733 :
1734 0 : chu_cond = 4.59988159904764*1.0d6 ! REPLACE WITH DATA FROM DEFS_BASIS
1735 0 : chu_cth = 1.078637439971599*1.0d4 ! REPLACE WITH DATA FROM DEFS_BASIS
1736 0 : femto = 1.0d-15
1737 :
1738 0 : tol_wtk = tol7/nkptirr/nband
1739 :
1740 0 : ABI_MALLOC(fermie ,(ntemper))
1741 0 : ABI_MALLOC(tmp_gkk_qpt ,(2,nbranch**2,nqpt))
1742 0 : ABI_MALLOC(tmp_gkk_rpt ,(2,nbranch**2,nrpt))
1743 0 : ABI_MALLOC(tmp_gkk_kpt ,(2,nbranch**2))
1744 0 : ABI_MALLOC(tmp_gkk_kpt2 ,(2,nbranch,nbranch))
1745 0 : ABI_MALLOC(gkk_kpt ,(2,nbranch,nbranch))
1746 0 : ABI_MALLOC(a2f_2d, (nene))
1747 0 : ABI_MALLOC(a2f_2d2, (nene))
1748 0 : ABI_MALLOC(inv_tau_k, (ntemper,nsppol,nkpt,nband))
1749 0 : ABI_MALLOC(tau_k, (ntemper,nsppol,nkpt,nband))
1750 0 : ABI_MALLOC(tmp_tau_k ,(ntemper,nsppol,new_nkptirr,nband))
1751 :
1752 0 : if (elph_ds%gkqwrite == 0) then
1753 0 : call wrtout(std_out,' get_tau_k : keeping gkq matrices in memory','COLL')
1754 0 : else if (elph_ds%gkqwrite == 1) then
1755 0 : fname=trim(elph_ds%elph_base_name) // '_GKKQ'
1756 0 : write (message,'(2a)')' get_tau_k : reading gkq matrices from file ',trim(fname)
1757 0 : call wrtout(std_out,message,'COLL')
1758 : else
1759 0 : write (message,'(a,i0)')' Wrong value for gkqwrite = ',elph_ds%gkqwrite
1760 0 : ABI_BUG(message)
1761 : end if
1762 :
1763 : !=========================================================
1764 : !Get equivalence between a kpt_phon pair and a qpt in qpt_full
1765 : !only works if the qpt grid is complete (identical to
1766 : !the kpt one, with a basic shift of (0,0,0)
1767 : !=========================================================
1768 :
1769 : !mapping of k + q onto k' for k and k' in full BZ
1770 : !for dense k grid
1771 0 : ABI_MALLOC(FSfullpktofull,(nkpt,nkpt))
1772 0 : ABI_MALLOC(mqtofull,(nkpt))
1773 :
1774 : !kpttokpt(itim,isym,iqpt) = kpoint index which transforms to ikpt under isym and with time reversal itim.
1775 0 : ABI_MALLOC(kpttokpt,(2,Cryst%nsym,nkpt))
1776 :
1777 0 : call wrtout(std_out,'get_tau_k: calling mkqptequiv to set up the FS kpoint set',"COLL")
1778 :
1779 0 : call mkqptequiv (FSfullpktofull,Cryst,elph_ds%k_phon%kpt,nkpt,nkpt,kpttokpt,elph_ds%k_phon%kpt,mqtofull)
1780 :
1781 : !=========================================================
1782 : !=========================================================
1783 :
1784 0 : omega_max = elph_ds%omega_max
1785 0 : omega_min = elph_ds%omega_min
1786 0 : domega = elph_ds%domega
1787 0 : enemax = maxval(eigenGS(elph_ds%maxFSband,:,:))
1788 0 : enemin = minval(eigenGS(elph_ds%minFSband,:,:))
1789 :
1790 0 : if (enemin < (elph_ds%fermie-0.2)) then
1791 0 : enemin = elph_ds%fermie-0.2
1792 : end if
1793 0 : if (enemax > (elph_ds%fermie+0.2)) then
1794 0 : enemax = elph_ds%fermie+0.2
1795 : end if
1796 :
1797 0 : nspline = elph_ds%ep_nspline
1798 0 : nene_all = INT((enemax-enemin+domega)/(nspline*domega)) + 1
1799 0 : deltaene = domega
1800 0 : write(std_out,*) 'E_min= ',enemin, 'E_max= ',enemax
1801 0 : write(std_out,*) 'Number of energy points= ',nene_all
1802 0 : write(std_out,'(a,I8)') 'scale factor for spline interpolation in RTA = ', elph_ds%ep_nspline
1803 0 : write(std_out,*) 'delta_ene before spline interpolation= ',deltaene*nspline
1804 0 : write(std_out,*) 'delta_ene after spline interpolation= ',deltaene
1805 0 : write(std_out,*) 'Omega_min= ',omega_min, 'Omega_max= ',omega_max
1806 0 : write(std_out,*) 'Number of phonon points= ',elph_ds%na2f
1807 0 : write(std_out,*) 'delta_omega= ',domega
1808 0 : write(std_out,*) 'number of bands= ', elph_ds%nband, nband
1809 :
1810 0 : ABI_MALLOC(tmp_wtk,(nband,nkpt,nsppol,nene_all))
1811 0 : ABI_MALLOC(tmp2_wtk,(nene_all))
1812 0 : ABI_MALLOC(ff2,(nene_all))
1813 0 : ABI_MALLOC(ene_pt,(nene_all))
1814 0 : ABI_MALLOC(ene_ptfine,(nene_all*nspline))
1815 0 : ABI_MALLOC(tmp_wtk1,(nene_all*nspline))
1816 0 : ABI_MALLOC(tmp_wtk2,(nene_all*nspline))
1817 0 : ABI_MALLOC(dos_e,(nsppol,nene_all))
1818 :
1819 : !Get energy points for spline interpolation
1820 0 : do iene = 1, nene_all
1821 0 : ene_pt(iene) = enemin + (iene-1)*nspline*deltaene
1822 : end do
1823 :
1824 0 : do iene = 1, nene_all*nspline
1825 0 : ene_ptfine(iene) = enemin + (iene-1)*deltaene
1826 : end do
1827 :
1828 0 : ABI_MALLOC(tmp_wtq,(elph_ds%nbranch, elph_ds%k_phon%nkpt, elph_ds%na2f+1))
1829 0 : ABI_MALLOC(wtq,(elph_ds%nbranch, elph_ds%k_phon%nkpt, elph_ds%na2f))
1830 0 : ABI_MALLOC(tmp2_wtq,(elph_ds%nbranch, elph_ds%na2f))
1831 :
1832 : !phonon
1833 0 : ABI_MALLOC(phfrq,(nbranch, nkptirr))
1834 0 : ABI_MALLOC(displ,(2, nbranch, nbranch, nkptirr))
1835 0 : ABI_MALLOC(pheigvec,(2*nbranch*nbranch, nkptirr))
1836 :
1837 0 : do iFSqpt = 1, nkptirr
1838 0 : call ifc%fourq(cryst,elph_ds%k_phon%kptirr(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt),out_eigvec=pheigvec(:,iFSqpt))
1839 : end do
1840 :
1841 0 : omega_min = omega_min - domega
1842 :
1843 : !bxu, obtain wtq for the q_fine, then condense to q_phon
1844 : call ep_ph_weights(phfrq,elph_ds%a2fsmear,omega_min,omega_max,elph_ds%na2f+1,Cryst%gprimd,elph_ds%kptrlatt, &
1845 0 : & elph_ds%nbranch,elph_ds%telphint,elph_ds%k_phon,tmp_wtq)
1846 0 : omega_min = omega_min + domega
1847 :
1848 0 : do iomega = 1, elph_ds%na2f
1849 0 : wtq(:,:,iomega) = tmp_wtq(:,:,iomega+1)
1850 : !write(1005,*) omega_min+(iomega-1)*domega, sum(tmp_wtq(:,:,iomega+1))/nkpt
1851 : end do
1852 0 : ABI_FREE(tmp_wtq)
1853 :
1854 : ! electron
1855 0 : tmp_wtk =zero
1856 0 : dos_e = zero
1857 : call ep_el_weights(elph_ds%ep_b_min, elph_ds%ep_b_max, eigenGS(elph_ds%minFSband:elph_ds%minFSband+nband-1,:,:), &
1858 : & elph_ds%elphsmear, &
1859 : & enemin, enemax, nene_all, Cryst%gprimd, elph_ds%k_phon%irredtoGS, elph_ds%kptrlatt, max_occ, &
1860 0 : & 1, nband, elph_ds%nFSband, nsppol, elph_ds%telphint, elph_ds%k_phon, tmp_wtk)
1861 : !& elph_ds%minFSband, elph_ds%nband, elph_ds%nFSband, nsppol, elph_ds%telphint, elph_ds%k_phon, tmp_wtk)
1862 :
1863 0 : do isppol = 1, nsppol
1864 0 : do iene = 1, nene_all
1865 0 : dos_e(isppol,iene) = sum(tmp_wtk(:,:,isppol,iene))/nkpt
1866 : end do
1867 : end do
1868 :
1869 0 : ABI_MALLOC(coskr1, (nqpt,nrpt))
1870 0 : ABI_MALLOC(sinkr1, (nqpt,nrpt))
1871 0 : call ftgam_init(ifc%gprim, nqpt, nrpt, elph_ds%k_phon%kpt, Ifc%rpt, coskr1, sinkr1)
1872 0 : ABI_MALLOC(coskr2, (nkptirr,nrpt))
1873 0 : ABI_MALLOC(sinkr2, (nkptirr,nrpt))
1874 0 : call ftgam_init(ifc%gprim, nkptirr, nrpt, elph_ds%k_phon%kpt, Ifc%rpt, coskr2, sinkr2)
1875 :
1876 : !get fermie for itemp
1877 0 : fermie = elph_ds%fermie
1878 0 : do itemp=1,ntemper ! runs over termperature in K
1879 0 : Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp)
1880 :
1881 0 : Bst%occopt = 3
1882 0 : Bst%tsmear = Temp*kb_HaK
1883 0 : call Bst%update_occ(-99.99_dp)
1884 0 : write(message,'(a,f12.6,a,E20.12)')'At T=',Temp,' Fermi level is:',Bst%fermie
1885 0 : call wrtout(std_out,message,'COLL')
1886 :
1887 0 : if (abs(elph_ds%fermie) < tol10) then
1888 0 : fermie(itemp) = Bst%fermie
1889 : end if
1890 : end do
1891 :
1892 0 : inv_tau_k = zero
1893 : !get a2f_2d = \sum_{q,nbranch,jband'} |gkk|^2*\delta(\epsilon_{k'j'}-\epsilon')*\delta(\omega_q-\omega)
1894 0 : do isppol=1,nsppol
1895 0 : write (std_out,*) '##############################################'
1896 0 : write (std_out,*) 'get_tau_k : Treating spin polarization ', isppol
1897 0 : write (std_out,*) '##############################################'
1898 :
1899 : ! do iFSkpt =1,nkpt
1900 0 : do ik_this_proc =1,elph_ds%k_phon%my_nkpt
1901 0 : iFSkpt = elph_ds%k_phon%my_ikpt(ik_this_proc)
1902 0 : write (std_out,*) 'get_tau_k : working on kpt # ', iFSkpt, '/', nkpt
1903 0 : do jband = 1, nband
1904 : ! write(*,*)'i am here 1 ', isppol,iFSkpt,jband
1905 0 : a2f_2d = zero
1906 0 : a2f_2d2 = zero
1907 :
1908 : !sum from here
1909 0 : nskip = 0
1910 0 : do jpband = 1, nband
1911 0 : jbeff = jpband+(jband-1)*nband
1912 :
1913 0 : if (elph_ds%gkqwrite == 0) then
1914 0 : tmp_gkk_qpt(:,:,:) = elph_ds%gkk_qpt(:,jbeff,:,ik_this_proc,isppol,:)
1915 0 : else if (elph_ds%gkqwrite == 1) then
1916 0 : irec = (ik_this_proc-1)*elph_ds%k_phon%my_nkpt + iqpt
1917 0 : if (iFSkpt == 1) then
1918 0 : write (std_out,*) ' get_tau_k read record ', irec
1919 : end if
1920 0 : read (elph_ds%unitgkq,REC=irec) tmp_gkk_qpt(:,:,iqpt_fullbz)
1921 : end if
1922 :
1923 : !FT to real space
1924 0 : call ftgam(Ifc%wghatm,tmp_gkk_qpt,tmp_gkk_rpt,natom,nqpt,nrpt,1,coskr1,sinkr1)
1925 :
1926 : !sum over irred q over k_phon, with corresponding weights
1927 0 : do iFSqpt = 1, nkptirr
1928 0 : iqpt_fullbz = elph_ds%k_phon%irredtoGS(iFSqpt)
1929 0 : ikpt_kpq = FSfullpktofull(iFSkpt,iqpt_fullbz)
1930 :
1931 0 : imqpt_fullbz = mqtofull(iqpt_fullbz)
1932 0 : ikpt_kmq = FSfullpktofull(iFSkpt,imqpt_fullbz)
1933 :
1934 : !Do FT from real-space gamma grid to 1 kpt in k_phon%new_kptirr
1935 0 : call ftgam(Ifc%wghatm,tmp_gkk_kpt,tmp_gkk_rpt,natom,1,nrpt,0,coskr2(iqpt_fullbz,:),sinkr2(iqpt_fullbz,:))
1936 : !tmp_gkk_kpt(:,:)=tmp_gkk_qpt(:,:,iFSqpt)
1937 :
1938 : !if ep_scalprod==0 we have to dot in the displacement vectors here
1939 0 : if (elph_ds%ep_scalprod==0) then
1940 :
1941 0 : call phdispl_cart2red(natom,Cryst%gprimd,displ(:,:,:,iFSqpt),displ_red)
1942 :
1943 0 : tmp_gkk_kpt2 = reshape (tmp_gkk_kpt(:,:), (/2,nbranch,nbranch/))
1944 0 : call gam_mult_displ(nbranch, displ_red, tmp_gkk_kpt2, gkk_kpt)
1945 :
1946 0 : do jbranch=1,nbranch
1947 0 : eigval(jbranch) = gkk_kpt(1, jbranch, jbranch)
1948 0 : imeigval(jbranch) = gkk_kpt(2, jbranch, jbranch)
1949 :
1950 0 : if (abs(imeigval(jbranch)) > tol10) then
1951 0 : write (message,'(a,i0,a,es16.8)')" real values branch = ",jbranch,' eigval = ',eigval(jbranch)
1952 0 : ABI_WARNING(message)
1953 0 : write (message,'(a,i0,a,es16.8)')" imaginary values branch = ",jbranch,' imeigval = ',imeigval(jbranch)
1954 0 : ABI_WARNING(message)
1955 : end if
1956 :
1957 : end do
1958 :
1959 : ! if ep_scalprod==1 we have to diagonalize the matrix we interpolated.
1960 0 : else if (elph_ds%ep_scalprod == 1) then
1961 :
1962 : ! MJV NOTE : gam_now is being recast as a (3*natom)**2 matrix here
1963 : call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, cone, tmp_gkk_kpt, 3*natom,&
1964 0 : & pheigvec(:,iFSqpt), 3*natom, czero, tmp_gkk_kpt2, 3*natom)
1965 :
1966 : call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, cone, pheigvec(:,iFSqpt), 3*natom,&
1967 0 : & tmp_gkk_kpt2, 3*natom, czero, gkk_kpt, 3*natom)
1968 :
1969 0 : diagerr = zero
1970 0 : do ibranch=1,nbranch
1971 0 : eigval(ibranch) = gkk_kpt(1,ibranch,ibranch)
1972 0 : do jbranch=1,ibranch-1
1973 0 : diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch))
1974 : end do
1975 0 : do jbranch=ibranch+1,nbranch
1976 0 : diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch))
1977 : end do
1978 : end do
1979 :
1980 0 : if (diagerr > tol12) then
1981 0 : write(message,'(a,es15.8)') 'get_tau_k: residual in diagonalization of gamma with phon eigenvectors: ', diagerr
1982 0 : ABI_WARNING(message)
1983 : end if
1984 :
1985 : else
1986 0 : write (message,'(a,i0)')' Wrong value for ep_scalprod = ',elph_ds%ep_scalprod
1987 0 : ABI_BUG(message)
1988 : end if ! end ep_scalprod if
1989 :
1990 : !For k'=k-q
1991 : !Do FT from real-space gamma grid to 1 kpt in k_phon%new_kptirr
1992 0 : call ftgam(Ifc%wghatm,tmp_gkk_kpt,tmp_gkk_rpt,natom,1,nrpt,0,coskr2(imqpt_fullbz,:),sinkr2(imqpt_fullbz,:))
1993 : !tmp_gkk_kpt(:,:)=tmp_gkk_qpt(:,:,iFSqpt)
1994 :
1995 : !if ep_scalprod==0 we have to dot in the displacement vectors here
1996 0 : if (elph_ds%ep_scalprod==0) then
1997 :
1998 0 : call phdispl_cart2red(natom,Cryst%gprimd,displ(:,:,:,iFSqpt),displ_red)
1999 :
2000 0 : tmp_gkk_kpt2 = reshape (tmp_gkk_kpt(:,:), (/2,nbranch,nbranch/))
2001 0 : call gam_mult_displ(nbranch, displ_red, tmp_gkk_kpt2, gkk_kpt)
2002 :
2003 0 : do jbranch=1,nbranch
2004 0 : eigval2(jbranch) = gkk_kpt(1, jbranch, jbranch)
2005 0 : imeigval(jbranch) = gkk_kpt(2, jbranch, jbranch)
2006 :
2007 0 : if (abs(imeigval(jbranch)) > tol10) then
2008 0 : write (message,'(a,i0,a,es16.8)')" real values branch = ",jbranch,' eigval = ',eigval2(jbranch)
2009 0 : ABI_WARNING(message)
2010 0 : write (message,'(a,i0,a,es16.8)')" imaginary values branch = ",jbranch,' imeigval = ',imeigval(jbranch)
2011 0 : ABI_WARNING(message)
2012 : end if
2013 :
2014 : end do
2015 :
2016 : ! if ep_scalprod==1 we have to diagonalize the matrix we interpolated.
2017 0 : else if (elph_ds%ep_scalprod == 1) then
2018 :
2019 : ! MJV NOTE : gam_now is being recast as a (3*natom)**2 matrix here
2020 : call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, cone, tmp_gkk_kpt, 3*natom,&
2021 0 : & pheigvec(:,iFSqpt), 3*natom, czero, tmp_gkk_kpt2, 3*natom)
2022 :
2023 : call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, cone, pheigvec(:,iFSqpt), 3*natom,&
2024 0 : & tmp_gkk_kpt2, 3*natom, czero, gkk_kpt, 3*natom)
2025 :
2026 0 : diagerr = zero
2027 0 : do ibranch=1,nbranch
2028 0 : eigval2(ibranch) = gkk_kpt(1,ibranch,ibranch)
2029 0 : do jbranch=1,ibranch-1
2030 0 : diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch))
2031 : end do
2032 0 : do jbranch=ibranch+1,nbranch
2033 0 : diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch))
2034 : end do
2035 : end do
2036 :
2037 0 : if (diagerr > tol12) then
2038 0 : write(message,'(a,es15.8)') 'get_tau_k: residual in diagonalization of gamma with phon eigenvectors: ', diagerr
2039 0 : ABI_WARNING(message)
2040 : end if
2041 :
2042 : else
2043 0 : write (message,'(a,i0)')' Wrong value for ep_scalprod = ',elph_ds%ep_scalprod
2044 0 : ABI_BUG(message)
2045 : end if ! end ep_scalprod if
2046 :
2047 0 : tmp2_wtk(:) = tmp_wtk(jpband,ikpt_kpq,isppol,:)
2048 0 : yp1 = (tmp2_wtk(2)-tmp2_wtk(1))/nspline/deltaene
2049 0 : ypn = (tmp2_wtk(nene_all)-tmp2_wtk(nene_all-1))/nspline/deltaene
2050 0 : call spline(ene_pt,tmp2_wtk,nene_all,yp1,ypn,ff2)
2051 0 : call splint(nene_all,ene_pt,tmp2_wtk,ff2,nene_all*nspline,ene_ptfine,tmp_wtk1)
2052 :
2053 0 : tmp2_wtk(:) = tmp_wtk(jpband,ikpt_kmq,isppol,:)
2054 0 : yp1 = (tmp2_wtk(2)-tmp2_wtk(1))/nspline/deltaene
2055 0 : ypn = (tmp2_wtk(nene_all)-tmp2_wtk(nene_all-1))/nspline/deltaene
2056 0 : call spline(ene_pt,tmp2_wtk,nene_all,yp1,ypn,ff2)
2057 0 : call splint(nene_all,ene_pt,tmp2_wtk,ff2,nene_all*nspline,ene_ptfine,tmp_wtk2)
2058 :
2059 0 : tmp2_wtq(:,:) = wtq(:,iFSqpt,:)
2060 0 : do iene=1,nene
2061 0 : e_k = eigenGS(elph_ds%minFSband+jband-1,iFSkpt,isppol)
2062 0 : ene = e_k - omega_max + (iene-1)*deltaene
2063 0 : if (ene<enemin .or. ene>enemax) cycle
2064 0 : iene_fine = NINT((ene-enemin+deltaene)/deltaene)
2065 0 : tmp_wtkpq = tmp_wtk1(iene_fine) * elph_ds%k_phon%wtkirr(iFSqpt)
2066 0 : tmp_wtkmq = tmp_wtk2(iene_fine) * elph_ds%k_phon%wtkirr(iFSqpt)
2067 :
2068 0 : if (tmp_wtkpq+tmp_wtkmq < tol_wtk ) then
2069 0 : nskip = nskip +1
2070 : cycle
2071 : end if
2072 :
2073 0 : do ibranch = 1, nbranch
2074 0 : if (abs(phfrq(ibranch,iFSqpt)) < tol7) cycle
2075 :
2076 0 : if (ene > e_k) then
2077 0 : omega = ene - e_k
2078 0 : if (abs(omega) < tol7 .or. abs(omega) > omega_max) cycle
2079 0 : iomega = NINT((omega-omega_min+domega)/domega)
2080 :
2081 : a2f_2d(iene) = a2f_2d(iene) +&
2082 : & eigval(ibranch)/phfrq(ibranch,iFSqpt)*&
2083 0 : & tmp_wtkpq * tmp2_wtq(ibranch,iomega)
2084 : end if
2085 :
2086 0 : if (ene < e_k) then
2087 0 : omega = e_k - ene
2088 0 : if (abs(omega) < tol7 .or. abs(omega) > omega_max) cycle
2089 0 : iomega = NINT((omega-omega_min+domega)/domega)
2090 :
2091 : a2f_2d2(iene) = a2f_2d2(iene) +&
2092 : & eigval(ibranch)/phfrq(ibranch,iFSqpt)*&
2093 0 : & tmp_wtkmq * tmp2_wtq(ibranch,iomega)
2094 : end if
2095 :
2096 : end do ! ibranch 3
2097 : end do ! nene 800
2098 : end do ! kptirr 216
2099 : end do ! j' band 3
2100 : ! print *, ' skipped ', nskip, ' energy points out of ', nene*nband*nkptirr
2101 :
2102 : ! get inv_tau_k
2103 0 : do itemp=1,ntemper ! runs over termperature in K
2104 0 : Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp)
2105 0 : do iene=1,nene
2106 0 : e_k = eigenGS(elph_ds%minFSband+jband-1,iFSkpt,isppol)
2107 0 : ene = e_k - omega_max + (iene-1)*deltaene
2108 0 : if (ene<enemin .or. ene>enemax) cycle
2109 :
2110 0 : xx=(ene-fermie(itemp))/(kb_HaK*Temp)
2111 0 : occ_e=1.0_dp/(exp(xx)+1.0_dp)
2112 0 : if (ene > e_k .and. (ene-e_k) .le. omega_max) then
2113 0 : omega = ene - e_k
2114 0 : if (abs(omega) < tol7) cycle
2115 0 : xx = omega/(kb_HaK*Temp)
2116 0 : occ_omega=1.0_dp/(exp(xx)-1.0_dp)
2117 :
2118 0 : therm_factor = occ_e + occ_omega
2119 :
2120 : inv_tau_k(itemp,isppol,iFSkpt,jband) = inv_tau_k(itemp,isppol,iFSkpt,jband) +&
2121 0 : a2f_2d(iene)*therm_factor*deltaene
2122 : end if
2123 0 : if (ene < e_k .and. (e_k-ene) .le. omega_max) then
2124 0 : omega = e_k - ene
2125 0 : if (abs(omega) < tol7) cycle
2126 0 : xx = omega/(kb_HaK*Temp)
2127 0 : occ_omega=1.0_dp/(exp(xx)-1.0_dp)
2128 :
2129 0 : therm_factor = 1 - occ_e + occ_omega
2130 :
2131 : inv_tau_k(itemp,isppol,iFSkpt,jband) = inv_tau_k(itemp,isppol,iFSkpt,jband) +&
2132 0 : a2f_2d2(iene)*therm_factor*deltaene
2133 : end if
2134 :
2135 : end do ! nene
2136 : end do ! Temp
2137 : ! write(*,*)'i am here 2 ', isppol,iFSkpt,jband
2138 : end do ! jband
2139 : end do ! kpt
2140 : end do ! nsppol
2141 :
2142 : !write (300+mpi_enreg%me,*) inv_tau_k
2143 0 : call xmpi_sum (inv_tau_k, xmpi_world, ierr)
2144 :
2145 0 : ABI_FREE(phfrq)
2146 0 : ABI_FREE(displ)
2147 0 : ABI_FREE(pheigvec)
2148 0 : ABI_FREE(tmp2_wtk)
2149 0 : ABI_FREE(ff2)
2150 0 : ABI_FREE(ene_pt)
2151 0 : ABI_FREE(ene_ptfine)
2152 0 : ABI_FREE(tmp_wtk1)
2153 0 : ABI_FREE(tmp_wtk2)
2154 0 : ABI_FREE(tmp2_wtq)
2155 0 : ABI_FREE(wtq)
2156 0 : ABI_FREE(coskr1)
2157 0 : ABI_FREE(sinkr1)
2158 0 : ABI_FREE(coskr2)
2159 0 : ABI_FREE(sinkr2)
2160 0 : ABI_FREE(kpttokpt)
2161 0 : ABI_FREE(FSfullpktofull)
2162 0 : ABI_FREE(mqtofull)
2163 0 : ABI_FREE(tmp_gkk_qpt)
2164 0 : ABI_FREE(tmp_gkk_rpt)
2165 0 : ABI_FREE(tmp_gkk_kpt)
2166 0 : ABI_FREE(tmp_gkk_kpt2)
2167 0 : ABI_FREE(gkk_kpt)
2168 0 : ABI_FREE(a2f_2d)
2169 0 : ABI_FREE(a2f_2d2)
2170 :
2171 : !output inv_tau_k and tau_k
2172 0 : fname = trim(elph_ds%elph_base_name) // '_INVTAUK'
2173 0 : if (open_file(fname,message,newunit=unit_invtau,status='unknown') /= 0) then
2174 0 : ABI_ERROR(message)
2175 : end if
2176 :
2177 : !print header to relaxation time file
2178 0 : write (unit_invtau,*) '# k-dep inverse of the relaxation time as a function of temperature.'
2179 0 : write (unit_invtau,*) '# '
2180 0 : write (unit_invtau,*) '# nkptirr= ', nkptirr, 'nband= ', nband
2181 0 : write (unit_invtau,*) '# number of temperatures= ', ntemper
2182 0 : write (unit_invtau,*) '# tau [femtosecond^-1] '
2183 :
2184 0 : fname = trim(elph_ds%elph_base_name) // '_TAUK'
2185 0 : if (open_file(fname,message,newunit=unit_tau,status='unknown') /= 0) then
2186 0 : ABI_ERROR(message)
2187 : end if
2188 :
2189 : !print header to relaxation time file
2190 0 : write (unit_tau,*) '# k-dep relaxation time as a function of temperature.'
2191 0 : write (unit_tau,*) '# '
2192 0 : write (unit_tau,*) '# nkptirr= ', nkptirr, 'nband= ', nband
2193 0 : write (unit_tau,*) '# number of temperatures= ', ntemper
2194 0 : write (unit_tau,*) '# tau [femtosecond] '
2195 :
2196 0 : tau_k = zero
2197 0 : do itemp=1,ntemper ! runs over termperature in K
2198 0 : Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp)
2199 0 : write(unit_invtau,'(a,f16.8)') '# Temperature = ', Temp
2200 0 : write(unit_tau,'(a,f16.8)') '# Temperature = ', Temp
2201 0 : do isppol=1,nsppol
2202 0 : write(unit_invtau,'(a,i6)') '# For isppol = ', isppol
2203 0 : write(unit_tau,'(a,i6)') '# For isppol = ', isppol
2204 0 : do iFSkpt = 1,nkpt
2205 : !FIXME: check when tau_k is too small, whether there should be a phonon
2206 : !scattering or not, and should tau_k be zero or not.
2207 0 : do jband = 1,nband
2208 0 : if (abs(inv_tau_k(itemp,isppol,iFSkpt,jband)) < tol9) then
2209 0 : inv_tau_k(itemp,isppol,iFSkpt,jband) = zero
2210 0 : tau_k(itemp,isppol,iFSkpt,jband) = zero
2211 : else
2212 : !no need to *nkpt due to wtkirr, as we need /nkpt for the sum
2213 : !no need to *two_pi due to the missing prefactor in gkk (see mka2f_tr_lova)
2214 0 : inv_tau_k(itemp,isppol,iFSkpt,jband) = inv_tau_k(itemp,isppol,iFSkpt,jband)*elph_ds%occ_factor
2215 0 : tau_k(itemp,isppol,iFSkpt,jband) = one/inv_tau_k(itemp,isppol,iFSkpt,jband)
2216 : end if
2217 : end do ! nband
2218 0 : write(unit_invtau,'(a,i8,a,3f12.6)') '# kpt# ', iFSkpt, ' kpt=', elph_ds%k_phon%kptirr(:,iFSkpt)
2219 0 : write(unit_invtau,'(100D16.8)') (inv_tau_k(itemp,isppol,iFSkpt,iband)*femto/Time_Sec,iband=1,nband)
2220 0 : write(unit_tau,'(a,i8,a,3f12.6)') '# kpt# ', iFSkpt, ' kpt=', elph_ds%k_phon%kptirr(:,iFSkpt)
2221 0 : write(unit_tau,'(100D16.8)') (tau_k(itemp,isppol,iFSkpt,iband)*Time_Sec/femto,iband=1,nband)
2222 : end do ! nkptirr
2223 0 : write(unit_invtau,*) ' '
2224 0 : write(unit_tau,*) ' '
2225 : end do ! nsppol
2226 0 : write(unit_invtau,*) ' '
2227 0 : write(unit_invtau,*) ' '
2228 0 : write(unit_tau,*) ' '
2229 0 : write(unit_tau,*) ' '
2230 : end do ! ntemper
2231 :
2232 : ! Only use the irred k for eigenGS and tau_k
2233 0 : ABI_MALLOC(tmp_eigenGS,(elph_ds%nband,elph_ds%k_phon%new_nkptirr,elph_ds%nsppol))
2234 :
2235 0 : do ikpt_irr = 1, new_nkptirr
2236 0 : tmp_eigenGS(:,ikpt_irr,:) = eigenGS(:,elph_ds%k_phon%new_irredtoGS(ikpt_irr),:)
2237 0 : tmp_tau_k(:,:,ikpt_irr,:) = tau_k(:,:,elph_ds%k_phon%new_irredtoGS(ikpt_irr),:)*Time_Sec
2238 : end do
2239 :
2240 : !BoltzTraP output files in SIESTA format
2241 0 : if (elph_ds%prtbltztrp == 1) then
2242 : !Prevent use in case occopt = 9
2243 0 : if (Bst%occopt==9) then
2244 0 : ABI_ERROR("Boltztrap outputting not possible with occopt = 9 at the moment")
2245 : end if
2246 : call ebands_prtbltztrp_tau_out (tmp_eigenGS(elph_ds%minFSband:elph_ds%maxFSband,:,:),&
2247 : & elph_ds%tempermin,elph_ds%temperinc,ntemper,fermie, &
2248 : & elph_ds%elph_base_name,elph_ds%k_phon%new_kptirr,nband,elph_ds%nelect,new_nkptirr, &
2249 0 : & elph_ds%nspinor,nsppol,Cryst%nsym,Cryst%rprimd,Cryst%symrel,tmp_tau_k)
2250 : end if !prtbltztrp
2251 0 : ABI_FREE(tmp_eigenGS)
2252 0 : ABI_FREE(tmp_tau_k)
2253 :
2254 : !Get the energy dependence of tau.
2255 : !Eq. (6) in Restrepo et al. Appl. Phys. Lett. 94, 212103 (2009) [[cite:Restrepo2009]]
2256 :
2257 0 : fname = trim(elph_ds%elph_base_name) // '_TAUE'
2258 0 : if (open_file(fname,message,newunit=unit_taue,status='unknown') /= 0) then
2259 0 : ABI_ERROR(message)
2260 : end if
2261 :
2262 : !print header to relaxation time file
2263 0 : write (unit_taue,*) '# Energy-dep relaxation time as a function of temperature.'
2264 0 : write (unit_taue,*) '# '
2265 0 : write (unit_taue,*) '# number of temperatures= ', ntemper
2266 0 : write (unit_taue,*) '# ene[Ha] tau [femtosecond] DOS[au] '
2267 :
2268 0 : fname = trim(elph_ds%elph_base_name) // '_MFP'
2269 0 : if (open_file(fname,message,newunit=unit_mfp,status='unknown') /= 0) then
2270 0 : ABI_ERROR(message)
2271 : end if
2272 :
2273 0 : write (unit_mfp,*) '# Energy-dep mean free path as a function of temperature.'
2274 0 : write (unit_mfp,*) '# '
2275 0 : write (unit_mfp,*) '# number of temperatures= ', ntemper
2276 0 : write (unit_mfp,*) '# ene[Ha] mfp [femtometer] '
2277 :
2278 0 : do itemp=1,ntemper ! runs over termperature in K
2279 0 : Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp)
2280 0 : write(unit_taue,'(a,f16.8)') '# Temperature = ', Temp
2281 0 : do isppol = 1, nsppol
2282 0 : write(unit_taue,*) '# Tau_e for isppol = ',isppol
2283 0 : do iene = 1, nene_all
2284 0 : rate_e = zero
2285 0 : do iFSkpt = 1, nkpt
2286 0 : do jband = 1, nband
2287 : rate_e = rate_e + inv_tau_k(itemp,isppol,iFSkpt,jband)* &
2288 0 : & tmp_wtk(jband,iFSkpt,isppol,iene)
2289 : end do ! jband
2290 : end do ! kpt
2291 0 : if (dabs(dos_e(isppol,iene)) < tol7) then
2292 : rate_e = zero
2293 : else
2294 0 : rate_e = rate_e/nkpt/dos_e(isppol,iene)
2295 : end if
2296 0 : write(unit_taue,"(3D16.8)") enemin+(iene-1)*deltaene*nspline, rate_e*femto/Time_Sec, dos_e(isppol,iene)
2297 : end do ! number of energies
2298 0 : write(unit_taue,*) ' '
2299 : end do ! nsppol
2300 0 : write(unit_taue,*) ' '
2301 : end do ! ntemperature
2302 :
2303 : ! calculate and output mean free path
2304 0 : do itemp=1,ntemper ! runs over termperature in K
2305 0 : Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp)
2306 0 : write(unit_mfp,'(a,f16.8)') '# Temperature = ', Temp
2307 0 : do isppol = 1, nsppol
2308 0 : do icomp = 1, 3
2309 0 : write(unit_mfp,*) '# Mean free path for isppol, icomp= ',isppol,icomp
2310 0 : do iene = 1, nene_all
2311 0 : mfp_e = zero
2312 0 : do iFSkpt = 1, nkptirr
2313 0 : do jband = 1, nband
2314 : mfp_e = mfp_e + tau_k(itemp,isppol,iFSkpt,jband)* &
2315 : & elph_tr_ds%el_veloc(iFSkpt,elph_ds%minFSband+jband-1,icomp,isppol)* &
2316 0 : & tmp_wtk(jband,iFSkpt,isppol,iene)
2317 : !& elph_ds%k_phon%new_wtkirr(iFSqpt)
2318 : end do ! jband
2319 : end do ! kpt
2320 0 : if (dabs(dos_e(isppol,iene)) < tol7) then
2321 : mfp_e = zero
2322 : else
2323 0 : mfp_e = mfp_e/nkptirr/dos_e(isppol,iene)
2324 : end if
2325 0 : write(unit_mfp,"(2D16.8)") enemin+(iene-1)*deltaene*nspline, mfp_e*Bohr_Ang*1.0d-10/femto
2326 : end do ! number of energies
2327 0 : write(unit_mfp,*) ' '
2328 : end do ! icomp
2329 0 : write(unit_mfp,*) ' '
2330 : end do ! nsppol
2331 0 : write(unit_mfp,*) ' '
2332 : end do ! ntemperature
2333 :
2334 0 : ABI_MALLOC(cond_e ,(ntemper,nsppol,nene_all,9))
2335 :
2336 : !get cond_e
2337 0 : cond_e = zero
2338 0 : do itemp=1,ntemper ! runs over termperature in K
2339 0 : do isppol = 1, nsppol
2340 0 : do iene = 1, nene_all
2341 : ! do iFSkpt =1,nkpt
2342 0 : do ik_this_proc =1,elph_ds%k_phon%my_nkpt
2343 0 : iFSkpt = elph_ds%k_phon%my_ikpt(ik_this_proc)
2344 0 : do jband = 1, nband
2345 0 : do icomp = 1, 3
2346 0 : do jcomp = 1, 3
2347 0 : itensor = (icomp-1)*3+jcomp
2348 : cond_e(itemp,isppol,iene,itensor) = cond_e(itemp,isppol,iene,itensor) + &
2349 : & tau_k(itemp,isppol,iFSkpt,jband)* &
2350 : & elph_tr_ds%el_veloc(iFSkpt,elph_ds%minFSband+jband-1,icomp,isppol)* &
2351 : & elph_tr_ds%el_veloc(iFSkpt,elph_ds%minFSband+jband-1,jcomp,isppol)* &
2352 0 : & tmp_wtk(jband,iFSkpt,isppol,iene)
2353 : end do
2354 : end do
2355 : end do ! jband
2356 : end do ! kpt
2357 : end do ! number of energies
2358 : end do ! nsppol
2359 : end do ! ntemperature
2360 :
2361 : ! MG FIXME: Why xmpi_world, besides only master should perform IO in the section below.
2362 0 : call xmpi_sum (cond_e, xmpi_world, ierr)
2363 :
2364 0 : cond_e = cond_e/nkpt
2365 :
2366 : !get transport coefficients
2367 :
2368 0 : fname = trim(elph_ds%elph_base_name) // '_COND'
2369 0 : if (open_file(fname,message,newunit=unit_cond,status='unknown') /= 0) then
2370 0 : ABI_ERROR(message)
2371 : end if
2372 :
2373 : !print header to conductivity file
2374 0 : write (unit_cond,*) '# Conductivity as a function of temperature.'
2375 0 : write (unit_cond,*) '# the formalism is isotropic, so non-cubic crystals may be wrong'
2376 0 : write (unit_cond,*) '# '
2377 0 : write (unit_cond,*) '# Columns are: '
2378 0 : write (unit_cond,*) '# temperature[K] cond[au] cond [SI] '
2379 0 : write (unit_cond,*) '# '
2380 :
2381 0 : fname = trim(elph_ds%elph_base_name) // '_CTH'
2382 0 : if (open_file(fname,message,newunit=unit_therm,status='unknown') /= 0) then
2383 0 : ABI_ERROR(message)
2384 : end if
2385 :
2386 : !print header to thermal conductivity file
2387 0 : write (unit_therm,'(a)') '# Thermal conductivity as a function of temperature.'
2388 0 : write (unit_therm,'(a)') '# the formalism is isotropic, so non-cubic crystals may be wrong'
2389 0 : write (unit_therm,'(a)') '# '
2390 0 : write (unit_therm,'(a)') '# Columns are: '
2391 0 : write (unit_therm,'(a)') '# temperature[K] thermal cond [au] thermal cond [SI]'
2392 0 : write (unit_therm,'(a)') '# '
2393 :
2394 0 : fname = trim(elph_ds%elph_base_name) // '_SBK'
2395 0 : if (open_file(fname,message,newunit=unit_sbk,status='unknown') /=0) then
2396 0 : ABI_ERROR(message)
2397 : end if
2398 :
2399 : !print header to relaxation time file
2400 0 : write (unit_sbk,*) '# Seebeck Coefficint as a function of temperature.'
2401 0 : write (unit_sbk,*) '# the formalism is isotropic, so non-cubic crystals may be wrong'
2402 0 : write (unit_sbk,*) '# '
2403 0 : write (unit_sbk,*) '# Columns are: '
2404 0 : write (unit_sbk,*) '# temperature[K] S [au] S [SI] '
2405 0 : write (unit_sbk,*) '# '
2406 :
2407 0 : ABI_MALLOC(cond ,(ntemper,nsppol,3,3))
2408 0 : ABI_MALLOC(cth ,(ntemper,nsppol,3,3))
2409 0 : ABI_MALLOC(sbk ,(ntemper,nsppol,3,3))
2410 0 : ABI_MALLOC(seebeck ,(ntemper,nsppol,3,3))
2411 :
2412 0 : cond = zero
2413 0 : cth = zero
2414 0 : sbk = zero
2415 0 : seebeck = zero
2416 0 : do isppol=1,nsppol
2417 0 : do icomp=1, 3
2418 0 : do jcomp=1, 3
2419 0 : itensor=(icomp-1)*3+jcomp
2420 0 : do itemp=1,ntemper
2421 0 : Temp=elph_ds%tempermin + elph_ds%temperinc*dble(itemp)
2422 0 : do iene = 1, nene_all
2423 0 : factor = (enemin+(iene-1)*deltaene*nspline - fermie(itemp))/(kb_HaK*Temp)
2424 0 : if (factor < -40.0d0) then
2425 : dfermide = zero
2426 0 : else if (factor > 40.0d0) then
2427 : dfermide = zero
2428 : else
2429 0 : dfermide = EXP(factor)/(kb_HaK*Temp*(EXP(factor)+one)**2.0d0)
2430 : end if
2431 : cond(itemp,isppol,icomp,jcomp) = cond(itemp,isppol,icomp,jcomp) + &
2432 0 : & cond_e(itemp,isppol,iene,itensor)*dfermide*deltaene*nspline
2433 : cth(itemp,isppol,icomp,jcomp) = cth(itemp,isppol,icomp,jcomp) + cond_e(itemp,isppol,iene,itensor)* &
2434 0 : & (enemin+(iene-1)*deltaene*nspline - fermie(itemp))**2.0d0*dfermide*deltaene*nspline
2435 : sbk(itemp,isppol,icomp,jcomp) = sbk(itemp,isppol,icomp,jcomp) + cond_e(itemp,isppol,iene,itensor)* &
2436 0 : & (enemin+(iene-1)*deltaene*nspline - fermie(itemp))*dfermide*deltaene*nspline
2437 : end do
2438 : end do ! temperature
2439 : end do ! jcomp
2440 : end do ! icomp
2441 : end do !end isppol
2442 :
2443 0 : do isppol=1,nsppol
2444 0 : do itemp=1,ntemper
2445 0 : cond_inv(:,:)=cond(itemp,isppol,:,:)
2446 0 : call matrginv(cond_inv,3,3)
2447 : call DGEMM('N','N',3,3,3,one,sbk(itemp,isppol,:,:),3,cond_inv,&
2448 0 : & 3,zero,seebeck(itemp,isppol,:,:),3)
2449 : end do
2450 : end do
2451 :
2452 0 : do isppol=1,nsppol
2453 0 : do icomp=1, 3
2454 0 : do jcomp=1, 3
2455 0 : itensor=(icomp-1)*3+jcomp
2456 0 : write(unit_cond,*) '# Conductivity for isppol, itrten= ',isppol,itensor
2457 0 : write(unit_therm,*) '# Thermal conductivity for isppol, itrten= ',isppol,itensor
2458 0 : write(unit_sbk,*) '# Seebeck coefficient for isppol, itrten= ',isppol,itensor
2459 0 : do itemp=1,ntemper
2460 0 : Temp=elph_ds%tempermin + elph_ds%temperinc*dble(itemp)
2461 :
2462 0 : seebeck(itemp,isppol,icomp,jcomp) = -1.0d0*seebeck(itemp,isppol,icomp,jcomp)/(kb_HaK*Temp)
2463 0 : cond(itemp,isppol,icomp,jcomp) = cond(itemp,isppol,icomp,jcomp)/cryst%ucvol
2464 0 : cth(itemp,isppol,icomp,jcomp) = cth(itemp,isppol,icomp,jcomp)/(kb_HaK*Temp)/cryst%ucvol
2465 0 : write(unit_cond,'(3D20.10)')Temp,cond(itemp,isppol,icomp,jcomp),cond(itemp,isppol,icomp,jcomp)*chu_cond
2466 0 : write(unit_therm,'(3D20.10)')Temp,cth(itemp,isppol,icomp,jcomp),cth(itemp,isppol,icomp,jcomp)*chu_cth
2467 0 : write(unit_sbk,'(3D20.10)')Temp,seebeck(itemp,isppol,icomp,jcomp),seebeck(itemp,isppol,icomp,jcomp)*kb_HaK
2468 : end do ! temperature
2469 0 : write(unit_cond,*)
2470 0 : write(unit_therm,*)
2471 0 : write(unit_sbk,*)
2472 : end do ! jcomp
2473 : end do ! icomp
2474 : end do !end isppol
2475 :
2476 :
2477 0 : ABI_FREE(inv_tau_k)
2478 0 : ABI_FREE(tau_k)
2479 0 : ABI_FREE(tmp_wtk)
2480 0 : ABI_FREE(dos_e)
2481 0 : ABI_FREE(cond_e)
2482 0 : ABI_FREE(fermie)
2483 0 : ABI_FREE(cond)
2484 0 : ABI_FREE(sbk)
2485 0 : ABI_FREE(cth)
2486 0 : ABI_FREE(seebeck)
2487 :
2488 0 : close (unit=unit_tau)
2489 0 : close (unit=unit_taue)
2490 0 : close (unit=unit_mfp)
2491 0 : close (unit=unit_invtau)
2492 0 : close (unit=unit_cond)
2493 0 : close (unit=unit_therm)
2494 0 : close (unit=unit_sbk)
2495 :
2496 0 : end subroutine get_tau_k
2497 : !!***
2498 :
2499 : end module m_a2ftr
2500 : !!***
|
