!
!     HUCKEL MO METHOD
!          CODED BY T. ISHIDA
!          REFERENCE) "BUNSHIKIDOUHOU," O. KIKUCHI, 1971.
!
      PROGRAM HUCKEL
!
      PARAMETER (NN=100,NC=12)
      DIMENSION H(NN,NN),C(NN,NN),OC(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.txt')
      OPEN (11,FILE='HUCKEL.OUT.txt',STATUS='REPLACE')
      WRITE (*,100)
      WRITE (11,100)
  100 FORMAT (/17H HUCKEL MO METHOD/) 
      CALL FILEINP(H,H0,N,NELEC,NN)
      CALL DIAG(H,C,N,NN)
      CALL REARR(H,C,N,NN)
      CALL DEGEN(H,N,NELEC,OC,NN,ND,NH)
      CALL ENECOE(H,C,N,OC,NN,NC)
      CALL BORDER(C,N,NELEC,OC,PRS,NN,NC,ND,NH)
      CALL RINDEX(C,H,N,NELEC,FRE,FRN,FRR,NN,ND,NH)
      CALL MCLAC(C,H,N,NELEC,OC,PRS,SPINH,SPINM,P,NN,ND,NH)
      WRITE (*,110)
      WRITE (11,110)
  110 FORMAT (//19H END OF CALCULATION/)
      CLOSE (10)
      CLOSE (11)
      STOP
      END PROGRAM HUCKEL
!
!     N:     number of AOs (MOs)
!     NELEC: number of es
!     NN:    matrix size limit
!     H:     Huckel matrix
!     E:     eigenvalue
!     C:     coefficient matrix
!     ND:    degeneracy around HOMO
!     NH:    next HOMO#
!     OC:    occupancy
!     loop threshold = 0.000000001 (in DIAG)
!     degeneracy criterion = 0.001 (in DEGEN,RINDEX)
!
!
      SUBROUTINE FILEINP(H,H0,N,NELEC,NN)
!
      DIMENSION H(NN,NN),H0(NN,NN),TITLE(12)
      READ (10,100) TITLE  
  100 FORMAT (12A4)     
      READ (10,*) N,NELEC
      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 (*,100) TITLE
      WRITE (11,100) TITLE
      WRITE (*,200)
      WRITE (11,200)
  200 FORMAT (/11H INPUT DATA)
      WRITE (*,210) N,NELEC
      WRITE (11,210) N,NELEC
  210 FORMAT (2I8)
      DO 20 I=1,N
           WRITE (*,220) (H(I,J),J=1,I)
           WRITE (11,220) (H(I,J),J=1,I)
  220      FORMAT (30F5.2/)
   20 CONTINUE
      WRITE (*,230)
      WRITE (11,230)
  230 FORMAT (1H )
      RETURN
      END
!
!
      SUBROUTINE DIAG(A,C,N,NN)
!     JACOBI METHOD
!
      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.0.000000001) 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
!                WRITE (*,*) '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
!                WRITE (*,*) 'COS ',COS,'    SIN ',SIN
   90 CONTINUE
      GO TO 20
      END
!
!
      SUBROUTINE REARR(E,C,N,NN)
!
      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
!
!
      SUBROUTINE DEGEN(E,N,NELEC,OC,NN,ND,NH)
!
      DIMENSION E(NN,NN),OC(NN)
      NOCC=(NELEC+1)/2
      ND=0
      NH=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
   10 CONTINUE
!      WRITE (*,*) 'DEGENERACY',ND,' NextHOMO#',NH
      DO 20 I=1,N
           OC(I)=0.
   20 CONTINUE
      DO 30 I=1,NH
           OC(I)=2.
   30 CONTINUE
      DO 40 I=NH+1,NH+ND  
           OC(I)=FLOAT(NELEC-NH*2)/FLOAT(ND)
   40 CONTINUE
      RETURN
      END
!
!
      SUBROUTINE ENECOE(E,C,N,OC,NN,NC)
!
      DIMENSION E(NN,NN),C(NN,NN),OC(NN)
      WRITE (*,100)
      WRITE (11,100)
  100 FORMAT (/53H ORBITAL ENERGIES, OCCUPANCIES, 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 (*,110) (E(K,K), K=I1,I2)
           WRITE (11,110) (E(K,K), K=I1,I2)
           WRITE (*,110) (OC(K), K=I1,I2)
           WRITE (11,110) (OC(K), K=I1,I2)
  110      FORMAT (/1H ,30F8.4)
           WRITE (*,120)
           WRITE (11,120)
  120      FORMAT (1H )
           DO 10 J=1,N
                WRITE (*,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
!
!
      SUBROUTINE BORDER(C,N,NELEC,OC,PRS,NN,NC,ND,NH)
!
      DIMENSION C(NN,NN),PRS(NN,NN),OC(NN)
      DO 10 I=1,N
           DO 10 J=1,N
                PRS(I,J)=0.
                DO 20 K=1,N
   20                PRS(I,J)=PRS(I,J)+OC(K)*C(I,K)*C(J,K)
   10 PRS(J,I)=PRS(I,J)
      WRITE (*,100)
      WRITE (11,100)
  100 FORMAT (/1H ,/35H Pi BOND ORDER AND ELECTRON DENSITY)
      INDEX=(N-1)/NC+1
      DO 30 I=1,INDEX
           WRITE (*,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 (*,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
!
!
      SUBROUTINE RINDEX(C,E,N,NELEC,FRE,FRN,FRR,NN,ND,NH)
!
      DIMENSION C(NN,NN),E(NN,NN),FRE(NN),FRN(NN),FRR(NN)
      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 (*,100)
      WRITE (11,100)
  100 FORMAT (//17H REACTION INDEXES //'     FR(E)    FR(N)    FR(R)  ')
      WRITE (*,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
!
!
      SUBROUTINE MCLAC(C,E,N,NELEC,OC,PRS,SPINH,SPINM,P,NN,ND,NH)
!
      DIMENSION E(NN,NN),C(NN,NN),OC(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 (*,100)
      WRITE (11,100)
  100 FORMAT (//13H    ELECTRON-,20H        SPIN-DENSITY/
&       12H     DENSITY,10H       HMO,15H      MCLACHLAN/)
      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 (*,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 
           DO 60 I=1,N
                SPINH(I)=0.
                DO 60 J=1,ND
   60                SPINH(I)=SPINH(I)+OC(NH+J)*C(I,NH+J)**2
           WRITE (*,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
      RETURN
      END
!
