1
2! Copyright (C) 2018 T. Mueller, J. K. Dewhurst, S. Sharma and E. K. U. Gross.
3! This file is distributed under the terms of the GNU General Public License.
4! See the file COPYING for license details.
5
6subroutine occupyulr
7use modmain
8use modulr
9implicit none
10! local variables
11integer, parameter :: maxit=1000
12integer ik0,ik,ist,it
13real(8) e0,e1,e
14real(8) chg,x,t1
15! external functions
16real(8), external :: stheta
17! find minimum and maximum eigenvalues
18e0=evalu(1,1)
19e1=e0
20do ik0=1,nkpt0
21  do ist=1,nstulr
22    e=evalu(ist,ik0)
23    if (e.lt.e0) e0=e
24    if (e.gt.e1) e1=e
25  end do
26end do
27if (e0.lt.e0min) then
28  write(*,*)
29  write(*,'("Warning(occupyulr): minimum eigenvalue less than minimum &
30   &linearisation energy : ",2G18.10)') e0,e0min
31  write(*,'(" for s.c. loop ",I5)') iscl
32end if
33t1=1.d0/swidth
34! determine the Fermi energy using the bisection method
35do it=1,maxit
36  efermi=0.5d0*(e0+e1)
37  chg=0.d0
38  do ik0=1,nkpt0
39! central k-point
40    ik=(ik0-1)*nkpa+1
41    do ist=1,nstulr
42      e=evalu(ist,ik0)
43      if (e.lt.e0min) then
44        occulr(ist,ik0)=0.d0
45      else
46        x=(efermi-e)*t1
47        occulr(ist,ik0)=occmax*stheta(stype,x)
48        chg=chg+wkpt(ik)*occulr(ist,ik0)
49      end if
50    end do
51  end do
52  if (chg.lt.chgval) then
53    e0=efermi
54  else
55    e1=efermi
56  end if
57  if ((e1-e0).lt.1.d-12) return
58end do
59write(*,*)
60write(*,'("Warning(occupyulr): could not find Fermi energy")')
61end subroutine
62
63