C ======================================================================
C
C
      SUBROUTINE EFERMI (NEL,NBANDS,DEL,NSPPTS,NDIM8,NDIM10,
     +                   WEIGHT,OCC,EF,EIGVAL,SORT)
C
C-----------------------------------------------------------------------
C
C     WRITTEN BY RICHARD NEEDS ON 9TH DECEMBER 1983
C     GIVEN THE EIGENVALUES IN EIGVAL AND THE WEIGHTS OF THE
C     K-POINTS IN WEIGHT THIS SUBROUTINE CALCULATES THE FERMI LEVEL
C     EF AND THE OCCUPANCY OF THE STATES OCC.
C
C     METHOD: C-L FU AND K-M HO, PHYS. REV. B 28, 5480 (1983)
C     GAUSSIAN SMEARING OF EIGENVALUES WHEN COMPUTING OCCUPATION
C     NOTE: FOR SUM OF BANDS WE DO NOT SMEAR EIGENVALUES,
C     AS WAS DONE BY FU AND HO. THE DIFFERENCE IS EASILY
C     CALCULATED.
C
C     NEL ..... NUMBER OF ELECTRONS PER UNIT CELL
C     NBANDS .. NUMBER OF BANDS FOR EACH K-POINT
C     DEL ..... WIDTH OF GAUSSIAN SMEARING FUNCTION
C     NSPPTS .. NUMBER OF K-POINTS
C     NDIM8 ... MAXIMUM NUMBER OF BANDS AT A K-POINT
C     WEIGHT .. THE WEIGHT OF EACH K-POINT
C     OCC ..... THE OCCUPANCY OF EACH STATE
C     EF ...... THE FERMI ENERGY
C     SORT .... THE EIGENVALUES ARE WRITTEN INTO SORT WHICH IS
C               THEN SORTED INTO ASCENDING NUMERICAL VALUE, FROM
C               WHICH BOUNDS ON EF CAN EASILY BE OBTAINED
C     EIGVAL .. CONTAINS THE BEST EIGENVALUES AVAILABLE
C---------------------------------------------------------------------
C     12-Mar-90 Obtained from B. Hammer 
C     12-MAR-90 NBANDS changed to be the same for all k-points XW
C
C
      REAL SORT(NDIM8*NDIM10)
      REAL OCC(NDIM8,NDIM10),WEIGHT(NDIM10),EIGVAL(NDIM8,NDIM10)
      EXTERNAL SERFC
C.....WARNINGS
      COMMON /WARN/ IWARN
      PARAMETER ( NCYCLE = 20, NDIV = 9 )
C--------------------------------------------------------------------
C
      DSOR = 1.0 /FLOAT(NDIV)
      Z    = FLOAT(NEL)
C
C     COPY EIGVAL INTO SORT
      NEIG = 0
      DO 10 ISPPT = 1,NSPPTS
        DO 20  J = 1, NBANDS
          NEIG = NEIG + 1
20        SORT(NEIG) = EIGVAL(J,ISPPT)
10      CONTINUE
C
C  SORT THE ARRAY INTO ASCENDING ORDER OF EIGENVALUE
C
C  SORT THE RECORDS WITH 'STRAIGHT INSERTION'-METHOD (NOT A TOO BAD
C  CHOICE OF SORT-ALGORITHM AS THE ENERGIES ARE ALLMOST ALREADY SORTED)
C=======================================================================
      DO 26 N=2,NSPPTS*NBANDS
        EN=SORT(N)
        DO 22 NN=N-1,1,-1
          IF (SORT(NN).LE.EN) GO TO 24
          SORT(NN+1)=SORT(NN)
  22    CONTINUE
        NN=0
  24    SORT(NN+1)=EN
  26  CONTINUE
C
C     FIND AN UPPER BOUND E2 AND A LOWER BOUND E1 ON THE
C     FERMI ENERGY
      INT = NEL*NSPPTS
      IF ( MOD(INT,2) .EQ. 1) THEN
        I1 = INT/2
        I2 = INT/2 + 2
      ELSE
        I1 = INT/2 - 1
        I2 = INT/2 + 1
        ENDIF
      I1 = MAX(I1, 1)
      I2 = MIN(I2, NEIG)
C
30    E1 =  SORT(I1)
C
      Z1 = 0.0 
      DO 40 ISPPT = 1,NSPPTS
        DO 50 J = 1,NBANDS
          X = (E1 - EIGVAL(J,ISPPT))/DEL
50        Z1 = Z1 + WEIGHT(ISPPT)*( 2.0 - SERFC(X) )
40      CONTINUE
C
      IF (Z1 .GT. Z-1.0E-6) THEN
        I1 = I1 - 1
        IF (I1 .GT. 0) THEN
          GOTO 30
        ELSE
          I1 = 1
          E1 = SORT(I1) - DEL
          WRITE (*,*) 'EFERMI *** LOWER BOUND ON EF SET TO',E1
          GOTO 70
          ENDIF
        ENDIF
C
70      E2 = SORT(I2)
C
      Z2 = 0.0 
      DO 60 ISPPT = 1, NSPPTS
        DO 80 J = 1,NBANDS
          X = (E2 - EIGVAL(J,ISPPT))/DEL
80        Z2 = Z2 + WEIGHT(ISPPT)*( 2.0 - SERFC(X) )
60      CONTINUE
C
      IF (Z2 .LT. Z + 1.E-6) THEN
        I2 = I2 + 1
        IF (I2 .LE. NEIG) THEN
          GOTO 70
        ELSE
          I2 = NEIG
          E2 = SORT(I2) + DEL
          WRITE (*,1100) E2
1100      FORMAT(' EFERMI *** WARNING ***'/
     +      ' FERMI LEVEL UPPER BOUND IS > LARGEST EIGENVALUE, E2=',
     +      F12.4,' eV')
          IWARN = 1
          ENDIF
        ENDIF
C
C     FIND FERMI ENERGY ENERGY BETWEEN BOUNDS E1 AND E2
      DO 90 ILOOP = 1,NCYCLE
        DIV = (E2-E1)*DSOR
        DO 100 I = 1,NDIV
          Z2 = 0.0
          EUP = E1 + FLOAT(I)*DIV
          DO 120 ISPPT = 1,NSPPTS
            DO 110 J = 1,NBANDS
              X = (EUP - EIGVAL(J,ISPPT))/DEL
110           Z2 = Z2 + WEIGHT(ISPPT)*( 2.0 - SERFC(X) )
120         CONTINUE
          IF (Z1 .GT. Z) GOTO 1000
          IF (Z2 .GT. Z) GOTO 130
C         THIS TEST IS NEEDED TO CATCH SEMICONDUCTORS
C          IF (ABS(Z2-Z) .LT. 1.0E-10) GOTO 150
          IF (ABS(Z2-Z) .LT. 1.0E-5) GOTO 150
          Z1 = Z2
100       CONTINUE
        GOTO 1010
C130     IF ( Z2-Z1 .LT. 1.0E-10) GOTO 160
130     IF ( Z2-Z1 .LT. 1.0E-5) GOTO 160
        E1 = EUP - DIV
90      E2 = EUP
      WRITE(*,*) '***WARNING*** FERMI ENERGY MAY NOT BE ACCURATE'
      WRITE(*,140) NCYCLE
140   FORMAT(' AFTER',I6,' CYCLES, REQUIRED CONVERGENCE NOT OBTAINED')
      GOTO 160
C
C     WRITE OUT FERMI ENERGY
150   EF = EUP
      GOTO 170
160   EF = EUP - 0.5*DIV
170   WRITE (*,180) EF
180   FORMAT (' FERMI ENERGY = ',F12.8,' EV')
C
C     FORM OCCUPATIONS OCC(NBDS,NSPPTS)
      DO 190 ISPPT = 1, NSPPTS
        DO 200 J = 1,NBANDS
          X = ( EF-EIGVAL(J,ISPPT))/DEL
200       OCC(J,ISPPT) = 2.0 - SERFC(X)
190     CONTINUE
C
C     TEST WHETHER OCCUPANCY ADDS UP TO Z
      TEST = 0.0
      DO 210 ISPPT = 1,NSPPTS
        DO 215 J = 1,NBANDS
215       TEST = TEST + WEIGHT(ISPPT)*OCC(J,ISPPT)
210     CONTINUE
      IF ( ABS(TEST-Z) .GT. 1.0E-9) THEN
        WRITE(*,*) '*** WARNING ***'
        WRITE(*,220) TEST,NEL
220     FORMAT(' SUM OF OCCUPANCIES =',F15.12 ,' BUT NEL =',I5)
      ELSE
        WRITE(*,230) TEST
230     FORMAT(' TOTAL CHARGE = ',F12.8)
        ENDIF
C
C     TEST WHETHER THE MATERIAL IS A SEMICONDUCTOR
      IF ( MOD( NEL, 2) .EQ. 1) RETURN
      INEL = NEL/2
      ELOW = EIGVAL(INEL+1,1)
      DO 310 ISPPT = 2,NSPPTS
        ELOW = AMIN1( ELOW, EIGVAL(INEL+1,ISPPT))
310     CONTINUE
      DO 320 ISPPT = 1,NSPPTS
        IF (ELOW .LT. EIGVAL(INEL,ISPPT)) RETURN
320     CONTINUE
      WRITE (*,*) 'MATERIAL MAY BE A SEMICONDUCTOR'
C
      RETURN
C----------------------------------------------------------------------
C      ERROR MESSAGES
C----------------------------------------------------------------------
1000  WRITE(*,*) '***** ERROR *****'
      WRITE(*,*) 'FERMI ENERGY LESS THAN LOWER SEARCH BOUND SET'
      WRITE(*,*) EIGVAL
      CALL EXIT
      RETURN
C
1010  WRITE(*,*) '***** ERROR *****'
      WRITE(*,*) 'FERMI ENERGY GREATER THAN UPPER SEARCH BOUND SET'
      WRITE(*,*) 'EigVal: ', EIGVAL
      write(*,*) 'Weight: ', WEIGHT
      WRITE(*,*) 'Z1, Z2, Z:', Z1,Z2,Z
      WRITE(*,*) 'E1,EUP:',  E1,EUP 
      CALL EXIT
      RETURN
C
      END
C
C =========================================================================
C
C     DOUBLE PRECISION FUNCTION ERFC(XX)
      FUNCTION SERFC(XX1)
      IMPLICIT DOUBLE PRECISION (A-H,O-W,Y-Z)
C
C     COMPLEMENTARY ERROR FUNCTION
C     FROM THE SANDIA MATHEMATICAL PROGRAM LIBRARY
C
C     XMAX IS THE VALUE BEYOND WHICH ERFC(X) = 0 .
C     IT IS COMPUTED AS SQRT(LOG(RMIN)), WHERE RMIN IS THE
C     SMALLEST REAL NUMBER REPRESENTABLE ON THE MACHINE.
C     IBM VALUE: (THE INTRINSIC ERFC COULD ALSO BE USED)
C     PARAMETER ( XMAX = 13.4 )
C     VAX VALUE: (XMAX = 9.3)
C -----------------------------------  
C     12-Mar-90  Obtained from B. Hammer 
C     12-MAR-90  Changed to single precision at the end XW
C                also XX1
      REAL XMAX
      PARAMETER ( XMAX = 9.3D0)
C
      DIMENSION P1(4),Q1(4),P2(6),Q2(6),P3(4),Q3(4)
C     REAL*4 XX1, SERFC
C
      DATA P1 /242.66 79552 30531 8D0 , 21.979 26161 82941 5D0 ,
     + 6.9963 83488 61913 6D0 , -3.5609 84370 18153 9D-2/
      DATA Q1 /215.05 88758 69861 2D0 , 91.164 90540 45149 0D0,
     + 15.082 79763 04077 9D0 , 1.0D0/
      DATA P2 /22.898 99285 1659D0 , 26.094 74695 6075D0 ,
     + 14.571 89859 6926D0 , 4.2677 20107 0898D0 ,
     + 0.56437 16068 6381D0 , -6.0858 15195 9688 D-6/
      DATA Q2 /22.898 98574 9891D0 , 51.933 57068 7552D0 ,
     + 50.273 20286 3803D0 , 26.288 79575 8761D0 ,
     + 7.5688 48229 3618D0 , 1.0D0/
      DATA P3 /-1.2130 82763 89978 D-2 , -0.11990 39552 68146 0D0 ,
     + -0.24391 10294 88626D0 , -3.2431 95192 77746 D-2/
      DATA Q3 /4.3002 66434 52770 D-2 , 0.48955 24419 61437D0 ,
     + 1.4377 12279 37118D0 , 1.0D0/
C     1/SQRT(PI)
      DATA SQPI /0.56418 95835 47756D0/
C
      XX=XX1
C----------------------------------------------------------------------
      IF (XX .GT.  XMAX)    GOTO 330
      IF (XX .LT. -XMAX)    GOTO 320
      X = ABS(XX)
      X2 = X*X
      IF (X .GT. 4.0D0)     GOTO 300
      IF (X .GT. 0.46875D0) GOTO 200
C
C     -46875 < X < 0.46875
      ERFC = X*(P1(1) + X2*(P1(2) + X2*(P1(3) + X2*P1(4))))
      ERFC = ERFC/(Q1(1) + X2*(Q1(2) + X2*(Q1(3) + X2*Q1(4))))
      IF (XX .LT. 0.0) ERFC = - ERFC
      ERFC = 1.0D0 - ERFC
      GOTO 9999
C
200   ERFC = EXP( -X2)*(P2(1) + X*(P2(2) + X*(P2(3) + X*(P2(4) +
     + X*(P2(5) + X*P2(6))))))
      ERFC = ERFC/(Q2(1) + X*(Q2(2) + X*(Q2(3) + X*(Q2(4) + X*(Q2(5) +
     + X*Q2(6))))))
      IF (XX .LE. 0.0) ERFC = 2.0D0 - ERFC
      GOTO 9999
C
300   XI2 = 1.0D0/X2
      ERFC = XI2*(P3(1) + XI2*(P3(2) + XI2*(P3(3) + XI2*P3(4))))/
     + (Q3(1) + XI2*(Q3(2) + XI2*(Q3(3) + XI2*Q3(4))))
      ERFC = EXP( -X2)*(SQPI + ERFC)/X
      IF (XX .LT. 0.0) ERFC = 2.0D0 - ERFC
      GOTO 9999
C
320   ERFC = 2.0D0
      GOTO 9999
330   ERFC = 0.0D0
C
9999  SERFC=ERFC
      RETURN
      END
