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