1 program fermi_int_1 2 ! 3 ! Written by Burak Himmetoglu (burakhmmtgl@gmail.com) 4 ! Uses some parts of the PW distribution 5 ! 6 ! This program computes the transport integrals at a given Fermi level using 7 ! constant scattering rate, and results in conductivity and Seebeck coefficients. 8 ! 9 ! Integrals that are computed: 10 ! 11 ! I0 = \sum_{n,k} (-df_nk/de) 12 ! I1(i,j) = \sum_{n,k} tau (-df_nk/de) v_nk,i v_nk,j 13 ! I2(i,j) = \sum_{n,k} tau (-df_nk/de) (e_nk - e_F) v_nk,i v_nk,j 14 ! 15 ! The difference between fermi_int_0 is that the calculation proceeds in a reduced 16 ! grid such that |e_nk-Ef| <= cut*kT. This allows faster integration times for a 17 ! given Fermi level (Ef). 18 ! 19 ! Description of the input card: 20 ! 21 ! fil_a2F : prefix.a2Fsave file 22 ! fil_info : file containing direct and reciprocal lattice vectors 23 ! T : Temperature in K 24 ! alat : lattice parameter in a.u. (Bohr) 25 ! vol : volume of lattice in a.u.^3 26 ! cut : Reduction parameter for integrals (i.e. Integration reduced to the region satisfying |e_nk-Ef| <= cut*kT) 27 ! invtau : inverse tau in eV 28 ! phband_i, phband_f : starting and ending indices of bands of interest. The bands must lie between efmin & efmax 29 ! lsoc : if .true. then the band structure is non-collinear 30 ! nthreads : Number of threads for OpenMP parallelization 31 ! lscissors : if .true. conduction bands are shifted by a constant energy up 32 ! shift : value of the shift applied to conduction bands 33 ! cbm_i : the initial conduction band for which the shift is applied ( all bands with index >= cbm_i are shifted up) 34 35!$ USE omp_lib 36 37 implicit none 38 39 integer :: i,j,k,nu,ik,ikk,nk1fit,nk2fit,nk3fit,nkfit, & 40 & nbnd, nksfit, npk, nsym, Nmu, imu, & 41 & s(3,3,48),ns,nrot,ibnd,io,phband_i, phband_f, & 42 & nphband, n, nn, jbnd, ibnd_ph, ind_k, cbm_i 43 ! 44 double precision :: wk, at(3,3), bg(3,3), efermi, alat, & 45 & T, wo(3), al(3), invtau,aa,cut,deg, & 46 & invT, fd, dfd, fac, vol, shift 47 ! 48 logical :: lsoc, lscissors 49 ! 50 double precision, allocatable :: xkfit(:,:),etfit(:,:),wkfit(:), & 51 & dfk(:,:,:),vk(:,:,:) 52 ! 53 double precision :: I0,I1(3,3),I2(3,3),sig(3,3),Se(3,3),inv_I1(3,3) 54 ! 55 integer, allocatable :: eqkfit(:), sfit(:), iflag(:,:),nkeff(:) 56 ! 57 ! OMP variables 58 double precision :: t0 59 ! 60 integer :: nthreads 61 ! 62 character*20 :: fil_a2F, fil_info 63 ! 64 double precision, PARAMETER :: Rytocm1 = 109737.57990d0, & 65 & RytoGHz = 3.289828D6, & 66 & RytoTHz = RytoGHz/1000.d0, & 67 & RytoeV = 13.60569253d0, & 68 & tpi = 6.283185307d0, & 69 & convsig = 2.89d5, & 70 & KtoRy = 1.d0/38.681648/300/RytoeV 71 ! convsig: conversion factor for 72 ! conductivity. From au to (Ohm cm)-1 73 ! 74 namelist /input/ fil_info, fil_a2F, T, phband_i, phband_f,cut, & 75 & efermi, invtau, alat, vol, nthreads, lsoc, & 76 & cbm_i, shift, lscissors 77 78 read(5,input) 79 80 !Set number of threads 81 !$ call omp_set_num_threads(nthreads) 82 83 t0 = 0.0 84 85 npk = 40000 86 ! 87 ! Convert from eV to Ryd 88 T = T * KtoRy 89 invtau = invtau / RytoeV 90 efermi = efermi / RytoeV 91 ! 92 ! Total number of bands of interest (usually the number of relevant conduction/valence bands) 93 nphband = phband_f - phband_i + 1 94 ! 95 ! Read the a2Fsave file 96 open(11,file=fil_a2F,status='unknown') 97 ! 98 read(11,*) nbnd, nksfit 99 ! 100 allocate(etfit(nbnd,nksfit), xkfit(3,nksfit), wkfit(nksfit)) 101 ! 102 ! Read band structure and k-point data 103 read(11,*) etfit 104 read(11,*) ((xkfit(i,ik), i=1,3), ik=1,nksfit) 105 read(11,*) wkfit 106 read(11,*) nk1fit, nk2fit, nk3fit 107 read(11,* ) nsym 108 do ns=1,nsym 109 read(11,*) ((s(i,j,ns),j=1,3),i=1,3) 110 enddo 111 ! 112 close(11) 113 ! 114 ! Regular grid size 115 nkfit=nk1fit*nk2fit*nk3fit 116 ! 117 ! Read info file on k-points (and lattice) 118 open(11,file=fil_info,status='unknown') 119 ! 120 read(11,*) 121 read(11,*) ((at(i,j),i=1,3),j=1,3) 122 ! 123 read(11,*) 124 read(11,*) 125 ! 126 read(11,*) ((bg(i,j),i=1,3),j=1,3) 127 ! 128 close(11) 129 ! 130 allocate (eqkfit(nkfit),sfit(nkfit)) 131 allocate (iflag(nphband,nkfit),nkeff(nphband)) 132 allocate (dfk(nphband,nkfit,3),vk(nphband,nkfit,3)) 133 ! 134 ! Find the map between IBZ and the full-grid 135 call lint ( nsym, s, .true., at, bg, npk, 0,0,0, & 136 & nk1fit,nk2fit,nk3fit, nksfit, xkfit, 1, nkfit, eqkfit, sfit) 137 ! 138 ! Deallocate unnecessary variables 139 deallocate(wkfit,sfit,xkfit) 140 ! 141 ! Weights for regular grid 142 if ( lsoc .eqv. .true.) then 143 wk = 1.0/nkfit 144 else 145 wk = 2.0/nkfit 146 end if 147 ! 148 ! If lscissors is true, then shift band energies (move conduction bands higher in energy) 149 if (lscissors) then 150 etfit(cbm_i:nbnd,:) = etfit(cbm_i:nbnd,:) + shift/RytoeV 151 !cbm = cbm + shift 152 end if 153 ! 154 ! Call band velocities and forward derivatives 155 call vband_ibz( nk1fit,nk2fit,nk3fit,nphband,nksfit,etfit(phband_i:phband_f,:),eqkfit,at, vk, dfk) 156 ! 157 ! Include the 2pi/a factor 158 vk = vk / tpi * alat 159 ! 160 ! Determine the reduced grid for each band 161 call reducegrid(nkfit,nksfit,nphband,etfit(phband_i:phband_f,:),eqkfit,efermi, & 162 & cut,T, nkeff,iflag) 163 ! 164 ! Initiate Fermi integrals 165 I0 = 0.0 166 I1 = 0.0 167 I2 = 0.0 168 invT = 1.0/T 169 ! 170 !$ t0 = omp_get_wtime() 171 ! Loop over bands and kpoints 172 !$omp parallel default(shared) & 173 !$omp private(ibnd,ibnd_ph,ik,ikk,ind_k,i,j,fd,dfd,fac) 174 do ibnd=phband_i,phband_f 175 ! 176 ibnd_ph = ibnd - phband_i + 1 177 ! 178 !$omp do reduction(+: I0, I1, I2) 179 do ik=1,nkeff(ibnd_ph) 180 ! 181 ikk = iflag(ibnd_ph,ik) ! ikk is in full-grid (just reduced) 182 ind_k = eqkfit(ikk) ! ind_k is in IBZ 183 ! 184 ! Fermi factors 185 fac = etfit(ibnd,ind_k) - efermi 186 fd = 1.0/( exp(fac*invT) + 1.0 ) 187 dfd = invT * fd * (1.0 - fd) 188 ! 189 ! Compute I0 (related to Thomas-Fermi screening) 190 I0 = I0 + wk * dfd 191 ! 192 ! Compute I1, I2 193 do i=1,3 194 do j=1,3 195 I1(i,j) = I1(i,j) + wk * dfd / invtau * vk(ibnd_ph,ikk,i) & 196 & * vk(ibnd_ph,ikk,j) 197 ! 198 I2(i,j) = I2(i,j) + wk * dfd * fac / invtau * vk(ibnd_ph,ikk,i) & 199 & * vk(ibnd_ph,ikk,j) 200 ! 201 end do ! j 202 end do ! i 203 ! 204 end do ! ik 205 !$omp end do 206 ! 207 end do ! ibnd 208 !$omp end parallel 209 ! 210 !Total integration time 211 !$ t0 = omp_get_wtime() - t0 212 ! 213 ! Conductivity (convert to units of Ohm^-1 cm^-1) 214 sig = I1 * convsig / vol 215 ! 216 !Compute Seebeck 217 Se = 0.0 218 call inv33(I1, inv_I1) 219 do i=1,3 220 do j=1,3 221 do k=1,3 222 Se(i,j) = Se(i,j) + I2(i,k)*inv_I1(k,j) 223 end do 224 end do 225 end do 226 ! Units (convert to units of V/K) 227 Se = Se * (-RytoeV * KtoRy / T) 228 ! 229 ! Write D(Ef), Conductivity and Seebeck into file 230 open(11,file='sig_1.out',status='unknown') 231 open(12,file='Se_1.out',status='unknown') 232 open(13,file='Def_1.out',status='unknown') 233 write(11,"(10e14.6)") efermi * RytoeV, ((sig(i,j),i=1,3),j=1,3) 234 write(12,"(10e14.6)") efermi * RytoeV, ((Se(i,j),i=1,3),j=1,3) 235 write(13,"(2e14.6)") efermi * RytoeV, I0 / RytoeV 236 close(11) 237 close(12) 238 close(13) 239 ! 240 ! 241 open(11,file='report_1',status='unknown') 242 if ( t0 > 0 ) then 243 write(11,"(A,I2)") "Number of threads = ", nthreads 244 write(11,"(A,e14.6)") "Integration time(s) =", t0 245 end if 246 ! 247 write(11,"(A,I5)") "Number of kpoints (regular grid) = ", nkfit 248 do ibnd=1,nphband 249 write(11,"(A,2I5)") "band, reduced grid size", ibnd, nkeff(ibnd) 250 end do 251 close(11) 252 ! Free memory 253 deallocate(etfit,eqkfit,dfk,vk,nkeff,iflag) 254 255 contains 256 ! 257 ! Determinant of a 3x3 matrix 258 double precision function det33 (a) 259 260 implicit none 261 double precision, intent(in) :: a(3,3) 262 263 ! Determinant of a 264 det33 = a(1,1)*( a(2,2)*a(3,3)-a(2,3)*a(3,2) ) - & 265 & a(1,2)*( a(1,1)*a(3,3)-a(1,3)*a(3,1) ) + & 266 & a(1,3)*( a(2,1)*a(3,2)-a(2,2)*a(3,1) ) 267 268 end function det33 269 ! 270 ! 271 ! Inverse of a 3x3 matrix (analytical) 272 subroutine inv33(a, inv_a) 273 274 implicit none 275 276 double precision, intent(in) :: a(3,3) 277 double precision, intent(out) :: inv_a(3,3) 278 279 double precision :: deta_inv 280 281 deta_inv = 1.0 / det33(a) 282 283 inv_a(1,1) = deta_inv * ( a(2,2)*a(3,3)-a(2,3)*a(3,2) ) 284 inv_a(1,2) = deta_inv * ( a(1,3)*a(3,2)-a(1,2)*a(3,3) ) 285 inv_a(1,3) = deta_inv * ( a(1,2)*a(2,3)-a(1,3)*a(2,2) ) 286 287 inv_a(2,1) = deta_inv * ( a(2,3)*a(3,1)-a(2,1)*a(3,3) ) 288 inv_a(2,2) = deta_inv * ( a(1,1)*a(3,3)-a(1,3)*a(3,1) ) 289 inv_a(2,3) = deta_inv * ( a(1,3)*a(2,1)-a(1,1)*a(2,3) ) 290 291 inv_a(3,1) = deta_inv * ( a(2,1)*a(3,2)-a(2,2)*a(3,1) ) 292 inv_a(3,2) = deta_inv * ( a(1,2)*a(3,1)-a(1,1)*a(3,2) ) 293 inv_a(3,3) = deta_inv * ( a(1,1)*a(2,2)-a(1,2)*a(2,1) ) 294 295 end subroutine inv33 296 297 298 end program fermi_int_1 299 300