1      SUBROUTINE DIAG(FAO,VECTOR,NOCC,EIG,MDIM,N)
2      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3      INCLUDE 'SIZES'
4      DIMENSION FAO(*),VECTOR(MDIM,*),EIG(*),WS(MAXORB)
5C***********************************************************************
6C
7C   "FAST" DIAGONALISATION PROCEDURE.
8C
9C    ON INPUT FAO CONTAINS THE LOWER HALF TRIANGLE OF THE MATRIX TO BE
10C                         DIAGONALISED, PACKED.
11C             VECTOR  CONTAINS THE OLD EIGENVECTORS ON INPUT, THE NEW
12C             VECTORS ON EXITING.
13C             NOCC = NUMBER OF OCCUPIED MOLECULAR ORBITALS.
14C             EIG  = EIGENVALUES FROM AN EXACT DIAGONALISATION
15C             MDIM = DECLARED SIZE OF MATRIX "C".
16C             N = NUMBER OF ATOMIC ORBITALS IN BASIS SET
17C
18C  DIAG IS A PSEUDO-DIAGONALISATION PROCEDURE, IN THAT THE VECTORS THAT
19C       ARE GENERATED BY IT ARE MORE NEARLY ABLE TO BLOCK-DIAGONALISE
20C       THE FOCK MATRIX OVER MOLECULAR ORBITALS THAN THE STARTING
21C       VECTORS. IT MUST BE CONSIDERED PSEUDO FOR SEVERAL REASONS:
22C       (A) IT DOES NOT GENERATE EIGENVECTORS - THE SECULAR DETERMINANT
23C           IS NOT DIAGONALISED, ONLY THE OCCUPIED-VIRTUAL INTERSECTION.
24C       (B) MANY SMALL ELEMENTS IN THE SEC.DET. ARE IGNORED AS BEING TOO
25C           SMALL COMPARED WITH THE LARGEST ELEMENT.
26C       (C) WHEN ELEMENTS ARE ELIMINATED BY ROTATION, THE REST OF THE
27C           SEC. DET. IS ASSUMED NOT TO CHANGE, I.E. ELEMENTS CREATED
28C           ARE IGNORED.
29C       (D) THE ROTATION REQUIRED TO ELIMINATE THOSE ELEMENTS CONSIDERED
30C           SIGNIFICANT IS APPROXIMATED TO USING THE EIGENVALUES OF THE
31C           EXACT DIAGONALISATION THROUGHOUT THE REST OF THE ITERATIVE
32C           PROCEDURE.
33C
34C  (NOTE:- IN AN ITERATIVE PROCEDURE ALL THE APPROXIMATIONS PRESENT IN
35C          DIAG BECOME VALID AT SELF-CONSISTENCY, SELF-CONSISTENCY IS
36C          NOT SLOWED DOWN BY USE OF THESE APPROXIMATIONS)
37C
38C    REFERENCE:
39C             "FAST SEMIEMPIRICAL CALCULATIONS",
40C             STEWART. J.J.P., CSASZAR, P., PULAY, P., J. COMP. CHEM.,
41C             3, 227, (1982)
42C
43C***********************************************************************
44      COMMON /SCRACH/ FMO(MORB2), XDUMY(MAXPAR**2-MORB2)
45C             FMO  IS A WORK-SPACE OF SIZE (N-NOCC)*NOCC, IT WILL HOLD
46C                  THE FOCK MOLECULAR ORBITAL INTERACTION MATRIX.
47C
48C  FIRST, CONSTRUCT THAT PART OF A SECULAR DETERMINANT OVER MOLECULAR
49C  ORBITALS WHICH CONNECTS THE OCCUPIED AND VIRTUAL SETS.
50C
51C***********************************************************************
52C
53      LOGICAL FIRST
54      DATA FIRST /.TRUE./
55      IF(FIRST)THEN
56         FIRST=.FALSE.
57C
58C   EPS IS THE SMALLEST NUMBER WHICH, WHEN ADDED TO 1.D0, IS NOT
59C   EQUAL TO 1.D0
60         CALL EPSETA(EPS,ETA)
61C
62C   INCREASE EPS TO ALLOW FOR A LOT OF ROUND-OFF
63C
64         BIGEPS=10.D0*SQRT(EPS)
65      ENDIF
66      TINY=0.D0
67      LUMO=NOCC+1
68      IJ=0
69C#      CALL TIMER('SQUARING')
70      DO 60 I=LUMO,N
71         KK=0
72         DO 30 J=1,N
73            SUM=0.D0
74            DO 10 K=1,J
75               KK=KK+1
76   10       SUM=SUM+FAO(KK)*VECTOR(K,I)
77            IF(J.EQ.N) GOTO 30
78            J1=J+1
79            K2=KK
80            DO 20 K=J1,N
81               K2=K2+K-1
82   20       SUM=SUM+FAO(K2)*VECTOR(K,I)
83   30    WS(J)=SUM
84         DO 50 J=1,NOCC
85            IJ=IJ+1
86            SUM=0.D0
87            DO 40 K=1,N
88   40       SUM=SUM+WS(K)*VECTOR(K,J)
89            IF(TINY.LT.ABS(SUM)) TINY=ABS(SUM)
90   50    FMO(IJ)=SUM
91   60 CONTINUE
92      TINY=0.05D0*TINY
93C***********************************************************************
94C
95C   NOW DO A CRUDE 2 BY 2 ROTATION TO "ELIMINATE" SIGNIFICANT ELEMENTS
96C
97C***********************************************************************
98C#      CALL TIMER('ROTATING')
99      IJ=0
100      DO 90 I=LUMO,N
101         DO 80 J=1,NOCC
102            IJ=IJ+1
103            IF(ABS(FMO(IJ)).LT.TINY) GOTO 80
104C
105C      BEGIN 2 X 2 ROTATIONS
106C
107            A=EIG(J)
108            B=EIG(I)
109            C=FMO(IJ)
110            D=A-B
111C
112C    USE BIGEPS TO DETERMINE WHETHER TO DO A 2 BY 2 ROTATION
113C
114            IF(ABS(C/D).LT.BIGEPS) GOTO 80
115C
116C  AT THIS POINT WE KNOW THAT
117            E=SIGN(SQRT(4.D0*C*C+D*D),D)
118            ALPHA=SQRT(0.5D0*(1.D0+D/E))
119            BETA=-SIGN(SQRT(1.D0-ALPHA*ALPHA),C)
120C
121C      ROTATION OF PSEUDO-EIGENVECTORS
122C
123            DO 70 M=1,N
124               A=VECTOR(M,J)
125               B=VECTOR(M,I)
126               VECTOR(M,J)=ALPHA*A+BETA*B
127               VECTOR(M,I)=ALPHA*B-BETA*A
128   70       CONTINUE
129   80    CONTINUE
130   90 CONTINUE
131C#      CALL TIMER('RETURNING')
132      RETURN
133      END
134