TABLE OF CONTENTS
ABINIT/get_tau_k [ Functions ]
NAME
get_tau_k
FUNCTION
Calculate the k-dependent relaxation time due to EPC. Impelementation based on derivation from Grmvall's book or OD Restrepo's paper (PRB 94 212103 (2009))
COPYRIGHT
Copyright (C) 2013-2018 ABINIT group (BXu) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
Cryst<crystal_t>=Info on the unit cell and on its symmetries. Ifc<ifc_type>=Object containing the interatomic force constants. elph_ds = elphon datastructure with data and dimensions eigenGS = Ground State eigenvalues max_occ = maximal occupancy for a band
OUTPUT
tau_k(nsppol,nkptirr,nband)=mode relaxation time due to electron phonono coupling rate_e(nene)= scattering rate due to electron phonono coupling vs. energy
PARENTS
elphon
CHILDREN
dgemm,ebands_prtbltztrp_tau_out,ebands_update_occ,ep_el_weights ep_ph_weights,ftgam,ftgam_init,gam_mult_displ,ifc_fourq,matrginv mkqptequiv,phdispl_cart2red,spline,splint,wrtout,xmpi_sum,zgemm
SOURCE
37 #if defined HAVE_CONFIG_H 38 #include "config.h" 39 #endif 40 41 #include "abi_common.h" 42 43 subroutine get_tau_k(Cryst,ifc,Bst,elph_ds,elph_tr_ds,eigenGS,max_occ) 44 45 use defs_basis 46 use defs_elphon 47 use defs_datatypes 48 use m_kptrank 49 use m_errors 50 use m_profiling_abi 51 use m_xmpi 52 use m_splines 53 use m_ifc 54 use m_ebands 55 56 use m_io_tools, only : open_file 57 use m_abilasi, only : matrginv 58 use m_geometry, only : phdispl_cart2red 59 use m_dynmat, only : ftgam_init, ftgam 60 use m_crystal, only : crystal_t 61 62 !This section has been created automatically by the script Abilint (TD). 63 !Do not modify the following lines by hand. 64 #undef ABI_FUNC 65 #define ABI_FUNC 'get_tau_k' 66 use interfaces_14_hidewrite 67 use interfaces_77_ddb, except_this_one => get_tau_k 68 !End of the abilint section 69 70 implicit none 71 72 !Arguments ------------------------------------ 73 type(crystal_t),intent(in) :: Cryst 74 type(ifc_type),intent(in) :: ifc 75 type(ebands_t),intent(inout) :: Bst 76 type(elph_type),intent(inout) :: elph_ds 77 type(elph_tr_type), intent(inout) :: elph_tr_ds 78 real(dp),intent(in) :: max_occ 79 real(dp),intent(in) :: eigenGS(elph_ds%nband,elph_ds%k_phon%nkpt,elph_ds%nsppol) 80 81 !Local variables------------------------------- 82 !scalars 83 character(len=500) :: message 84 character(len=fnlen) :: fname 85 integer :: ntemper,nsppol,nbranch,nband,natom 86 integer :: nkpt,nqpt,nkptirr,nqptirr,new_nkptirr 87 integer :: isppol,iFSkpt,iFSqpt,iqpt,iqpt_fullbz,imqpt_fullbz,ikpt_kpq,ikpt_kmq 88 integer :: iband,jband,jpband,jbeff,ibranch,jbranch,itemp 89 integer :: irec,ierr,nrpt,ik_this_proc 90 integer :: unit_tau,unit_invtau 91 integer :: nene,nene_all,iene,iene_fine,unit_taue,unit_mfp 92 integer :: icomp,jcomp,itensor 93 integer :: ikpt_irr,iomega,unit_cond,unit_therm,unit_sbk 94 integer :: nskip,nspline 95 real(dp) :: occ_omega,occ_e 96 real(dp) :: xx,Temp,therm_factor 97 real(dp) :: factor,dfermide 98 real(dp) :: e_k,chu_tau,rate_e,mfp_e 99 real(dp) :: ene,enemin,enemax,deltaene 100 real(dp) :: omega,omega_min,omega_max,domega 101 real(dp) :: diagerr 102 real(dp) :: chu_mfp,chu_cond,chu_cth,chu_sbk,femto 103 real(dp) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch) 104 real(dp) :: eigval(elph_ds%nbranch),eigval2(elph_ds%nbranch) 105 real(dp) :: imeigval(elph_ds%nbranch) 106 real(dp) :: tmp_wtkpq, tmp_wtkmq, tol_wtk 107 real(dp) :: yp1,ypn 108 !arrays 109 integer,allocatable :: FSfullpktofull(:,:),mqtofull(:) 110 integer,allocatable :: kpttokpt(:,:,:) 111 real(dp) :: cond_inv(3,3) 112 real(dp),allocatable :: fermie(:) 113 real(dp),allocatable :: tmp_eigenGS(:,:,:) 114 real(dp),allocatable :: tmp_gkk_qpt(:,:,:),tmp_gkk_rpt(:,:,:),tmp_gkk_kpt(:,:) 115 real(dp),allocatable :: tmp_gkk_kpt2(:,:,:), gkk_kpt(:,:,:) 116 real(dp),allocatable :: tau_k(:,:,:,:),inv_tau_k(:,:,:,:),tmp_tau_k(:,:,:,:) 117 real(dp),allocatable :: phfrq(:,:),pheigvec(:,:) 118 real(dp),allocatable :: displ(:,:,:,:) 119 real(dp),allocatable :: a2f_2d(:),a2f_2d2(:) 120 real(dp),allocatable :: tmp_wtk(:,:,:,:),tmp2_wtk(:),tmp_wtk1(:),tmp_wtk2(:) 121 real(dp),allocatable :: ene_pt(:),ene_ptfine(:),ff2(:) 122 real(dp),allocatable :: wtq(:,:,:),tmp_wtq(:,:,:),tmp2_wtq(:,:) 123 real(dp),allocatable :: dos_e(:,:) 124 real(dp),allocatable :: coskr1(:,:),sinkr1(:,:) 125 real(dp),allocatable :: coskr2(:,:),sinkr2(:,:) 126 real(dp),allocatable :: cond_e(:,:,:,:),cond(:,:,:,:),sbk(:,:,:,:),seebeck(:,:,:,:),cth(:,:,:,:) 127 128 ! ************************************************************************* 129 130 write(std_out,*) 'get_tau_k : enter ' 131 132 nrpt = ifc%nrpt 133 natom = cryst%natom 134 135 nsppol = elph_ds%nsppol 136 nbranch = elph_ds%nbranch 137 nband = elph_ds%ngkkband 138 nkpt = elph_ds%k_phon%nkpt 139 nqpt = elph_ds%nqpt_full 140 nkptirr = elph_ds%k_phon%nkptirr 141 new_nkptirr = elph_ds%k_phon%new_nkptirr 142 nqptirr = elph_ds%nqptirred 143 ntemper = elph_ds%ntemper 144 nene = 2*elph_ds%na2f-1 ! only need e_k +- omega_max range, take deltaene=delta_oemga 145 146 chu_tau = 2.4188843265*1.0d-17 147 chu_mfp = 5.291772*1.0d-11 148 chu_cond = 4.59988159904764*1.0d6 149 chu_cth = 1.078637439971599*1.0d4 150 chu_sbk = 8.617343101*1.0d-5 151 femto = 1.0d-15 152 153 tol_wtk = tol7/nkptirr/nband 154 155 ABI_ALLOCATE(fermie ,(ntemper)) 156 ABI_ALLOCATE(tmp_gkk_qpt ,(2,nbranch**2,nqpt)) 157 ABI_ALLOCATE(tmp_gkk_rpt ,(2,nbranch**2,nrpt)) 158 ABI_ALLOCATE(tmp_gkk_kpt ,(2,nbranch**2)) 159 ABI_ALLOCATE(tmp_gkk_kpt2 ,(2,nbranch,nbranch)) 160 ABI_ALLOCATE(gkk_kpt ,(2,nbranch,nbranch)) 161 ABI_ALLOCATE(a2f_2d, (nene)) 162 ABI_ALLOCATE(a2f_2d2, (nene)) 163 ABI_ALLOCATE(inv_tau_k, (ntemper,nsppol,nkpt,nband)) 164 ABI_ALLOCATE(tau_k, (ntemper,nsppol,nkpt,nband)) 165 ABI_ALLOCATE(tmp_tau_k ,(ntemper,nsppol,new_nkptirr,nband)) 166 167 if (elph_ds%gkqwrite == 0) then 168 call wrtout(std_out,' get_tau_k : keeping gkq matrices in memory','COLL') 169 else if (elph_ds%gkqwrite == 1) then 170 fname=trim(elph_ds%elph_base_name) // '_GKKQ' 171 write (message,'(2a)')' get_tau_k : reading gkq matrices from file ',trim(fname) 172 call wrtout(std_out,message,'COLL') 173 else 174 write (message,'(a,i0)')' Wrong value for gkqwrite = ',elph_ds%gkqwrite 175 MSG_BUG(message) 176 end if 177 178 !========================================================= 179 !Get equivalence between a kpt_phon pair and a qpt in qpt_full 180 !only works if the qpt grid is complete (identical to 181 !the kpt one, with a basic shift of (0,0,0) 182 !========================================================= 183 184 !mapping of k + q onto k' for k and k' in full BZ 185 !for dense k grid 186 ABI_ALLOCATE(FSfullpktofull,(nkpt,nkpt)) 187 ABI_ALLOCATE(mqtofull,(nkpt)) 188 189 !kpttokpt(itim,isym,iqpt) = kpoint index which transforms to ikpt under isym and with time reversal itim. 190 ABI_ALLOCATE(kpttokpt,(2,Cryst%nsym,nkpt)) 191 192 call wrtout(std_out,'get_tau_k: calling mkqptequiv to set up the FS kpoint set',"COLL") 193 194 call mkqptequiv (FSfullpktofull,Cryst,elph_ds%k_phon%kpt,nkpt,nkpt,kpttokpt,elph_ds%k_phon%kpt,mqtofull) 195 196 !========================================================= 197 !========================================================= 198 199 omega_max = elph_ds%omega_max 200 omega_min = elph_ds%omega_min 201 domega = elph_ds%domega 202 enemax = maxval(eigenGS(elph_ds%maxFSband,:,:)) 203 enemin = minval(eigenGS(elph_ds%minFSband,:,:)) 204 205 if (enemin < (elph_ds%fermie-0.2)) then 206 enemin = elph_ds%fermie-0.2 207 end if 208 if (enemax > (elph_ds%fermie+0.2)) then 209 enemax = elph_ds%fermie+0.2 210 end if 211 212 nspline = elph_ds%ep_nspline 213 nene_all = INT((enemax-enemin+domega)/(nspline*domega)) + 1 214 deltaene = domega 215 write(std_out,*) 'E_min= ',enemin, 'E_max= ',enemax 216 write(std_out,*) 'Number of energy points= ',nene_all 217 write(std_out,'(a,I8)') 'scale factor for spline interpolation in RTA = ', elph_ds%ep_nspline 218 write(std_out,*) 'delta_ene before spline interpolation= ',deltaene*nspline 219 write(std_out,*) 'delta_ene after spline interpolation= ',deltaene 220 write(std_out,*) 'Omega_min= ',omega_min, 'Omega_max= ',omega_max 221 write(std_out,*) 'Number of phonon points= ',elph_ds%na2f 222 write(std_out,*) 'delta_omega= ',domega 223 write(std_out,*) 'number of bands= ', elph_ds%nband, nband 224 225 ABI_ALLOCATE(tmp_wtk,(nband,nkpt,nsppol,nene_all)) 226 ABI_ALLOCATE(tmp2_wtk,(nene_all)) 227 ABI_ALLOCATE(ff2,(nene_all)) 228 ABI_ALLOCATE(ene_pt,(nene_all)) 229 ABI_ALLOCATE(ene_ptfine,(nene_all*nspline)) 230 ABI_ALLOCATE(tmp_wtk1,(nene_all*nspline)) 231 ABI_ALLOCATE(tmp_wtk2,(nene_all*nspline)) 232 ABI_ALLOCATE(dos_e,(nsppol,nene_all)) 233 234 !Get energy points for spline interpolation 235 do iene = 1, nene_all 236 ene_pt(iene) = enemin + (iene-1)*nspline*deltaene 237 end do 238 239 do iene = 1, nene_all*nspline 240 ene_ptfine(iene) = enemin + (iene-1)*deltaene 241 end do 242 243 ABI_ALLOCATE(tmp_wtq,(elph_ds%nbranch, elph_ds%k_phon%nkpt, elph_ds%na2f+1)) 244 ABI_ALLOCATE(wtq,(elph_ds%nbranch, elph_ds%k_phon%nkpt, elph_ds%na2f)) 245 ABI_ALLOCATE(tmp2_wtq,(elph_ds%nbranch, elph_ds%na2f)) 246 247 !phonon 248 ABI_ALLOCATE(phfrq,(nbranch, nkptirr)) 249 ABI_ALLOCATE(displ,(2, nbranch, nbranch, nkptirr)) 250 ABI_ALLOCATE(pheigvec,(2*nbranch*nbranch, nkptirr)) 251 252 do iFSqpt = 1, nkptirr 253 call ifc_fourq(ifc,cryst,elph_ds%k_phon%kptirr(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt),out_eigvec=pheigvec(:,iFSqpt)) 254 end do 255 256 omega_min = omega_min - domega 257 258 !bxu, obtain wtq for the q_fine, then condense to q_phon 259 call ep_ph_weights(phfrq,elph_ds%a2fsmear,omega_min,omega_max,elph_ds%na2f+1,Cryst%gprimd,elph_ds%kptrlatt, & 260 & elph_ds%nbranch,elph_ds%telphint,elph_ds%k_phon,tmp_wtq) 261 omega_min = omega_min + domega 262 263 do iomega = 1, elph_ds%na2f 264 wtq(:,:,iomega) = tmp_wtq(:,:,iomega+1) 265 !write(1005,*) omega_min+(iomega-1)*domega, sum(tmp_wtq(:,:,iomega+1))/nkpt 266 end do 267 ABI_DEALLOCATE(tmp_wtq) 268 269 ! electron 270 tmp_wtk =zero 271 dos_e = zero 272 call ep_el_weights(elph_ds%ep_b_min, elph_ds%ep_b_max, eigenGS(elph_ds%minFSband:elph_ds%minFSband+nband-1,:,:), & 273 & elph_ds%elphsmear, & 274 & enemin, enemax, nene_all, Cryst%gprimd, elph_ds%k_phon%irredtoGS, elph_ds%kptrlatt, max_occ, & 275 & 1, nband, elph_ds%nFSband, nsppol, elph_ds%telphint, elph_ds%k_phon, tmp_wtk) 276 !& elph_ds%minFSband, elph_ds%nband, elph_ds%nFSband, nsppol, elph_ds%telphint, elph_ds%k_phon, tmp_wtk) 277 278 do isppol = 1, nsppol 279 do iene = 1, nene_all 280 dos_e(isppol,iene) = sum(tmp_wtk(:,:,isppol,iene))/nkpt 281 end do 282 end do 283 284 ABI_ALLOCATE(coskr1, (nqpt,nrpt)) 285 ABI_ALLOCATE(sinkr1, (nqpt,nrpt)) 286 call ftgam_init(ifc%gprim, nqpt, nrpt, elph_ds%k_phon%kpt, Ifc%rpt, coskr1, sinkr1) 287 ABI_ALLOCATE(coskr2, (nkptirr,nrpt)) 288 ABI_ALLOCATE(sinkr2, (nkptirr,nrpt)) 289 call ftgam_init(ifc%gprim, nkptirr, nrpt, elph_ds%k_phon%kpt, Ifc%rpt, coskr2, sinkr2) 290 291 !get fermie for itemp 292 fermie = elph_ds%fermie 293 do itemp=1,ntemper ! runs over termperature in K 294 Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp) 295 296 Bst%occopt = 3 297 Bst%tsmear = Temp*kb_HaK 298 call ebands_update_occ(Bst,-99.99_dp) 299 write(message,'(a,f12.6,a,E20.12)')'At T=',Temp,' Fermi level is:',Bst%fermie 300 call wrtout(std_out,message,'COLL') 301 302 if (abs(elph_ds%fermie) < tol10) then 303 fermie(itemp) = Bst%fermie 304 end if 305 end do 306 307 inv_tau_k = zero 308 !get a2f_2d = \sum_{q,nbranch,jband'} |gkk|^2*\delta(\epsilon_{k'j'}-\epsilon')*\delta(\omega_q-\omega) 309 do isppol=1,nsppol 310 write (std_out,*) '##############################################' 311 write (std_out,*) 'get_tau_k : Treating spin polarization ', isppol 312 write (std_out,*) '##############################################' 313 314 ! do iFSkpt =1,nkpt 315 do ik_this_proc =1,elph_ds%k_phon%my_nkpt 316 iFSkpt = elph_ds%k_phon%my_ikpt(ik_this_proc) 317 write (std_out,*) 'get_tau_k : working on kpt # ', iFSkpt, '/', nkpt 318 do jband = 1, nband 319 ! write(*,*)'i am here 1 ', isppol,iFSkpt,jband 320 a2f_2d = zero 321 a2f_2d2 = zero 322 323 !sum from here 324 nskip = 0 325 do jpband = 1, nband 326 jbeff = jpband+(jband-1)*nband 327 328 if (elph_ds%gkqwrite == 0) then 329 tmp_gkk_qpt(:,:,:) = elph_ds%gkk_qpt(:,jbeff,:,ik_this_proc,isppol,:) 330 else if (elph_ds%gkqwrite == 1) then 331 irec = (ik_this_proc-1)*elph_ds%k_phon%my_nkpt + iqpt 332 if (iFSkpt == 1) then 333 write (std_out,*) ' get_tau_k read record ', irec 334 end if 335 read (elph_ds%unitgkq,REC=irec) tmp_gkk_qpt(:,:,iqpt_fullbz) 336 end if 337 338 !FT to real space 339 call ftgam(Ifc%wghatm,tmp_gkk_qpt,tmp_gkk_rpt,natom,nqpt,nrpt,1,coskr1,sinkr1) 340 341 !sum over irred q over k_phon, with corresponding weights 342 do iFSqpt = 1, nkptirr 343 iqpt_fullbz = elph_ds%k_phon%irredtoGS(iFSqpt) 344 ikpt_kpq = FSfullpktofull(iFSkpt,iqpt_fullbz) 345 346 imqpt_fullbz = mqtofull(iqpt_fullbz) 347 ikpt_kmq = FSfullpktofull(iFSkpt,imqpt_fullbz) 348 349 !Do FT from real-space gamma grid to 1 kpt in k_phon%new_kptirr 350 call ftgam(Ifc%wghatm,tmp_gkk_kpt,tmp_gkk_rpt,natom,1,nrpt,0,coskr2(iqpt_fullbz,:),sinkr2(iqpt_fullbz,:)) 351 !tmp_gkk_kpt(:,:)=tmp_gkk_qpt(:,:,iFSqpt) 352 353 !if ep_scalprod==0 we have to dot in the displacement vectors here 354 if (elph_ds%ep_scalprod==0) then 355 356 call phdispl_cart2red(natom,Cryst%gprimd,displ(:,:,:,iFSqpt),displ_red) 357 358 tmp_gkk_kpt2 = reshape (tmp_gkk_kpt(:,:), (/2,nbranch,nbranch/)) 359 call gam_mult_displ(nbranch, displ_red, tmp_gkk_kpt2, gkk_kpt) 360 361 do jbranch=1,nbranch 362 eigval(jbranch) = gkk_kpt(1, jbranch, jbranch) 363 imeigval(jbranch) = gkk_kpt(2, jbranch, jbranch) 364 365 if (abs(imeigval(jbranch)) > tol10) then 366 write (message,'(a,i0,a,es16.8)')" real values branch = ",jbranch,' eigval = ',eigval(jbranch) 367 MSG_WARNING(message) 368 write (message,'(a,i0,a,es16.8)')" imaginary values branch = ",jbranch,' imeigval = ',imeigval(jbranch) 369 MSG_WARNING(message) 370 end if 371 372 end do 373 374 ! if ep_scalprod==1 we have to diagonalize the matrix we interpolated. 375 else if (elph_ds%ep_scalprod == 1) then 376 377 ! MJV NOTE : gam_now is being recast as a (3*natom)**2 matrix here 378 call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, cone, tmp_gkk_kpt, 3*natom,& 379 & pheigvec(:,iFSqpt), 3*natom, czero, tmp_gkk_kpt2, 3*natom) 380 381 call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, cone, pheigvec(:,iFSqpt), 3*natom,& 382 & tmp_gkk_kpt2, 3*natom, czero, gkk_kpt, 3*natom) 383 384 diagerr = zero 385 do ibranch=1,nbranch 386 eigval(ibranch) = gkk_kpt(1,ibranch,ibranch) 387 do jbranch=1,ibranch-1 388 diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch)) 389 end do 390 do jbranch=ibranch+1,nbranch 391 diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch)) 392 end do 393 end do 394 395 if (diagerr > tol12) then 396 write(message,'(a,es15.8)') 'get_tau_k: residual in diagonalization of gamma with phon eigenvectors: ', diagerr 397 MSG_WARNING(message) 398 end if 399 400 else 401 write (message,'(a,i0)')' Wrong value for ep_scalprod = ',elph_ds%ep_scalprod 402 MSG_BUG(message) 403 end if ! end ep_scalprod if 404 405 !For k'=k-q 406 !Do FT from real-space gamma grid to 1 kpt in k_phon%new_kptirr 407 call ftgam(Ifc%wghatm,tmp_gkk_kpt,tmp_gkk_rpt,natom,1,nrpt,0,coskr2(imqpt_fullbz,:),sinkr2(imqpt_fullbz,:)) 408 !tmp_gkk_kpt(:,:)=tmp_gkk_qpt(:,:,iFSqpt) 409 410 !if ep_scalprod==0 we have to dot in the displacement vectors here 411 if (elph_ds%ep_scalprod==0) then 412 413 call phdispl_cart2red(natom,Cryst%gprimd,displ(:,:,:,iFSqpt),displ_red) 414 415 tmp_gkk_kpt2 = reshape (tmp_gkk_kpt(:,:), (/2,nbranch,nbranch/)) 416 call gam_mult_displ(nbranch, displ_red, tmp_gkk_kpt2, gkk_kpt) 417 418 do jbranch=1,nbranch 419 eigval2(jbranch) = gkk_kpt(1, jbranch, jbranch) 420 imeigval(jbranch) = gkk_kpt(2, jbranch, jbranch) 421 422 if (abs(imeigval(jbranch)) > tol10) then 423 write (message,'(a,i0,a,es16.8)')" real values branch = ",jbranch,' eigval = ',eigval2(jbranch) 424 MSG_WARNING(message) 425 write (message,'(a,i0,a,es16.8)')" imaginary values branch = ",jbranch,' imeigval = ',imeigval(jbranch) 426 MSG_WARNING(message) 427 end if 428 429 end do 430 431 ! if ep_scalprod==1 we have to diagonalize the matrix we interpolated. 432 else if (elph_ds%ep_scalprod == 1) then 433 434 ! MJV NOTE : gam_now is being recast as a (3*natom)**2 matrix here 435 call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, cone, tmp_gkk_kpt, 3*natom,& 436 & pheigvec(:,iFSqpt), 3*natom, czero, tmp_gkk_kpt2, 3*natom) 437 438 call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, cone, pheigvec(:,iFSqpt), 3*natom,& 439 & tmp_gkk_kpt2, 3*natom, czero, gkk_kpt, 3*natom) 440 441 diagerr = zero 442 do ibranch=1,nbranch 443 eigval2(ibranch) = gkk_kpt(1,ibranch,ibranch) 444 do jbranch=1,ibranch-1 445 diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch)) 446 end do 447 do jbranch=ibranch+1,nbranch 448 diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch)) 449 end do 450 end do 451 452 if (diagerr > tol12) then 453 write(message,'(a,es15.8)') 'get_tau_k: residual in diagonalization of gamma with phon eigenvectors: ', diagerr 454 MSG_WARNING(message) 455 end if 456 457 else 458 write (message,'(a,i0)')' Wrong value for ep_scalprod = ',elph_ds%ep_scalprod 459 MSG_BUG(message) 460 end if ! end ep_scalprod if 461 462 tmp2_wtk(:) = tmp_wtk(jpband,ikpt_kpq,isppol,:) 463 yp1 = (tmp2_wtk(2)-tmp2_wtk(1))/nspline/deltaene 464 ypn = (tmp2_wtk(nene_all)-tmp2_wtk(nene_all-1))/nspline/deltaene 465 call spline(ene_pt,tmp2_wtk,nene_all,yp1,ypn,ff2) 466 call splint(nene_all,ene_pt,tmp2_wtk,ff2,nene_all*nspline,ene_ptfine,tmp_wtk1) 467 468 tmp2_wtk(:) = tmp_wtk(jpband,ikpt_kmq,isppol,:) 469 yp1 = (tmp2_wtk(2)-tmp2_wtk(1))/nspline/deltaene 470 ypn = (tmp2_wtk(nene_all)-tmp2_wtk(nene_all-1))/nspline/deltaene 471 call spline(ene_pt,tmp2_wtk,nene_all,yp1,ypn,ff2) 472 call splint(nene_all,ene_pt,tmp2_wtk,ff2,nene_all*nspline,ene_ptfine,tmp_wtk2) 473 474 tmp2_wtq(:,:) = wtq(:,iFSqpt,:) 475 do iene=1,nene 476 e_k = eigenGS(elph_ds%minFSband+jband-1,iFSkpt,isppol) 477 ene = e_k - omega_max + (iene-1)*deltaene 478 if (ene<enemin .or. ene>enemax) cycle 479 iene_fine = NINT((ene-enemin+deltaene)/deltaene) 480 tmp_wtkpq = tmp_wtk1(iene_fine) * elph_ds%k_phon%wtkirr(iFSqpt) 481 tmp_wtkmq = tmp_wtk2(iene_fine) * elph_ds%k_phon%wtkirr(iFSqpt) 482 483 if (tmp_wtkpq+tmp_wtkmq < tol_wtk ) then 484 nskip = nskip +1 485 cycle 486 end if 487 488 do ibranch = 1, nbranch 489 if (abs(phfrq(ibranch,iFSqpt)) < tol7) cycle 490 491 if (ene > e_k) then 492 omega = ene - e_k 493 if (abs(omega) < tol7 .or. abs(omega) > omega_max) cycle 494 iomega = NINT((omega-omega_min+domega)/domega) 495 496 a2f_2d(iene) = a2f_2d(iene) +& 497 & eigval(ibranch)/phfrq(ibranch,iFSqpt)*& 498 & tmp_wtkpq * tmp2_wtq(ibranch,iomega) 499 end if 500 501 if (ene < e_k) then 502 omega = e_k - ene 503 if (abs(omega) < tol7 .or. abs(omega) > omega_max) cycle 504 iomega = NINT((omega-omega_min+domega)/domega) 505 506 a2f_2d2(iene) = a2f_2d2(iene) +& 507 & eigval(ibranch)/phfrq(ibranch,iFSqpt)*& 508 & tmp_wtkmq * tmp2_wtq(ibranch,iomega) 509 end if 510 511 end do ! ibranch 3 512 end do ! nene 800 513 end do ! kptirr 216 514 end do ! j' band 3 515 ! print *, ' skipped ', nskip, ' energy points out of ', nene*nband*nkptirr 516 517 ! get inv_tau_k 518 do itemp=1,ntemper ! runs over termperature in K 519 Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp) 520 do iene=1,nene 521 e_k = eigenGS(elph_ds%minFSband+jband-1,iFSkpt,isppol) 522 ene = e_k - omega_max + (iene-1)*deltaene 523 if (ene<enemin .or. ene>enemax) cycle 524 525 xx=(ene-fermie(itemp))/(kb_HaK*Temp) 526 occ_e=1.0_dp/(exp(xx)+1.0_dp) 527 if (ene > e_k .and. (ene-e_k) .le. omega_max) then 528 omega = ene - e_k 529 if (abs(omega) < tol7) cycle 530 xx = omega/(kb_HaK*Temp) 531 occ_omega=1.0_dp/(exp(xx)-1.0_dp) 532 533 therm_factor = occ_e + occ_omega 534 535 inv_tau_k(itemp,isppol,iFSkpt,jband) = inv_tau_k(itemp,isppol,iFSkpt,jband) +& 536 a2f_2d(iene)*therm_factor*deltaene 537 end if 538 if (ene < e_k .and. (e_k-ene) .le. omega_max) then 539 omega = e_k - ene 540 if (abs(omega) < tol7) cycle 541 xx = omega/(kb_HaK*Temp) 542 occ_omega=1.0_dp/(exp(xx)-1.0_dp) 543 544 therm_factor = 1 - occ_e + occ_omega 545 546 inv_tau_k(itemp,isppol,iFSkpt,jband) = inv_tau_k(itemp,isppol,iFSkpt,jband) +& 547 a2f_2d2(iene)*therm_factor*deltaene 548 end if 549 550 end do ! nene 551 end do ! Temp 552 ! write(*,*)'i am here 2 ', isppol,iFSkpt,jband 553 end do ! jband 554 end do ! kpt 555 end do ! nsppol 556 557 !write (300+mpi_enreg%me,*) inv_tau_k 558 call xmpi_sum (inv_tau_k, xmpi_world, ierr) 559 560 ABI_DEALLOCATE(phfrq) 561 ABI_DEALLOCATE(displ) 562 ABI_DEALLOCATE(pheigvec) 563 ABI_DEALLOCATE(tmp2_wtk) 564 ABI_DEALLOCATE(ff2) 565 ABI_DEALLOCATE(ene_pt) 566 ABI_DEALLOCATE(ene_ptfine) 567 ABI_DEALLOCATE(tmp_wtk1) 568 ABI_DEALLOCATE(tmp_wtk2) 569 ABI_DEALLOCATE(tmp2_wtq) 570 ABI_DEALLOCATE(wtq) 571 ABI_DEALLOCATE(coskr1) 572 ABI_DEALLOCATE(sinkr1) 573 ABI_DEALLOCATE(coskr2) 574 ABI_DEALLOCATE(sinkr2) 575 ABI_DEALLOCATE(kpttokpt) 576 ABI_DEALLOCATE(FSfullpktofull) 577 ABI_DEALLOCATE(mqtofull) 578 ABI_DEALLOCATE(tmp_gkk_qpt) 579 ABI_DEALLOCATE(tmp_gkk_rpt) 580 ABI_DEALLOCATE(tmp_gkk_kpt) 581 ABI_DEALLOCATE(tmp_gkk_kpt2) 582 ABI_DEALLOCATE(gkk_kpt) 583 ABI_DEALLOCATE(a2f_2d) 584 ABI_DEALLOCATE(a2f_2d2) 585 586 !output inv_tau_k and tau_k 587 fname = trim(elph_ds%elph_base_name) // '_INVTAUK' 588 if (open_file(fname,message,newunit=unit_invtau,status='unknown') /= 0) then 589 MSG_ERROR(message) 590 end if 591 592 !print header to relaxation time file 593 write (unit_invtau,*) '# k-dep inverse of the relaxation time as a function of temperature.' 594 write (unit_invtau,*) '# ' 595 write (unit_invtau,*) '# nkptirr= ', nkptirr, 'nband= ', nband 596 write (unit_invtau,*) '# number of temperatures= ', ntemper 597 write (unit_invtau,*) '# tau [femtosecond^-1] ' 598 599 fname = trim(elph_ds%elph_base_name) // '_TAUK' 600 if (open_file(fname,message,newunit=unit_tau,status='unknown') /= 0) then 601 MSG_ERROR(message) 602 end if 603 604 !print header to relaxation time file 605 write (unit_tau,*) '# k-dep relaxation time as a function of temperature.' 606 write (unit_tau,*) '# ' 607 write (unit_tau,*) '# nkptirr= ', nkptirr, 'nband= ', nband 608 write (unit_tau,*) '# number of temperatures= ', ntemper 609 write (unit_tau,*) '# tau [femtosecond] ' 610 611 tau_k = zero 612 do itemp=1,ntemper ! runs over termperature in K 613 Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp) 614 write(unit_invtau,'(a,f16.8)') '# Temperature = ', Temp 615 write(unit_tau,'(a,f16.8)') '# Temperature = ', Temp 616 do isppol=1,nsppol 617 write(unit_invtau,'(a,i6)') '# For isppol = ', isppol 618 write(unit_tau,'(a,i6)') '# For isppol = ', isppol 619 do iFSkpt = 1,nkpt 620 !FIXME: check when tau_k is too small, whether there should be a phonon 621 !scattering or not, and should tau_k be zero or not. 622 do jband = 1,nband 623 if (abs(inv_tau_k(itemp,isppol,iFSkpt,jband)) < tol9) then 624 inv_tau_k(itemp,isppol,iFSkpt,jband) = zero 625 tau_k(itemp,isppol,iFSkpt,jband) = zero 626 else 627 !no need to *nkpt due to wtkirr, as we need /nkpt for the sum 628 !no need to *two_pi due to the missing prefactor in gkk (see mka2f_tr_lova) 629 inv_tau_k(itemp,isppol,iFSkpt,jband) = inv_tau_k(itemp,isppol,iFSkpt,jband)*elph_ds%occ_factor 630 tau_k(itemp,isppol,iFSkpt,jband) = one/inv_tau_k(itemp,isppol,iFSkpt,jband) 631 end if 632 end do ! nband 633 write(unit_invtau,'(a,i8,a,3f12.6)') '# kpt# ', iFSkpt, ' kpt=', elph_ds%k_phon%kptirr(:,iFSkpt) 634 write(unit_invtau,'(100D16.8)') (inv_tau_k(itemp,isppol,iFSkpt,iband)*femto/chu_tau,iband=1,nband) 635 write(unit_tau,'(a,i8,a,3f12.6)') '# kpt# ', iFSkpt, ' kpt=', elph_ds%k_phon%kptirr(:,iFSkpt) 636 write(unit_tau,'(100D16.8)') (tau_k(itemp,isppol,iFSkpt,iband)*chu_tau/femto,iband=1,nband) 637 end do ! nkptirr 638 write(unit_invtau,*) ' ' 639 write(unit_tau,*) ' ' 640 end do ! nsppol 641 write(unit_invtau,*) ' ' 642 write(unit_invtau,*) ' ' 643 write(unit_tau,*) ' ' 644 write(unit_tau,*) ' ' 645 end do ! ntemper 646 647 ! Only use the irred k for eigenGS and tau_k 648 ABI_ALLOCATE(tmp_eigenGS,(elph_ds%nband,elph_ds%k_phon%new_nkptirr,elph_ds%nsppol)) 649 650 do ikpt_irr = 1, new_nkptirr 651 tmp_eigenGS(:,ikpt_irr,:) = eigenGS(:,elph_ds%k_phon%new_irredtoGS(ikpt_irr),:) 652 tmp_tau_k(:,:,ikpt_irr,:) = tau_k(:,:,elph_ds%k_phon%new_irredtoGS(ikpt_irr),:)*chu_tau 653 end do 654 655 !BoltzTraP output files in SIESTA format 656 if (elph_ds%prtbltztrp == 1) then 657 call ebands_prtbltztrp_tau_out (tmp_eigenGS(elph_ds%minFSband:elph_ds%maxFSband,:,:),& 658 & elph_ds%tempermin,elph_ds%temperinc,ntemper,fermie, & 659 & elph_ds%elph_base_name,elph_ds%k_phon%new_kptirr,nband,elph_ds%nelect,new_nkptirr, & 660 & elph_ds%nspinor,nsppol,Cryst%nsym,Cryst%rprimd,Cryst%symrel,tmp_tau_k) 661 end if !prtbltztrp 662 ABI_DEALLOCATE(tmp_eigenGS) 663 ABI_DEALLOCATE(tmp_tau_k) 664 665 !Get the energy dependence of tau. 666 !Eq. (6) in Restrepo et al. Appl. Phys. Lett. 94, 212103 (2009) 667 668 fname = trim(elph_ds%elph_base_name) // '_TAUE' 669 if (open_file(fname,message,newunit=unit_taue,status='unknown') /= 0) then 670 MSG_ERROR(message) 671 end if 672 673 !print header to relaxation time file 674 write (unit_taue,*) '# Energy-dep relaxation time as a function of temperature.' 675 write (unit_taue,*) '# ' 676 write (unit_taue,*) '# number of temperatures= ', ntemper 677 write (unit_taue,*) '# ene[Ha] tau [femtosecond] DOS[au] ' 678 679 fname = trim(elph_ds%elph_base_name) // '_MFP' 680 if (open_file(fname,message,newunit=unit_mfp,status='unknown') /= 0) then 681 MSG_ERROR(message) 682 end if 683 684 write (unit_mfp,*) '# Energy-dep mean free path as a function of temperature.' 685 write (unit_mfp,*) '# ' 686 write (unit_mfp,*) '# number of temperatures= ', ntemper 687 write (unit_mfp,*) '# ene[Ha] mfp [femtometer] ' 688 689 do itemp=1,ntemper ! runs over termperature in K 690 Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp) 691 write(unit_taue,'(a,f16.8)') '# Temperature = ', Temp 692 do isppol = 1, nsppol 693 write(unit_taue,*) '# Tau_e for isppol = ',isppol 694 do iene = 1, nene_all 695 rate_e = zero 696 do iFSkpt = 1, nkpt 697 do jband = 1, nband 698 rate_e = rate_e + inv_tau_k(itemp,isppol,iFSkpt,jband)* & 699 & tmp_wtk(jband,iFSkpt,isppol,iene) 700 end do ! jband 701 end do ! kpt 702 if (dabs(dos_e(isppol,iene)) < tol7) then 703 rate_e = zero 704 else 705 rate_e = rate_e/nkpt/dos_e(isppol,iene) 706 end if 707 write(unit_taue,"(3D16.8)") enemin+(iene-1)*deltaene*nspline, rate_e*femto/chu_tau, dos_e(isppol,iene) 708 end do ! number of energies 709 write(unit_taue,*) ' ' 710 end do ! nsppol 711 write(unit_taue,*) ' ' 712 end do ! ntemperature 713 714 ! calculate and output mean free path 715 do itemp=1,ntemper ! runs over termperature in K 716 Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp) 717 write(unit_mfp,'(a,f16.8)') '# Temperature = ', Temp 718 do isppol = 1, nsppol 719 do icomp = 1, 3 720 write(unit_mfp,*) '# Mean free path for isppol, icomp= ',isppol,icomp 721 do iene = 1, nene_all 722 mfp_e = zero 723 do iFSkpt = 1, nkptirr 724 do jband = 1, nband 725 mfp_e = mfp_e + tau_k(itemp,isppol,iFSkpt,jband)* & 726 & elph_tr_ds%el_veloc(iFSkpt,elph_ds%minFSband+jband-1,icomp,isppol)* & 727 & tmp_wtk(jband,iFSkpt,isppol,iene) 728 !& elph_ds%k_phon%new_wtkirr(iFSqpt) 729 end do ! jband 730 end do ! kpt 731 if (dabs(dos_e(isppol,iene)) < tol7) then 732 mfp_e = zero 733 else 734 mfp_e = mfp_e/nkptirr/dos_e(isppol,iene) 735 end if 736 write(unit_mfp,"(2D16.8)") enemin+(iene-1)*deltaene*nspline, mfp_e*chu_mfp/femto 737 end do ! number of energies 738 write(unit_mfp,*) ' ' 739 end do ! icomp 740 write(unit_mfp,*) ' ' 741 end do ! nsppol 742 write(unit_mfp,*) ' ' 743 end do ! ntemperature 744 745 ABI_ALLOCATE(cond_e ,(ntemper,nsppol,nene_all,9)) 746 747 !get cond_e 748 cond_e = zero 749 do itemp=1,ntemper ! runs over termperature in K 750 do isppol = 1, nsppol 751 do iene = 1, nene_all 752 ! do iFSkpt =1,nkpt 753 do ik_this_proc =1,elph_ds%k_phon%my_nkpt 754 iFSkpt = elph_ds%k_phon%my_ikpt(ik_this_proc) 755 do jband = 1, nband 756 do icomp = 1, 3 757 do jcomp = 1, 3 758 itensor = (icomp-1)*3+jcomp 759 cond_e(itemp,isppol,iene,itensor) = cond_e(itemp,isppol,iene,itensor) + & 760 & tau_k(itemp,isppol,iFSkpt,jband)* & 761 & elph_tr_ds%el_veloc(iFSkpt,elph_ds%minFSband+jband-1,icomp,isppol)* & 762 & elph_tr_ds%el_veloc(iFSkpt,elph_ds%minFSband+jband-1,jcomp,isppol)* & 763 & tmp_wtk(jband,iFSkpt,isppol,iene) 764 end do 765 end do 766 end do ! jband 767 end do ! kpt 768 end do ! number of energies 769 end do ! nsppol 770 end do ! ntemperature 771 772 ! MG FIXME: Why xmpi_world, besides only master should perform IO in the section below. 773 call xmpi_sum (cond_e, xmpi_world, ierr) 774 775 cond_e = cond_e/nkpt 776 777 !get transport coefficients 778 779 fname = trim(elph_ds%elph_base_name) // '_COND' 780 if (open_file(fname,message,newunit=unit_cond,status='unknown') /= 0) then 781 MSG_ERROR(message) 782 end if 783 784 !print header to conductivity file 785 write (unit_cond,*) '# Conductivity as a function of temperature.' 786 write (unit_cond,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 787 write (unit_cond,*) '# ' 788 write (unit_cond,*) '# Columns are: ' 789 write (unit_cond,*) '# temperature[K] cond[au] cond [SI] ' 790 write (unit_cond,*) '# ' 791 792 fname = trim(elph_ds%elph_base_name) // '_CTH' 793 if (open_file(fname,message,newunit=unit_therm,status='unknown') /= 0) then 794 MSG_ERROR(message) 795 end if 796 797 !print header to thermal conductivity file 798 write (unit_therm,'(a)') '# Thermal conductivity as a function of temperature.' 799 write (unit_therm,'(a)') '# the formalism is isotropic, so non-cubic crystals may be wrong' 800 write (unit_therm,'(a)') '# ' 801 write (unit_therm,'(a)') '# Columns are: ' 802 write (unit_therm,'(a)') '# temperature[K] thermal cond [au] thermal cond [SI]' 803 write (unit_therm,'(a)') '# ' 804 805 fname = trim(elph_ds%elph_base_name) // '_SBK' 806 if (open_file(fname,message,newunit=unit_sbk,status='unknown') /=0) then 807 MSG_ERROR(message) 808 end if 809 810 !print header to relaxation time file 811 write (unit_sbk,*) '# Seebeck Coefficint as a function of temperature.' 812 write (unit_sbk,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 813 write (unit_sbk,*) '# ' 814 write (unit_sbk,*) '# Columns are: ' 815 write (unit_sbk,*) '# temperature[K] S [au] S [SI] ' 816 write (unit_sbk,*) '# ' 817 818 ABI_ALLOCATE(cond ,(ntemper,nsppol,3,3)) 819 ABI_ALLOCATE(cth ,(ntemper,nsppol,3,3)) 820 ABI_ALLOCATE(sbk ,(ntemper,nsppol,3,3)) 821 ABI_ALLOCATE(seebeck ,(ntemper,nsppol,3,3)) 822 823 cond = zero 824 cth = zero 825 sbk = zero 826 seebeck = zero 827 do isppol=1,nsppol 828 do icomp=1, 3 829 do jcomp=1, 3 830 itensor=(icomp-1)*3+jcomp 831 do itemp=1,ntemper 832 Temp=elph_ds%tempermin + elph_ds%temperinc*dble(itemp) 833 do iene = 1, nene_all 834 factor = (enemin+(iene-1)*deltaene*nspline - fermie(itemp))/(kb_HaK*Temp) 835 if (factor < -40.0d0) then 836 dfermide = zero 837 else if (factor > 40.0d0) then 838 dfermide = zero 839 else 840 dfermide = EXP(factor)/(kb_HaK*Temp*(EXP(factor)+one)**2.0d0) 841 end if 842 cond(itemp,isppol,icomp,jcomp) = cond(itemp,isppol,icomp,jcomp) + & 843 & cond_e(itemp,isppol,iene,itensor)*dfermide*deltaene*nspline 844 cth(itemp,isppol,icomp,jcomp) = cth(itemp,isppol,icomp,jcomp) + cond_e(itemp,isppol,iene,itensor)* & 845 & (enemin+(iene-1)*deltaene*nspline - fermie(itemp))**2.0d0*dfermide*deltaene*nspline 846 sbk(itemp,isppol,icomp,jcomp) = sbk(itemp,isppol,icomp,jcomp) + cond_e(itemp,isppol,iene,itensor)* & 847 & (enemin+(iene-1)*deltaene*nspline - fermie(itemp))*dfermide*deltaene*nspline 848 end do 849 end do ! temperature 850 end do ! jcomp 851 end do ! icomp 852 end do !end isppol 853 854 do isppol=1,nsppol 855 do itemp=1,ntemper 856 cond_inv(:,:)=cond(itemp,isppol,:,:) 857 call matrginv(cond_inv,3,3) 858 call DGEMM('N','N',3,3,3,one,sbk(itemp,isppol,:,:),3,cond_inv,& 859 & 3,zero,seebeck(itemp,isppol,:,:),3) 860 end do 861 end do 862 863 do isppol=1,nsppol 864 do icomp=1, 3 865 do jcomp=1, 3 866 itensor=(icomp-1)*3+jcomp 867 write(unit_cond,*) '# Conductivity for isppol, itrten= ',isppol,itensor 868 write(unit_therm,*) '# Thermal conductivity for isppol, itrten= ',isppol,itensor 869 write(unit_sbk,*) '# Seebeck coefficient for isppol, itrten= ',isppol,itensor 870 do itemp=1,ntemper 871 Temp=elph_ds%tempermin + elph_ds%temperinc*dble(itemp) 872 873 seebeck(itemp,isppol,icomp,jcomp) = -1.0d0*seebeck(itemp,isppol,icomp,jcomp)/(kb_HaK*Temp) 874 cond(itemp,isppol,icomp,jcomp) = cond(itemp,isppol,icomp,jcomp)/cryst%ucvol 875 cth(itemp,isppol,icomp,jcomp) = cth(itemp,isppol,icomp,jcomp)/(kb_HaK*Temp)/cryst%ucvol 876 write(unit_cond,'(3D20.10)')Temp,cond(itemp,isppol,icomp,jcomp),cond(itemp,isppol,icomp,jcomp)*chu_cond 877 write(unit_therm,'(3D20.10)')Temp,cth(itemp,isppol,icomp,jcomp),cth(itemp,isppol,icomp,jcomp)*chu_cth 878 write(unit_sbk,'(3D20.10)')Temp,seebeck(itemp,isppol,icomp,jcomp),seebeck(itemp,isppol,icomp,jcomp)*chu_sbk 879 end do ! temperature 880 write(unit_cond,*) 881 write(unit_therm,*) 882 write(unit_sbk,*) 883 end do ! jcomp 884 end do ! icomp 885 end do !end isppol 886 887 888 ABI_DEALLOCATE(inv_tau_k) 889 ABI_DEALLOCATE(tau_k) 890 ABI_DEALLOCATE(tmp_wtk) 891 ABI_DEALLOCATE(dos_e) 892 ABI_DEALLOCATE(cond_e) 893 ABI_DEALLOCATE(fermie) 894 ABI_DEALLOCATE(cond) 895 ABI_DEALLOCATE(sbk) 896 ABI_DEALLOCATE(cth) 897 ABI_DEALLOCATE(seebeck) 898 899 close (unit=unit_tau) 900 close (unit=unit_taue) 901 close (unit=unit_mfp) 902 close (unit=unit_invtau) 903 close (unit=unit_cond) 904 close (unit=unit_therm) 905 close (unit=unit_sbk) 906 907 end subroutine get_tau_k