C
C $Id: doscf.F,v 1.9 1998/07/16 16:39:41 jjv5 Exp $
C
C------------------------------------------------------------------------
      SUBROUTINE DOSCF(EELECT1,ECORE1,ETOT1,IERROR)
C
C
C     DRIVER ROUTINE TO CARRY OUT THE ITERATIVE SELF-CONSISTENT FIELD
C     CALCULATION.  ENERGIES ARE RETURNED IN eV.
C
C     EELECT1 = ELECTRONIC ENERGY
C     ECORE1  = CORE-CORE REPULSIONS
C     ETOT1   = EELECT1 + ECORE1
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "divcon.dim"
#include "divcon.h"
C
      LOGICAL FIRST,PSEUDO,DOSHFT,GOODF,pvecnow,lastit
      DATA FIRST /.TRUE./
      SAVE FIRST,EEOLD,USHIFT,DESCF,DPSCF,ETOLD

      logical frst_call

      pvecnow = .false.
      lastit = .false.
      if (wrtscr) then
         tcalls(1) = 0.0d0
         tcalls(2) = 0.0d0
         ncalls(1) = 0
         ncalls(2) = 0
         icall = 0
      endif

      iIDump = index(keywrd,'PDUMP=')
      if (iIDump.ne.0) then
         call rdnum (keywrd, iIDump, Dump, ierror)
         iNDump = Dump
      endif


      frst_call = first
      tfock = 0.0
C
      IF(FIRST)THEN
        EELECT1 = 99999.0D0
        etot1 = 999999.0D0
        USHIFT = 0.0D0
        ISHIFT = INDEX(KEYWRD,'SHIFT=')
        IF(ISHIFT.NE.0) CALL RDNUM(KEYWRD,ISHIFT,USHIFT,IERROR)
        DESCF = -1.0D0
        DPSCF = -1.0D0
C
C       USER-SPECIFIED CONVERGENCE CRITERIA:
C
        IDESCF = INDEX(KEYWRD,'DESCF=')
        IDPSCF = INDEX(KEYWRD,'DPSCF=')
        IF(IDESCF.NE.0) CALL RDNUM(KEYWRD,IDESCF,DESCF,IERROR)
        IF(IDPSCF.NE.0) CALL RDNUM(KEYWRD,IDPSCF,DPSCF,IERROR)
      ENDIF
      ASHIFT = USHIFT
      DOSHFT = USHIFT.GT.1.0D-6
C
      IERROR = 0
      IIMAX = IIMAT(NATOMS+1)-1
      IPMAX = IP1(NATOMS+1)
      IJMAX = IJMAT(IPMAX)-1
C
C     ENERGY (eV) AND DENSITY MATRIX CONVERGENCE CRITERIA:
C
      IF(DESCF.GT.0.0D0) EECRT = DESCF
      IF(DPSCF.GT.0.0D0) DENCRT = DPSCF
C
C     OVERRIDE ABOVE CRITERIA IF DENSITY MATRIX HAS VERY GOOD
C     CONVERGENCE.
C
      EEAUX = 10.0D0*EECRT
      DENAUX = DENCRT/10.0D0
      IF(DESCF.GT.0.0D0) EEAUX = DESCF
      IF(DPSCF.GT.0.0D0) DENAUX = DPSCF
C
C     RELAX CONVERGENCE CRITERIA IF THE GRADIENT NORM IS LARGE AND THE
C     STRUCTURE IS NOT TINY.
C
!      IF(GNORM.GT.5.0D0.AND.IIMAX.GT.30)THEN
!        EECRT = EECRT*GNORM*0.2D0
!        DENCRT = 0.25D0*DSQRT(EECRT)
!        EECRT = MIN(EECRT,0.0001D0)
!        DENCRT = MIN(DENCRT,0.002D0)
!        EEAUX = EECRT
!        DENAUX = DENCRT
!      ENDIF
C
      DEPREV = -1.0D0
      DPPREV = 1.0D0
      DPMAX = 1.0D0
C
      PSEUDO = .FALSE.
C
C     DO SCF ITERATIONS.
C
      MAXIT = maxitscf

      NOSCIL = 0
      ITER = 0

      call etimer(tstart)

      ndiag = 0
      DO 500 K=1,MAXIT

c       if(k.gt.2) call frmchk
C
C       IDIAG WILL BE A FLAG TO KEEP TRACK OF WHEN DENSITY MATRIX HAS
C       BEEN INITIALIZED TO DIAGONAL FORM.  IF IDIAG.EQ.1 THEN MOSUB
C       WILL BE TOLD THAT THERE IS NO GOOD ESTIMATE OF THE FERMI ENERGY.
C       ALSO, IF FSHIFT IS CALLED, THEN NO SHIFT WILL BE APPLIED BECAUSE
C       P DOESN'T EVEN COME CLOSE TO COMMUTING WITH F.
C
 10     IF(FIRST.AND.K.EQ.1)THEN
          IDIAG = 1
        ELSE
          IDIAG = 0
        ENDIF

        if (wrtscr) then
           if (setch) then
              do i=1,ncores
C-RDC                  write(iscr,'(" efermi(",i4,")=",F20.10)') i, efermi(i)
              enddo
           else
C-RDC               write(iscr,'(" efermi=",F20.10)') efermi(1)
           endif
        endif

C
C       RESET DENSITY MATRIX AND INVOKE SHIFT (IF IT'S NOT ALREADY IN
C       EFFECT) IF WE'RE UNDERGOING LARGE-SCALE OSCILLATIONS AFTER 10
C       ITERATIONS.
C
        if ((.not.fdmx).and.(.not.pvecnow)) then
C     .    don't mess with the density matrix when a frozen
C     .    calculation on a guess density matrix is performed
           IF(ITER.GT.10)THEN
C
C     .       VG 11/9/98
C     .       Original DIVCON statement:
C     .       IF(EELECT1.GT.(EEOLD+1.0D0).AND.ABS(DPMAX).GT.0.1D0)THEN
C
C     .       In SCRF cycles ECORE1 changes EACH cycle because it contains surface
C     .       charge contributions which are recalculated at each iteration, thus
C     .       the SCF convergence criterion must be done on total energy(ETOT1)
C     .       and NOT on the electronic energy (EELECT1), because it leads to
C     .       excesive solute polarization.
C     .       New statement:
              IF(ETOT1.GT.(ETOLD+1.0D0).AND.ABS(DPMAX).GT.0.1D0)THEN


                 if (screen) then
                    write(*,'(" resetting density matrix, ",
     &                   "invoking shift")')
                 endif

                 ITER = 0
                 INIT = 0
                 CALL INITP(INIT, ierror)
                 if (ierror.ne.0) return
                 IDIAG = 1
                 EELECT1 = 99999.0D0
                 etot1 = 999999.0D0
                 DOSHFT = .TRUE.
                 ASHIFT = MAX(2.0*ASHIFT,USHIFT,10.0D0)
                 DEPREV = 1.0D0
                 DPPREV = 1.0D0
                 NOSCIL = 0
              ENDIF
           ENDIF
        endif 
             
        if ((k.eq.maxit).and.pvec) then
           pvecnow = .true.
           call setprtvec
        endif

        ITER = ITER + 1
C
C       COMPUTE GLOBAL FOCK MATRIX, AND, IF IT'S THE FIRST SCF ITERATION,
C       COMPUTE CORE-CORE REPULSTIONS.
C

        call etimer(tf1)
        if (pmeqm) then
         CALL FOCKPME(K,ECORE1)
        else
         CALL FOCK(K,ECORE1)
        endif
       
        call etimer(tf2)
        tfock = tfock + tf2-tf1

        if (wrtscr) then
C-RDC          write(iscr,'(" time for FOCK =",F11.3)') tf2-tf1
      endif

C
C       GET ELECTRONIC ENERGY USING THE CURRENT FOCK MATRIX AND ORBITALS.
C
        EEOLD = EELECT1
        etold = etot1
        CALL ESCF(EELECT1)
        ETOT1 = EELECT1 + ECORE1
      
   
         DEE = EELECT1 - EEOLD
c        dee = etot1 - etold
! 05/2008  convergence criteria relative to the total energy 
!       dee = (etot1 - etold) / etot1
! 06/2011 search for machine precision
        call deeprec(etot1,eecrt,iout)
        eeaux = 10.0d0*eecrt

!       if (wrtscr) then
!           write(0,'("k=",i4,"   iter=",i4,/" eelect=",f25.8,
!    .          "   etot=",f25.8,
!    .          /" dee=",f20.8,5x,"dpmax=",f15.8)')
!    .          k,iter,eelect1,etot1, DEE,DPMAX
!       endif

        if (last_scrf) return

C     . we're done when the eigenvectors have been printed
        if (lastit) goto 1000

C     . we're done if only one scf cycle is requested
        if (onescf) goto 1000

C
C       CHECK FOR CONVERGENCE.
C
        IF(ABS(DEE).LT.EECRT.AND.ABS(DPMAX).LT.DENCRT)THEN
           ECRT2 = EECRT*2.0D0
           DCRT2 = DENCRT*2.0D0
           IF(ABS(DEPREV).LT.ECRT2.AND.ABS(DPPREV).LT.DCRT2)THEN

!             if (wrtscr) then
!                write(iscr,'(" converged, dee=",f20.8,
!    $                " dpmax=",F20.8)') dee,dpmax
!             endif
              if (pvec) then
                 if (.not.pvecnow) then
!                  if(wrtscr)
!    $                   write(iscr,'(
!    $                   /"will perform one extra iteration ",
!    $                   "to write out eigenvectors"/)')
                    pvecnow = .true.
                    call setprtvec
                 endif
              else
                 goto 1000
              endif
           ENDIF
        ELSEIF(ABS(DEE).LT.EEAUX.AND.ABS(DPMAX).LT.DENAUX)THEN
           EAUX2 = EEAUX*2.0D0
           DAUX2 = DENAUX*2.0D0
           IF(ABS(DEPREV).LT.EAUX2.AND.ABS(DPPREV).LT.DAUX2)THEN

!             if (wrtscr) then
!                write(iscr,'(" converged, dee=",f20.8,
!    $                " dpmax=",F20.8)') dee,dpmax
!             endif
              if (pvec) then
                 if (.not.pvecnow) then
!                  if (wrtscr)
!    $                   write(iscr,'(
!    $                   /"will perform one extra iteration ",
!    $                   "to write out eigenvectors"/)')
                    pvecnow = .true.
                    call setprtvec
                 endif
              else
                 goto 1000
              endif
          ENDIF
       ENDIF
       DEPREV = DEE
       DPPREV = DPMAX
C
C       SEE IF A PSEUDO-DIAGONALIZATION CAN BE DONE INSTEAD OF
C       A FULL DIAGONALIZATION.
C
         IF(STAND.AND..NOT.FULLSCF.AND.ITER.GT.1)THEN
c          PSEUDO = ABS(DEE/EELECT1).LT.0.05D0.AND.ABS(DPMAX).LT.0.1D0
           PSEUDO = ABS(DEE/ETOT1).LT.0.05D0.AND.ABS(DPMAX).LT.0.1D0
         ENDIF
C
C       DO A FULL DIAGONALIZATION EVERY 5 OR 10 ITERATIONS, DEPENDING
C       ON WHETHER A DYNAMIC LEVEL SHIFT IS BEING USED.
C
         IF(PSEUDO.AND.ITER.GT.10)THEN
           IF(MOD(ITER,10).EQ.0)THEN
             PSEUDO = .FALSE.
           ELSEIF(DOSHFT.AND.MOD(ITER,5).EQ.0)THEN
             PSEUDO = .FALSE.
           ENDIF
         ENDIF
C
C       KEEP TRACK OF ENERGY OSCILLATIONS.
C
        IF(DEE*DEPREV.LT.0.0D0.AND.ITER.GT.5) NOSCIL = NOSCIL + 1
C
C       FOCK MATRIX MIXING SCHEME TO ACCELERATE CONVERGENCE:
C
        if ((.not.fdmx).and.(.not.pvecnow)) then
C     .    don't mess with the Fock and density matrix when a
C     .    frozen calc. on the guess density matrix is performed
C           CALL FMIX(ITER,ICNVRG)
C
!           IF(ICNVRG.EQ.1)THEN

!              if (wrtscr) then
C-RDC                  write(iscr,'(" converged on FMIX DE")')
!              endif
!              if (pvec) then
!                 if (.not.pvecnow) then
C-RDC                    if (wrtscr)
C-RDC      $                   write(iscr,'(
C-RDC      $                   /"will perform one extra iteration ",
C-RDC      $                   "to write out eigenvectors"/)')
!                    pvecnow = .true.
!                    call setprtvec
!                 endif
!              else
!                 goto 1000
!              endif

!           ENDIF

C
C     .    DYNAMIC OR STATIC LEVEL SHIFT:
C
            IF(DOSHFT)THEN
               IF(ITER.EQ.1.OR..NOT.STAND)THEN
                  EGAP = 0.0D0
               ELSE
C
C     .          USE PREVIOUS HOMO-LUMO GAP IF WE'RE DOING A
C     .          PSEUDO-DIAGONALIZATION.
C
                  IF(.NOT.PSEUDO)THEN
                     NOCC = NELEC/2
                     EGAP = EVAL(NOCC+1) - EVAL(NOCC)
                  ENDIF
               ENDIF
               IF(EGAP.LT.ASHIFT)THEN
                  ESHIFT = 0.5D0*(EGAP - ASHIFT)
C
C     .          FSHIFT WILL NOT ACTUALLY DO A SHIFT IF THE DENSITY MATRIX
C     .          HAS JUST BEEN INITIALIZED TO DIAGONAL FORM (IDIAG.EQ.1).
C
                  CALL FSHIFT(IDIAG,ESHIFT)
               ENDIF
            ENDIF
        endif
C
C       MOSUB WILL BUILD NEW DENSITY MATRIX, SO SAVE OLD ONE TO CHECK
C       FOR CONVERGENCE.
C
        DO 290 II=1,IIMAX
           PIIOLD(II) = PDIAG(II)
 290    CONTINUE
        DO 300 IJ=1,IJMAX
           PIJOLD(IJ) = PDIAT(IJ)
 300    CONTINUE
C
C     DETERMINE MO'S AND DENSITY MATRIX.
C

        if (wrtscr) then
           icall = 1
           if(pseudo) icall = 2
        endif

C
C       A GOOD ESTIMATE OF FERMI ENERGY EXISTS IF IDIAG.NE.1
C
        GOODF = IDIAG.NE.1
        if (frozenmc) then
           CALL MOSUBfdm(GOODF,ITER,PSEUDO,IERROR)
        elseif (fdmx) then
           CALL MOSUBfdmx(GOODF,ITER,PSEUDO,pvecnow,IERROR) 
        else
           CALL MOSUB(GOODF,ITER,PSEUDO,pvecnow,IERROR)
        endif
        if (pvecnow) then
C-RDC            if (wrtscr) write(iscr,'("wrote eigenvectors")')
           lastit = .true.
           goto 10
        endif

        ndiag = ndiag + 1
        IF(IERROR.NE.0) RETURN
C
C       CHECK FOR CONVERGENCE IN DENSITY MATRIX AND ENERGY.
C

        DPOLD = DPMAX
        DPMAX = 0.0D0
        DO 370 II=1,IIMAX
           DII = PDIAG(II) - PIIOLD(II)
           IF(ABS(DII).GT.ABS(DPMAX)) DPMAX = DII
 370    CONTINUE
        DO 380 IJ=1,IJMAX
           DIJ = PDIAT(IJ) - PIJOLD(IJ)
           IF(ABS(DIJ).GT.ABS(DPMAX)) DPMAX = DIJ
 380    CONTINUE
C
C     . DENSITY MATRIX MIXING SCHEME TO ACCELERATE CONVERGENCE:
C

        if (.not.fdmx) then
C     .    don't mess with the density matrix when a frozen calculation
C     .    on guess density matrix is performed
           CALL PMIX(ITER)
C
C     .    USE A 50-50 DENSITY MATRIX MIXTURE ON FIRST ITERATION TO
C     .    AVOID OSCILLATIONS.
C
           IF(DOSHFT.AND.ITER.EQ.1)THEN

              if (wrtscr) then
C-RDC                  write(iscr,'(" fractional mix")')
              endif

              POLD = 0.5D0
              PNEW = 0.5D0
              DO 390 II=1,IIMAX
                 PDIAG(II) = PNEW*PDIAG(II) + POLD*PIIOLD(II)
 390          CONTINUE
              DO 400 IJ=1,IJMAX
                 PDIAT(IJ) = PNEW*PDIAT(IJ) + POLD*PIJOLD(IJ)
 400          CONTINUE
           ENDIF
        endif
        if (iIDump.ne.0) then
c          if(mod(k,iNDump).eq.0) call wrtdmx (ierror)
           if (ierror.ne.0) return
        endif

 500  CONTINUE
c      IERROR = 1
      WRITE(IOUT,'(/" NO CONVERGENCE IN SCF CALCULATION")')
 1000 FIRST = .FALSE.

 
      call etimer(tstop)
      if (wrtedmx) then
c        call wrtdmx(ierror)
         if (ierror.ne.0) return
      endif
        
c     if(frst_call.and.(.not.STAND)) ndiag = ndiag + 1
      tscf = tstop-tstart
      tscfav = tscf/ndiag
      tfockav = tfock/k

!     if (wrtscr) then
!        titer = (tstop-tstart)/ndiag
!        write(iscr,'(" ndiag=",i10,
!    $        /" average SCF cycle=",F11.3," seconds")')
!    $        ndiag, titer
!        if(STAND)then
!           write(iscr,'(" average full diag=",F11.3,
!    $           /" average pseudo diag=",f11.3)')
!    $           tcalls(1)/ncalls(1),tcalls(2)/ncalls(2)
!        endif
!     endif

      RETURN
      END