File: occupyulr.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (63 lines) | stat: -rw-r--r-- 1,397 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

! Copyright (C) 2018 T. Mueller, J. K. Dewhurst, S. Sharma and E. K. U. Gross.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

subroutine occupyulr
use modmain
use modulr
implicit none
! local variables
integer, parameter :: maxit=1000
integer ik,ist,it
real(8) e0,e1,e
real(8) chg,x,t1
! external functions
real(8) stheta
external stheta
! find minimum and maximum eigenvalues
e0=evalu(1,1)
e1=e0
do ik=1,nkpt0
  do ist=1,nstulr
    e=evalu(ist,ik)
    if (e.lt.e0) e0=e
    if (e.gt.e1) e1=e
  end do
end do
if (e0.lt.e0min) then
  write(*,*)
  write(*,'("Warning(occupyulr): minimum eigenvalue less than minimum &
   &linearisation energy : ",2G18.10)') e0,e0min
  write(*,'(" for s.c. loop ",I5)') iscl
end if
t1=1.d0/swidth
! determine the Fermi energy using the bisection method
do it=1,maxit
  efermi=0.5d0*(e0+e1)
  chg=0.d0
  do ik=1,nkpt0
    do ist=1,nstulr
      e=evalu(ist,ik)
      if (e.lt.e0min) then
        occulr(ist,ik)=0.d0
      else
        x=(efermi-e)*t1
        occulr(ist,ik)=occmax*stheta(stype,x)
        chg=chg+wkpt0(ik)*occulr(ist,ik)
      end if
    end do
  end do
  if (chg.lt.chgvalu) then
    e0=efermi
  else
    e1=efermi
  end if
  if ((e1-e0).lt.1.d-12) goto 10
end do
write(*,*)
write(*,'("Warning(occupyulr): could not find Fermi energy")')
10 continue
return
end subroutine