C
C $Id: doferm.F,v 1.3 1998/07/16 16:39:39 jjv5 Exp $
C
C------------------------------------------------------------------------

C     new code: NOT PARALLELIZED!!!!!
C     more than one fermi energies when setch == .true.

      SUBROUTINE DOFERM(IERROR)
C
C     USES AN ITERATIVE BISECTION TECHNIQUE TO DETERMINE THE FERMI
C     ENERGY EFERMI.  ERROR FLAG IS SET IF MORE THAN 100 ITERATIONS
C     ARE NECESSARY.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "divcon.dim"
#include "divcon.h"

C
C     BOLTZMANN CONSTANT IN eV/Kelvin:
C
      DATA BOLTZ /8.617335408D-5/
      SAVE BOLTZ
C
      BETAH = 1.0D0/(BOLTZ*TEMPK)
      ETOLER = 1.0D-8
      MAXIT = 100

      if (setch) then

         do 10 icr=1,ncores

            i0 = isubend(icr-1)
            i1 = i0+1
            i2 = isubend(icr)
            n = i2-i0
            iorbmn = iorbpt(i1)
            iorbmx = iorbpt(i2+1)-1
C
C     .     DETERMINE THE OBSERVED RANGE OF EIGENVALUES FOR ALL SUBSYSTEMS.
C
            imax = iorbpt(i1+1)-1
            emin = eval(iorbmn)
            emax = eval(imax)
            IF(N.GT.1)THEN
               DO K=i1,i2
                  IMIN = IORBPT(K)
                  IMAX = IORBPT(K+1)-1
                  EMIN = MIN(EMIN,EVAL(IMIN))
                  EMAX = MAX(EMAX,EVAL(IMAX))
               enddo
            ENDIF
C
C     .     USE A BISECTION TECHNIQUE TO DETERMINE THE FERMI ENERGY.
C
            ELECS = NELECEF(icr)
            NITER = 0
 2          EFERMI(icr) = (EMAX + EMIN)*0.5
C
C     .     SET FERMI OCCUPATION NUMBERS AND USE THE SUM OF WEIGHTED SQUARED
C     .     EIGENVECTOR COEFFICIENTS EVECSQ TO GET THE NUMBER OF ELECTRONS.
C
            NITER = NITER + 1
            ETOTAL = 0.0D0
            DO I=iorbmn,IORBMX
               DELTAE = EVAL(I) - EFERMI(icr)
               ARG = BETAH*DELTAE
               IF(ARG.GT.50.0D0)THEN
                  FERMI(I) = 0.0D0
               ELSEIF(ARG.LT.-50.0D0)THEN
                  FERMI(I) = 1.0D0
               ELSE
                  FERMI(I) = 1.0D0/(1.0D0 + EXP(ARG))
               ENDIF
               ETOTAL = ETOTAL + FERMI(I)*EVECSQ(I)
            enddo
C
C     .     WE'RE DONE IF THE COMPUTED NUMBER OF ELECTRONS IS CLOSE ENOUGH
C     .     TO THE TARGET.
C
            DIFF = ABS(ETOTAL - ELECS)
            IF(DIFF.LT.ETOLER) then
               goto 10
            endif
C
C     .     THROW AWAY LOWER OR UPPER HALF OF INTERVAL, DEPENDING UPON WHETHER
C     .     THE COMPUTED NUMBER OF ELECTRONS IS TOO LOW OR TOO HIGH.  NOTE THAT
C     .     RAISING THE FERMI ENERGY INCREASES THE NUMBER OF ELECTRONS BECAUSE
C     .     MORE MO'S BECOME FULLY OCCUPIED.
C
            IF(ETOTAL.LT.ELECS)THEN
               EMIN = EFERMI(icr)
            ELSE
               EMAX = EFERMI(icr)
            ENDIF

            IF(NITER.EQ.MAXIT)THEN
               IERROR = 1
               RETURN
            ENDIF
            GO TO 2

 10      continue

      else
         IORBMX = IORBPT(NSUB+1)-1
C
C     DETERMINE THE OBSERVED RANGE OF EIGENVALUES FOR ALL SUBSYSTEMS.
C
         EMIN = EVAL(1)
         EMAX = EVAL(IORBPT(2)-1)
         IF(NSUB.GT.1)THEN
            DO 20 K=1,NSUB
               IMIN = IORBPT(K)
               IMAX = IORBPT(K+1)-1
               EMIN = MIN(EMIN,EVAL(IMIN))
               EMAX = MAX(EMAX,EVAL(IMAX))
 20         CONTINUE
         ENDIF

C
C     USE A BISECTION TECHNIQUE TO DETERMINE THE FERMI ENERGY.
C
         ELECS = NELEC
         NITER = 0
 100     EFERMI(1) = (EMAX + EMIN)/2.0D0
C
C     SET FERMI OCCUPATION NUMBERS AND USE THE SUM OF WEIGHTED SQUARED
C     EIGENVECTOR COEFFICIENTS EVECSQ TO GET THE NUMBER OF ELECTRONS.
C
         NITER = NITER + 1
         ETOTAL = 0.0D0
         DO 200 I=1,IORBMX
            DELTAE = EVAL(I) - EFERMI(1)
            ARG = BETAH*DELTAE
            IF(ARG.GT.50.0D0)THEN
               FERMI(I) = 0.0D0
            ELSEIF(ARG.LT.-50.0D0)THEN
               FERMI(I) = 1.0D0
            ELSE
               FERMI(I) = 1.0D0/(1.0D0 + EXP(ARG))
            ENDIF
            ETOTAL = ETOTAL + FERMI(I)*EVECSQ(I)
 200     CONTINUE
C
C     WE'RE DONE IF THE COMPUTED NUMBER OF ELECTRONS IS CLOSE ENOUGH
C     TO THE TARGET.
C
         DIFF = ABS(ETOTAL - ELECS)

         IF(DIFF.LT.ETOLER) RETURN
C
C     THROW AWAY LOWER OR UPPER HALF OF INTERVAL, DEPENDING UPON WHETHER
C     THE COMPUTED NUMBER OF ELECTRONS IS TOO LOW OR TOO HIGH.  NOTE THAT
C     RAISING THE FERMI ENERGY INCREASES THE NUMBER OF ELECTRONS BECAUSE
C     MORE MO'S BECOME FULLY OCCUPIED.
C

         IF(ETOTAL.LT.ELECS)THEN
            EMIN = EFERMI(1)
         ELSE
            EMAX = EFERMI(1)
         ENDIF
         IF(NITER.EQ.MAXIT)THEN
            IERROR = 1
            RETURN
         ENDIF
         GO TO 100

      endif
      END