LCOV - code coverage report
Current view: top level - fermi - ef_newton.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 40 54 74.1 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             :       MODULE m_efnewton
       2             :       use m_juDFT
       3             :       CONTAINS
       4         122 :       RECURSIVE SUBROUTINE ef_newton(
       5             :      >     n,irank,
       6         122 :      >     inkem,nocst,index,tkb,e,
       7         122 :      X     w_near_ef,ef,we,recursion_level)
       8             : 
       9             : c***********************************************************************
      10             : c     
      11             : c     This subroutine adjusts the Fermi-Energy to obtain charge
      12             : c     neutrality. This is done by solving
      13             : c     
      14             : c     ( sum ff(e) * we ) - w_near_ef = 0
      15             : c     e
      16             : c     
      17             : c     using the Newton-Method.
      18             : c     Here the sum is over the eigenvalues between ef-8kt and ef+8kt, 
      19             : c     ff is the Fermi-Function, we is the weight of the according 
      20             : c     k-point and w_near_ef is the weight required between ef-8kt 
      21             : c     and ef+8kt to obtain neutrality.
      22             : c     
      23             : c***********************************************************************
      24             : 
      25             :       IMPLICIT NONE
      26             : 
      27             : C     .. Scalar Arguments ..
      28             :       INTEGER, INTENT (IN)    :: n
      29             :       INTEGER, INTENT (IN)    :: inkem,nocst,irank
      30             :       REAL,    INTENT (IN)    :: tkb
      31             :       REAL,    INTENT (INOUT) :: ef,w_near_ef
      32             :       INTEGER,INTENT(IN),OPTIONAL:: recursion_level
      33             : C     ..
      34             : C     .. Array Arguments ..
      35             :       INTEGER, INTENT (IN)    :: index(n)
      36             :       REAL,    INTENT (IN)    :: e(:) !(nkptd*neigd*jspd)
      37             :       REAL,    INTENT (INOUT) :: we(:) !(nkptd*neigd*jspd)
      38             : C     ..
      39             : C     .. Local Scalars ..
      40             :       REAL dff,expo,sdff,sff,efmin,efmax 
      41             :       INTEGER icnt,idim,rec_level
      42             :       LOGICAL efminok,efmaxok
      43             :       
      44             : C     ..
      45             : C     .. Local Arrays ..
      46         244 :       REAL ff(size(e))
      47             : C     ..
      48             : C     .. Parameters ..
      49             :       REAL,PARAMETER:: eps=1.0e-10
      50             : C     ..
      51             : c***********************************************************************
      52             : c     ABBREVIATIONS
      53             : c     
      54             : c     e          : linear list of the eigenvalues within the highest
      55             : c     energy-window
      56             : c     we         : list of weights of the eigenvalues in e
      57             : c     ef         : fermi energy determined to obtain charge neutrality
      58             : c     tkb        : value of temperature (kt) broadening around fermi
      59             : c     energy in htr units
      60             : c     index      : index list, e(index(n)) is the list of the 
      61             : c     eigenvalues in ascending order
      62             : c     ikem       : number of eigenvalues below ef-8kt
      63             : c     nocst      : number of eigenvalues below ef+8kt
      64             : c     w_near_ef  : weight (charge/spindg) between ef-8kt and ef+8kt
      65             : c     needed to obtain charge neutrality
      66             : c     
      67             : c***********************************************************************
      68             : 
      69         122 :       IF (PRESENT(recursion_level)) THEN
      70           0 :          rec_level=recursion_level+1
      71           0 :          IF (rec_level>20) CALL juDFT_error
      72             :      +     ("Determination of fermi-level did not converge",hint=
      73             :      + 'change temp. broad., set gauss = T, or use finer k mesh',
      74           0 :      + calledby="ef_newton")
      75             :       ELSE
      76         122 :          rec_level=0
      77         122 :          IF ( irank == 0 ) WRITE (6,FMT='(/,5x,''EF_NEWTON:  '',
      78         122 :      +''Adjust Fermi-Energy by Newton-Method.'',/)')
      79             :       ENDIF
      80             : c     
      81             : c     --->    NEWTON CYCLE
      82             : c
      83         122 :       efminok= .false. 
      84         122 :       efmaxok= .false. 
      85         333 :       DO icnt = 1,50
      86         333 :          sff = 0.0 
      87         333 :          sdff = 0.0
      88        2643 :          DO idim = inkem + 1,nocst
      89             : c     
      90             : c --->    COMPUTE FERMI-FUNCTION ff AND THE DERIVATIVE dff
      91             : c
      92        2310 :             expo= -ABS(e(index(idim))-ef)/tkb 
      93        2310 :             expo= EXP(expo)
      94        2310 :             IF (e(index(idim))<ef) THEN
      95        1220 :                ff(idim) = 1./ (expo+1.)
      96             :             ELSE
      97        1090 :                ff(idim)= expo/ (expo+1.)
      98             :             ENDIF 
      99        2310 :             dff= (1.+expo)**2 
     100        2310 :             dff= expo/(dff*tkb)
     101             : c     
     102             : c     --->    MULTIPLY WITH THE WEIGHT
     103             : c     
     104        2310 :             ff(idim) = we(index(idim))*ff(idim)
     105        2310 :             dff = we(index(idim))*dff
     106             : c     
     107             : c     --->    AND SUM IT UP
     108             : c     
     109        2310 :             sff = sff + ff(idim)
     110        2643 :             sdff = sdff + dff
     111             :          END DO
     112         333 :          sff = sff - w_near_ef
     113         333 :          IF (abs(sff).LT.eps) THEN
     114             :             !Converged, so do some output and return
     115         122 :             w_near_ef = sff + w_near_ef
     116         122 :             IF ( irank == 0 ) WRITE (6,FMT=8010) icnt,sff,-sff/sdff
     117             : 
     118        1696 :             DO idim = inkem + 1,nocst
     119         909 :                we(index(idim)) = ff(idim)
     120             :             END DO
     121             : 
     122             :  8000     FORMAT (15x,'ef_newton failed after      :',i3,'iterations.',/,
     123             :      +           15x,'The error in the weight is  : ',e12.5,/,
     124             :      +           15x,'The error in ef is          : ',e12.5,' htr',/)
     125             :  8010       FORMAT (15x,'Number of iterations needed  : ',i3,/,
     126             :      +           15x,'The error in the weight is   : ',e12.5,/,
     127             :      +           15x,'The error in ef is           : ',e12.5,' htr',/)
     128             :             RETURN 
     129             :          ENDIF
     130             : 
     131         211 :          IF (abs(sdff).LT.1e-29) THEN
     132           0 :             if (irank==0) THEN
     133           0 :                write(6,*) "Instability in determination of fermi-level,"
     134           0 :                write(6,*) "doubled temperature broading to continue"
     135           0 :                write(*,*) "Instability in determination of fermi-level,"
     136           0 :                write(*,*) "doubled temperature broading to continue"
     137             :             ENDIF
     138             :             CALL  ef_newton(
     139             :      >             n,irank,
     140             :      >             inkem,nocst,index,tkb*2.0,e,
     141           0 :      X             w_near_ef,ef,we,rec_level)
     142           0 :             RETURN
     143             :          ENDIF
     144         211 :          IF (sff > 0.) THEN
     145          94 :            efmax= ef 
     146          94 :            efmaxok= .true.  
     147             :          ELSE 
     148         117 :            efmin= ef 
     149         117 :            efminok= .true. 
     150             :          ENDIF
     151         211 :          ef = ef - sff/sdff
     152             :          ! if the Newton method is not stable, nested intervals are used 
     153         211 :          IF ( efminok .and. efmaxok ) THEN 
     154          14 :            IF ( (ef<efmin) .or. (ef>efmax) ) THEN 
     155           0 :              ef= (efmin+efmax)/2.
     156             :            ENDIF 
     157             :          ENDIF 
     158             :       END DO
     159             : c     
     160             : c---  > NOT CONVERGED AFTER 50 ITERATIONS
     161             : c     
     162           0 :       IF ( irank == 0 ) WRITE (6,FMT=8000) icnt,sff,-sff/sdff
     163           0 :       ef=ef+0.001
     164             :       CALL  ef_newton(
     165             :      >     n,irank,
     166             :      >     inkem,nocst,index,tkb*2.0,e,
     167           0 :      X     w_near_ef,ef,we,rec_level)
     168             :     
     169             :  
     170             :       END SUBROUTINE ef_newton
     171             :       END MODULE m_efnewton

Generated by: LCOV version 1.13