C
C $Id: diagp.F,v 1.2 1998/07/16 16:39:35 jjv5 Exp arjan $
C
C------------------------------------------------------------------------
      SUBROUTINE DIAGP(NORBS,NOCC,F,E,FOV,V,C)
C
C     PSEUDO-DIAGONALIZATION ROUTINE FOR THE MATRIX F.  THIS ROUTINE
C     IS USED WHEN EIGENVALUES IN SCF CALCULATION HAVE SUFFICIENTLY
C     CONVERGED.  THEIR VALUES ARE USED TO APPROXIMATELY ELIMINATE
C     SIGNIFICANT OFF-DIAGONAL ELEMENTS IN THE OCCUPIED-VIRTUAL
C     INTERSECTION OF THE MO'S WITH F.  SEE J.J.P. STEWART, C. CSASZAR,
C     AND P. PULAY, J. COMP. CHEM, VOL. 3, 227-228 (1982).
C
C     NORBS = THE TOTAL NUMBER OF MOLECULAR ORBITALS.
C
C      NOCC = THE NUMBER OF OCCUPIED MOLECULAR ORBITALS.
C
C         F = MATRIX TO BE PSEUDODIAGONALIZED.  THE STORAGE IS LINEAR
C             AND REQUIRES AT LEAST NORBS*NORBS SLOTS.  THE 2-D TO 1-D
C             CORRESPONDENCE IS (I,J) --> (I+(J-1)*NORBS).  ONLY THE
C             LOWER TRIANGLE (I.GE.J) IS NEEDED AND ONLY THE UPPER
C             TRIANGLE IS ALTERED.
C
C         E = EIGENVALUES OF F.  REQUIRES NORBS SLOTS.
C
C       FOV = TEMPORARY STORAGE FOR PARTIALLY DIAGONALIZED VERSION OF F.
C             REQUIRES AS MUCH STORAGE AS F.
C
C         V = TEMPORARY STORAGE.  REQUIRES NORBS SLOTS.
C
C         C = EIGENVECTORS FOR F.  STORAGE IS THE SAME AS F.  THIS ROUTINE
C             ROTATES THE EIGENVECTORS TO HELP ELIMINATE THE OCCUPIED-
C             VIRTUAL INTERSECTION.
C
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION F(*),E(*),FOV(*),V(*),C(*)
C

      EPSLON = 1.0D-8
      TOLERA = 0.05D0
C
C     COPY LOWER TRIANGLE OF F TO UPPER TRIANGLE FOR QUICKER ACCESS.
C
      ITERM = 0
      DO 20 I=1,NORBS
        JTERM = 0
        DO 10 J=1,I
          F(J+ITERM) = F(I+JTERM)
          JTERM = JTERM + NORBS
 10     CONTINUE
        ITERM = ITERM + NORBS
 20   CONTINUE
C
C     FORM THE MATRIX FOV FROM THE OCCUPIED-VIRTUAL INTERSECTION OF
C     THE MOLECULAR ORBITALS:
C
C                  FOV = C(OCC)'*F*C(VIRT)
C
      FIJMAX = 0.0D0
      JTERM = NOCC*NORBS
      DO 180 J=NOCC+1,NORBS
        ITERM = 0
        DO 120 I=1,NORBS
          VI = 0.0D0
          DO 100 K=1,NORBS
            VI = VI + F(K+ITERM)*C(K+JTERM)
 100      CONTINUE
          V(I) = VI
          ITERM = ITERM + NORBS
 120    CONTINUE
        ITERM = 0
        DO 160 I=1,NOCC
          FIJ = 0.0D0
          DO 140 K=1,NORBS
            FIJ = FIJ + C(K+ITERM)*V(K)
 140      CONTINUE
          IF(ABS(FIJ).GT.FIJMAX) FIJMAX = ABS(FIJ)
          FOV(I+JTERM) = FIJ
          ITERM = ITERM + NORBS
 160    CONTINUE
        JTERM = JTERM + NORBS
 180  CONTINUE
      THRESH = TOLERA*FIJMAX
C
C     DO A PSEUDO-JACOBI ROTATION TO APPROXIMATELY ELIMINATE SIGNIFICANT
C     ELEMENTS OF FOV.
C
      ITERM = 0
      DO 240 I=1,NOCC
        JTERM = NOCC*NORBS
        DO 220 J=NOCC+1,NORBS
          IF(ABS(FOV(I+JTERM)).LT.THRESH) GO TO 210
          EOCC = E(I)
          EVIRT = E(J)
          FIJ = FOV(I+JTERM)
          D = EOCC - EVIRT
          IF(D.NE.0.0D0)THEN
            IF(ABS(FIJ/D).LT.EPSLON) GO TO 210
          ENDIF
          D = D*0.5D0
          SIGD = SIGN(1.0D0,D)
          W = -SIGD*FIJ/DSQRT(FIJ**2 + D**2)
          SINE = W/DSQRT(2.0D0*(1.0D0 + DSQRT(1.0D0 - W*W)))
          COSINE = DSQRT(1.0D0 - SINE**2)
C
C         ROTATE PSEUDO-EIGENVECTORS.
C
          DO 200 K=1,NORBS
            KI = K + ITERM
            KJ = K + JTERM
            CKI = C(KI)
            CKJ = C(KJ)
            C(KI) = CKI*COSINE - CKJ*SINE
            C(KJ) = CKI*SINE + CKJ*COSINE
 200      CONTINUE
 210      JTERM = JTERM + NORBS
 220    CONTINUE
        ITERM = ITERM + NORBS
 240  CONTINUE
      RETURN
      END