1 program efermi 2 !--------------------------------------------------------------------------------- 3 ! Written by Burak Himmetoglu (burakhmmtgl@gmail.com) 4 ! Uses some parts of the PW distribution 5 ! 6 ! This is a simple code that reads the electron concentration and band 7 ! structure (from a2Fsave file) and then computes the Fermi level. 8 ! The system must be insulating. 9 ! 10 ! Description of the input card: 11 ! 12 ! fil_a2F : prefix.a2Fsave file 13 ! fil_info : file containing direct and reciprocal lattice vectors 14 ! T : Temperature in K 15 ! vol : Volume in au^3 16 ! cbm, vbm : Conduction band min. and Valence band max. (in eV) 17 ! doping : Electron concentration (in cm^-3). For hole doping, use negative value 18 ! ndop : Number of doping levels considered 19 ! nthreads : Number of threads for OpenMP parallelization 20 ! lscissors : if .true. conduction bands are shifted by a constant energy up 21 ! shift : value of the shift applied to conduction bands 22 ! cbm_i : the initial conduction band for which the shift is applied ( all bands with index >= cbm_i are shifted up) 23 24!$ use omp_lib 25 26 implicit none 27 ! 28 integer :: i,j,iq,ik,imu,itemp,nu,nqtot,nsig,nat,nk1fit,nk2fit, & 29 & nk3fit, nkfit, nksfit_real, nbnd, nksfit, npk, nsym, & 30 & s(3,3,48),ns,nrot,ibnd,ndop,cbm_i 31 ! 32 ! OMP 33 integer :: TID, nthreads 34 double precision :: t0 35 36 ! 37 double precision :: at(3,3), bg(3,3),vol,T, fd, nelec1, & 38 & cbm, vbm, ef_mid, doping(10), nelec, & 39 & Ef_IBZ(10), shift 40 ! 41 double precision, allocatable :: xkfit(:,:), etfit(:,:), wkfit(:),& 42 & wk(:) 43 ! 44 integer, allocatable :: eqkfit(:), sfit(:) 45 ! 46 logical :: lscissors 47 ! 48 character*20 :: fil_kp, fil_a2F, fil_info 49 ! 50 double precision, PARAMETER :: Rytocm1 = 109737.57990d0, & 51 & RytoGHz = 3.289828D6, & 52 & RytoTHz = RytoGHz/1000.d0, & 53 & RytoeV = 13.60569253d0, & 54 & tpi = 6.283185307d0, & 55 & convsig = 2.89d5, & 56 & KtoRy = 1.d0/38.681648/300/RytoeV, & 57 & autocm = 5.2917721092d-9, & 58 & convRH = 9.2522431d-13 59 ! convsig: conversion factor for 60 ! conductivity. From au to (Ohm cm)-1 61 62 namelist /input/ fil_info, fil_a2F, vol, cbm, vbm, T, & 63 & doping, ndop, nthreads, lscissors, shift, cbm_i 64 65 read(5,input) 66 67 t0 = 0.0 68 69 !$ call omp_set_num_threads(nthreads) 70 71 !$ t0 = omp_get_wtime() 72 npk = 40000 73 ! 74 ! Convert temperature from K to Ryd 75 T = T * KtoRy 76 ! 77 ! Read a2F 78 ! 79 open(11,file=fil_a2F,status='unknown') 80 ! 81 read(11,*) nbnd, nksfit 82 ! 83 allocate(etfit(nbnd,nksfit), xkfit(3,nksfit), wkfit(nksfit)) 84 ! 85 ! Read band structure and k-point data 86 read(11,*) etfit 87 read(11,*) ((xkfit(i,ik), i=1,3), ik=1,nksfit) 88 read(11,*) wkfit 89 read(11,*) nk1fit, nk2fit, nk3fit 90 read(11,* ) nsym 91 do ns=1,nsym 92 read(11,*) ((s(i,j,ns),j=1,3),i=1,3) 93 enddo 94 ! 95 close(11) 96 ! 97 ! Regular grid size 98 nkfit=nk1fit*nk2fit*nk3fit 99 ! Read info file on k-points (and lattice) 100 ! 101 open(11,file=fil_info,status='unknown') 102 ! 103 read(11,*) 104 read(11,*) ((at(i,j),i=1,3),j=1,3) 105 ! 106 read(11,*) 107 read(11,*) 108 ! 109 read(11,*) ((bg(i,j),i=1,3),j=1,3) 110 ! 111 close(11) 112 ! 113 allocate (eqkfit(nkfit), sfit(nkfit),wk(nkfit)) 114 ! eqkfit : pointers to band energies in uniform grid 115 ! sfit : pointers to symmetries 116 ! 117 call lint ( nsym, s, .true., at, bg, npk, 0,0,0, & 118 & nk1fit,nk2fit,nk3fit, nksfit, xkfit, 1, nkfit, eqkfit, sfit) 119 ! 120 ! eqkfit(nk) : maps IBZ to full grid. The full grid is in crystal coords 121 ! 122 ! If lscissors is true, then shift band energies (move conduction bands higher in energy) 123 if (lscissors) then 124 etfit(cbm_i:nbnd,:) = etfit(cbm_i:nbnd,:) + shift/RytoeV 125 cbm = cbm + shift 126 end if 127 ! 128 ! Determine the number of electrons (neutral, insulating system) 129 ef_mid = (cbm + vbm)/2.0/RytoeV ! Fermi level at the middle of the gap 130 ! 131 nelec1 = 0.d0 132 ! 133 ! Irreducable BZ 134 nelec1 = sumk (etfit,nbnd,nksfit,wkfit,T,1,ef_mid) 135 ! 136 do i=1,ndop 137 ! Determine number of electrons from doping levels (given in cm-3) 138 nelec = doping(i) * vol * autocm**3 + nelec1 139 ! 140 ! Compute efermi using the bisection method 141 Ef_IBZ(i) = fermi_en (etfit,wkfit,nbnd,nksfit,nelec,T,1) 142 ! 143 end do 144 ! 145 ! Memory clean 146 ! 147 deallocate (xkfit,etfit,wkfit,wk,eqkfit,sfit) 148 ! 149 !$t0 = omp_get_wtime() - t0 150 ! 151 if (t0 >0) write(6,"(A,e14.6)") 'Walltime= ', t0 152 ! 153 do i=1,ndop 154 write(6,"(A,e14.6,f14.6)") 'doping, Efermi', doping(i), Ef_IBZ(i) * RytoeV 155 end do 156 ! 157 ! 158 contains 159 ! ----------------------------------------------- 160 ! FUNCTIONS 161 ! ----------------------------------------------- 162 163 double precision function sumk (et,nbnd,nks,wk,degauss,ngauss,ee) 164 165!$ use omp_lib 166 use kinds 167 168 implicit none 169 170 integer, intent(in) :: nbnd, nks, ngauss 171 double precision, intent(in) :: ee,degauss 172 ! wk: weight of kpt 173 ! ee: E 174 ! degauss: temperature/degauss 175 double precision, intent(in) :: et(nbnd,nks), wk(nks) 176 177 double precision :: fd, x, arg, a, hp, hd, sqrtpm1 178 179 integer :: i, j, ik, n, ni 180 181 real(dp), external :: qe_erf 182 183 sqrtpm1 = 1.0d0/1.77245385090551602729d0 184 185 sumk = 0.d0 186 187 n = 1 ! First order Gauss-Hermite 188 189 if (ngauss .eq. 1) then ! FD 190 !$omp parallel do & 191 !$omp default(shared) & 192 !$omp private(i,ik,fd) & 193 !$omp reduction(+ : sumk) 194 do i=1,nbnd 195 do ik=1,nks 196 fd = 1.d0 / ( 1.d0+exp( (et(i,ik)-ee)/degauss ) ) 197 sumk = sumk + wk(ik) * fd 198 end do 199 end do 200 !$omp end parallel do 201 else if (ngauss .eq. 0) then ! Gaussian (for zero temperature) 202 !$omp parallel do & 203 !$omp default(shared) & 204 !$omp private(i,ik,fd,x) & 205 !$omp reduction(+ : sumk) 206 do i=1,nbnd 207 do ik=1,nks 208 x = (ee-et(i,ik))/degauss/dsqrt(2.d0) 209 fd = 0.5d0 * (1.d0 + qe_erf(x) ) 210 sumk = sumk + wk(ik) * fd 211 end do 212 end do 213 !$omp end parallel do 214 else if ( ngauss .eq. 2 ) then ! MP 215 !$omp parallel do & 216 !$omp default(shared) & 217 !$omp private(i,j,ik,fd,x) & 218 !$omp reduction(+ : sumk) 219 do i=1,nbnd 220 do ik=1,nks 221 ! 222 x = (ee-et(i,ik))/degauss/dsqrt(2.d0) 223 fd = 0.5d0 * (1.d0 + qe_erf(x) ) 224 hd = 0.d0 225 arg = min (200.d0, x**2) 226 hp = exp ( - arg) 227 ni = 0 228 a = sqrtpm1 229 ! 230 do j = 1, n 231 hd = 2.0d0 * x * hp - 2.0d0 * DBLE (ni) * hd 232 ni = ni + 1 233 a = - a / (DBLE (j) * 4.0d0) 234 fd = fd - a * hd 235 hp = 2.0d0 * x * hd-2.0d0 * DBLE (ni) * hp 236 ni = ni + 1 237 end do 238 ! 239 sumk = sumk + wk(ik) * fd 240 end do 241 end do 242 !$omp end parallel do 243 244 end if 245 246 end function sumk 247 ! 248 double precision function fermi_en (eig,wk,nbnd,nktot,nelec,temperature,ngauss) 249 250 implicit none 251 252 integer, intent(in) :: nbnd, nktot, ngauss 253 double precision, intent(in) :: eig(nbnd,nktot),wk(nktot), nelec, & 254 & temperature 255 256 ! Local variables 257 integer :: i, ik 258 double precision :: Elw, Eup, Ef, sumkup, sumklw, sumkmid 259 integer, PARAMETER :: maxiter = 300 260 double precision, PARAMETER :: eps = 1d-10 261 262 Elw = eig(1,1) 263 Eup = eig(nbnd,1) 264 do ik = 2, nktot 265 Elw = min ( Elw, eig (1, ik) ) 266 Eup = max ( Eup, eig (nbnd, ik) ) 267 end do 268 Eup = Eup + 2*temperature 269 Elw = Elw - 2*temperature 270 271 sumkup = sumk(eig, nbnd, nktot, wk, temperature, ngauss, Eup) 272 sumklw = sumk(eig, nbnd, nktot, wk, temperature, ngauss, Elw) 273 ! 274 if ( (sumkup - nelec) < -eps .or. (sumklw - nelec) > eps ) then 275 fermi_en = -99 276 return 277 end if 278 ! 279 do 100 i=1, maxiter 280 Ef = 0.5d0 * (Eup + Elw) 281 sumkmid = sumk(eig, nbnd, nktot, wk, temperature, ngauss, Ef) 282 if ( abs (sumkmid-nelec) < eps ) then 283 fermi_en = Ef 284 exit 285 else if ( (sumkmid-nelec) < -eps) then 286 Elw = Ef 287 else 288 Eup = Ef 289 end if 290100 continue 291 292 end function fermi_en 293 ! 294 end program efermi 295