1 2! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl. 3! This file is distributed under the terms of the GNU General Public License. 4! See the file COPYING for license details. 5 6!BOP 7! !ROUTINE: occupy 8! !INTERFACE: 9subroutine occupy 10! !USES: 11use modmain 12use modtest 13! !DESCRIPTION: 14! Finds the Fermi energy and sets the occupation numbers for the 15! second-variational states using the routine {\tt fermi}. 16! 17! !REVISION HISTORY: 18! Created February 2004 (JKD) 19! Added gap estimation, November 2009 (F. Cricchio) 20! Added adaptive smearing width, April 2010 (T. Bjorkman) 21!EOP 22!BOC 23implicit none 24! local variables 25integer, parameter :: maxit=1000 26integer ik,ist,it 27real(8) e0,e1,e 28real(8) chg,w,x,t1 29! external functions 30real(8), external :: sdelta,stheta 31! determine the smearing width automatically if required 32if ((autoswidth).and.(iscl.gt.1)) call findswidth 33! find minimum and maximum eigenvalues 34e0=evalsv(1,1) 35e1=e0 36do ik=1,nkpt 37 do ist=1,nstsv 38 e=evalsv(ist,ik) 39 if (e.lt.e0) e0=e 40 if (e.gt.e1) e1=e 41 end do 42end do 43if (e0.lt.e0min) then 44 write(*,*) 45 write(*,'("Warning(occupy): minimum eigenvalue less than minimum & 46 &linearisation energy : ",2G18.10)') e0,e0min 47 write(*,'(" for s.c. loop ",I5)') iscl 48end if 49t1=1.d0/swidth 50! determine the Fermi energy using the bisection method 51do it=1,maxit 52 efermi=0.5d0*(e0+e1) 53 chg=0.d0 54 do ik=1,nkpt 55 w=wkpt(ik) 56 do ist=1,nstsv 57 e=evalsv(ist,ik) 58 if (e.lt.e0min) then 59 occsv(ist,ik)=0.d0 60 else 61 x=(efermi-e)*t1 62 occsv(ist,ik)=occmax*stheta(stype,x) 63 chg=chg+w*occsv(ist,ik) 64 end if 65 end do 66 end do 67 if (chg.lt.chgval) then 68 e0=efermi 69 else 70 e1=efermi 71 end if 72 if ((e1-e0).lt.1.d-12) goto 10 73end do 74write(*,*) 75write(*,'("Warning(occupy): could not find Fermi energy")') 7610 continue 77! find the density of states at the Fermi surface in units of 78! states/Hartree/unit cell 79fermidos=0.d0 80do ik=1,nkpt 81 w=wkpt(ik) 82 do ist=1,nstsv 83 x=(evalsv(ist,ik)-efermi)*t1 84 fermidos=fermidos+w*sdelta(stype,x)*t1 85 end do 86 if (abs(occsv(nstsv,ik)).gt.epsocc) then 87 write(*,*) 88 write(*,'("Warning(occupy): not enough empty states for k-point ",I6)') ik 89 write(*,'(" and s.c. loop ",I5)') iscl 90 end if 91end do 92fermidos=fermidos*occmax 93! write Fermi density of states to test file 94call writetest(500,'DOS at Fermi energy',tol=5.d-3,rv=fermidos) 95! estimate the indirect band gap (FC) 96e0=-1.d8 97e1=1.d8 98ikgap(1)=1 99ikgap(2)=1 100! these loops are inefficiently ordered to fix a bug in versions 17 and 18 of 101! the Intel compiler 102do ist=1,nstsv 103 do ik=1,nkpt 104 e=evalsv(ist,ik) 105 if (e.lt.efermi) then 106 if (e.gt.e0) then 107 e0=e 108 ikgap(1)=ik 109 end if 110 else 111 if (e.lt.e1) then 112 e1=e 113 ikgap(2)=ik 114 end if 115 end if 116 end do 117end do 118bandgap(1)=e1-e0 119! write band gap to test file 120call writetest(510,'estimated indirect band gap',tol=2.d-2,rv=bandgap(1)) 121! estimate the direct band gap 122e=1.d8 123ikgap(3)=1 124do ik=1,nkpt 125 e0=-1.d8 126 e1=1.d8 127 do ist=1,nstsv 128 t1=evalsv(ist,ik) 129 if (t1.le.efermi) then 130 if (t1.gt.e0) e0=t1 131 else 132 if (t1.lt.e1) e1=t1 133 end if 134 end do 135 t1=e1-e0 136 if (t1.lt.e) then 137 e=t1 138 ikgap(3)=ik 139 end if 140end do 141bandgap(2)=e 142end subroutine 143!EOC 144 145