TABLE OF CONTENTS
ABINIT/mka2f_tr [ Functions ]
NAME
mka2f_tr
FUNCTION
calculates the FS averaged Transport alpha^2F_tr function calculates and outputs the associated electrical conductivity, relaxation time, and Seebeck coefficient and thermal conductivities for the first task : copied from mka2F
COPYRIGHT
Copyright (C) 2004-2018 ABINIT group (JPC, MJV, BXU) This file is distributed under the terms of the GNU General Public Licence, see ~abinit/COPYINGS= or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
crystal<crystal_t>=data type gathering info on the crystalline structure. Ifc<ifc_type>=Object containing the interatomic force constants. elph_ds elph_ds%gkk2 = gkk2 matrix elements on full FS grid for each phonon mode elph_ds%nbranch = number of phonon branches = 3*natom elph_ds%nFSband = number of bands included in the FS integration elph_ds%k_phon%nkpt = number of kpts included in the FS integration elph_ds%k_phon%wtk = integration weights on the FS delph_ds%n0 = DOS at the Fermi level calculated from the k_phon integration weights elph_ds%k_phon%kpt = coordinates of all FS kpoints mustar = coulomb pseudopotential parameter eventually for 2 spin channels natom = number of atoms ntemper = number of temperature points to calculate, from tempermin to tempermin+ntemper*temperinc tempermin = minimum temperature at which resistivity etc are calculated (in K) temperinc = interval for temperature grid on which resistivity etc are calculated (in K) elph_tr_ds%dos_n0 = DOS at the Fermi level calculated from the k_phon integration weights, but has a temperature dependence elph_tr_ds%dos_n = DOS at varied energy levels around Fermi level elph_tr_ds%veloc_sq0 = Fermi velocity square with T dependence
OUTPUT
elph_ds
PARENTS
elphon
CHILDREN
d2c_wtq,dgemm,ep_ph_weights,ftgam,ftgam_init,gam_mult_displ,ifc_fourq matrginv,simpson_int,wrtout,xmpi_sum,zgemm
NOTES
copied from ftiaf9.f
SOURCE
57 #if defined HAVE_CONFIG_H 58 #include "config.h" 59 #endif 60 61 #include "abi_common.h" 62 63 64 subroutine mka2f_tr(crystal,ifc,elph_ds,ntemper,tempermin,temperinc,pair2red,elph_tr_ds) 65 66 use defs_basis 67 use defs_elphon 68 use m_errors 69 use m_profiling_abi 70 use m_xmpi 71 72 use m_io_tools, only : open_file 73 use m_numeric_tools, only : simpson_int 74 use m_abilasi, only : matrginv 75 use m_crystal, only : crystal_t 76 use m_ifc, only : ifc_type, ifc_fourq 77 use m_dynmat, only : ftgam_init, ftgam 78 79 !This section has been created automatically by the script Abilint (TD). 80 !Do not modify the following lines by hand. 81 #undef ABI_FUNC 82 #define ABI_FUNC 'mka2f_tr' 83 use interfaces_14_hidewrite 84 use interfaces_77_ddb, except_this_one => mka2f_tr 85 !End of the abilint section 86 87 implicit none 88 89 !Arguments ------------------------------------ 90 !scalars 91 integer,intent(in) :: ntemper 92 real(dp),intent(in) :: tempermin,temperinc 93 type(ifc_type),intent(in) :: ifc 94 type(crystal_t),intent(in) :: crystal 95 type(elph_tr_type),intent(inout) :: elph_tr_ds 96 type(elph_type),intent(inout) :: elph_ds 97 !arrays 98 integer,intent(in) :: pair2red(elph_ds%nenergy,elph_ds%nenergy) 99 100 !Local variables ------------------------- 101 !x =w/(2kbT) 102 !scalars 103 integer :: iFSqpt,ibranch,iomega,isppol,jbranch,nerr 104 integer :: unit_a2f_tr,natom,ii,jj 105 integer :: idir, iatom, k1, kdir 106 integer :: unit_lor,unit_rho,unit_tau,unit_sbk, unit_therm 107 integer :: itemp, tmp_nenergy 108 integer :: itrtensor, icomp, jcomp!, kcomp 109 integer :: ie, ie_1, ie2, ie_2, ie1, ie_tmp, ssp, s1(4), s2(4) 110 integer :: ie2_left, ie2_right 111 integer :: ik_this_proc, ierr,nrpt 112 logical,parameter :: debug=.False. 113 real(dp) :: Temp,chgu,chtu,chsu,chwu,diagerr,ucvol 114 real(dp) :: a2fprefactor, gtemp 115 real(dp) :: lambda_tr,lor0,lorentz,maxerr,omega 116 real(dp) :: rho,tau,wtherm,xtr 117 real(dp) :: lambda_tr_trace 118 real(dp) :: domega, omega_min, omega_max 119 real(dp) :: gaussval, gaussprefactor, gaussfactor, gaussmaxarg, xx 120 real(dp) :: qnorm2, tmp_fermie 121 real(dp) :: e1, e2, diff, xe 122 real(dp) :: occ_omega, occ_e1, occ_e2 123 real(dp) :: nv1, nv2, sigma1, sigma2 124 real(dp) :: dos_n_e2, veloc_sq_icomp, veloc_sq_jcomp 125 real(dp) :: tointegq00_1, tointegq00_2, tointegq01_1, tointegq01_2,tointegq11_1, tointegq11_2 126 real(dp) :: j00, j01, j11 127 real(dp) :: tointegq00,tointegq01,tointegq11 128 real(dp) :: pref_s, pref_w, tmp_veloc_sq0, tmp_veloc_sq1, tmp_veloc_sq2 129 character(len=500) :: message 130 character(len=fnlen) :: fname 131 !arrays 132 real(dp),parameter :: c0(2)=(/0.d0,0.d0/),c1(2)=(/1.d0,0.d0/) 133 real(dp) :: gprimd(3,3) 134 real(dp) :: eigval(elph_ds%nbranch) 135 real(dp) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch) 136 real(dp) :: gam_now(2,elph_ds%nbranch*elph_ds%nbranch) 137 real(dp) :: tmpa2f(elph_ds%na2f) 138 real(dp) :: tmpgam1(2,elph_ds%nbranch,elph_ds%nbranch) 139 real(dp) :: tmpgam2(2,elph_ds%nbranch,elph_ds%nbranch) 140 real(dp) :: q11_inv(3,3) 141 !real(dp) :: fullq(6,6) 142 real(dp),allocatable :: phfrq(:,:) 143 real(dp),allocatable :: tmp_a2f_1d_tr(:,:,:,:,:) 144 real(dp),allocatable :: displ(:,:,:,:) 145 real(dp),allocatable :: pheigvec(:,:) 146 real(dp),allocatable :: tmp_wtq(:,:,:) 147 real(dp),allocatable :: integrho(:), tointegrho(:) 148 real(dp),allocatable :: integrand_q00(:),integrand_q01(:),integrand_q11(:) 149 real(dp),allocatable :: q00(:,:,:,:), q01(:,:,:,:),q11(:,:,:,:) 150 real(dp),allocatable :: seebeck(:,:,:,:)!, rho_nm(:,:,:,:) 151 real(dp),allocatable :: rho_T(:) 152 real(dp),allocatable :: coskr(:,:), sinkr(:,:) 153 154 ! ********************************************************************* 155 !calculate a2f_tr for frequencies between 0 and omega_max 156 157 158 write(std_out,*) 'mka2f_tr : enter ' 159 160 ucvol = crystal%ucvol 161 natom = crystal%natom 162 gprimd = crystal%gprimd 163 164 ! number of real-space points for FT interpolation 165 nrpt = ifc%nrpt 166 ! 167 !MG: the step should be calculated locally using nomega and the extrema of the spectrum. 168 !One should not rely on previous calls for the setup of elph_ds%domega 169 !I will remove elph_ds%domega since mka2f.F90 will become a method of gamma_t 170 domega =elph_ds%domega 171 omega_min = elph_ds%omega_min 172 omega_max = elph_ds%omega_max 173 174 if (elph_ds%ep_lova .eq. 1) then 175 tmp_nenergy = 1 176 else if (elph_ds%ep_lova .eq. 0) then 177 tmp_nenergy = elph_ds%nenergy 178 end if 179 180 !! defaults for number of temperature steps and max T (all in Kelvin...) 181 !ntemper=1000 182 !tempermin=zero 183 !temperinc=one 184 185 ABI_ALLOCATE(rho_T,(ntemper)) 186 187 188 gaussprefactor = sqrt(piinv) / elph_ds%a2fsmear 189 gaussfactor = one / elph_ds%a2fsmear 190 gaussmaxarg = sqrt(-log(1.d-90)) 191 !lor0=(pi*kb_HaK)**2/3. 192 lor0=pi**2/3.0_dp 193 194 !maximum value of frequency (a grid has to be chosen for the representation of alpha^2 F) 195 !WARNING! supposes this value has been set in mkelph_linwid. 196 197 !ENDMG 198 199 maxerr=0. 200 nerr=0 201 202 ABI_ALLOCATE(tmp_wtq,(elph_ds%nbranch, elph_ds%k_fine%nkpt, elph_ds%na2f+1)) 203 ABI_ALLOCATE(elph_ds%k_fine%wtq,(elph_ds%nbranch, elph_ds%k_fine%nkpt, elph_ds%na2f)) 204 ABI_ALLOCATE(elph_ds%k_phon%wtq,(elph_ds%nbranch, elph_ds%k_phon%nkpt, elph_ds%na2f)) 205 206 ABI_ALLOCATE(phfrq,(elph_ds%nbranch, elph_ds%k_fine%nkpt)) 207 ABI_ALLOCATE(displ,(2, elph_ds%nbranch, elph_ds%nbranch, elph_ds%k_fine%nkpt)) 208 ABI_ALLOCATE(pheigvec,(2*elph_ds%nbranch*elph_ds%nbranch, elph_ds%k_fine%nkpt)) 209 210 do iFSqpt=1,elph_ds%k_fine%nkpt 211 call ifc_fourq(ifc,crystal,elph_ds%k_fine%kpt(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt),out_eigvec=pheigvec(:,iFSqpt)) 212 end do 213 omega_min = omega_min - domega 214 215 !bxu, obtain wtq for the q_fine, then condense to q_phon 216 call ep_ph_weights(phfrq,elph_ds%a2fsmear,omega_min,omega_max,elph_ds%na2f+1,gprimd,elph_ds%kptrlatt_fine, & 217 & elph_ds%nbranch,elph_ds%telphint,elph_ds%k_fine,tmp_wtq) 218 !call ep_ph_weights(phfrq,elph_ds%a2fsmear,omega_min,omega_max,elph_ds%na2f+1,gprimd,elph_ds%kptrlatt_fine, & 219 !& elph_ds%nbranch,1,elph_ds%k_fine,tmp_wtq) 220 omega_min = omega_min + domega 221 222 do iomega = 1, elph_ds%na2f 223 elph_ds%k_fine%wtq(:,:,iomega) = tmp_wtq(:,:,iomega+1) 224 end do 225 ABI_DEALLOCATE(tmp_wtq) 226 227 if (elph_ds%use_k_fine == 1) then 228 call d2c_wtq(elph_ds) 229 end if 230 231 ABI_DEALLOCATE(phfrq) 232 ABI_DEALLOCATE(displ) 233 ABI_DEALLOCATE(pheigvec) 234 235 !reduce the dimension from fine to phon for phfrq and pheigvec 236 ABI_ALLOCATE(phfrq,(elph_ds%nbranch, elph_ds%k_phon%nkpt)) 237 ABI_ALLOCATE(displ,(2, elph_ds%nbranch, elph_ds%nbranch, elph_ds%k_phon%nkpt)) 238 ABI_ALLOCATE(pheigvec,(2*elph_ds%nbranch*elph_ds%nbranch, elph_ds%k_phon%nkpt)) 239 240 do iFSqpt=1,elph_ds%k_phon%nkpt 241 call ifc_fourq(ifc,crystal,elph_ds%k_phon%kpt(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt),out_eigvec=pheigvec(:,iFSqpt)) 242 end do 243 244 ABI_ALLOCATE(elph_tr_ds%a2f_1d_tr,(elph_ds%na2f,9,elph_ds%nsppol,4,tmp_nenergy**2,ntemper)) 245 ABI_ALLOCATE(tmp_a2f_1d_tr,(elph_ds%na2f,9,elph_ds%nsppol,4,tmp_nenergy**2)) 246 247 ! prepare phase factors 248 ABI_ALLOCATE(coskr, (elph_ds%k_phon%nkpt, nrpt)) 249 ABI_ALLOCATE(sinkr, (elph_ds%k_phon%nkpt, nrpt)) 250 call ftgam_init(ifc%gprim, elph_ds%k_phon%nkpt, nrpt, elph_ds%k_phon%kpt, ifc%rpt, coskr, sinkr) 251 252 elph_tr_ds%a2f_1d_tr = zero 253 tmp_a2f_1d_tr = zero 254 255 256 do ie = 1, elph_ds%n_pair 257 do ssp = 1,4 258 do isppol = 1, elph_ds%nsppol 259 260 ! loop over qpoint in full kpt grid (presumably dense) 261 do ik_this_proc =1,elph_ds%k_phon%my_nkpt 262 iFSqpt = elph_ds%k_phon%my_ikpt(ik_this_proc) 263 264 qnorm2 = sum(elph_ds%k_phon%kpt(:,iFSqpt)**2) 265 ! if (flag_to_exclude_soft_modes = .false.) qnorm2 = zero 266 do itrtensor=1,9 267 268 if (elph_ds%ep_int_gkk == 1) then 269 gam_now(:,:) = elph_tr_ds%gamma_qpt_tr(:,itrtensor,:,isppol,iFSqpt) 270 else 271 ! Do FT from real-space gamma grid to 1 qpt. 272 call ftgam(ifc%wghatm,gam_now,elph_tr_ds%gamma_rpt_tr(:,itrtensor,:,isppol,:,ssp,ie),natom,1,nrpt,0,& 273 & coskr(iFSqpt,:), sinkr(iFSqpt,:)) 274 end if 275 276 ! Diagonalize gamma matrix at this qpoint (complex matrix). 277 278 ! if ep_scalprod==0 we have to dot in the displacement vectors here 279 if (elph_ds%ep_scalprod==0) then 280 281 displ_red(:,:,:) = zero 282 do jbranch=1,elph_ds%nbranch 283 do iatom=1,natom 284 do idir=1,3 285 ibranch=idir+3*(iatom-1) 286 do kdir=1,3 287 k1 = kdir+3*(iatom-1) 288 displ_red(1,ibranch,jbranch) = displ_red(1,ibranch,jbranch) + & 289 & gprimd(kdir,idir)*displ(1,k1,jbranch,iFSqpt) 290 displ_red(2,ibranch,jbranch) = displ_red(2,ibranch,jbranch) + & 291 & gprimd(kdir,idir)*displ(2,k1,jbranch,iFSqpt) 292 end do 293 end do 294 end do 295 end do 296 297 tmpgam2 = reshape (gam_now, (/2,elph_ds%nbranch,elph_ds%nbranch/)) 298 call gam_mult_displ(elph_ds%nbranch, displ_red, tmpgam2, tmpgam1) 299 do jbranch=1,elph_ds%nbranch 300 eigval(jbranch) = tmpgam1(1, jbranch, jbranch) 301 end do 302 303 else if (elph_ds%ep_scalprod == 1) then 304 305 306 ! NOTE: in these calls gam_now and pheigvec do not have the right rank, but blas usually does not care 307 308 call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, c1, gam_now, 3*natom,& 309 & pheigvec(:,iFSqpt), 3*natom, c0, tmpgam1, 3*natom) 310 call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, c1, pheigvec(:,iFSqpt), 3*natom,& 311 & tmpgam1, 3*natom, c0, tmpgam2, 3*natom) 312 diagerr = zero 313 do ibranch=1,elph_ds%nbranch 314 eigval(ibranch) = tmpgam2(1,ibranch,ibranch) 315 do jbranch=1,ibranch-1 316 diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch)) 317 end do 318 do jbranch=ibranch+1,elph_ds%nbranch 319 diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch)) 320 end do 321 end do 322 if (diagerr > tol12) then 323 nerr=nerr+1 324 maxerr=max(diagerr, maxerr) 325 end if 326 end if ! end ep_scalprod if 327 328 ! Add all contributions from the phonon modes at this qpoint to a2f and the phonon dos. 329 do ibranch=1,elph_ds%nbranch 330 ! if (abs(phfrq(ibranch,iFSqpt)) < tol10) then 331 if ( abs(phfrq(ibranch,iFSqpt)) < tol7 .or. & 332 & (phfrq(ibranch,iFSqpt) < tol4 .and. qnorm2 > 0.03 )) then 333 ! note: this should depend on the velocity of sound, to accept acoustic modes! 334 a2fprefactor = zero 335 else 336 ! a2fprefactor = eigval (ibranch)/(two_pi*abs(phfrq(ibranch,iFSqpt))*elph_ds%n0(isppol)) 337 ! Use the dos_n0 at the lowest input temperature, assuming to be low 338 a2fprefactor = eigval (ibranch)/(two_pi*abs(phfrq(ibranch,iFSqpt))) 339 end if 340 341 omega = omega_min 342 tmpa2f(:) = zero 343 do iomega=1,elph_ds%na2f 344 xx = (omega-phfrq(ibranch,iFSqpt))*gaussfactor 345 omega = omega+domega 346 if (abs(xx) > gaussmaxarg) cycle 347 348 gaussval = gaussprefactor*exp(-xx*xx) 349 gtemp = gaussval*a2fprefactor 350 351 if (dabs(gtemp) < 1.0d-50) gtemp = zero 352 tmpa2f(iomega) = tmpa2f(iomega) + gtemp 353 end do 354 355 ! tmpa2f(:) = zero 356 ! do iomega=1,elph_ds%na2f 357 ! gtemp = a2fprefactor*elph_ds%k_phon%wtq(ibranch,iFSqpt,iomega) 358 ! if (dabs(gtemp) < 1.0d-50) gtemp = zero 359 ! tmpa2f(iomega) = tmpa2f(iomega) + gtemp 360 ! end do 361 362 tmp_a2f_1d_tr (:,itrtensor,isppol,ssp,ie) = tmp_a2f_1d_tr (:,itrtensor,isppol,ssp,ie) + tmpa2f(:) 363 364 end do ! end ibranch 365 end do ! end itrtensor 366 end do ! end iFSqpt - loop done in parallel 367 end do ! end isppol 368 end do ! ss' 369 end do ! n_pair 370 371 ! MG: FIXME: Why xmpi_world? besides only one CPU should perform IO (see below) 372 ! Likely this routine is never executed in parallel 373 call xmpi_sum (tmp_a2f_1d_tr, xmpi_world, ierr) 374 375 ABI_DEALLOCATE(coskr) 376 ABI_DEALLOCATE(sinkr) 377 378 do itemp=1,ntemper ! runs over termperature in K 379 do isppol=1,elph_ds%nsppol 380 do jj=1,tmp_nenergy**2 381 do ii=1,4 382 elph_tr_ds%a2f_1d_tr(:,:,isppol,ii,jj,itemp) = tmp_a2f_1d_tr(:,:,isppol,ii,jj)/elph_tr_ds%dos_n0(itemp,isppol) 383 end do 384 end do 385 end do 386 end do 387 388 ABI_DEALLOCATE(tmp_a2f_1d_tr) 389 390 !second 1 / elph_ds%k_phon%nkpt factor for the integration weights 391 elph_tr_ds%a2f_1d_tr = elph_tr_ds%a2f_1d_tr / elph_ds%k_phon%nkpt 392 393 if (elph_ds%ep_scalprod == 1) then 394 write(std_out,*) 'mka2f_tr: errors in diagonalization of gamma_tr with phon eigenvectors: ', nerr,maxerr 395 end if 396 397 !output the elph_tr_ds%a2f_1d_tr 398 fname = trim(elph_ds%elph_base_name) // '_A2F_TR' 399 if (open_file(fname,message,newunit=unit_a2f_tr,status='unknown') /= 0) then 400 MSG_ERROR(message) 401 end if 402 403 write (unit_a2f_tr,'(a)') '#' 404 write (unit_a2f_tr,'(a)') '# ABINIT package : a2f_tr file' 405 write (unit_a2f_tr,'(a)') '#' 406 write (unit_a2f_tr,'(a)') '# a2f_tr function integrated over the FS. omega in a.u.' 407 write (unit_a2f_tr,'(a,I10)') '# number of kpoints integrated over : ', elph_ds%k_phon%nkpt 408 write (unit_a2f_tr,'(a,I10)') '# number of energy points : ',elph_ds%na2f 409 write (unit_a2f_tr,'(a,E16.6,a,E16.6,a)') '# between omega_min = ', omega_min,' Ha and omega_max = ', omega_max, ' Ha' 410 write (unit_a2f_tr,'(a,E16.6)') '# and the smearing width for gaussians is ', elph_ds%a2fsmear 411 write (unit_a2f_tr,'(a)') '#' 412 413 !done with header 414 do isppol=1,elph_ds%nsppol 415 write (unit_a2f_tr,'(a,E16.6)') '# The DOS at Fermi level is ', elph_tr_ds%dos_n0(1,isppol) 416 omega = omega_min 417 do iomega=1,elph_ds%na2f 418 ! bxu, at which eps and eps' should I save it 419 ! better to save them all, but could be too many 420 write (unit_a2f_tr, '(10D16.6)') omega, elph_tr_ds%a2f_1d_tr(iomega,:,isppol,1,INT(elph_ds%n_pair/2)+1,1) 421 omega=omega+domega 422 end do 423 write (unit_a2f_tr,*) 424 end do !isppol 425 426 close (unit=unit_a2f_tr) 427 428 !calculation of transport properties 429 ABI_ALLOCATE(integrho,(elph_ds%na2f)) 430 ABI_ALLOCATE(tointegrho,(elph_ds%na2f)) 431 ABI_ALLOCATE(integrand_q00,(elph_ds%na2f)) 432 ABI_ALLOCATE(integrand_q01,(elph_ds%na2f)) 433 ABI_ALLOCATE(integrand_q11,(elph_ds%na2f)) 434 ABI_ALLOCATE(q00,(ntemper,3,3,elph_ds%nsppol)) 435 ABI_ALLOCATE(q01,(ntemper,3,3,elph_ds%nsppol)) 436 ABI_ALLOCATE(q11,(ntemper,3,3,elph_ds%nsppol)) 437 ABI_ALLOCATE(seebeck,(elph_ds%nsppol,ntemper,3,3)) 438 !ABI_ALLOCATE(rho_nm,(elph_ds%nsppol,ntemper,3,3)) 439 440 fname = trim(elph_ds%elph_base_name) // '_RHO' 441 if (open_file(fname,message,newunit=unit_rho,status='unknown') /= 0) then 442 MSG_ERROR(message) 443 end if 444 !print header to resistivity file 445 write (unit_rho,*) '# Resistivity as a function of temperature.' 446 write (unit_rho,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 447 write (unit_rho,*) '# ' 448 write (unit_rho,*) '# Columns are: ' 449 write (unit_rho,*) '# temperature[K] rho[au] rho [SI] rho/temp [au]' 450 write (unit_rho,*) '# ' 451 452 fname = trim(elph_ds%elph_base_name) // '_TAU' 453 if (open_file(fname,message,newunit=unit_tau,status='unknown') /= 0) then 454 MSG_ERROR(message) 455 end if 456 !print header to relaxation time file 457 write (unit_tau,*) '# Relaxation time as a function of temperature.' 458 write (unit_tau,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 459 write (unit_tau,*) '# ' 460 write (unit_tau,*) '# Columns are: ' 461 write (unit_tau,*) '# temperature[K] tau[au] tau [SI] ' 462 write (unit_tau,*) '# ' 463 464 fname = trim(elph_ds%elph_base_name) // '_SBK' 465 if (open_file(fname,message,newunit=unit_sbk,status='unknown') /= 0) then 466 MSG_ERROR(message) 467 end if 468 !print header to relaxation time file 469 write (unit_sbk,*) '# Seebeck Coefficint as a function of temperature.' 470 write (unit_sbk,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 471 write (unit_sbk,*) '# ' 472 write (unit_sbk,*) '# Columns are: ' 473 write (unit_sbk,*) '# temperature[K] S [au] S [SI] ' 474 write (unit_sbk,*) '# ' 475 476 fname = trim(elph_ds%elph_base_name) // '_WTH' 477 if (open_file(fname,message,newunit=unit_therm,status='unknown') /= 0) then 478 MSG_ERROR(message) 479 end if 480 481 !print header to thermal conductivity file 482 write (unit_therm,'(a)') '# Thermal conductivity/resistivity as a function of temperature.' 483 write (unit_therm,'(a)') '# the formalism is isotropic, so non-cubic crystals may be wrong' 484 write (unit_therm,'(a)') '# ' 485 write (unit_therm,'(a)') '# Columns are: ' 486 write (unit_therm,'(a)') '# temperature[K] thermal rho[au] thermal cond [au] thermal rho [SI] thermal cond [SI]' 487 write (unit_therm,'(a)') '# ' 488 489 fname = trim(elph_ds%elph_base_name) // '_LOR' 490 if (open_file(fname,message,newunit=unit_lor,status='unknown') /= 0) then 491 MSG_ERROR(message) 492 end if 493 494 !print header to lorentz file 495 write (unit_lor,*) '# Lorentz number as a function of temperature.' 496 write (unit_lor,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 497 write (unit_lor,*) '# ' 498 write (unit_lor,*) '# Columns are: ' 499 write (unit_lor,*) '# temperature[K] Lorentz number[au] Lorentz quantum = (pi*kb_HaK)**2/3' 500 write (unit_lor,*) '# ' 501 502 do isppol=1,elph_ds%nsppol 503 lambda_tr_trace = zero 504 do itrtensor=1,9 505 omega = omega_min 506 tointegrho = zero 507 do iomega=1,elph_ds%na2f 508 if(omega<=0) then 509 omega=omega+domega 510 cycle 511 end if 512 ! bxu, agian, which eps and eps' to use? 513 ! sometimes Ef is in the gap 514 tointegrho(iomega)=two*elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1)/omega 515 omega=omega+domega 516 end do 517 518 integrho = zero 519 call simpson_int(elph_ds%na2f,domega,tointegrho,integrho) 520 lambda_tr=integrho(elph_ds%na2f) 521 write (message, '(a,2i3,a,es16.6)' )& 522 & ' mka2f_tr: TRANSPORT lambda for isppol itrtensor', isppol, itrtensor, ' = ', lambda_tr 523 call wrtout(std_out,message,'COLL') 524 if (itrtensor == 1 .or. itrtensor == 5 .or. itrtensor == 9) lambda_tr_trace = lambda_tr_trace + lambda_tr 525 end do !end itrtensor do 526 527 lambda_tr_trace = lambda_tr_trace / three 528 write (message, '(a,i3,a,es16.6)' )& 529 & ' mka2f_tr: 1/3 trace of TRANSPORT lambda for isppol ', isppol, ' = ', lambda_tr_trace 530 call wrtout(std_out,message,'COLL') 531 call wrtout(ab_out,message,'COLL') 532 end do !end isppol do 533 534 !constant to change units of rho from au to SI 535 chgu=2.173969*1.0d-7 536 chtu=2.4188843265*1.0d-17 537 chsu=8.617343101*1.0d-5 ! au to J/C/K 538 chwu=9.270955772*1.0d-5 ! au to mK/W 539 540 !change the fermi level to zero, as required for q01 to vanish. 541 tmp_fermie = elph_ds%fermie 542 !Get Q00, Q01, Q11, and derive rho, tau 543 q00 = zero 544 q01 = zero 545 q11 = zero 546 ! prepare s1 and s2 arrays 547 s1 = (/1, 1, -1, -1/) 548 s2 = (/1, -1, 1, -1/) 549 550 do isppol=1,elph_ds%nsppol 551 do icomp=1, 3 552 do jcomp=1, 3 553 itrtensor=(icomp-1)*3+jcomp 554 555 write(unit_rho,*) '# Rho for isppol, itrten = ', isppol, itrtensor 556 write(unit_tau,*) '# Tau for isppol, itrten = ', isppol, itrtensor 557 558 do itemp=1,ntemper ! runs over termperature in K 559 Temp=tempermin+temperinc*dble(itemp) 560 tmp_veloc_sq0 = sqrt(elph_tr_ds%veloc_sq0(itemp,icomp,isppol)*elph_tr_ds%veloc_sq0(itemp,jcomp,isppol)) 561 562 integrand_q00 = zero 563 integrand_q01 = zero 564 integrand_q11 = zero 565 566 omega = omega_min 567 do iomega=1,elph_ds%na2f 568 if(omega .le. 0) then 569 omega=omega+domega 570 cycle 571 end if 572 xtr=omega/(kb_HaK*Temp) 573 occ_omega=1.0_dp/(exp(xtr)-1.0_dp) 574 575 tmp_veloc_sq1 = zero 576 tmp_veloc_sq2 = zero 577 do ie1=1,elph_ds%nenergy 578 e1 = elph_tr_ds%en_all(isppol,ie1) 579 580 !! BXU, the tolerance here needs to be used with caution 581 !! which depends on the dimensions of the system 582 !! E.g. in 2D system, DOS can be much smaller 583 if (elph_tr_ds%dos_n(ie1,isppol)/natom .lt. 0.1d0) cycle ! energy in the gap forbidden 584 585 xtr=(e1-tmp_fermie)/(kb_HaK*Temp) 586 occ_e1=1.0_dp/(exp(xtr)+1.0_dp) 587 588 e2 = e1 + omega 589 xtr=(e2-tmp_fermie)/(kb_HaK*Temp) 590 occ_e2=1.0_dp/(exp(xtr)+1.0_dp) 591 ! Do we need to change the fermie to the one with T dependence? 592 ! find which ie2 give the closest energy 593 if (e2 .gt. elph_tr_ds%en_all(isppol,elph_ds%nenergy)) then 594 ie2 = 0 595 cycle 596 else 597 ie_tmp = 1 598 diff = dabs(e2-elph_tr_ds%en_all(isppol,1)) 599 do ie2 = 2, elph_ds%nenergy 600 if (dabs(e2-elph_tr_ds%en_all(isppol,ie2)) .lt. diff) then 601 diff = dabs(e2-elph_tr_ds%en_all(isppol,ie2)) 602 ie_tmp = ie2 603 end if 604 end do 605 ie2 = ie_tmp 606 607 if (e2 < elph_tr_ds%en_all(isppol,ie2)) then 608 ie2_right = ie2 609 ie2_left = ie2-1 610 else 611 ie2_right = ie2+1 612 ie2_left = ie2 613 end if 614 615 end if 616 617 if (elph_tr_ds%dos_n(ie2,isppol)/natom .lt. 0.1d0) cycle 618 619 tointegq00 = zero 620 tointegq01 = zero 621 tointegq11 = zero 622 623 ! BXU linear interpolation of volec_sq and dos_n 624 xe=(e2-elph_tr_ds%en_all(isppol,ie2_left))/ & 625 & (elph_tr_ds%en_all(isppol,ie2_right)-elph_tr_ds%en_all(isppol,ie2_left)) 626 veloc_sq_icomp = elph_tr_ds%veloc_sq(icomp,isppol,ie2_left)*(1.0d0-xe) + & 627 & elph_tr_ds%veloc_sq(icomp,isppol,ie2_right)*xe 628 veloc_sq_jcomp = elph_tr_ds%veloc_sq(jcomp,isppol,ie2_left)*(1.0d0-xe) + & 629 & elph_tr_ds%veloc_sq(jcomp,isppol,ie2_right)*xe 630 dos_n_e2 = elph_tr_ds%dos_n(ie2_left,isppol)*(1.0d0-xe) + & 631 & elph_tr_ds%dos_n(ie2_right,isppol)*xe 632 633 tmp_veloc_sq1 = sqrt(elph_tr_ds%veloc_sq(icomp,isppol,ie1)*elph_tr_ds%veloc_sq(jcomp,isppol,ie1)) 634 ! tmp_veloc_sq2 = sqrt(elph_tr_ds%veloc_sq(icomp,isppol,ie2)*elph_tr_ds%veloc_sq(jcomp,isppol,ie2)) 635 tmp_veloc_sq2 = sqrt(veloc_sq_icomp*veloc_sq_jcomp) 636 637 ! ie_1 = (ie1-1)*elph_ds%nenergy + ie2 638 ! ie_2 = (ie2-1)*elph_ds%nenergy + ie1 639 ie_1 = pair2red(ie1,ie2) 640 ie_2 = pair2red(ie2,ie1) 641 if (ie_1 .eq. 0 .or. ie_2 .eq. 0) then 642 MSG_BUG('CHECK pair2red!') 643 end if 644 645 do ssp=1,4 ! (s,s'=+/-1, condense the indices) 646 647 nv1 = 1.0_dp/(elph_tr_ds%dos_n(ie1,isppol)*sqrt(tmp_veloc_sq1)) 648 sigma1 = sqrt(3.0_dp)*(e1-tmp_fermie)/(pi*Temp*kb_HaK) 649 !DEBUG 650 if (elph_ds%ep_lova .eq. 1) then 651 nv1 = 1.0_dp/(elph_tr_ds%dos_n0(itemp,isppol)*sqrt(tmp_veloc_sq0)) 652 sigma1 = sqrt(3.0_dp)*(e1-tmp_fermie)/(pi*Temp*kb_HaK) 653 end if 654 !ENDDEBUG 655 656 tointegq00_1 = zero 657 tointegq01_1 = zero 658 tointegq11_1 = zero 659 660 !DEBUG 661 if (elph_ds%ep_lova .eq. 1) then 662 nv2 = 1.0_dp/(elph_tr_ds%dos_n0(itemp,isppol)*sqrt(tmp_veloc_sq0)) 663 sigma2 = sqrt(3.0_dp)*(e2-tmp_fermie)/(pi*Temp*kb_HaK) 664 j00 = (nv1 + s1(ssp)*nv2)*(nv1 + s2(ssp)*nv2)/4.0_dp 665 j01 = (nv1 + s1(ssp)*nv2)*(nv1*sigma1 + s2(ssp)*nv2*sigma2)/4.0_dp 666 j11 = (nv1*sigma1 + s1(ssp)*nv2*sigma2)*(nv1*sigma1 + s2(ssp)*nv2*sigma2)/4.0_dp 667 tointegq00_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* & 668 & occ_e1*(1.0_dp-occ_e2)*j00*occ_omega 669 tointegq01_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* & 670 & occ_e1*(1.0_dp-occ_e2)*j01*occ_omega 671 tointegq11_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* & 672 & occ_e1*(1.0_dp-occ_e2)*j11*occ_omega 673 !END DEBUG 674 else if (elph_ds%ep_lova .eq. 0) then 675 nv2 = 1.0_dp/(dos_n_e2*sqrt(tmp_veloc_sq2)) 676 sigma2 = sqrt(3.0_dp)*(e2-tmp_fermie)/(pi*Temp*kb_HaK) 677 j00 = (nv1 + s1(ssp)*nv2)*(nv1 + s2(ssp)*nv2)/4.0_dp 678 j01 = (nv1 + s1(ssp)*nv2)*(nv1*sigma1 + s2(ssp)*nv2*sigma2)/4.0_dp 679 j11 = (nv1*sigma1 + s1(ssp)*nv2*sigma2)*(nv1*sigma1 + s2(ssp)*nv2*sigma2)/4.0_dp 680 ! bxu TEST 681 if (debug) then 682 if (ssp .eq. 1 .and. itrtensor .eq. 1) then 683 write(21,"(3i5,4E20.12)") iomega, ie1, ie2, & 684 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp), j01, & 685 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01, & 686 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01*occ_e1*(1.0_dp-occ_e2)*occ_omega 687 end if 688 if (ssp .eq. 2 .and. itrtensor .eq. 1) then 689 write(22,"(3i5,4E20.12)") iomega, ie1, ie2, & 690 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp), j01, & 691 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01, & 692 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01*occ_e1*(1.0_dp-occ_e2)*occ_omega 693 end if 694 if (ssp .eq. 3 .and. itrtensor .eq. 1) then 695 write(23,"(3i5,4E20.12)") iomega, ie1, ie2, & 696 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp), j01, & 697 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01, & 698 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01*occ_e1*(1.0_dp-occ_e2)*occ_omega 699 end if 700 if (ssp .eq. 4 .and. itrtensor .eq. 1) then 701 write(24,"(3i5,4E20.12)") iomega, ie1, ie2, & 702 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp), j01, & 703 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01, & 704 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01*occ_e1*(1.0_dp-occ_e2)*occ_omega 705 end if 706 end if 707 tointegq00_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)* & 708 & occ_e1*(1.0_dp-occ_e2)*j00*occ_omega 709 tointegq01_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)* & 710 & occ_e1*(1.0_dp-occ_e2)*j01*occ_omega 711 tointegq11_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)* & 712 & occ_e1*(1.0_dp-occ_e2)*j11*occ_omega 713 end if 714 715 tointegq00_2 = zero 716 tointegq01_2 = zero 717 tointegq11_2 = zero 718 719 !DEBUG 720 if (elph_ds%ep_lova .eq. 1) then 721 nv2 = 1.0_dp/(elph_tr_ds%dos_n0(itemp,isppol)*sqrt(tmp_veloc_sq0)) 722 sigma2 = sqrt(3.0_dp)*(e2-tmp_fermie)/(pi*Temp*kb_HaK) 723 j00 = (nv2 + s1(ssp)*nv1)*(nv2 + s2(ssp)*nv1)/4.0_dp 724 j01 = (nv2 + s1(ssp)*nv1)*(nv2*sigma2 + s2(ssp)*nv1*sigma1)/4.0_dp 725 j11 = (nv2*sigma2 + s1(ssp)*nv1*sigma1)*(nv2*sigma2 + s2(ssp)*nv1*sigma1)/4.0_dp 726 tointegq00_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* & 727 & occ_e1*(1.0_dp-occ_e2)*j00*(occ_omega+1) 728 tointegq01_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* & 729 & occ_e1*(1.0_dp-occ_e2)*j01*(occ_omega+1) 730 tointegq11_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* & 731 & occ_e1*(1.0_dp-occ_e2)*j11*(occ_omega+1) 732 !END DEBUG 733 else if (elph_ds%ep_lova .eq. 0) then 734 nv2 = 1.0_dp/(dos_n_e2*sqrt(tmp_veloc_sq2)) 735 sigma2 = sqrt(3.0_dp)*(e2-tmp_fermie)/(pi*Temp*kb_HaK) 736 j00 = (nv2 + s1(ssp)*nv1)*(nv2 + s2(ssp)*nv1)/4.0_dp 737 j01 = (nv2 + s1(ssp)*nv1)*(nv2*sigma2 + s2(ssp)*nv1*sigma1)/4.0_dp 738 j11 = (nv2*sigma2 + s1(ssp)*nv1*sigma1)*(nv2*sigma2 + s2(ssp)*nv1*sigma1)/4.0_dp 739 !DEBUG bxu TEST 740 if (debug) then 741 if (ssp .eq. 1 .and. itrtensor .eq. 1) then 742 write(31,"(3i5,4E20.12)") iomega, ie2, ie1, & 743 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp), j01, & 744 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01, & 745 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01*occ_e2*(1.0_dp-occ_e1)*(occ_omega+1) 746 end if 747 if (ssp .eq. 2 .and. itrtensor .eq. 1) then 748 write(32,"(3i5,4E20.12)") iomega, ie2, ie1, & 749 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp), j01, & 750 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01, & 751 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01*occ_e2*(1.0_dp-occ_e1)*(occ_omega+1) 752 end if 753 if (ssp .eq. 3 .and. itrtensor .eq. 1) then 754 write(33,"(3i5,4E20.12)") iomega, ie2, ie1, & 755 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp), j01, & 756 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01, & 757 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01*occ_e2*(1.0_dp-occ_e1)*(occ_omega+1) 758 end if 759 if (ssp .eq. 4 .and. itrtensor .eq. 1) then 760 write(34,"(3i5,4E20.12)") iomega, ie2, ie1, & 761 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp), j01, & 762 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01, & 763 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01*occ_e2*(1.0_dp-occ_e1)*(occ_omega+1) 764 end if 765 end if 766 !ENDDEBUG 767 tointegq00_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)* & 768 & occ_e2*(1.0_dp-occ_e1)*j00*(occ_omega+1) 769 tointegq01_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)* & 770 & occ_e2*(1.0_dp-occ_e1)*j01*(occ_omega+1) 771 tointegq11_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)* & 772 & occ_e2*(1.0_dp-occ_e1)*j11*(occ_omega+1) 773 end if ! elph_ds%ep_lova 774 775 tointegq00 = tointegq00 + tointegq00_1 + tointegq00_2 776 tointegq01 = tointegq01 + tointegq01_1 + tointegq01_2 777 tointegq11 = tointegq11 + tointegq11_1 + tointegq11_2 778 779 end do ! ss' = 4 780 integrand_q00(iomega) = integrand_q00(iomega) + elph_tr_ds%de_all(isppol,ie1)*tointegq00 781 integrand_q01(iomega) = integrand_q01(iomega) + elph_tr_ds%de_all(isppol,ie1)*tointegq01 782 integrand_q11(iomega) = integrand_q11(iomega) + elph_tr_ds%de_all(isppol,ie1)*tointegq11 783 end do ! ie1 ~ 20 784 omega=omega+domega 785 q00(itemp,icomp,jcomp,isppol) = q00(itemp,icomp,jcomp,isppol) +& 786 & domega*integrand_q00(iomega) 787 q01(itemp,icomp,jcomp,isppol) = q01(itemp,icomp,jcomp,isppol) +& 788 & domega*integrand_q01(iomega) 789 q11(itemp,icomp,jcomp,isppol) = q11(itemp,icomp,jcomp,isppol) +& 790 & domega*integrand_q11(iomega) 791 end do ! omega ~ 400 792 793 q00(itemp,icomp,jcomp,isppol)=q00(itemp,icomp,jcomp,isppol)* & 794 & ucvol*two_pi*elph_tr_ds%dos_n0(itemp,isppol)/(kb_HaK*Temp) 795 q01(itemp,icomp,jcomp,isppol)=q01(itemp,icomp,jcomp,isppol)* & 796 & ucvol*two_pi*elph_tr_ds%dos_n0(itemp,isppol)/(kb_HaK*Temp) 797 q11(itemp,icomp,jcomp,isppol)=q11(itemp,icomp,jcomp,isppol)* & 798 & ucvol*two_pi*elph_tr_ds%dos_n0(itemp,isppol)/(kb_HaK*Temp) 799 800 rho = 0.5_dp*q00(itemp,icomp,jcomp,isppol) 801 ! is tau energy dependent? 802 tau = 2.0d0*ucvol/(q00(itemp,icomp,jcomp,isppol)*elph_tr_ds%dos_n0(itemp,isppol)*tmp_veloc_sq0) 803 write(unit_rho,'(4D20.10)')temp,rho,rho*chgu,rho/temp 804 write(unit_tau,'(3D20.10)')temp,tau,tau*chtu 805 rho_T(itemp)=rho 806 end do ! temperature = 1? 807 write(unit_rho,*) 808 write(unit_tau,*) 809 810 end do ! jcomp = 3 811 end do ! icomp = 3 812 end do ! isppol = 2 813 814 !----------------------------- 815 816 seebeck = zero 817 !rho_nm = zero 818 819 !do isppol=1,elph_ds%nsppol 820 !do itemp=1,ntemper 821 !q11_inv(:,:)=q11(itemp,:,:,isppol) 822 !call matrginv(q11_inv,3,3) 823 !do icomp=1,3 824 !do jcomp=1,3 825 !do kcomp=1,3 826 !seebeck(isppol,itemp,icomp,jcomp) = seebeck(isppol,itemp,icomp,jcomp) + & 827 !& q01(itemp,icomp,kcomp,isppol)*q11_inv(kcomp,jcomp) 828 !end do 829 !end do 830 !end do 831 !end do 832 !end do 833 834 do isppol=1,elph_ds%nsppol 835 do itemp=1,ntemper 836 q11_inv(:,:)=q11(itemp,:,:,isppol) 837 call matrginv(q11_inv,3,3) 838 call DGEMM('N','N',3,3,3,one,q01(itemp,:,:,isppol),3,q11_inv,& 839 & 3,zero,seebeck(isppol,itemp,:,:),3) 840 ! call DGEMM('N','N',3,3,3,one,seebeck(isppol,itemp,:,:),3,q01(itemp,:,:,isppol),& 841 ! & 3,zero,rho_nm(isppol,itemp,:,:),3) 842 end do 843 end do 844 pref_s = pi/sqrt(3.0_dp) 845 seebeck=pref_s*seebeck 846 847 !fullq = zero 848 !do icomp=1,3 849 !do jcomp=1,3 850 !fullq(icomp,jcomp) = q00(1,icomp,jcomp,1) 851 !end do 852 !end do 853 !do icomp=1,3 854 !do jcomp=4,6 855 !fullq(icomp,jcomp) = q01(1,icomp,jcomp-3,1) 856 !end do 857 !end do 858 !do icomp=4,6 859 !do jcomp=1,3 860 !fullq(icomp,jcomp) = q01(1,icomp-3,jcomp,1) 861 !end do 862 !end do 863 !do icomp=4,6 864 !do jcomp=4,6 865 !fullq(icomp,jcomp) = q11(1,icomp-3,jcomp-3,1) 866 !end do 867 !end do 868 !write(*,*)' fullq' 869 !write(*,"(6E20.12)") (fullq(1,jcomp),jcomp=1,6) 870 !write(*,"(6E20.12)") (fullq(2,jcomp),jcomp=1,6) 871 !write(*,"(6E20.12)") (fullq(3,jcomp),jcomp=1,6) 872 !write(*,"(6E20.12)") (fullq(4,jcomp),jcomp=1,6) 873 !write(*,"(6E20.12)") (fullq(5,jcomp),jcomp=1,6) 874 !write(*,"(6E20.12)") (fullq(6,jcomp),jcomp=1,6) 875 write(message,'(a)') 'q00:' 876 call wrtout(std_out,message,'COLL') 877 write(message,'(3E20.12)') (q00(1,1,jcomp,1),jcomp=1,3) 878 call wrtout(std_out,message,'COLL') 879 write(message,'(3E20.12)') (q00(1,2,jcomp,1),jcomp=1,3) 880 call wrtout(std_out,message,'COLL') 881 write(message,'(3E20.12)') (q00(1,3,jcomp,1),jcomp=1,3) 882 call wrtout(std_out,message,'COLL') 883 write(message,'(a)') 'q01:' 884 call wrtout(std_out,message,'COLL') 885 write(message,'(3E20.12)') (q01(1,1,jcomp,1),jcomp=1,3) 886 call wrtout(std_out,message,'COLL') 887 write(message,'(3E20.12)') (q01(1,2,jcomp,1),jcomp=1,3) 888 call wrtout(std_out,message,'COLL') 889 write(message,'(3E20.12)') (q01(1,3,jcomp,1),jcomp=1,3) 890 call wrtout(std_out,message,'COLL') 891 write(message,'(a)') 'q11, q11_inv:' 892 call wrtout(std_out,message,'COLL') 893 write(message,'(6E20.12)') (q11(1,1,jcomp,1),jcomp=1,3),(q11_inv(1,jcomp),jcomp=1,3) 894 call wrtout(std_out,message,'COLL') 895 write(message,'(6E20.12)') (q11(1,2,jcomp,1),jcomp=1,3),(q11_inv(2,jcomp),jcomp=1,3) 896 call wrtout(std_out,message,'COLL') 897 write(message,'(6E20.12)') (q11(1,3,jcomp,1),jcomp=1,3),(q11_inv(3,jcomp),jcomp=1,3) 898 call wrtout(std_out,message,'COLL') 899 !q11_inv = zero 900 !do icomp = 1, 3 901 !q11_inv(icomp,icomp) = 2.0_dp 902 !end do 903 904 !call matrginv(fullq,6,6) 905 906 !do isppol=1,elph_ds%nsppol 907 !do itemp=1,ntemper 908 !rho_nm(isppol,itemp,:,:) = q00(itemp,:,:,isppol) - rho_nm(isppol,itemp,:,:) 909 !end do 910 !end do 911 !rho_nm = 0.5_dp*rho_nm 912 913 !Output of Seebeck coefficient 914 do isppol=1,elph_ds%nsppol 915 do icomp=1,3 916 do jcomp=1,3 917 itrtensor=(icomp-1)*3+jcomp 918 write(unit_sbk,*) '# Seebeck for isppol, itrten = ', isppol, itrtensor 919 ! write(88,*) '# Rho for isppol, itrten = ', isppol, itrtensor 920 ! write(89,*) '# Rho for isppol, itrten = ', isppol, itrtensor 921 do itemp=1,ntemper 922 Temp=tempermin+temperinc*dble(itemp) 923 write(unit_sbk,'(3D20.10)')temp, seebeck(isppol,itemp,icomp,jcomp), seebeck(isppol,itemp,icomp,jcomp)*chsu 924 ! write(88,'(3D20.10)')temp, rho_nm(isppol,itemp,icomp,jcomp), rho_nm(isppol,itemp,icomp,jcomp)*chgu 925 ! write(89,'(3D20.10)')temp, 0.5_dp/fullq(1,1), 0.5_dp*chgu/fullq(1,1) 926 end do ! temperature 927 write(unit_sbk,*) 928 ! write(88,*) 929 ! write(89,*) 930 end do ! jcomp 931 end do ! icomp 932 end do ! isppol 933 934 !Get thermal resistivity, based on eqn. (52) in Allen's PRB 17, 3725 (1978) 935 !WARNING: before 6.13.1 the thermal resistivity and Lorentz number were not in 936 !atomic units, BUT the SI units are good. 937 pref_w = 3.0_dp/(2.0_dp*pi**2.0d0) 938 do isppol=1,elph_ds%nsppol 939 do icomp=1, 3 940 do jcomp=1, 3 941 itrtensor=(icomp-1)*3+jcomp 942 943 write(unit_therm,*) '# Thermal resistivity for isppol, itrten= ', isppol 944 write(unit_lor,*) '# Lorentz coefficient for isppol, itrten= ', isppol 945 946 do itemp=1,ntemper 947 948 Temp=tempermin + temperinc*dble(itemp) 949 950 wtherm = pref_w*q11(itemp,icomp,jcomp,isppol)/(kb_HaK*Temp) 951 952 ! write(unit_therm,'(5D20.10)')temp,wtherm,1./wtherm,wtherm/3.4057d9,1./(wtherm) *3.4057d9 953 write(unit_therm,'(5D20.10)')temp,wtherm,1.0_dp/wtherm,wtherm*chwu,1.0_dp/(wtherm*chwu) 954 955 lorentz=rho_T(itemp)/(wtherm*kb_HaK*Temp) 956 write(unit_lor,*)temp,lorentz,lor0 957 958 end do 959 write(unit_therm,*) 960 write(unit_lor,*) 961 end do ! jcomp 962 end do ! icomp 963 end do ! isppol 964 965 966 ABI_DEALLOCATE(phfrq) 967 ABI_DEALLOCATE(displ) 968 ABI_DEALLOCATE(pheigvec) 969 ABI_DEALLOCATE(integrand_q00) 970 ABI_DEALLOCATE(integrand_q01) 971 ABI_DEALLOCATE(integrand_q11) 972 ABI_DEALLOCATE(q00) 973 ABI_DEALLOCATE(q01) 974 ABI_DEALLOCATE(q11) 975 ABI_DEALLOCATE(seebeck) 976 ABI_DEALLOCATE(rho_T) 977 ABI_DEALLOCATE(integrho) 978 ABI_DEALLOCATE(tointegrho) 979 980 close (unit=unit_lor) 981 close (unit=unit_rho) 982 close (unit=unit_tau) 983 close (unit=unit_sbk) 984 close (unit=unit_therm) 985 986 ABI_DEALLOCATE(elph_ds%k_fine%wtq) 987 ABI_DEALLOCATE(elph_ds%k_phon%wtq) 988 989 ABI_DEALLOCATE(elph_tr_ds%a2f_1d_tr) 990 991 ABI_DEALLOCATE(elph_tr_ds%gamma_qpt_tr) 992 ABI_DEALLOCATE(elph_tr_ds%gamma_rpt_tr) 993 write(std_out,*) ' mka2f_tr : end ' 994 995 end subroutine mka2f_tr