LCOV - code coverage report
Current view: top level - fermi - ferhis.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 73 80 91.2 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_ferhis
       8             : CONTAINS
       9         501 :   SUBROUTINE ferhis(input,kpts,mpi, index,idxeig,idxkpt,idxjsp,n,&
      10         501 :        nstef,ws,spindg,weight, e,ne,we, noco,cell,ef,seigv,w_iks,results)
      11             :     !***********************************************************************
      12             :     !
      13             :     !     This subroutine determines the fermi energy and the sum of the
      14             :     !     single particle eigenvalues by histogram method.
      15             :     !
      16             :     !
      17             :     !     Theory :   zelec(nwd) = sum{ sum{ sum{ we * f(e) } } }
      18             :     !                             sp    e    k
      19             :     !
      20             :     !
      21             :     !                seigv = sum{ sum{ sum{ w(k) * f(e) * e }
      22             :     !                         sp   e    k
      23             :     !
      24             :     !                the weights w(k) are normalized: sum{w(k)} = 1
      25             :     !                                                  k                -6
      26             :     !                         a) 1                           for kT < 10
      27             :     !                we    = {                           -1             -6
      28             :     !                         b){ 1+exp(e(k,nu) -ef)/kt) }   for kt >=10
      29             :     !
      30             :     !                in case a) we choose the Fermi energy the highest
      31             :     !                           valence state
      32             :     !                        b) we choose as Fermi energy the arithmetric
      33             :     !                           mean between the highest occupied and lowest
      34             :     !                           unoccupied state if this energy difference
      35             :     !                           Delta E <= kT, otherwise as a).
      36             :     !
      37             :     !                                      stefan bl"ugel , kfa , oct 1991
      38             :     !
      39             :     !               free energy and extrapolation T -> 0  implemented
      40             :     !                         (see M.J.Gillan, J.Phys.: Condens. Matter 1,
      41             :     !                          (1989) 689-711 )
      42             :     !
      43             :     !                                      peter richard, kfa, jun. 1995
      44             :     !
      45             :     !               adapted to flapw7
      46             :     !
      47             :     !                                      philipp kurz, kfa, oct. 1995
      48             :     !               entropy calculation changed
      49             :     !                     
      50             :     !                                      r.pentcheva, kfa, may  1996
      51             :     !
      52             :     !***********************************************************************
      53             :     USE m_efnewton
      54             :     USE m_types
      55             :     USE m_xmlOutput
      56             :     USE m_constants
      57             :     IMPLICIT NONE
      58             :     TYPE(t_results),INTENT(INOUT)   :: results
      59             :     TYPE(t_mpi),INTENT(IN)          :: mpi
      60             :     TYPE(t_input),INTENT(IN)        :: input
      61             :     TYPE(t_kpts),INTENT(IN)         :: kpts
      62             :     TYPE(t_noco),INTENT(IN),OPTIONAL         :: noco
      63             :     TYPE(t_cell),INTENT(IN),OPTIONAL         :: cell
      64             :     !     ..
      65             :     !     .. Scalar Arguments ..
      66             :     INTEGER,INTENT(IN)  ::  n ,nstef
      67             :     REAL,INTENT(IN)     ::  spindg,ws,weight
      68             :     REAL,INTENT(INOUT)  ::  ef,seigv
      69             :     REAL,INTENT(OUT)    ::  w_iks(:,:,:)
      70             :     !     ..
      71             :     !     .. Array Arguments ..
      72             :     INTEGER, INTENT (IN) :: idxeig(:)!(dimension%neigd*kpts%nkpt*dimension%jspd)
      73             :     INTEGER, INTENT (IN) :: idxjsp(:)!(dimension%neigd*kpts%nkpt*dimension%jspd)
      74             :     INTEGER, INTENT (IN) :: idxkpt(:)!(dimension%neigd*kpts%nkpt*dimension%jspd)
      75             :     INTEGER, INTENT (IN) ::  INDEX(:)!(dimension%neigd*kpts%nkpt*dimension%jspd)
      76             :     INTEGER, INTENT (IN) ::     ne(:,:)!(kpts%nkpt,dimension%jspd)
      77             :     REAL,    INTENT (IN) ::      e(:)!(kpts%nkpt*dimension%neigd*dimension%jspd)
      78             :     REAL,    INTENT (INOUT) ::  we(:)!(kpts%nkpt*dimension%neigd*dimension%jspd)
      79             : 
      80             :     !--- J constants
      81             :     !--- J constants
      82             : 
      83             :     !     ..
      84             :     !     .. Local Scalars ..
      85             :     REAL,PARAMETER:: del=1.e-6
      86             :     REAL :: efermi,emax,emin,entropy,fermikn,gap,&
      87             :               wfermi,wvals,w_below_emin,w_near_ef,tkb
      88             :     INTEGER ink,inkem,j,js,k,kpt,nocc,nocst,i,nspins
      89             : 
      90             :     !     .. Local Arrays ..      
      91             :     REAL :: qc(3)
      92             :     CHARACTER(LEN=20)    :: attributes(2)
      93             : 
      94             :     !     ..
      95             :     !***********************************************************************
      96             :     !------->          ABBREVIATIONS
      97             :     !
      98             :     !     eig        : array of eigenvalues within all energy-windows
      99             :     !     wtkpt      : list of the weights of each k-point (from inp-file)
     100             :     !     e          : linear list of the eigenvalues within the highest
     101             :     !                  energy-window
     102             :     !     we         : list of weights of the eigenvalues in e
     103             :     !     w          : array of weights (output, needed to calculate the
     104             :     !                  new charge-density)
     105             :     !     zelec      : number of electrons in a window
     106             :     !     spindg     : spindegeneracy (2 in nonmagnetic calculations)
     107             :     !     seigv      : weighted sum of the occupied valence eigenvalues
     108             :     !     seigsc     : weighted sum of the semi-core eigenvalues
     109             :     !     seigscv    : sum of seigv and seigsc
     110             :     !     ts         : entropy contribution to the free energy
     111             :     !     tkb        : value of temperature (kt) broadening around fermi
     112             :     !                  energy in htr units
     113             :     !     ef         : fermi energy determined to obtain charge neutrality
     114             :     !     wfermi     : weight of the occupation number for states at the
     115             :     !                  fermi energy.
     116             :     !     fd         : fermi dirac distribution
     117             :     !     fermikn    : fermi dirac distribution for the k-th point 
     118             :     !                  and n-th state
     119             :     !**********************************************************************
     120             :     !     ..
     121         167 :     nspins=input%jspins
     122         167 :     if (noco%l_noco) nspins=1
     123         167 :     tkb=input%tkb !might be modified if we have an insulator
     124         167 :     IF ( mpi%irank == 0 ) THEN
     125         167 :        WRITE (6,FMT='(/)')
     126         167 :        WRITE (6,FMT='(''FERHIS:  Fermi-Energy by histogram:'')')
     127             :     END IF
     128             : 
     129         167 :     efermi = ef
     130         167 :     IF (nstef.LT.n) THEN
     131         167 :        gap = e(INDEX(nstef+1)) - ef
     132         167 :        results%bandgap = gap*hartree_to_ev_const
     133         167 :        IF ( mpi%irank == 0 ) THEN
     134         501 :           attributes = ''
     135         167 :           WRITE(attributes(1),'(f20.10)') gap*hartree_to_ev_const
     136         167 :           WRITE(attributes(2),'(a)') 'eV'
     137         167 :           CALL writeXMLElement('bandgap',(/'value','units'/),attributes)
     138         167 :           WRITE (6,FMT=8050) gap
     139             :        END IF
     140             :     END IF
     141         167 :     IF ( mpi%irank == 0 ) THEN
     142         167 :        WRITE ( 6,FMT=8010) spindg* (ws-weight)
     143             :     END IF
     144             :     !
     145             :     !---> DETERMINE OCCUPATION AT THE FERMI LEVEL
     146             :     !
     147         167 :     wfermi = ws - weight
     148             :     !                                          -6
     149             :     !======> DETERMINE FERMI ENERGY for kT >= 10
     150             :     !
     151             :     !
     152         167 :     IF (tkb.GE.del) THEN
     153             :        !
     154             :        !---> TEMPERATURE BROADENING
     155             :        !
     156         167 :        IF (nstef.LT.n) THEN
     157             :           !
     158             :           !--->    STATES ABOVE EF AVAILABLE           
     159             :           !
     160         167 :           ef = 0.5* (e(INDEX(nstef+1))+ef)
     161         167 :           emax = ef + 8.0*tkb
     162         167 :           emin = ef - 8.0*tkb
     163         167 :           w_near_ef = 0.0
     164         167 :           w_below_emin = 0.0
     165         167 :           inkem = 0
     166       15942 :           ink_loop: DO ink = 1,n
     167             : 
     168       16109 :              IF (e(INDEX(ink)).LT.emin) THEN
     169       14988 :                 inkem = ink
     170       14988 :                 w_below_emin = w_below_emin + we(INDEX(ink))
     171         954 :              ELSE IF (e(INDEX(ink)).GT.emax) THEN
     172             :                 EXIT ink_loop
     173             :              END IF
     174             : 
     175             :           ENDDO ink_loop
     176         167 :           IF (ink>n) THEN
     177           0 :              IF ( mpi%irank == 0 ) THEN
     178           0 :                 WRITE (6,*) 'CAUTION!!!  All calculated eigenvalues ', 'are below ef + 8kt.'
     179             :              END IF
     180             :           ENDIF
     181             : 
     182         167 :           w_near_ef = weight - w_below_emin
     183             : 
     184         167 :           IF (w_near_ef.GT.del) THEN
     185             :              !
     186             :              !--->       STATES BETWEEN EF-8kt AND EF+8kt AVAILABLE:
     187             :              !--->            ADJUST FERMI-ENERGY BY NEWTON-METHOD
     188             :              !
     189         122 :              nocst = ink - 1
     190         122 :              CALL ef_newton(n,mpi%irank, inkem,nocst,index,tkb,e, w_near_ef,ef,we)
     191             :              !
     192         122 :              IF ( mpi%irank == 0 ) THEN
     193         122 :                 WRITE (6,FMT=8030) ef,spindg*weight, spindg*w_below_emin,spindg* (w_below_emin+w_near_ef)
     194             :              END IF
     195             : 
     196             :           ELSE
     197             :              !
     198             :              !--->       NO STATES BETWEEN EF-8kt AND EF+8kt AVAILABLE
     199             :              !
     200          45 :              IF ( mpi%irank == 0 ) WRITE (6,FMT=8020)
     201          45 :              nocst = nstef
     202          45 :              we(INDEX(nocst)) = we(INDEX(nocst)) - wfermi
     203          45 :              ef = efermi
     204          45 :              tkb = 0.0
     205             :           END IF
     206             :        ELSE
     207             :           !
     208             :           !--->    NO STATES ABOVE EF AVAILABLE
     209             :           !
     210           0 :           tkb = 0.0
     211           0 :           nocst = nstef
     212           0 :           we(INDEX(nocst)) = we(INDEX(nocst)) - wfermi
     213             :        END IF
     214             : 
     215             :     ELSE
     216             :        !
     217             :        !---> NO TEMPERATURE BROADENING IF tkb < del
     218             :        !
     219           0 :        nocst = nstef
     220           0 :        we(INDEX(nocst)) = we(INDEX(nocst)) - wfermi
     221             :     END IF
     222             :     !
     223             :     !      write(6,*) nocst,'    nocst in ferhis'
     224             :     !      do  ink = 1,nocst
     225             :     !         write(6,*) ink,index(ink),we(index(ink)),
     226             :     !     +      '    ink,index(ink),we(index(ink)): weights for eigenvalues'
     227             :     !      end do
     228             :     !
     229             :     !
     230             :     !=======>   DETERMINE OCCUPATION NUMBER AND WEIGHT OF EIGENVALUES
     231             :     !                     FOR EACH K_POINT
     232             :     !
     233         344 :     w_iks(:,:,:) = 0.0
     234             : 
     235         167 :     IF ( mpi%irank == 0 ) WRITE (6,FMT=8080) nocst
     236       15942 :     DO i=1,nocst
     237       15942 :        w_iks(idxeig(INDEX(i)),idxkpt(INDEX(i)),idxjsp(INDEX(i))) = we(INDEX(i))
     238             :     ENDDO
     239             :     !
     240             :     !======>   CHECK SUM OF VALENCE WEIGHTS
     241             :     !
     242             : 
     243         167 :     wvals = 0.0
     244         344 :     DO  js = 1,nspins
     245        1030 :        DO  k = 1,kpts%nkpt
     246         863 :           wvals = wvals + SUM(w_iks(:ne(k,js),k,js))
     247             :        ENDDO
     248             :     ENDDO
     249             : 
     250         167 :     IF ( mpi%irank == 0 ) WRITE (6,FMT=8070) wvals
     251             :     !
     252             :     !
     253             :     !=======>   DETERMINE ENTROPY
     254             :     !
     255             :     !
     256             :     ! --->   formula for entropy:
     257             :     !
     258             :     !        entropy = - two * sum wtkpt(kpt) * 
     259             :     !                          kpt
     260             :     !                       { sum ( f(e(kpt,nu))*log(f(e(kpt,nu)))
     261             :     !                          n
     262             :     !                              +(1-f(e(kpt,nu)))*log(1-f(e(kpt,nu))) )  }
     263             :     !
     264             :     !        here we have   w(n,kpt,js)= spindg*wghtkp(kpt)*f(e(kpt,n))
     265             :     !
     266         167 :     entropy = 0.0
     267         344 :     DO js = 1,nspins
     268        1030 :        DO kpt = 1 , kpts%nkpt
     269       29359 :           DO nocc=1,ne(kpt,js) 
     270       28496 :              fermikn = w_iks(nocc,kpt,js)/kpts%wtkpt(kpt)
     271       28496 :              IF ( fermikn .GT. 0.0 .AND. fermikn .LT. 1.0 ) &
     272        1473 :                   entropy = entropy + kpts%wtkpt(kpt) * ( fermikn * LOG( fermikn) + ( 1.0 - fermikn) * LOG( 1.0 - fermikn) )
     273             :           END DO
     274             :        END DO
     275             :     ENDDO
     276         167 :     entropy = -spindg*entropy
     277         167 :     results%ts = tkb*entropy
     278         167 :     IF ( mpi%irank == 0 ) WRITE (6,FMT=8060) entropy,entropy*3.0553e-6 !: boltzmann constant in htr/k
     279             : 
     280             : 
     281             : 
     282             :     !
     283             :     !=======>   DETERMINE SINGLE PARTICLE ENERGY
     284             :     !
     285             :     !
     286             : 
     287         167 :     seigv = seigv+spindg*DOT_PRODUCT(e(INDEX(:nocst)),we(INDEX(:nocst)))
     288         167 :     IF (mpi%irank == 0) THEN
     289         501 :        attributes = ''
     290         167 :        WRITE(attributes(1),'(f20.10)') seigv
     291         167 :        WRITE(attributes(2),'(a)') 'Htr'
     292         167 :        CALL writeXMLElement('sumValenceSingleParticleEnergies',(/'value','units'/),attributes)
     293         167 :        WRITE (6,FMT=8040) seigv
     294             :     END IF
     295             : 
     296             : 
     297             :     !
     298             :     ! 7.12.95 r.pentcheva   seigscv = seigsc + seigv   will be
     299             :     ! calculated in fermie
     300             :     !
     301             : 8000 FORMAT (/,10x,'==>efrmhi: not enough wavefunctions.',i10,2e20.10)
     302             : 8010 FORMAT (10x,'charge neutrality (T=0)     :',f11.6,'    (zero if ',&
     303             :          &       'the highest occ. eigenvalue is "entirely" occupied)')
     304             : 8020 FORMAT (/,10x,'no eigenvalues within 8 tkb of ef',&
     305             :          &       ' reverts to the t=0 k method.')
     306             : 8030 FORMAT (/,5x,'-->  new fermi energy            :',f11.6,' htr',&
     307             :          &       /,10x,'valence charge              :',f11.6,' e ',/,10x,&
     308             :          &       'actual charge blw ef-8kt    :',f11.6,' e ',/,10x,&
     309             :          &       'actual charge blw ef+8kt    :',f11.6,' e ')
     310             : 8040 FORMAT (/,10x,'sum of val. single particle energies: ',f20.10,&
     311             :          &       ' htr',/)
     312             : 8050 FORMAT (/,10x,'bandgap                     :',f11.6,' htr')
     313             : 8060 FORMAT (10x,'entropy         :',f11.6,' *kb htr/K =',&
     314             :          &       f10.5,' htr/K')
     315             : 8070 FORMAT (10x,'sum of the valence weights  :',f11.6)
     316             : 8080 FORMAT (10x,'number of occ. states       :',i10)
     317             : 
     318         167 :   END SUBROUTINE ferhis
     319             : END MODULE m_ferhis

Generated by: LCOV version 1.13