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