TABLE OF CONTENTS
ABINIT/conducti_nc [ Functions ]
NAME
conducti_nc
FUNCTION
This program computes the elements of the optical frequency dependent conductivity tensor and the conductivity along the three principal axes from the Kubo-Greenwood formula.
COPYRIGHT
Copyright (C) 2002-2018 ABINIT group (VRecoules, PGhosh) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
(main routine)
OUTPUT
(main routine)
NOTES
bantot doccde(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy. dom=frequency range eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree). eigen11(2*mband*mband*nkpt_rbz*nsppol)=first-order eigenvalues (hartree) in reciprocal direction 100 eigen12(2*mband*mband*nkpt_rbz*nsppol)=first-order eigenvalues (hartree) in reciprocal direction 010 eigen13(2*mband*mband*nkpt_rbz*nsppol)=first-order eigenvalues (hartree) in reciprocal direction 001 ecut=kinetic energy planewave cutoff (hartree). entropy= entropy associated with the smearing (adimensional) fermie= fermi energy (Hartree) gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{2}$). gmet_inv(3,3)=inverse of reciprocal space metric ($\textrm{bohr}^{2}$). gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1). kin11= Onsager kinetic coeficient=optical conductivity kin12= Onsager kinetic coeficient kin21= Onsager kinetic coeficient kin22= Onsager kinetic coeficient Kth=thermal conductivity mom=number of frequency for conductivity computation mband=maximum number of bands. natom = number of atoms in the unit cell. nband(nkpt*nsppol)=number of bands at each RF k point for each spin. nelect=number of electrons per unit cell nkpt=number of k points in the IBZ for this perturbation ngfft(3)=integer fft box dimensions. nspinor=number of spinorial components of the wavefunctions. nsppol=1 for unpolarized, 2 for spin-polarized. ntypat = number of atom types. occ(mband*nkpt*nsppol)=occupation number for each band and k. occopt==option for occupancies rmet(3,3)=real space metric ($\textrm{bohr}^{2}$). rprimd(3,3)=real space primitive translations. of primitive translations. Sth=thermopower tsmear=smearing width (or temperature) in Hartree ucvol=unit cell volume in ($\textrm{bohr}^{3}$). wind=frequency windows for computations of sigma wtk(nkpt)=weight assigned to each k point. znucl(natom)=atomic number of atoms np_sum=noziere-pines sumrule cond_kg(mom)=kubo-greenwood conductivity
PARENTS
conducti
CHILDREN
getnel,hdr_free,jacobi,matr3inv,metric,msig,nctk_fort_or_ncfile wfk_close,wfk_open_read,wfk_read_eigk
SOURCE
79 #if defined HAVE_CONFIG_H 80 #include "config.h" 81 #endif 82 83 #include "abi_common.h" 84 85 subroutine conducti_nc(filnam,filnam_out,mpi_enreg) 86 87 use defs_basis 88 use defs_abitypes 89 use m_xmpi 90 use m_profiling_abi 91 use m_errors 92 use m_wfk 93 use m_nctk 94 use m_hdr 95 96 use m_io_tools, only : open_file, get_unit 97 use m_geometry, only : metric 98 use m_fstrings, only : sjoin 99 use m_abilasi, only : jacobi 100 use m_occ, only : getnel 101 102 !This section has been created automatically by the script Abilint (TD). 103 !Do not modify the following lines by hand. 104 #undef ABI_FUNC 105 #define ABI_FUNC 'conducti_nc' 106 use interfaces_32_util 107 use interfaces_67_common, except_this_one => conducti_nc 108 !End of the abilint section 109 110 implicit none 111 112 !Arguments ----------------------------------- 113 !scalars 114 character(len=fnlen) :: filnam,filnam_out 115 type(MPI_type),intent(in) :: mpi_enreg 116 117 !Local variables------------------------------- 118 !scalars 119 integer,parameter :: formeig0=0,formeig1=1 120 integer :: bantot,bd2tot_index,bdtot0_index,bdtot_index 121 integer :: headform,iband,ii,jj,ikpt,iunt 122 integer :: index_1,iom,isppol,jband,l1,l2,mband,mom,natom,nband1 123 integer :: nrot,iomode 124 integer :: nband_k,nkpt,nlign,nrest,nspinor,nsppol,ntypat 125 integer :: occopt,comm 126 integer :: tens_unt,lij_unt,sig_unt,kth_unt,ocond_unt 127 real(dp) :: deltae,dosdeltae,diff_occ,dom,ecut,entropy,fermie,maxocc 128 real(dp) :: nelect,np_sum,np_sum_k1,np_sum_k2,omin,oml,socc,socc_k,sig 129 real(dp) :: tphysel,tsmear,ucvol,wind,Tatm 130 character(len=fnlen) :: filnam0,filnam1,filnam2,filnam3 131 character(len=500) :: msg 132 type(hdr_type) :: hdr 133 type(wfk_t) :: gswfk,ddk1,ddk2,ddk3 134 !arrays 135 integer,allocatable :: nband(:) 136 real(dp) :: gmet(3,3),gmet_inv(3,3),gprimd(3,3),gprimd_inv(3,3),rmet(3,3),rprimd(3,3) 137 real(dp),allocatable :: cond_kg(:,:,:),cond_kg_cart(:,:,:),cond_nd(:,:,:),dhdk2_r(:,:,:,:),dhdk2_g(:,:) 138 real(dp),allocatable :: doccde(:),doccde_k(:),cond_kg_xx(:),cond_kg_yy(:),cond_kg_zz(:),trace(:) 139 real(dp),allocatable :: eig0_k(:),eig0tmp(:),eig1_k(:,:),eigen0(:),eigen11(:) 140 real(dp),allocatable :: eigen12(:),eigtmp(:) 141 real(dp),allocatable :: eigen13(:),occ(:),occ_k(:),wtk(:),cond_tot(:),oml1(:) 142 real(dp),allocatable :: kin11(:),kin12(:),kin21(:),kin22(:) 143 real(dp),allocatable :: kin11_k(:),kin12_k(:),kin21_k(:),kin22_k(:),Kth(:),Stp(:) 144 real(dp) :: cond_kg_w(3,3),z(3,3) 145 real(dp) :: eig_cond(3) 146 147 ! ********************************************************************************* 148 149 ABI_UNUSED(mpi_enreg%paral_kgb) 150 151 !Read data file 152 if (open_file(filnam,msg,newunit=iunt,form='formatted', status="old") /= 0) then 153 MSG_ERROR(msg) 154 end if 155 156 rewind(iunt) 157 read(iunt,*) 158 read(iunt,'(a)')filnam1 ! first ddk file 159 read(iunt,'(a)')filnam2 ! second ddk file 160 read(iunt,'(a)')filnam3 ! third ddk file 161 read(iunt,'(a)')filnam0 ! ground-state data 162 163 !Open the GS Wavefunction file and the 3 DDK files. 164 165 ! TODO: one should perform basic consistency tests for the GS WFK and the DDK files, e.g. 166 ! k-points and their order, spins, number of bands could differ in the four files. 167 ! Note indeed that we are assuming the same numer of bands in all the files. 168 comm = xmpi_comm_self 169 call nctk_fort_or_ncfile(filnam0, iomode, msg) 170 if (len_trim(msg) /= 0) MSG_ERROR(msg) 171 call wfk_open_read(gswfk,filnam0, formeig0, iomode, get_unit(), comm) 172 173 call nctk_fort_or_ncfile(filnam1, iomode, msg) 174 if (len_trim(msg) /= 0) MSG_ERROR(msg) 175 call wfk_open_read(ddk1,filnam1, formeig1, iomode, get_unit(), comm, hdr_out=hdr) 176 177 call nctk_fort_or_ncfile(filnam2, iomode, msg) 178 if (len_trim(msg) /= 0) MSG_ERROR(msg) 179 call wfk_open_read(ddk2,filnam2, formeig1, iomode, get_unit(), comm) 180 181 call nctk_fort_or_ncfile(filnam3, iomode, msg) 182 if (len_trim(msg) /= 0) MSG_ERROR(msg) 183 call wfk_open_read(ddk3,filnam3, formeig1, iomode, get_unit(), comm) 184 185 if (wfk_compare(ddk1, ddk2) /= 0) then 186 MSG_ERROR("ddk1 and ddk2 are not consistent. see above messages") 187 end if 188 if (wfk_compare(ddk1, ddk3) /= 0) then 189 MSG_ERROR("ddk1 and ddk3 are not consistent. see above messages") 190 end if 191 192 !Extract params from the header of the first ddk file (might have been the GS file ?) 193 194 !Extract info from the header 195 headform=hdr%headform 196 bantot=hdr%bantot 197 ecut=hdr%ecut_eff 198 natom=hdr%natom 199 nkpt=hdr%nkpt 200 nspinor=hdr%nspinor 201 nsppol=hdr%nsppol 202 ntypat=hdr%ntypat 203 occopt=hdr%occopt 204 rprimd(:,:)=hdr%rprimd(:,:) 205 ABI_ALLOCATE(nband,(nkpt*nsppol)) 206 ABI_ALLOCATE(occ,(bantot)) 207 fermie=hdr%fermie 208 occ(1:bantot)=hdr%occ(1:bantot) 209 nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol) 210 211 !Get mband, as the maximum value of nband(nkpt) 212 mband=maxval(nband(:)) 213 214 write(std_out,*) 215 write(std_out,'(a,3f10.5,a)' )' rprimd(bohr) =',rprimd(1:3,1) 216 write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,2) 217 write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,3) 218 write(std_out,'(a,i8)') ' natom =',natom 219 write(std_out,'(a,2i8)') ' nkpt,mband =',nkpt,mband 220 write(std_out,'(a, f10.5,a)' ) ' ecut =',ecut,' Ha' 221 write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie =',fermie,' Ha',fermie*Ha_eV,' eV' 222 223 !Prepare the reading of ddk Wff files 224 ABI_ALLOCATE(eigtmp,(2*mband*mband)) 225 ABI_ALLOCATE(eig0tmp,(mband)) 226 227 !Read the eigenvalues of ground-state and ddk files 228 ABI_ALLOCATE(eigen0,(mband*nkpt*nsppol)) 229 ABI_ALLOCATE(eigen11,(2*mband*mband*nkpt*nsppol)) 230 ABI_ALLOCATE(eigen12,(2*mband*mband*nkpt*nsppol)) 231 ABI_ALLOCATE(eigen13,(2*mband*mband*nkpt*nsppol)) 232 bdtot0_index=0 ; bdtot_index=0 233 do isppol=1,nsppol 234 do ikpt=1,nkpt 235 nband1=nband(ikpt+(isppol-1)*nkpt) 236 call wfk_read_eigk(gswfk,ikpt,isppol,xmpio_single,eig0tmp) 237 eigen0(1+bdtot0_index:nband1+bdtot0_index)=eig0tmp(1:nband1) 238 239 call wfk_read_eigk(ddk1,ikpt,isppol,xmpio_single,eigtmp) 240 eigen11(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2) 241 242 call wfk_read_eigk(ddk2,ikpt,isppol,xmpio_single,eigtmp) 243 eigen12(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2) 244 245 call wfk_read_eigk(ddk3,ikpt,isppol,xmpio_single,eigtmp) 246 eigen13(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2) 247 248 bdtot0_index=bdtot0_index+nband1 249 bdtot_index=bdtot_index+2*nband1**2 250 end do 251 end do 252 253 ! Close files 254 call wfk_close(gswfk) 255 call wfk_close(ddk1) 256 call wfk_close(ddk2) 257 call wfk_close(ddk3) 258 259 ABI_DEALLOCATE(eigtmp) 260 ABI_DEALLOCATE(eig0tmp) 261 262 !--------------------------------------------------------------------------------- 263 !gmet inversion 264 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 265 call matr3inv(gmet,gmet_inv) 266 call matr3inv(gprimd,gprimd_inv) 267 268 !--------------------------------------------------------------------------------- 269 !derivative of occupation wrt the energy. 270 ABI_ALLOCATE(doccde,(mband*nkpt*nsppol)) 271 ABI_ALLOCATE(wtk,(nkpt)) 272 273 read(iunt,*)tsmear 274 Tatm=tsmear*Ha_K 275 write(std_out,'(a,f12.5,a,f12.5,a)') ' Temp =',tsmear,' Ha ',Tatm,' Kelvin' 276 ! 277 nlign=nkpt/6 278 nrest=nkpt-6*nlign 279 index_1=0 280 do ii=1,nlign 281 read(iunt,*)wtk(1+index_1:6+index_1) 282 index_1=index_1+6 283 end do 284 if (nrest/=0) then 285 read(iunt,*)wtk(6*nlign+1:nkpt) 286 end if 287 ! 288 if (occopt==1) then 289 write(std_out,'(a,i4)') ' occopt =',occopt 290 doccde=0.0d0 291 else 292 tphysel=zero 293 maxocc=two/(nsppol*nspinor) 294 dosdeltae=zero 295 call getnel(doccde,dosdeltae,eigen0,entropy,fermie,maxocc,mband,nband,& 296 & nelect,nkpt,nsppol,occ,occopt,1,tphysel,tsmear,11,wtk) 297 ! DEBUG 298 ! write(std_out,'(a,f10.5)')' getnel : nelect =',nelect 299 ! ENDDEBUG 300 end if 301 !--------------------------------------------------------------------------------- 302 !size of the frequency range 303 read(iunt,*)dom,wind 304 close(iunt) 305 mom=dint(wind/dom) 306 ABI_ALLOCATE(oml1,(mom)) 307 do iom=1,mom 308 oml1(iom)=tol10*1000.0d0+dble(iom)*dom 309 end do 310 311 ABI_ALLOCATE(cond_nd,(mom,3,3)) 312 ABI_ALLOCATE(cond_kg,(mom,3,3)) 313 ABI_ALLOCATE(cond_kg_cart,(mom,3,3)) 314 ABI_ALLOCATE(cond_kg_xx,(mom)) 315 ABI_ALLOCATE(cond_kg_yy,(mom)) 316 ABI_ALLOCATE(trace,(mom)) 317 ABI_ALLOCATE(cond_kg_zz,(mom)) 318 ABI_ALLOCATE(cond_tot,(mom)) 319 ABI_ALLOCATE(kin11,(mom)) 320 ABI_ALLOCATE(kin12,(mom)) 321 ABI_ALLOCATE(kin21,(mom)) 322 ABI_ALLOCATE(kin22,(mom)) 323 ABI_ALLOCATE(kin11_k,(mom)) 324 ABI_ALLOCATE(kin12_k,(mom)) 325 ABI_ALLOCATE(kin21_k,(mom)) 326 ABI_ALLOCATE(kin22_k,(mom)) 327 ABI_ALLOCATE(Kth,(mom)) 328 ABI_ALLOCATE(Stp,(mom)) 329 write(std_out,'(a,i8,2f10.5,a)')' mom,wind,dom =',mom,wind,dom,' Ha' 330 331 !--------------------------------------------------------------------------------- 332 333 kin11 = 0.0d0 334 kin12 = 0.0d0 335 kin21 = 0.0d0 336 kin22 = 0.0d0 337 np_sum = 0.0d0 338 socc = 0.0d0 339 cond_kg = 0.0d0 340 341 342 !LOOP OVER SPINS 343 do isppol=1,nsppol 344 ! 345 bdtot_index = 0 346 bd2tot_index = 0 347 ! 348 deltae = 0.0d0 349 ! 350 ! BIG FAT k POINT LOOP 351 ! 352 do ikpt=1,nkpt 353 ! 354 nband_k=nband(ikpt+(isppol-1)*nkpt) 355 ! 356 ABI_ALLOCATE(eig0_k,(nband_k)) 357 ABI_ALLOCATE(eig1_k,(2*nband_k**2,3)) 358 ABI_ALLOCATE(occ_k,(nband_k)) 359 ABI_ALLOCATE(doccde_k,(nband_k)) 360 ABI_ALLOCATE(dhdk2_r,(nband_k,nband_k,3,3)) 361 ABI_ALLOCATE(dhdk2_g,(nband_k,nband_k)) 362 363 cond_nd = 0.0d0 364 kin11_k = 0.0d0 365 kin12_k = 0.0d0 366 kin21_k = 0.0d0 367 kin22_k = 0.0d0 368 np_sum_k1 = 0.0d0 369 np_sum_k2 = 0.0d0 370 socc_k = 0.0d0 371 dhdk2_r = 0.0d0 372 dhdk2_g = 0.0d0 373 ! 374 ! eigenvalue for k-point 375 eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index) 376 ! first derivative eigenvalues for k-point 377 eig1_k(:,1)=eigen11(1+bd2tot_index:2*nband_k**2+bd2tot_index) 378 eig1_k(:,2)=eigen12(1+bd2tot_index:2*nband_k**2+bd2tot_index) 379 eig1_k(:,3)=eigen13(1+bd2tot_index:2*nband_k**2+bd2tot_index) 380 ! occupation numbers for k-point 381 occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index) 382 ! derivative of occupation number for k-point 383 doccde_k(:)=doccde(1+bdtot_index:nband_k+bdtot_index) 384 ! 385 ! DEBUG 386 ! write(16,*) 387 ! write(16,*)' conducti : ikpt=',ikpt 388 ! do iband=1,nband_k 389 ! write(16, '(i4,4es22.12)' )iband,wtk(ikpt),occ_k(iband),& 390 ! & doccde_k(iband),eig0_k(iband) 391 ! end do 392 ! write(16,*) 393 ! ENDDEBUG 394 ! 395 ! LOOP OVER BAND 396 do iband=1,nband_k 397 do jband=1,nband_k 398 ! 399 do l1=1,3 400 do l2=1,3 401 do ii=1,3 402 do jj=1,3 403 dhdk2_r(iband,jband,l1,l2)=dhdk2_r(iband,jband,l1,l2)+(rprimd(l1,ii)& 404 & *eig1_k(2*iband-1+(jband-1)*2*nband_k,ii)*& 405 & rprimd(l2,jj)*eig1_k(2*iband-1+(jband-1)*2*nband_k,jj)& 406 & +rprimd(l1,ii)*eig1_k(2*iband +(jband-1)*2*nband_k,ii)*& 407 & rprimd(l2,jj)*eig1_k(2*iband+(jband-1)*2*nband_k,jj)) 408 end do 409 end do 410 end do 411 end do 412 413 do l1=1,3 414 do l2=1,3 415 dhdk2_r(iband,jband,l1,l2)=dhdk2_r(iband,jband,l1,l2)/two_pi/two_pi 416 end do 417 end do 418 ! 419 do l1=1,3 420 do l2=1,3 421 dhdk2_g(iband,jband)=dhdk2_g(iband,jband)+gmet_inv(l1,l2)*( & 422 & eig1_k(2*iband-1+(jband-1)*2*nband_k,l1)*& 423 & eig1_k(2*iband-1+(jband-1)*2*nband_k,l2) & 424 & +eig1_k(2*iband +(jband-1)*2*nband_k,l1)*& 425 & eig1_k(2*iband +(jband-1)*2*nband_k,l2)) 426 end do 427 end do 428 dhdk2_g(iband,jband)=dhdk2_g(iband,jband)/two_pi/two_pi 429 ! 430 diff_occ = occ_k(iband)-occ_k(jband) 431 ! if (dabs(diff_occ)>=tol8) then 432 ! 433 ! Conductivity for each omega 434 omin = 0.0d0 435 do iom=1,mom 436 oml=oml1(iom) 437 if (jband>iband) then 438 sig= dhdk2_g(iband,jband)& 439 & *(diff_occ)/oml*(dexp(-((eig0_k(jband)-eig0_k(iband)-oml)/dom)**2)& 440 & -dexp(-((eig0_k(iband)-eig0_k(jband)-oml)/dom)**2)) 441 kin11_k(iom)=kin11_k(iom)+sig 442 kin12_k(iom)=kin12_k(iom)-sig*(eig0_k(jband)-fermie) 443 kin21_k(iom)=kin21_k(iom)-sig*(eig0_k(iband)-fermie) 444 kin22_k(iom)=kin22_k(iom) + & 445 & sig*(eig0_k(iband)-fermie)*(eig0_k(jband)-fermie) 446 end if 447 do l1=1,3 448 do l2=1,3 449 cond_nd(iom,l1,l2)=cond_nd(iom,l1,l2) +dhdk2_r(iband,jband,l1,l2)& 450 & *(diff_occ)/oml*dexp(-((eig0_k(jband)-eig0_k(iband)-oml)/dom)**2) 451 end do 452 end do 453 454 end do 455 ! 456 ! Sumrule start 457 if (dabs(eig0_k(iband)-eig0_k(jband))>=tol10) then 458 np_sum_k1=np_sum_k1 -dhdk2_g(iband,jband)& 459 & *(diff_occ)/(eig0_k(iband)-eig0_k(jband)) 460 else 461 np_sum_k2=np_sum_k2 - doccde_k(iband)*dhdk2_g(iband,jband) 462 end if 463 ! 464 465 ! end loop over band 466 ! end if 467 end do 468 socc_k=socc_k+occ_k(iband) 469 end do 470 ! 471 do iom=1,mom 472 kin11(iom)=kin11(iom)+wtk(ikpt)*kin11_k(iom) 473 kin12(iom)=kin12(iom)+wtk(ikpt)*kin12_k(iom) 474 kin21(iom)=kin21(iom)+wtk(ikpt)*kin21_k(iom) 475 kin22(iom)=kin22(iom)+wtk(ikpt)*kin22_k(iom) 476 do l1=1,3 477 do l2=1,3 478 cond_kg(iom,l1,l2)=cond_kg(iom,l1,l2)+wtk(ikpt)*cond_nd(iom,l1,l2) 479 end do 480 end do 481 end do 482 483 np_sum=np_sum + wtk(ikpt)*(np_sum_k1+np_sum_k2) 484 socc=socc+wtk(ikpt)*socc_k 485 ! 486 ! validity limit 487 deltae=deltae+(eig0_k(nband_k)-fermie) 488 489 bd2tot_index=bd2tot_index+2*nband_k**2 490 bdtot_index=bdtot_index+nband_k 491 ABI_DEALLOCATE(eig0_k) 492 ABI_DEALLOCATE(eig1_k) 493 ABI_DEALLOCATE(occ_k) 494 ABI_DEALLOCATE(doccde_k) 495 ABI_DEALLOCATE(dhdk2_r) 496 ABI_DEALLOCATE(dhdk2_g) 497 ! End loop over k 498 end do 499 500 write(std_out,'(a,3f10.5)')' sumrule =',np_sum/socc/three,socc 501 write(std_out,'(a,f10.5,a,f10.5,a)')& 502 & ' Emax-Efermi =',deltae/dble(nkpt),' Ha',deltae/dble(nkpt)*Ha_eV,' eV' 503 504 ! End loop over spins 505 end do 506 507 cond_kg=cond_kg*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 508 509 510 !Check that new output file does NOT exist 511 !Keep this line : prevent silly (compiler ?) bug on HP 8000 512 write(std_out,*)' conducti : call isfile ' 513 ! 514 if (open_file(trim(filnam_out)//'_tens',msg,newunit=tens_unt,form='formatted',action="write") /= 0) then 515 MSG_ERROR(msg) 516 end if 517 if (open_file(trim(filnam_out)//'_Lij',msg,newunit=lij_unt,form='formatted',action="write") /= 0) then 518 MSG_ERROR(msg) 519 end if 520 write(lij_unt,'(a)')' # omega(ua) L12 L21 L22 L22' 521 522 if (open_file(trim(filnam_out)//'_sig',msg,newunit=sig_unt,form='formatted',action="write") /= 0) then 523 MSG_ERROR(msg) 524 end if 525 write(sig_unt,'(a)')' # omega(ua) hbar*omega(eV) cond(ua) cond(ohm.cm)-1' 526 527 if (open_file(trim(filnam_out)//'_Kth',msg,newunit=kth_unt,form='formatted',action="write") /= 0) then 528 MSG_ERROR(msg) 529 end if 530 write(kth_unt,'(a)')& 531 & ' #omega(ua) hbar*omega(eV) thermal cond(ua) Kth(W/m/K) thermopower(ua) Stp(microohm/K)' 532 533 if (open_file(trim(filnam_out)//'.out',msg,newunit=ocond_unt,form='formatted',action="write") /= 0) then 534 MSG_ERROR(msg) 535 end if 536 write(ocond_unt,'(a)' )' Conducti output file:' 537 write(ocond_unt,'(a)' )' Contains all results produced by conducti utility' 538 write(ocond_unt,'(a)' )' ' 539 write(ocond_unt,'(a)')' # omega(ua) cond(ua) thermal cond(ua) thermopower(ua)' 540 ! 541 !call isfile(filnam_out,'new') 542 543 !Keep this line : prevent silly (compiler ?) bug on HP 8000 544 write(std_out,*)' conducti : after call isfile ' 545 ! 546 !Compute thermal conductivity and thermopower 547 do iom=1,mom 548 oml=oml1(iom) 549 kin11(iom)=kin11(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 550 kin21(iom)=kin21(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 551 kin12(iom)=kin12(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 552 kin22(iom)=kin22(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 553 if (dabs(kin11(iom))<10.0d-20) kin11(iom)=zero 554 Kth(iom)=kin22(iom) 555 Stp(iom)=zero 556 if(kin11(iom)/=zero) then 557 Kth(iom)=Kth(iom)-(kin12(iom)*kin21(iom)/kin11(iom)) 558 Stp(iom)=kin12(iom)/(kin11(iom)*Tatm) 559 end if 560 if (dabs(Kth(iom))<10.0d-20) Kth(iom)=0.0d0 561 if (dabs(Stp(iom))<10.0d-20) Stp(iom)=0.0d0 562 if (dabs(kin12(iom))<10.0d-20) kin12(iom)=zero 563 if (dabs(kin21(iom))<10.0d-20) kin21(iom)=zero 564 if (dabs(kin22(iom))<10.0d-20) kin22(iom)=zero 565 566 write(lij_unt,'(f12.5,4es22.12)')oml,kin12(iom),kin21(iom),kin22(iom),kin22(iom)/Tatm*3.4057d9 567 write(sig_unt,'(2f12.5,2es22.12)') oml,oml*Ha_eV,kin11(iom),kin11(iom)*Ohmcm 568 write(kth_unt,'(2f12.5,4es22.12)') oml,oml*Ha_eV,Kth(iom),Kth(iom)*3.4057d9/Tatm,Stp(iom),Stp(iom)*3.6753d-2 569 write(ocond_unt,'(1f12.5,3es22.12)') oml,kin11(iom),Kth(iom),Stp(iom) 570 end do 571 ! 572 573 write(tens_unt,'(a)' )' Conductivity file ' 574 write(tens_unt,'(a)' )' ----------------- ' 575 write(tens_unt,'(a)' )' Contain first the full conductivity tensor, for the desired set of energies,' 576 write(tens_unt,'(a)' )' then, the three principal values, for the desired set of energies' 577 write(tens_unt,'(a)' )' (note that eigenvalues are not directly associated with xx,yy,zz)' 578 write(tens_unt,'(a)' )' ' 579 580 write(ocond_unt,'(a)' )' ' 581 write(ocond_unt,'(a)' )' full conductivity tensor, for the desired set of energies' 582 write(ocond_unt,'(a)' )' then, the three principal values, for the desired set of energies:' 583 584 do iom=1,mom 585 oml=oml1(iom)*Ha_eV 586 write(tens_unt, '(a,es16.6,a)' ) ' energy (in eV) =',oml,', conductivity tensor (in Ohm.cm-1) follows :' 587 write(ocond_unt, '(a,es16.6,a)' ) ' energy (in eV) =',oml,', conductivity tensor (in Ohm.cm-1) follows :' 588 do l1=1,3 589 write(tens_unt,"(3f25.15)") (cond_kg(iom,l1,l2)*Ohmcm,l2=1,3) 590 write(ocond_unt,"(3f25.15)") (cond_kg(iom,l1,l2)*Ohmcm,l2=1,3) 591 end do 592 end do 593 594 !Diagonalizing the conductivity matrix for sigma_xx,sigma_yy,sigma_zz 595 cond_kg_xx=0d0 596 cond_kg_yy=0d0 597 cond_kg_zz=0d0 598 !trace=0d0 used for checking with the original version of the code 599 do iom=1,mom 600 oml=oml1(iom)*Ha_eV 601 cond_kg_w=0d0 602 do l1=1,3 603 do l2=1,3 604 cond_kg_w(l1,l2)=cond_kg(iom,l1,l2) 605 end do 606 end do 607 call jacobi(cond_kg_w,3,3,eig_cond,z,nrot) 608 609 ! When the value is too small, set it to zero before printing 610 if(abs(eig_cond(1))<tol10)eig_cond(1)=zero 611 if(abs(eig_cond(2))<tol10)eig_cond(2)=zero 612 if(abs(eig_cond(3))<tol10)eig_cond(3)=zero 613 614 cond_kg_xx(iom)=eig_cond(1) 615 cond_kg_yy(iom)=eig_cond(2) 616 cond_kg_zz(iom)=eig_cond(3) 617 ! trace(iom)=cond_kg_xx(iom)+cond_kg_yy(iom)+cond_kg_zz(iom) 618 end do 619 620 !DEBUG Keep this line : prevent silly (compiler ?) bug on HP 8000 621 !write(std_out,*)' conducti : after open ' 622 !ENDDEBUG 623 624 write(tens_unt,'(a,a)')ch10,' Now, print principal values of the conductivity tensor.' 625 write(tens_unt,'(a)')' ' 626 write(tens_unt,'(a)')' #omega(ua) cond_1(ua) cond_2(ua) cond_3(ua) cond_tot(ua)' 627 628 write(ocond_unt,'(a)')' ' 629 write(ocond_unt,'(a,a)')ch10,' Now, print principal values of the conductivity tensor.' 630 write(ocond_unt,'(a)')' ' 631 write(ocond_unt,'(a)')' #omega(ua) cond_1(ua) cond_2(ua) cond_3(ua) cond_tot(ua)' 632 633 634 do iom=1,mom 635 cond_tot(iom)=cond_kg_xx(iom)+cond_kg_yy(iom)+cond_kg_zz(iom) 636 write(tens_unt,'(f12.5,4es22.12)')oml1(iom),cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom) 637 write(ocond_unt,'(f12.5,4es22.12)')oml1(iom),cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom) 638 end do 639 640 write(tens_unt,*) 641 write(tens_unt,'(a)')' #hbar*omega(eV) cond_1(ohm.cm)-1 cond_2(ohm.cm)-1 cond_3(ohm.cm)-1 cond_t(ohm.cm)-1' 642 write(ocond_unt,*) 643 write(ocond_unt,'(a)')' #hbar*omega(eV) cond_1(ohm.cm)-1 cond_2(ohm.cm)-1 cond_3(ohm.cm)-1 cond_t(ohm.cm)-1' 644 645 do iom=1,mom 646 oml=oml1(iom)*Ha_eV 647 cond_tot(iom)=cond_tot(iom)*Ohmcm 648 cond_kg_xx(iom)=cond_kg_xx(iom)*Ohmcm 649 cond_kg_yy(iom)=cond_kg_yy(iom)*Ohmcm 650 cond_kg_zz(iom)=cond_kg_zz(iom)*Ohmcm 651 write(tens_unt,'(f12.5,4es22.12)')oml,cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom) 652 write(ocond_unt,'(f12.5,4es22.12)')oml,cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom) 653 end do 654 !Calculate the imaginary part of the conductivity (principal value) 655 !+derived optical properties. 656 657 call msig (kin11,mom,oml1,filnam_out) 658 659 close(tens_unt) 660 close(lij_unt) 661 close(sig_unt) 662 close(kth_unt) 663 close(ocond_unt) 664 665 ABI_DEALLOCATE(nband) 666 ABI_DEALLOCATE(oml1) 667 ABI_DEALLOCATE(occ) 668 ABI_DEALLOCATE(eigen11) 669 ABI_DEALLOCATE(eigen12) 670 ABI_DEALLOCATE(eigen13) 671 ABI_DEALLOCATE(eigen0) 672 ABI_DEALLOCATE(doccde) 673 ABI_DEALLOCATE(wtk) 674 ABI_DEALLOCATE(cond_nd) 675 ABI_DEALLOCATE(cond_kg) 676 ABI_DEALLOCATE(cond_kg_cart) 677 ABI_DEALLOCATE(cond_kg_xx) 678 ABI_DEALLOCATE(cond_kg_yy) 679 ABI_DEALLOCATE(trace) 680 ABI_DEALLOCATE(cond_kg_zz) 681 ABI_DEALLOCATE(cond_tot) 682 ABI_DEALLOCATE(kin11) 683 ABI_DEALLOCATE(kin22) 684 ABI_DEALLOCATE(kin12) 685 ABI_DEALLOCATE(kin21) 686 ABI_DEALLOCATE(kin11_k) 687 ABI_DEALLOCATE(kin22_k) 688 ABI_DEALLOCATE(kin12_k) 689 ABI_DEALLOCATE(kin21_k) 690 ABI_DEALLOCATE(Stp) 691 ABI_DEALLOCATE(Kth) 692 693 call hdr_free(hdr) 694 695 end subroutine conducti_nc