The Elk Code
occupy.f90
Go to the documentation of this file.
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:
9 subroutine occupy
10 ! !USES:
11 use modmain
12 use modtest
13 ! !DESCRIPTION:
14 ! Finds the Fermi energy and sets the occupation numbers for the
15 ! second-variational states using the routine {\tt fermi}. Also determines
16 ! the density of states at the Fermi surface as well as the direct and
17 ! indirect band gaps.
18 !
19 ! !REVISION HISTORY:
20 ! Created February 2004 (JKD)
21 ! Added gap estimation, November 2009 (F. Cricchio)
22 ! Added adaptive smearing width, April 2010 (T. Bjorkman)
23 !EOP
24 !BOC
25 implicit none
26 ! local variables
27 integer, parameter :: maxit=1000
28 integer ik,ist,it
29 real(8) e0,e1,e,de
30 real(8) ei0,ei1,ed0,ed1
31 real(8) chg,w,x,t0
32 ! external functions
33 real(8), external :: sdelta,stheta
34 ! determine the smearing width automatically if required
35 if ((autoswidth).and.(iscl > 1)) call findswidth
36 ! find minimum and maximum eigenvalues
37 e0=evalsv(1,1)
38 e1=e0
39 do ik=1,nkpt
40  do ist=1,nstsv
41  e=evalsv(ist,ik)
42  if (e < e0) e0=e
43  if (e > e1) e1=e
44  end do
45 end do
46 if (e0 < e0min) then
47  write(*,*)
48  write(*,'("Warning(occupy): minimum eigenvalue less than minimum &
49  &linearisation energy : ",2G18.10)') e0,e0min
50  write(*,'(" for s.c. loop ",I5)') iscl
51 end if
52 t0=1.d0/swidth
53 ! determine the Fermi energy using the bisection method
54 do it=1,maxit
55  efermi=0.5d0*(e0+e1)
56  chg=0.d0
57  do ik=1,nkpt
58  w=wkpt(ik)
59  do ist=1,nstsv
60  e=evalsv(ist,ik)
61  if (e < e0min) then
62  occsv(ist,ik)=0.d0
63  else
64  x=(efermi-e)*t0
65  occsv(ist,ik)=occmax*stheta(stype,x)
66  chg=chg+w*occsv(ist,ik)
67  end if
68  end do
69  end do
70  if (chg < chgval) then
71  e0=efermi
72  else
73  e1=efermi
74  end if
75  if ((e1-e0) < 1.d-12) goto 10
76 end do
77 write(*,*)
78 write(*,'("Warning(occupy): could not find Fermi energy")')
79 10 continue
80 if (any(abs(occsv(nstsv,1:nkpt)) > epsocc)) then
81  write(*,*)
82  write(*,'("Warning(occupy): not enough empty states for s.c. loop ",I5)') iscl
83 end if
84 ! find the density of states at the Fermi surface in units of
85 ! states/Hartree/unit cell
86 fermidos=0.d0
87 do ik=1,nkpt
88  w=wkpt(ik)
89  do ist=1,nstsv
90  x=(evalsv(ist,ik)-efermi)*t0
91  fermidos=fermidos+w*sdelta(stype,x)
92  end do
93 end do
95 ! write Fermi density of states to test file
96 call writetest(500,'DOS at Fermi energy',tol=5.d-3,rv=fermidos)
97 ! estimate the indirect and direct band gaps (FC)
98 ei0=-1.d8; ei1=1.d8
99 de=1.d8
100 ikgap(1:3)=1
101 do ik=1,nkpt
102  ed0=-1.d8; ed1=1.d8
103  do ist=1,nstsv
104  e=evalsv(ist,ik)
105  if (e <= efermi) then
106  if (e > ed0) ed0=e
107  if (e > ei0) then
108 ! transfer is a workaround for a bug in Intel Fortran versions 17 and 18
109  ikgap(1)=transfer(ik,ik)
110  ei0=e
111  end if
112  else
113  if (e < ed1) ed1=e
114  if (e < ei1) then
115  ikgap(2)=ik
116  ei1=e
117  end if
118  end if
119  end do
120  e=ed1-ed0
121  if (e < de) then
122  ikgap(3)=ik
123  de=e
124  end if
125 end do
126 bandgap(1)=ei1-ei0
127 bandgap(2)=de
128 ! automatically find the difference between the fixed linearisation and Fermi
129 ! energies if required
130 if (autodlefe) call finddlefe
131 ! write band gap to test file
132 call writetest(510,'estimated indirect band gap',tol=2.d-2,rv=bandgap(1))
133 end subroutine
134 !EOC
135 
real(8) efermi
Definition: modmain.f90:907
subroutine writetest(id, descr, nv, iv, iva, tol, rv, rva, zv, zva)
Definition: modtest.f90:16
subroutine finddlefe
Definition: finddlefe.f90:7
integer, dimension(3) ikgap
Definition: modmain.f90:917
real(8), dimension(:,:), allocatable evalsv
Definition: modmain.f90:921
subroutine occupy
Definition: occupy.f90:10
integer nkpt
Definition: modmain.f90:461
integer iscl
Definition: modmain.f90:1051
subroutine findswidth
Definition: findswidth.f90:10
real(8) swidth
Definition: modmain.f90:895
real(8) epsocc
Definition: modmain.f90:903
real(8) fermidos
Definition: modmain.f90:913
real(8) e0min
Definition: modmain.f90:825
integer nstsv
Definition: modmain.f90:889
real(8), dimension(:), allocatable wkpt
Definition: modmain.f90:475
real(8) occmax
Definition: modmain.f90:901
real(8), dimension(:,:), allocatable occsv
Definition: modmain.f90:905
integer stype
Definition: modmain.f90:891
real(8) chgval
Definition: modmain.f90:722
logical autodlefe
Definition: modmain.f90:833
logical autoswidth
Definition: modmain.f90:897
real(8), dimension(2) bandgap
Definition: modmain.f90:915