TABLE OF CONTENTS
ABINIT/mka2f_tr_lova [ Functions ]
NAME
mka2f_tr_lova
FUNCTION
calculates the FS averaged Transport alpha^2F_tr alpha^2F_trout alpha^2F_trin functions calculates and outputs the associated electrical and thermal conductivities for the first task: copied from mka2F
COPYRIGHT
Copyright (C) 2004-2018 ABINIT group (JPC, MJV) 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_fine%nkpt = number of kpts included in the FS integration elph_ds%k_fine%wtk = integration weights on the FS delph_ds%n0 = DOS at the Fermi level calculated from the k_fine integration weights elph_ds%k_fine%kpt = coordinates of all FS kpoints mustar = coulomb pseudopotential parameter eventually for 2 spin channels 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)
OUTPUT
elph_ds
PARENTS
elphon
CHILDREN
ftgam,ftgam_init,gam_mult_displ,ifc_fourq,simpson_int,wrtout,zgemm
NOTES
copied from ftiaf9.f
SOURCE
50 #if defined HAVE_CONFIG_H 51 #include "config.h" 52 #endif 53 54 #include "abi_common.h" 55 56 subroutine mka2f_tr_lova(crystal,ifc,elph_ds,ntemper,tempermin,temperinc,elph_tr_ds) 57 58 use defs_basis 59 use defs_elphon 60 use m_profiling_abi 61 use m_errors 62 63 use m_io_tools, only : open_file 64 use m_numeric_tools, only : simpson_int 65 use m_crystal, only : crystal_t 66 use m_ifc, only : ifc_type, ifc_fourq 67 use m_dynmat, only : ftgam_init, ftgam 68 69 !This section has been created automatically by the script Abilint (TD). 70 !Do not modify the following lines by hand. 71 #undef ABI_FUNC 72 #define ABI_FUNC 'mka2f_tr_lova' 73 use interfaces_14_hidewrite 74 use interfaces_77_ddb, except_this_one => mka2f_tr_lova 75 !End of the abilint section 76 77 implicit none 78 79 !Arguments ------------------------------------ 80 !scalars 81 integer,intent(in) :: ntemper 82 real(dp),intent(in) :: tempermin,temperinc 83 type(crystal_t),intent(in) :: crystal 84 type(ifc_type),intent(in) :: ifc 85 type(elph_tr_type),intent(inout) :: elph_tr_ds 86 type(elph_type),intent(inout) :: elph_ds 87 88 !Local variables ------------------------- 89 !x =w/(2kbT) 90 !scalars 91 integer :: iFSqpt,ibranch,iomega,isppol,jbranch,nerr 92 integer :: unit_a2f_tr, unit_a2f_trout, unit_a2f_trin, natom 93 integer :: idir, iatom, k1, kdir,unit_lor,unit_rho,unit_tau,unit_therm 94 integer :: itemp,nrpt,itrtensor, icomp, jcomp 95 real(dp) :: Temp,chgu,chtu,femto,diagerr,firh,firhT,gaussfactor,domega 96 real(dp) :: firh_tau,firhT_tau ! added by BX to get Tau 97 real(dp) :: a2fprefactor_in, temp_in 98 real(dp) :: a2fprefactor_out, temp_out 99 real(dp) :: gaussprefactor,gaussval,lambda_tr,lor0,lorentz,maxerr,maxx,omega 100 real(dp) :: rho,tau,tolexp,wtherm,xtr,xx 101 real(dp) :: lambda_tr_trace,omega_min, omega_max,qnorm2,spinfact 102 character(len=500) :: message 103 character(len=fnlen) :: fname 104 !arrays 105 real(dp),parameter :: c0(2)=(/0.d0,0.d0/),c1(2)=(/1.d0,0.d0/) 106 real(dp) :: gprimd(3,3) 107 real(dp) :: eigval_in(elph_ds%nbranch) 108 real(dp) :: eigval_out(elph_ds%nbranch) 109 real(dp) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch) 110 real(dp) :: gam_now_in (2,elph_ds%nbranch*elph_ds%nbranch) 111 real(dp) :: gam_now_out(2,elph_ds%nbranch*elph_ds%nbranch) 112 real(dp) :: tmpa2f_in (elph_ds%na2f) 113 real(dp) :: tmpa2f_out(elph_ds%na2f) 114 real(dp) :: tmpgam1(2,elph_ds%nbranch,elph_ds%nbranch) 115 real(dp) :: tmpgam2(2,elph_ds%nbranch,elph_ds%nbranch) 116 real(dp),allocatable :: phfrq(:,:) 117 real(dp),allocatable :: displ(:,:,:,:) 118 real(dp),allocatable :: pheigvec(:,:) 119 real(dp),allocatable :: integrho(:),integtau(:),tointegrho(:),tointega2f(:),tointegtau(:) 120 real(dp),allocatable :: rho_T(:),tau_T(:) 121 real(dp),allocatable :: coskr(:,:) 122 real(dp),allocatable :: sinkr(:,:) 123 124 ! ********************************************************************* 125 126 !calculate a2f_tr for frequencies between 0 and omega_max 127 write(std_out,*) 'mka2f_tr_lova : enter ' 128 ! 129 !MG: the step should be calculated locally using nomega and the extrema of the spectrum. 130 !One should not rely on previous calls for the setup of elph_ds%domega 131 !I will remove elph_ds%domega since mka2f.F90 will become a method of gamma_t 132 domega =elph_ds%domega 133 134 ! Number of points for FFT interpolation 135 nrpt = ifc%nrpt 136 natom = crystal%natom 137 gprimd = crystal%gprimd 138 139 ABI_ALLOCATE(elph_tr_ds%a2f_1d_tr,(elph_ds%na2f,9,elph_ds%nsppol,1,1,1)) 140 ABI_ALLOCATE(elph_tr_ds%a2f_1d_trin,(elph_ds%na2f,9,elph_ds%nsppol)) 141 ABI_ALLOCATE(elph_tr_ds%a2f_1d_trout,(elph_ds%na2f,9,elph_ds%nsppol)) 142 143 !! defaults for number of temperature steps and max T (all in Kelvin...) 144 !ntemper=1000 145 !tempermin=zero 146 !temperinc=one 147 ABI_ALLOCATE(rho_T,(ntemper)) 148 ABI_ALLOCATE(tau_T,(ntemper)) 149 150 151 !tolerance on gaussian being = 0 152 tolexp = 1.d-100 153 maxx = sqrt(-log(tolexp)) 154 lor0=(pi*kb_HaK)**2/3. 155 156 !maximum value of frequency (a grid has to be chosen for the representation of alpha^2 F) 157 !WARNING! supposes this value has been set in mkelph_linwid. 158 159 gaussprefactor = sqrt(piinv) / elph_ds%a2fsmear 160 gaussfactor = one / elph_ds%a2fsmear 161 162 !spinfact should be 1 for a normal non sppol calculation without spinorbit 163 !for spinors it should also be 1 as bands are twice as numerous but n0 has been divided by 2 164 !for sppol 2 it should be 0.5 as we have 2 spin channels to sum 165 spinfact = one / elph_ds%nsppol !/ elph_ds%nspinor 166 167 !ENDMG 168 169 elph_tr_ds%a2f_1d_tr = zero 170 elph_tr_ds%a2f_1d_trin = zero 171 elph_tr_ds%a2f_1d_trout = zero 172 173 maxerr=0. 174 nerr=0 175 176 ABI_ALLOCATE(phfrq,(elph_ds%nbranch, elph_ds%k_fine%nkpt)) 177 ABI_ALLOCATE(displ,(2, elph_ds%nbranch, elph_ds%nbranch, elph_ds%k_fine%nkpt)) 178 ABI_ALLOCATE(pheigvec,(2*elph_ds%nbranch*elph_ds%nbranch, elph_ds%k_fine%nkpt)) 179 180 do iFSqpt=1,elph_ds%k_fine%nkpt 181 call ifc_fourq(ifc,crystal,elph_ds%k_fine%kpt(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt),out_eigvec=pheigvec(:,iFSqpt)) 182 end do 183 184 omega_min = minval(phfrq) 185 omega_max = maxval(phfrq) 186 187 ABI_ALLOCATE(coskr, (elph_ds%k_fine%nkpt,nrpt)) 188 ABI_ALLOCATE(sinkr, (elph_ds%k_fine%nkpt,nrpt)) 189 call ftgam_init(Ifc%gprim, elph_ds%k_fine%nkpt, nrpt, elph_ds%k_fine%kpt, Ifc%rpt, coskr, sinkr) 190 191 do isppol=1,elph_ds%nsppol 192 193 ! loop over qpoint in full kpt grid (presumably dense) 194 do iFSqpt=1,elph_ds%k_fine%nkpt 195 qnorm2 = sum(elph_ds%k_fine%kpt(:,iFSqpt)**2) 196 ! if (flag_to_exclude_soft_modes = .false.) qnorm2 = zero 197 do itrtensor=1,9 198 ! Do FT from real-space gamma grid to 1 qpt. 199 200 if (elph_ds%ep_int_gkk == 1) then 201 gam_now_in(:,:) = elph_tr_ds%gamma_qpt_trin(:,itrtensor,:,isppol,iFSqpt) 202 gam_now_out(:,:) = elph_tr_ds%gamma_qpt_trout(:,itrtensor,:,isppol,iFSqpt) 203 else 204 call ftgam(Ifc%wghatm,gam_now_in, elph_tr_ds%gamma_rpt_trin(:,itrtensor,:,isppol,:),natom,1,nrpt,0,& 205 & coskr(iFSqpt,:), sinkr(iFSqpt,:)) 206 call ftgam(Ifc%wghatm,gam_now_out,elph_tr_ds%gamma_rpt_trout(:,itrtensor,:,isppol,:),natom,1,nrpt,0,& 207 & coskr(iFSqpt,:), sinkr(iFSqpt,:)) 208 end if 209 210 ! Diagonalize gamma matrix at this qpoint (complex matrix). 211 212 ! if ep_scalprod==0 we have to dot in the displacement vectors here 213 if (elph_ds%ep_scalprod==0) then 214 215 displ_red(:,:,:) = zero 216 do jbranch=1,elph_ds%nbranch 217 do iatom=1,natom 218 do idir=1,3 219 ibranch=idir+3*(iatom-1) 220 do kdir=1,3 221 k1 = kdir+3*(iatom-1) 222 displ_red(1,ibranch,jbranch) = displ_red(1,ibranch,jbranch) + & 223 & gprimd(kdir,idir)*displ(1,k1,jbranch,iFSqpt) 224 displ_red(2,ibranch,jbranch) = displ_red(2,ibranch,jbranch) + & 225 & gprimd(kdir,idir)*displ(2,k1,jbranch,iFSqpt) 226 end do 227 end do 228 end do 229 end do 230 231 tmpgam2 = reshape (gam_now_in, (/2,elph_ds%nbranch,elph_ds%nbranch/)) 232 call gam_mult_displ(elph_ds%nbranch, displ_red, tmpgam2, tmpgam1) 233 do jbranch=1,elph_ds%nbranch 234 eigval_in(jbranch) = tmpgam1(1, jbranch, jbranch) 235 end do 236 237 tmpgam2 = reshape (gam_now_out, (/2,elph_ds%nbranch,elph_ds%nbranch/)) 238 call gam_mult_displ(elph_ds%nbranch, displ_red, tmpgam2, tmpgam1) 239 do jbranch=1,elph_ds%nbranch 240 eigval_out(jbranch) = tmpgam1(1, jbranch, jbranch) 241 end do 242 243 else if (elph_ds%ep_scalprod == 1) then 244 245 ! 246 ! NOTE: in these calls gam_now and pheigvec do not have the right rank, but blas usually does not care 247 call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, c1, gam_now_in, 3*natom,& 248 & pheigvec(:,iFSqpt), 3*natom, c0, tmpgam1, 3*natom) 249 call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, c1, pheigvec(:,iFSqpt), 3*natom,& 250 & tmpgam1, 3*natom, c0, tmpgam2, 3*natom) 251 diagerr = zero 252 253 do ibranch=1,elph_ds%nbranch 254 eigval_in(ibranch) = tmpgam2(1,ibranch,ibranch) 255 do jbranch=1,ibranch-1 256 diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch)) 257 end do 258 do jbranch=ibranch+1,elph_ds%nbranch 259 diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch)) 260 end do 261 end do 262 if (diagerr > tol12) then 263 nerr=nerr+1 264 maxerr=max(diagerr, maxerr) 265 end if 266 267 call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, c1, gam_now_out, 3*natom,& 268 & pheigvec(:,iFSqpt), 3*natom, c0, tmpgam1, 3*natom) 269 call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, c1, pheigvec(:,iFSqpt), 3*natom,& 270 & tmpgam1, 3*natom, c0, tmpgam2, 3*natom) 271 diagerr = zero 272 273 do ibranch=1,elph_ds%nbranch 274 eigval_out(ibranch) = tmpgam2(1,ibranch,ibranch) 275 do jbranch=1,ibranch-1 276 diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch)) 277 end do 278 do jbranch=ibranch+1,elph_ds%nbranch 279 diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch)) 280 end do 281 end do 282 if (diagerr > tol12) then 283 nerr=nerr+1 284 maxerr=max(diagerr, maxerr) 285 end if 286 end if 287 ! end ep_scalprod if 288 289 ! Add all contributions from the phonon modes at this qpoint to 290 ! a2f and the phonon dos. 291 do ibranch=1,elph_ds%nbranch 292 ! if (abs(phfrq(ibranch,iFSqpt)) < tol10) then 293 if ( abs(phfrq(ibranch,iFSqpt)) < tol7 .or. & 294 & (phfrq(ibranch,iFSqpt) < tol4 .and. qnorm2 > 0.03 )) then ! 295 ! note: this should depend on the velocity of sound, to accept acoustic 296 ! modes! 297 a2fprefactor_in = zero 298 a2fprefactor_out= zero 299 else 300 a2fprefactor_in = eigval_in (ibranch)/(two_pi*abs(phfrq(ibranch,iFSqpt))*elph_ds%n0(isppol)) 301 a2fprefactor_out = eigval_out(ibranch)/(two_pi*abs(phfrq(ibranch,iFSqpt))*elph_ds%n0(isppol)) 302 end if 303 304 omega = omega_min 305 tmpa2f_in (:) = zero 306 tmpa2f_out(:) = zero 307 do iomega=1,elph_ds%na2f 308 xx = (omega-phfrq(ibranch,iFSqpt))*gaussfactor 309 gaussval = gaussprefactor*exp(-xx*xx) 310 311 temp_in = gaussval*a2fprefactor_in 312 temp_out = gaussval*a2fprefactor_out 313 314 if (dabs(temp_in) < 1.0d-50) temp_in = zero 315 if (dabs(temp_out) < 1.0d-50) temp_out = zero 316 tmpa2f_in (iomega) = tmpa2f_in (iomega) + temp_in 317 tmpa2f_out(iomega) = tmpa2f_out(iomega) + temp_out 318 omega = omega+domega 319 end do 320 321 elph_tr_ds%a2f_1d_trin (:,itrtensor,isppol) = elph_tr_ds%a2f_1d_trin (:,itrtensor,isppol) + tmpa2f_in(:) 322 elph_tr_ds%a2f_1d_trout(:,itrtensor,isppol) = elph_tr_ds%a2f_1d_trout(:,itrtensor,isppol) + tmpa2f_out(:) 323 324 end do ! end ibranch do 325 end do ! end itrtensor do 326 end do ! end iFSqpt do 327 end do ! end isppol 328 329 ABI_DEALLOCATE(coskr) 330 ABI_DEALLOCATE(sinkr) 331 332 !second 1 / elph_ds%k_fine%nkpt factor for the integration weights 333 elph_tr_ds%a2f_1d_trin = elph_tr_ds%a2f_1d_trin / elph_ds%k_fine%nkpt 334 elph_tr_ds%a2f_1d_trout = elph_tr_ds%a2f_1d_trout / elph_ds%k_fine%nkpt 335 336 if (elph_ds%ep_scalprod == 1) then 337 write(std_out,*) 'mka2f_tr_lova: errors in diagonalization of gamma_tr with phon eigenvectors: ', nerr,maxerr 338 end if 339 340 elph_tr_ds%a2f_1d_tr(:,:,:,1,1,1) = elph_tr_ds%a2f_1d_trout(:,:,:) - elph_tr_ds%a2f_1d_trin(:,:,:) 341 342 !output the elph_tr_ds%a2f_1d_tr 343 fname = trim(elph_ds%elph_base_name) // '_A2F_TR' 344 if (open_file (fname,message,newunit=unit_a2f_tr,status='unknown') /= 0) then 345 MSG_ERROR(message) 346 end if 347 348 fname = trim(elph_ds%elph_base_name) // '_A2F_TRIN' 349 if (open_file(fname,message,newunit=unit_a2f_trin,status='unknown') /= 0) then 350 MSG_ERROR(message) 351 end if 352 353 fname = trim(elph_ds%elph_base_name) // '_A2F_TROUT' 354 if (open_file (fname,message,newunit=unit_a2f_trout,status='unknown') /=0) then 355 MSG_ERROR(message) 356 end if 357 358 write (unit_a2f_tr,'(a)') '#' 359 write (unit_a2f_tr,'(a)') '# ABINIT package : a2f_tr file' 360 write (unit_a2f_tr,'(a)') '#' 361 write (unit_a2f_tr,'(a)') '# a2f_tr function integrated over the FS. omega in a.u.' 362 write (unit_a2f_tr,'(a,I10)') '# number of kpoints integrated over : ', elph_ds%k_fine%nkpt 363 write (unit_a2f_tr,'(a,I10)') '# number of energy points : ',elph_ds%na2f 364 write (unit_a2f_tr,'(a,E16.6,a,E16.6,a)') '# between omega_min = ', omega_min,' Ha and omega_max = ', omega_max, ' Ha' 365 write (unit_a2f_tr,'(a,E16.6)') '# and the smearing width for gaussians is ', elph_ds%a2fsmear 366 write (unit_a2f_tr,'(a)') '#' 367 368 write (unit_a2f_trin,'(a)') '#' 369 write (unit_a2f_trin,'(a)') '# ABINIT package : a2f_trin file' 370 write (unit_a2f_trin,'(a)') '#' 371 write (unit_a2f_trin,'(a)') '# a2f_trin function integrated over the FS. omega in a.u.' 372 write (unit_a2f_trin,'(a,I10)') '# number of kpoints integrated over : ', elph_ds%k_fine%nkpt 373 write (unit_a2f_trin,'(a,I10)') '# number of energy points : ',elph_ds%na2f 374 write (unit_a2f_trin,'(a,E16.6,a,E16.6,a)') '# between omega_min = ', omega_min,' Ha and omega_max = ', omega_max, ' Ha' 375 write (unit_a2f_trin,'(a,E16.6)') '# and the smearing width for gaussians is ', elph_ds%a2fsmear 376 write (unit_a2f_trin,'(a)') '#' 377 378 write (unit_a2f_trout,'(a)') '#' 379 write (unit_a2f_trout,'(a)') '# ABINIT package : a2f_trout file' 380 write (unit_a2f_trout,'(a)') '#' 381 write (unit_a2f_trout,'(a)') '# a2f_trout function integrated over the FS. omega in a.u.' 382 write (unit_a2f_trout,'(a,I10)') '# number of kpoints integrated over : ', elph_ds%k_fine%nkpt 383 write (unit_a2f_trout,'(a,I10)') '# number of energy points : ',elph_ds%na2f 384 write (unit_a2f_trout,'(a,E16.6,a,E16.6,a)') '# between omega_min = ', omega_min,' Ha and omega_max = ', omega_max, ' Ha' 385 write (unit_a2f_trout,'(a,E16.6)') '# and the smearing width for gaussians is ', elph_ds%a2fsmear 386 write (unit_a2f_trout,'(a)') '#' 387 388 !done with header 389 do isppol=1,elph_ds%nsppol 390 write (unit_a2f_tr,'(a,E16.6)') '# The DOS at Fermi level is ', elph_ds%n0(isppol) 391 write (unit_a2f_trin,'(a,E16.6)') '# The DOS at Fermi level is ', elph_ds%n0(isppol) 392 write (unit_a2f_trout,'(a,E16.6)') '# The DOS at Fermi level is ', elph_ds%n0(isppol) 393 ! omega = zero 394 omega = omega_min 395 do iomega=1,elph_ds%na2f 396 write (unit_a2f_tr, '(10D16.6)') omega, elph_tr_ds%a2f_1d_tr (iomega,:,isppol,1,1,1) 397 write (unit_a2f_trin, '(10D16.6)') omega, elph_tr_ds%a2f_1d_trin (iomega,:,isppol) 398 write (unit_a2f_trout,'(10D16.6)') omega, elph_tr_ds%a2f_1d_trout(iomega,:,isppol) 399 omega=omega+domega 400 end do 401 write (unit_a2f_tr,*) 402 write (unit_a2f_trin,*) 403 write (unit_a2f_trout,*) 404 end do !isppol 405 406 close (unit=unit_a2f_tr) 407 close (unit=unit_a2f_trin) 408 close (unit=unit_a2f_trout) 409 410 !calculation of transport properties 411 ABI_ALLOCATE(integrho,(elph_ds%na2f)) 412 ABI_ALLOCATE(tointegrho,(elph_ds%na2f)) 413 ABI_ALLOCATE(tointega2f,(elph_ds%na2f)) 414 ABI_ALLOCATE(integtau,(elph_ds%na2f)) 415 ABI_ALLOCATE(tointegtau,(elph_ds%na2f)) 416 417 fname = trim(elph_ds%elph_base_name) // '_RHO' 418 if (open_file(fname,message,newunit=unit_rho,status='unknown') /= 0) then 419 MSG_ERROR(message) 420 end if 421 422 !print header to resistivity file 423 write (unit_rho,*) '# Resistivity as a function of temperature.' 424 write (unit_rho,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 425 write (unit_rho,*) '# ' 426 write (unit_rho,*) '# Columns are: ' 427 write (unit_rho,*) '# temperature[K] rho[au] rho [SI] rho/temp [au]' 428 write (unit_rho,*) '# ' 429 430 fname = trim(elph_ds%elph_base_name) // '_TAU' 431 if (open_file(fname,message,newunit=unit_tau,status='unknown') /= 0) then 432 MSG_ERROR(message) 433 end if 434 435 !print header to relaxation time file 436 write (unit_tau,*) '# Relaxation time as a function of temperature.' 437 write (unit_tau,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 438 write (unit_tau,*) '# ' 439 write (unit_tau,*) '# Columns are: ' 440 write (unit_tau,*) '# temperature[K] tau[au] tau [femtosecond] ' 441 write (unit_tau,*) '# ' 442 443 fname = trim(elph_ds%elph_base_name) // '_WTH' 444 if (open_file(fname,message,newunit=unit_therm,status='unknown') /= 0) then 445 MSG_ERROR(message) 446 end if 447 448 !print header to thermal conductivity file 449 write (unit_therm,'(a)') '# Thermal conductivity/resistivity as a function of temperature.' 450 write (unit_therm,'(a)') '# the formalism is isotropic, so non-cubic crystals may be wrong' 451 write (unit_therm,'(a)') '# ' 452 write (unit_therm,'(a)') '# Columns are: ' 453 write (unit_therm,'(a)') '# temperature[K] thermal rho[au] thermal cond [au] thermal rho [SI] thermal cond [SI]' 454 write (unit_therm,'(a)') '# ' 455 456 fname = trim(elph_ds%elph_base_name) // '_LOR' 457 if (open_file(fname,message,newunit=unit_lor,status='unknown') /= 0) then 458 MSG_ERROR(message) 459 end if 460 461 !print header to lorentz file 462 write (unit_lor,*) '# Lorentz number as a function of temperature.' 463 write (unit_lor,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 464 write (unit_lor,*) '# ' 465 write (unit_lor,*) '# Columns are: ' 466 write (unit_lor,*) '# temperature[K] Lorentz number[au] Lorentz quantum = (pi*kb_HaK)**2/3' 467 write (unit_lor,*) '# ' 468 469 do isppol=1,elph_ds%nsppol 470 lambda_tr_trace = zero 471 do itrtensor=1,9 472 omega = omega_min 473 tointega2f = zero 474 do iomega=1,elph_ds%na2f 475 if(omega<=0) then 476 omega=omega+domega 477 cycle 478 end if 479 tointega2f(iomega)=elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1)/omega 480 omega=omega+domega 481 end do 482 483 integrho = zero 484 call simpson_int(elph_ds%na2f,domega,tointega2f,integrho) 485 lambda_tr = two * spinfact * integrho(elph_ds%na2f) 486 write (message, '(a,2i3,a,es16.6)' )& 487 & ' mka2f_tr_lova : TRANSPORT lambda for isppol itrtensor', isppol, itrtensor, ' = ', lambda_tr 488 call wrtout(std_out,message,'COLL') 489 if (itrtensor == 1 .or. itrtensor == 5 .or. itrtensor == 9) lambda_tr_trace = lambda_tr_trace + lambda_tr 490 end do !end itrtensor do 491 492 lambda_tr_trace = lambda_tr_trace / three 493 write (message, '(a,i3,a,es16.6)' )& 494 & ' mka2f_tr_lova: 1/3 trace of TRANSPORT lambda for isppol ', isppol, ' = ', lambda_tr_trace 495 call wrtout(std_out,message,'COLL') 496 call wrtout(ab_out,message,'COLL') 497 end do !end isppol do 498 499 !constant to change units of rho from au to SI 500 chgu=2.173969d-7 501 chtu=2.4188843265d-17 502 femto=1.0d-15 503 504 do isppol=1,elph_ds%nsppol 505 do icomp=1, 3 506 do jcomp=1, 3 507 itrtensor=(icomp-1)*3+jcomp 508 509 ! prefactor for resistivity integral 510 ! firh=6.d0*pi*crystal%ucvol*kb_HaK/(elph_ds%n0(isppol)*elph_tr_ds%FSelecveloc_sq(isppol)) 511 ! FIXME: check factor of 2 which is different from Savrasov paper. 6 below for thermal conductivity is correct. 512 firh=2.d0*pi*crystal%ucvol*kb_HaK/elph_ds%n0(isppol)/& 513 & sqrt(elph_tr_ds%FSelecveloc_sq(icomp,isppol)*elph_tr_ds%FSelecveloc_sq(jcomp,isppol)) 514 515 ! Add by BX to get Tau_elph 516 firh_tau = 2.0d0*pi*kb_HaK 517 ! End Adding 518 519 write(unit_rho,*) '# Rho for isppol, itrten = ', isppol, itrtensor 520 write(unit_tau,*) '# Tau for isppol, itrten = ', isppol, itrtensor 521 522 ! jmb 523 tointegtau(:)=0. 524 tointegrho(:)=0. 525 do itemp=1,ntemper ! runs over termperature in K 526 Temp=tempermin+temperinc*dble(itemp) 527 firhT=firh*Temp 528 firhT_tau=firh_tau*Temp 529 omega = omega_min 530 do iomega=1,elph_ds%na2f 531 if(omega<=0) then 532 omega=omega+domega 533 cycle 534 end if 535 xtr=omega/(2*kb_HaK*Temp) 536 if(xtr < log(huge(zero)*tol16)/2)then 537 tointegrho(iomega)=spinfact*firhT*omega*elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1) & 538 & /(((2*Temp*kb_HaK)**2)*((exp(xtr)-exp(-xtr))/2)**2) 539 ! Add by BX to get Tau 540 tointegtau(iomega)=spinfact*firhT_tau*omega*elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1) & 541 & /(((2*Temp*kb_HaK)**2)*((exp(xtr)-exp(-xtr))/2)**2) 542 else 543 tointegrho(iomega)=zero 544 tointegtau(iomega)=zero 545 end if 546 omega=omega+domega 547 end do 548 549 call simpson_int(elph_ds%na2f,domega,tointegrho,integrho) 550 call simpson_int(elph_ds%na2f,domega,tointegtau,integtau) 551 rho=integrho(elph_ds%na2f) 552 tau=1.0d99 553 if(dabs(integtau(elph_ds%na2f)) < tol7) then 554 write(message,'(a)') ' Cannot get a physical relaxation time ' 555 MSG_WARNING(message) 556 else 557 tau=1.0d0/integtau(elph_ds%na2f) 558 end if 559 ! if(elph_ds%na2f < 350.0) then 560 ! tau=1.0d0/integtau(elph_ds%na2f) 561 ! end if 562 write(unit_rho,'(4D20.10)')temp,rho,rho*chgu,rho/temp 563 write(unit_tau,'(3D20.10)')temp,tau,tau*chtu/femto 564 rho_T(itemp)=rho 565 tau_T(itemp)=tau 566 end do ! temperature 567 write(unit_rho,*) 568 write(unit_tau,*) 569 570 end do ! jcomp 571 end do ! icomp 572 end do ! isppol 573 574 !----------------------------- 575 576 577 do isppol=1,elph_ds%nsppol 578 do icomp=1, 3 579 do jcomp=1, 3 580 itrtensor=(icomp-1)*3+jcomp 581 ! prefactor for integral of thermal conductivity 582 ! firh=(18.*crystal%ucvol)/(pi*kb_HaK*elph_ds%n0(isppol)*elph_tr_ds%FSelecveloc_sq(isppol)) 583 firh=(6.d0*crystal%ucvol)/(pi*kb_HaK*elph_ds%n0(isppol))/ & 584 & sqrt(elph_tr_ds%FSelecveloc_sq(icomp,isppol)*elph_tr_ds%FSelecveloc_sq(jcomp,isppol)) 585 586 587 write(unit_therm,*) '# Thermal resistivity for isppol, itrten= ', isppol 588 write(unit_lor,*) '# Lorentz coefficient for isppol, itrten= ', isppol 589 590 tointegrho(:)=0. 591 do itemp=1,ntemper 592 593 Temp=tempermin + temperinc*dble(itemp) 594 omega = omega_min 595 do iomega=1,elph_ds%na2f 596 if(omega<=0) then 597 omega=omega+domega 598 cycle 599 end if 600 xtr=omega/(2*kb_HaK*Temp) 601 if(xtr < log(huge(zero)*tol16)/2)then 602 tointegrho(iomega) = spinfact*xtr**2/omega*& 603 & ( elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1)+& 604 & 4*xtr**2*elph_tr_ds%a2f_1d_trout(iomega,itrtensor,isppol)/pi**2+ & 605 & 2*xtr**2*elph_tr_ds%a2f_1d_trin(iomega,itrtensor,isppol)/pi**2) & 606 & /(((exp(xtr)-exp(-xtr))/2)**2) 607 else 608 tointegrho(iomega) = zero 609 end if 610 omega=omega+domega 611 end do 612 613 call simpson_int(elph_ds%na2f,domega,tointegrho,integrho) 614 wtherm=integrho(elph_ds%na2f)*firh 615 616 if(abs(wtherm) > tol12)then 617 write(unit_therm,'(5D20.10)')temp,wtherm,1./wtherm,wtherm/3.4057d9,1./(wtherm) *3.4057d9 618 619 lorentz=rho_T(itemp)/(wtherm*temp) 620 write(unit_lor,*)temp,lorentz,lor0 621 else 622 write(unit_therm,'(5D20.10)')temp,zero,huge(one),zero,huge(one) 623 write(unit_lor,*)temp,huge(one),lor0 624 end if 625 626 end do 627 write(unit_therm,*) 628 write(unit_lor,*) 629 end do ! jcomp 630 end do ! icomp 631 end do !end isppol do 632 633 634 ABI_DEALLOCATE(phfrq) 635 ABI_DEALLOCATE(displ) 636 ABI_DEALLOCATE(pheigvec) 637 ABI_DEALLOCATE(rho_T) 638 ABI_DEALLOCATE(tau_T) 639 640 close (unit=unit_lor) 641 close (unit=unit_rho) 642 close (unit=unit_tau) 643 close (unit=unit_therm) 644 645 ABI_DEALLOCATE(integrho) 646 ABI_DEALLOCATE(integtau) 647 ABI_DEALLOCATE(tointega2f) 648 ABI_DEALLOCATE(tointegrho) 649 ABI_DEALLOCATE(tointegtau) 650 ABI_DEALLOCATE(elph_tr_ds%a2f_1d_tr) 651 ABI_DEALLOCATE(elph_tr_ds%a2f_1d_trin) 652 ABI_DEALLOCATE(elph_tr_ds%a2f_1d_trout) 653 654 ABI_DEALLOCATE(elph_tr_ds%gamma_qpt_trin) 655 ABI_DEALLOCATE(elph_tr_ds%gamma_qpt_trout) 656 ABI_DEALLOCATE(elph_tr_ds%gamma_rpt_trin) 657 ABI_DEALLOCATE(elph_tr_ds%gamma_rpt_trout) 658 659 !DEBUG 660 write(std_out,*) ' mka2f_tr_lova : end ' 661 !ENDDEBUG 662 663 end subroutine mka2f_tr_lova