PROGRAM HMO DIMENSION P(1275), C(1275), ROOT(50), AH(50,50) DIMENSION VECT(50,50), FV1(50), FV2(50) LOGICAL DEGEN, IPFLAG CHARACTER ITITLE(2)*10,IPROM*3,INEW*3,IPRINT*3,ICHANG*3 C C HUCKEL MOLECULAR ORBITAL PROGRAM C VERSION DATED APRIL 1984. C updated for VAX VMS October, 1990. C updated for Unix November, 1992. C C H. S. RZEPA, C DEPARTMENT OF CHEMISTRY, C IMPERIAL COLLEGE, C LONDON, SW7 2AT, UK. C C TEL. 071 225 8339 E-mail: rzepa@ic.ac.uk C C C C DEFINE I/O BUFFERS C OPEN(5,FILE='INPUT') C OPEN(6,FILE='OUTPUT') C C ZERO HUCKEL MATRIX DO 9 J9=1,1275 9 C(J9) = 0.0 10 IPFLAG = .FALSE. 3001 WRITE(*,3000) 3000 FORMAT(' TYPE Y IF YOU WANT DETAILED PROMPTS ISSUED DURING THE PRO 1GRAM, N IF YOU'/' ONLY WANT SHORT PROMPTS OR PRESS CTRL+d TO EXIT 2(=N)') READ (*,9932,END=140,ERR=3001) IPROM IF(IPROM.EQ.'Y' .OR. IPROM.EQ.'YES') IPFLAG = .TRUE. IF(IPROM.EQ.'y' .OR. IPROM.EQ.'yes') IPFLAG = .TRUE. IF(IPFLAG) WRITE(*,999) 999 FORMAT(' YOU WILL BE ASKED VARIOUS QUESTIONS ABOUT YOUR MOLECULE,W 1HICH'/' MOSTLY REQUIRE YOU TO TYPE Y (FOR YES) OR N (FOR NO).IF IN 2STEAD YOU ENTER'/' CTRL+d YOU WILL EITHER EXIT FROM THE PROGRAM, O 3R A DEFAULT OPTION WILL'/' BE SELECTED (SHOWN IN THE PRECEEDING PR 4OMPT) OR THE PROMPT WILL BE REPEATED.') C INITIALLY, READ IN A TITLE C IEXCIT = 0 DEGEN = .FALSE. ITITLE(1) = 'HMO CALCUL' ITITLE(2) = 'ATION ' WRITE (*,1000) READ (*,1015,END=140) (ITITLE(J),J=1,2) C C FIRST, READ IN THE ORDER OF YOUR SECULAR MATRIX. 726 WRITE (*,1010) IF(IPFLAG) WRITE(*,9056) 9056 FORMAT(' TYPE TWO INTEGERS, SEPARATING THE VALUES BY A COMMA, OR A 1T LEAST ONE SPACE.') 928 READ (*,*,END=925,ERR=925) IN,IEL IF(IN.GT.50) WRITE(*,725) 725 FORMAT(' THE PROGRAM CANNOT HANDLE MORE THAN 50 ATOMS') IF(IN.GT.50)GO TO 726 GO TO 926 925 WRITE(*,927) 927 FORMAT(' TYPE VALUES FOR THE NUMBER OF ATOMS AND ELECTRONS.') C CLOSE(5) C OPEN(5,FILE='INPUT') GO TO 928 C NEXT, READ IN HOW MANY ELECTRONS YOU HAVE. C NOCCO REPRESENTS THE NUMBER OF CLOSED SHELL FILLED ORBITALS 926 NOCCO = IEL/2 IOD = NOCCO*2 C IF THERE IS AN UNPAIRED ELECTRON, REDEFINE NOCCO TO C REPRESENT THE NUMBER OF OCCUPIED ORBITALS. IF (IEL.GT.IOD) NOCCO = NOCCO+1 IF (IEL.EQ.IOD) NCLOSE = NOCCO IF (IEL.GT.IOD) NCLOSE = NOCCO-1 IF (IEL.GT.IOD) IODD = NOCCO C C OPTION FOR CALCULATING AN EXCITED STATE. C IF (IEL.EQ.IOD) WRITE (*,1030) IF (IEL.EQ.IOD) READ (*,*,END=1031,ERR=1031) IEXCIT 1031 CONTINUE C 1031 CLOSE(5) C OPEN(5,FILE='INPUT') IF (IEXCIT.NE.0) NCLOSE = NOCCO-1 IF (IEXCIT.NE.0) NOCCO = NOCCO+1 IF (IEXCIT.NE.0) IODD = NOCCO-1 C C IF HMO MATRIX ALREADY IN, SKIP THE READ. C WRITE (*,1012) IF(IPFLAG) WRITE(*,2101) 2101 FORMAT(' TYPE N (DEFAULT) IF YOU WANT TO REPEAT YOUR PREVIOUS CALC 1ULATION CHANGING JUST'/ 2' ONE OR TWO ELEMENTS, OR IF YOU WANT TO ENTER A MATRIX OF ZEROS') INEW = 'N' C CALCULATE THE NUMBER OF MATRIX ELEMENTS THAT WILL HAVE C TO BE READ IN. NI = IN*(IN+1)/2 READ (*,9932,END=226,ERR=226) INEW 9932 FORMAT(A3) IF(INEW.EQ.'NO'. OR. INEW.EQ.'N') GO TO 226 IF(INEW.EQ.'no'. OR. INEW.EQ.'n') GO TO 226 C READ IN THE BOTTOM TRIANGLE OF YOUR MATRIX ELEMENTS. WRITE (*,1040) NI IF(IPFLAG) WRITE (*,1041) 1041 FORMAT(' IF YOU PRESS THE RETURN KEY AT THIS STAGE, A MATRIX CONSI 1STING ENTIRELY OF ZEROS'/' WILL BE SET UP. INDIVIDUAL ELEMENTS IN 2THIS CAN BE ENTERED SEPARATELY LATER. THIS WOULD BE USED'/' FOR LA 3RGE MATRICES FOR EXAMPLE WHERE IT IS EASIER TO TYPE INDIVIDUAL ELE 4MENTS'/' OR WHERE HETEROATOM VALUES ARE TO BE ENTERED.') READ (*,*,END=228,ERR=26) (C(J),J=1,NI) GO TO 20 26 WRITE(*,1055) 228 CONTINUE C 228 CLOSE(5) C OPEN(5,FILE='INPUT') IPRINT= 'YES' 20 WRITE (*,1069) 1069 FORMAT(' PRINT THE HUCKEL MATRIX (=N).') READ(*,9932,END=111,ERR=111) IPRINT WRITE(*,'(A3)') IPRINT GOTO 112 111 IPRINT='Y' C CLOSE(5) C OPEN(5,FILE='INPUT') 112 IF(IPRINT.EQ.'Y' .OR. IPRINT.EQ.'YES') WRITE (*,1060) (ITITLE(J),J 1=1,2) IF(IPRINT.EQ.'y' .OR. IPRINT.EQ.'yes') WRITE (*,1060) 1 (ITITLE(J),J=1,2) IF(IPRINT.EQ.'Y' .OR. IPRINT.EQ.'YES') CALL VECPRT (C,210,IN) IF(IPRINT.EQ.'y'.OR.IPRINT.EQ.'yes')CALL VECPRT (C,210,IN) C CONVERT VECTOR TO MATRIX IX1 = 0 DO 101 I1 = 1,IN DO 101 J1 = 1,I1 IX1 = IX1 + 1 AH(I1,J1) = C(IX1) AH(J1,I1) = C(IX1) 101 CONTINUE C C CHANGE MATRIX ELEMENTS IF DESIRED 230 CONTINUE C 230 CLOSE(5) C OPEN(5,FILE='INPUT') WRITE(*,*) IPRINT WRITE(*,201) 201 FORMAT(' CHANGE ANY ELEMENTS OF THE MATRIX (Y/N,=N)') IF(IPFLAG) WRITE(*,209) 209 FORMAT(' THIS OPTION IS USED IF YOU MADE AN ERROR ABOVE, OR IF YO 1U WANT TO CHANGE JUST ONE'/' MATRIX ELEMENT IN THE PREVIOUS MATRI 2X, IE TO INSERT DIFFERENT HETEROATOM VALUES ETC.') ICHANG = 'N' READ (*,9932,END=230,ERR=230) ICHANG IF(ICHANG.EQ.'YES' .OR. ICHANG.EQ.'Y') GOTO 223 IF(ICHANG.EQ.'yes' .OR. ICHANG.EQ.'y') GOTO 223 GOTO 224 223 WRITE(*,225) 225 FORMAT(' TYPE ROW AND COLUMN NUMBERS AND NEW VALUE OF MATRIX ELEME 1NT'/' (CTRL+d WHEN FINISHED).') IF(IPFLAG) WRITE(*,332) 332 FORMAT(' SEPARATING EACH BY AT LEAST A COMMA OR A SPACE, IE'/ 1'3,2,0.5'/' WHEN YOU HAVE MADE ALL THE CHANGES YOU WANT, ENTER A C 1TRL+Z IN RESPONSE TO THIS QUESTION.') READ(*,*,END=226) I2,J2, AH(I2,J2) AH(J2,I2) = AH(I2,J2) GOTO 223 226 CONTINUE C 226 CLOSE(5) C OPEN(5,FILE='INPUT') IX2 = 0 DO 227 I3 = 1,IN DO 227 J3 = 1,I3 IX2 = IX2 + 1 C(IX2) = AH(I3,J3) 227 CONTINUE GOTO 228 224 CONTINUE C SINCE THE RESONANCE INTEGRAL BETA IS A -VE QUANTITY, C MULTIPLY ALL MATRIX ELEMENTS BY - TO ORDER THE C EIGENVALUES CORRECTLY, IE MOST -VE ARE BONDING. C DO 30 KK = 1,NI C(KK) = -C(KK) 30 CONTINUE C WRITE OUT THE SECULAR MATRIX. C C DIAGONALIZE THE SECULAR MATRIX. C CALL SECOND (TIME1) CALL RSP (50,IN,NI,C,ROOT,1,VECT,FV1,FV2,IERR) C CALL SECOND (TIME2) DO 31 KK = 1,NI C(KK) = - C(KK) 31 CONTINUE TIME= 0.0 C TIME = TIME2-TIME1 WRITE (*,1120) TIME C WRITE OUT THE EIGENVALUES AND EIGENVECTORS. CALL MATOUT (VECT,ROOT,IN,IN,50) C WRITE OUT THE NUMBER OF DOUBLY AND SINGLY OCCUPIED ORBITALS WRITE (*,1070) NOCCO,NCLOSE C C TEST FOR DEGENERACY OF HOMO/LUMO. C IF (NOCCO.EQ.IN) GO TO 35 IF (ABS(ROOT(NOCCO)-ROOT(NOCCO+1)).LT.1.E-04) DEGEN = .TRUE. IF (DEGEN) WRITE (*,1080) C C NOW COMPUTE THE BOND ORDER MATRIX C LOWER TRIANGLE ONLY, SINCE THE MATRIX IS SYMMETRIC. 35 IE = 0 DO 90 J = 1,IN DO 80 I = 1,J IE = IE+1 AA = 0. C SUM OVER THE DOUBLY OCCUPIED ORBITALS C IF (NCLOSE.EQ.0) GO TO 50 DO 40 K = 1,NCLOSE IF (DEGEN.AND.K.EQ.NCLOSE.AND.NCLOSE.EQ.NOCCO) GO TO 40 AA = AA+VECT(I,K)*VECT(J,K) 40 CONTINUE C C IF HOMO/LUMO DEGENERATE, TAKE LINEAR COMBINATION OF BOTH. C IF (DEGEN.AND.NCLOSE.EQ.NOCCO) AA = AA+0.5*(VECT(I,NCLOSE)*VEC 1 T(J,NCLOSE)+VECT(I,NCLOSE+1)*VECT(J,NCLOSE+1)) AA = 2.*AA 50 IF (IEL.EQ.IOD.AND.IEXCIT.EQ.0) GO TO 70 C SUM OVER THE SINGLY OCCUPIED ORBITALS DO 60 K = IODD,NOCCO IF (DEGEN.AND.K.EQ.NOCCO) GO TO 60 AA = AA+VECT(I,K)*VECT(J,K) 60 CONTINUE C C IF HOMO/LUMO DEGENERATE, TAKE LINEAR COMBINATION OF BOTH. C IF (DEGEN) AA = AA+0.5*(VECT(I,NOCCO)*VECT(J,NOCCO)+VECT(I,NOC 1 CO+1)*VECT(J,NOCCO+1)) 70 CONTINUE P(IE) = AA 80 CONTINUE 90 CONTINUE C C FIND PI BOND ENERGY. C PIENA = 0. PIEN = 0. IF (NCLOSE.EQ.0) GO TO 110 DO 100 K = 1,NCLOSE PIEN = PIEN-2.*ROOT(K) PIENA = PIENA+2. 100 CONTINUE 110 IF (IEL.EQ.IOD.AND.IEXCIT.EQ.0) GO TO 130 DO 120 J = IODD,NOCCO PIENA = PIENA+1. 120 PIEN = PIEN-ROOT(J) 130 CONTINUE C C WRITE OUT THE PI ELECTRON ENERGY. C WRITE (*,1090) PIENA,PIEN C C C WRITE OUT THE DENSITY (BOND ORDER) MATRIX. C WRITE (*,1100) CALL VECPRT (P,210,IN) WRITE (*,1110) GO TO 10 140 STOP C C C 1000 FORMAT (' READ A TITLE (20 CHARACTERS) OR PRESS CTRL+d TO EXIT.') 1010 FORMAT (' READ ORDER OF SECULAR MATRIX (IE NUMBER OF ATOMS EXCLUDI 1NG HYDROGEN)'/ 1' AND NUMBER OF ELECTRONS CONJUGATED IN PI SYSTEM.') 1012 FORMAT (' READ A NEW HUCKEL MATRIX (= Y)') 1015 FORMAT (2A10) 1030 FORMAT (' TYPE 0 IF GROUND ELECTRONIC STATE OR 1 IF 1ST EXCITED ST 1ATE.') 1040 FORMAT (' READ BOTTOM TRIANGLE OF YOUR HUCKEL MATRIX, A ROW AT A T 1IME,'/' SEPARATING EACH VALUE BY AT LEAST A SPACE OR A COMMA. A TO 3TAL OF ',I3,' NUMBERS MUST BE TYPED'/' THEY CAN ALL BE ENTERED ON 2 ONE LINE, IE BUTADIENE:'/' 0 1 0 0 1 0 0 0 1 0') 1050 FORMAT (' YOU HAVE NOT READ IN SUFFICIENT DATA, START AGAIN.') 1055 FORMAT (' YOU HAVE NOT TYPED NUMBERS --- A MATRIX OF ZEROS HAS BEE 1N SET UP') 1060 FORMAT (/,' THE HMO SECULAR MATRIX FOR ',2A10,' IS:'/) 1070 FORMAT (/,' NUMBER OF OCCUPIED ORBITALS ',I3/'NUMBER OF CLOSED SHE 1LL ORBITALS ',I3/) 1080 FORMAT (/' HOMO/LUMO ASSUMED DEGENERATE'/) 1090 FORMAT (/' PI ELECTRON ENERGY ',F4.1,' ALPHA + ',F7.4,' BETA'/) 1100 FORMAT (/,' BOND ORDER (DENSITY) MATRIX: '/) 1110 FORMAT (///) 1120 FORMAT (/' TIME FOR DIAGONALIZATION',5X,F10.5,' SECONDS'//'EIGENVA 1LUES AND EIGENVECTORS: '//) END SUBROUTINE RSP (NM,N,NV,A,W,MATZ,Z,FV1,FV2,IERR) C INTEGER I,J,N,NM,NV,IERR,MATZ REAL A(NV), W(N), Z(NM,N), FV1(N), FV2(N) C C THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF C SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) C TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) C OF A REAL SYMMETRIC PACKED MATRIX. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRIX A, C C NV IS AN INTEGER VARIABLE SET EQUAL TO THE C DIMENSION OF THE ARRAY A AS SPECIFIED FOR C A IN THE CALLING PROGRAM. NV MUST NOT BE C LESS THAN N*(N+1)/2, C C A CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC C PACKED MATRIX STORED ROW-WISE, C C MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF C ONLY EIGENVALUES ARE DESIRED, OTHERWISE IT IS SET TO C ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. C C ON OUTPUT- C C W CONTAINS THE EIGENVALUES IN ASCENDING ORDER, C C Z CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO, C C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN C ERROR COMPLETION CODE DESCRIBED IN SECTION 2B OF THE C DOCUMENTATION. THE NORMAL COMPLETION CODE IS ZERO, C C FV1 AND FV2 ARE TEMPORARY STORAGE ARRAYS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ----------------------------------------------------------------- C IF (N.LE.NM) GO TO 10 IERR = 10*N GO TO 60 10 IF (NV.GE.(N*(N+1))/2) GO TO 20 IERR = 20*N GO TO 60 C 20 CALL TRED3 (N,NV,A,W,FV1,FV2) IF (MATZ.NE.0) GO TO 30 C ********** FIND EIGENVALUES ONLY ********** C CALL TQLRAT (N,W,FV2,IERR) GO TO 60 C ********** FIND BOTH EIGENVALUES AND EIGENVECTORS ********** 30 DO 50 I = 1,N C DO 40 J = 1,N Z(J,I) = 0.0 40 CONTINUE C Z(I,I) = 1.0 50 CONTINUE C CALL TQL2 (NM,N,W,FV1,Z,IERR) IF (IERR.NE.0) GO TO 60 CALL TRBAK3 (NM,N,NV,A,N,Z) 60 RETURN C ********** LAST CARD OF RSP ********** END C C ----------------------------------------------------------------- C C C ----------------------------------------------------------------- C SUBROUTINE TQL2 (NM,N,D,E,Z,IERR) C INTEGER I,J,K,L,M,N,II,L1,NM,MML,IERR REAL D(N), E(N), Z(NM,N) DIMENSION CC(50), SS(50) REAL B,C,F,G,H,P,R,S,MACHEP C REAL SQRT,ABS,SIGN C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2, C NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND C WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS C OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD. C THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO C BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS C FULL MATRIX TO TRIDIAGONAL FORM. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRIX, C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX, C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY, C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE C REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS C OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN C THE IDENTITY MATRIX. C C ON OUTPUT- C C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT C UNORDERED FOR INDICES 1,2,...,IERR-1, C C E HAS BEEN DESTROYED, C C Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC C TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE, C Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED C EIGENVALUES, C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ----------------------------------------------------------------- C C ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC. C C ********** MACHEP = 2.**(-47) C IERR = 0 IF (N.EQ.1) GO TO 170 C DO 10 I = 2,N 10 E(I-1) = E(I) C F = 0.0 B = 0.0 E(N) = 0.0 C DO 120 L = 1,N J = 0 H = MACHEP*(ABS(D(L))+ABS(E(L))) IF (B.LT.H) B = H C ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT ********** DO 20 M = L,N IF (ABS(E(M)).LE.B) GO TO 30 C ********** E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP ********** 20 CONTINUE C 30 IF (M.EQ.L) GO TO 110 40 IF (J.EQ.30) GO TO 160 J = J+1 C ********** FORM SHIFT ********** L1 = L+1 G = D(L) P = (D(L1)-G)/(2.0*E(L)) R = SQRT(P*P+1.0) D(L) = E(L)/(P+SIGN(R,P)) H = G-D(L) C DO 50 I = L1,N 50 D(I) = D(I)-H C F = F+H C ********** QL TRANSFORMATION ********** P = D(M) C = 1.0 S = 0.0 MML = M-L C ********** FOR I=M-1 STEP -1 UNTIL L DO -- ********** C MOD 2 C MODIFICATION BY K.Z. C OLD CODE READ DO 200 II = 1, MML DO 80 II = 1,MML I = M-II G = C*E(I) H = C*P IF (ABS(P).LT.ABS(E(I))) GO TO 60 C = E(I)/P R = SQRT(C*C+1.0) E(I+1) = S*P*R S = C/R C = 1.0/R GO TO 70 60 C = P/E(I) R = SQRT(C*C+1.0) E(I+1) = S*E(I)*R S = 1.0/R C = C*S 70 P = C*D(I)-S*G D(I+1) = H+S*(C*G+S*D(I)) C CODING INSERTED BY K.Z. C CC(I) = C SS(I) = S 80 CONTINUE C ********** FORM VECTOR ********** C MODIFICATION BY K.Z. C OLD CODE: C DO 180 K = 1,N C H = Z(K,I+1) C Z(K,I+1) = S * Z(K,I) + C*H C Z(K,I) = C * Z(K,I) - S*H C 180 CONTINUE C DO 100 K = 1,N T1 = Z(K,M) DO 90 II = 1,MML I = M-II T2 = Z(K,I) Z(K,I+1) = SS(I)*T2+CC(I)*T1 T1 = CC(I)*T2-SS(I)*T1 90 CONTINUE Z(K,L) = T1 100 CONTINUE C E(L) = S*P D(L) = C*P IF (ABS(E(L)).GT.B) GO TO 40 110 D(L) = D(L)+F 120 CONTINUE C ********** ORDER EIGENVALUES AND EIGENVECTORS ********** DO 150 II = 2,N I = II-1 K = I P = D(I) C DO 130 J = II,N IF (D(J).GE.P) GO TO 130 K = J P = D(J) 130 CONTINUE C IF (K.EQ.I) GO TO 150 D(K) = D(I) D(I) = P C DO 140 J = 1,N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 140 CONTINUE C 150 CONTINUE C GO TO 170 C ********** SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS ********** 160 IERR = L 170 RETURN C ********** LAST CARD OF TQL2 ********** END C C ----------------------------------------------------------------- C SUBROUTINE TRBAK3 (NM,N,NV,A,M,Z) C INTEGER I,J,K,L,M,N,IK,IZ,NM,NV REAL A(NV), Z(NM,M) REAL H,S C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK3, C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING C SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY TRED3. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRIX, C C NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A C AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT, C C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANSFORMATIONS C USED IN THE REDUCTION BY TRED3 IN ITS FIRST C N*(N+1)/2 POSITIONS, C C M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED, C C Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED C IN ITS FIRST M COLUMNS. C C ON OUTPUT- C C Z CONTAINS THE TRANSFORMED EIGENVECTORS C IN ITS FIRST M COLUMNS. C C NOTE THAT TRBAK3 PRESERVES VECTOR EUCLIDEAN NORMS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ----------------------------------------------------------------- C IF (M.EQ.0) GO TO 50 IF (N.EQ.1) GO TO 50 C DO 40 I = 2,N L = I-1 IZ = (I*L)/2 IK = IZ+I H = A(IK) IF (H.EQ.0.0) GO TO 40 C DO 30 J = 1,M S = 0.0 IK = IZ C DO 10 K = 1,L IK = IK+1 S = S+A(IK)*Z(K,J) 10 CONTINUE C ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW ********** S = (S/H)/H IK = IZ C DO 20 K = 1,L IK = IK+1 Z(K,J) = Z(K,J)-S*A(IK) 20 CONTINUE C 30 CONTINUE C 40 CONTINUE C 50 RETURN C ********** LAST CARD OF TRBAK3 ********** END C C ----------------------------------------------------------------- SUBROUTINE TRED3 (N,NV,A,D,E,E2) C INTEGER I,J,K,L,N,II,IZ,JK,NV REAL A(NV), D(N), E(N), E2(N) REAL F,G,H,HH,SCALE C REAL SQRT,ABS,SIGN C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3, C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX, STORED AS C A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX C USING ORTHOGONAL SIMILARITY TRANSFORMATIONS. C C ON INPUT- C C N IS THE ORDER OF THE MATRIX, C C NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A C AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT, C C A CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC C INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL C ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS. C C ON OUTPUT- C C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL C TRANSFORMATIONS USED IN THE REDUCTION, C C D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX, C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO, C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ----------------------------------------------------------------- C C ********** FOR I=N STEP -1 UNTIL 1 DO -- ********** DO 110 II = 1,N I = N+1-II L = I-1 IZ = (I*L)/2 H = 0.0 SCALE = 0.0 IF (L.LT.1) GO TO 20 C ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) ********** DO 10 K = 1,L IZ = IZ+1 D(K) = A(IZ) SCALE = SCALE+ABS(D(K)) 10 CONTINUE C IF (SCALE.NE.0.0) GO TO 30 20 E(I) = 0.0 E2(I) = 0.0 GO TO 100 C 30 DO 40 K = 1,L D(K) = D(K)/SCALE H = H+D(K)*D(K) 40 CONTINUE C E2(I) = SCALE*SCALE*H F = D(L) G = -SIGN(SQRT(H),F) E(I) = SCALE*G H = H-F*G D(L) = F-G A(IZ) = SCALE*D(L) IF (L.EQ.1) GO TO 100 F = 0.0 C DO 80 J = 1,L G = 0.0 JK = (J*(J-1))/2 C ********** FORM ELEMENT OF A*U ********** C C THE ORIGINAL FORTRAN CODE AT THIS POINT READS C DO 180 K = 1, L C JK = JK + 1 C IF (K .GT. J) JK = JK + K - 2 C G = G + A(JK) * D(K) C 180 CONTINUE C C A BETTER (AND EQUIVALENT) CODING IS C C DO 180 K=1,J C JK=JK+1 C G=G+A(JK)*D(K) C 180 CONTINUE C IF(J.EQ.L) GO TO 182 C J1= J+1 C DO 181 K=J1,L C JK= JK+K-1 C G= G+A(JK)*D(K) C 181 CONTINUE C 182 CONTINUE C C END OF NEW CODE C C THE ABOVE CODE WAS DUE TO ROBERT K. DEWAR C THE CODE THAT FOLLOW IS DUE TO K.ZINK K = 0 50 K = K+1 JK = JK+1 G = G+A(JK)*D(K) IF (K.LT.J) GO TO 50 IF (K.EQ.L) GO TO 70 60 JK = JK+K K = K+1 G = G+A(JK)*D(K) IF (K.LT.L) GO TO 60 70 CONTINUE C END OF NEW CODE DUE TO K.Z. C C ********** FORM ELEMENT OF P ********** E(J) = G/H F = F+E(J)*D(J) 80 CONTINUE C HH = F/(H+H) JK = 0 C ********** FORM REDUCED A ********** DO 90 J = 1,L F = D(J) G = E(J)-HH*F E(J) = G C DO 90 K = 1,J JK = JK+1 A(JK) = A(JK)-F*E(K)-G*D(K) 90 CONTINUE C 100 D(I) = A(IZ+1) A(IZ+1) = SCALE*SQRT(H) 110 CONTINUE C RETURN C ********** LAST CARD OF TRED3 ********** END C SUBROUTINE MATOUT (A,B,NC,NR,NDIM) DIMENSION A(NDIM,NC), B(NC),BX(50) CHARACTER LINE(21)*6 C DO 555 I = 1,21 555 LINE(I) = '------' C *** WISH TO OUTPUT THE MATRIX A(NR,NC) C DO 666 NN = 1,NC 666 BX(NN) = - B(NN) KA = 1 KC = 6 10 KB = MIN0(KC,NC) WRITE (*,1010) (BX(I),I=KA,KB) WRITE (*,1000) LA = 1 LC = 48 20 LB = MIN0(LC,NR) N = 0 DO 30 I = LA,LB WRITE (*,1030) I,(A(I,J),J=KA,KB) N = N+1 IF (N.LT.6) GO TO 30 WRITE (*,1020) N = 0 30 CONTINUE IF (LB.EQ.NR) GO TO 40 LA = LC+1 LC = LC+48 WRITE (*,1040) GO TO 20 40 IF (KB.EQ.NC) RETURN KA = KC+1 KC = KC+6 IF (NR.GT.30) WRITE (*,1040) GO TO 10 C C C C 1000 FORMAT (/3X,'LCAO COEFFICIENT') 1010 FORMAT (/3X,22H ORBITAL ENERGIES,5H A+ ,6(F6.3,2HB;,' A+ ')) 1020 FORMAT (2H ) 1030 FORMAT (20H PZ ORBITAL ON ATOM ,I5,6F12.5) 1040 FORMAT (1H1) END SUBROUTINE VECPRT (A,ISIZE,NUMB) DIMENSION A(ISIZE) CHARACTER LINE(21)*6 DO 10 I = 1,21 10 LINE(I) = '------' LIMIT = (NUMB*(NUMB+1))/2 KK = 8 NA = 1 20 LL = 0 M = MIN0((NUMB+1-NA),5) MA = 2*M+1 M = NA+M-1 WRITE (*,1000) (I,I=NA,M) WRITE (*,1010) (LINE(K),K=1,MA) DO 40 I = NA,NUMB LL = LL+1 K = (I*(I-1))/2 L = MIN0((K+M),(K+I)) K = K+NA IF ((KK+LL).LE.50) GO TO 30 WRITE (*,1020) WRITE (*,1000) (N,N=NA,M) WRITE (*,1010) (LINE(N),N=1,MA) KK = 4 LL = 0 30 WRITE (*,1030) I,(A(N),N=K,L) 40 CONTINUE IF (L.GE.LIMIT) GO TO 50 KK = KK+LL+4 NA = M+1 IF ((KK+NUMB+1-NA).LE.50) GO TO 20 KK = 4 WRITE (*,1020) GO TO 20 50 RETURN C C C 1000 FORMAT (1H0/8X,7HATOM ,5(I5,7X)) 1010 FORMAT (1H ,6X,21A6) 1020 FORMAT (1H1) 1030 FORMAT (1H ,6H ATOM ,I3,2X,10F12.6) END