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