fergwt.f90 4.78 KB
Newer Older
1 2 3 4 5 6 7 8
MODULE  m_fergwt
  USE m_juDFT
  !****************************************************************
  !     determines the fermi energy and weights for the k-space
  !     integration using gaussing-smearing method.
  !                                               c.l.fu
  !*****************************************************************
CONTAINS
9
  SUBROUTINE fergwt(kpts,input,mpi, ne,eig, ef,w_iks,seigv)
10

11
    USE m_constants
12 13 14 15 16 17
    USE m_types
    IMPLICIT NONE

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

  END SUBROUTINE fergwt
END MODULE m_fergwt