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