ef_newton.f 5.56 KB
Newer Older
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
      MODULE m_efnewton
      use m_juDFT
      CONTAINS
      RECURSIVE SUBROUTINE ef_newton(
     >     n,irank,
     >     inkem,nocst,index,tkb,e,
     X     w_near_ef,ef,we,recursion_level)

c***********************************************************************
c     
c     This subroutine adjusts the Fermi-Energy to obtain charge
c     neutrality. This is done by solving
c     
c     ( sum ff(e) * we ) - w_near_ef = 0
c     e
c     
c     using the Newton-Method.
c     Here the sum is over the eigenvalues between ef-8kt and ef+8kt, 
c     ff is the Fermi-Function, we is the weight of the according 
c     k-point and w_near_ef is the weight required between ef-8kt 
c     and ef+8kt to obtain neutrality.
c     
c***********************************************************************

      IMPLICIT NONE

C     .. Scalar Arguments ..
      INTEGER, INTENT (IN)    :: n
      INTEGER, INTENT (IN)    :: inkem,nocst,irank
      REAL,    INTENT (IN)    :: tkb
      REAL,    INTENT (INOUT) :: ef,w_near_ef
      INTEGER,INTENT(IN),OPTIONAL:: recursion_level
C     ..
C     .. Array Arguments ..
      INTEGER, INTENT (IN)    :: index(n)
      REAL,    INTENT (IN)    :: e(:) !(nkptd*neigd*jspd)
      REAL,    INTENT (INOUT) :: we(:) !(nkptd*neigd*jspd)
C     ..
C     .. Local Scalars ..
40
      REAL dff,expo,sdff,sff,efmin,efmax 
41
      INTEGER icnt,idim,rec_level
42 43
      LOGICAL efminok,efmaxok
      
44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
C     ..
C     .. Local Arrays ..
      REAL ff(size(e))
C     ..
C     .. Parameters ..
      REAL,PARAMETER:: eps=1.0e-10
C     ..
c***********************************************************************
c     ABBREVIATIONS
c     
c     e          : linear list of the eigenvalues within the highest
c     energy-window
c     we         : list of weights of the eigenvalues in e
c     ef         : fermi energy determined to obtain charge neutrality
c     tkb        : value of temperature (kt) broadening around fermi
c     energy in htr units
c     index      : index list, e(index(n)) is the list of the 
c     eigenvalues in ascending order
c     ikem       : number of eigenvalues below ef-8kt
c     nocst      : number of eigenvalues below ef+8kt
c     w_near_ef  : weight (charge/spindg) between ef-8kt and ef+8kt
c     needed to obtain charge neutrality
c     
c***********************************************************************

      IF (PRESENT(recursion_level)) THEN
         rec_level=recursion_level+1
         IF (rec_level>20) CALL juDFT_error
     +     ("Determination of fermi-level did not converge",hint
     +        ='change temperature broad. or set gauss = T',calledby
     +        ="ef_newton")
      ELSE
         rec_level=0
         IF ( irank == 0 ) WRITE (6,FMT='(/,5x,''EF_NEWTON:  '',
     +''Adjust Fermi-Energy by Newton-Method.'',/)')
      ENDIF
c     
c     --->    NEWTON CYCLE
82 83 84
c
      efminok= .false. 
      efmaxok= .false. 
85 86 87 88 89
      DO icnt = 1,50
         sff = 0.0 
         sdff = 0.0
         DO idim = inkem + 1,nocst
c     
90 91 92 93 94 95 96 97 98 99 100
c --->    COMPUTE FERMI-FUNCTION ff AND THE DERIVATIVE dff
c
            expo= -ABS(e(index(idim))-ef)/tkb 
            expo= EXP(expo)
            IF (e(index(idim))<ef) THEN
               ff(idim) = 1./ (expo+1.)
            ELSE
               ff(idim)= expo/ (expo+1.)
            ENDIF 
            dff= (1.+expo)**2 
            dff= expo/(dff*tkb)
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
c     
c     --->    MULTIPLY WITH THE WEIGHT
c     
            ff(idim) = we(index(idim))*ff(idim)
            dff = we(index(idim))*dff
c     
c     --->    AND SUM IT UP
c     
            sff = sff + ff(idim)
            sdff = sdff + dff
         END DO
         sff = sff - w_near_ef
         IF (abs(sff).LT.eps) THEN
            !Converged, so do some output and return
            w_near_ef = sff + w_near_ef
            IF ( irank == 0 ) WRITE (6,FMT=8010) icnt,sff,-sff/sdff

            DO idim = inkem + 1,nocst
               we(index(idim)) = ff(idim)
            END DO

122
 8000     FORMAT (15x,'ef_newton failed after      :',i3,'iterations.',/,
123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143
     +           15x,'The error in the weight is  : ',e12.5,/,
     +           15x,'The error in ef is          : ',e12.5,' htr',/)
 8010       FORMAT (15x,'Number of iterations needed  : ',i3,/,
     +           15x,'The error in the weight is   : ',e12.5,/,
     +           15x,'The error in ef is           : ',e12.5,' htr',/)
            RETURN 
         ENDIF

         IF (abs(sdff).LT.1e-29) THEN
            if (irank==0) THEN
               write(6,*) "Instability in determination of fermi-level,"
               write(6,*) "doubled temperature broading to continue"
               write(*,*) "Instability in determination of fermi-level,"
               write(*,*) "doubled temperature broading to continue"
            ENDIF
            CALL  ef_newton(
     >             n,irank,
     >             inkem,nocst,index,tkb*2.0,e,
     X             w_near_ef,ef,we,rec_level)
            RETURN
         ENDIF
144 145 146 147 148 149 150
         IF (sff > 0.) THEN
           efmax= ef 
           efmaxok= .true.  
         ELSE 
           efmin= ef 
           efminok= .true. 
         ENDIF
151
         ef = ef - sff/sdff
152 153 154 155 156 157
         ! if the Newton method is not stable, nested intervals are used 
         IF ( efminok .and. efmaxok ) THEN 
           IF ( (ef<efmin) .or. (ef>efmax) ) THEN 
             ef= (efmin+efmax)/2.
           ENDIF 
         ENDIF 
158 159 160 161 162 163 164 165 166 167 168 169 170 171
      END DO
c     
c---  > NOT CONVERGED AFTER 50 ITERATIONS
c     
      IF ( irank == 0 ) WRITE (6,FMT=8000) icnt,sff,-sff/sdff
      ef=ef+0.001
      CALL  ef_newton(
     >     n,irank,
     >     inkem,nocst,index,tkb*2.0,e,
     X     w_near_ef,ef,we,rec_level)
    
 
      END SUBROUTINE ef_newton
      END MODULE m_efnewton