      PROGRAM LCOAO
c----GC-C-MCD--06.OCT.1983--Update-Dec.2001/Dec.2005--------JS-L  
c                                                                       
c    ===========================================================              
c       Linear Combination of Orthogonalized Atomic Orbitals              
c    ===========================================================               
c
c    Jens Spanget-Larsen: "The alternant hydrocarbon pairing 
c    theorem and all-valence electrons theory. An approximate LCOAO
c    theory for the electronic absorption and MCD spectra of
c    conjugated organic compounds"
c    Part 1. Croat. Chem. Acta 59, 711-717 (1986).
c    Part 2. Theor. Chem. Acc. 98, 137-153 (1997).
c                                                                                                          
c  - Restricted Open Shell Formalism                                    
c  - Grand Canonical/Canonical Monoexcited CI for Open Shell Systems            
c                                                                      
c    Development version, not user friendly!
c
c    Array Dimensions: 
c		  100 = Max. Number of atoms
c		  300 = Max. Number of orbitals
c		  101 = 100+1
c		  201 = 2*100+1
c		  301 = 3*100+1
c		  400 = 4*100
c   		10000 = 100*100   
c		 5050 = 100(100+1)/2
C		90000 = 300*300
c
c    INPUT FILE - 
c
c 1) FORMAT(18(A4)): text heading
c 2) FORMAT(5I3):
c      NE       charge of the molecule
c      NPE      nr. of electrons in partly occupied MOs
c               (defaults: NPE=0 for even, NPE=1 for odd 
c                number of electrons)
c      NPMO     nr. of partly occupied MOs
c               (defaults: NPMO=0 for even, NPMO=1 for odd 
c                number of electrons)
c      IPI      If IPI>0, only PI-PI* cfgs. included in CI
c               (molecular plane assumed to be X,Y) 
c      IPRNT    print code (IPRNT=0-4 generates min-max printout)
c      NMOP     nr. of HOMOs and nr. of LUMOs printed (default 32)
c 3) One line per atom, FORMAT(I2,F13.6,2F15.6):
c      IQ       atomic number
c      X,Y,Z    atomic coordinates in Angstrom
c 4) One blank line, terminating input of atomic coordinates
c 5) FORMAT(6I3,F10.5):
c      IGCCI    CI calculation performed if IGCCI=1
c      NB       nr. of highest occ. valence MOs incl. in CI (default 10)         
c      NA       nr. of lowest unocc. valence MOs incl. in CI (default 10)
c               (note that NB and NA define an all-valence MO 'window', 
c               and in general, NB*NA singly excited cfgs. are generated; 
c               but if IPI>0, only those cfgs. that are of PI-PI* type are 
c               included in the CI)
c      ITRIP    triplet states computed if ITRIP=1 
c      IDOB     dummy in this version of the program
c      MCD      MCD B-terms for PI-PI* states computed if MCD=>0 (only
c               if IPI>0; MCD=1-5 generates min-max printout)
c      EMAX     excited cfg. energy cutoff limit in eV (default 15)
c  
c----------------------------------------------------------------JS-L  
      implicit real*8(A-H,O-Z)
      character*14 nin,nout                                               
      COMMON /A/ N,NH,NE,NQ,ECORE,ICOUNT,JCOUNT,EMAX,IPRNT,IGCCI,       
     1           EATOM,ISCF,ITGT,ICI,GAH,ITRIP,IPEN,ETOT,ETOT1,IDOB,    
     2           NB,NAB                                                 
      COMMON /B/ E(300),EE(300),ZCORE(300),ET(300),SUM(300),C(400),          
     1           GAX(300),NTIT(20)                                       
      COMMON /C/ S(300,300),F(300,300)                                      
      COMMON /MOM/ SX(300),SY(300),SZ(300)                                 
      COMMON/E/NCYCL,EPREV,ITERR,N5,NTOTAL,N2S,N2P,NTN,NG,DAMP          
      COMMON /OCCDAT/ NHOMO,NLUMO,NPE,NPMO,OCC(300)                      
      COMMON /D/ V1(300),V2(300),V3(300),V4(300),V5(300),                    
     1  INDEX(2,300),INZ(2,300),DIAG1(300),DIAG(400)                       
      COMMON /F/ PENSS,PENSP,PENPP,XG,FPPPI,FPPSIG,FPS,FSS,             
     .           DPI,DSIG,DDD(300),FPPPIM,FDM                            
      COMMON /PIMCD/ MCD,NC,NOCC,NST,LPI(300),
     .  XP(100),YP(100),BETA(5050),CPI(10000), 
     .  CIPI(10000),IFROM(100),ITO(100),EP(100),POCC(100)          
      COMMON /SBIAS/ BIAS,IPI                                           
  105  FORMAT ((5X,16F8.4))                                              
  110  FORMAT (' SINGLET ENERGIES',/)                                    
  111  FORMAT (' TRIPLET ENERGIES',/)                                    
  112  FORMAT (//' CONFIGURATION INTERACTION:'//' IGCCI=',I2,            
     .  ' NB=',I2,' NAB=',I2,' ITRIP=',I2,' IDOB=',I2,' MCD=',I2,       
     .  '   EMAX=',F6.2/)                                               
  113 FORMAT (//1X,18A4)                                          
  130 FORMAT (6I3,F10.5)                                                
 2002 FORMAT (/' SINGLET MATRIX',/)                                      
 2003 FORMAT (/' TRIPLET MATRIX',/)                                      
 2004 FORMAT (//' SINGLET STATE EIGENVECTORS (VERTICAL)')                 
 2005 FORMAT (//' TRIPLET STATE EIGENVECTORS (VERTICAL)')                 
 9999 format(4f20.10)
      write(*,9)
    9 format(' Name of input file: ')
      read(*,11) NIN
      write(*,12)
   12 format(' Name of output file: ')
      read(*,11) NOUT
   11 format(a)
      open(5,file=nin)
      open(6,file=nout,status='new')
      open(1,file='f1',status='unknown',form='formatted')
      open(2,file='f2',status='unknown',form='formatted')
      open(3,file='f3',status='unknown',form='formatted')
      open(4,file='f4',status='unknown',form='formatted')
    1 CALL READER  
      N2=N+N                                                            
      N4=N2+N2                                                          
      NT=N+NH                                                           
      CALL GMATR(1)                                                     
      CALL SMATR                                                        
      CALL F0MATR                                                       
      CALL LOWDIN                                                                                                         
      WRITE(6,3010)                                                
 3010 FORMAT(2X/10X,'G-, S-, AND F0-MATRICES'/)             
      CALL SCF                                                                                                           
      WRITE(6,3020)                                             
 3020 FORMAT(2X/10X,'SCF-CALCULATION')                      
      DO 8 I=1,N5                                                       
    8 ET(I)=SUM(I)                                                      
c-----CI-CALCULATION;  READ PARAMETERS, ASSIGN DEFAULTS:                
   80 READ(5,130) IGCCI,NB,NAB,ITRIP,IDOB,MCD,  EMAX                    
      IF(IGCCI.LE.0) GOTO 1                                             
      DO 90 I=1,N5                                                      
   90 SUM(I)=ET(I)                                                      
      IF(NB.LE.0)  NB =10                                               
      IF(NAB.LE.0) NAB=10                                               
      IF(IDOB.EQ.0) IDOB=2                                              
      IF(NHOMO.GE.NLUMO) IDOB=-1                                        
      IF(MCD.EQ.0) MCD=1                                                
      MCD=MIN0(MCD,5)                                                   
      IF(IPI.LE.0) MCD=-1                                               
      IF(MCD.GT.0)IDOB=-1                                               
      IF(EMAX.LE..0) EMAX=15.0                                          
      WRITE(6,113) (NTIT(II),II=1,18)                                   
      WRITE(6,112) IGCCI,NB,NAB,ITRIP,IDOB,MCD,EMAX                     
      CALL GCCI                                                         
      IF(IPRNT.GE.2) WRITE(6,2002)                                                     
      IF(IPRNT.GE.2) CALL PRINT(S,JCOUNT,0,0,0,0,0)                                                                        
      WRITE(6,3040)                                                
 3040 FORMAT(2X/10X,'SINGLET(+TRIPLET) SCI-MATRIX')         
      NTEMP=N                                                           
      N=JCOUNT                                                          
      ICI=1                                                             
      CALL DIAGON                                                       
      WRITE(6,110)                                                      
      WRITE(6,105) (SUM(I),I=1,JCOUNT)                                  
      WRITE(6,2004)                                                     
      CALL PRINT(F,JCOUNT,JCOUNT,1,32,0,0)                                                                                                                             
 3050 FORMAT(2X/10X,'DIAG. SINGLET SCI-MATRIX')             
      CALL DMOM (SX,SY,SZ,V3,V4,V5,F,ZCORE,SUM,JCOUNT)                  
      CALL SPRT1(0,DUMMY)                                               
c-MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD-B
      IR=0                                                              
      DO 1000 J=1,NST                                                   
      EP(J)=SUM(J)                                                      
      DO 1000 I=1,NST                                                   
      IR=IR+1                                                           
 1000 CIPI(IR)=F(I,J) *  SQRT( POCC(IFROM(I)) - POCC(ITO(I)) )          
      CIPI(1)=F(1,1)                                                    
      CALL MCDB(1)                                                      
c-MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD-B
      N=NTEMP    
      JCOUNT=1                                                                                                      
      IF (ITRIP.LE.0) GO TO 10                                                                                             
      REWIND 3                                                          
      JCOUNT=JCOUNT-1                                                   
      JCO=JCOUNT-1                                                      
      N=JCOUNT                                                          
      DO 5 I=1,JCO                                                      
      IL=I+1                                                            
      DO 5 J=IL,JCOUNT                                                  
    5 READ(3,9999) S(J,I)                                                    
      DO 60 I=1,JCOUNT                                                  
   60 S(I,I)=EE(I+1)                                                    
      IF(IPRNT.GE.2) WRITE(6,2003)                                                     
      IF(IPRNT.GE.2) CALL PRINT(S,JCOUNT,0,0,0,0,0)                     
      CALL DIAGON                                                       
      N=NTEMP                                                           
      IF(IPRNT.GE.2) WRITE(6,111)                                                      
      IF(IPRNT.GE.2) WRITE(6,105) (SUM(I) ,I=1,JCOUNT)                                 
      WRITE(6,2005)                                                     
      CALL PRINT(F,JCOUNT,JCOUNT,1,32,0,0)                              
      WRITE(6,3100)                                              
 3100 FORMAT(2X/10X,'LOAD AND DIAG. OF TRIPLET SCI-MATRIX') 
   10 CONTINUE                                                         
c     GO TO 80                                                          
      END                                                               
c#####################################################################  
      SUBROUTINE SHIFT (A,N)                                            
      implicit real*8(A-H,O-Z)
      DIMENSION A(300,300)                                                
      DO 1 I=1,N                                                        
      DO 1 J=I,N                                                        
      AH=A(I,J)                                                         
      A(I,J)=A(J,I)                                                     
    1 A(J,I)=AH                                                         
      RETURN                                                            
      END                                                               
c#####################################################################  
      SUBROUTINE STORE(A,IUN,K)                                         
      implicit real*8(A-H,O-Z)
c-----STORE ARRAY A ON UNIT IUN    /JS-L                                
      DIMENSION A(300,300)                                                
      REWIND IUN                                                        
      WRITE(IUN,9999) ((A(I,J),J=1,K),I=1,K)                                 
 9999 format(4f20.10)
c     WRITE(6,100) K,K,IUN                                              
c 100 FORMAT(2X/10X,I2,' BY ',I2,' ARRAY STORED ON UNIT',I2/)           
      RETURN                                                            
      END                                                               
c#####################################################################  
      SUBROUTINE LOAD(A,IUN,K)                                          
      implicit real*8(A-H,O-Z)
c-----LOAD UNIT IUN INTO ARRAY A   /JS-L                                
      DIMENSION A(300,300)                                                
      REWIND IUN                                                        
      READ(IUN,9999) ((A(I,J),J=1,K),I=1,K)                                  
 9999 format(4f20.10)
c     WRITE(6,100) K,K,IUN                                              
c 100 FORMAT(2X/10X,I2,' BY ',I2,' ARRAY LOADED FROM UNIT',I2/)         
      RETURN                                                            
      END                                                               
c#####################################################################  
      SUBROUTINE PRINT(H,NI,MI,MINI,MAXI,K1,K2)                         
      implicit real*8(A-H,O-Z)
c     PROGRAMMED BY P.A.STRAUB AND H.BAUMANN                            
c-----MODIFIED BY JS-L                                                  
c      THIS SUBROUTINE PRINTS A NI*MI MATRIX                             
c     FOR MI=0 ONLY THE SYMMETRIC PART BELOW DIAGONAL                   
c     FOR NI=0 ONLY THE SYMMETRIC PART ABOVE DIAGONAL                   
      DIMENSION H(300,300)                                                
      DATA LINE/4H----/                                                 
  100 FORMAT(1X)                                                        
  103 FORMAT(1X,4H----,16(2A4))                                         
  101 FORMAT(/,16(5X,I3))                                             
  102 FORMAT(I3,1X,16F8.4)                                              
      MIN=MINI                                                          
      MAX=MAXI                                                          
      L1=K1                                                             
      L2=K2                                                             
      IF(K1.EQ.0) L1=100000                                             
      IF(K2.EQ.0) L2=1000                                               
      M=MIN0(MI,NI)                                                     
      N=MAX0(MI,NI)                                                     
      M1 = M                                                            
      IF(M1.LE.0) M1=N                                                  
      IF(MIN.LE.0) MIN=1                                                
      IF(MAX.GT.0) M1=MIN0(M1,MAX)                                      
      DO 11 J=MIN,M1,16                                                 
      JEND =MIN0(M1,J+15)                                               
      WRITE(6,101) (K,K=J,JEND)                                         
      WRITE(6,103) ((LINE,LL=1,2),K=J,JEND)                             
      IA=1                                                              
      IF(M.EQ.0) IA=J                                                   
      DO 10 I=IA,N                                                      
      IL2=I-1-L2                                                        
      IF(IL2.LT.0) GOTO 8                                               
      IF(FLOAT(IL2)/FLOAT(L1)-FLOAT(IL2/L1).LT..000001) WRITE(6,100)    
    8 JIEND = JEND                                                      
      IF (M.NE.0) GO TO 9                                               
      IF (JEND.GE.I) JIEND=I                                            
      IF (NI.EQ.0) GO TO 12                                             
    9 WRITE(6,102) I,(H(I,K),K=J,JIEND)                                 
      GO TO 10                                                          
   12 WRITE(6,102) I,(H(K,I),K=J,JIEND)                                 
   10 CONTINUE                                                          
   11 CONTINUE                                                          
      RETURN                                                            
      END                                                               

c#####################################################################  
      SUBROUTINE DMOM (SX,SY,SZ,SCIX,SCIY,SCIZ,U,O,E,JC)                
      implicit real*8(A-H,O-Z)
c-----PRINTS EXC. ENERGIES, TRANSITION MOMENTS AND OSC. STRENGTHS /JS-L 
      COMMON /D/ XJSL(5,300),IQ(2,300)                                    
      COMMON /OCCDAT/ NHOMO,NLUMO,NPE,NPMO,OCC(300)                      
      DIMENSION SX(300),SY(300),SZ(300),SCIX(300),SCIY(300),
     1      SCIZ(300),O(300),U(300,300),E(300)                                              
  101 FORMAT(1X,102(1H-))                                               
  103 FORMAT (I5,F12.6,F12.4,F12.1,5F12.6)                              
      IF(ABS(U(1,1)).LT.SQRT(.5).AND.E(1).LT..0) WRITE(6,104)           
  104 FORMAT(//' ----CAUTION----'//' SCF GROUND STATE UNSTABLE'//       
     .     ' TRANSITION MOMENTS AND OSC. STRENGTHS UNRELIABLE'//)       
      WRITE(6,102)                                                      
  102 FORMAT (/,3X,2HNO,7X,2HEV,9X,2HKK,12X,2HNM,9X,2HSX,10X,2HSY,      
     .        10X,2HSZ,10X,3HOSC,7X,6HLOGEPS)
      WRITE(6,101)                           
      IF(NHOMO.GE.NLUMO) WRITE(6,105)                                   
  105 FORMAT(105X,'HOLE STATE CHARACTER')                           
      DO 10 J=2,JC                                                                                                 
      E(J)=E(J)-E(1)                                                    
      SCIX(J)=0.0                                                       
      SCIY(J)=0.0                                                       
      SCIZ(J)=0.0                                                       
      HOLE=0.0                                                          
      DO 5 I=2,JC                                                       
      UIJ=U(I,J)                                                        
      SCIX(J)=SCIX(J)+SX(I)*UIJ                                         
      SCIY(J)=SCIY(J)+SY(I)*UIJ                                         
      SCIZ(J)=SCIZ(J)+SZ(I)*UIJ                                         
      IF(IQ(2,I).LE.NHOMO) HOLE=HOLE+UIJ*UIJ                            
    5 CONTINUE                                                          
      EKK = E(J)*8.066                                                  
      ENM = 10000./EKK                                                  
      SCI2=SCIX(J)*SCIX(J)+SCIY(J)*SCIY(J)+SCIZ(J)*SCIZ(J)              
      O(J)=.00379277*E(J)*SCI2                                          
      ELOG=ABS(O(J))*5.43658E+4                                         
      IF (ELOG.LE.1.) ELOG=0.                                           
      IF (ELOG.NE.0.) ELOG=ALOG10(ELOG)                                 
      WRITE(6,103) J,E(J),EKK,ENM,SCIX(J),SCIY(J),SCIZ(J),O(J),ELOG     
      IF(NHOMO.GE.NLUMO) WRITE(6,106) HOLE,J                            
  106 FORMAT(1H+,110X,F8.5,I5)                                          
   10 CONTINUE                                                          
      RETURN                                                            
      END                                                               
c#####################################################################  
      FUNCTION H(LX,IX,JX,MX)                                           
      implicit real*8(A-H,O-Z)
c     INTEGRAL (LX,IX/JX,MX) IS COMPUTED BY THE GAMMA INTEGRALS STORED  
c     IN THE LOWER PART OF G AND THE VECTORS STORED IN U                
c     PROGRAMMED BY H. BAUMANN                                          
c-----MODIFIED BY JS-L                                                  
      COMMON /B/ ES(300),ER(300),ZCORE(300),ET(300),SUM(300),
     1           C(400),GAX(300)
      COMMON /E/ NCYCL,EPREV,ITERR,N5                                   
      COMMON /C/ G(90000),U(90000)                                      
      DIMENSION W(300)                                                   
      N4=N5-1                                                           
      I1=300*(LX-1)                                                      
      J1=300*(JX-1)                                                      
      K1=300*(IX-1)                                                      
      L1=300*(MX-1)                                                      
      DO 93 NZ=1,N5                                                     
      ES(NZ)=0.                                                         
   93 W(NZ)=U(I1+NZ)*U(K1+NZ)                                           
      H=W(N5)*GAX(N5)*U(J1+N5)*U(L1+N5)                                 
      NZ1=-300                                                           
      DO 90 NZ=1,N4                                                     
      NZ1=NZ1+300                                                        
      WN=W(NZ)                                                          
      MI=NZ+1                                                           
      H1=WN*GAX(NZ)                                                     
      DO 94 MZ=MI,N5                                                    
      ES(MZ)=ES(MZ)+WN*G(NZ1+MZ)                                        
   94 H1=W(MZ)*G(NZ1+MZ)+H1                                             
   90 H=H+U(J1+NZ)*U(L1+NZ)*(H1+ES(NZ))                                 
      H=H+U(J1+N5)*U(L1+N5)*ES(N5)                                      
      RETURN                                                            
      END                                                               
c#####################################################################  
      SUBROUTINE READER                                                 
      implicit real*8(A-H,O-Z)
c      READ INPUT DATA                                                   
c     PROGRAMMED BY H. BAUMANN                                          
c-----NUMEROUS MODIFICATIONS BY JS-L                                    
      DIMENSION SDAT(36),ASDAT(36),APDAT(36),BDAT(36),ZDAT(36),GA(36),  
     . X(100),Y(100),Z(100),ZS(100),IQ(300),BSDAT(36),BPDAT(36),             
     . USDAT(36),UPDAT(36)                                              
      COMMON /A/ N,NH,NE,NQ,ECORE,ICOUNT,JCOUNT,EMAX,IPRNT,IGCCI,       
     .           EATOM,ISCF,ITGT,ICI,GAH,ITRIP,IPEN,ETOT,ETOT1,IDOB     
      COMMON /B/ E(300),EE(300),ZCORE(300),ET(300),SUM(300),C(400),          
     .           GAX(300),NTIT(20)                                       
      COMMON /E/ NCYCL,EPREV,ITERR,N5,NTOTAL,N2S,N2P,NTN,NG,DAMP        
      COMMON /SPTERM/ SP(100)                                            
      COMMON /OCCDAT/ NHOMO,NLUMO,NPE,NPMO,OCC(300)                      
      COMMON /F/ PENSS,PENSP,PENPP,XG,FPPPI,FPPSIG,FPS,FSS,             
     .           DPI,DSIG,DDD(300),FPPPIM,FDM                            
      COMMON /PIMCD/ MCD,NC,NOCC,NST,LPI(300),
     .  XP(100),YP(100),BETA(5050),CPI(10000),CIPI(10000),
     .  IFROM(100),ITO(100),EP(100),POCC(100)          
      COMMON /SBIAS/ BIAS,IPI       
      COMMON /NMOP/ NMOP                                    
      EQUIVALENCE (C(1),X(1)),(C(101),Y(1)),(C(201),Z(1)),
     . 		(C(301),ZS(1)),(SUM(1),IQ(1))                                        
      DATA NQS /300/,IBLANK/4H    /                                      
      DATA SDAT(1)  ,ASDAT(1)         ,BDAT(1),ZDAT(1),GA(1)            
     .    /1.2,   -7.175,        -12. ,1.   ,12.85/                     
c -----ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF   
      DATA SDAT(2),ASDAT(2),APDAT(2),BDAT(2),ZDAT(2),GA(2)              
     .   /  10.0  ,  99999., 99999. ,  0.01 , 0.000 , 12.85/            
c-----ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF   
      DATA SDAT(3),ASDAT(3),APDAT(3),BDAT(3),ZDAT(3),GA(3)              
     .   /.65,   -3.105, -2.05,  -3.,  1.,  2.98/                       
      DATA SDAT(4),ASDAT(4),APDAT(4),BDAT(4),ZDAT(4),GA(4)              
     1    /.975,-6.55,-3.035,-4.,2.,5.85/                               
      DATA SDAT(5),ASDAT(5),APDAT(5),BDAT(5),ZDAT(5),GA(5)              
     1    /1.3,-10.305,-4.37,-5.,3.,8.10/                               
      DATA SDAT(6),ASDAT(6) ,APDAT(6),BDAT(6),ZDAT(6),GA(6)             
     1    /1.625,-14.960,-5.805,-17.5,4.   ,10.93/                      
      DATA SDAT(7),ASDAT(7) ,APDAT(7),BDAT(7),ZDAT(7),GA(7)             
     1    /1.950,-20.485,-8.480,-26.0,5.   ,13.10/                      
      DATA SDAT(8),ASDAT(8) ,APDAT(8),BDAT(8),ZDAT(8),GA(8)             
     1    /2.275,-27.255,-10.965,-45. ,6.   ,15.27/                     
      DATA SDAT(9),ASDAT(9) ,APDAT(9),BDAT(9),ZDAT(9),GA(9)             
     1    /2.600,-28.48,-12.18,-50.,7.   ,17.36/                        
      DATA SDAT(11),ASDAT(11),APDAT(11),BDAT(11),ZDAT(11),GA(11)        
     1    /1.1,-2.805,-1.565,-.5,1.,2.95/                               
      DATA SDAT(12),ASDAT(12),APDAT(12),BDAT(12),ZDAT(12),GA(12)        
     1    /1.425,-5.875,-2.29,-1.,2.,4.46/                              
      DATA SDAT(13),ASDAT(13),APDAT(13),BDAT(13),ZDAT(13),GA(13)        
     1    /  1.75,-8.595,-3.92,-1.5,3.,5.1/                             
      DATA SDAT(14),ASDAT(14),APDAT(14),BDAT(14),ZDAT(14),GA(14)        
     1    /2.075,-12.125,-6.005,-5.25,4.,6.37/                          
      DATA SDAT(15),ASDAT(15),APDAT(15),BDAT(15),ZDAT(15),GA(15)        
     1    /2.4,-14.34,-7.235,-7.8,5.,9.31/                              
      DATA SDAT(16),ASDAT(16),APDAT(16),BDAT(16),ZDAT(16),GA(16)        
     1    /1.817,-15.81,-7.385,-18.0,6.,10.01/                          
      DATA SDAT(17),ASDAT(17),APDAT(17),BDAT(17),ZDAT(17),GA(17)        
     1    /3.05,-17.5,-9.38,-15.,7.,11.3/                               
      DATA SDAT(34),ASDAT(34),APDAT(34),BDAT(34),ZDAT(34),GA(34)        
     1    /1.565,-15.425,-7.10,-14.0,6.,9.16/                           
   99  FORMAT (18(A4))                                                   
  100 FORMAT (1X,18(A4),//)                                                
  101 FORMAT (6I3,12F5.2)                                               
  104 FORMAT (I2,F13.6,2F15.6,5F5.2)                                    
  105 FORMAT(
     .   '                    >-----------<'/
     .   '                      L C O A O  '/
     .   '                    >-----------<'//              
     .   '"Linear Combination of Orthogonalized Atomic Orbitals"'/
     .   '  J. Spanget-Larsen: Theor. Chem. Acc. 59, 137 (1997) '//
     .   '       100 atoms - 300 orbitals developement-vs.'/
     .   '                       Dec. 2005'/   
     .   '------------------------------------------------BKVH--'//)
  106 FORMAT (' N=',I2,'  NH=',I2,'  NE=',I3,                           
     1      '  NPE=',I2,'  NPMO=',I2,'  IPI=',I2, 
     2       '  IPRNT=',I2,'  NMOP=',I2)
  107 FORMAT(1X/' PP(PI)=',F6.3,'   PP(SIG)=',F6.3,'   PS=',F6.3,       
     .       '   SS=',F6.3,'  D(PI)=',F6.3,'  D(SIG)=',F6.3/           
     .' XG=',F6.3,'  PEN(SS)=',F6.3,' PEN(SP)=',F6.3,' PEN(PP)=',F6.3/
     .' PRM. FOR MAGN. DIPOLE INTEGRALS:  PP(PI)=',F6.3,'  FD=',F6.3/
     .' SIGMA ORBITAL ''BIAS'':  ',F6.3,'    ZDAT(2): ',F6.3)           
  698 FORMAT (I5,3F15.6,I11,2F10.3,10X,F10.3,10X,1F10.3)                
  699 FORMAT (I5,3F15.6,I11,6F10.3)                                     
 1000 FORMAT(1H )                                                       
 7499 FORMAT(//4X,1HI,11X,1HX,14X,1HY,14X,1HZ,11X,4HN(Z),5X,3HEXP,6X,  
     .    4HX(S),6X,4HX(P),5X,7HBETA(S),3X,7HBETA(P),3X,5HGAMMA/        
     .    1X,130(1H-))                                                  
      NQ=NQS                                                            
    1 READ(5,99) (NTIT(I),I=1,18)                                       
      WRITE(6,1000)                                                     
      WRITE(6,105)                                                      
      WRITE(6,100) (NTIT(I),I=1,18)                                     
      IF (NTIT(1).EQ.IBLANK) STOP                                       
      READ(5,101) NE,NPE,NPMO,IPI,IPRNT,NMOP,FPPPI,FPPSIG,FPS,FSS,           
     .            DPI,DSIG,XG,PENSS,PENSP,PENPP,BIAS,DAMP               
c-----DEFAULTS                                                          
                        ISCF   =  30                                    
      IF(DAMP.LE..01)   DAMP   =  0.5                                   
      IF(PENSS.LE..0)   PENSS  =  0.5                                   
      IF(PENSP.LE..0)   PENSP  =  0.5                                   
      IF(PENPP.LE..0)   PENPP  =  1.5                                   
      IF(XG.LE..0)      XG     =  1.0                                   
      IF(DPI.EQ..0)     DPI    =  0.23                                  
      IF(DSIG.EQ..0)    DSIG   =  0.23                                  
      IF(FPPPI.LE..0)   FPPPI  =  0.75                                  
      IF(FPPSIG.LE..0)  FPPSIG =  1.00                                  
      IF(FPS.LE..0)     FPS    =  1.00                                  
      IF(FSS.LE..0)     FSS    =  1.00  
      IF(NMOP.EQ.0)     NMOP   =  32                                
                        FDM    =  0.0                                   
                        FPPPIM =  1.50/EXP(-FDM*1.397*1.625/.529167)    
      DO 990 I=1,36                                                     
      BSDAT(I)=ASDAT(I)                                                 
  990 BPDAT(I)=APDAT(I)                                                 
      BSDAT(2)=0.001                                                    
      BPDAT(2)=0.001                                                    
      N=0                                                               
      NH=0                                                              
      DO 13 I=1,1000                                                    
c-----ASD, APD, BSD, BPD, SD ARE NEW ATOMIC PARAMETERS    /JS-L         
      READ(5,104) IQ(I),X(I),Y(I),Z(I),ASD,APD,BSD,BPD,SD               
      IQT=IQ(I)                                                         
      IF (IQT.LT.1) GO TO 14                                            
c-----ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF   
      IF (IQT.EQ.2) ZDAT(2) = ASD                                       
      IF (IQT.EQ.2) GOTO 8888                                           
c-----ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF   
      IF(ABS(ASD).GT..0) ASDAT(IQT)=-ABS(ASD)                           
      IF(ABS(APD).GT..0) APDAT(IQT)=-ABS(APD)                           
 8888 IF(ABS(BSD).GT..0) BSDAT(IQT)=-ABS(BSD)                           
      IF(ABS(BPD).GT..0) BPDAT(IQT)=-ABS(BPD)                           
      IF(ABS(SD).GT..0)   SDAT(IQT)= ABS(SD)                            
      IF (IQT.LE.2) NH=NH+1                                             
      IF (IQT.GT.2) N=N+1                                               
   13 CONTINUE                                                          
   14 N2=N+N                                                            
      N4=N2+N2                                                          
      N5=N4+NH                                                          
      NTOTAL=NH+N                                                       
      DO 15 I=1,NTOTAL                                                  
      IQT=IQ(I)                                                         
      DO 16 J=I,NTOTAL                                                  
      IF (IQ(J).GE.IQT) GO TO 16                                        
      XT=X(I)                                                           
      YT=Y(I)                                                           
      ZT=Z(I)                                                           
      IQ(I)=IQ(J)                                                       
      X(I)=X(J)                                                         
      Y(I)=Y(J)                                                         
      Z(I)=Z(J)                                                         
      IQ(J)=IQT                                                         
      X(J)=XT                                                           
      Y(J)=YT                                                           
      Z(J)=ZT                                                           
      IQT=IQ(I)                                                         
   16 CONTINUE                                                          
   15 CONTINUE                                                          
c-MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD- 
      NPI=NH+1                                                          
      NAT=NH+N                                                          
      IR=0                                                              
      DO 2000 IA=NPI,NAT                                                
      IR=IR+1                                                           
      XP(IR)=X(IA)                                                      
 2000 YP(IR)=Y(IA)                                                      
c-MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD- 
      N2S=NH+1                                                          
      N2P=NTOTAL+1                                                      
      NTN=NTOTAL+N                                                      
      N4=0                                                              
      TWOSQ3=2.*SQRT(3.)                                                
      DO 52 I1=1,NTOTAL                                                 
      I2=I1+N                                                           
      I3=I2+N                                                           
      I4=I3+N                                                           
      IKI=IQ(I1)                                                        
c-----SP-TERMS FOR TRANSITION MOMENTS (SEE TRAMOM)      /JS-L           
      SP(I1)=.0                                                         
      IF(IKI.LE.2) GO TO 999                                            
      FN=5.                                                             
      IF(IKI.GT.10) FN=7.                                               
      SP(I1)=FN/(SDAT(IKI)*TWOSQ3)                                      
  999 ET(I1)=BSDAT(IKI)                                                 
      DDD(I1)=ASDAT(IKI)                                                
      ZCORE(I1)=ZDAT(IKI)                                               
c-----ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF   
      NESF = INT(ZDAT(IKI))                                             
      IF(IKI.EQ.2) NESF = 0                                             
c     N4=N4+INT(ZCORE(I1))               ( ORIGINAL VERSION )           
      N4=N4+NESF                                                        
c-----ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF---ESF   
      ZS(I1)=SDAT(IKI)                                                  
      GAX(I1)=GA(IKI)                                                   
      IF (IKI.GT.2) GOTO 51                                             
      GO TO 52                                                          
   51 DDD(I2)=APDAT(IKI)                                                
      DDD(I3)=APDAT(IKI)                                                
      DDD(I4)=APDAT(IKI)                                                
      ET(I2)=BPDAT(IKI)                                                 
      ET(I3)=BPDAT(IKI)                                                 
      ET(I4)=BPDAT(IKI)                                                 
      GAX(I2)=GA(IKI)                                                   
      GAX(I3)=GA(IKI)                                                   
      GAX(I4)=GA(IKI)                                                   
   52 CONTINUE                                                          
      GAH=GA(1)                                                         
      NE=N4-NE                                                          
      WRITE(6,106) N,NH,NE,NPE,NPMO,IPI,IPRNT,NMOP                           
      WRITE(6,107) FPPPI,FPPSIG,FPS,FSS,DPI,DSIG,                       
     .             XG,PENSS,PENSP,PENPP,FPPPIM,FDM,BIAS                 
     .             ,ZDAT(2)                                             
      WRITE(6,7499)                                                     
      IF(NH.EQ.0) GO TO 33                                              
      DO 30 I=1,NH                                                      
   30 WRITE(6,698) I,X(I),Y(I),Z(I),IQ(I),ZS(I),DDD(I),ET(I),GAX(I)     
   33 IF(N.EQ.0) GO TO 34                                               
      DO 50 I=1,N                                                       
      J=NH+I                                                            
      JJ=J+N                                                            
      IQJ=IQ(J)                                                         
   50 WRITE(6,699) J,X(J),Y(J),Z(J),IQJ,ZS(J),DDD(J),DDD(JJ),           
     .             ET(J),ET(JJ),GAX(J)                                  
   34 RETURN                                                            
      END                                                               
c#####################################################################
      SUBROUTINE AUFBAU(NE,NMO,E,LOOK)                                  
      implicit real*8(A-H,O-Z)
      DIMENSION E(300)                                                   
      COMMON /OCCDAT/ NHOMO,NLUMO,NPE,NPMO,OCC(300)                      
      IF(LOOK) 110,50,1                                                 
c-----DETERMINE OPEN SHELL DEGENERACY (NPMO) AND TOTAL OCCUPANCY (NPE)  
    1 NPMO=0                                                            
      NPE=0                                                             
      NFMO=NE/2                                                         
      I=(NE+1)/2                                                        
   10 IF(I.GE.NMO) GOTO 20                                              
      IF(ABS(E(I+1)-E(I)).GT..0001) GOTO 20                             
      I=I+1                                                             
      GOTO 10                                                           
   20 IF(NE.EQ.2*I) GOTO 60                                             
      J=(NE+2)/2                                                        
   30 IF(J.EQ.1) GOTO 40                                                
      IF(ABS(E(J)-E(J-1)).GT..0001) GOTO 40                             
      J=J-1                                                             
      GOTO 30                                                           
   40 NFMO=J-1                                                          
      NPMO=I-NFMO                                                       
      NPE=NE-2*NFMO                                                     
      GOTO 60                                                           
   50 LOOK=-1                                                           
      NFMO=(NE-NPE)/2                                                   
      IF(2*NFMO.NE.NE-NPE) NPE=NPE+1                                    
      NPMO=MAX0(NPMO,(NPE+1)/2)                                         
c-----ASSIGN GRAND CANONICALLY AVERAGED OCCUPATION NUMBERS              
   60 NFMO1=NFMO+1                                                      
      DO 70 I=NFMO1,NMO                                                 
   70 OCC(I)=0.0                                                        
      DO 80 I=1,NFMO                                                    
   80 OCC(I)=1.0                                                        
      IF(NPMO.EQ.0) GOTO 100                                            
      OCCP=0.5*FLOAT(NPE)/FLOAT(NPMO)                                   
      DO 90 I=1,NPMO                                                    
   90 OCC(NFMO+I)=OCCP                                                  
  100 NHOMO=NFMO+NPMO                                                   
      NLUMO=NFMO+1                                                      
  110 RETURN                                                            
      END                                                               
c#######################################################################
      SUBROUTINE SPRT1(IFLAG,BTERM)                                     
      implicit real*8(A-H,O-Z)
c     A PRINTPLOT OF THE ELECTRONIC TRANSITIONS IS PRODUCED             
c     PROGRAMMED BY J.KELEMEN UND H.BAUMANN                             
c-----MODIFIED TO INDICATE POL. DIR.S AND MCD B-SIGNS BY JS-L           
      COMMON /A/ NT                                                     
      COMMON /B/ ET(300),EE(300),OS(300),SO(300),E(300),C(400),
     1           GX(300),LT(20) 
      COMMON /D/ V1(300),V2(300),SCIX(300),SCIY(300),SCIZ(300)               
      DIMENSION IPLT(100),BTERM(1)                                      
      DATA HX/.5/,IBR/1H-/,IPL/1H+/,IPNT/1H:/,IBLNK/1H /                
      DATA IX/1HX/,IY/1HY/,IZ/1HZ/,I0/1H0/                              
    2 FORMAT (1H ,49X,21HL O G   E P S I L O N )                        
    3 FORMAT (1H , 8X,3H1.0,7X,3H1.5, 7X,3H2.0, 7X,3H2.5, 7X, 3H3.0,7X, 
     .3H3.5, 7X, 3H4.0, 7X,3H4.5, 7X, 3H5.0, 7X, 3H5.5, 7X,3H6.0)       
    4 FORMAT (1H ,4X, 6HKK   :, 10(10H . . . . :),4X, 2HNM )            
    5 FORMAT (1H ,9X,101(1H-))                                          
    6 FORMAT (1H /)                                                     
    7 FORMAT (1H )                                                      
    8 FORMAT (1H ,F7.3,3H  :,99A1,1H:,F7.0)                             
    9 FORMAT (1H )                                                      
      N=NT-1                                                            
      WRITE(6,9)                                                        
      WRITE(6,1000) (LT(I),I=1,18)                                      
 1000 FORMAT(//11X,18(A4))                                           
      EP=1000.                                                          
      OMAX=0.                                                           
      DO 1 J=2,N                                                        
      E(J)=ABS(E(J))*8.066                                              
      IF (OS(J).LE.1.E-5) OS(J)=1.E-5                                   
      IF(IFLAG.NE.1) OS(J)=ALOG10(ABS(OS(J))*5.43658E+4)                
      EP=AMIN1(E(J),EP)                                                 
      IF (E(J).LE.100.) OMAX=AMAX1(OS(J),OMAX)                           
    1 CONTINUE                                                          
      SP=AMAX1(EP,57.)                                                  
      EP=AINT(EP)-5.                                                    
      IF (EP.LE.0.) EP=1.                                               
      SP=AINT(SP)+5.                                                    
      LIM=(SP-EP)/HX                                                    
      WW=.36*BL**2                                                      
      WRITE(6,6)                                                        
      WRITE(6,2)                                                        
      WRITE(6,7)                                                        
      WRITE(6,3)                                                        
      WRITE(6,4)                                                        
      WRITE(6,5)                                                        
      WN=SP+HX                                                          
      DO 10 LL=1,LIM                                                    
      WN=WN-HX                                                          
      A=WN-0.5*HX                                                       
      B=WN+0.5*HX                                                       
      WL=(1.0E+4)/WN                                                    
      DO 20 K=1,100                                                     
   20 IPLT(K)=IBLNK                                                     
      IPLT(20)=IPNT                                                     
      IPLT(40)=IPNT                                                     
      IPLT(60)=IPNT                                                     
      IPLT(80)=IPNT                                                     
      DO 30 J=1,N                                                       
      IF (E(J).LT.A.OR.E(J).GE.B) GO TO 30                              
      R=AMAX1(AMAX1(ABS(SCIX(J)),ABS(SCIY(J))),ABS(SCIZ(J)))            
      IR=IX                                                             
      IF(ABS(SCIY(J)).GE.R) IR=IY                                       
      IF(ABS(SCIZ(J)).GE.R) IR=IZ                                       
      IF(R.LT..00001) IR=I0                                             
      M=1                                                               
      IF (OS(J).GT.1.) M=(OS(J)-1.)*20.                                 
      IF(M.GT.100) M=100                                                
      DO 40 K=1,M                                                       
      ISYM=IBR                                                          
      IF( IFLAG.EQ.1 .AND. BTERM(J).GT..0 ) ISYM=IPL                    
      IP=IPLT(K)                                                        
      IF(IP.NE.IX.AND.IP.NE.IY.AND.IP.NE.IZ.AND.IP.NE.I0) IPLT(K)=ISYM  
   40 CONTINUE                                                          
      IPLT(M)=IR                                                        
   30 CONTINUE                                                          
      WRITE(6,8) WN,(IPLT(J),J=1,99),WL                                 
   10 CONTINUE                                                          
      WRITE(6,5)                                                        
      WRITE(6,4)                                                        
      WRITE(6,3)                                                        
      WRITE(6,7)                                                        
      WRITE(6,2)                                                        
      WRITE(6,9)                                                        
      DO 50 J=1,N                                                       
      E(J)=E(J)/8.066                                                   
   50 CONTINUE                                                          
      RETURN                                                            
      END                                                               
c#####################################################################  
      SUBROUTINE GMATR(IFLAG)                                           
      implicit real*8(A-H,O-Z)
c      COMPUTES THE GAMMA INTEGRAL PARAMETERS                            
c-----PROGRAMMED BY H. BAUMANN, STRONGLY MODIFIED BY JS-L               
      DIMENSION X(100),Y(100),Z(100),ZS(100)                                
      COMMON /A/ N,NH,NE,NQ,ECORE,ICOUNT,JCOUNT,EMAX,IPRNT,             
     1           IGCCI,EATOM,ISCF,ITGT,ICI,GAH                               
      COMMON /B/ E(300),EE(300),ZCORE(300),ET(300),SUM(300),
     1           C(400),GAX(300)   
      COMMON /C/ S(300,300),G(300,300)                                      
      COMMON /E/ NCYCL,EPREV,ITERR,N5,NTOTAL,N2S,N2P,NTN,NG             
      COMMON /F/ PENSS,PENSP,PENPP,XG,FPPPI,FPPSIG,FPS,FSS,             
     .           DPI,DSIG,DDD(300),FPPPIM,FDM                            
      EQUIVALENCE (C(1),X(1)),(C(101),Y(1)),(C(201),Z(1)),
     .            (C(301),ZS(1)) 
  101 FORMAT (//' DIATOMIC ELECTRON REPULSION INTEGRALS (EV)')            
  102 FORMAT (//' INTERATOMIC DISTANCES (ANGSTROM)')                      
c     GAMMA-MATRIX IN UPPER HALF OF G                                   
c     DISTANCE-MATRIX IN LOWER HALF OF G                                
      XG=ABS(XG)                                                        
      XGINV=-1./XG                                                      
      DO 2 I=1,NTOTAL                                                   
      DO 2 J=1,NTOTAL                                                   
    2 G(I,J)=0.                                                         
      IF (NH.LT.1) GO TO 8                                              
      DO 7 I=1,NH                                                       
      DO 5 J=I,NH                                                       
      IF (I.EQ.J) GO TO 5                                               
      D=SQRT((X(I)-X(J))**2+(Y(I)-Y(J))**2+(Z(I)-Z(J))**2)              
      G(J,I)=D                                                          
      G(I,J)=14.3942*(D**XG+(14.3942/GAH)**XG)**XGINV                   
    5 CONTINUE                                                          
      IF (N.LT.1) GO TO 7                                               
      DO 6 J=1,N                                                        
      L=J+NH                                                            
      D=SQRT((X(I)-X(L))**2+(Y(I)-Y(L))**2+(Z(I)-Z(L))**2)              
      G(L,I)=D                                                          
      GG = 1./GAX(L) + 1./GAH                                           
      G(I,L)=14.3942 * (D**XG+(7.1971*GG)**XG) **XGINV                  
    6 CONTINUE                                                          
    7 CONTINUE                                                          
    8 IF (N.LT.1) GO TO 18                                              
      DO 10 I=1,N                                                       
      II=I+NH                                                           
      DO 10 J=I,N                                                       
      JJ=J+NH                                                           
      IF (I.EQ.J) GO TO 10                                              
      D=SQRT((X(II)-X(JJ))**2+(Y(II)-Y(JJ))**2+(Z(II)-Z(JJ))**2)        
      G(JJ,II)=D                                                        
      GG = 1./GAX(II) + 1./GAX(JJ)                                      
      G(II,JJ)=14.3942 * (D**XG+(7.1971*GG)**XG) **XGINV                
   10 CONTINUE                                                          
      IF(IFLAG.EQ.1) WRITE(6,102)                                       
      IF(IFLAG.EQ.1) CALL PRINT(G,NTOTAL,0,0,0,0,0)                     
      DO 4 I=1,NTOTAL                                                   
    4 G(I,I)=GAX(I)                                                     
      IF(IPRNT.LT.1) GOTO 24                                            
      WRITE(6,101)                                                      
      CALL PRINT(G,0,NTOTAL,0,0,0,0)                                    
   24 DO 11 J=N2S,NTOTAL                                                
      DO 11 L=2,4                                                       
      I=J+N*(L-1)                                                       
      G(I,I)=G(J,J)                                                     
   11 CONTINUE                                                          
      DO 14 I=1,NH                                                      
      DO 13 J=N2S,NTOTAL                                                
      DO 12 JJ=2,4                                                      
      L=J+N*(JJ-1)                                                      
   12 G(I,L)=G(I,J)                                                     
   13 CONTINUE                                                          
   14 CONTINUE                                                          
      DO 17 K=N2S,NTOTAL                                                
      DO 17 M=N2S,NTOTAL                                                
      IG=MAX0(K,M)                                                      
      IK=MIN0(K,M)                                                      
      DO 16 II=1,4                                                      
      I=IK+N*(II-1)                                                     
      DO 15 JJ=1,4                                                      
      J=IG+N*(JJ-1)                                                     
      IMA=MAX0(I,J)                                                     
      IMI=MIN0(I,J)                                                     
      IF (J.NE.I) G(IMI,IMA)=G(IK,IG)                                   
   15 CONTINUE                                                          
   16 CONTINUE                                                          
   17 CONTINUE                                                          
   18 ECORE=0.                                                          
      NH1=NH-1                                                          
      DO 20 I=1,NH1                                                     
      I1=I+1                                                            
      DO 20 J=I1,NH                                                     
      ECORE=ECORE+ZCORE(I)*ZCORE(J)*G(I,J)                              
   20 CONTINUE                                                          
      IF (N.LT.1) GO TO 23                                              
      DO 21 I=1,NH                                                      
      DO 21 JJ=N2S,NTOTAL                                               
      ECORE=ECORE+ZCORE(I)*ZCORE(JJ)*G(I,JJ)                            
   21 CONTINUE                                                          
      IF (N.LE.1) GO TO 23                                              
      N1=NTOTAL-1                                                       
      DO 22 II=N2S,N1                                                   
      I1=II+1                                                           
      DO 22 JJ=I1,NTOTAL                                                
      ECORE=ECORE+ZCORE(II)*ZCORE(JJ)*G(II,JJ)                          
   22 CONTINUE                                                          
   23 CONTINUE                                                          
      CALL STORE(G,2,N5)                                                
      RETURN                                                            
      END                                                               
c#####################################################################  
      SUBROUTINE SMATR                                                  
      implicit real*8(A-H,O-Z)
c-----COMPUTE OVERLAP MATRIX, SCALED OVERLAP MATRICES, ETC.     /JS-L   
c     BASED ON:    QCPE 228, M.J.S. DEWAR AND P.J. STUDENT              
      COMMON/A/ NX,NH,NE,NQ,ECORE,ICOUNT,JCOUNT,EMAX,IPRNT,             
     .     IGCCI,EATOM,ISCF,ITGT                                        
      COMMON /B/ EK(300),EE(300),ZCORE(300),ET(300),SUM(300),C(400)          
      COMMON /C/ S(300,300),U(300,300)                                      
      COMMON /E/ NCYCL,EPREV,ITERR,N5,NUMAT,N2S                         
      COMMON /F/ PENSS,PENSP,PENPP,XG,FPPPI,FPPSIG,FPS,FSS,             
     .           DPI,DSIG,DDD(300),FPPPIM,FDM                            
      COMMON /PIMCD/ MCD,NC,NOCC,NST,LPI(300),XP(100),YP(100),BET(5050),   
     .  BETAM(100,100),CIPI(10000),IFROM(100),ITO(100),EP(100),POCC(100)       
      DIMENSION ZA(100),Z(5),A(7),B(7),T(3),X(100),
     .    Y(100),Z1(100),IQ(300)    
      EQUIVALENCE (C(1),X(1)),(C(101),Y(1)),(C(201),Z1(1)),
     . 		(C(301),ZA(1)),(IQ(1),SUM(1))                                        
  101 FORMAT (//'OVERLAP MATRIX  ')                                      
      DO 9999 I=1,N5                                                    
      DO 9997 J=1,I                                                     
      U(J,I)=0.0                                                        
      S(I,J)=0.0                                                        
 9997 S(J,I)=0.0                                                        
 9999 S(I,I)=1.0                                                        
      DO 16 I1=1,NUMAT                                                  
      NA=IQ(I1)                                                         
c-----NEXT STATEMENT CHANGED FROM  IPLUS=NS2  TO AVOID REPETITION /JS-L 
      IPLUS=MIN0(NUMAT,I1+1)                                            
      DO 15 J1=IPLUS,NUMAT                                              
      IF (I1.EQ.J1) GO TO 15                                            
      IG=MAX0(I1,J1)                                                    
      IK=MIN0(I1,J1)                                                    
      NB=IQ(J1)                                                         
c     II=1 FIRST-FIRST                                                  
c     II=2 FIRST-SECOND                                                 
c     II=3 FIRST-THIRD                                                  
c     II=4 SECOND-SECOND                                                
c     II=5 SECOND-THIRD                                                 
c     II=6 THIRD-THIRD                                                  
      II=1                                                              
      IF ((NA.LT.3).AND.(NB.LT.3)) GO TO 14                             
      II=6                                                              
      IF ((NA.GT.10).AND.(NB.GT.10)) GO TO 14                           
      IF ((NA.GT.10).OR.(NB.GT.10)) GO TO 13                            
      II=4                                                              
      IF ((NA.LT.3).OR.(NB.LT.3)) II=2                                  
      GO TO 14                                                          
   13 II=3                                                              
c-----NEXT STATEMENT INSERTED BY JS-L                                   
      IF ((NA.LT.3).OR.(NB.LT.3)) GOTO 14                               
      IF ((NA.GT.2).OR.(NB.GT.2)) II=5                                  
   14 RAB=U(IG,IK)/.529167                                              
      IF(NA.GT.NB) GO TO 1                                              
      ISP=2                                                             
      IPS=3                                                             
      FACTOR = +1.0                                                     
      SA=ZA(I1)                                                         
      SB=ZA(J1)                                                         
      GO TO 2                                                           
    1 ISP=3                                                             
      IPS=2                                                             
      FACTOR = -1.0                                                     
      SA=ZA(J1)                                                         
      SB=ZA(I1)                                                         
    2 J=II+2                                                            
      IF(II.GT.3) J=J-1                                                 
      ALPHA=0.5*RAB*(SA+SB)                                             
      BETA=0.5*RAB*(SB-SA)                                              
c-----PREPARE PARAMETERS FOR CONSTRUCTION OF F(0) MATRIX                
      XPI =EXP(DPI *ALPHA)                                              
      XSIG=EXP(DSIG*ALPHA)                                              
      PPPI =(1.+FPPPI * XPI  )/2.                                       
      PPSIG=(1.+FPPSIG* XSIG )/2.                                       
      PS   =(1.+FPS   * XSIG )/2.                                       
      SS   =(1.+FSS   * XSIG )/2.                                       
      DO 3  I=1,J                                                       
      N=I-1                                                             
      A(I)=OA(ALPHA,N)                                                  
      B(I)=OB(BETA ,N)                                                  
    3 CONTINUE                                                          
      DO 4  I=1,5                                                       
    4 Z(I)=0.0                                                          
      GO TO (5,6,7,8,9,10),II                                           
c     ------------------------------------------------------------------
c *** THE ORDERING OF THE ELEMENTS WITHIN Z IS                          
c *** Z(1)=(S(B)/S(A))   Z(2)=(P-SIGMA(B)/S(A))   Z(3)=(S(B)/P-SIGMA(A))
c *** Z(4)=(P-SIGMA(B)/P-SIGMA(A))                Z(5)=(P-PI(B)/P-PI(A))
c     ------------------------------------------------------------------
c *** FIRST ROW - FIRST ROW OVERLAPS                                    
    5 Z(1)=0.25*SQRT((SA*SB*RAB**2)**3)*(A(3)*B(1)-B(3)*A(1))           
      GO TO 11                                                          
c *** FIRST ROW - SECOND ROW OVERLAPS                                   
    6 W=SQRT((SA**3)*(SB**5))*(RAB**4)*0.125                            
      Z(1)=W*SQRT(1.0/3.0)*(A(4)*B(1)-B(4)*A(1)+A(3)*B(2)-B(3)*A(2))    
      Z(ISP)=W*(A(3)*B(1)-B(3)*A(1)+A(4)*B(2)-B(4)*A(2))*FACTOR         
      GO TO 11                                                          
c *** FIRST ROW - THIRD ROW OVERLAPS                                    
    7 W=SQRT((SA**3)*(SB**7)/7.5)*(RAB**5)*0.0625                       
      Z(1)=W*(A(5)*B(1)-B(5)*A(1)+2.0*(A(4)*B(2)-B(4)*A(2)))/SQRT(3.0)  
      Z(ISP)=W*FACTOR*(A(4)*(B(1)+B(3))-B(4)*(A(1)+A(3))                
     1                +B(2)*(A(3)+A(5))-A(2)*(B(3)+B(5)))               
      GO TO 11                                                          
c *** SECOND ROW - SECOND ROW OVERLAPS                                  
    8 W=SQRT((SA*SB)**5)*(RAB**5)*0.0625                                
      RT3=1.0/SQRT(3.0)                                                 
      Z(1)=W*(A(5)*B(1)+B(5)*A(1)-2.0*A(3)*B(3))/3.0                    
      D=A(4)*(B(1)-B(3))-A(2)*(B(3)-B(5))                               
      E=B(4)*(A(1)-A(3))-B(2)*(A(3)-A(5))                               
      Z(ISP)=W*RT3*FACTOR*(D+E)                                         
      Z(IPS)=W*RT3*FACTOR*(E-D)                                         
      Z(4)=W*(B(3)*(A(5)+A(1))-A(3)*(B(5)+B(1)))                        
      Z(5)=0.5*W*(A(5)*(B(1)-B(3))-B(5)*(A(1)-A(3))-A(3)*B(1)+B(3)*A(1))
      GO TO 11                                                          
c *** SECOND ROW - THIRD ROW OVERLAPS                                   
    9 W=SQRT((SA**5)*(SB**7)/7.5)*(RAB**6)*0.03125                      
      RT3=1.0/SQRT(3.0)                                                 
      Z(1)=W*(A(6)*B(1)+A(5)*B(2)-2.0*(A(4)*B(3)+A(3)*B(4))+A(2)*B(5)   
     1       +A(1)*B(6))/3.0                                            
      Z(ISP)=W*RT3*FACTOR*(A(6)*B(2)+A(5)*B(1)-2.0*(A(4)*B(4)+A(3)*B(3))
     1                    +A(2)*B(6)+A(1)*B(5))                         
      Z(IPS)=W*RT3*FACTOR*(A(5)*(2.0*B(3)-B(1))-B(5)*(2.0*A(3)-A(1))    
     1                    -A(2)*(B(6)-2.0*B(4))+B(2)*(A(6)-2.0*A(4)))   
      Z(4)=W*(B(4)*(A(1)+A(5))-A(4)*(B(1)+B(5))+B(3)*(A(2)+A(6))        
     1      -A(3)*(B(2)+B(6)))                                          
      Z(5)=0.5*W*(A(6)*(B(1)-B(3))-B(6)*(A(1)-A(3))+A(5)*(B(2)-B(4))    
     1       -B(5)*(A(2)-A(4))-A(4)*B(1)+B(4)*A(1)-A(3)*B(2)+B(3)*A(2)) 
      GO TO 11                                                          
c *** THIRD ROW - THIRD ROW OVERLAPS                                    
   10 W=SQRT((SA*SB*RAB**2)**7)/480.0                                   
      RT3=1.0/SQRT(3.0)                                                 
      Z(1)=W*(A(7)*B(1)-3.0*(A(5)*B(3)-A(3)*B(5))-A(1)*B(7))/3.0        
      D=A(6)*(B(1)-B(3))-2.0*A(4)*(B(3)-B(5))+A(2)*(B(5)-B(7))          
      E=B(6)*(A(1)-A(3))-2.0*B(4)*(A(3)-A(5))+B(2)*(A(5)-A(7))          
      Z(ISP)=W*RT3*FACTOR*(D-E)                                         
      Z(IPS)=W*RT3*FACTOR*(-D-E)                                        
      Z(4)=W*(A(3)*(B(7)+2.0*B(3))-A(5)*(B(1)+2.0*B(5))-B(5)*A(1)       
     1      +A(7)*B(3))                                                 
      Z(5)=0.5*W*(A(7)*(B(1)-B(3))+B(7)*(A(1)-A(3))                     
     1   +A(5)*(B(5)-B(3)-B(1))+B(5)*(A(5)-A(3)-A(1))+2.0*A(3)*B(3))    
   11 U(IK,IG)= Z(1) * PENSS                                            
      S(IK,IG)= Z(1) * SS                                               
      S(IG,IK)= Z(1)                                                    
      IF (II.EQ.1) GO TO 15                                             
c      COMPUTE THE BOND-DIRECTION COSINES                                
      IG=MAX0(I1,J1)                                                    
      IK=MIN0(I1,J1)                                                    
      T(1)=(X(I1)-X(J1))/U(IG,IK)                                       
      T(2)=(Y(I1)-Y(J1))/U(IG,IK)                                       
      T(3)=(Z1(I1)-Z1(J1))/U(IG,IK)                                     
      IF (NA.LT.3) GO TO 17                                             
c     S-PX, S-PY, S-PZ                                                  
c-----ERROR CORRECTED: J2-->I2, J1-->I1           /JS-L                 
      I2=I1                                                             
      DO 18 K=1,3                                                       
      I2=I2+NX                                                          
      IG=MAX0(J1,I2)                                                    
      IK=MIN0(J1,I2)                                                    
      SH=T(K)*Z(3)                                                      
      U(IK,IG)= SH * PENSP                                              
      S(IK,IG)= SH * PS                                                 
   18 S(IG,IK)= SH                                                      
   17 IF (NB.LT.3) GO TO 19                                             
c     PX-S, PY-S, PZ-S                                                  
      J2=J1                                                             
      DO 20 K=1,3                                                       
      J2=J2+NX                                                          
      IG=MAX0(I1,J2)                                                    
      IK=MIN0(I1,J2)                                                    
      SH=T(K)*Z(2)                                                      
      U(IK,IG) = SH * PENSP                                             
      S(IK,IG) = SH * PS                                                
   20 S(IG,IK) = SH                                                     
   19 IF (II.LT.4) GO TO 15                                             
c       XX-, XY-, XZ-, YY-, YZ-, ZZ-WECHSELWIRKUNG                       
      I2=I1                                                             
      DO 21 K=1,3                                                       
      I2=I2+NX                                                          
      J2=J1                                                             
      DO 21 L=1,3                                                       
      J2=J2+NX                                                          
      IG=MAX0(I2,J2)                                                    
      IK=MIN0(I2,J2)                                                    
      FF=0.                                                             
      IF (K.EQ.L) FF=1.                                                 
      SH1=T(K)*T(L)*Z(4)                                                
      SH2=(FF-T(K)*T(L))*Z(5)                                           
      U(IK,IG)= (SH1 + SH2) * PENPP                                     
      S(IK,IG)= SH1 * PPSIG + SH2 * PPPI                                
      S(IG,IK)= SH1         + SH2                                       
c-MCD- -MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD- 
      IF(K*L.NE.9)GO TO 21                                              
      BETAM(I1-NH,J1-NH)= SH2 * ( 1.+ FPPPIM * EXP(FDM*ALPHA) ) / 2.    
      BETAM(J1-NH,I1-NH)=BETAM(I1-NH,J1-NH)                             
c-MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD- 
   21 CONTINUE                                                          
   15 CONTINUE                                                          
   16 CONTINUE                                                          
      IF (IPRNT.LT.3) RETURN                                            
      WRITE(6,101)                                                      
      CALL PRINT(S,N5,0,0,0,NX,NH)                                      
      RETURN                                                            
      END                                                               
c#####################################################################  
      FUNCTION OA(A,K)                                                  
      implicit real*8(A-H,O-Z)
c     THIS QCPE 228 FUNCTION EVALUATES THE A INTEGRALS                  
      B=1./A                                                            
      S=1.                                                              
      OA=1.                                                             
      IF (K.LT.1) GO TO 2                                               
      DO 1 M=1,K                                                        
      L=K-M+1                                                           
      S=L*S*B                                                           
    1 OA=OA+S                                                           
    2 OA=OA*B*EXP(-A)                                                   
      RETURN                                                            
      END                                                               
c#####################################################################  
      FUNCTION OB(BETA,N)                                               
      implicit real*8(A-H,O-Z)
c      THIS QCPE 228 FUNCTION EVALUATES THE B INTEGRALS                  
      B=BETA**2                                                         
      FN=N                                                              
      J=MOD(N,2)                                                        
      IF (J.EQ.0) GO TO 2                                               
      FNUMER=FN+2.                                                      
      SUM=BETA/FNUMER                                                   
      FACTOR=-2.                                                        
      FI=3.0                                                            
      GO TO 3                                                           
    2 FNUMER=FN+1.                                                      
      SUM=1./FNUMER                                                     
      FACTOR=2.                                                         
      FI=2.                                                             
    3 T=SUM                                                             
    4 DENOM=FNUMER+2.                                                   
      T=(T/FI)*(B/(FI-1.))*(FNUMER/DENOM)                               
      IF (T.LE.1.E-10) GO TO 5                                          
      SUM=T+SUM                                                         
      FI=FI+2.                                                          
      FNUMER=DENOM                                                      
      GO TO 4                                                           
    5 OB=FACTOR*SUM                                                     
      RETURN                                                            
      END                                                               
c#####################################################################  
      SUBROUTINE F0MATR                                                 
      implicit real*8(A-H,O-Z)
c-----COMPUTE F0 IN AO BASIS; DIAGONAL IN DDD(*)                        
c     ON ENTRY:  F0( S ) SF ) , SP( R ) SP )
c     J. Spanget-Larsen, Theor. Chem. Acc. 98, 137-153 (1997)                        
      COMMON /A/ N,NH,NE,NQ,ECORE,ICOUNT,JCOUNT,EMAX,IPRNT              
      COMMON /B/ E(300),EE(300),ZCORE(300),ET(300),SUM(300)                  
      COMMON /C/ F0(300,300),SP(300,300)                                    
      COMMON /E/ NCYCL,EPREV,ITERR,N5,NTOTAL                            
      COMMON /F/ PENSS,PENSP,PENPP,XG,FPPPI,FPPSIG,FPS,FSS,             
     .           DPI,DSIG,DDD(300),FPPPIM,FDM                            
      COMMON /PIMCD/ MCD,NC,NOCC,NST,LPI(300),XP(100),YP(100),BET(5050),   
     1  BETAM(100,100),CIPI(10000),IFROM(100),ITO(100),EP(100),POCC(100)       
      DO 10 I=1,N5                                                      
      DO 10 J=1,I                                                       
   10 SP(I,J)=0.0                                                       
c-----COMPUTE PENETRATION TERMS FOR ONE-CENTER F0 ELEMENTS              
      DO 40 I=1,NTOTAL                                                  
      J=1                                                               
      IF(I.GT.NH) J=4                                                   
      DO 40 I1=1,J                                                      
      IG=I+(I1-1)*N                                                     
      DO 40 I2=1,I1                                                     
      IK=I+(I2-1)*N                                                     
      TEMP=0.0                                                          
      DO 30 L=1,N5                                                      
      IG1=MAX0(IK,L)                                                    
      IK1=MIN0(IK,L)                                                    
      IG2=MAX0(L,IG)                                                    
      IK2=MIN0(L,IG)                                                    
   30 TEMP=TEMP + SP(IK1,IG1) * DDD(L) * SP(IK2,IG2)                    
      IF(IK.LT.IG) SP(IG,IK)=TEMP                                       
   40 IF(IK.EQ.IG) SUM(IK)  =TEMP                                       
c-----OFF-DIAGONAL F0 ELEMENTS (IN UPPER HALF OF F0)                    
      DO 50 IG=1,N5                                                     
      DO 50 IK=1,IG                                                     
   50 IF(IK.LT.IG) F0(IK,IG) = SP(IG,IK) + F0(IK,IG)*(ET(IK)+ET(IG))    
c-----DIAGONAL F0 ELEMENTS (IN DDD)                                     
      DO 60 I=1,N5                                                      
   60 DDD(I)=DDD(I)+SUM(I)                                              
c-MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD- 
      DO 6000 I=1,N                                                     
      BETAM(I,I)=DDD(N+NH+I)-SUM(N+NH+I)                                
      DO 6000 J=1,I                                                     
      IF(I.NE.J) BETAM(I,J)=BETAM(I,J)*(BETAM(I,I)+BETAM(J,J))          
 6000 BETAM(J,I)=BETAM(I,J)                                             
c-MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD- 
      WRITE(6,1000) (SUM(I),I=1,N5)                                     
 1000 FORMAT(//' DIAGONAL PENETRATION TERMS'//(10F13.8))               
      RETURN                                                            
c     ON EXIT:  F0( S ) F0* ), SP( - ! SP )                             
      END                                                               
c#####################################################################
      SUBROUTINE LOWDIN                                                 
      implicit real*8(A-H,O-Z)
c-----LOWDIN TRANSFORM F0 MATRIX                        /JS-L           
      COMMON /A/ N,NH,NE,NQ,ECORE,ICOUNT,JCOUNT,EMAX,IPRNT,IGCCI,       
     1           EATOM,ISCF,ITGT,ICI,GAH,ITRIP,IPEN,ETOT,ETOT1,IDOB     
      COMMON /B/ E(300),EE(300),ZCORE(300),ET(300),SUM(300)                  
      COMMON /C/ H(300,300),W(300,300)                                      
      COMMON /E/ NCYCL,EPREV,ITERR,N5,NTOTAL,N2S,N2P,NTN                
      COMMON /F/ PENSS,PENSP,PENPP,XG,FPPPI,FPPSIG,FPS,FSS,             
     .     DPI,DSIG,DDD(300),FPPPIM,FDM                                  
      COMMON /PIMCD/ MCD,NC,NOCC,NST,LPI(300),XP(100),YP(100),
     1               BETA(5050),BETAM(100,100),CIPI(10000),
     2               IFROM(100),ITO(100),EP(100),POCC(100)      
      COMMON /SBIAS/ BIAS,IPI                                           
c-----COMPUTE INVERSE SQUARE ROOT OF OVERLAP-MATRIX                     
   20 NTEMP=N                                                           
      N=N5                                                              
      HBIAS = BIAS                                                      
      BIAS = 0.0                                                        
      CALL DIAGON                                                       
      BIAS = HBIAS                                                      
      N=NTEMP                                                           
      IF(IPRNT.LT.4) GOTO 30                                            
      WRITE(6,1000) (SUM(I),I=1,N5)                                     
 1000 FORMAT(//' EIGENVALUES OF OVERLAP-MATRIX'//(10F13.8))            
      WRITE(6,1100)                                                     
 1100 FORMAT(//' EIGENVECTORS OF OVERLAP-MATRIX')                      
      CALL PRINT(W,N5,N5,0,0,N,NH)                                      
   30 DO 40 I=1,N5                                                      
      IF(SUM(I).GE..0001) GOTO 40                                       
      WRITE(6,1200) I,SUM(I)                                            
 1200 FORMAT(//I4,'. EIGENVALUE OF OVERLAP-MATRIX =',F13.8,            
     .  //' MATRIX NOT POSITIVE DEFINITE, CALCULATION INTERRUPTED')     
      STOP                                                              
   40 SUM(I)=1./SQRT(SUM(I))                                            
      DO 60 I=1,N5                                                      
      DO 60 J=1,I                                                       
      TEMP=0.0                                                          
      DO 50 K=1,N5                                                      
   50 TEMP=TEMP+W(I,K)*SUM(K)*W(J,K)                                    
   60 H(I,J)=TEMP                                                       
      IF(IPRNT.LT.3) GOTO 70                                            
      WRITE(6,1300)                                                     
 1300 FORMAT(//' INVERSE SQUARE ROOT OF OVERLAP-MATRIX')               
      CALL PRINT(H,N5,0,0,0,N,NH)                                       
c-  - -REARRANGE                                                         
   70 DO 90 I=1,N5                                                      
      DO 80 J=1,I                                                       
      W(I,J)=H(I,J)                                                     
   80 W(J,I)=W(I,J)                                                     
   90 H(I,I)=DDD(I)                                                     
c- - -PRINT                                                             
      IF(IPRNT.GT.2) WRITE(6,1350)                                      
 1350 FORMAT(//'F(0) MATRIX IN AO BASIS (EV)')       
      IF(IPRNT.GT.2) CALL PRINT(H,0,N5,0,0,N,NH)                        
c-----LOWDIN TRANSFORM F(0) MATRIX                                      
      DO 150 I=1,N5                                                     
      DO 130 L=1,N5                                                     
      TEMP=0.0                                                          
      IF(L.EQ.1) GOTO 110                                               
      L1=L-1                                                            
      DO 100 K=1,L1                                                     
  100 TEMP=TEMP+H(K,L)*W(I,K)                                           
  110 DO 120 K=L,N5                                                     
  120 TEMP=TEMP+H(L,K)*W(I,K)                                           
  130 SUM(L)=TEMP                                                       
      DO 150 J=1,I                                                      
      TEMP=0.0                                                          
      DO 140 L=1,N5                                                     
  140 TEMP=TEMP+SUM(L)*W(L,J)                                           
      IF(I.EQ.J) EE(I) =TEMP                                            
      IF(I.NE.J) H(I,J)=TEMP                                            
  150 CONTINUE                                                          
      DO 160 I=1,N5                                                     
  160 H(I,I)=EE(I)                                                      
c- - -PRINT                                                             
  180 IF(IPRNT.GT.1) WRITE(6,1400)                                      
 1400 FORMAT(//'F(0) MATRIX IN LOWDIN BASIS (EV)')   
      IF(IPRNT.GT.1) CALL PRINT(H,N5,0,0,0,N,NH)                        
c- - -STORE                                                             
  190 CALL STORE(H,1,N5)                                                
c-MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD- 
      IR=0                                                              
      NP=NH+3*N                                                         
      DO 3002 I=1,N                                                     
      DO 3001 L=1,N                                                     
      SUM(L)=.0                                                         
      DO 3001 K=1,N                                                     
 3001 SUM(L)=SUM(L)+BETAM(K,L)*W(NP+I,NP+K)                             
      DO 3002 J=1,I                                                     
      IR=IR+1                                                           
      BETA(IR)=.0                                                       
      DO 3002 L=1,N                                                     
 3002 BETA(IR)=BETA(IR)+SUM(L)*W(NP+L,NP+J)                             
c-MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD- 
      RETURN                                                            
      END                                                               
c#####################################################################  
      SUBROUTINE SCF                                                    
      implicit real*8(A-H,O-Z)
c-----SCF PROCEDURE  /JS-L                                              
      COMMON /A/ N,NH,NE,NQ,ECORE,ICOUNT,JCOUNT,EMAX,IPRNT,IGCCI,       
     1       EATOM,ISCF,ITGT,ICI,GAH,ITRIP,IPEN,ETOT,ETOT1,IDOB,NB,NAB  
      COMMON /B/ E(300),EE(300),ZCORE(300),ET(300),SUM(300),
     1           C(400),GAX(300),NTIT(20)                                               
      COMMON /C/ S(300,300),U(300,300)                                      
      COMMON /E/ NCYCL,EPREV,ITERR,N5,NTOTAL,N2S,N2P,NTN,NG,DAMP        
      COMMON /MOM/ SX(300),SY(300)                                        
      COMMON /OCCDAT/ NHOMO,NLUMO,NPE,NPMO,OCC(300)                      
      COMMON /F/ PENSS,PENSP,PENPP,XG,FPPPI,FPPSIG,FPS,FSS,             
     .           DPI,DSIG,DDD(300),FPPPIM,FDM                            
      COMMON /SBIAS/ BIAS,IPI
      COMMON /NMOP/ NMOP                                           
   72 FORMAT (I4,5X,'D=',F4.2,5X,'DE =',F15.8,' EV',5X,                 
     .        'DP=',F12.8,5X)                               
   73 FORMAT (I4,5X,'D=',F4.2,5X,'DE =',F15.8,' EV',5X,                 
     .        'DP=',F12.8,5X,'  ENERGIES SATISFIED')                               
   74 FORMAT (I4,5X,'D=',F4.2,5X,'DE =',F15.8,' EV',5X,                 
     .        'DP=',F12.8,5X,'  DENSITIES SATISFIED')                               
   75 FORMAT (I4,5X,'D=',F4.2,5X,'DE =',F15.8,' EV',5X,                 
     .        'DP=',F12.8,5X,'  ENERGIES AND DENSITIES SATISFIED')                               
  820 FORMAT (10X,I3,': ',F12.6,F12.6)                                  
 1961 FORMAT (//' MO ENERGIES (EV) AND HALF OCCUPATION NUMBERS'/)       
 3002 FORMAT (//' MO COEFFICIENTS STORED COLUMNWISE (LOWDIN BASIS)')      
 3003 FORMAT (//' DENSITY MATRIX (LOWDIN BASIS)')                         
      NTEMP=N                                                           
      LOOK=1                                                            
      IF(NPMO+NPE.NE.0) LOOK=0                                          
      ETOT1=0.                                                          
      CALL LOAD(S,1,N5)                                                 
      DO 1 I=1,N5                                                       
      ET(I)=0.0                                                         
    1 SX(I)=0.0                                                         
      DAMP=ALOG(DAMP)                                                   
      DD=1.                                                             
c-----SCF LOOP                                                          
      WRITE(6,3009) (NTIT(I),I=1,18)                                     
 3009 FORMAT(1H //1X,18(A4)//' INITIAL DENSITIES AND CHARGES:'/)      
      DO 10 NCYCL=1,ISCF                                                
      N=N5                                                              
      CALL DIAGON                                                       
      N=NTEMP                                                           
      CALL AUFBAU(NE,N5,SUM,LOOK)                                       
      IF(ITERR.NE.0) GOTO 30                                            
      CALL FMATR(DE,DP,NCYCL,DD)                                        
      IF(NCYCL.EQ.1) WRITE(6,3050)                                      
 3050 FORMAT(/' SCF PROCEDURE:'/)                                       
	  IF(DP.LT.1.E-4 .AND. DE.LT.5.E-4) WRITE(6,75) NCYCL,DD,DE,DP
	  IF(DP.LT.1.E-4 .AND. DE.LT.5.E-4) GO TO 20
      IF(DP.LT.1.E-4) WRITE(6,74) NCYCL,DD,DE,DP                                     
	  IF(DP.LT.1.E-4) GO TO 10
      IF(DE.LT.5.E-4) WRITE(6,73) NCYCL,DD,DE,DP 
      IF(DE.LT.5.E-4) GO TO 10                                    
      WRITE(6,72) NCYCL,DD,DE,DP                                   
   10 CONTINUE                                                          
      WRITE(6,3006)                                                      
 3006 FORMAT(//' SCF PROCEDURE INTERRUPTED, CONVERGENCE NOT ATTAINED'//)
c-----TOTAL ENERGY                                                      
   20 N=N5                                                              
      CALL DIAGON                                                       
      N=NTEMP                                                           
      BIAS=0.0                                                          
      IF(IPRNT.LT.2) GOTO 30                                            
      WRITE(6,3001)                                                     
 3001 FORMAT(//' GROUND STATE FOCK MATRIX IN LOWDIN BASIS (EV)')        
      CALL PRINT(S,N5,0,0,0,N,NH)                                       
   30 CALL LOAD(S,2,N5)                                                 
      CALL SHIFT(S,N5)                                                  
      CALL ENERGY(1,ECORE,SUM)                                          
      WRITE(6,1961)                                                     
c-----PRINT MO ENERGIES AND OCCUPATION NUMBERS                          
      DO 40 I=1,N5                                                      
      J=N5+1-I                                                          
      WRITE(6,820) J,SUM(J),OCC(J)                                      
      IF(J.EQ.1) GOTO 40                                                
      IF(ABS(OCC(J)-OCC(J-1)).GT.0.001) WRITE(6,3010)                   
 3010 FORMAT(1X)                                                        
   40 CONTINUE                                                          
      WRITE(6,3002)                                                     
      MIN=MAX0(1,NLUMO-NMOP)                                              
      MAX=MIN0(N5,NHOMO+NMOP)                                             
      CALL PRINT(U,N5,N5,MIN,MAX,N,NH)                                  
      CALL STORE(U,4,N5)                                                
      IF (IPRNT.LT.2) GO TO 83                                          
      WRITE(6,3003)                                                     
      CALL PMATR (S,U,N5)                                               
      CALL PRINT(S,N5,0,0,0,N,NH)                                       
   83 WRITE(6,3005)                                                      
 3005 FORMAT(//' ATOMIC POPULATIONS'/)                                  
      CALL PRATP(SX,N,NH,ZCORE)                                         
      CALL PI                                                           
      RETURN                                                            
      END                                                               
c###################################################################    
      SUBROUTINE PI                                                     
      implicit real*8(A-H,O-Z)
      COMMON /A/ NC,NH                                                  
      COMMON /C/ S(300,300),U(300,300)                                      
      COMMON /PIMCD/ MCD,N,NOCC,NST,LPI(300),XP(100),YP(100),BETA(5050),   
     1  CPI(10000),CIPI(10000),IFROM(100),ITO(100),EP(100),POCC(100)          
      COMMON /MOM/ Q(300)                                                
      COMMON /OCCDAT/ NHOMO,NLUMO,NPE,NPMO,OCC(300)                      
      COMMON /SBIAS/ BIAS,IPI                                           
      REAL NEL                                                          
      N=NC                                                              
      NPI=NH+3*NC+1                                                     
      NORB=NH+4*NC                                                      
      IP=0                                                              
      IR=0                                                              
      NEL=0.0                                                           
      DO 20 MO=1,NORB                                                   
      LPI(MO)=1                                                         
      IF(IPI.LE.0) GO TO 20                                             
      SPI=0.0                                                           
      IX=0                                                              
      DO 10 MY=NPI,NORB                                                 
      IX=IX+1                                                           
   10  SPI=SPI+ABS(U(MY,MO))                                             
      IF(SPI.LT.0.1) LPI(MO)=0                                          
      IF(SPI.LT.0.1) GO TO 20                                           
      IP=IP+1                                                           
      POCC(IP)=OCC(MO)                                                  
      IF(OCC(MO).GT.0.01)  NOCC=IP                                      
      LPI(MO)=IP                                                        
      NEL=NEL + 2.0 * OCC(MO)                                           
      IX=0                                                              
      DO 15 MY=NPI,NORB                                                 
      IX=IX+1                                                           
      C=0.0                                                             
      C=U(MY,MO)                                                        
      IR=IR+1                                                           
   15  CPI(IR)=C                                                         
   20 CONTINUE                                                          
      IF(IPI.LE.0) RETURN                                               
c-----DETERMINE CENTER OF CHARGE                                        
      IF( NEL .GT. 0.1 ) GO TO 25                                       
      WRITE(6,500) NEL                                                  
  500 FORMAT(//1X,'   ERROR   NO. OF PI ELECTRONS',F6.2//)              
      NEL=1.0                                                           
   25 X0=.0                                                             
      Y0=.0                                                             
      IR=0                                                              
      DO 30 MY=NPI,NORB                                                 
      IR=IR+1                                                           
      X0=X0+XP(IR)*Q(MY)/NEL                                            
      Y0=Y0+YP(IR)*Q(MY)/NEL                                            
  30  CONTINUE                                                          
      WRITE(6,1000) X0,Y0                                               
 1000 FORMAT(' PI ELECTRON CENTER OF CHARGE:   X,Y  =',F8.4,' ,',F8.4) 
c-----SHIFT TO CHARGE CENTERED COORDINATE SYSTEM                        
      WRITE(6,2000)                                                     
 2000 FORMAT(/1X,'CHARGE CENTERED X,Y-COORDINATES FOR PI SYSTEM:'/)     
      DO 40 MY=1,N                                                      
      XP(MY)=XP(MY)-X0                                                  
      YP(MY)=YP(MY)-Y0                                                  
   40 WRITE(6,3000) MY,XP(MY),YP(MY)                                     
 3000 FORMAT(I5,2F9.4)                                                  
      RETURN                                                            
      END                                                               
c###################################################################    
      SUBROUTINE PRATP(PAO,NC,NH,ZCORE)                                 
      implicit real*8(A-H,O-Z)
      DIMENSION PAO(1),ZCORE(1),PA(300),CHARGE(300)                       
c-----COMPUTE ATOMIC POPULATIONS AND NET CHARGES                        
      SPA=0.0                                                           
      SCH=0.0                                                           
      NT=NC+NH                                                          
      DO 20 I=1,NT                                                      
      PA(I)=PAO(I)                                                      
      IF(I.LE.NH) GOTO 10                                               
      PA(I)=PAO(I)+PAO(I+NC)+PAO(I+2*NC)+PAO(I+3*NC)                    
   10 CHARGE(I)=ZCORE(I)-PA(I)                                          
      SPA=SPA+PA(I)                                                     
   20 SCH=SCH+CHARGE(I)                                                 
c-----PRINT                                                             
      WRITE(6,1000)                                                     
 1000 FORMAT(' POSITION',11X,'AO POPULATIONS',16X,'ATOMIC      NET',    
     .      /13X,'S',8X,'PX',7X,'PY',7X,'PZ',6X,'POPULATION   CHARGE',  
     .      /1X,67('-'))                                                
      IF(NH.EQ.0) GOTO 40                                               
      DO 30 I=1,NH                                                      
   30 WRITE(6,2000) I,PAO(I),PA(I),CHARGE(I)                            
 2000 FORMAT(I5,4X,F9.4,27X,2F11.4)                                     
   40 IF(NC.EQ.0) GOTO 60                                               
      DO 50 J=1,NC                                                      
      I=NH+J                                                            
   50 WRITE(6,3000) I,PAO(I),PAO(I+NC),PAO(I+2*NC),PAO(I+3*NC),PA(I),   
     .              CHARGE(I)                                           
 3000 FORMAT(I5,4X,4F9.4,2F11.4)                                        
   60 WRITE(6,4000) SPA,SCH                                             
 4000 FORMAT(/42X,'SUM',2F11.4/1X)                                      
      RETURN                                                            
      END                                                               
c###################################################################    
      SUBROUTINE ENERGY(IFLAG,ECORE,EMO)                                
      implicit real*8(A-H,O-Z)
      DIMENSION EMO(1)                                                  
      COMMON /OCCDAT/ NHOMO,NLUMO,NPE,NPMO,OCC(300)                      
c-----HUECKEL ENERGY                                                    
      EHUECK=0.0                                                        
      DO 10 I=1,NHOMO                                                   
   10 EHUECK = EHUECK + 2.0 * OCC(I) * EMO(I)                           
      WRITE(6,1000) EHUECK                                              
 1000 FORMAT(//' HUECKEL ENERGY                  ',F15.6,' EV')        
      IF(IFLAG.LT.1) RETURN                                             
c-----ROOTHAAN CLOSED SHELL ELECTRONIC ENERGY                           
      ER=EHUECK                                                         
c- - -DIAGONAL CONTRIBUTIONS                                            
      DO 20 I=1,NHOMO                                                   
   20 ER = ER - OCC(I)**2 * H(I,I,I,I)                                  
c- - -OFF-DIAGONAL CONTRIBUTIONS                                        
      NHOMO1=NHOMO-1                                                    
      DO 30 I=1,NHOMO1                                                  
      I1=I+1                                                            
      DO 30 J=I1,NHOMO                                                  
   30 ER = ER - 4.0 * OCC(I) * OCC(J) * ( H(I,I,J,J) - H(I,J,J,I)/2.0 ) 
      E = ER + ECORE                                                    
      WRITE(6,2000) ER,ECORE,E                                          
 2000 FORMAT(  /' ROOTHAAN CLOSED SHELL ENERGY    ',F15.6,' EV',        
     .        //' CORE REPULSION                  ',F15.6,' EV',        
     .        //' TOTAL VALENCE SHELL ENERGY      ',F15.6,' EV')        
      IF(NHOMO.LT.NLUMO) RETURN                                         
c-----PROJECTION OF OPEN SHELL STATES                                   
      WRITE(6,3000)                                                      
 3000 FORMAT(//' PROJECTION OF OPEN SHELL STATES:')                     
      IF(NPMO-2) 40,50,80                                               
c- - -A SINGLE UNPAIRED ELECTRON IN ONE SPACE-ORBITAL                   
   40 E = E - H(NHOMO,NHOMO,NHOMO,NHOMO)/4.                             
      WRITE(6,4000) E                                                    
 4000 FORMAT(  /' DOUBLET                         ',F15.6,' EV')        
      RETURN                                                            
   50 A = H(NHOMO,NHOMO,NHOMO,NHOMO)                                    
      B = H(NLUMO,NLUMO,NLUMO,NLUMO)                                    
      C = H(NHOMO,NHOMO,NLUMO,NLUMO)                                    
      D = H(NHOMO,NLUMO,NLUMO,NHOMO)                                    
      IF(NPE-2) 70,60,70                                                
c- - -TWO ELECTRONS IN TWO DEGENERATE SPACE-ORBITALS                    
   60 E1 = E - A/4. - B/4.     - 0.5*D                                  
      E2 = E - A/4. - B/4.     + 1.5*D                                  
      E3 = E + A/4. + B/4. - C - 0.5*D                                  
      E4 = E + A/4. + B/4. - C + 1.5*D                                  
      X1 = 0.                                                           
      X2 = E2 - E1                                                      
      X3 = E3 - E1                                                      
      X4 = E4 - E1                                                      
      WRITE(6,5000) E1,X1,E2,X2,E3,X3,E4,X4                             
 5000 FORMAT(  /' TRIPLET        U V              ',2F15.6,' EV',       
     .        //'             (  U V              ',2F15.6,' EV',       
     .         /'             (',                                       
     .         /' SINGLETS    (  U U - V V        ',2F15.6,' EV',       
     .         /'             (',                                       
     .         /'             (  U U + V V        ',2F15.6,' EV')       
      RETURN                                                            
c- - -ONE OR THREE ELECTRONS IN TWO DEGENERATE SPACE-ORBITALS           
   70 E = E - A/16. - B/16. - C/4. + D/8.                               
      WRITE(6,4000) E                                                   
      RETURN                                                            
   80 WRITE(6,6000)                                                     
 6000 FORMAT(//' (PRESENT VERSION DOES NOT TREAT SPACE DEGENERACIES',  
     .          ' LARGER THAN TWO)')                                    
      RETURN                                                            
      END                                                               
c###################################################################    
      SUBROUTINE PMATR (W,U,N)                                          
      implicit real*8(A-H,O-Z)
c-----RETURNS BOND ORDERS IN LOWER HALF OF W                            
      COMMON /OCCDAT/ NHOMO,NLUMO,NPE,NPMO,OCC(300)                      
      DIMENSION W(300,300),U(300,300)                                       
      DO 20 I=1,N                                                       
      DO 20 J=1,I                                                       
      TEMP=0.0                                                          
      DO 10 K=1,NHOMO                                                   
   10 TEMP=TEMP+OCC(K)*U(I,K)*U(J,K)                                    
   20 W(I,J)=TEMP+TEMP                                                  
      RETURN                                                            
      END                                                               
c###################################################################    
      SUBROUTINE FMATR(DE,DP,NCYC,DD)                                   
      implicit real*8(A-H,O-Z)
c-----CONSTRUCT CLOSED SHELL OR PSEUDO-CLOSED SHELL CNDO FOCK MATRIX    
      COMMON /A/ N,NH,NE,NQ,ECORE,ICOUNT,JCOUNT,EMAX,IPRNT,IGCCI,       
     1       EATOM,ISCF,ITGT,ICI,GAH,ITRIP,IPEN,ETOT,ETOT1              
      COMMON /B/ E(300),EE(300),ZCORE(300),ET(300),SUM(300),
     1           C(400),GAX(300)   
      COMMON /C/ S(300,300),U(300,300)                                      
      COMMON /E/ NCYCL,EPREV,ITERR,N5,NTOTAL,N2S,N2P,NTN,NG,DAMP        
      COMMON /MOM/ SX(300),SY(300),SZ(300)                                 
c-----COMPUTE BOND ORDER MATRIX                                         
      CALL PMATR(S,U,N5)                                                
c- - -DAMPING                                   /JS-L                   
      IF(DD.GT..05) GOTO 2                                              
      DD=0.0                                                            
      GOTO 4                                                            
    2  IF(NCYC.NE.1) CALL LOAD(U,3,N5)                                   
      CALL STORE(S,3,N5)                                                
      IF(NCYC.EQ.1) GOTO 4                                              
      DD = EXP( DAMP * (NCYC-1) )                                       
      DD1= 1.0 - DD                                                     
      DO 3 I=1,N5                                                       
      DO 3 J=1,I                                                        
    3  S(I,J) = DD1 * S(I,J)  +  DD * U(I,J)                             
c-----EVALUATE CONVERGENCE PARAMETERS DE AND DP                         
    4 DE=0.0                                                            
      DP=0.0                                                            
      DO 10 I=1,N5                                                      
      DE=AMAX1(DE,ABS(ET(I)-SUM(I)))                                    
      ET(I)=SUM(I)                                                      
      DP=AMAX1(DP,ABS(SX(I)-S(I,I)))                                    
   10 SX(I)=S(I,I)                                                      
      IF(NCYC.EQ.1) CALL PRATP(SX,N,NH,ZCORE)                           
c-----COMPUTE GROSS ATOMIC ORBITAL AND ATOMIC CHARGES                   
      DO 30 I=1,NTOTAL                                                  
      L=1                                                               
      IF(I.GT.NH) L=4                                                   
      ZC=ZCORE(I)/FLOAT(L)                                              
      TEMP=0.0                                                          
      DO 20 J=1,L                                                       
      K=I+(J-1)*N                                                       
      SY(K)=SX(K)-ZC                                                    
      TEMP=TEMP+SY(K)                                                   
c- - - GROSS AO CHARGE REDEFINED AS *P(K,K)-1* FOR USE OF TRAD. PRM*S:   
   20 SY(K)=SY(K)+ZC-1.                                                 
c  20 CONTINUE                                                          
   30 SZ(I)=TEMP                                                        
c- - -ENTER F(0)                                                        
      CALL LOAD(U,1,N5)                                                 
      DO 40 I=1,N5                                                      
      DO 40 J=1,I                                                       
   40 S(J,I)=U(I,J)                                                     
c- - -ENTER GAMMA MATRIX                                                
      CALL LOAD(U,2,N5)                                                 
c-- ---CONSTRUCT DIAGONAL CNDO FOCK MATRIX ELEMENTS                      
      N3=N*3                                                            
      DO 60 I=1,NTOTAL                                                  
      TEMP=0.0                                                          
      DO 50 J=1,NTOTAL                                                  
      MIN=MIN0(I,J)                                                     
      MAX=MAX0(I,J)                                                     
   50 IF(I.NE.J)  TEMP = TEMP + SZ(J) * U(MIN,MAX)                      
   60 SUM(I)=TEMP                                                       
      DO 80 I=1,NTOTAL                                                  
      U(I,I) = S(I,I) + ( SZ(I) - SY(I)/2. ) * GAX(I)  +  SUM(I)        
      IF(I.LE.NH) GOTO 80                                               
      DO 70 K=N,N3,N                                                    
      L=I+K                                                             
   70 U(L,L) = S(L,L) + ( SZ(I) - SY(L)/2. ) * GAX(I)  +  SUM(I)        
   80 CONTINUE                                                          
c-----CONSTRUCT OFF-DIAGONAL CNDO FOCK MATRIX ELEMENTS                  
      DO 90 I=2,N5                                                      
      I1=I-1                                                            
      DO 90 J=1,I1                                                      
   90 U(I,J) = S(J,I) - 0.5 * S(I,J) * U(J,I)                           
c      PRESENT STORAGE:   U( FOCK ) GAMMA ) ,   S( P ( HCORE )           
c- ----RETURN FOCK MATRIX IN LOWER HALF OF S                             
      DO 110 I=1,N5                                                     
      DO 110 J=1,I                                                      
  110 S(I,J)=U(I,J)                                                     
      RETURN                                                            
      END                                                               
c###################################################################    
      SUBROUTINE DIAGON                                                 
      implicit real*8(A-H,O-Z)
c     EISPACK PROGRAM TRED2 CALCULATES OF THE SYMMETRIC MATRIX A        
c     THE TRIDIAGONAL MATRIX. VECTOR D CONTAINS THE DIAGONAL ELEMENTS   
c     OF THE TRIDIAGONAL MATRIX. VECTOR E CONTAINS IN HIS LAST N-1      
c     POSITIONS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX.     
c     THE ELEMENT E(1) IS SET TO ZERO. Z CONTAINS THE ORTHOGONAL        
c     TRANSFORMATION MATRIX, PRODUCED IN THE REDUCTION TO THE           
c     TRIDIAGONAL FORM.                                                 
        COMMON /A/ N,NH,NE,NM                                             
       COMMON /B/ EL(300),EE(300),ZCORE(300),ET(300),D(300)                   
       COMMON /C/ A(300,300),Z(300,300)                                      
      COMMON /D/ E(300)                                                  
      DO 100 I = 1,N                                                    
      DO 100 J=1,I                                                      
      Z(I,J) = A(I,J)                                                   
  100 CONTINUE                                                          
      IF (N.EQ.1) GO TO 320                                             
c-----FOR I=N STEP -1 UNTIL 2 DO--                                      
      DO 300 II = 2,N                                                   
      I = N+2-II                                                        
      L=I-1                                                             
      H=0.0                                                             
      SCALE = 0.0                                                       
      IF (L.LT.2) GO TO 130                                             
c-----SCALE ROW                                                         
      DO 120 K=1,L                                                      
  120 SCALE = SCALE + ABS(Z(I,K))                                       
      IF (SCALE.NE.0.0) GO TO 140                                       
  130 E(I) = Z(I,L)                                                     
      GO TO 290                                                         
  140 DO 150 K = 1, L                                                   
      Z(I,K) = Z(I,K)/SCALE                                             
      H = H + Z(I,K)*Z(I,K)                                             
  150 CONTINUE                                                          
      F = Z(I,L)                                                        
      G = -SIGN(SQRT(H),F)                                              
      E(I) = SCALE*G                                                    
      H = H-F*G                                                         
      Z(I,L) = F-G                                                      
      F = 0.0                                                           
      DO 240 J = 1,L                                                    
      Z(J,I) = Z(I,J)/(SCALE*H)                                         
      G = 0.0                                                           
c-----FORM ELEMENT OF A*U                                               
      DO 180 K = 1,J                                                    
  180 G = G+Z(J,K)*Z(I,K)                                               
      JP1 = J+1                                                         
      IF (L.LT.JP1) GO TO 220                                           
      DO 200 K = JP1, L                                                 
  200 G=G+Z(K,J)*Z(I,K)                                                 
c-----FORM ELEMENT OF P                                                 
  220 E(J) = G/H                                                        
      F = F + E(J)*Z(I,J)                                               
  240 CONTINUE                                                          
      HH = F/(H+H)                                                      
c-----FORM REDUCED A                                                    
      DO 260 J = 1,L                                                    
      F = Z(I,J)                                                        
      G = E(J) - HH*F                                                   
      E(J) = G                                                          
      DO 260 K = 1,J                                                    
      Z(J,K) = Z(J,K) -F*E(K)-G*Z(I,K)                                  
  260 CONTINUE                                                          
      DO 280 K = 1,L                                                    
  280 Z(I,K) = SCALE*Z(I,K)                                             
  290 D(I) = H                                                          
  300 CONTINUE                                                          
  320 D(1) = 0.0                                                        
      E(1) = 0.0                                                        
c-----ACCUMULATION OF TRANSFORMATION MATRICES                           
      DO 500 I=1,N                                                      
      L = I-1                                                           
      IF (D(I).EQ.0.0) GO TO 380                                        
      DO 360 J = 1,L                                                    
      G = 0.0                                                           
      DO 340 K = 1,L                                                    
  340 G = G + Z(I,K)*Z(K,J)                                             
      DO 360 K = 1,L                                                    
      Z(K,J) = Z(K,J) - G*Z(K,I)                                        
  360 CONTINUE                                                          
  380 D(I) = Z(I,I)                                                     
      Z(I,I) = 1.0                                                      
      IF (L.LT.1) GO TO 500                                             
      DO 400 J = 1,L                                                    
      Z(I,J) = 0.0                                                      
      Z(J,I) = 0.0                                                      
  400 CONTINUE                                                          
  500 CONTINUE                                                          
      CALL IMTQL2                                                       
      RETURN                                                            
      END                                                               
c###################################################################    
      SUBROUTINE IMTQL2                                                 
      implicit real*8(A-H,O-Z)
c     EISPACK PROGRAM IMTQL2 CALCULATESOF A TRIDIAGONAL MATRIX THE EIGEN
c     AND EIGENVALUES. THE DIAGONAL OF THE TRIDIAGONAL MATRIX           
c     IS GIVEN IN D AND THE OFF-DIAGONAL IN E. ON ENTRY Z CONTAINS      
c     THE TRANSFORMATION-MATRIX. AT THE END THE VECTORS AND D           
c     THE EIGENVALUES ARE GIVEN IN D IN ASCENDING ORDER.                
      COMMON /A/ N,NH,NE,NM                                             
      COMMON /B/ EL(300),EE(300),ZCORE(300),ET(300),D(300)                   
      COMMON /C/ A(300,300),Z(300,300)                                      
      COMMON /D/ E(300)                                                  
      COMMON /E/ NCYCL,EPREV,IERR                                       
c-----ACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING                 
c     THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.              
      ACHEP=2.**(-47)                                                   
      IERR=0                                                            
      IF (N.EQ.1) GO TO 1001                                            
      DO 100 I=2,N                                                      
  100 E(I-1)=E(I)                                                       
      E(N)=0.                                                           
      DO 240 L=1,N                                                      
      J=0                                                               
c-----LOOK FOR SMALL SUB-DIAGONAL ELEMENT                               
  105 DO 110 MX=L,N                                                     
      M=MX                                                              
      IF (M.EQ.N) GO TO 120                                             
      IF (ABS(E(M)) .LE. ACHEP*(ABS(D(M))+ABS(D(M+1))))                 
     .         GO TO 120                                                
  110 CONTINUE                                                          
  120 P=D(L)                                                            
      IF (M.EQ.L) GO TO 240                                             
      IF (J.EQ.30) GO TO 1000                                           
      J=J+1                                                             
c-----FORM SHIFT                                                        
      G=(D(L+1)-P)/(2.0*E(L))                                           
      R=SQRT(G**2+1.)                                                   
      G=D(M)-P+E(L)/(G+SIGN(R,G))                                       
      S=1.                                                              
      C=1.                                                              
      P=0.                                                              
      MML=M-L                                                           
c-----FOR I=M-1 STEP -1 UNTIL L DO --                                   
      DO 200 II=1,MML                                                   
      I=M-II                                                            
      F=S*E(I)                                                          
      B=C*E(I)                                                          
      IF (ABS(F).LT.ABS(G)) GO TO 150                                   
      C=G/F                                                             
      R=SQRT(C**2+1.)                                                   
      E(I+1)=F*R                                                        
      S=1./R                                                            
      C=C*S                                                             
      GO TO 160                                                         
  150 S=F/G                                                             
      R=SQRT(S**2+1.)                                                   
      E(I+1)=G*R                                                        
      C=1./R                                                            
      S=S*C                                                             
  160 G=D(I+1)-P                                                         
      R=(D(I)-G)*S+2.*C*B                                               
      P=S*R                                                             
      D(I+1)=G+P                                                        
      G=C*R-B                                                           
c-----FORM VECTOR                                                       
      DO 180 K=1,N                                                      
      F=Z(K,I+1)                                                        
      Z(K,I+1)=S*Z(K,I)+C*F                                             
      Z(K,I)=C*Z(K,I)-S*F                                               
  180 CONTINUE                                                          
  200 CONTINUE                                                          
      D(L)=D(L)-P                                                       
      E(L)=G                                                            
      E(M)=0.                                                           
      GO TO 105                                                         
  240 CONTINUE                                                          
c-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-N
c    SIGMA BIAS FOR SPECIAL PURPOSES:                                  
      CALL SIGMA                                                        
c-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-NB-N
c-----ORDER EIGENVALUES AND EIGENVECTORS                                
      DO 300 II=2,N                                                     
      I=II-1                                                            
      K=I                                                               
      P=D(I)                                                            
      DO 260 J=II,N                                                     
      IF (D(J).GT.P) GO TO 260                                          
      K=J                                                               
      P=D(J)                                                            
  260 CONTINUE                                                          
      IF (K.EQ.I) GO TO 300                                             
      D(K)=D(I)                                                         
      D(I)=P                                                            
      DO 280 J=1,N                                                      
      P=Z(J,I)                                                          
      Z(J,I)=Z(J,K)                                                     
      Z(J,K)=P                                                          
  280 CONTINUE                                                          
  300 CONTINUE                                                          
      GO TO 1001                                                        
c-----SET ERROR FLAG:  NO CONVERGENCE AFTER 30 ITERATIONS               
 1000 IERR=L                                                            
      WRITE(6,101) IERR                                                 
  101 FORMAT (' MORE THAN',I3,' IMTQL2 ITERATIONS')                     
 1001 CONTINUE                                                          
      RETURN                                                            
      END                                                               
c###################################################################    
      SUBROUTINE SIGMA                                                  
      implicit real*8(A-H,O-Z)
c-----SHIFT ALL SIGMA MO ENERGIES BY AN AMOUNT *BIAS*                   
      COMMON /A/ N,NH                                                   
      COMMON /B/ EP(300),EE(300),ZC(300),ET(300),D(300)                      
      COMMON /C/ A(300,300),Z(300,300)                                      
      COMMON /SBIAS/ BIAS,IPI                                           
      IF(ABS(BIAS).LT.0.1) RETURN                                       
      NPI=N-(N-NH)/4+1                                                  
      DO 20 MO=1,N                                                      
      SPI=0.0                                                           
      DO 10 MY=NPI,N                                                    
   10 SPI=SPI+ABS(Z(MY,MO))                                             
      IF(SPI.LT.0.1) D(MO)=D(MO)-BIAS                                   
   20 CONTINUE                                                          
      RETURN                                                            
      END                                                               
c###################################################################    
      SUBROUTINE GCCI                                                   
      implicit real*8(A-H,O-Z)
c     SINGLY EXCITED CONFIGURATIONS                                     
c     PROGRAMMED BY H.BAUMANN                                           
      DIMENSION X(100),Y(100),Z(100),ZS(100),INDIZ(5,300)                    
      EQUIVALENCE (C(1),X(1)),(C(101),Y(1)),(C(201),Z(1)),
     1		(C(301),ZS(1)),(XJSL(1,1),INDIZ(1,1))                                
      COMMON /OCCDAT/ NHOMO,NLUMO,NPE,NPMO,OCC(300)                      
      COMMON /A/ N,NH,NE,NQ,ECORE,ICOUNT,JCOUNT,EMAX,IPRNT,IGCCI,       
     1           EATOM,ISCF,ITGT,ICI,GAH,ITRIP,IPEN,ETOT,ETOT1,IDOB,    
     2           NB,NAB                                                 
      COMMON /ZE/ N4                                                    
      COMMON /B/ ES(300),ER(300),EQ(300),SUM(300),ET(300),C(400)             
     1 ,GAX(300)                                                         
      COMMON /C/ ESING(300,300),U(300,300)                                  
      COMMON /D/ XJSL(5,300),IQ(2,300),INZ(2,300),W(300)                    
      COMMON /MOM/ S1(300),S2(300),S3(300)                                 
      COMMON /E/ NCYCL,EPREV,ITERR,N5                                   
      COMMON /PIMCD/ MCD,NC,NOCC,NST,LPI(300),XP(100),
     1  YP(100),BETA(5050),CPI(10000),
     2  CIPI(10000),IFROM(100),ITO(100),EP(100),POCC(100)          
  101 FORMAT (//38X,' CONFIGURATION ENERGY',37X,'TRANSITION MOMENT')    
  102 FORMAT ('  OCC.ORB.',8X,9HVIRT.ORB.,18X,5HIN EV,20X,5HCONF.)      
  120 FORMAT (5X,3HNO.,13X,3HNO.,14X,7HSINGLET,6X,7HTRIPLET,14X,3HNO.,  
     .       18X,2HSX,8X,2HSY,8X,2HSZ,7X,3HOSC/1X,132(1H-))             
  121 FORMAT (4X,I3,1X,11H   ---->   ,1X,I3,14X,F8.4,5X,F8.4,13X,I3,13X,
     .       4F10.5)                                                    
 9999 format(4f20.10)
c-----LOAD GROUND CFG. INTO IQ(*,1)                                     
      IQ(1,1)=0                                                         
      IQ(2,1)=0                                                         
      DO 8 I=1,110                                                       
      DO 8 J=I,110                                                       
      ESING(I,J)=0.                                                     
    8 ESING(J,I)=0.                                                     
      CALL LOAD(ESING,2,N5)                                             
      CALL LOAD(U,4,N5)                                                 
      JCOUNT=1                                                          
      JCO=0                                                             
      IF (NLUMO.GT.N5.OR.NHOMO.LT.1) GO TO 21                           
      WRITE(6,101)                                                      
      WRITE(6,102)                                                      
      WRITE(6,120)                                                      
      IF(NHOMO.GE.NLUMO) WRITE(6,1000)                                  
 1000 FORMAT(/' GRAND CANONICAL ENSEMBLE AVERAGING PROCEDURE:'/)        
      NMIN=MAX0(1,NLUMO-NB)                                             
      NMAX=MIN0(N5,NHOMO+NAB)                                           
      NB=NB/IDOB                                                        
      NAB=NAB/IDOB                                                      
      NMI=MAX0(2,NLUMO-NB)                                              
      NMA=MIN0(N5-1,NHOMO+NAB)                                          
      TEST=ABS(SUM(NMA)-SUM(NMA+1))                                     
      IF (TEST.LT.1.E-4) NAB=NAB-1                                      
      TEST=ABS(SUM(NMI)-SUM(NMI-1))                                     
      IF (TEST.LT.1.E-4) NB=NB-1                                        
      IF (NMAX.EQ.N5) GO TO 5                                           
      TEST=ABS(SUM(NMAX)-SUM(NMAX+1))                                   
      IF (TEST.LT.1.E-4) NMAX=NMAX-1                                    
    5 IF (NMIN.EQ.1) GO TO 6                                            
      TEST=ABS(SUM(NMIN)-SUM(NMIN-1))                                   
      IF (TEST.LT.1.E-4) NMIN=NMIN+1                                    
    6 CALL SHIFT(ESING,N5)                                             
      DO 3 I=NLUMO,NMAX                                                 
      IF(LPI(I).LT.1) GO TO 3                                           
      JJJ=0                                                             
      DO 12 JJ=NMIN,NHOMO                                               
      J=NHOMO-JJJ                                                       
      JJJ=JJJ+1                                                         
      IF(LPI(J).LT.1) GO TO 12                                          
      IF(I.LE.J) GOTO 12                                                
      DOCC=OCC(J)-OCC(I)                                                
      IF(ABS(DOCC).LE..00001) GOTO 12                                   
      ERI=SUM(I)-SUM(J)-DOCC*H(I,I,J,J)                                 
      ESI=ERI+DOCC*2.*H(I,J,J,I)                                        
      JCO=JCO+1                                                         
      IF (ESI.GE.EMAX) GO TO 12                                         
      CALL TRAMOM (SX,SY,SZ,I,J,U,X,Y,Z)                                
   16 JCOUNT=JCOUNT+1                                                   
      ESING(JCOUNT,JCOUNT)=ESI                                          
      IQ(1,JCOUNT)=J                                                    
      IQ(2,JCOUNT)=I                                                    
      O=.00379277*ESI*(SX*SX+SY*SY+SZ*SZ)                               
      S1(JCOUNT)=SX                                                     
      S2(JCOUNT)=SY                                                     
      S3(JCOUNT)=SZ                                                     
      ER(JCOUNT)=ERI                                                    
      WRITE(6,121) J,I,ESI,ERI,JCOUNT,SX,SY,SZ,O                        
      IF (JCOUNT.EQ.NQ) GO TO 14                                        
   12 CONTINUE                                                          
    3 CONTINUE                                                          
  14  CONTINUE                                                          
      IF (JCOUNT.LE.0) GO TO 21                                         
      REWIND 3                                                          
      JCO=JCOUNT-1                                                      
      DO 1 I=1,JCOUNT                                                   
    1 ESING(1,I)=0.                                                     
      DO 41 I=2,JCO                                                     
      IL=I+1                                                            
      IX=IQ(1,I)                                                        
      JX=IQ(2,I)                                                        
      DOCC=OCC(IX)-OCC(JX)                                              
      DO 4 J=IL,JCOUNT                                                  
      LX=IQ(1,J)                                                        
      MX=IQ(2,J)                                                        
      DDOCC=SQRT(DOCC*(OCC(LX)-OCC(MX)))                                
      IF (IX.NE.LX) GO TO 7                                             
      IF (JX.EQ.MX) GO TO 4                                             
    7 ESIN=-H(LX,IX,JX,MX)*DDOCC                                        
      IF (ITRIP.GT.0) WRITE (3,9999) ESIN                                    
      HILF=H(LX,MX,JX,IX)*DDOCC                                         
      ESING(I,J)=HILF+HILF+ESIN                                         
    4 CONTINUE                                                          
   41 CONTINUE                                                          
   21 CONTINUE                                                          
      IF(IGCCI.EQ.1) CALL CCI(ESING,IQ,JCOUNT)                          
      NGRO=MAX0(JCOUNT,N5)                                              
      CALL STORE(ESING,1,NGRO)                                          
      CALL SHIFT (ESING,NGRO)                                           
c-MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD- 
      NST=MIN0(JCOUNT,100)                                               
      DO 2000 I=1,NST                                                   
      IFROM(I)=LPI(IQ(1,I))                                             
 2000 ITO  (I)=LPI(IQ(2,I))                                             
c-MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD--MCD- 
      RETURN                                                            
      END                                                               
c###################################################################    
      SUBROUTINE CCI(E,IQ,N)                                            
      implicit real*8(A-H,O-Z)
c-----CONV. SINGLET GRAND CANONICAL CI-MATRIX TO CANONICAL CI-MATRIX    
      DIMENSION E(300,300),IQ(2,300)                                       
      INTEGER O                                                         
      COMMON /OCCDAT/ NHOMO,NLUMO,NPE,NPMO,OCC(300)                      
      COMMON /MOM/    S1(300),S2(300),S3(300)                              
c-----RETURN IF NOT THE CASE OF A SINGLE HALF FILLED ORBITAL            
      IF(NHOMO.NE.NLUMO) RETURN                                         
      WRITE(6,1000)                                                     
 1000 FORMAT(/' CANONICAL ENSEMBLE AVERAGING PROCEDURE:'/)              
      O=NHOMO                                                           
      DO 20 J=2,N                                                        
      L =IQ(1,J)                                                        
      K =IQ(2,J)                                                        
      OCCLK=OCC(L)-OCC(K)                                               
      DO 10 I=2,J                                                       
      L1=IQ(1,I)                                                        
      K1=IQ(2,I)                                                        
      C = 0.0                                                           
      IF( L .EQ.O .AND. L1.EQ.O )  C = C - H(K1,O,O,K)                  
      IF( K .EQ.O .AND. K1.EQ.O )  C = C - H(L1,O,O,L)                  
      IF( L1.EQ.O .AND. K .EQ.O )  C = C + H(L,O,K1,O)                  
      IF( K1.EQ.O .AND. L .EQ.O )  C = C + H(K,O,L1,O)                  
   10 E(I,J)=E(I,J)+0.25*C/SQRT(OCCLK*(OCC(L1)-OCC(K1)))                
      OSC=0.00379277*E(J,J)*(S1(J)**2+S2(J)**2+S3(J)**2)                
   20 WRITE(6,2000) L,K,E(J,J),J,S1(J),S2(J),S3(J),OSC                  
 2000 FORMAT(4X,I3,'    ---->    ',I3,14X,F8.4,26X,I3,13X,4F10.5)       
      RETURN                                                            
      END                                                               
c###################################################################    
      SUBROUTINE TRAMOM (SX,SY,SZ,I,J,U,X,Y,Z)                          
      implicit real*8(A-H,O-Z)
c     PROGRAMMED BY H.BAUMANN                                           
c     TRANSITION MOMENT BETWEEN GROUND CFG. AND SINGLY EXC. CFG. J->I   
c-----MODIFIED TO INCLUDE SP-TERMS BY JS-L                              
      DIMENSION X(100),Y(100),Z(100),U(300,300)                              
      COMMON /A/ N,NH,NE,NQ,ECORE,ICOUNT,JCOUNT,EMAX,IPRNT,IGCCI,       
     1           EATOM,ISCF,ITGT,ICI,GAH,ITRIP,IPEN,ETOT,ETOT1,IDOB     
      COMMON /E/ NCYCL,EPREV,ITERR,N5,NT,N2S                            
      COMMON /SPTERM/ SP(100)                                            
      COMMON /OCCDAT/ NHOMO,NLUMO,NPE,NPMO,OCC(300)                      
      FACTOR=4.8033*SQRT(2.)                                            
      SX=0.                                                             
      SY=0.                                                             
      SZ=0.                                                             
      IF (1.GT.NH) GO TO 16                                             
      DO 14 K=1,NH                                                      
      UIJ=U(K,I)*U(K,J)                                                 
      SX=SX+UIJ*X(K)                                                    
      SY=SY+UIJ*Y(K)                                                    
      SZ=SZ+UIJ*Z(K)                                                    
   14 CONTINUE                                                          
   16 IF (1.GT.N) GO TO 17                                              
      DO 15 NMULT=1,4                                                   
      DO 15 K=N2S,NT                                                    
      KK=K+N*(NMULT-1)                                                  
      UIJ=U(KK,I)*U(KK,J)                                               
      SX=SX+UIJ*X(K)                                                    
      SY=SY+UIJ*Y(K)                                                    
      SZ=SZ+UIJ*Z(K)                                                    
   15 CONTINUE                                                          
c-----ADD SP-CONTRIBUTIONS                                              
      DO 99 IN=1,N                                                      
      NS=NH+IN                                                          
      NPX=NS+N                                                          
      NPY=NPX+N                                                         
      NPZ=NPY+N                                                         
      USI=U(NS,I)                                                       
      USJ=U(NS,J)                                                       
      SPN=SP(NS)                                                        
      SX=SX + (USI*U(NPX,J)+USJ*U(NPX,I)) * SPN                         
      SY=SY + (USI*U(NPY,J)+USJ*U(NPY,I)) * SPN                         
   99 SZ=SZ + (USI*U(NPZ,J)+USJ*U(NPZ,I)) * SPN                         
   17 FAC=SQRT(OCC(J)-OCC(I))*FACTOR                                    
      SX=SX*FAC                                                         
      SY=SY*FAC                                                         
      SZ=SZ*FAC                                                         
      RETURN                                                            
      END                                                                           
c#######################################################################
      SUBROUTINE MCDB(IFLAG)                                            
      implicit real*8(A-H,O-Z)
      COMMON /PIMCD/ MCD,N,NOCC,NST,LPI(300),X(100),Y(100),BETA(5050),     
     1   C(10000),CI(10000),FROM(100),TO(100),E(100),OCC(100)                 
      LOGICAL P(5)                                                      
      INTEGER FROM,TO                                                   
      EXTERNAL MATPRT                                                   
      IF(MCD.LT.0) RETURN                                               
      DO 10 I=1,5                                                       
   10 P(I)=.FALSE.                                                      
      IF(IFLAG.EQ.2) GOTO 100                                           
      DO 15 I=1,MCD                                                     
   15 P(6-I)=.TRUE.                                                     
      WRITE(6,1000)                                                     
 1000 FORMAT(/' -------------------------------------------------',                          
     .       /' CALCULATION OF MCD B TERMS FOR PI-PI* TRANSITIONS',       
     .       /' -------------------------------------------------',      
     .     //' PI MOLECULAR ORBITAL COEFFICIENTS')                      
      CALL MATPRT(C,N,0)                                                
      WRITE(6,1500) (OCC(I),I=1,N)                                      
 1500 FORMAT(1X//' OCCUPATION NUMBERS'//5X,16F8.4)                      
      WRITE(6,2000)                                                     
 2000 FORMAT(1X//' PI TYPE BETA MATRIX FOR MAGNETIC DIPOLE INTEGRALS')  
      CALL MATPRT(BETA,N,1)                                             
  100 CALL BTERM(N,NOCC,NST,X,Y,BETA,C,CI,FROM,TO,E,MATPRT,P,OCC)       
      RETURN                                                            
      END                                                               
c#######################################################################
      SUBROUTINE BTERM(N,NOCC,NST,X,Y,BETA,C,CI,FROM,TO,E,MATPRT,P,OCC) 
c     Adapted from Computer Program BTERM, QCPE Bull. 1, 37 (1981).
c     This routine was originally developed for use with a 
c     PPP program operating with one-dimensional arrays.   (JS-L)
      implicit real*8(A-H,O-Z)
      LOGICAL P(5)                                                      
      INTEGER FROM(1),TO(1),F,FI,FF,GF                                  
      DIMENSION X(1),Y(1),BETA(1),C(1),CI(1),E(1),OCC(1),               
     1 DEXMO(1),DEYMO(1),DMZMO(1),DEXCFG(1),DEYCFG(1),DMZAO(1),         
     1 AXI(1),AYI(1)                                                    
      COMMON/C/ INDEX(100,100),BG(10000),BF(10000),B(100),AI(100),            
     1  AZI(100),SBG(100),SBF(100),DEXST(5050),DEYST(5050),                
     1 DMZST(5050),DMZCFG(5050)                                         
      EQUIVALENCE (DEXST(1),DEXMO(1)),(DEYST(1),DEYMO(1)),              
     .        (DMZST(1),DMZMO(1)),(BG(1),DEXCFG(1)),(BF(1),DEYCFG(1)),  
     .        (DMZCFG(1),DMZAO(1)),(SBG(1),AXI(1)),(SBF(1),AYI(1))      
      DATA F1/0.13124/,F2/4.8033/,F3/8.0658/,F4/0.00379277/             
      S(I,J)=SIGN(1.,FLOAT(I-J))
      PI180 = 3.141592654/180.0                                
      IF(N.GT.100.OR.NST.GT.100)GOTO 260                                  
      MAX=MAX0(N,NST)                                                   
      DO 10 J=1,MAX                                                     
      DO 10 I=1,J                                                       
      INDEX(I,J)=I+(J*(J-1))/2                                          
   10 INDEX(J,I)=INDEX(I,J)                                             
      SQR2=SQRT(2.0)                                                    
c-----------------------------------------------------------------------
c   ELEMENTS OF MAGNETIC DIPOLE MOMENT OPERATOR OVER AO*S, Z-COMPONENT  
c   (ACCORDING TO LINDERBERG*S FORMULA)                                 
c-----------------------------------------------------------------------
      MYNY=0                                                            
      DO 20 MY=1,N                                                      
            DO 20 NY=1,MY                                               
                  MYNY=MYNY+1                                           
   20             DMZAO(MYNY)=-BETA(MYNY)*(X(MY)*Y(NY)-X(NY)*Y(MY))*F1  
      IF(P(1)) WRITE(6,1000)                                            
 1000 FORMAT(//'Z-COMPONENT OF MAGNETIC DIPOLE OPERATOR',             
     .         /' IN AO BASIS  (BOHR-MAGNETONS)')                       
      IF(P(1)) CALL MATPRT(DMZAO,N,1)                                   
c-----------------------------------------------------------------------
c   ELEMENTS OF DIPOLE OPERATORS OVER MO*S                              
c   (ELECTRIC DIPOLE OPERATOR ACCORDING TO DIPOLE LENGTH EXPRESSION)    
c-----------------------------------------------------------------------
      DO 70 I=1,N                                                       
            II=N*(I-1)                                                  
            DO 40 NY=1,N                                                
                  A=0.0                                                 
                  DO 30 MY=1,N                                          
                       MYNY=INDEX(MY,NY)                               
   30                   A=A+DMZAO(MYNY)*S(MY,NY)*C(MY+II)               
   40             AI(NY)=A                                              
            DO 60 J=1,I                                                 
                  JJ=N*(J-1)                                            
                  AX=0.0                                                
                  AY=0.0                                                
                  AZ=0.0                                                
                  DO 50 MY=1,N                                          
                        CC=C(MY+II)*C(MY+JJ)                            
                        AX=AX-CC*X(MY)                                  
                        AY=AY-CC*Y(MY)                                  
   50                   AZ=AZ+C(MY+JJ)*AI(MY)                           
                 JI=INDEX(J,I)                                         
                 DEXMO(JI)=AX * F2                                     
                 DEYMO(JI)=AY * F2                                     
   60            DMZMO(JI)=AZ                                          
   70       CONTINUE                                                    
      IF(P(2)) WRITE(6,1100)                                            
 1100 FORMAT(//'X-COMPONENT OF ELECTRIC DIPOLE OPERATOR',             
     .         /' IN MO BASIS  (DEBYES)')                               
      IF(P(2)) CALL MATPRT(DEXMO,N,1)                                   
      IF(P(2)) WRITE(6,1200)                                            
 1200 FORMAT(//' Y-COMPONENT OF ELECTRIC DIPOLE OPERATOR',            
     .          /' IN MO BASIS  (DEBYES)')                              
      IF(P(2)) CALL MATPRT(DEYMO,N,1)                                   
      IF(P(2)) WRITE(6,1300)                                            
 1300 FORMAT(//' Z-COMPONENT OF MAGNETIC DIPOLE OPERATOR'             
     .          /' IN MO BASIS  (BOHR-MAGNETONS)')                      
      IF(P(2)) CALL MATPRT(DMZMO,N,1)                                   
c-----------------------------------------------------------------------
c   ELEMENTS OF DIPOLE OPERATORS OVER GROUND AND SINGLY EXCITED         
c   CONFIGURATIONS (AS DEFINED BY THE LISTS *FROM* AND *TO*)            
c-----------------------------------------------------------------------
      DO 150 I=1,NST                                                    
             I1=FROM(I)                                                 
            I2=TO(I)                                                   
             DO 150 J=1,I                                               
                   IF(J.LT.I) GO TO 110                                
c-----DIAGONAL ELEMENTS                                                 
                   IF(J.GT.1) GO TO 90                                 
c- - -GROUND CONFIGURATION                                              
                   DEXCFG(1)=0.0                                       
                   DEYCFG(1)=0.0                                       
                   DMZCFG(1)=0.0                                       
                   DO 80 K=1,NOCC                                      
                         KK=INDEX(K,K)                                 
                         DEXCFG(1)=DEXCFG(1)+2.0*DEXMO(KK)*OCC(K)      
   80                    DEYCFG(1)=DEYCFG(1)+2.0*DEYMO(KK)*OCC(K)      
                   GO TO 150                                           
c- - -EXCITED CONFIGURATION                                             
   90              II=INDEX(I,I)                                       
                   I1I1=INDEX(I1,I1)                                   
                   I2I2=INDEX(I2,I2)                                   
                   DEXCFG(II)=DEXCFG(1)-DEXMO(I1I1)+DEXMO(I2I2)        
                   DEYCFG(II)=DEYCFG(1)-DEYMO(I1I1)+DEYMO(I2I2)        
  100              DMZCFG(II)=.0                                       
                   GO TO 150                                           
c-----OFF-DIAGONAL ELEMENTS                                             
  110               IF(J.GT.1) GO TO 120                                
c- - -EXCITED CONFIGURATION / GROUND CONFIGURATION                      
                   I1I2=INDEX(I1,I2)                                   
                   AX=SQR2*DEXMO(I1I2)                                 
                   AY=SQR2*DEYMO(I1I2)                                 
                   AZ=SQR2*DMZMO(I1I2)                                 
                   GO TO 140                                           
c- - -EXCITED CONFIGURATION / EXCITED CONFIGURATION                     
  120              J1=FROM(J)                                          
                   J2=TO(J)                                            
                   AX=.0                                               
                   AY=.0                                               
                   AZ=.0                                               
                   IF(I1.NE.J1) GO TO 130                              
                   I2J2=INDEX(I2,J2)                                   
                   AX=DEXMO(I2J2)                                      
                   AY=DEYMO(I2J2)                                      
                   AZ=DMZMO(I2J2)*S(I2,J2)                             
                   GO TO 140                                           
  130               IF(I2.NE.J2) GO TO 140                              
                   J1I1=INDEX(J1,I1)                                   
                   AX=-DEXMO(J1I1)                                     
                   AY=-DEYMO(J1I1)                                     
                   AZ=-DMZMO(J1I1)*S(J1,I1)                            
  140               IJ=INDEX(I,J)                                       
                   DEXCFG(IJ)=AX                                       
                   DEYCFG(IJ)=AY                                       
                   DMZCFG(IJ)=AZ                                       
  150               CONTINUE                                            
      IF(P(3)) WRITE(6,1400)                                            
 1400 FORMAT(//'X-COMPONENT OF ELECTRIC DIPOLE OPERATOR',              
     .        /' IN CONFIGURATION BASIS  (DEBYES)')                     
      IF(P(3)) CALL MATPRT(DEXCFG,NST,1)                                
      IF(P(3)) WRITE(6,1500)                                            
 1500 FORMAT(//' Y-COMPONENT OF ELECTRIC DIPOLE OPERATOR',            
     .          /' IN CONFIGURATION BASIS  (DEBYES)')                   
      IF(P(3)) CALL MATPRT(DEYCFG,NST,1)                                
      IF(P(3)) WRITE(6,1600)                                            
 1600 FORMAT(//' Z-COMPONENT OF MAGNETIC DIPOLE OPERATOR',            
     .          /' IN CONFIGURATION BASIS  (BOHR-MAGNETONS)')           
      IF(P(3)) CALL MATPRT(DMZCFG,NST,1)                                
c-----------------------------------------------------------------------
c   ELEMENTS OF DIPOLE OPERATORS OVER CI-STATES                         
c-----------------------------------------------------------------------
      DO 190 I=1,NST                                                    
            II=NST*(I-1)                                               
            DO 170 L=1,NST                                             
                   AX=.0                                               
                   AY=.0                                               
                   AZ=.0                                               
                   DO 160 K=1,NST                                      
                          KL=INDEX(K,L)                                
                          AX=AX+DEXCFG(KL)*CI(K+II)                    
                          AY=AY+DEYCFG(KL)*CI(K+II)                    
  160                     AZ=AZ+DMZCFG(KL)*S(K,L)*CI(K+II)             
                   AXI(L)=AX                                           
                   AYI(L)=AY                                           
  170              AZI(L)=AZ                                           
            DO 190 J=1,I                                               
                   JJ=NST*(J-1)                                        
                   AX=.0                                               
                   AY=.0                                               
                   AZ=.0                                               
                   DO 180 L=1,NST                                      
                          AX=AX+CI(L+JJ)*AXI(L)                        
                          AY=AY+CI(L+JJ)*AYI(L)                        
  180                     AZ=AZ+CI(L+JJ)*AZI(L)                        
                   IJ=INDEX(I,J)                                       
                   DEXST(IJ)=AX                                        
                   DEYST(IJ)=AY                                        
                   DMZST(IJ)=AZ                                        
  190               CONTINUE                                            
      IF(P(4)) WRITE(6,1700)                                            
 1700 FORMAT(//'X-COMPONENT OF ELECTRIC DIPOLE OPERATOR',             
     .         /' IN STATE BASIS  (DEBYES)')                            
      IF(P(4)) CALL MATPRT(DEXST,NST,1)                                 
      IF(P(4)) WRITE(6,1800)                                            
 1800 FORMAT(//' Y-COMPONENT OF ELECTRIC DIPOLE OPERATOR',            
     .          /' IN STATE BASIS  (DEBYES)')                           
      IF(P(4)) CALL MATPRT(DEYST,NST,1)                                 
      IF(P(4)) WRITE(6,1900)                                            
 1900 FORMAT(//' Z-COMPONENT OF MAGNETIC DIPOLE OPERATOR',            
     .          /' IN STATE BASIS  (BOHR-MAGNETONS)')                   
      IF(P(4)) CALL MATPRT(DMZST,NST,1)                                 
c-----------------------------------------------------------------------
c   EXPANSION OF B-TERMS  B(F) = SUMMA(I) ( BG(I,F) + BF(I,F) )         
c-----------------------------------------------------------------------
      DO 240 F=1,NST                                                    
            FF=NST*(F-1)                                               
            GF=INDEX(1,F)                                              
            SBG(F)=.0                                                  
            SBF(F)=.0                                                  
            DO 230 I=1,NST                                             
                   FI=INDEX(F,I)                                       
                   IG=INDEX(I,1)                                       
c- - -CONTRIBUTION DUE TO MIXING OF I*TH STATE WITH GROUND STATE        
                   IF(I.NE.1) GO TO 200                                
                   BG(1+FF)=.0                                         
                   GO TO 210                                           
  200               EIE1=E(I)-E(1)                                      
                   IF(ABS(EIE1).LT.0.000001) GOTO 270                  
                   BG(I+FF)= DMZST(IG)                                 
     .                   * (DEXST(GF)*DEYST(FI)-DEXST(FI)*DEYST(GF))    
     .                   / ( EI E1  * F3)                               
                   SBG(F)=SBG(F)+BG(I+FF)                              
c- - -CONTRIBUTION DUE TO MIXING OF I*TH STATE WITH EXC.STATE F         
  210               IF(I.NE.F) GO TO 220                                
                   BF(F+FF)=.0                                         
                   GO TO 230                                           
  220               EIEF=E(I)-E(F)                                      
                   IF(ABS(EIEF).LT.0.000001) GOTO 270                  
                   BF(I+FF)= DMZST(FI)*S(F,I)                          
     .                   * (DEXST(GF)*DEYST(IG)-DEXST(IG)*DEYST(GF))    
     .                   / ( EI EF * F3)                                
                   SBF(F)=SBF(F)+BF(I+FF)                              
  230               CONTINUE                                            
  240        B(F)= SBG(F) + SBF(F)                                      
      IF(P(5)) WRITE(6,2000)                                            
 2000 FORMAT(//' CONTRIBUTIONS TO MCD B-TERMS B(F) DUE TO MIXING',     
     .          ' OF I''TH STATE WITH GROUNDSTATE'/' I/F')              
      IF(P(5)) CALL MATPRT(BG,NST,0)                                    
      IF(P(5)) WRITE(6,2100)                                            
 2100 FORMAT(//' CONTRIBUTIONS TO MCD B-TERMS B(F) DUE TO MIXING',   
     .            ' OF I''TH STATE WITH F''TH STATE'/' I/F')            
      IF(P(5)) CALL MATPRT(BF,NST,0)                                    
c-----PRINT B-TERMS                                                     
      WRITE(6,2200)                                                     
 2200 FORMAT(//' WAVENUMBERS IN KK = 1000CM-1,'/                       
     .            ' MCD B-TERMS IN (BOHR-MAGNETON DEBYE**2)/KK,'/                   
     .            ' ELECTRIC TRANSITION MOMENTS IN DEBYES'/
     .            ' TRANSITION MOMENT ANGLES IN DEG'/               
     .            ' OSCILLATOR STRENGTHS  OSC'//                            
     .'   F    W         BG     +    BF     =     B         M(X)',      
     .'        M(Y)     PHI(X-->Y)     OSC'/1X,100('-'))                              
      DO 250 F=2,NST                                                    
             GF=INDEX(1,F)                                            
             W= F3 * (E(F)-E(1)) 
             DE2= DEXST(GF)**2 + DEYST(GF)**2                        
             O= F4 * (E(F)-E(1)) * DE2
             PHI= 0.0
             IF ( SQRT(DE2).LE.0.00001 ) GOTO 250
             PHI= 90.0
             IF ( ABS(DEXST(GF)).LE.0.0000001 ) GOTO 250
             PHI= ATAN(DEYST(GF)/DEXST(GF)) / PI180
  250        WRITE(6,2500)F,W,SBG(F),SBF(F),B(F),DEXST(GF),
     .       DEYST(GF),PHI,O
 2500 FORMAT(1X,I3,F8.3,65F12.6)                                        
c-----PLOT SPECTRUM WITH INDICATION OF MCD B-SIGNS                      
      CALL SPRT1(1,B)                                                   
      RETURN                                                            
c-----------------------------------------------------------------------
c   TROUBLE SECTION                                                     
c-----------------------------------------------------------------------
  260 WRITE(6,3000) N,NST                                               
 3000 FORMAT(1X//' N =',I3,', NST=',I3,                                
     .           /' DIMENSIONS TOO LARGE FOR SUBROUTINE BTERM',         
     .           /' RECOMPILE WITH APPROPRIATE ARRAY DIMENSIONS'/1H1)   
      RETURN                                                            
  270 WRITE(6,4000) I,F                                                 
 4000 FORMAT(1X//' TERMINATION DUE TO NEAR DEGENERACY OF STATES NO.',  
     .           I3,' AND',I3/1H1)                                     
      RETURN                                                            
      END                                                               
c#######################################################################
       SUBROUTINE MATPRT(A,N,MS)                                         
      implicit real*8(A-H,O-Z)
c ---------------------------------------------------------------------- 
c    PRINT N * N MATRIX  A, STORAGE MODE MS ( 0 OR 1 )                   
c---------------------------------------------------------------------- 
      DIMENSION A(1)                                                    
      DATA LINE/4H----/                                                 
      DO 20 J1=1,N,16                                                   
           J2=MIN0(N,J1+15)                                            
           WRITE(6,1000) (J,J=J1,J2)                                   
           WRITE(6,1500) ((LINE,L=1,2),J=J1,J2)                        
           I1=1                                                        
           IF(MS.EQ.1) I1=J1                                           
           DO 20 I=I1,N                                                
                 IF(MS.EQ.1) GO TO 10                                  
                 K1=N*(J1-1)+I                                         
                 K2=N*(J2-1)+I                                         
                 WRITE(6,2000) I,(A(K),K=K1,K2,N)                      
                 GO TO 20                                              
   10             II=I*(I-1)/2                                          
                 K1=II+I1                                              
                 K2=II+MIN0(I,J2)                                      
                 WRITE(6,2000) I,(A(K),K=K1,K2)                        
   20             CONTINUE                                              
 1000 FORMAT(/16(5X,I3))                                             
 1500 FORMAT(1X,4H----,16(2A4))                                         
 2000 FORMAT(I4,16F8.4)                                                 
      RETURN                                                            
      END                                                               






