      SUBROUTINE ITER  (H, W, WJ, WK, EE, FULSCF,RAND)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DOUBLE PRECISION MECI
      DIMENSION H(*), W(*), WJ(*), WK(*)
      COMMON /FOKMAT/ F(MPACK), FB(MPACK)
      COMMON /DENSTY/ P(MPACK), PA(MPACK), PB(MPACK)
      COMMON /VECTOR/ C(MORB2),EIGS(MAXORB),CBETA(MORB2),EIGB(MAXORB)
      COMMON /GRADNT/ DUMY(MAXPAR),GNORM
      COMMON /LAST  / LAST
      COMMON /MESAGE/ IFLEPO,IITER
      COMMON /ATHEAT/ ATHEAT
      COMMON /ENUCLR/ ENUCLR
      COMMON /CITERM/ XI,XJ,XK
      COMMON /PATH  / LATOM,LPARAM,REACT(200)
      COMMON /NUMCAL/ NUMCAL
      COMMON /SCFTYP/ EMIN, LIMSCF
      COMMON /TIME  / TIME0
      LOGICAL FULSCF, RAND, LIMSCF
      DOUBLE PRECISION WJ, WK
C***********************************************************************
C
C     ITER GENERATES A SCF FIELD AND RETURNS THE ENERGY IN "ENERGY"
C
C THE MAIN ARRAYS USED IN ITER ARE:
C            P      ONLY EVER CONTAINS THE TOTAL DENSITY MATRIX
C            PA     ONLY EVER CONTAINS THE ALPHA DENSITY MATRIX
C            PB     ONLY EVER CONTAINS THE BETA DENSITY MATRIX
C            C      ONLY EVER CONTAINS THE EIGENVECTORS
C            H      ONLY EVER CONTAINS THE ONE-ELECTRON MATRIX
C            F      STARTS OFF CONTAINING THE ONE-ELECTRON MATRIX,
C                   AND IS USED TO HOLD THE FOCK MATRIX
C            W      ONLY EVER CONTAINS THE TWO-ELECTRON MATRIX
C
C THE MAIN INTEGERS CONSTANTS IN ITER ARE:
C
C            LINEAR SIZE OF PACKED TRIANGLE = NORBS*(NORBS+1)/2
C
C THE MAIN INTEGER VARIABLES ARE
C            NITER  NUMBER OF ITERATIONS EXECUTED
C
C  PRINCIPAL REFERENCES:
C
C   ON MNDO: "GROUND STATES OF MOLECULES. 38. THE MNDO METHOD.
C             APPROXIMATIONS AND PARAMETERS."
C             DEWAR, M.J.S., THIEL,W., J. AM. CHEM. SOC.,99,4899,(1977).
C   ON SHIFT: "THE DYNAMIC 'LEVEL SHIFT' METHOD FOR IMPROVING THE
C             CONVERGENCE OF THE SCF PROCEDURE", A. V. MITIN, J. COMP.
C             CHEM. 9, 107-110 (1988)
C   ON HALF-ELECTRON: "MINDO/3 COMPARISON OF THE GENERALIZED S.C.F.
C             COUPLING OPERATOR AND "HALF-ELECTRON" METHODS FOR
C             CALCULATING THE ENERGIES AND GEOMETRIES OF OPEN SHELL
C             SYSTEMS"
C             DEWAR, M.J.S., OLIVELLA, S., J. CHEM. SOC. FARA. II,
C             75,829,(1979).
C   ON PULAY'S CONVERGER: "CONVERGANCE ACCELERATION OF ITERATIVE
C             SEQUENCES. THE CASE OF SCF ITERATION", PULAY, P.,
C             CHEM. PHYS. LETT, 73, 393, (1980).
C   ON CNVG:  IT ENCORPORATES THE IMPROVED ITERATION SCHEME (IIS) BY
C             PIOTR BADZIAG & FRITZ SOLMS. ACCEPTED FOR PUBLISHING
C             IN COMPUTERS & CHEMISTRY
C   ON PSEUDODIAGONALISATION: "FAST SEMIEMPIRICAL CALCULATIONS",
C             STEWART. J.J.P., CSASZAR, P., PULAY, P., J. COMP. CHEM.,
C             3, 227, (1982)
C
C***********************************************************************
      DIMENSION POLD(MPACK), POLD2(MPACK), POLD3(MAXORB+400)
      DIMENSION PBOLD(MPACK), PBOLD2(MPACK), PBOLD3(MAXORB+400)
************************************************************************
*                                                                      *
*   PACK ALL THE ARRAYS USED BY PULAY INTO A COMMON BLOCK SO THAT THEY *
*   CAN BE USED BY THE C.I. DERIVATIVE, IF NEEDED                      *
*                                                                      *
************************************************************************
      COMMON /WORK3/POLD,POLD2,PBOLD,PBOLD2
      COMMON /WORK1/ AR1,AR2,AR3,AR4,BR1,BR2,BR3,BR4
      DIMENSION  AR1(2*NPULAY), AR2(2*NPULAY), AR3(2*NPULAY),
     1 AR4(2*NPULAY)
      DIMENSION  BR1(2*NPULAY), BR2(2*NPULAY), BR3(2*NPULAY),
     1 BR4(8*NPULAY)
      DIMENSION ESCF0(10)
      COMMON /PRECI / SELCON
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,
     2                NALPHA, NBETA, NCLOSE, NOPEN, NDUMY, FRACT
     3      /MOLORB/ DUMMY(MAXORB),PDIAG(MAXORB)
     4      /KEYWRD/ KEYWRD
     5      /NUMSCF/ NSCF
      SAVE ICALCN, DEBUG, PRTFOK, PRTEIG, PRTDEN, PRTENG, PRT1EL, ABPRT
      SAVE LINEAR, MINPRT, NEWDG, SCFCRT, PRTPL, PRTVEC, PL
      SAVE BSHIFT, PLTEST, ITRMAX, NA2EL, NA1EL, NB2EL,NB1EL
      SAVE IFILL, CAMKIN, CI, OKPULY, UHF, SCF1, OKNEWD, EXCITD, TIMES
      SAVE FORCE, ALLCON, TRANS, HALFE, W1, W2, RANDOM
      CHARACTER KEYWRD*241, ABPRT(3)*5
      LOGICAL PRTFOK,PRTEIG,PRTDEN, DEBUG, PRTENG, TIMES, CI
     1,UHF, NEWDG, SCF1, HALFE, FORCE, PRT1EL,PRTPL, OKNEWD
     2,EXCITD, MINPRT, FRST, BFRST, OKPULY, READY, PRTVEC,
     3CAMKIN, ALLCON, MAKEA, MAKEB, INCITR, CAPPS, TIMITR
      DATA ICALCN/0/, DEBUG/.FALSE./, PRTFOK/.FALSE./
      DATA PRTEIG/.FALSE./,PRTDEN,PRTENG/.FALSE.,.FALSE./
      DATA PRT1EL/.FALSE./
      DATA ABPRT/'     ','ALPHA',' BETA'/
C
C  INITIALIZE
C
      IFILL=0
      IHOMO=MAX(1,NCLOSE+NALPHA)
      IHOMOB=MAX(1,NCLOSE+NBETA)
      EOLD=1.D2
      READY=.FALSE.
      IF (ICALCN.NE.NUMCAL) THEN
         CALL EPSETA(EPS,ETA)
C
C  ULTIMATE SCF CRITERION: HEAT OF FORMATION CONVERGED WITHIN A FACTOR
C  OF 10 OF THE LIMITING PRECISION OF THE COMPUTER
C
         EPS=23.061D0*EPS*10.D0
         IRRR=5
         SHIFT=0.D0
         ICALCN=NUMCAL
         SHFMAX=20.D0
         LINEAR=(NORBS*(NORBS+1))/2
C
C    DEBUG KEY-WORDS WORKED OUT
C
         DEBUG=( INDEX(KEYWRD,'DEBUG') .NE. 0 )
         MINPRT=(INDEX(KEYWRD,'SADDLE')+
     1      LATOM .EQ.0 .OR. DEBUG)
         PRTEIG=( INDEX(KEYWRD,'EIGS') .NE. 0 )
         PRTPL =( INDEX(KEYWRD,' PL ')  .NE.0 )
         PRT1EL=( INDEX(KEYWRD,'1ELE') .NE.0 .AND. DEBUG)
         PRTDEN=( INDEX(KEYWRD,'DENS').NE.0 .AND. DEBUG)
         PRTFOK=( INDEX(KEYWRD,'FOCK') .NE. 0  .AND. DEBUG)
         PRTVEC=( INDEX(KEYWRD,'VECT') .NE. 0  .AND. DEBUG)
         DEBUG=( INDEX(KEYWRD,'ITER') .NE. 0 )
C
C INITIALIZE SOME LOGICALS AND CONSTANTS
C
         NEWDG =.FALSE.
         PLCHEK=0.005D0
         PL    =1.D0
         BSHIFT=0.D0
         SHIFT=1.D0
*
* SCFCRT IS MACHINE-PRECISION DEPENDENT
*
         SCFCRT=1.D-4
         ITRMAX = 200
         NA2EL=NCLOSE
         NA1EL=NALPHA+NOPEN
         NB2EL=0
         NB1EL=NBETA+NOPEN
C
C  USE KEY-WORDS TO ASSIGN VARIOUS CONSTANTS
C
         IF(INDEX(KEYWRD,'FILL').NE.0)
     1      IFILL=-READA(KEYWRD,INDEX(KEYWRD,'FILL'))
         IF(INDEX(KEYWRD,'SHIFT').NE.0)
     1      BSHIFT=-READA(KEYWRD,INDEX(KEYWRD,'SHIFT'))
         IF(BSHIFT.NE.0)TEN=BSHIFT
         IF(INDEX(KEYWRD,'ITRY').NE.0)
     1      ITRMAX=READA(KEYWRD,INDEX(KEYWRD,'ITRY'))
         CAMKIN=(INDEX(KEYWRD,'KING')+INDEX(KEYWRD,'CAMP') .NE. 0)
         CI=(INDEX(KEYWRD,'MICROS')+INDEX(KEYWRD,'C.I.') .NE. 0)
         OKPULY=.FALSE.
         OKPULY=(INDEX(KEYWRD,'PULAY').NE.0)
         UHF=(INDEX(KEYWRD,'UHF') .NE. 0)
         SCF1=(INDEX(KEYWRD,'1SCF') .NE. 0)
         OKNEWD=ABS(BSHIFT) .LT. 0.001D0
         IF(CAMKIN.AND.ABS(BSHIFT).GT.1.D-5) BSHIFT=4.44D0
         EXCITD=(INDEX(KEYWRD,'EXCI') .NE. 0)
         TIMES=(INDEX(KEYWRD,'TIMES') .NE. 0)
         TIMITR=(TIMES.AND.DEBUG)
         FORCE=(INDEX(KEYWRD,'FORCE') .NE. 0)
         ALLCON=(OKPULY.OR.CAMKIN)
C
C   DO WE NEED A CAPPED ATOM CORRECTION?
C
         J=0
         DO 10 I=1,NUMAT
   10    IF(NAT(I).EQ.102)J=J+1
         CAPPS=(J.GT.0)
         IITER=1
         TRANS=0.1D0
         IF(INDEX(KEYWRD,'RESTART')+INDEX(KEYWRD,'OLDENS')
     1      .NE. 0) THEN
            IF(INDEX(KEYWRD,'OLDENS').NE.0)
     1   OPEN(UNIT=10,FILE='FOR010',STATUS='UNKNOWN',FORM='UNFORMATTED')
            REWIND 10
            READ(10)(PA(I),I=1,LINEAR)
            IF( UHF) THEN
               READ(10)(PB(I),I=1,LINEAR)
               DO 20 I=1,LINEAR
                  POLD(I)=PA(I)
                  PBOLD(I)=PB(I)
   20          P(I)=PA(I)+PB(I)
            ELSE
               DO 30 I=1,LINEAR
                  PB(I)=PA(I)
                  PBOLD(I)=PA(I)
                  POLD(I)=PA(I)
   30          P(I)=PA(I)*2.D0
            ENDIF
         ELSE
            NSCF=0
            DO 40 I=1,LINEAR
               P(I)=0.D0
               PA(I)=0.D0
   40       PB(I)=0.D0
            W1=NA1EL/(NA1EL+1.D-6+NB1EL)
            W2=1.D0-W1
            IF(W1.LT.1.D-6)W1=0.5D0
            IF(W2.LT.1.D-6)W2=0.5D0
C
C  SLIGHTLY PERTURB THE DENSITY MATRIX IN CASE THE SYSTEM IS
C  TRAPPED IN A S**2 = 0 STATE.
C
            RANDOM=1.0D0
            IF(UHF.AND.NA1EL.EQ.NB1EL) RANDOM=1.1D0
            DO 50 I=1,NORBS
               J=(I*(I+1))/2
               P(J)=PDIAG(I)
               PA(J)=P(J)*W1*RANDOM
               RANDOM=1.D0/RANDOM
   50       PB(J)=P(J)*W2*RANDOM
            DO 60 I=1,LINEAR
               PBOLD(I)=PB(I)
   60       POLD(I)=PA(I)
         ENDIF
         HALFE=(NOPEN .NE. NCLOSE.AND.FRACT.NE.2.D0.AND.FRACT.NE.0.D0)
C
C   DETERMINE THE SELF-CONSISTENCY CRITERION
C
         IF(INDEX(KEYWRD,'PREC') .NE. 0)
     1                               SCFCRT=SCFCRT*0.01D0
         IF( INDEX(KEYWRD,'POLAR') + INDEX(KEYWRD,'NLLSQ') +
     1 INDEX(KEYWRD,'SIGMA') .NE. 0) SCFCRT=SCFCRT*0.001D0
         IF(FORCE)                   SCFCRT=SCFCRT*0.0001D0
         IF(NOPEN-NCLOSE.GT.4)       SCFCRT=SCFCRT*0.1D0
         SCFCRT=MAX(SCFCRT,1.D-12)
         IF(INDEX(KEYWRD,'POLAR').NE.0)SCFCRT=1.D-11
C
C  THE USER CAN STATE THE SCF CRITERION, IF DESIRED.
C
         I=INDEX(KEYWRD,'SCFCRT')
         IF(I.NE.0) THEN
            SCFCRT=READA(KEYWRD,I)
            WRITE(6,'(''  SCF CRITERION ='',G14.4)')SCFCRT
            IF(SCFCRT.LT.1.D-12)
     1 WRITE(6,'(//2X,'' THERE IS A RISK OF INFINITE LOOPING WITH'',
     2'' THE SCFCRT LESS THAN 1.D-12'')')
         ELSE
            IF(DEBUG)WRITE(6,'(''  SCF CRITERION ='',G14.4)')SCFCRT
         ENDIF
         IF(.NOT.SCF1)LAST=0
C
C   END OF INITIALIZATION SECTION.
C
      ELSEIF(FORCE.AND.NSCF.GT.0.AND..NOT.UHF)THEN
C
C   RESET THE DENSITY MATRIX IF MECI HAS FORMED AN EXCITED STATE.  THIS
C   PREVENTS THE SCF GETTING TRAPPED ON AN EXCITED STATE, PARTICULARLY
C   IF THE PULAY CONVERGER IS USED.
C
         DO 70 I=1,LINEAR
   70    P(I)=2.D0*PA(I)
      ENDIF
C
C   INITIALIZATION OPERATIONS DONE EVERY TIME ITER IS CALLED
C
      MAKEA=.TRUE.
      MAKEB=.TRUE.
      IEMIN=0
      IEMAX=0
C
C  TURN OFF SHIFT IF NOT A FULL SCF.
C
      IF(.NOT.FULSCF) SHIFT=0.D0
      IF(NEWDG) NEWDG=(ABS(BSHIFT).LT.0.001D0)
      IF(LAST.EQ.1) NEWDG=.FALSE.
C
C   SELF-CONSISTENCY CRITERIA: SELCON IS IN KCAL/MOL, PLTEST IS
C   A LESS IMPORTANT TEST TO MAKE SURE THAT THE SELCON TEST IS NOT
C   PASSED 'BY ACCIDENT'
C                              IF GNORM IS LARGE, MAKE SELCON BIGGER
C
      SELCON=SCFCRT
      IF(.NOT. FORCE .AND. .NOT. HALFE) THEN
         IF(GNORM.GT.5.D0) SELCON=SCFCRT*GNORM*0.2D0
         IF(GNORM.GT.200.D0) SELCON=SCFCRT*40.D0
      ENDIF
      PLTEST=0.05D0*SQRT(SELCON)
C
C  SOMETIMES HEAT GOES SCF BUT DENSITY IS STILL FLUCTUATING IN UHF
C  IN WHICH CASE PAY LESS ATTENTION TO DENSITY MATRIX
C
      IF(NALPHA.NE.NBETA.AND.UHF)PLTEST=0.001D0
      IF(DEBUG)THEN
         WRITE(6,'(''  SELCON, PLTEST'',3G16.7)')SELCON, PLTEST
      ENDIF
      TITER1=SECOND()
      IF(PRT1EL) THEN
         WRITE(6,'(//10X,''ONE-ELECTRON MATRIX AT ENTRANCE TO ITER'')')
         CALL VECPRT(H,NORBS)
      ENDIF
      IREDY=1
   80 NITER=0
      TIME1=SECOND()
      FRST=.TRUE.
      IF(CAMKIN) THEN
         MODEA=1
         MODEB=1
      ELSE
         MODEA=0
         MODEB=0
      ENDIF
      BFRST=.TRUE.
**********************************************************************
*                                                                    *
*                                                                    *
*                START THE SCF LOOP HERE                             *
*                                                                    *
*                                                                    *
**********************************************************************
      INCITR=.TRUE.
   90 INCITR=(MODEA.NE.3.AND.MODEB.NE.3)
      IF(INCITR)NITER=NITER+1
      IF(TIMITR)THEN
         TITER=SECOND()
         WRITE(6,*)
         WRITE(6,'(A,F7.2)')'     TIME FOR ITERATION:', TITER-TITER0
         TITER0=TITER
      ENDIF
      IF(NITER.GT.ITRMAX-10.AND..NOT.ALLCON) THEN
************************************************************************
*                                                                      *
*                   SWITCH ON ALL CONVERGERS                           *
*                                                                      *
************************************************************************
         WRITE(6,'(//,'' ALL CONVERGERS ARE NOW FORCED ON'',/
     1          '' SHIFT=10, PULAY ON, CAMP-KING ON'',/
     2          '' AND ITERATION COUNTER RESET'',//)')
         ALLCON=.TRUE.
         BSHIFT=4.44D0
         IREDY=-4
         EOLD=100.D0
         OKPULY=.TRUE.
         NEWDG=.FALSE.
         CAMKIN=(.NOT.HALFE)
         GOTO 80
      ENDIF
************************************************************************
*                                                                      *
*                        MAKE THE ALPHA FOCK MATRIX                    *
*                                                                      *
************************************************************************
      IF(ABS(SHIFT).GT.1.D-10.AND.BSHIFT .NE. 0.D0) THEN
         L=0
         IF(NITER.GT.1)THEN
            IF(NEWDG.AND..NOT.(HALFE.OR.CAMKIN))THEN
C
C  SHIFT WILL APPLY TO THE VIRTUAL ENERGY LEVELS USED IN THE
C  PSEUDODIAGONALIIZATION. IF DIFF IS -VE, GOOD, THEN LOWER THE
C  HOMO-LUMO GAP BY 0.1EV, OTHERWISE INCREASE IT.
               IF(DIFF.GT.0)THEN
                  SHIFT=1.D0
C
C IF THE PSEUDODIAGONALIZATION APPROXIMATION -- THAT THE WAVEFUNCTION
C IS ALMOST STABLE -- IS INVALID, TURN OFF NEWDG
                  IF(DIFF.GT.1)NEWDG=.FALSE.
               ELSE
                  SHIFT=-0.1D0
               ENDIF
            ELSE
               SHIFT=TEN+EIGS(IHOMO+1)-EIGS(IHOMO)+SHIFT
            ENDIF
            IF(DIFF.GT.0.D0) THEN
               IF(SHIFT.GT.4.D0)SHFMAX=4.5D0
               IF(SHIFT.GT.SHFMAX)SHFMAX=MAX(SHFMAX-0.5D0,0.D0)
            ENDIF
C
C   IF SYSTEM GOES UNSTABLE, LIMIT SHIFT TO THE RANGE -INFINITY - SHFMAX
C   BUT IF SYSTEM IS STABLE, LIMIT SHIFT TO THE RANGE -INFINITY - +20
C
            SHIFT=MAX(-20.D0,MIN(SHFMAX,SHIFT))
            IF(ABS(SHIFT-SHFMAX).LT.1.D-5)SHFMAX=SHFMAX+0.01D0
C
C  THE CAMP-KING AND PULAY CONVERGES NEED A CONSTANT SHIFT.
C  IF THE SHIFT IS ALLOWED TO VARY, THESE CONVERGERS WILL NOT
C  WORK PROPERLY.
C
            IF(OKPULY.OR.ABS(BSHIFT-4.44D0).LT.1.D-5)THEN
               SHIFT=-8.D0
               IF(NEWDG) SHIFT=0
            ENDIF
            IF(UHF)THEN
               IF(NEWDG.AND..NOT.(HALFE.OR.CAMKIN))THEN
                  SHIFTB=TEN-TENOLD
               ELSE
                  SHIFTB=TEN+EIGS(IHOMOB+1)-EIGS(IHOMOB)+SHIFTB
               ENDIF
               IF(DIFF.GT.0.D0)SHIFTB=MIN(4.D0,SHIFTB)
               SHIFTB=MAX(-20.D0,MIN(SHFMAX,SHIFTB))
               IF(OKPULY.OR.ABS(BSHIFT-4.44D0).LT.1.D-5)THEN
                  SHIFTB=-8.D0
                  IF(NEWDG)SHIFT=0
               ENDIF
               DO 100 I=IHOMOB+1,NORBS
  100          EIGB(I)=EIGB(I)+SHIFTB
            ELSE
            ENDIF
         ENDIF
         TENOLD=TEN
         IF(PL.GT.PLCHEK)THEN
            SHFTBO=SHIFTB
            SHFTO=SHIFT
         ELSE
            SHIFTB=SHFTBO
            SHIFT=SHFTO
         ENDIF
         DO 110 I=IHOMO+1,NORBS
  110    EIGS(I)=EIGS(I)+SHIFT
         DO 130 I=1,NORBS
            DO 120 J=1,I
               L=L+1
  120       F(L)=H(L)+SHIFT*PA(L)
  130    F(L)=F(L)-SHIFT
      ELSEIF (I.EQ.77.AND.LAST.EQ.0.AND.NITER.LT.2.AND.FULSCF)THEN
C
C  SLIGHTLY PERTURB THE FOCK MATRIX IN CASE THE SYSTEM IS
C  TRAPPED IN A METASTABLE EXCITED ELECTRONIC STATE
C
         RANDOM=0.001D0
         DO 140 I=1,LINEAR
            RANDOM=-RANDOM
  140    F(I)=H(I)+RANDOM
      ELSE
         DO 150 I=1,LINEAR
  150    F(I)=H(I)
      ENDIF
  160 CONTINUE
      IF(TIMITR)THEN
         T0=SECOND()
         WRITE(6,'(A,F7.2)')' LOAD FOCK MAT. INTEGRAL',T0-TITER0
      ENDIF
C#      CALL TIMER('BEFORE FOCK2')
      CALL FOCK2(F,P,PA,W, WJ, WK,NUMAT,NFIRST,NMIDLE,NLAST)
C#      CALL TIMER('AFTER FOCK2')
C#      CALL TIMER('BEFORE FOCK1')
      CALL FOCK1(F,P,PA,PB)
C#      CALL TIMER('AFTER FOCK1')
      IF(TIMITR)THEN
         T0=SECOND()
         TF1=TF1+T0-T1
         WRITE(6,'(2(A,F7.2))')'  FOCK1:',T0-T1,'INTEGRAL:',T0-TITER0
      ENDIF
************************************************************************
*                                                                      *
*                        MAKE THE BETA FOCK MATRIX                     *
*                                                                      *
************************************************************************
      IF (UHF) THEN
         IF(SHIFTB .NE. 0.D0) THEN
            L=0
            DO 180 I=1,NORBS
               DO 170 J=1,I
                  L=L+1
  170          FB(L)=H(L)+SHIFTB*PB(L)
  180       FB(L)=FB(L)-SHIFTB
         ELSEIF (RAND.AND.LAST.EQ.0.AND.NITER.LT.2.AND.FULSCF)THEN
            RANDOM=0.001
            DO 190 I=1,LINEAR
               RANDOM=-RANDOM
  190       FB(I)=H(I)+RANDOM
         ELSE
            DO 200 I=1,LINEAR
  200       FB(I)=H(I)
         ENDIF
         CALL FOCK2(FB,P,PB,W, WJ, WK,NUMAT,NFIRST,NMIDLE,NLAST)
         CALL FOCK1(FB,P,PB,PA)
      ENDIF
      IF( .NOT. FULSCF) GOTO 380
      IF(PRTFOK) THEN
         WRITE(6,210)NITER
  210    FORMAT('   FOCK MATRIX ON ITERATION',I3)
         CALL VECPRT (F,NORBS)
      ENDIF
C
C   CODE THE FOLLOWING LINE IN PROPERLY SOMETIME
C   THIS OPERATION IS BELIEVED TO GIVE RISE TO A BETTER FOCK MATRIX
C   THAN THE CONVENTIONAL GUESS.
C
      IF(IRRR.EQ.0)THEN
         DO 220 I=1,NORBS
  220    F((I*(I+1))/2)=F((I*(I+1))/2)*0.5D0
         IRRR=2
      ENDIF
************************************************************************
*                                                                      *
*                        CALCULATE THE ENERGY IN KCAL/MOLE             *
*                                                                      *
************************************************************************
      IF (NITER .GE. ITRMAX) THEN
         IF(DIFF.LT.1.D-3.AND.PL.LT.1.D-4.AND..NOT.FORCE)THEN
            WRITE(6,'('' """""""""""""""UNABLE TO ACHIEVE SELF-CONSISTEN
     1CE, JOB CONTINUING'')')
            GOTO 380
         ENDIF
         IF(MINPRT)WRITE (6,230)
  230    FORMAT (//10X,'"""""""""""""UNABLE TO ACHIEVE SELF-CONSISTENCE'
     1,/)
         WRITE (6,240) DIFF,PL
  240    FORMAT (//,10X,'DELTAE= ',E12.4,5X,'DELTAP= ',E12.4,///)
         IFLEPO=9
         IITER=2
         CALL WRITE (TIME0,ESCF)
         STOP
      ENDIF
      EE=HELECT(NORBS,PA,H,F)
      IF(UHF)THEN
         EE=EE+HELECT(NORBS,PB,H,FB)
      ELSE
         EE=EE*2.D0
      ENDIF
      IF(CAPPS)EE=EE+CAPCOR(NAT,NFIRST,NLAST,NUMAT,P,H)
      IF(BSHIFT.NE.0.D0)
     1SCORR=SHIFT*(NOPEN-NCLOSE)*23.061D0*0.25D0*(FRACT*(2.D0-FRACT))
      ESCF=(EE+ENUCLR)*23.061D0+ATHEAT+SCORR
      IF(INCITR)THEN
         DIFF=ESCF-EOLD
         IF(DIFF.GT.0)THEN
            TEN=TEN-1.D0
         ELSE
            TEN=TEN*0.975D0+0.05D0
         ENDIF
C
C MAKE SURE SELF-CONSISTENCY TEST IS NOT MORE STRINGENT THAN THE
C COMPUTER CAN HANDLE
C
         SELLIM=MAX(SELCON,EPS*MAX(ABS(EE),1.D0))
C
C SCF TEST:  CHANGE IN HEAT OF FORMATION IN KCAL/MOL SHOULD BE
C            LESS THAN SELLIM.  THE OTHER TESTS ARE SAFETY MEASURES
C
         IF(.NOT.(NITER.GT.4.AND.(PL.EQ.0.D0.OR.PL.LT.PLTEST.AND.
     1   ABS(DIFF).LT.SELLIM) .AND. READY)) GOTO 270
************************************************************************
*                                                                      *
*          SELF-CONSISTENCY TEST, EXIT MODE FROM ITERATIONS            *
*                                                                      *
************************************************************************
  250    IF (ABS(SHIFT) .LT. 1.D-10) GOTO 380
         SHIFT=0.D0
         SHIFTB=0.D0
         DO 260 I=1,LINEAR
  260    F(I)=H(I)
         MAKEA=.TRUE.
         MAKEB=.TRUE.
         GOTO 160
  270    CONTINUE
***********************************************************************
***********************************************************************
         IF(LIMSCF.AND.EMIN.NE.0.D0.AND..NOT.(CI.OR.HALFE))THEN
C
C  THE FOLLOWING TESTS ARE INTENDED TO ALLOW A FAST EXIT FROM ITER
C  IF THE RESULT IS 'GOOD ENOUGH' FOR THE CURRENT STEP IN THE GEOMETRY
C  OPTIMIZATION
C
            IF(ESCF.LT.EMIN)THEN
C
C  THE ENERGY IS LOWER THAN THE PREVIOUS MINIMUM.  NOW CHECK THAT
C  IT IS CONSISTENTLY LOWER.
C
               IEMAX=0
               IEMIN=MIN(5,IEMIN+1)
               DO 280 I=2,IEMIN
  280          ESCF0(I-1)=ESCF0(I)
               ESCF0(IEMIN)=ESCF
C
C  IS THE DIFFERENCE IN ENERGY BETWEEN TWO ITERATIONS LESS THAN 5%
C  OF THE ENERGY GAIN FOR THIS GEOMETRY RELATIVE TO THE PREVIOUS
C  MINIMUM.
C
               IF(IEMIN.GT.3)THEN
                  DO 290 I=2,IEMIN
                     IF(ABS(ESCF0(I)-ESCF0(I-1)).GT.0.05D0*(EMIN-ESCF))
     1GOTO 320
  290             CONTINUE
C
C IS GOOD ENOUGH -- RAPID EXIT
C
                  IF(DEBUG) WRITE(6,*)
     1' RAPID EXIT BECAUSE ENERGY IS CONSISTENTLY LOWER'
                  GOTO 250
               ENDIF
            ELSE
C
C  THE ENERGY HAS RISEN ABOVE THAT OF THE PREVIOUS MINIMUM.
C  WE NEED TO CHECK WHETHER THIS IS A FLUKE OR IS THIS REALLY
C  A BAD GEOMETRY.
C
               IEMIN=0
               IEMAX=MIN(5,IEMAX+1)
               DO 300 I=2,IEMAX
  300          ESCF0(I-1)=ESCF0(I)
               ESCF0(IEMAX)=ESCF
C
C  IS THE DIFFERENCE IN ENERGY BETWEEN TWO ITERATIONS LESS THAN 5%
C  OF THE ENERGY LOST FOR THIS GEOMETRY RELATIVE TO THE PREVIOUS
C  MINIMUM.
C
               IF(IEMAX.GT.3)THEN
                  DO 310 I=2,IEMAX
                     IF(ABS(ESCF0(I)-ESCF0(I-1)).GT.0.05D0*(ESCF-EMIN))
     1GOTO 320
  310             CONTINUE
C
C IS GOOD ENOUGH -- RAPID EXIT
C
                  IF(DEBUG) WRITE(6,*)
     1' RAPID EXIT BECAUSE ENERGY IS CONSISTENTLY HIGHER'
                  GOTO 250
               ENDIF
            ENDIF
         ENDIF
  320    READY=(IREDY.GT.0.AND.(ABS(DIFF).LT.SELLIM*10.D0.OR.PL.EQ.0.D0)
     1)
         IREDY=IREDY+1
      ENDIF
      IF(PRTPL.OR.DEBUG.AND.NITER.GT.ITRMAX-20) THEN
         IF(ABS(ESCF).GT.99999.D0) ESCF=SIGN(9999.D0,ESCF)
         IF(ABS(DIFF).GT.9999.D0)DIFF=0.D0
         IF(INCITR)
     1    WRITE(6,'('' ITERATION'',I3,'' PLS='',2E10.3,'' ENERGY  '',
     2F14.7,'' DELTAE'',F13.7)')NITER,PL,PLB,ESCF,DIFF
      ENDIF
      IF(INCITR)EOLD=ESCF
************************************************************************
*                                                                      *
*                        INVOKE THE CAMP-KING CONVERGER                *
*                                                                      *
************************************************************************
      IF(NITER.GT.2 .AND. CAMKIN .AND. MAKEA)
     1CALL INTERP(NORBS,NA1EL,NORBS-NA1EL, MODEA, ESCF/23.061D0,
     2F, C, AR1, AR2, AR3, AR4, AR1)
      MAKEB=.FALSE.
      IF(MODEA.EQ.3)GOTO 340
      MAKEB=.TRUE.
      IF(TIMITR)THEN
         T0=SECOND()
         WRITE(6,'(2(A,F7.2))')' ADJUST DAMPER  INTEGRAL',T0-TITER0
      ENDIF
      IF( NEWDG ) THEN
************************************************************************
*                                                                      *
*                        INVOKE PULAY'S CONVERGER                      *
*                                                                      *
************************************************************************
         IF(OKPULY.AND.MAKEA.AND.IREDY.GT.1)
     1CALL PULAY(F,PA,NORBS,POLD,POLD2,POLD3,JALP,IALP,MPACK,FRST,PL)
************************************************************************
*                                                                      *
*           DIAGONALIZE THE ALPHA OR RHF SECULAR DETERMINANT           *
* WHERE POSSIBLE, USE THE PULAY-STEWART METHOD, OTHERWISE USE BEPPU'S  *
*                                                                      *
************************************************************************
         IF (HALFE.OR.CAMKIN) THEN
            CALL HQRII(F,NORBS,NORBS,EIGS,C)
         ELSE
C#      CALL TIMER('BEFORE DIAG')
            CALL DIAG (F,C,NA1EL,EIGS,NORBS,NORBS)
C#      CALL TIMER('AFTER DIAG')
         ENDIF
      ELSE
C#      CALL TIMER('BEFORE HQRII')
         CALL HQRII(F,NORBS,NORBS,EIGS,C)
C#      CALL TIMER('AFTER HQRII')
         IF(TIMITR)THEN
            T1=SECOND()
            WRITE(6,'(2(A,F7.2))')'  HQRII:',T1-T0,' INTEGRAL',T1-TITER0
         ENDIF
      ENDIF
      J=1
      IF(PRTVEC) THEN
         J=1
         IF(UHF)J=2
         WRITE(6,'(//10X,A,
     1'' EIGENVECTORS AND EIGENVALUES ON ITERATION'',I3)')
     2   ABPRT(J),NITER
         CALL MATOUT(C,EIGS,NORBS,NORBS,NORBS)
      ELSE
         IF (PRTEIG) WRITE(6,330)ABPRT(J),NITER,(EIGS(I),I=1,NORBS)
      ENDIF
  330 FORMAT(10X,A,'  EIGENVALUES ON ITERATION',I3,/10(6G13.6,/))
  340 IF(IFILL.NE.0)CALL SWAP(C,NORBS,NORBS,NA2EL,IFILL)
************************************************************************
*                                                                      *
*            CALCULATE THE ALPHA OR RHF DENSITY MATRIX                 *
*                                                                      *
************************************************************************
      IF(UHF)THEN
         CALL DENSIT( C,NORBS, NORBS, NA2EL,NA1EL, FRACT, PA, 1)
         IF(MODEA.NE.3.AND..NOT. (NEWDG.AND.OKPULY))
     1    CALL CNVG(PA, POLD, POLD2, NORBS, NITER, PL)
      ELSE
C#      CALL TIMER('BEFORE DENSIT')
         CALL DENSIT( C,NORBS, NORBS, NA2EL,NA1EL, FRACT, P, 1)
C#      CALL TIMER('AFTER DENSIT')
         IF(MODEA.NE.3.AND..NOT. (NEWDG.AND.OKPULY))THEN
C#      CALL TIMER('BEFORE CNVG')
            CALL CNVG(P, POLD, POLD2, NORBS, NITER, PL)
C#      CALL TIMER('AFTER CNVG')
         ENDIF
      ENDIF
************************************************************************
*                                                                      *
*                       UHF-SPECIFIC CODE                              *
*                                                                      *
************************************************************************
      IF( UHF )THEN
************************************************************************
*                                                                      *
*                        INVOKE THE CAMP-KING CONVERGER                *
*                                                                      *
************************************************************************
         IF(NITER.GT.2 .AND. CAMKIN .AND. MAKEB )
     1CALL INTERP(NORBS,NB1EL,NORBS-NB1EL, MODEB, ESCF/23.061D0,
     2FB, CBETA, BR1, BR2, BR3, BR4, BR1)
         MAKEA=.FALSE.
         IF(MODEB.EQ.3) GOTO 350
         MAKEA=.TRUE.
         IF( NEWDG ) THEN
************************************************************************
*                                                                      *
*                        INVOKE PULAY'S CONVERGER                      *
*                                                                      *
************************************************************************
            IF( OKPULY.AND.MAKEB.AND.IREDY.GT.1)
     1CALL PULAY(FB,PB,NORBS,PBOLD,PBOLD2,
     2PBOLD3,JBET,IBET,MPACK,BFRST,PLB)
************************************************************************
*                                                                      *
*           DIAGONALIZE THE ALPHA OR RHF SECULAR DETERMINANT           *
* WHERE POSSIBLE, USE THE PULAY-STEWART METHOD, OTHERWISE USE BEPPU'S  *
*                                                                      *
************************************************************************
            IF (HALFE.OR.CAMKIN) THEN
               CALL HQRII(FB,NORBS,NORBS,EIGB,CBETA)
            ELSE
               CALL DIAG (FB,CBETA,NB1EL,EIGB,NORBS,NORBS)
            ENDIF
         ELSE
            CALL HQRII(FB,NORBS,NORBS,EIGB,CBETA)
         ENDIF
         IF(PRTVEC) THEN
            WRITE(6,'(//10X,A,'' EIGENVECTORS AND EIGENVALUES ON '',
     1''ITERATION'',I3)')ABPRT(3),NITER
            CALL MATOUT(CBETA,EIGB,NORBS,NORBS,NORBS)
         ELSE
            IF (PRTEIG) WRITE(6,330)ABPRT(3),NITER,(EIGB(I),I=1,NORBS)
         ENDIF
************************************************************************
*                                                                      *
*                CALCULATE THE BETA DENSITY MATRIX                     *
*                                                                      *
************************************************************************
  350    CALL DENSIT( CBETA,NORBS, NORBS, NB2EL, NB1EL, FRACT, PB, 1)
         IF( .NOT. (NEWDG.AND.OKPULY))
     1CALL CNVG(PB, PBOLD, PBOLD2, NORBS, NITER, PLB)
      ENDIF
************************************************************************
*                                                                      *
*                   CALCULATE THE TOTAL DENSITY MATRIX                 *
*                                                                      *
************************************************************************
      IF(UHF) THEN
         DO 360 I=1,LINEAR
  360    P(I)=PA(I)+PB(I)
      ELSE
         DO 370 I=1,LINEAR
            PA(I)=P(I)*0.5D0
  370    PB(I)=PA(I)
      ENDIF
      IF(PRTDEN) THEN
         WRITE(6,'('' DENSITY MATRIX ON ITERATION'',I4)')NITER
         CALL VECPRT (P,NORBS)
      ENDIF
      OKNEWD=(PL.LT.SELLIM .OR. OKNEWD)
      NEWDG=(PL.LT.TRANS .AND. OKNEWD .OR. NEWDG)
      IF(PL.LT.TRANS*0.3333D0)OKNEWD=.TRUE.
      GO TO 90
**********************************************************************
*                                                                    *
*                                                                    *
*                      END THE SCF LOOP HERE                         *
*                NOW CALCULATE THE ELECTRONIC ENERGY                 *
*                                                                    *
*                                                                    *
**********************************************************************
*          SELF-CONSISTENCE ACHEIVED.
*
  380 EE=HELECT(NORBS,PA,H,F)
      IF(UHF) THEN
         EE=EE+HELECT(NORBS,PB,H,FB)
      ELSE
         EE=EE*2.D0 +
     1SHIFT*(NOPEN-NCLOSE)*23.061D0*0.25D0*(FRACT*(2.D0-FRACT))
      ENDIF
      IF(CAPPS)EE=EE+CAPCOR(NAT,NFIRST,NLAST,NUMAT,P,H)
C
C   NORMALLY THE EIGENVALUES ARE INCORRECT BECAUSE THE
C   PSEUDODIAGONALIZATION HAS BEEN USED.  IF THIS
C   IS THE LAST SCF, THEN DO AN EXACT DIAGONALIZATION
      IF( NSCF.EQ.0 .OR. LAST.EQ.1 .OR. CI .OR. HALFE ) THEN
C
C  PUT F AND FB INTO POLD IN ORDER TO NOT DESTROY F AND FB
C  AND DO EXACT DIAGONALISATIONS
         DO 390 I=1,LINEAR
  390    POLD(I)=F(I)
         CALL HQRII(POLD,NORBS,NORBS,EIGS,C)
         IF(UHF) THEN
            DO 400 I=1,LINEAR
  400       POLD(I)=FB(I)
            CALL HQRII(POLD,NORBS,NORBS,EIGB,CBETA)
            DO 410 I=1,LINEAR
  410       POLD(I)=PA(I)
         ELSE
            DO 420 I=1,LINEAR
  420       POLD(I)=P(I)
         ENDIF
         IF(CI.OR.HALFE) THEN
C#        CALL TIMER('BEFORE MECI')
            SUM=MECI(EIGS,C)
C#        CALL TIMER('AFTER MECI')
            EE=EE+SUM
            IF(PRTPL)THEN
               ESCF=(EE+ENUCLR)*23.061 +ATHEAT
               WRITE(6,'(27X,''AFTER MECI, ENERGY  '',F14.7)')ESCF
            ENDIF
         ENDIF
      ENDIF
      NSCF=NSCF+1
      TITER2=SECOND()
      IF(DEBUG)WRITE(6,'('' NO. OF ITERATIONS ='',I6)')NITER
C            IF(FORCE)  SCFCRT=1.D-5
      IF(ALLCON.AND.ABS(BSHIFT-4.44D0).LT.1.D-7)THEN
         CAMKIN=.FALSE.
         ALLCON=.FALSE.
         NEWDG=.FALSE.
         BSHIFT=-10.D0
         OKPULY=.FALSE.
      ENDIF
      SHIFT=1.D0
      IF(EMIN.EQ.0.D0)THEN
         EMIN=ESCF
      ELSE
         EMIN=MIN(EMIN,ESCF)
      ENDIF
      RETURN
      END
