LCOV - code coverage report
Current view: top level - fermi - fergwt.f90 (source / functions) Hit Total Coverage
Test: FLEUR test coverage Lines: 0 74 0.0 %
Date: 2024-03-28 04:22:06 Functions: 0 1 0.0 %

          Line data    Source code
       1             : MODULE  m_fergwt
       2             :   USE m_juDFT
       3             :   !****************************************************************
       4             :   !     determines the fermi energy and weights for the k-space
       5             :   !     integration using gaussing-smearing method.
       6             :   !                                               c.l.fu
       7             :   !*****************************************************************
       8             : CONTAINS
       9           0 :   SUBROUTINE fergwt(kpts,input,fmpi, ne,eig, ef,w_iks,seigv,l_output)
      10             : 
      11             :     USE m_constants
      12             :     USE m_types
      13             :     IMPLICIT NONE
      14             : 
      15             :     TYPE(t_mpi),INTENT(IN)       :: fmpi
      16             :     TYPE(t_input),INTENT(IN)     :: input
      17             :     TYPE(t_kpts),INTENT(IN)      :: kpts
      18             :     REAL,        INTENT(INOUT)   :: ef,seigv
      19             :     REAL,INTENT(INOUT)           :: w_iks(:,:,:)
      20             :     !     ..
      21             :     !     ..
      22             :     !     .. Array Arguments ..
      23             :     INTEGER, INTENT (IN) :: ne(:,:)    !(kpts%nkpt,dimension%jspd)
      24             :     REAL,    INTENT (IN) :: eig(:,:,:) !input%neig,kpts%nkpt,dimension%jspd)
      25             :     !     .. Logical Arguments ..
      26             :     LOGICAL,INTENT(IN) :: l_output
      27             :     !     ..
      28             :     !     .. Local Scalars ..
      29             :     REAL chmom,de,ef0,ef1,elow,en,eps,eup,fac,fact1,s,s0,s1,s2,&
      30             :          workf,wt,wtk,zcdiff,zero,seigv1
      31             :     INTEGER i,ifl,it,jspin,k,nbnd
      32             :     !     ..
      33             :     !     .. External Functions ..
      34             :     REAL  erf
      35             :     !     ..
      36             :     !     .. Data statements ..
      37             :     DATA zero/0.e0/,eps/1.e-5/,eup/3.0e0/,elow/-3.0e0/
      38             :     !     ..
      39           0 :     fact1 = input%tkb/SQRT(pi_const)
      40             :     !     ---> determines ef
      41           0 :     ifl = 0
      42             :     conv_loop:DO WHILE (.TRUE.)
      43           0 :        DO  it = 1,50
      44           0 :           s = 0.
      45           0 :           DO  jspin = 1,input%jspins
      46           0 :              DO  k = 1,kpts%nkpt
      47           0 :                 wtk = kpts%wtkpt(k)
      48           0 :                 nbnd = ne(k,jspin)
      49           0 :                 DO  i = 1,nbnd
      50           0 :                    en = eig(i,k,jspin)
      51           0 :                    de = (en-ef)/input%tkb
      52           0 :                    wt = 2.0
      53           0 :                    IF (de.GT.eup) wt = 0.0
      54           0 :                    IF (de.GE.elow .AND. de.LE.eup) THEN
      55           0 :                       IF (de.LT.zero) THEN
      56           0 :                          wt = 1. + ERF(-de)
      57             :                       ELSE
      58           0 :                          wt = 1. - ERF(de)
      59             :                       END IF
      60             :                    END IF
      61           0 :                    s = s + wt*wtk
      62           0 :                    w_iks(i,k,jspin) = wt/2.
      63             :                 ENDDO
      64             :              ENDDO
      65             :           ENDDO
      66           0 :           s = s/REAL(input%jspins)
      67           0 :           zcdiff = input%zelec - s
      68           0 :           IF (ABS(zcdiff).LT.eps) EXIT conv_loop
      69           0 :           IF (ifl.EQ.0) THEN
      70           0 :              ifl = 1
      71           0 :              ef0 = ef
      72           0 :              ef = ef + 0.003
      73           0 :              s0 = s
      74             :           ELSE
      75           0 :              fac = (s0-s)/ (input%zelec-s)
      76           0 :              IF (ABS(fac).LT.1.0e-1) THEN
      77           0 :                 ef0 = ef
      78           0 :                 s0 = s
      79           0 :                 IF (zcdiff.GE.zero) THEN
      80           0 :                    ef = ef + 0.003
      81             :                 ELSE
      82           0 :                    ef = ef - 0.003
      83             :                 END IF
      84             :              ELSE
      85           0 :                 ef1 = ef
      86           0 :                 ef = ef + (ef0-ef)/fac
      87           0 :                 ef0 = ef1
      88           0 :                 s0 = s
      89             :              END IF
      90             :           END IF
      91             :        ENDDO
      92           0 :        eps = 1.25*eps
      93           0 :        IF ( fmpi%irank == 0 .and. l_output) WRITE (oUnit,FMT=8000) eps
      94             : 8000   FORMAT (10x,'warning: eps has been increased to',e12.5)
      95             :     ENDDO conv_loop
      96           0 :     workf = -hartree_to_ev_const*ef
      97           0 :     IF ( fmpi%irank == 0 .and. l_output) THEN
      98           0 :        WRITE (oUnit,FMT=8010) ef,workf,s
      99             :     END IF
     100             : 8010 FORMAT (/,10x,'fermi energy=',f10.5,' har',3x,'work function=',&
     101             :                f10.5,' ev',/,10x,'number of valence electrons=',f10.5)
     102           0 :     IF (ABS(zcdiff).GT.5.0e-4) THEN
     103             :        CALL juDFT_error('Fermi-level determination did not converge',&
     104           0 :                         hint ="change temperature or set input = F", calledby ="fergwt")
     105             :     ENDIF
     106           0 :     DO  jspin = 1,input%jspins
     107           0 :        IF ( fmpi%irank == 0 .and. l_output) WRITE (oUnit,FMT=8020) jspin
     108             : 8020   FORMAT (/,/,5x,'band-weighting factor for spin=',i5)
     109           0 :        DO  k = 1,kpts%nkpt
     110           0 :           nbnd = ne(k,jspin)
     111           0 :           IF ( fmpi%irank == 0 .and. l_output) WRITE (oUnit,FMT=8030) k
     112             : 8030      FORMAT (/,5x,'k-point=',i5,/)
     113           0 :           w_iks(:,k,jspin) = kpts%wtkpt(k)*w_iks(:,k,jspin)
     114           0 :           IF ( fmpi%irank == 0 .and. l_output) THEN
     115           0 :              WRITE (oUnit,FMT=8040) (w_iks(i,k,jspin),i=1,nbnd)
     116             :           END IF
     117             : 8040      FORMAT (5x,16f6.3)
     118             :        ENDDO
     119             :     ENDDO
     120           0 :     s1 = 0.
     121           0 :     s2 = 0.
     122           0 :     DO  jspin = 1,input%jspins
     123           0 :        s = 0.
     124           0 :        DO  k = 1,kpts%nkpt
     125           0 :           DO  i = 1,ne(k,jspin)
     126           0 :              s = s + w_iks(i,k,jspin)
     127           0 :              seigv = seigv + w_iks(i,k,jspin)*eig(i,k,jspin)
     128           0 :              en = eig(i,k,jspin)
     129           0 :              de = (en-ef)/input%tkb
     130             :              !     ---> correction term
     131           0 :              IF (ABS(de).LT.3.) THEN
     132           0 :                 de = de*de
     133           0 :                 s2 = s2 + EXP(-de)*kpts%wtkpt(k)
     134             :              END IF
     135             :           ENDDO
     136             :        ENDDO
     137           0 :        s1 = s1 + s
     138             :     ENDDO
     139           0 :     seigv = (2/input%jspins)*seigv
     140           0 :     seigv1 = (1/input%jspins)*fact1*s2
     141           0 :     chmom = s1 - input%jspins*s
     142           0 :     IF ( fmpi%irank == 0 .and. l_output) THEN
     143           0 :        WRITE (oUnit,FMT=8050) seigv - seigv1,s1,chmom
     144             :     END IF
     145             : 8050 FORMAT (/,10x,'sum of eigenvalues-correction=',f12.5,/,10x,&
     146             :           'sum of weight                =',f12.5,/,10x,&
     147             :           'total moment                 =',f12.5,/)
     148             : 
     149           0 :   END SUBROUTINE fergwt
     150             : END MODULE m_fergwt

Generated by: LCOV version 1.14