mopac606_nbo/src/diag.f

   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)
   5 C***********************************************************************
   6 C
   7 C   "FAST" DIAGONALISATION PROCEDURE.
   8 C
   9 C    ON INPUT FAO CONTAINS THE LOWER HALF TRIANGLE OF THE MATRIX TO BE
  10 C                         DIAGONALISED, PACKED.
  11 C             VECTOR  CONTAINS THE OLD EIGENVECTORS ON INPUT, THE NEW
  12 C             VECTORS ON EXITING.
  13 C             NOCC = NUMBER OF OCCUPIED MOLECULAR ORBITALS.
  14 C             EIG  = EIGENVALUES FROM AN EXACT DIAGONALISATION
  15 C             MDIM = DECLARED SIZE OF MATRIX "C".
  16 C             N = NUMBER OF ATOMIC ORBITALS IN BASIS SET
  17 C
  18 C  DIAG IS A PSEUDO-DIAGONALISATION PROCEDURE, IN THAT THE VECTORS THAT
  19 C       ARE GENERATED BY IT ARE MORE NEARLY ABLE TO BLOCK-DIAGONALISE
  20 C       THE FOCK MATRIX OVER MOLECULAR ORBITALS THAN THE STARTING
  21 C       VECTORS. IT MUST BE CONSIDERED PSEUDO FOR SEVERAL REASONS:
  22 C       (A) IT DOES NOT GENERATE EIGENVECTORS - THE SECULAR DETERMINANT
  23 C           IS NOT DIAGONALISED, ONLY THE OCCUPIED-VIRTUAL INTERSECTION.
  24 C       (B) MANY SMALL ELEMENTS IN THE SEC.DET. ARE IGNORED AS BEING TOO
  25 C           SMALL COMPARED WITH THE LARGEST ELEMENT.
  26 C       (C) WHEN ELEMENTS ARE ELIMINATED BY ROTATION, THE REST OF THE
  27 C           SEC. DET. IS ASSUMED NOT TO CHANGE, I.E. ELEMENTS CREATED
  28 C           ARE IGNORED.
  29 C       (D) THE ROTATION REQUIRED TO ELIMINATE THOSE ELEMENTS CONSIDERED
  30 C           SIGNIFICANT IS APPROXIMATED TO USING THE EIGENVALUES OF THE
  31 C           EXACT DIAGONALISATION THROUGHOUT THE REST OF THE ITERATIVE
  32 C           PROCEDURE.
  33 C
  34 C  (NOTE:- IN AN ITERATIVE PROCEDURE ALL THE APPROXIMATIONS PRESENT IN
  35 C          DIAG BECOME VALID AT SELF-CONSISTENCY, SELF-CONSISTENCY IS
  36 C          NOT SLOWED DOWN BY USE OF THESE APPROXIMATIONS)
  37 C
  38 C    REFERENCE:
  39 C             "FAST SEMIEMPIRICAL CALCULATIONS",
  40 C             STEWART. J.J.P., CSASZAR, P., PULAY, P., J. COMP. CHEM.,
  41 C             3, 227, (1982)
  42 C
  43 C***********************************************************************
  44       COMMON /SCRACH/ FMO(MORB2), XDUMY(MAXPAR**2-MORB2)
  45 C             FMO  IS A WORK-SPACE OF SIZE (N-NOCC)*NOCC, IT WILL HOLD
  46 C                  THE FOCK MOLECULAR ORBITAL INTERACTION MATRIX.
  47 C
  48 C  FIRST, CONSTRUCT THAT PART OF A SECULAR DETERMINANT OVER MOLECULAR
  49 C  ORBITALS WHICH CONNECTS THE OCCUPIED AND VIRTUAL SETS.
  50 C
  51 C***********************************************************************
  52 C
  53       LOGICAL FIRST
  54       DATA FIRST /.TRUE./
  55       IF(FIRST)THEN
  56          FIRST=.FALSE.
  57 C
  58 C   EPS IS THE SMALLEST NUMBER WHICH, WHEN ADDED TO 1.D0, IS NOT
  59 C   EQUAL TO 1.D0
  60          CALL EPSETA(EPS,ETA)
  61 C
  62 C   INCREASE EPS TO ALLOW FOR A LOT OF ROUND-OFF
  63 C
  64          BIGEPS=10.D0*SQRT(EPS)
  65       ENDIF
  66       TINY=0.D0
  67       LUMO=NOCC+1
  68       IJ=0
  69 C#      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
  93 C***********************************************************************
  94 C
  95 C   NOW DO A CRUDE 2 BY 2 ROTATION TO "ELIMINATE" SIGNIFICANT ELEMENTS
  96 C
  97 C***********************************************************************
  98 C#      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
 104 C
 105 C      BEGIN 2 X 2 ROTATIONS
 106 C
 107             A=EIG(J)
 108             B=EIG(I)
 109             C=FMO(IJ)
 110             D=A-B
 111 C
 112 C    USE BIGEPS TO DETERMINE WHETHER TO DO A 2 BY 2 ROTATION
 113 C
 114             IF(ABS(C/D).LT.BIGEPS) GOTO 80
 115 C
 116 C  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)
 120 C
 121 C      ROTATION OF PSEUDO-EIGENVECTORS
 122 C
 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
 131 C#      CALL TIMER('RETURNING')
 132       RETURN
 133       END