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

Generated by: LCOV version 1.13