C=385    INTERACTIVE HUCKEL PGM FOR THE DEC 10 
      PROGRAM HUCKEL     
C  
C  
C AN INTERACTIVE HUCKEL PROGRAM FOR THE DECSYSTEM 10      
C HUCKEL ENERGIES, MO COEFFICIENTS, BOND ORDER MATRICES,  
C NET CHARGES, FREE VALENCES, AND SPIN DENSITIES ARE CALC-D          
C THE GROUND STATE IS CALC-D FIRST, AND EXCITED STATES MAY BE        
C DONE AS WELL.          
C  
C HETERO ATOMS (B,N,O,F,SI,P,S,CL) MAY BE INCLUDED.       
C  
      PARAMETER NMAX=80  
      DIMENSION H(NMAX,NMAX),V(NMAX,NMAX),P(NMAX,NMAX),E(NMAX),JC(NMAX),         

     1NOC(NMAX),A(NMAX,NMAX),NTP(NMAX),ZC(NMAX),Q(NMAX),FV(NMAX),    
     2RHO(NMAX),FVO(13),NDEG(NMAX)  
      EQUIVALENCE (JC,NDEG)         
      COMMON/LABEL/ NAME(13)        
      DATA FVO/1.73205,1.70526,1.39262,1.58280,0.90934,0.94208,0.17914,          

     11.73205,1.40886,1.66578,0.96197,1.22930,0.32111/    
      DATA NAME/'C ','B ','N1','N2','O1','O2','F ','SI','P1','P2','S1',          

     1'S2','CL'/         
      IDH=NMAX
C  
C  
C THIS IS HOUSEKEEPING, CLEANING OUT THE ARRAYS.          
C  
C  
    1 DO 2 I=1,IDH       
      NTP(I)=0
      E(I)=0. 
      DO 2 J=1,IDH       
      A(I,J)=0.          
      H(I,J)=0.          
      V(I,J)=0.          
    2 P(I,J)=0.          
C  
C  
C THE PROGRAM INTRODUCES ITSELF AND QUERIES THE USER. FOR BEGINNERS  
C A BRIEF 'HELP' FILE MAY BE TYPED OUT UPON REQUEST.      
C  
C  
      WRITE(5,3)IDH      
    3 FORMAT(//5X,'HUCKEL MOLECULAR ORBITAL CALC-N, MAXIMUM SIZE = ',I4)         

    4 WRITE(5,5)         
    5 FORMAT(1H$,5X,'DO YOU WANT TO RUN>(TYPE Y OR N, H FOR HELP)    
     1[Y]    ')          
      READ(5,8)ANS       
      IF(ANS.EQ.'H')GO TO 67        
      IF (ANS.EQ.'N')GO TO 66       
C  
C  
C THIS SECTION CHOOSES THE METHOD OF CONSTRUCTION OF THE HUCKEL      
C MATRIX. *EASY* BUILDS THE SYSTEM BY CROSSLINKING A LINEAR OR       
C MONOCYCLIC POLYENE. *HARD* WORKS FROM A SERIES OF CONNECTED-ATOM   
C LISTS.      
C  
C  
    6 WRITE(5,7)         
    7 FORMAT(/1H$,5X,'BUILD FROM LINEAR OR MONOCYCLIC SYSTEM>        
     1(Y OR N)[Y] ')     
      READ(5,8)ANS       
    8 FORMAT(A1)         
      IF(ANS.EQ.'N')GO TO 10        
    9 CALL EASY(H,N,IDH) 
      GO TO 11
   10 CALL HARD(H,N,IDH) 
C  
C  
C THIS SECTION ASSIGNS ALL ATOMS AS CARBONS.   
C IF THERE ARE HETERO ATOMS, *HETERO* MAKES THE CHANGES REQUIRED     
C BASED ON INTERACTIVE INPUT.       
C  
C  
   11 WRITE(5,12)        
   12 FORMAT(1H$,5X,'ARE THERE ANY HETEROATOMS>(Y OR N)[N]   ')      
      READ(5,8)ANS       
      DO 13 I=1,N        
      NTP(I)=1
      ZC(I)=1.
   13 CONTINUE
      IF(ANS.NE.'Y')GO TO 15        
   14 CALL HETERO(H,NTP,ZC,N,IDH)   
   15 CONTINUE
C  
C  
C THIS IS AN INACTIVE SECTION FOR PRINTING OUT THE HUCKEL MATRIX     
C IT MAY BE REACTIVATED BY REMOVING THE LEADING 'C-S'.    
C  
C  
C  
C  
C     WRITE(5,7777)      
C7777 FORMAT(/10X,'HUCKEL MATRIX BY COLUMNS')  
C     CALL PPRNT(H,N,IDH)
C  
C  
   16 WRITE(5,17)        
   17 FORMAT(/1H$,5X,'PI IONIC CHARGE = [0]    ')         
      READ(5,18)NCHG     
   18 FORMAT(I)          
      NE=0    
      DO 19 I=1,N        
      NE=NE+ZC(I)        
   19 CONTINUE
      NE=NE-NCHG         
      WRITE(5,20)        
   20 FORMAT(1H$,5X,'CHECK DATA; ARE YOU READY TO RUN(Y OR N)>[Y]   ')           

      READ(5,8) ANS      
      IF(ANS.NE.'N') GO TO 120      
      CALL CHANGE(H,N,IDH)          
  120 DO 121 I=1,N       
      DO 121 J=1,N       
  121 A(I,J)=H(I,J)      
   21 CALL MXGIV(A,V,E,P,-1.,N,IDH) 
C  
C  
C THIS SECTION ASSIGNS THE ORBITAL OCCUPATION NUMBERS. ORBITALS      
C WHOSE ENERGIES DIFFER BY 0.00005 OR LESS ARE CONSIDERED TO BE      
C DEGENERATE;;
C  
C  
      NLEV=0  
      NORB=0  
      GO TO 24
   23 NORB=NORB+NDEG(NLEV)          
      IF(NORB.GE.N) GO TO 26        
   24 NLEV=NLEV+1        
      I=NORB+1
      NDEG(NLEV)=1       
      IF(I.GE.N) GO TO 26
      I1=I+1  
      DO 25 K=I1,N       
      DE=ABS(E(I)-E(K))  
      IF(DE.GT.5.0E-05) GO TO 25    
      NDEG(NLEV)=NDEG(NLEV)+1       
   25 CONTINUE
      GO TO 23
   26 DO 27 I=1,N        
   27 NOC(I)=0
      NORB=0  
      DO 32 ILEV=1,NLEV  
      IF(ILEV.EQ.1) GO TO 28        
      NORB=NORB+NDEG(ILEV-1)        
      IF(NORB.GE.N) GO TO 33        
   28 NMULT=NDEG(ILEV)   
      IF(NMULT.EQ.1) GO TO 31       
      DO 30 ITIME=1,2    
      DO 29 K=1,NMULT    
      IF(NE.LE.0) GO TO 33          
      I=NORB+K
      NOC(I)=NOC(I)+1    
      NE=NE-1 
   29 CONTINUE
   30 CONTINUE
      GO TO 32
   31 IF(NE.LE.0) GO TO 33          
      I=NORB+1
      INCR=AMIN0(2,NE)   
      NOC(I)=INCR        
      NE=NE-INCR         
   32 CONTINUE
   33 CONTINUE
C  
C  
C EIGENVALUE-EIGENVECTOR PRINTOUT.  
C FOR ZERO ELECTRONS SKIP BOND ORDER, ETC., CALC-N        
C  
C  
   34 CALL EVPRNT(E,V,NTP,N,IDH)    
C  
C  
C CALC HUCKEL ENERGY     
C PRINT OUT OCCUPANCY, ENERGY       
C  
C  
   35 ETOT=0. 
      DO 36 I=1,N        
      ETOT=ETOT+NOC(I)*E(I)         
   36 CONTINUE
      WRITE(5,37)        
   37 FORMAT(//10X,'ORBITAL OCCUPANCY')        
      DO 38 I=1,N,7      
      K=I+6   
      IF(K.GT.N) K=N     
   38 WRITE(5,39)((J,NOC(J)),J=I,K) 
   39 FORMAT(2X,7('MO',I2,'-'I1,4X))
      WRITE(5,40)ETOT    
   40 FORMAT(/10X,'HUCKEL ENERGY = ',F10.5)    
C  
C  
C CALC BOND ORDER-CHARGE DENSITY MATRIX AND PRINT IT.     
C  
C  
      WRITE(5,41)        
   41 FORMAT(//1H$,5X,'CALCULATE BOND-ORDER MATRIX (Y OR N)>[Y]   ') 
      READ(5,8)ANS       
      IF(ANS.EQ.'N') GO TO 53       
   42 WRITE(5,43)        
   43 FORMAT(//10X,'BOND-ORDER MATRIX BY COLUMNS')        
      DO 45 IP=1,N       
      DO 45 IQ=IP,N      
      SUM=0.  
      DO 44 I=1,N        
   44 SUM=SUM+NOC(I)*V(IP,I)*V(IQ,I)
      P(IP,IQ)=SUM       
      P(IQ,IP)=SUM       
   45 CONTINUE
      CALL PPRNT(P,N,IDH)
C  
C  
C CALC NET CHARGE, FREE VALENCE, AND SPIN DENSITY         
C PRINT OUT THE ABOVE    
C  
C  
      WRITE(5,46)        
      WRITE(5,47)        
   46 FORMAT(///18X,'CHARGE',9X,'NET',8X,'FREE',8X,'SPIN')
   47 FORMAT(17X,'DENSITY',6X,'CHARGE',5X,'VALENCE',5X,'DENSITY'//)  
      DO 50 IP=1,N       
      Q(IP)=ZC(IP)-P(IP,IP)         
      RHO(IP)=0.         
      SUM=0.  
      DO 48 IQ=1,N       
      IF(IP.EQ.IQ)GO TO 48          
      IF(H(IP,IQ).EQ.0.) GO TO 48   
      SUM=SUM+P(IP,IQ)   
   48 CONTINUE
      IT=NTP(IP)         
      FV(IP)=FVO(IT)-SUM 
      DO 49 I=1,N        
      IF(NOC(I).EQ.1) RHO(IP)=RHO(IP)+V(IP,I)**2          
   49 CONTINUE
   50 WRITE(5,51)NAME(IT),IP,P(IP,IP),Q(IP),FV(IP),RHO(IP)
   51 FORMAT(3X,A2,I4,3X,4F12.5)    
      WRITE(5,52)        
   52 FORMAT(/10X,'SPIN DENSITY IS FOR HIGHEST MULTIPLICITY')        
C  
C SECTION FOR CHANGING ORBITAL OCCUPANCY. BOND-ORDERS, ENERGIES, ETC 
C FOR EXCITED STATES CAN BE CALC-D. 
C  
C  
   53 WRITE(5,54)        
   54 FORMAT(//1H$,10X,'NEW ELECTRONIC CONFIGURATION (Y OR N)>[N]   ')           

      READ(5,8)ANS       
      IF(ANS.NE.'Y') GO TO 1        
      WRITE(5,55)        
   55 FORMAT( 
     1/10X,'REFER TO *LATEST* OCCUPANCY LIST.' 
     2/ 10X,'TRANSFER 1 ELECTRON FROM ORBITAL I TO ORBITAL J.'       
     3/10X,'I=0 MEANS INCREASE # ELECTRONS(ADD 1 TO J).'  
     4/10X,' J=0 MEANS DECREASE # ELECTRONS(TAKE 1 FROM I).'         
     5/10X,'BOTH I AND J ZERO MEANS NO MORE TRANSFERS.'   
     6/10X,' I OR J OUT OF RANGE = NO MORE TRANSFERS.'/)  
   56 WRITE(5,57)        
   57 FORMAT(1H$,5X,'ENTER I,J   ') 
      READ(5,58) I,J     
   58 FORMAT(2I)         
      IF(I.LT.0.OR.I.GT.N) GO TO 35 
      IF(J.LT.0.OR.J.GT.N) GO TO 35 
      IF(I.EQ.0.AND.J.EQ.0) GO TO 35
      IF(I.EQ.0) GO TO 62
      IF(NOC(I).LE.0) GO TO 60      
   59 NOC(I)=NOC(I)-1    
      GO TO 62
   60 WRITE(5,61) I      
   61 FORMAT(/10X,'** ORBITAL',I3,' IS EMPTY;;;')         
      GO TO 56
   62 IF(J.EQ.0) GO TO 56
      IF(NOC(J).GE.2)GO TO 64       
   63 NOC(J)=NOC(J)+1    
      GO TO 56
   64 WRITE(5,65) J      
   65 FORMAT(/10X,'** ORBITAL',I3,' IS ALREADY FULL;;;')  
      GO TO 56
   66 STOP 'END OF CALC-N'          
C  
C  
C THIS SECTION TYPES OUT THE *HELP* FILE FOR THE PROGRAM  
C  
C  
   67 WRITE(5,68)        
   68 FORMAT(//15X,'THE PROGRAM BUILDS A CARBON SKELETON FIRST,'     
     1/10X,'INTRODUCING HETEROATOMS IN A LATER SECTION.'  
     2//15X,'MOST SYSTEMS CAN BE BUILT FROM A LINEAR OR MONOCYCLIC'  
     3/10X,'HYDROCARBON BY INTRODUCING CROSSLINKS FOR THE EXTRA BONDS.'          

     4// 5X,'STYRENE = LINEAR, N=8, 1-6 CROSSLINK'        
     5/5X,'ANTHRACENE = MONOCYCLIC, N=14, 1-6, 8-13 CROSSLINKS'      
     6//15X,'IF YOU CHOOSE THIS OPTION YOU NEED TO SUPPLY THE SIZE'  
     7/10X,'OF THE SYSTEM, A Y OR N ANSWER TO THE CYCLIC QUERY, A'   
     8/10X,'Y OR N ANSWER TO THE CROSSLINK QUERY, THE # OF CROSS-'   
     9/10X,'LINKS, AND THE #-S OF ATOMS CONNECTED BY EACH CROSSLINK.'
     1/)      
      WRITE(5,69)        
   69 FORMAT(1H$,10X,' CAN YOU BUILD YOUR SYSTEM THIS WAY(TYPE Y OR N)>          

     1[N]  ') 
      READ(5,8)ANS       
      IF(ANS.EQ.'Y')GO TO 9         
      WRITE(5,70)        
   70 FORMAT(15X,'ANY OTHER SYSTEM MUST BE BUILT USING A SERIES OF'  
     1/10X,'CONNECTED-ATOM LISTS. EACH LIST MAY CONTAIN UP TO 20'    
     2/10X,'ATOM #-S. ANY NUMBER OF LISTS MAY BE USED. SOME EXAMPLES'
     3/10X,'ARE GIVEN BELOW0'       
     4//5X,'STYRENE = 2 LISTS'      
     5/15X,'1 2 3 4 5 6 1'          
     6/15X,'6 7 8'       
     7//5X,'TRIMETHYLENE CYCLOPROPANE = 3 LISTS'          
     8/15X,'4 1 2 5'     
     9/15X,'1 3 2'       
     1/15X,'3 6'         
     2//15X,'THE REQUIRED INPUT IS THE SIZE OF THE SYSTEM, THE NUMBER'           

     3/10X,'OF CONNECTED-ATOM LISTS, AND THE NUMBER OF ATOMS IN EACH'
     4/10X,'LIST ALONG WITH THE ATOM #-S FOR THAT LIST.'  
     5/)      
      WRITE(5,71)        
   71 FORMAT(1H$,10X,'IS YOUR DATA IN THIS FORM(TYPE Y OR N)>[N]   ')
      READ(5,8)ANS       
      IF(ANS.EQ.'Y') GO TO 10       
      WRITE(5,72)        
   72 FORMAT(10X,'YOU HAVE BOMBED OUT;;; GET ORGANIZED;;;;;')        
      GO TO 66
      END     
      SUBROUTINE PPRNT(P,N,ID)      
C  
C  
C A ROUTINE FOR PRINTING OUT A SQUARE MATRIX BY COLUMNS   
C THE FORMAT IS FOR A 72-COLUMN WIDE DEVICE.   
C  
C  
      DIMENSION P(ID,N)  
      IGO=1   
      NEND=0  
    1 NBEG=NEND+1        
      NEND=NEND+6        
      IF(NEND.LT.N) GO TO 3         
    2 NEND=N  
      IGO=2   
    3 WRITE(5,4)((K),K=NBEG,NEND)   
    4 FORMAT(1H0,10X,6(5X,I3,2X))   
      WRITE(5,5)         
    5 FORMAT(1H )        
      DO 6 I=1,N         
    6 WRITE(5,7)I,(P(I,J),J=NBEG,NEND)         
    7 FORMAT(1H ,I8,2X,6F10.5)      
      GO TO (1,8),IGO    
    8 RETURN  
      END     
      SUBROUTINE EASY(H,N,ID)       
      DIMENSION H(ID,ID) 
C  
C A SUBROUTINE FOR BUILDING A HUCKEL MATRIX FOR ANY SYSTEM WHICH MAY BE          

C BE DESCRIBED AS A CROSSLINKED LINEAR OR MONOCYCLIC POLYENE.        
C  
      WRITE(5,1)         
    1 FORMAT(1H$,5X,'SIZE OF HUCKEL SYSTEM (0 OR NEG  = STOP)[0]     
     1    ')  
    2 READ(5,3)N         
      IF(N.LE.0) GO TO 19
      IF(N.GT.ID)GO TO 20
    3 FORMAT(I)          
      NO=N-1  
      DO 4 I=1,NO        
      J=I+1   
      H(I,J)=1.          
    4 H(J,I)=1.          
      WRITE(5,5)         
    5 FORMAT(1H$,5X,'IS THE SYSTEM CYCLIC(Y OR N)>[N]    ')          
      READ(5,6)ANS       
    6 FORMAT(A1)         
      IF(ANS.NE.'Y') GO TO 8        
    7 H(1,N)=1.          
      H(N,1)=1.          
    8 WRITE(5,9)         
    9 FORMAT(1H$,5X,'HOW MANY CROSSLINKS ARE PRESENT>[0]   ')        
      READ(5,3)NC        
      IF(NC.LT.0) NC=0   
      IF(NC.EQ.0) GO TO 18          
   10 DO 16 IC=1,NC      
   11 WRITE(5,12)IC      
   12 FORMAT(1H$,5X,'CROSSLINK',I3,' CONNECTS ATOMS   ')  
      READ(5,17)I,J      
      IF(I.LE.0.OR.J.LE.0) GO TO 13 
      IF(I.GT.N.OR.J.GT.N) GO TO 13 
      GO TO 15
   13 WRITE(5,14)        
   14 FORMAT(5X,'EITHER I OR J  OUT OF RANGE'/)
      GO TO 11
   15 H(I,J)=1.          
   16 H(J,I)=1.          
   17 FORMAT(2I)         
   18 RETURN  
   19 STOP 'N LE 0'      
   20 WRITE(5,21)ID      
   21 FORMAT(1H$,10X,'TOO BIG',I4,' MAX, ENTER NEW N   ') 
      GO TO 2 
      END     
      SUBROUTINE HARD(H,N,ID)       
      DIMENSION JL(30)   
      DIMENSION H(ID,ID) 
C  
C A SUBROUTINE FOR CONSTRUCTING A HUCKEL MATRIX FOR A COMPLEX HYDROCARBO         

C FROM CONNECTED-ATOM LISTS.        
C  
      WRITE(5,1)         
    1 FORMAT(1H$,/5X,'BUILDING WITH CONNECTED-ATOM LISTS;'//         
     15X,'HOW MANY ATOMS IN THE SYSTEM(0 OR NEG = STOP)>[0]   ')     
    2 READ(5,3)N         
      IF(N.LE.0)GO TO 13 
      IF(N.GT.ID)GO TO 14
    3 FORMAT(I)          
      WRITE(5,4)         
    4 FORMAT(1H$,        
     1/5X,'HOW MANY CONNECTED-ATOM LISTS>'     
     2/15X,'(20 ENTRIES MAX PER LIST)[0]   ')  
      READ(5,3)NC        
      IF(NC.LE.0)GO TO 13
      DO 12 IC=1,NC      
    5 WRITE(5,6)IC       
    6 FORMAT(1H$,5X,'HOW MANY ENTRIES IN LIST',I3,' (20 MAX)>[0]   ')
      READ(5,3)NL        
      IF(NL.LE.0)GO TO 13
      WRITE(5,7)IC       
    7 FORMAT(5X,'ENTER LIST',I4)    
      READ(5,11)(JL(K),K=1,NL)      
      DO 9 K=1,NL        
      IF(JL(K).GT.0.AND.JL(K).LE.N) GO TO 9    
      WRITE(5,8) JL(K)   
    8 FORMAT(5X,'LIST ENTRY OUT OF RANGE;',I5) 
      GO TO 5 
    9 CONTINUE
      NLO=NL-1
      DO 10 IF=1,NLO     
      JF=IF+1 
      I=JL(IF)
      J=JL(JF)
      H(I,J)=1.          
   10 H(J,I)=1.          
   11 FORMAT(20I)        
   12 CONTINUE
      RETURN  
   13 STOP 'ZERO OR NEG. SIZE OR LENGTH OF ATOM LIST'     
   14 WRITE(5,15)ID      
   15 FORMAT(1H$,10X,'TOO BIG',I4,'MAX, ENTER NEW N     ')
      GO TO 2 
      END     
      SUBROUTINE EVPRNT(E,V,NTP,N,ID)          
C  
C  
C A ROUTINE FOR PRINTING OUT COLUMNWISE EIGENVALUES AND THEIR        
C RESPECTIVE EIGENVECTORS, ON A 72-COLUMN DEVICE.         
C  
C  
      DIMENSION E(N),V(ID,N),NTP(N) 
      COMMON/LABEL/NAME(13)         
      IGO=1   
      NEND=0  
    1 NBEG=NEND+1        
      NEND=NEND+6        
      IF(NEND.LT.N) GO TO 3         
    2 NEND=N  
      IGO=2   
    3 WRITE(5,4)((K),K=NBEG,NEND)   
    4 FORMAT(1H0,10X,6(4X,'MO-',I3))
      WRITE(5,5)(E(I),I=NBEG,NEND)  
    5 FORMAT(1H0,2X,'E/BETA',2X,6F10.5)        
      WRITE(5,6)         
    6 FORMAT(/)          
      DO 7 IP=1,N        
      IT=NTP(IP)         
    7 WRITE(5,8)NAME(IT),IP,(V(IP,K),K=NBEG,NEND)         
    8 FORMAT(1H ,2X,A2,'-',I3,2X,6F10.6)       
      GO TO(1,9),IGO     
    9 RETURN  
      END     
      SUBROUTINE HETERO(H,NTP,ZC,N,ID)         
C  
C  
C A ROUTINE FOR INTRODUCING HETEROATOMS INTO AN EXISTING HUCKEL      
C SYSTEM. ONE-CENTER TERMS ARE PUT IN. TWO-CENTER TERMS ARE          
C ONLY SUBSTITUTED FOR EXISTING TERMS. NO NEW INTERACTIONS
C MAY BE INCLUDED.       
C THE PARAMETERS ARE BASED ON HINZE-BEVERIDGE PPP CALC-NS.
C  
C  
      DIMENSION H(ID,N),NTP(N),ZC(N)
      DIMENSION A(13),B(13,13),ZCORE(13)       
      DATA ZCORE/1.,0.,1.,2.,1.,2.,2.,1.,1.,2.,1.,2.,2./  
      DATA A/0.,-0.45,0.51,1.37,0.97,2.09,2.71,0.,0.19,0.75,0.46,1.11,           

     11.48/   
      DATA B/ 
     11.00,0.73,1.02,0.89,1.06,0.66,0.52,0.75,0.77,0.76,0.81,0.69,0.62,          

     20.73,0.87,0.66,0.53,0.60,0.35,0.26,0.57,0.53,0.54,0.51,0.44,0.41,          

     31.02,0.66,1.09,0.99,1.14,0.80,0.65,0.72,0.78,0.81,0.83,0.78,0.77,          

     40.89,0.53,0.99,0.98,1.13,0.89,0.77,0.43,0.55,0.64,0.68,0.73,0.80,          

     51.06,0.60,1.14,1.13,1.26,1.02,0.92,0.65,0.75,0.82,0.84,0.85,0.88,          

     60.66,0.35,0.80,0.89,1.02,0.95,0.94,0.24,0.31,0.39,0.43,0.54,0.70,          

     70.52,0.26,0.65,0.77,0.92,0.94,1.04,0.17,0.21,0.22,0.28,0.32,0.51,          

     80.75,0.57,0.72,0.43,0.65,0.24,0.17,0.64,0.62,0.52,0.61,0.40,0.34,          

     90.77,0.53,0.78,0.55,0.75,0.31,0.21,0.62,0.63,0.58,0.65,0.48,0.35,          

     10.76,0.54,0.81,0.64,0.82,0.39,0.22,0.52,0.58,0.63,0.65,0.60,0.55,          

     20.81,0.51,0.83,0.68,0.84,0.43,0.28,0.61,0.65,0.65,0.68,0.58,0.52,          

     30.69,0.44,0.78,0.73,0.85,0.54,0.32,0.40,0.48,0.60,0.58,0.63,0.59,          

     40.62,0.41,0.77,0.80,0.88,0.70,0.51,0.34,0.35,0.55,0.52,0.59,0.68/          

      WRITE(5,1)         
    1 FORMAT(10X,'THE FOLLOWING HETERO-ATOMS HAVE THE ASSIGNED TYPE #.'          

     1//15X,' 1 C  - SP2 CARBON'    
     2/15X,' 2 B  - SP2 BORON'      
     3/15X,' 3 N1 - IMINE NITROGEN (-N=)' 
     4/15X,' 4 N2 - PYRROLE NITROGEN (-N:-)'          
     5/15X,' 5 O1 - CARBONYL OXYGEN'
     6/15X,' 6 O2 - FURAN OXYGEN'   
     7/15X,' 7 F  - FLUORINE'       
     8/15X,' 8 SI - SP2 SILICON'    
     9/15X,' 9 P1 - PHOSPHOROUS, LIKE N1'      
     1/15X,'10 P2 - PHOSPHOROUS. LIKE N2'      
     2/15X,'11 S1 - THIONE SULFUR'  
     3/15X,'12 S2 - THIOPHENE SULFUR'          
     4/15X,'13 CL - CHLORINE'/)     
      WRITE(5,2)         
    2 FORMAT(1H$,10X,'ARE YOUR HETEROATOMS LISTED(Y OR N)>[Y] ')     
      READ(5,3) ANS      
    3 FORMAT(A1)         
      IF( ANS.EQ.'N') GO TO 19      
    4 WRITE(5,5)         
    5 FORMAT(1H$,10X,'HOW MANY HETEROATOMS>[0]  ')        
      READ(5,6)NH        
    6 FORMAT(I)          
      IF(NH.LT.0) NH=0   
      IF(NH.GT.0.AND.NH.LE.N) GO TO 8          
      IF(NH.EQ.0) GO TO 21          
      WRITE(5,7) NH      
    7 FORMAT(5X,'# HETEROATOMS OUTOF RANGE;',I5)          
      GO TO 4 
    8 DO 16 IH=1,NH      
    9 WRITE(5,10)IH      
   10 FORMAT(1H$,10X,'HETEROATOM',I3,', ENTER POSITION # AND TYPE #  ')          

      READ(5,15)IP,IT    
      IF(IP.GT.0.AND.IP.LE.N)GO TO 12          
      WRITE(5,11)IP      
   11 FORMAT(5X,'ATOM POSITION # OUT OF RANGE;',I5)       
      GO TO 9 
   12 IF(IT.GT.0.AND.IT.LE.13) GO TO 14        
      WRITE(5,13)IT      
   13 FORMAT(5X,'HETEROATOM TYPE # OUT OF RANGE;',I5)     
      GO TO 9 
   14 NTP(IP)=IT         
   15 FORMAT(2I)         
   16 CONTINUE
      DO 18 I=1,N        
      IT=NTP(I)          
      IF(IT.EQ.1) GO TO 18          
      H(I,I)=A(IT)       
      ZC(I)=ZCORE(IT)    
      DO 17 J=1,N        
      IF(I.EQ.J) GO TO 17
      IF(H(I,J).EQ.0.) GO TO 17     
      JT=NTP(J)          
      H(I,J)=B(IT,JT)    
      H(J,I)=H(I,J)      
   17 CONTINUE
   18 CONTINUE
      RETURN  
   19 WRITE(5,20)        
   20 FORMAT(5X,'NO PARAMETERS FOR HETEROATOM')
      STOP    
   21 WRITE(5,22)        
   22 FORMAT(5X,'NO HETEROATOMS REQUESTED>>>') 
      WRITE(5,23)        
   23 FORMAT(1H$,5X,'CONTINUE CALC-N FOR HYDROCARBON(Y OR N)>[N]  ') 
      READ(5,3)ANS       
      IF(ANS.EQ.'Y') RETURN         
      STOP    
      END     
      SUBROUTINE CHANGE(H,N,IDIM)   
      DIMENSION H(IDIM,N)
C  
C  
C THIS ROUTINE IS USED TO CHANGE SELECTED ELEMENTS OF A HUCKEL MATRIX
C  
C  
      NTIME=1 
C  
C          PRINT THE CURRENT HUCKEL MATRIX     
C  
    1 WRITE(5,9999)      
 9999 FORMAT(1H0,'THIS IS THE CURRENT HUCKEL MATRIX')     
      CALL PPRNT(H,N,IDIM)          
      WRITE(5,9998)      
      READ(5,8999) ANS   
 9998 FORMAT(/1H$,10X,'IS THE MATRIX OK(Y OR N)>[Y]  ')   
 8999 FORMAT(A1)         
      IF(ANS.NE.'N') GO TO 900      
C  
C  
C ***** CHANGE SECTION *****        
C  
C  
      GO TO (5,10),NTIME 
    5 WRITE(5,9997)      
 9997 FORMAT(1H0,10X,'ENTER ROW(I) AND COL(J) FOR EACH ELEMENT TO BE',/,
     1 5X,'CHANGED, ONE AT A TIME. ZERO I OR J WILL END CHANGES.')    
      NTIME=2 
C  
C          GET ROW AND COLUMN INDICES          
C  
   10 WRITE(5,9996)      
 9996 FORMAT(/1H$,10X,'ENTER I,J   ')          
      READ(5,8998)I,J    
 8998 FORMAT(2I)         
C  
C          CHECK TO MAKE SURE INDICES ARE IN THE RANGE 1 TO N        
C  
      IF(I.LE.N.AND.J.LE.N)GO TO 15 
      WRITE(5,9995)      
 9995 FORMAT(1H ,10X,'I OR J GTR. THAN # AO-S;; TRY AGAIN;')         
      GO TO 10
   15 IF(I.LE.0.OR.J.LE.0) GO TO 1  
C  
C         PRINT OLD VALUE, ACCEPT NEW VALUE, SYMMETRIZE.  
C         NULL VALUE IS THE SAME AS ZERO;;;    
C  
      WRITE(5,9994)H(I,J)
 9994 FORMAT(1H$,10X,'OLD VALUE = ',F6.2,',NEW VAUE=>[0]   ')        
      READ(5,8997)H(I,J) 
 8997 FORMAT(F)          
      H(J,I)=H(I,J)      
      GO TO 10
  900 RETURN  
      END     
      SUBROUTINE MXGIV(A,V,E,TMP,XSIGN,N,ID)   
C  
C A ROUTINE FOR DIAGONALIZING A REAL SYMMETRIC MATRIX     
C USING A GIVENS ROUTINES.          
C  
C      A - THE ARRAY(2-D) HOLDING THE MATRIX   
C  
C      V - THE ARRAY(2-D) THAT WILL HOLD THE EIGENVECTORS,
C          STORED BY COLUMNS.       
C  
C      E - THE 1-D ARRAY HOLDING THE EIGENVALUES.         
C  
C    TMP - A TEMPORARY STORAGE ARRAY USED DURING THE CALC-N.         
C          FOR MOST MO PROBLEMS IT IS CONVENIENT TO USE THE          
C          DENSITY MATRIX FOR THIS PURPOSE.    
C  
C  XSIGN - A CONSTANT THAT CONTROLS THE ORDERING OF THE EIGENVALUES  
C          XSIGN = -1, E(1) IS THE MOST POSITIVE EIGENVALUE          
C          XSIGN = +1, E(1) IS THE MOST NEGATIVE EIGENVALUE          
C  
C      N - THE ORDER OF THE PROBLEM.
C  
C     ID - THE FIRST DIMENSION OF THE ARRAYS A AND V.     
C  
      DIMENSION A(ID,1),E(N),V(ID,N),TMP(1)    
C  
C  
C TRAP THE TRIVIAL CASE FOR N=1     
C  
C  
      IF(N.GT.1) GO TO 1 
      E(1)=A(1,1)        
      V(1,1)=1.          
      RETURN  
    1 CONTINUE
C  
C PACK THE ARRAY A (UPPER TRIANGULAR) BY COLUMNS INTO TMP.
C  
C  
      IB=0    
      DO 2 IC=1,N        
      DO 2 IR=1,IC       
       IB=IB+1
    2 TMP(IB)=A(IR,IC)*XSIGN        
C  
C  
C THIS SHOULD BE THE SINGLE PRECISION VERSION [GIVENS.FOR]
C  
C  
      CALL GIVENS(N,N,ID,TMP,A,E,V) 
      DO 5 I=1,N         
    5 E(I)=E(I)*XSIGN    
      RETURN  
  900 J=IER-128          
      WRITE(5,9997)J     
 9997 FORMAT(1H0,10X,'CONVERGENCE FAILURE AT EIGENVALUE',I5)         
      STOP    
      END     
      SUBROUTINE GIVENS (NX,NROOTX,NJX,A,B,ROOT,VECT)     
C 62.3  GIVENS  -EIGENVALUES AND EIGENVECTORS BY THE GIVENS METHOD.  
C      BY FRANKLIN PROSSER, INDIANA UNIVERSITY.
C      SEPTEMBER, 1967   
C      CALCULATES EIGENVALUES AND EIGENVECTORS OF REAL SYMMETRIC MATRIX          

C      STORED IN PACKED UPPER TRIANGULAR FORM. 
C  
C      THANKS ARE DUE TO F. E. HARRIS (STANFORD UNIVERSITY) AND H. H.
C      MICHELS (UNITED AIRCRAFT RESEARCH LABORATORIES) FOR EXCELLENT 
C      WORK ON NUMERICAL DIFFICULTIES WITH EARLIER VERSIONS OF THIS  
C      PROGRAM.          
C  
C      THE PARAMETERS FOR THE ROUTINE ARE...   
C          NX     ORDER OF MATRIX   
C          NROOTX NUMBER OF ROOTS WANTED.  THE NROOTX SMALLEST (MOST 
C       NEGATIVE) ROOTS WILL BE CALCULATED.  IF NO VECTORS
C       ARE WANTED, MAKE THIS NUMBER NEGATIVE. 
C          NJX    ROW DIMENSION OF VECT ARRAY.  SEE 'VECT' BELOW.    
C       NJX MUST BE NOT LESS THAN NX.          
C          A      MATRIX STORED BY COLUMNS IN PACKED UPPER TRIANGULAR
C      FORM, I.E. OCCUPYING NX*(NX+1)/2 CONSECUTIVE       
C      LOCATIONS.        
C          B      SCRATCH ARRAY USED BY GIVENS.  MUST BE AT LEAST    
C       NX*5 CELLS.      
C          ROOT   ARRAY TO HOLD THE EIGENVALUES.  MUST BE AT LEAST   
C      NROOTX CELLS LONG.  THE NROOTX SMALLEST ROOTS ARE  
C       ORDERED LARGEST FIRST IN THIS ARRAY.   
C          VECT   EIGENVECTOR ARRAY.  EACH COLUMN WILL HOLD AN       
C       EIGENVECTOR FOR THE CORRESPONDING ROOT.  MUST BE  
C       DIMENSIONED WITH 'NJX' ROWS AND AT LEAST 'NROOTX' 
C       COLUMNS, UNLESS NO VECTORS  
C       ARE REQUESTED (NEGATIVE NROOTX).  IN THIS LATTER  
C       CASE, THE ARGUMENT VECT IS JUST A DUMMY, AND THE  
C       STORAGE IS NOT USED.        
C  
C      THE ARRAYS A AND B ARE DESTROYED BY THE COMPUTATION.  THE RESULTS         

C      APPEAR IN ROOT AND VECT.     
C  
C      FOR PROPER FUNCTIONING OF THIS ROUTINE, THE RESULT OF A FLOATING          

C      POINT UNDERFLOW SHOULD BE A ZERO.       
C  
C      TO CONVERT THIS ROUTINE TO DOUBLE PRECISION (E.G. ON IBM 360  
C      MACHINES), BE SURE THAT ALL REAL VARIABLES AND FUNCTION       
C      REFERENCES ARE PROPERLY MADE DOUBLE PRECISION.     
C  
C      THE ORIGINAL REFERENCE TO THE GIVENS TECHNIQUE IS IN OAK RIDGE
C      REPORT NUMBER ORNL 1574 (PHYSICS), BY WALLACE GIVENS.         
C      THE METHOD AS PRESENTED IN THIS PROGRAM CONSISTS OF FOUR STEPS,           

C      ALL MODIFICATIONS OF THE ORIGINAL METHOD...        
C      FIRST, THE INPUT MATRIX IS REDUCED TO TRIDIAGONAL FORM BY THE 
C      HOUSEHOLDER TECHNIQUE (J. H. WILKINSON, COMP. J. 3, 23 (1960)).           

C      THE ROOTS ARE THEN LOCATED BY THE STURM SEQUENCE METHOD (J. M.
C      ORTEGA (SEE REFERENCE BELOW).  THE VECTORS OF THE TRIDIAGONAL 
C      FORM ARE THEN EVALUATED (J. H. WILKINSON, COMP. J. 1, 90 (1958)),         

C      AND LAST THE TRIDIAGONAL VECTORS ARE ROTATED TO VECTORS OF THE
C      ORIGINAL ARRAY (FIRST REFERENCE).       
C      VECTORS FOR DEGENERATE (OR NEAR-DEGENERATE) ROOTS ARE FORCED  
C      TO BE ORTHOGONAL, USING A METHOD SUGGESTED BY B. GARBOW, ARGONNE          

C      NATIONAL LABS (PRIVATE COMMUNICATION, 1964).  THE GRAM-SCHMIDT
C      PROCESS IS USED FOR THE ORTHOGONALIZATION.         
C  
C      AN EXCELLENT PRESENTATION OF THE GIVENS TECHNIQUE IS FOUND IN 
C      J. M. ORTEGA'S ARTICLE IN 'MATHEMATICS FOR DIGITAL COMPUTERS,'
C      VOLUME 2, ED. BY RALSTON AND WILF, WILEY (1967), PAGE 94.     
C  
      DIMENSION B(NX,5),A(1),ROOT(NROOTX),VECT(NJX,NROOTX)
C  
C ** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *         

C **   USERS PLEASE NOTE...         
C **   THE FOLLOWING TWO PARAMETERS, ETA AND THETA, SHOULD BE ADJUSTED           

C **   BY THE USER FOR HIS PARTICULAR MACHINE. 
C **   ETA IS AN INDICATION OF THE PRECISION OF THE FLOATING POINT   
C **   REPRESENTATION ON THE COMPUTER BEING USED (ROUGHLY 10**(-M),  
C **   WHERE M IS THE NUMBER OF DECIMALS OF PRECISION ).  
C **   THETA IS AN INDICATION OF THE RANGE OF NUMBERS THAT CAN BE    
C **   EXPRESSED IN THE FLOATING POINT REPRESENTATION (ROUGHLY THE   
C **   LARGEST POSITIVE NUMBER OR THE RECIPROCAL OF THE SMALLEST     
C **   POSITIVE NUMBER, WHICHEVER IS SMALLER). 
C **   SOME RECOMMENDED VALUES FOLLOW.         
C **   FOR IBM 7094, UNIVAC 1108, ETC. (27-BIT BINARY FRACTION, 8-BIT
C **   BINARY EXPONENT), ETA=1.E-8, THETA=1.E37.          
C **   FOR CONTROL DATA 3600 (36-BIT BINARY FRACTION, 11-BIT BINARY  
C **   EXPONENT), ETA=1.E-11, THETA=1.E307.    
C **   FOR CONTROL DATA 6600 (48-BIT BINARY FRACTION, 11-BIT BINARY  
C **   EXPONENT), ETA=1.E-14, THETA=1.E293.    
C **   FOR IBM 360/50 AND 360/65 DOUBLE PRECISION (56-BIT HEXADECIMAL
C **   FRACTION, 7-BIT HEXADECIMAL EXPONENT), ETA=1.E-16, THETA=1.E75.           

C **          
C  
C      THE FOLLOWING ARE THE VAX 11/780 VALUES 
C  
      ETA=1.E-7          
      THETA=1.E38        
C ** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *         

C  
      DEL1 = ETA/100.    
      DELTA = ETA**2*100.
      SMALL = ETA**2/100.
      DELBIG = THETA*DELTA/1000.    
      THETA1 = 1000./THETA          
C      TOLER  IS A FACTOR USED TO DETERMINE IF TWO ROOTS ARE CLOSE   
C      ENOUGH TO BE CONSIDERED DEGENERATE FOR PURPOSES OF ORTHOGONALI-           

C      ZING THEIR VECTORS.  FOR THE MATRIX NORMED TO UNITY, IF THE   
C      DIFFERENCE BETWEEN TWO ROOTS IS LESS THAN TOLER, THEN         
C      ORTHOGONALIZATION WILL OCCUR.
      TOLER = ETA*100.   
C  
C      INITIAL VALUE FOR PSEUDORANDOM NUMBER GENERATOR... (2**23)-3  
      RPOWER = 8388608.  
      RPOW1 = RPOWER/2.  
      RAND1 = RPOWER - 3.
C  
      N = NX  
      NROOT = IABS(NROOTX)          
      IF (NROOT.EQ.0) GO TO 78      
      IF (N-1) 78,1,2    
    1 ROOT(1) = A(1)     
      IF (NROOTX.GT.0) VECT(1,1) = 1.0         
      GO TO 78
    2 CONTINUE
C     NSIZE    NUMBER OF ELEMENTS IN THE PACKED ARRAY     
      NSIZE = (N*(N+1))/2
      NM1 = N-1          
      NM2 = N-2          
C  
C     SCALE MATRIX TO EUCLIDEAN NORM OF 1.  SCALE FACTOR IS ANORM.   
      FACTOR = 0.        
      DO 3 I=1,NSIZE     
    3 FACTOR = AMAX1(FACTOR,ABS(A(I)))         
      IF (FACTOR.NE.0.) GO TO 6     
C     NULL MATRIX.  FIX UP ROOTS AND VECTORS, THEN EXIT.  
      DO 5 I=1,NROOT     
      IF (NROOTX.LT.0) GO TO 5      
      DO 4 J=1,N         
    4 VECT(J,I) = 0.     
      VECT(I,I) = 1.0    
    5 ROOT(I) = 0.       
      GO TO 78
C  
    6 ANORM = 0.         
      J = 1   
      K = 1   
    7 DO 9 I=1,NSIZE     
      IF (I.NE.J) GO TO 8
      ANORM = ANORM + (A(I)/FACTOR)**2/2.      
      K = K+1 
      J = J+K 
      GO TO 9 
    8 ANORM = ANORM + (A(I)/FACTOR)**2         
    9 CONTINUE
   10 ANORM = SQRT(ANORM*2.)*FACTOR 
      DO 11 I=1,NSIZE    
   11 A(I) = A(I)/ANORM  
      ALIMIT = 1.0       
C  
C      TRIDIA SECTION.   
C      TRIDIAGONALIZATION OF SYMMETRIC MATRIX  
      ID = 0  
      IA = 1  
      IF (NM2.EQ.0) GO TO 24        
      DO 23  J=1,NM2     
C      J       COUNTS ROW  OF A-MATRIX TO BE DIAGONALIZED 
C      IA      START OF NON-CODIAGONAL ELEMENTS IN THE ROW
C      ID      INDEX OF CODIAGONAL ELEMENT ON ROW BEING CODIAGONALIZED.          

      IA = IA+J+2        
      ID = ID + J + 1    
      JP2 = J+2          
C      SUM SQUARES OF NON-CODIAGONAL ELEMENTS IN ROW J    
      II = IA 
      SUM = 0.0          
      DO 12 I=JP2,N      
      SUM = SUM + A(II)**2          
   12 II = II + I        
      TEMP = A(ID)       
      IF (SUM.GT.SMALL) GO TO 14    
C      NO TRANSFORMATION NECESSARY IF ALL THE NON-CODIAGONAL         
C      ELEMENTS ARE TINY.
   13 B(J,1) = TEMP      
      A(ID) = 0.         
      GO TO 23
C      NOW COMPLETE THE SUM OF OFF-DIAGONAL SQUARES       
   14 SUM = SQRT(SUM + TEMP**2)     
C      NEW CODIAGONAL ELEMENT       
      B(J,1) = -SIGN(SUM,TEMP)      
C      FIRST NON-ZERO ELEMENT OF THIS W-VECTOR 
      B(J+1,2) = SQRT((1.0 + ABS(TEMP)/SUM)/2.0)          
C      FORM REST OF THE W-VECTOR ELEMENTS      
      TEMP = SIGN(0.5/(B(J+1,2)*SUM),TEMP)     
      II = IA 
      DO 15 I=JP2,N      
      B(I,2) = A(II)*TEMP
   15 II = II + I        
C      FORM P-VECTOR AND SCALAR.  P-VECTOR = A-MATRIX*W-VECTOR.      
C      SCALAR = W-VECTOR*P-VECTOR.  
      AK = 0.0
C      IC      LOCATION OF NEXT DIAGONAL ELEMENT          
      IC = ID + 1        
      J1 = J + 1         
      DO 19  I=J1,N      
      JJ = IC 
      TEMP = 0.          
      DO 18  II=J1,N     
C      I       RUNS OVER THE NON-ZERO P-ELEMENTS          
C      II      RUNS OVER ELEMENTS OF W-VECTOR  
      TEMP = TEMP + B(II,2)*A(JJ)   
C      CHANGE INCREMENTING MODE AT THE DIAGONAL ELEMENTS. 
      IF (II.LT.I) GO TO 17         
   16 JJ = JJ + II       
      GO TO 18
   17 JJ = JJ + 1        
   18 CONTINUE
C      BUILD UP THE K-SCALAR (AK)   
      AK = AK + TEMP*B(I,2)         
      B(I,1) = TEMP      
C      MOVE IC TO TOP OF NEXT A-MATRIX 'ROW'   
   19 IC = IC + I        
C      FORM THE Q-VECTOR 
      DO 20  I=J1,N      
   20 B(I,1) = B(I,1) - AK*B(I,2)   
C      TRANSFORM THE REST OF THE A-MATRIX      
C      JJ      START-1 OF THE REST OF THE A-MATRIX        
      JJ = ID 
C      MOVE W-VECTOR INTO THE OLD A-MATRIX LOCATIONS TO SAVE SPACE   
C      I       RUNS OVER THE SIGNIFICANT ELEMENTS OF THE W-VECTOR    
      DO 22  I=J1,N      
      A(JJ) = B(I,2)     
      DO 21  II=J1,I     
      JJ = JJ + 1        
   21 A(JJ) = A(JJ) - 2.0*(B(I,1)*B(II,2) + B(I,2)*B(II,1))          
   22 JJ = JJ + J        
   23 CONTINUE
C      MOVE LAST CODIAGONAL ELEMENT OUT INTO ITS PROPER PLACE        
   24 CONTINUE
      B(NM1,1) = A(NSIZE-1)         
      A(NSIZE-1) = 0.    
      B(N,1)=0.          
C  
C     STURM SECTION.     
C     STURM SEQUENCE ITERATION TO OBTAIN ROOTS OF TRIDIAGONAL FORM.  
C     MOVE DIAGONAL ELEMENTS INTO SECOND N ELEMENTS OF B-VECTOR.     
C     THIS IS A MORE CONVENIENT INDEXING POSITION.        
C     ALSO, PUT SQUARE OF CODIAGONAL ELEMENTS IN THIRD N ELEMENTS.   
      JUMP=1  
      DO 25 J=1,N        
      B(J,2)=A(JUMP)     
      B(J,3) = B(J,1)**2 
   25 JUMP = JUMP+J+1    
      DO 26 I=1,NROOT    
   26 ROOT(I) = +ALIMIT  
      ROOTL = -ALIMIT    
C     ISOLATE THE ROOTS.  THE NROOT LOWEST ROOTS ARE FOUND, LOWEST FIRST         

      DO 35 I=1,NROOT    
C     FIND CURRENT 'BEST' UPPER BOUND          
      ROOTX = +ALIMIT    
      DO 27 J=I,NROOT    
   27 ROOTX = AMIN1(ROOTX,ROOT(J))  
      ROOT(I) = ROOTX    
C     GET IMPROVED TRIAL ROOT       
   28 TRIAL = (ROOTL + ROOT(I))*0.5 
      IF (TRIAL.EQ.ROOTL.OR.TRIAL.EQ.ROOT(I)) GO TO 35    
C     FORM STURM SEQUENCE RATIOS, USING ORTEGA'S ALGORITHM (MODIFIED).           

C     NOMTCH IS THE NUMBER OF ROOTS LESS THAN THE TRIAL VALUE.       
   29 CONTINUE
      NOMTCH = N         
      J = 1   
   30 F0 = B(J,2) - TRIAL
   31 CONTINUE
      IF (ABS(F0).LT.THETA1) GO TO 32          
      IF (F0.GE.0.) NOMTCH = NOMTCH - 1        
      J = J + 1          
      IF (J.GT.N) GO TO 33          
C     SINCE MATRIX IS NORMED TO UNITY, MAGNITUDE OF B(J,3) IS LESS THAN          

C     ONE, SO OVERFLOW IS NOT POSSIBLE AT THE DIVISION STEP, SINCE   
C     F0 IS GREATER THAN THETA1.    
C     F0 = B(J,2) - TRIAL - B(J-1,3)/F0        
      BB = B(J-1,3)/F0   
      IF (F0.EQ.0) BB=0. 
      F0 = B(J,2) - TRIAL - BB      
      GO TO 31
   32 J = J + 2          
      NOMTCH = NOMTCH - 1
      IF (J.LE.N) GO TO 30          
   33 CONTINUE
C     FIX NEW BOUNDS ON ROOTS       
      IF (NOMTCH.GE.I) GO TO 34     
      ROOTL = TRIAL      
      GO TO 28
   34 ROOT(I) = TRIAL    
      NOM = MIN0(NROOT,NOMTCH)      
      ROOT(NOM) = TRIAL  
      GO TO 28
   35 CONTINUE
C  
C     TRIVEC SECTION.    
C     EIGENVECTORS OF CODIAGONAL FORM          
   37 CONTINUE
C     QUIT NOW IF NO VECTORS WERE REQUESTED.   
      IF (NROOTX.LT.0) GO TO 76     
C     INITIALIZE VECTOR ARRAY.      
      DO 38 I=1,N        
      DO 38 J=1,NROOT    
   38 VECT(I,J) = 1.0    
      DO 70 I=1,NROOT    
      AROOT = ROOT(I)    
C     ORTHOGONALIZE IF ROOTS ARE CLOSE.        
      IF (I.EQ.1) GO TO 40          
C     THE ABSOLUTE VALUE IN THE NEXT TEST IS TO ASSURE THAT THE TRIVEC           

C     SECTION IS INDEPENDENT OF THE ORDER OF THE EIGENVALUES.                    

   39 IF (ABS(ROOT(I-1)-AROOT).LT.TOLER) GO TO 41                                

   40 IA = -1                                                                    

   41 IA = IA + 1                                                                

      ELIM1 = A(1) - AROOT                                                       

      ELIM2 = B(1,1)                                                             

      JUMP = 1                                                                   

      DO 44  J=1,NM1                                                             

      JUMP = JUMP+J+1                                                            

C     GET THE CORRECT PIVOT EQUATION FOR THIS STEP.                              

      IF (ABS(ELIM1).LE.ABS(B(J,1))) GO TO 42                                    

C     FIRST (ELIM1) EQUATION IS THE PIVOT THIS TIME.  CASE 1.                    

      B(J,2) = ELIM1                                                             

      B(J,3) = ELIM2                                                             

      B(J,4) = 0.                                                                

      TEMP = B(J,1)/ELIM1                                                        

      ELIM1 = A(JUMP) - AROOT - TEMP*ELIM2                                       

      ELIM2 = B(J+1,1)                                                           

      GO TO 43                                                                   

C     SECOND EQUATION IS THE PIVOT THIS TIME.  CASE 2.                           

   42 B(J,2) = B(J,1)                                                            

      B(J,3) = A(JUMP) - AROOT                                                   

      B(J,4) = B(J+1,1)                                                          

      TEMP = 1.0                                                                 

      IF (ABS(B(J,1)).GT.THETA1) TEMP = ELIM1/B(J,1)                             

      ELIM1 = ELIM2 - TEMP*B(J,3)                                                

      ELIM2 = -TEMP*B(J+1,1)                                                     

C     SAVE FACTOR FOR THE SECOND ITERATION.                                      

   43 B(J,5) = TEMP                                                              

   44 CONTINUE                                                                   

      B(N,2) = ELIM1                                                             

      B(N,3) = 0.                                                                

      B(N,4) = 0.                                                                

      B(NM1,4) = 0.                                                              

      ITER = 1                                                                   

      IF (IA.NE.0) GO TO 59                                                      

C     BACK SUBSTITUTE TO GET THIS VECTOR.                                        

   45 L = N + 1                                                                  

      DO 53 J=1,N                                                                

      L = L - 1                                                                  

   46 CONTINUE                                                                   

      IF (L-NM1) 49,48,47                                                        

   47 ELIM1=VECT(L,I)                                                            

      GO TO 50                                                                   

   48 ELIM1=VECT(L,I)-VECT(L+1,I)*B(L,3)                                         

      GO TO 50                                                                   

   49 ELIM1=VECT(L,I)-VECT(L+1,I)*B(L,3)-VECT(L+2,I)*B(L,4)                      

   50 CONTINUE                                                                   

C     IF OVERFLOW IS CONCEIVABLE, SCALE THE VECTOR DOWN.                         

C     THIS APPROACH IS USED TO AVOID MACHINE-DEPENDENT AND SYSTEM-               

C     DEPENDENT CALLS TO OVERFLOW ROUTINES.                                      

      IF (ABS(ELIM1).GT.DELBIG) GO TO 51                                         

      TEMP = B(L,2)                                                              

      IF (ABS(B(L,2)).LT.DELTA) TEMP = DELTA                                     

      VECT(L,I) = ELIM1/TEMP                                                     

      GO TO 53                                                                   

C     VECTOR IS TOO BIG.  SCALE IT DOWN.                                         

   51 DO 52 K=1,N                                                                

   52 VECT(K,I) = VECT(K,I)/DELBIG                                               

      GO TO 46                                                                   

   53 CONTINUE                                                                   

      GO TO (54,61), ITER                                                        

C     SECOND ITERATION.  (BOTH ITERATIONS FOR REPEATED-ROOT VECTORS).            

   54 ITER = ITER + 1                                                            

   55 ELIM1 = VECT(1,I)                                                          

      DO 58 J=1,NM1                                                              

      IF (B(J,2).EQ.B(J,1)) GO TO 57                                             

C     CASE ONE.                                                                  

   56 VECT(J,I) = ELIM1                                                          

      ELIM1 = VECT(J+1,I) - ELIM1*B(J,5)                                         

      GO TO 58                                                                   

C     CASE TWO.                                                                  

   57 VECT(J,I) = VECT(J+1,I)                                                    

      ELIM1 = ELIM1 - VECT(J+1,I)*B(J,5)                                         

   58 CONTINUE                                                                   

      VECT(N,I) = ELIM1                                                          

      GO TO 45                                                                   

C     PRODUCE A RANDOM VECTOR                                                    

   59 CONTINUE                                                                   

      DO 60 J=1,N                                                                

C     GENERATE PSEUDORANDOM NUMBERS WITH UNIFORM DISTRIBUTION IN (-1,1).         

C     THIS RANDOM NUMBER SCHEME IS OF THE FORM...                                

C     RAND1 = AMOD((2**12+3)*RAND1,2**23)                                        

C     IT HAS A PERIOD OF 2**21 NUMBERS.                                          

      RAND1 = AMOD(4099.*RAND1,RPOWER)                                           

   60 VECT(J,I) = RAND1/RPOW1 - 1.                                               

      GO TO 45                                                                   

C                                                                                

C     ORTHOGONALIZE THIS REPEATED-ROOT VECTOR TO OTHERS WITH THIS ROOT.          

   61 IF (IA.EQ.0) GO TO 65                                                      

      DO 64 J1=1,IA                                                              

      K = I - J1                                                                 

      TEMP = 0.                                                                  

      DO 62 J=1,N                                                                

   62 TEMP = TEMP + VECT(J,I)*VECT(J,K)                                          

      DO 63 J=1,N                                                                

   63 VECT(J,I) = VECT(J,I) - TEMP*VECT(J,K)                                     

   64 CONTINUE                                                                   

   65 GO TO (55,66), ITER                                                        

C     NORMALIZE THE VECTOR                                                       

   66 ELIM1 = 0.                                                                 

      DO 67 J=1,N                                                                

   67 ELIM1 = AMAX1(ABS(VECT(J,I)),ELIM1)                                        

      TEMP=0.                                                                    

      DO 68 J=1,N                                                                

      ELIM2=VECT(J,I)/ELIM1                                                      

   68 TEMP = TEMP + ELIM2**2                                                     

      TEMP=1.0/(SQRT(TEMP)*ELIM1)                                                

      DO 69 J=1,N                                                                

      VECT(J,I) = VECT(J,I)*TEMP                                                 

      IF (ABS(VECT(J,I)).LT.DEL1) VECT(J,I) = 0.                                 

   69 CONTINUE                                                                   

   70 CONTINUE                                                                   

C                                                                                

C      SIMVEC SECTION.                                                           

C      ROTATE CODIAGONAL VECTORS INTO VECTORS OF ORIGINAL ARRAY                  

C      LOOP OVER ALL THE TRANSFORMATION VECTORS                                  

      IF (NM2.EQ.0) GO TO 76                                                     

      JUMP = NSIZE - (N+1)                                                       

      IM = NM1                                                                   

      DO 75  I=1,NM2                                                             

      J1 = JUMP                                                                  

C      MOVE A TRANSFORMATION VECTOR OUT INTO BETTER INDEXING POSITION.           

      DO 71  J=IM,N                                                              

      B(J,2) = A(J1)                                                             

   71 J1 = J1 + J                                                                

C      MODIFY ALL REQUESTED VECTORS.                                             

      DO 74  K=1,NROOT                                                           

      TEMP = 0.                                                                  

C      FORM SCALAR PRODUCT OF TRANSFORMATION VECTOR WITH EIGENVECTOR             

      DO 72  J=IM,N                                                              

   72 TEMP = TEMP + B(J,2)*VECT(J,K)                                             

      TEMP = TEMP + TEMP                                                         

      DO 73  J=IM,N                                                              

   73 VECT(J,K) = VECT(J,K) - TEMP*B(J,2)                                        

   74 CONTINUE                                                                   

      JUMP = JUMP - IM                                                           

   75 IM = IM - 1                                                                

   76 CONTINUE                                                                   

C      RESTORE ROOTS TO THEIR PROPER SIZE.                                       

      DO 77 I=1,NROOT                                                            

   77 ROOT(I) = ROOT(I)*ANORM                                                    

   78 RETURN                                                                     

      END