C C HUCKEL MO METHOD C CODED BY ISHIDA C REFERENCE) "BUNSHIKIDOUHOU," O.KIKUCHI, 1971. C PROGRAM HUCKEL C PARAMETER (NN=100,NC=12) DIMENSION H(NN,NN),C(NN,NN),PRS(NN,NN),H0(NN,NN), + FRE(NN),FRN(NN),FRR(NN),SPINH(NN),SPINM(NN), + P(NN,NN) OPEN (10,FILE='HUCKEL.INP',STATUS='OLD') OPEN (11,FILE='HUCKEL.OUT',STATUS='NEW') WRITE (9,100) WRITE (11,100) 100 FORMAT (/17H HUCKEL MO METHOD/) CALL FILEINP(H,H0,N,NELEC,D,NN) CALL DIAG(H,C,N,D,NN) CALL REARR(H,C,N,NN) CALL ENECOE(H,C,N,NN,NC) CALL DEGEN(H,N,NELEC,NN,ND,NH) CALL BORDER(C,N,NELEC,PRS,NN,NC,ND,NH) CALL RINDEX(C,H,N,NELEC,FRE,FRN,FRR,NN,ND,NH) CALL MCLAC(C,H,N,NELEC,PRS,SPINH,SPINM,P,NN,ND,NH) C CALL OMEGA1(.....) WRITE (9,110) WRITE (11,110) 110 FORMAT (//19H END OF CALCULATION/) CLOSE (10) CLOSE (11) PAUSE STOP END C C SUBROUTINE FILEINP(H,H0,N,NELEC,D,NN) C DIMENSION H(NN,NN),H0(NN,NN),TITLE(12) READ (10,100) TITLE 100 FORMAT (12A4) READ (10,*) N,D IF (D.GE.1.0) THEN NELEC=INT(D) READ (10,*) D ELSE NELEC=N END IF READ (10,*) ((H(I,J),J=1,I),I=1,N) DO 10 I=1,N DO 10 J=1,I H(J,I)=H(I,J) H0(I,J)=H(I,J) H0(J,I)=H(I,J) 10 CONTINUE WRITE (9,100) TITLE WRITE (11,100) TITLE WRITE (9,200) WRITE (11,200) 200 FORMAT (/11H INPUT DATA) WRITE (9,210) N,NELEC,D WRITE (11,210) N,NELEC,D 210 FORMAT (/2I8/1F13.10/) DO 20 I=1,N WRITE (9,220) (H(I,J),J=1,I) WRITE (11,220) (H(I,J),J=1,I) 220 FORMAT (30F5.2/) 20 CONTINUE WRITE (9,230) WRITE (11,230) 230 FORMAT (1H ,/) RETURN END C C SUBROUTINE DIAG(A,C,N,D,NN) C JACOBI METHOD C DIMENSION A(NN,NN),C(NN,NN) DO 10 I=1,N DO 10 J=1,N C(I,J)=0. C(I,I)=1. 10 CONTINUE 20 AMAX=ABS(A(2,1)) DO 30 I=2,N I1=I-1 DO 30 J=1,I1 IF (ABS(A(I,J)).GT.AMAX) AMAX=ABS(A(I,J)) 30 CONTINUE IF (AMAX.LE.D) RETURN SML=0.1*AMAX DO 90 J=2,N J1=J-1 DO 90 K=1,J1 IF (ABS(A(J,K)).LE.SML) GO TO 90 C WRITE (9,*) 'PIVOT ',A(J,K) AKK=A(K,K) AJJ=A(J,J) TANB=(AJJ-AKK)/(2.*A(K,J)) IF (TANB.LE.0.) THEN TAN=TANB-SQRT(1.+TANB*TANB) ELSE TAN=TANB+SQRT(1.+TANB*TANB) END IF COS=1./SQRT(1.+TAN*TAN) SIN=COS*TAN DO 80 L=1,N IF (L.EQ.K) THEN A(K,K)=AKK*COS**2+A(J,K)*2.*SIN*COS + +AJJ*SIN**2 A(J,J)=AKK+AJJ-A(K,K) ELSE IF (L.NE.J) THEN AKL=A(K,L) AJL=A(J,L) A(K,L)=AKL*COS+AJL*SIN A(J,L)=AJL*COS-AKL*SIN A(L,K)=A(K,L) A(L,J)=A(J,L) END IF END IF 80 CONTINUE A(J,K)=A(J,K)*(COS*COS-SIN*SIN)+(AJJ-AKK)*SIN*COS A(K,J)=A(J,K) DO 85 I=1,N CIJ=C(I,J) CIK=C(I,K) C(I,J)=CIJ*COS-CIK*SIN C(I,K)=CIJ*SIN+CIK*COS 85 CONTINUE C WRITE (9,*) 'COS ',COS,' SIN ',SIN 90 CONTINUE C CALL PRTMTX(A,N,NN,NC) C CALL PRTMTX(C,N,NN,NC) C PAUSE GO TO 20 END C C SUBROUTINE REARR(E,C,N,NN) C DIMENSION E(NN,NN),C(NN,NN) DO 10 I=1,N DO 10 K=I,N IF (E(K,K).GT.E(I,I)) THEN X=E(I,I) E(I,I)=E(K,K) E(K,K)=X DO 20 J=1,N X=C(J,I) C(J,I)=C(J,K) C(J,K)=X 20 CONTINUE END IF 10 CONTINUE RETURN END C C SUBROUTINE ENECOE(E,C,N,NN,NC) C DIMENSION E(NN,NN),C(NN,NN) WRITE (9,100) WRITE (11,100) 100 FORMAT (/39H ORBITAL ENERGIES AND LCAO COEFFICIENTS) INDEX=(N-1)/NC+1 DO 10 I=1,INDEX I1=NC*(I-1)+1 I2=NC*I IF (N.LT.I2) I2=N WRITE (9,110) (E(K,K), K=I1,I2) WRITE (11,110) (E(K,K), K=I1,I2) 110 FORMAT (//1H ,30F8.4) WRITE (9,120) WRITE (11,120) 120 FORMAT (1H ) DO 10 J=1,N WRITE (9,130) (C(J,K), K=I1,I2) WRITE (11,130) (C(J,K), K=I1,I2) 130 FORMAT (1H ,30F8.4) 10 CONTINUE RETURN END C C SUBROUTINE PRTMTX(A,N,NN,NC) C DIMENSION A(NN,NN) WRITE (9,200) 200 FORMAT (1H ,/16H MATRIX ELEMENTS/) INDEX=(N-1)/NC+1 DO 10 I=1,INDEX I1=NC*(I-1)+1 I2=NC*I IF (N.LT.I2) I2=N DO 10 J=1,N WRITE (9,100) (A(J,K),K=I1,I2) 100 FORMAT (30F8.3) 10 CONTINUE RETURN END C C SUBROUTINE BORDER(C,N,NELEC,PRS,NN,NC,ND,NH) C DIMENSION C(NN,NN),PRS(NN,NN) C NOCC=N/2 FACTOR=FLOAT((NELEC-NH*2))/FLOAT(ND) C WRITE (9,*) FACTOR DO 10 I=1,N DO 10 J=1,N PRS(I,J)=0. DO 20 K=1,NH C DO 20 K=1,NOCC 20 PRS(I,J)=PRS(I,J)+2.*C(I,K)*C(J,K) DO 25 K=1,ND 25 PRS(I,J)=PRS(I,J)+FACTOR*C(I,NH+K)*C(J,NH+K) 10 PRS(J,I)=PRS(I,J) WRITE (9,100) WRITE (11,100) 100 FORMAT (//32H BOND ORDER AND ELECTRON DENSITY) INDEX=(N-1)/NC+1 DO 30 I=1,INDEX WRITE (9,110) WRITE (11,110) 110 FORMAT (1H /) I1=NC*(I-1)+1 I2=NC*I IF (N.LT.I2) I2=N DO 30 J=1,N WRITE (9,120) (PRS(J,K),K=I1,I2) WRITE (11,120) (PRS(J,K),K=I1,I2) 120 FORMAT (1H ,30F8.4) 30 CONTINUE RETURN END C C SUBROUTINE RINDEX(C,E,N,NELEC,FRE,FRN,FRR,NN,ND,NH) C DIMENSION C(NN,NN),E(NN,NN),FRE(NN),FRN(NN),FRR(NN) C WRITE (9,*) NELEC, ' DEG. ',ND,' NHOMO',NH IF((NELEC.EQ.NH*2).OR.(NELEC.EQ.((NH+ND)*2))) THEN NOCC=NELEC/2 NOCC1=NOCC+1 INDEX=0 DO 60 I=1,NOCC IF (ABS(E(NOCC,NOCC)-E(I,I)).LT. 0.001) INDEX=INDEX+1 60 CONTINUE DO 70 I=1,N FRE(I)=0. DO 65 J=1,INDEX FRE(I)=FRE(I)+2.*C(I,NOCC1-J)**2 65 CONTINUE FRE(I)=FRE(I)/FLOAT(INDEX) 70 CONTINUE INDEX=0 DO 80 I=NOCC1,N IF (ABS(E(NOCC1,NOCC1)-E(I,I)).LT. 0.001) INDEX=INDEX+1 80 CONTINUE DO 90 I=1,N FRN(I)=0. DO 85 J=1,INDEX FRN(I)=FRN(I)+2.*C(I,NOCC+J)**2 85 CONTINUE FRN(I)=FRN(I)/FLOAT(INDEX) FRR(I)=.5*(FRE(I)+FRN(I)) 90 CONTINUE WRITE (9,100) WRITE (11,100) 100 FORMAT (//17H REACTION INDEXES // + ' FR(E) FR(N) FR(R) ' /) WRITE (9,110) (FRE(I),FRN(I),FRR(I),I=1,N) WRITE (11,110) (FRE(I),FRN(I),FRR(I),I=1,N) 110 FORMAT (1H ,3F9.4) END IF RETURN END C C SUBROUTINE DEGEN(E,N,NELEC,NN,ND,NH) C DIMENSION E(NN,NN) NOCC=(NELEC+1)/2 C WRITE (9,*) 'NELEC',NELEC, ' NOCC',NOCC ND=0 DO 10 I=1,N DE=E(NOCC,NOCC)-E(I,I) IF (DE.LT.-0.001) NH=I IF (ABS(DE).LT.0.001) ND=ND+1 C WRITE (9,*) DE,NH 10 CONTINUE C WRITE (9,*) 'DEG.',ND,' NHOMO',NH RETURN END C C SUBROUTINE MCLAC(C,E,N,NELEC,PRS,SPINH,SPINM,P,NN,ND,NH) C DIMENSION E(NN,NN),C(NN,NN),PRS(NN,NN), + SPINH(NN),SPINM(NN),P(NN,NN) DATA RAMDA/1.2/ IF (NELEC.NE.NH*2).AND.(NELEC.NE.(NH+ND)*2) THEN WRITE (9,100) WRITE (11,100) 100 FORMAT (//4X9HELECTRON-,8X12HSPIN-DENSITY/ + 5X7HDENSITY,7X3HHMO,6X9HMCLACHLAN/) IF((NELEC.EQ.(NH*2+1)).AND.(ND.EQ.1)) THEN DO 10 I=1,N SPINH(I)=C(I,NH+1)**2 SPINM(I)=0. 10 CONTINUE NOCC=NH+1 DO 20 K=1,N DO 20 L=K,N P(K,L)=0. DO 30 I=1,NH DO 30 J=NOCC,N 30 P(K,L)=P(K,L)+4.*C(K,I)*C(K,J)*C(L,I)*C(L,J) + /(E(J,J)-E(I,I)) 20 P(L,K)=P(K,L) DO 40 I=1,N DO 50 J=1,N 50 SPINM(I)=SPINM(I)+P(I,J)*SPINH(J) 40 SPINM(I)=SPINH(I)-RAMDA*SPINM(I) WRITE (9,110) (PRS(I,I),SPINH(I),SPINM(I),I=1,N) WRITE (11,110) (PRS(I,I),SPINH(I),SPINM(I),I=1,N) 110 FORMAT (3F12.5) ELSE FACTOR=FLOAT(NELEC-NH*2)/FLOAT(ND) DO 60 I=1,N SPINH(I)=0. DO 60 J=1,ND 60 SPINH(I)=SPINH(I)+FACTOR*C(I,NH+J)**2 WRITE (9,120) (PRS(I,I),SPINH(I),I=1,N) WRITE (11,120) (PRS(I,I),SPINH(I),I=1,N) 120 FORMAT (2F12.5) END IF END IF C WRITE (9,*) 'STOP' RETURN END C C