TABLE OF CONTENTS
ABINIT/conducti_paw [ Functions ]
NAME
conducti_paw
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 for PAW formalism
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/Infos/contributors .
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,nkpt,mband,mband,nsppol)=first-order eigenvalues (hartree) in direction x eigen12(2,nkpt,mband,mband,nsppol)=first-order eigenvalues (hartree) in direction y eigen13(2,nkpt,mband,mband,nsppol)=first-order eigenvalues (hartree) in direction z ecut=kinetic energy planewave cutoff (hartree). fermie= fermi energy (Hartree) gmet(3,3)=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. 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}$).sigx(mom,nphicor)) 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
PARENTS
conducti
CHILDREN
hdr_free,hdr_io,hdr_read_from_fname,metric,msig,wffclose,wffopen
SOURCE
74 #if defined HAVE_CONFIG_H 75 #include "config.h" 76 #endif 77 78 #include "abi_common.h" 79 80 subroutine conducti_paw(filnam,filnam_out,mpi_enreg) 81 82 use defs_basis 83 use defs_abitypes 84 use m_errors 85 use m_xmpi 86 use m_wffile 87 use m_profiling_abi 88 use m_hdr 89 90 use m_io_tools, only : open_file, get_unit 91 use m_geometry, only : metric 92 93 !This section has been created automatically by the script Abilint (TD). 94 !Do not modify the following lines by hand. 95 #undef ABI_FUNC 96 #define ABI_FUNC 'conducti_paw' 97 use interfaces_67_common, except_this_one => conducti_paw 98 !End of the abilint section 99 100 implicit none 101 102 !Arguments ----------------------------------- 103 !scalars 104 character(len=fnlen) :: filnam,filnam_out 105 type(MPI_type),intent(in) :: mpi_enreg 106 107 !Local variables------------------------------- 108 !scalars 109 integer,parameter :: master=0 110 integer :: iomode,bantot,bdtot_index 111 integer :: fform1,headform,iband,ierr,ikpt 112 integer :: iom,isppol,jband,l1,l2,mband,me,mom 113 integer :: natom,nband_k,nkpt,nspinor,nsppol,ntypat 114 integer :: occopt,rdwr,spaceComm,iunt,opt_unt 115 integer :: lij_unt,sig_unt,kth_unt,ocond_unt 116 real(dp) :: del,deltae,diff_occ,ecut,fermie,maxocc 117 real(dp) :: np_sum,np_sum_k1,np_sum_k2,omin,omax,dom,oml,sig,socc,socc_k 118 real(dp) :: Tatm,tphysel,tsmear,ucvol 119 character(len=fnlen) :: filnam1,filnam_gen 120 character(len=500) :: msg 121 type(hdr_type) :: hdr 122 type(wffile_type) :: wff1 123 !arrays 124 integer,allocatable :: nband(:) 125 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3) 126 real(dp),allocatable :: cond_nd(:,:,:),dhdk2_r(:,:,:,:),dhdk2_g(:,:,:) 127 real(dp),allocatable :: doccde(:),doccde_k(:),eig0_k(:),eigen0(:) 128 real(dp),allocatable :: occ(:),occ_k(:),wtk(:),oml1(:) 129 real(dp),allocatable :: kin11(:,:),kin12(:),kin21(:),kin22(:) 130 real(dp),allocatable :: kin11_k(:),kin12_k(:),kin21_k(:),kin22_k(:),Kth(:),Stp(:) 131 real(dp),allocatable :: psinablapsi(:,:,:,:),sig_abs(:) 132 133 ! ********************************************************************************* 134 ABI_UNUSED(mpi_enreg%paral_kgb) 135 136 !write(std_out,'(a)')' The name of the output file is :',trim(filnam_out) 137 !Read data file 138 if (open_file(filnam, msg, newunit=iunt, form='formatted', action="read", status="old") /=0 ) then 139 MSG_ERROR(msg) 140 end if 141 rewind(iunt) 142 read(iunt,*) 143 read(iunt,'(a)')filnam_gen ! generic name for the files 144 filnam1=trim(filnam_gen)//'_OPT' 145 !Read size of the frequency range 146 read(iunt,*) dom,omin,omax,mom 147 close(iunt) 148 write(std_out,'(a,i8,3f10.5,a)')' npts,omin,omax,width =',mom,omin,omax,dom,' Ha' 149 150 !These default values are typical of sequential use 151 iomode=IO_MODE_FORTRAN ; spaceComm=xmpi_comm_self; me=0 152 153 ! Read the header of the optic files 154 call hdr_read_from_fname(hdr, filnam1, fform1, spaceComm) 155 call hdr_free(hdr) 156 if (fform1 /= 610) then 157 MSG_ERROR("Abinit8 requires an OPT file with fform = 610") 158 end if 159 160 !Open the conducti and/or optic files 161 opt_unt = get_unit() 162 call WffOpen(iomode,spaceComm,filnam1,ierr,wff1,master,me,opt_unt) 163 164 !Read the header from Ground state file 165 rdwr=1 166 call hdr_io(fform1,hdr,rdwr,wff1) 167 ABI_CHECK(fform1/=0,"Error opening wff1") 168 169 !Extract info from the header 170 headform=hdr%headform 171 bantot=hdr%bantot 172 ecut=hdr%ecut_eff 173 natom=hdr%natom 174 nkpt=hdr%nkpt 175 nspinor=hdr%nspinor 176 nsppol=hdr%nsppol 177 ntypat=hdr%ntypat 178 occopt=hdr%occopt 179 rprimd(:,:)=hdr%rprimd(:,:) 180 ABI_ALLOCATE(nband,(nkpt*nsppol)) 181 ABI_ALLOCATE(occ,(bantot)) 182 ABI_ALLOCATE(wtk,(nkpt)) 183 fermie=hdr%fermie 184 tsmear=hdr%tsmear 185 occ(1:bantot)=hdr%occ(1:bantot) 186 wtk(1:nkpt)=hdr%wtk(1:nkpt) 187 nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol) 188 189 !Get mband, as the maximum value of nband(nkpt) 190 mband=maxval(nband(:)) 191 192 write(std_out,*) 193 write(std_out,'(a,3f10.5,a)' )' rprimd(bohr) =',rprimd(1:3,1) 194 write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,2) 195 write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,3) 196 write(std_out,'(a,i8)') ' natom =',natom 197 write(std_out,'(a,3i8)') ' nkpt,mband,nsppol =',nkpt,mband,nsppol 198 write(std_out, '(a, f10.5,a)' ) ' ecut =',ecut,' Ha' 199 write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie =',fermie,' Ha',fermie*Ha_eV,' eV' 200 Tatm=tsmear*Ha_K 201 write(std_out,'(a,f12.5,a,f12.5,a)') ' Temp =',tsmear,' Ha ',Tatm,' Kelvin' 202 203 ABI_ALLOCATE(eigen0,(mband*nkpt*nsppol)) 204 read(opt_unt)(eigen0(iband),iband=1,mband*nkpt*nsppol) 205 ! 206 ! 207 !--------------------------------------------------------------------------------- 208 !gmet inversion to get ucvol 209 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 210 211 !--------------------------------------------------------------------------------- 212 !derivative of occupation wrt the energy. 213 ABI_ALLOCATE(doccde,(mband*nkpt*nsppol)) 214 if (occopt==1) then 215 write(std_out,'(a,i4)') ' occopt =',occopt 216 doccde=zero 217 else 218 tphysel=zero 219 maxocc=two/(nsppol*nspinor) 220 end if 221 !--------------------------------------------------------------------------------- 222 !size of the frequency range 223 del=(omax-omin)/(mom-1) 224 ABI_ALLOCATE(oml1,(mom)) 225 do iom=1,mom 226 oml1(iom)=omin+dble(iom-1)*del 227 end do 228 ABI_ALLOCATE(cond_nd,(mom,3,3)) 229 ABI_ALLOCATE(kin11,(mom,nsppol)) 230 ABI_ALLOCATE(kin12,(mom)) 231 ABI_ALLOCATE(kin21,(mom)) 232 ABI_ALLOCATE(kin22,(mom)) 233 ABI_ALLOCATE(sig_abs,(mom)) 234 ABI_ALLOCATE(kin11_k,(mom)) 235 ABI_ALLOCATE(kin12_k,(mom)) 236 ABI_ALLOCATE(kin21_k,(mom)) 237 ABI_ALLOCATE(kin22_k,(mom)) 238 ABI_ALLOCATE(Kth,(mom)) 239 ABI_ALLOCATE(Stp,(mom)) 240 241 !--------------------------------------------------------------------------------- 242 !Conductivity ------- 243 ! 244 ABI_ALLOCATE(psinablapsi,(2,3,mband,mband)) 245 kin11 = zero 246 kin12 = zero 247 kin21 = zero 248 kin22 = zero 249 np_sum = zero 250 socc = zero 251 sig_abs = zero 252 253 bdtot_index = 0 254 255 !LOOP OVER SPINS/K 256 deltae = zero 257 do isppol=1,nsppol 258 do ikpt=1,nkpt 259 nband_k=nband(ikpt+(isppol-1)*nkpt) 260 ABI_ALLOCATE(eig0_k,(nband_k)) 261 ABI_ALLOCATE(occ_k,(nband_k)) 262 ABI_ALLOCATE(doccde_k,(nband_k)) 263 ABI_ALLOCATE(dhdk2_r,(nband_k,nband_k,3,3)) 264 ABI_ALLOCATE(dhdk2_g,(natom,nband_k,nband_k)) 265 266 cond_nd = zero 267 kin11_k = zero 268 kin12_k = zero 269 kin21_k = zero 270 kin22_k = zero 271 np_sum_k1 = zero 272 np_sum_k2 = zero 273 socc_k = zero 274 dhdk2_r = zero 275 dhdk2_g = zero 276 277 ! eigenvalue for k-point 278 eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index) 279 ! first derivative eigenvalues for k-point 280 psinablapsi=zero 281 read(opt_unt)((psinablapsi(1:2,1,iband,jband),iband=1,nband_k),jband=1,nband_k) 282 read(opt_unt)((psinablapsi(1:2,2,iband,jband),iband=1,nband_k),jband=1,nband_k) 283 read(opt_unt)((psinablapsi(1:2,3,iband,jband),iband=1,nband_k),jband=1,nband_k) 284 ! DEBUG 285 ! write(963,*)isppol,ikpt,((psinablapsi(1:2,1,iband,jband),iband=1,nband_k),jband=1,nband_k) 286 ! write(963,*)isppol,ikpt,((psinablapsi(1:2,2,iband,jband),iband=1,nband_k),jband=1,nband_k) 287 ! write(963,*)isppol,ikpt,((psinablapsi(1:2,3,iband,jband),iband=1,nband_k),jband=1,nband_k) 288 ! ENDDEBUG 289 290 ! occupation numbers for k-point 291 occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index) 292 ! derivative of occupation number for k-point 293 doccde_k(:)=doccde(1+bdtot_index:nband_k+bdtot_index) 294 295 ! LOOP OVER BANDS 296 do iband=1,nband_k 297 do jband=1,nband_k 298 do l1=1,3 299 do l2=1,3 300 dhdk2_r(iband,jband,l1,l2)=dhdk2_r(iband,jband,l1,l2)+(& 301 & psinablapsi(1,l1,iband,jband)*psinablapsi(1,l2,iband,jband)& 302 & +psinablapsi(2,l1,iband,jband)*psinablapsi(2,l2,iband,jband)) 303 end do 304 end do 305 306 do l1=1,3 307 dhdk2_g(1,iband,jband)=dhdk2_g(1,iband,jband)+( & 308 & psinablapsi(1,l1,iband,jband)*psinablapsi(1,l1,iband,jband) & 309 & +psinablapsi(2,l1,iband,jband)*psinablapsi(2,l1,iband,jband)) 310 end do 311 312 diff_occ = occ_k(iband)-occ_k(jband) 313 if (dabs(diff_occ)>=tol8) then 314 315 ! Conductivity for each omega 316 ! omin = zero 317 do iom=1,mom 318 oml=oml1(iom) 319 if (jband>iband) then 320 sig= dhdk2_g(1,iband,jband)& 321 & *(diff_occ)/oml*(dexp(-((eig0_k(jband)-eig0_k(iband)-oml)/dom)**2)& 322 & -dexp(-((eig0_k(iband)-eig0_k(jband)-oml)/dom)**2)) 323 kin11_k(iom)=kin11_k(iom)+sig 324 kin12_k(iom)=kin12_k(iom)-sig*(eig0_k(jband)-fermie) 325 kin21_k(iom)=kin21_k(iom)-sig*(eig0_k(iband)-fermie) 326 kin22_k(iom)=kin22_k(iom) + & 327 & sig*(eig0_k(iband)-fermie)*(eig0_k(jband)-fermie) 328 end if 329 do l1=1,3 330 do l2=1,3 331 cond_nd(iom,l1,l2)=cond_nd(iom,l1,l2) +dhdk2_r(iband,jband,l1,l2)& 332 & *(diff_occ)/oml*dexp(-((eig0_k(jband)-eig0_k(iband)-oml)/dom)**2) 333 end do 334 end do 335 end do 336 337 ! Sumrule start 338 if (dabs(eig0_k(iband)-eig0_k(jband))>=tol10) then 339 np_sum_k1=np_sum_k1 -dhdk2_g(1,iband,jband)& 340 & *(diff_occ)/(eig0_k(iband)-eig0_k(jband)) 341 else 342 np_sum_k2=np_sum_k2 - doccde_k(iband)*dhdk2_g(1,iband,jband) 343 end if 344 345 ! end loop over band 346 end if 347 end do 348 socc_k=socc_k+occ_k(iband) 349 end do 350 do iom=1,mom 351 kin11(iom,isppol)=kin11(iom,isppol)+wtk(ikpt)*kin11_k(iom) 352 kin12(iom)=kin12(iom)+wtk(ikpt)*kin12_k(iom) 353 kin21(iom)=kin21(iom)+wtk(ikpt)*kin21_k(iom) 354 kin22(iom)=kin22(iom)+wtk(ikpt)*kin22_k(iom) 355 end do 356 np_sum=np_sum + wtk(ikpt)*(np_sum_k1+np_sum_k2) 357 socc=socc+wtk(ikpt)*socc_k 358 359 ! Validity limit 360 deltae=deltae+(eig0_k(nband_k)-fermie) 361 362 bdtot_index=bdtot_index+nband_k 363 ABI_DEALLOCATE(eig0_k) 364 ABI_DEALLOCATE(occ_k) 365 ABI_DEALLOCATE(doccde_k) 366 ABI_DEALLOCATE(dhdk2_r) 367 ABI_DEALLOCATE(dhdk2_g) 368 ! End loop over k 369 end do 370 ! End loop over Spin 371 end do 372 373 write(std_out,'(a,3f10.5)')' sumrule =',np_sum/socc/three/dble(nsppol),socc 374 write(std_out,'(a,f10.5,a,f10.5,a)')& 375 & ' Emax-Efermi =',deltae/dble(nkpt*nsppol),' Ha',deltae/dble(nkpt*nsppol)*Ha_eV,' eV' 376 377 if (open_file(trim(filnam_out)//'_Lij',msg, newunit=lij_unt, form='formatted', action="write") /= 0) then 378 MSG_ERROR(msg) 379 end if 380 write(lij_unt,'(a)')' # omega(ua) L12 L21 L22 L22' 381 382 if (open_file(trim(filnam_out)//'_sig', msg, newunit=sig_unt, form='formatted', action="write") /= 0) then 383 MSG_ERROR(msg) 384 end if 385 if (nsppol==1) then 386 write(sig_unt,'(a)')' # omega(ua) hbar*omega(eV) cond(ua) cond(ohm.cm)-1' 387 else 388 write(sig_unt,'(2a)')' # omega(ua) hbar*omega(eV) cond(ua) cond(ohm.cm)-1',& 389 & ' cond(ohm.cm)-1 UP cond(ohm.cm)-1 DN' 390 end if 391 392 if (open_file(trim(filnam_out)//'_Kth', msg, newunit=kth_unt, form='formatted', action="write") /=0) then 393 MSG_ERROR(msg) 394 end if 395 write(kth_unt,'(a)')& 396 & ' #omega(ua) hbar*omega(eV) thermal cond(ua) Kth(W/m/K) thermopower(ua) Stp(microohm/K)' 397 398 if (open_file(trim(filnam_out)//'.out', msg, newunit=ocond_unt, form='formatted', action="write") /= 0) then 399 MSG_ERROR(msg) 400 end if 401 write(ocond_unt,'(a)' )' #Conducti output file:' 402 write(ocond_unt,'(a)' )' #Contains all results produced by conducti utility' 403 write(ocond_unt,'(a)' )' ' 404 write(ocond_unt,'(a)')' # omega(ua) cond(ua) thermal cond(ua) thermopower(ua)' 405 406 !call isfile(filnam_out,'new') 407 408 !Compute thermal conductivity and thermopower 409 do iom=1,mom 410 oml=oml1(iom) 411 do isppol=1,nsppol 412 kin11(iom,isppol)=kin11(iom,isppol)*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 413 if (dabs(kin11(iom,isppol))<10.0d-20) kin11(iom,isppol)=zero 414 sig_abs(iom)=sig_abs(iom)+kin11(iom,isppol) 415 end do 416 kin21(iom)=kin21(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 417 kin12(iom)=kin12(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 418 kin22(iom)=kin22(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 419 Kth(iom)=kin22(iom) 420 Stp(iom)=zero 421 if(sig_abs(iom)/=zero) then 422 Kth(iom)=Kth(iom)-(kin12(iom)*kin21(iom)/sig_abs(iom)) 423 Stp(iom)=kin12(iom)/(sig_abs(iom)*Tatm) 424 end if 425 if (dabs(Kth(iom))<10.0d-20) Kth(iom)=zero 426 if (dabs(Stp(iom))<10.0d-20) Stp(iom)=zero 427 if (abs(kin12(iom))<10.0d-80) kin12=zero 428 if (abs(kin21(iom))<10.0d-80) kin21=zero 429 if (abs(kin22(iom))<10.0d-80) kin22=zero 430 431 write(lij_unt,'(f12.5,4es22.12)')oml,kin12(iom),kin21(iom),kin22(iom),kin22(iom)/Tatm*3.4057d9 432 433 if (nsppol==1) then 434 write(sig_unt,'(2f12.5,2es22.12)') oml,oml*Ha_eV,sig_abs(iom),sig_abs(iom)*Ohmcm 435 else 436 write(sig_unt,'(2f12.5,4es22.12)') oml,oml*Ha_eV,sig_abs(iom),sig_abs(iom)*Ohmcm,& 437 & kin11(iom,1)*Ohmcm,kin11(iom,2)*Ohmcm 438 end if 439 write(kth_unt,'(2f12.5,4es22.12)') oml,oml*Ha_eV,Kth(iom),Kth(iom)*3.4057d9/Tatm,& 440 & Stp(iom),Stp(iom)*3.6753d-2 441 write(ocond_unt,'(1f12.5,3es22.12)') oml,sig_abs(iom),Kth(iom),Stp(iom) 442 end do 443 444 !Calculate the imaginary part of the conductivity (principal value) 445 !+derived optical properties. 446 call msig (sig_abs,mom,oml1,filnam_out) 447 448 close(lij_unt) 449 close(sig_unt) 450 close(kth_unt) 451 close(ocond_unt) 452 453 write(std_out,'(2a)')ch10,'OUTPUT' 454 write(std_out,'(a)')trim(filnam_out)//'_Lij : Onsager kinetic coefficients' 455 write(std_out,'(a)')trim(filnam_out)//'_sig : Optical conductivity' 456 write(std_out,'(a)')trim(filnam_out)//'_Kth : Thermal conductivity and thermopower' 457 write(std_out,'(a)')trim(filnam_out)//'_eps : Dielectric fonction' 458 write(std_out,'(a)')trim(filnam_out)//'_abs : n, k, reflectivity, absorption' 459 460 call WffClose(wff1,ierr) 461 462 ABI_DEALLOCATE(psinablapsi) 463 ABI_DEALLOCATE(kin11) 464 ABI_DEALLOCATE(kin22) 465 ABI_DEALLOCATE(kin12) 466 ABI_DEALLOCATE(kin21) 467 ABI_DEALLOCATE(kin11_k) 468 ABI_DEALLOCATE(kin22_k) 469 ABI_DEALLOCATE(kin12_k) 470 ABI_DEALLOCATE(kin21_k) 471 ABI_DEALLOCATE(Stp) 472 ABI_DEALLOCATE(Kth) 473 ABI_DEALLOCATE(cond_nd) 474 ABI_DEALLOCATE(sig_abs) 475 ABI_DEALLOCATE(eigen0) 476 ABI_DEALLOCATE(nband) 477 ABI_DEALLOCATE(oml1) 478 ABI_DEALLOCATE(occ) 479 ABI_DEALLOCATE(doccde) 480 ABI_DEALLOCATE(wtk) 481 482 call hdr_free(hdr) 483 484 end subroutine conducti_paw