      SUBROUTINE DIAG(FAO,VECTOR,NOCC,EIG,MDIM,N)
      implicit double precision(A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension FAO(*),VECTOR(MDIM,*),EIG(*),WS(MAXORB)
C***********************************************************************
C
C   "FAST" DIAGONALISATION PROCEDURE.
C
C    ON INPUT FAO CONTAINS THE LOWER HALF TRIANGLE OF THE MATRIX TO BE
C                         DIAGONALISED, PACKED.
C             VECTOR  CONTAINS THE OLD EIGENVECTORS ON INPUT, THE NEW
C             VECTORS ON EXITING.
C             NOCC = NUMBER OF OCCUPIED MOLECULAR ORBITALS.
C             EIG  = EIGENVALUES FROM AN EXACT DIAGONALISATION
C             MDIM = DECLARED SIZE OF MATRIX "C".
C             N = NUMBER OF ATOMIC ORBITALS IN BASIS SET
C
C  DIAG IS A PSEUDO-DIAGONALISATION PROCEDURE, IN THAT THE VECTORS THAT
C       ARE GENERATED BY IT ARE MORE NEARLY ABLE TO BLOCK-DIAGONALISE
C       THE FOCK MATRIX OVER MOLECULAR ORBITALS THAN THE STARTING
C       VECTORS. IT MUST BE CONSIDERED PSEUdo FOR SEVERAL REASONS:
C       (A) IT DOES NOT GENERATE EIGENVECTORS - THE SECULAR DETERMINANT
C           IS NOT DIAGONALISED, ONLY THE OCCUPIED-VIRTUAL INTERSECTION.
C       (B) MANY SMALL ELEMENTS IN THE SEC.DET. ARE IGNORED AS BEING TOO
C           SMALL COMPARED WITH THE LARGEST ELEMENT.
C       (C) WHEN ELEMENTS ARE ELIMINATED BY ROTATION, THE REST OF THE
C           SEC. DET. IS ASSUMED NOT TO CHANGE, I.E. ELEMENTS CREATED
C           ARE IGNORED.
C       (D) THE ROTATION REQUIRED TO ELIMINATE THOSE ELEMENTS CONSIDERED
C           SIGNIFICANT IS APPROXIMATED TO USING THE EIGENVALUES OF THE
C           EXACT DIAGONALISATION THROUGHOUT THE REST OF THE ITERATIVE
C           PROCEDURE.
C
C  (NOTE:- IN AN ITERATIVE PROCEDURE ALL THE APPROXIMATIONS PRESENT IN
C          DIAG BECOME VALID AT SELF-CONSISTENCY, SELF-CONSISTENCY IS
C          NOT SLOWED DOWN BY USE OF THESE APPROXIMATIONS)
C
C    REFERENCE:
C             "FAST SEMIEMPIRICAL CALCULATIONS",
C             STEWART. J.J.P., CSASZAR, P., PULAY, P., J. COMP. CHEM.,
C             3, 227, (1982)
C
C***********************************************************************
      common /SCRACH/ FMO(MORB2), XDUMY(MAXPAR**2-MORB2)
C             FMO  IS A WORK-SPACE OF SIZE (N-NOCC)*NOCC, IT WILL HOLD
C                  THE FOCK MOLECULAR ORBITAL INTERACTION MATRIX.
C
C  FIRST, CONSTRUCT THAT PART OF A SECULAR DETERMINANT OVER MOLECULAR
C  ORBITALS WHICH CONNECTS THE OCCUPIED AND VIRTUAL SETS.
C
C***********************************************************************
C
      logical FIRST
      data FIRST /.TRUE./
      if (FIRST) then
         FIRST=.FALSE.
C
C   EPS IS THE SMALLEST NUMBER WHICH, WHEN ADDED TO 1.D0, IS NOT
C   EQUAL TO 1.D0
         call EPSETA(EPS,ETA)
C
C   INCREASE EPS TO ALLOW FOR A LOT OF ROUND-OFF
C
         BIGEPS=10.D0*SQRT(EPS)
      end if
      TINY=0.D0
      LUMO=NOCC+1
      IJ=0
C#      call TIMER('SQUARING')
      do 60 I=LUMO,N
         KK=0
         do 30 J=1,N
            SUM=0.D0
            do 10 K=1,J
               KK=KK+1
   10       SUM=SUM+FAO(KK)*VECTOR(K,I)
            if (J.EQ.N) go to 30
            J1=J+1
            K2=KK
            do 20 K=J1,N
               K2=K2+K-1
   20       SUM=SUM+FAO(K2)*VECTOR(K,I)
   30    WS(J)=SUM
         do 50 J=1,NOCC
            IJ=IJ+1
            SUM=0.D0
            do 40 K=1,N
   40       SUM=SUM+WS(K)*VECTOR(K,J)
            if (TINY.LT.ABS(SUM)) TINY=ABS(SUM)
   50    FMO(IJ)=SUM
   60 continue
      TINY=0.05D0*TINY
C***********************************************************************
C
C   NOW do A CRUDE 2 BY 2 ROTATION TO "ELIMINATE" SIGNifICANT ELEMENTS
C
C***********************************************************************
C#      call TIMER('ROTATING')
      IJ=0
      do 90 I=LUMO,N
         do 80 J=1,NOCC
            IJ=IJ+1
            if (ABS(FMO(IJ)).LT.TINY) go to 80
C
C      BEGIN 2 X 2 ROTATIONS
C
            A=EIG(J)
            B=EIG(I)
            C=FMO(IJ)
            D=A-B
C
C    USE BIGEPS TO DETERMINE WHETHER TO do A 2 BY 2 ROTATION
C
            if (ABS(C/D).LT.BIGEPS) go to 80
C
C  AT THIS POINT WE KNOW THAT
            E=SIGN(SQRT(4.D0*C*C+D*D),D)
            ALPHA=SQRT(0.5D0*(1.D0+D/E))
            BETA=-SIGN(SQRT(1.D0-ALPHA*ALPHA),C)
C
C      ROTATION OF PSEUDO-EIGENVECTORS
C
            do 70 M=1,N
               A=VECTOR(M,J)
               B=VECTOR(M,I)
               VECTOR(M,J)=ALPHA*A+BETA*B
               VECTOR(M,I)=ALPHA*B-BETA*A
   70       continue
   80    continue
   90 continue
C#      call TIMER('returnING')
      return
      end
