PROGRAM FCI
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /CFULLR/ CORE,THR,THRRES
      COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
     >               ,INTOFF(8),ITYPEA(256)
     >               ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
     >               ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
     >               ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
     >               ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
      COMMON /BIG/ Q(1)
      INTEGER IQ(1)
      EQUIVALENCE (Q(1),IQ(1))
      COMMON /CHAMIL/ ISYM,NNCI
      COMMON /TAPES / INP,IOUT
      DIMENSION MULTT(8,8)
      DATA MULTT /1,2,3,4,5,6,7,8, 2,1,4,3,6,5,8,7,
     2            3,4,1,2,7,8,5,6, 4,3,2,1,8,7,6,5,
     4            5,6,7,8,1,2,3,4, 6,5,8,7,2,1,4,3,
     6            7,8,5,6,3,4,1,2, 8,7,6,5,4,3,2,1/
      INP = 5
      IOUT = 6
      WRITE (IOUT,35)
35    FORMAT(' PROGRAM * FCI  (Determinant based Full CI)',
     >       5X,'Author: P.J. Knowles, 1984')
C...  SET UP DEFAULT OR INITIAL VALUES
      MEMORY = 100 000
      IPRTIM = -1
      IPRINT = -1
      IPRDIA = -1
      IPRHAM = -1
      IOPTIO = 0
      MAXIT = 25
      THR = 1D-5
      THRRES = 0.05
      NROOT = 1
      NACT = 0
      NELEC = 0
      MS2 = 0
      ISYM = 1
      DO 2 I=1,8
      NOCC(I) = 0
C...  D2H MULTIPLICATION TABLE
      DO 2 J=1,8
2     MULT(J,I) = MULTT(J,I)
      DO 3 I=1,256
3     ITYPEA(I)=1
C
C...  READ INPUT
      CALL INPDAT
C
C...  INITIALIZE MEMORY HANDLING
      I = INICOR(MEMORY)
C
C...  NUMBERS OF ALPHA, BETA ELECTRONS
      NA = (NELEC+MS2)/2
      NB = (NELEC-MS2)/2
      DO 36 I=1,NACT
36    NOCC(ITYPEA(I)) = NOCC(ITYPEA(I))+1
      CALL INITC
C
      WRITE (IOUT,25) ISYM
25    FORMAT(/' Space symmetry:',T30,I3)
      WRITE (IOUT,21) MAXIT,THR,THRRES,NROOT
21    FORMAT(
     4' Maximum iterations:',   T30,I3/
     5' Convergence threshold:',T30,F12.7/
     8' Output threshold:',     T30,F12.7/
     A' Number of roots:',      T30,I3)
      NNCI = NCI(ISYM)
      IV1 = ICORR(NNCI)
      IV2 = ICORR(NNCI)
      CALL DAVIDS (Q(IV1),Q(IV2),IOPTIO)
      CALL CORLSR(IV1)
C
      CALL ENDCOR
      END
      SUBROUTINE FCMXM (D,Z,E,NK,NIJ,NOPCNT)
CSUBR MATRIX MULTIPLICATION KERNEL FOR FULL CI PROGRAM
CSUBR  E = D * Z
CSUBR  Z IS INTEGRAL MATRIX AND NOT SPARSE, DIMENSION NIJ=NUMBER OF
CSUBR            ORBITAL PAIRS OF GIVEN SYMMETRY
CSUBR  D,E LEADING DIMENSION NK=NUMBER OF DETERMINANTS OF A SYMMETRY
CSUBR            BLOCK, OR SUB-BATCH OF THIS
CSUBR  D IS VERY SPARSE IN 1ST ITERATION OF CI, GRADUALLY RISES TO
CSUBR             USUALLY 50-80% POPULATED
CSUBR OPTIMUM FORM OF THIS ROUTINE IS DEPENDENT ON MACHINE ARCHITECTURE
CSUBR USUALLY NK>NIJ, SO BETTER TO VECTORISE ALONG NK, BUT SOMETIMES
CSUBR BETTER TO EXPLOIT SPARSITY OF D
CSUBR ROUTINE SHOULD RETURN NOPCNT -- NUMBER OF FLOATING POINT OPS --
CSUBR IF ADDITIONAL TIMING INFORMATION IS REQUIRED
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NHALF=30)
      DIMENSION D(NK,NIJ),Z(NIJ,NIJ),E(NK,NIJ)
c     pmwg replaced all code by a dgemm in Jan 2014
      call dgemm('n','n',nk,nij,nij,1d0,d,nk,z,nij,0d0,e,nk)
c      DO 1 J=1,NK*NIJ
c1     E(J,1) = 0D0
C...  EVALUATE SPARSITY OF D
c      IDS = 0
c      DO 7 J=1,NK*NIJ
c7     IF (D(J,1).NE.0D0) IDS=IDS+1
c      IZS = NIJ**2
c      IF (IDS*(NIJ+NHALF).LT.IZS*(NK+NHALF)) THEN
C...    USE SPARSITY OF D
c        NOPCNT = NOPCNT + 2*IDS*NIJ
C        CALL MXMA (Z,1,NIJ, D,NK,1, E,NK,1, NIJ,NIJ,NK)
c        DO 2 K=1,NK
c        DO 2 J=1,NIJ
c        IF (D(K,J).EQ.0D0) GOTO 2
c        DO 21 I=1,NIJ
c21      E(K,I) = E(K,I) + D(K,J)*Z(J,I)
c2       CONTINUE
c      ELSE
C...    VECTORISE ALONG D
c        NOPCNT = NOPCNT + 2*IZS*NK
C        CALL MXMA (D,1,NK, Z,NIJ,1, E,1,NK, NK,NIJ,NIJ)
c        DO 3 I=1,NIJ
c        DO 3 J=1,NIJ
C        IF (Z(J,I).EQ.0D0) GOTO 3
c        DO 31 K=1,NK
c31      E(K,I) = E(K,I) + D(K,J)*Z(J,I)
c3       CONTINUE
c      END IF
      RETURN
      END
      SUBROUTINE INPDAT
CSUBR DATA INPUT FOR PORTABLE FULL CI PROGRAM
CSUBR NAMELIST INPUT MAY REQUIRE ADJUSTMENT ON SOME MACHINES
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /CFULLR/ CORE,THR,THRRES
      COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
     >               ,INTOFF(8),ITYPEA(256)
     >               ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
     >               ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
     >               ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
     >               ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
      COMMON /BIG/ Q(1)
      INTEGER IQ(1)
      EQUIVALENCE (Q(1),IQ(1))
      COMMON /CHAMIL/ ISYM,NNCI
      COMMON /TAPES / INP,IOUT
      INTEGER ORBSYM(255)
      EQUIVALENCE (ORBSYM(1),ITYPEA(1))
      NAMELIST /FCI   / NORB,NELEC,MS2,ISYM,ORBSYM
     >                 ,IPRINT,IPRDIA,IPRHAM,IPRTIM
     >                         ,INT,MEMORY,CORE
     >                 ,MAXIT,THR,THRRES,NROOT
      INT = -1
      READ (INP,FCI)
      IF (INT.NE.-1) INP=INT
      NACT = NORB
      RETURN
      END
      SUBROUTINE INPINT (Z1,Z2)
CSUBR INTEGRAL INPUT FOR PORTABLE FULL CI PROGRAM
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /CFULLR/ CORE,THR,THRRES
      COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
     >               ,INTOFF(8),ITYPEA(256)
     >               ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
     >               ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
     >               ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
     >               ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
      COMMON /BIG/ Q(1)
      INTEGER IQ(1)
      EQUIVALENCE (Q(1),IQ(1))
      COMMON /TAPES / INP,IOUT
      DIMENSION Z2(*),Z1(*)
cpmg
c         rs = 2.0d0
c         R = sqrt(dfloat(na+nb))/2 * rs
         R = 0.25d0
         write(*,*) 'Scaling integrals for R =',R
cpmg
      CALL FZERO (Z1,NPAIR(1))
      CALL FZERO (Z2,NINT2)
1     READ (INP,*,END=99) Z,I,J,K,L
      IF (MAX(I,J,K,L).GT.NACT) GOTO 1
      IF (K.NE.0) THEN
C...    TWO ELECTRON INTEGRAL
cpmg
           z = z/R
cpmg
        ISYMIJ = MULT(ITYPEA(I),ITYPEA(J))
        IJ = IQ(INTADT+I-1+(J-1)*NACT)
        KL = IQ(INTADT+K-1+(L-1)*NACT)
        Z2(IJ+(KL-1)*NPAIR(ISYMIJ)+INTOFF(ISYMIJ)) = Z
        Z2(KL+(IJ-1)*NPAIR(ISYMIJ)+INTOFF(ISYMIJ)) = Z
      ELSE IF (I.NE.0) THEN
C...    ONE ELECTRON INTEGRAL
cpmg
           z = z/R/R
cpmg
        IJ = IQ(INTADT+I-1+(J-1)*NACT)
        Z1(IJ) = Z
      ELSE
C...    CORE REPULSION ENERGY
        CORE = Z
      END IF
      GOTO 1
99    RETURN
      END
      SUBROUTINE DAVIDS (V1,V2,IOPTIO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /TAPES / INP,IOUT
C
C...  DAVIDSON DIAGONALISER, SINGLE ROOT
C...  WARNING, NOT NECESSARILY SAFE WITH DEGENERACIES, ETC.
      PARAMETER (IFILD=92,IFIL1=93,IFIL2=94,MAXITR=200)
      COMMON /CFULLR/ CORE,THR,THRRES
      COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
     >               ,INTOFF(8),ITYPEA(256)
     >               ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
     >               ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
     >               ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
     >               ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
      COMMON /CHAMIL/ ISYMH,N
      DIMENSION V1(N),V2(N)
      DIMENSION ZM(MAXITR,MAXITR),ZV(MAXITR,MAXITR),R(MAXITR)
C
      TRIAL = 0
      RESULT = 0
      MAXITS = MIN(MAXIT,MAXITR,N)
      CALL FCDIAG (V2)
      IF (IOPTIO.GE.5) CALL OUTVEC (V2,N,'DIAGONAL ELEMENTS')
C
      OPEN (IFILD,STATUS='SCRATCH',FORM='UNFORMATTED')
      OPEN (IFIL1,STATUS='SCRATCH',FORM='UNFORMATTED')
      OPEN (IFIL2,STATUS='SCRATCH',FORM='UNFORMATTED')
      REWIND IFILD
      REWIND IFIL1
      REWIND IFIL2
C
      WRITE (IFILD) V2
      REWIND IFILD
      CALL TRIALV (V1,V2)
C
      WRITE (IOUT,'(/'' It  Tr    CPU'',2X,''Convergence'',6X,
     >''Energy''/)')
      DO 190 IT=1,MAXITS
      WRITE (IFIL1) V1
C...  CALL THE CI PROGRAM V2 = H . V1  + V2
      DO 30 I=1,N
30    V2(I)=0.0D0
      CALL FSIGMA (V1,V2)
      WRITE (IFIL2) V2
C
      REWIND IFIL1
      DO 50 JT=1,IT
      READ (IFIL1) V1
      ZM(JT,IT) = 0.0D0
      DO 40 I=1,N
40    ZM(JT,IT) = ZM(JT,IT) + V1(I)*V2(I)
50    ZM(IT,JT) = ZM(JT,IT)
      REWIND IFIL1
      REWIND IFIL2
C
C...  SMALL MATRIX DIAGONALISER
      DO 60 JT=1,IT
      DO 60 KT=1,IT
60    ZV(KT,JT) = ZM(KT,JT)
      CALL DIAG2 (MAXITR,IT,R,ZV)
C
      TEST = 0
      DO 61 IROOT=1,MIN(IT,NROOT)
      TES = ABS(ZV(IT,IROOT))
      IF (TES.LT.TEST) GOTO 61
      ITRACK = IROOT
      TEST = TES
61    CONTINUE
      IF (IT.EQ.1) EREF=R(1)
      WRITE (IOUT,70) IT,ITRACK,SECOND(),TEST
     >,(R(IROOT)+CORE,IROOT=1,MIN(NROOT,IT))
70    FORMAT(I3,I4,F7.1,F12.8,10F20.9)
      IF (TEST.LT.THR) GOTO 210
C
C...  MAKE RESIDUAL
      DO 80 I=1,N
80    V1(I) = 0.0D0
      DO 100 JT=1,IT
      READ (IFIL1) V2
      DO 90 I=1,N
90    V1(I) = V1(I) + (-R(ITRACK)*ZV(JT,ITRACK))*V2(I)
      READ (IFIL2) V2
      DO 100 I=1,N
100   V1(I) = V1(I) + ZV(JT,ITRACK) * V2(I)
      READ (IFILD) V2
      IF (IT.EQ.1) THEN
      DO 110 JJ = 1,N
      IF (ABS(V2(JJ)-EREF).LT.1D-9) V2(JJ)=1D20
110   CONTINUE
      END IF
      REWIND IFILD
      DO 130 I=1,N
130   V1(I) = V1(I) / (V2(I)-R(1)+1.0D-10)
C
C...  ORTHOGONALISE
140   REWIND IFIL1
      DO 160 JT=1,IT
      READ (IFIL1) V2
      ZZ = 0.0D0
      DO 150 I=1,N
150   ZZ = ZZ + V1(I)*V2(I)
      DO 160 I=1,N
160   V1(I) = V1(I) + (-ZZ) * V2(I)
      ZZ = 0.0D0
      DO 170 I=1,N
170   ZZ = ZZ + V1(I)**2
      IF (ZZ.LT.1D-20) GOTO 210
      ZZ = 1.0 / SQRT(ZZ)
      DO 180 I=1,N
180   V1(I) = V1(I) * ZZ
C...  REORTHOGONALISE IF A LOT HAS BEEN ANNIHILATED
      IF (ZZ.GT.1D3) GOTO 140
C
190   CONTINUE
      IT = MAXITS
      WRITE (IOUT,200)
200   FORMAT(/1X,'*** Convergence not achieved in max iterations')
C
210   DO 280 IROOT=1,NROOT
      REWIND IFIL1
      DO 220 I=1,N
220   V1(I) = 0.0D0
      DO 230 JT=1,IT
      READ (IFIL1) V2
      DO 230 I=1,N
230   V1(I) = V1(I) + ZV(JT,IROOT) * V2(I)
      WRITE (IOUT,'('' State'',I2,5X,''Energy'',F23.12)')
     > IROOT,R(IROOT)+CORE
      IF (IROOT.EQ.1) WRITE (IOUT,'('' Correlation energy'',F23.12)')
     > R(1)-EREF
      WRITE (IOUT,'(/'' Final CI vector''/)')
      CALL FCVECO (V1,ISYMH,THRRES)
280   CONTINUE
C
      CLOSE (IFILD,STATUS='DELETE')
      CLOSE (IFIL1,STATUS='DELETE')
      CLOSE (IFIL2,STATUS='DELETE')
C
      RETURN
      END
      SUBROUTINE TRIALV (V1,V2)
CSUBR SUBROUTINE TO GENERATE TRIAL CI VECTOR
CSUBR COULD BE REPLACED BY E.G. READING VECTOR IF REQUIRED
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /TAPES / INP,IOUT
      COMMON /CFULLR/ CORE,THR,THRRES
      COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
     >               ,INTOFF(8),ITYPEA(256)
     >               ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
     >               ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
     >               ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
     >               ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
      COMMON /BIG/ Q(1)
      INTEGER IQ(1)
      EQUIVALENCE (Q(1),IQ(1))
      COMMON /CHAMIL/ ISYMH,N
      DIMENSION IMIN(64)
      DIMENSION V1(N),V2(N)
      TRIAL = 0
      IF (NINT(TRIAL).EQ.0) THEN
C...  FIND THE DETERMINANT WITH LOWEST ENERGY
      VMIN = 1D20
      DO 10 I=1,N
      IF (V2(I).GT.VMIN) GOTO 10
      IMIN(1) = I
      VMIN = V2(I)
10    CONTINUE
C...  CONSTRUCT A SPIN EIGENFUNCTION CONTAINING THIS DETERMINANT
      CALL FCSPAD (V1,IMIN(1))
      WRITE (IOUT,20)
20    FORMAT(/1X,'Initial configuration generated:')
      NREF = 0
      DO 40 I=1,N
      IF (V1(I).EQ.0.0D0) GOTO 40
      WRITE (IOUT,30) I,V1(I),V2(I)+CORE
      NREF = NREF+1
      IMIN(NREF) = I
30    FORMAT(I8,2F18.7)
40    CONTINUE
      ELSE
C... READ FROM FILE
      CONTINUE
      END IF
      RETURN
      END
      SUBROUTINE FCVECO (V,ISYM,THRESH)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /TAPES / INP,IOUT
      COMMON /CFULLR/ CORE,THR,THRRES
      COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
     >               ,INTOFF(8),ITYPEA(256)
     >               ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
     >               ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
     >               ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
     >               ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
      COMMON /BIG/ Q(1)
      INTEGER IQ(1)
      EQUIVALENCE (Q(1),IQ(1))
      CHARACTER*130 TEST
      DIMENSION V(*)
      COMMON /CLOCAL/ ICGA(32),ICGB(32)
      N = NCI(ISYM)
      IBASE = ICORR(0)
      INTERA = ICORI(NA*NACT)
      CALL INITI (IQ(INTERA),NA,NACT)
      INTERB = ICORI(NB*NACT)
      CALL INITI (IQ(INTERB),NA,NACT)
      IC = ICORI(NSTRAA)
      ITZ = ICORI(NSTRBB)
      CALL FCMIC (IQ(IC),IQ(ITZ),ISYM)
      NFROZ=0
      DO 60 II=1,NCI(ISYM)
      IF (ABS(V(II)).LT.THRESH) GOTO 60
      CALL FCSTRG (II,ISYM,ICGA,ICGB)
      WRITE (IOUT,'(1X,F20.12,64I3)') V(II)
     1,(ICGA(I)+NFROZ,I=1,NA)
     2,(ICGB(I)+NFROZ,I=1,NB)
60    CONTINUE
      WRITE (IOUT,'(1X,''/EOF'')')
      CALL CORLSR(IBASE)
      RETURN
      END
      SUBROUTINE INITC
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /CFULLR/ CORE,THR,THRRES
      COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
     >               ,INTOFF(8),ITYPEA(256)
     >               ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
     >               ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
     >               ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
     >               ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
      COMMON /TAPES / INP,IOUT
      COMMON /BIG/ Q(1)
      INTEGER IQ(1)
      EQUIVALENCE (Q(1),IQ(1))
      DIMENSION ICG(32)
C...  LENGTHS OF CI EXPANSIONS
C...  NSTRA, NSTRB ARE NUMBERS OF ALPHA,BETA STRINGS OF EACH SYMMETRY
C...  NCI ARE NUMBERS OF DETERMINANTS OF EACH SYMMETRY
C...  NPAIR ARE NUMBERS OF ORBITAL PAIRS OF EACH SYMMETRY
C.... INTOFF OFFSET FOR SYMMETRY BLOCKS OF TWO ELECTRON INTEGRALS
      NSTRAA = 0
      NSTRBB = 0
      MAXAA  = 0
      MAXBB  = 0
      DO 100 ISYM=1,8
      ICG(1) = 0
      NSTRA(ISYM) = -1
80    NSTRA(ISYM) = NSTRA(ISYM)+1
      CALL FCSTRS (NA,ICG,ISYM,ITYPEA)
      IF (ICG(1).NE.0) GOTO 80
      NSTRAA = NSTRAA + NSTRA(ISYM)
      MAXAA = MAX(MAXAA,NSTRA(ISYM))
      ICG(1) = 0
      NSTRB(ISYM) = -1
90    NSTRB(ISYM) = NSTRB(ISYM)+1
      CALL FCSTRS (NB,ICG,ISYM,ITYPEA)
      IF (ICG(1).NE.0) GOTO 90
      NSTRBB = NSTRBB + NSTRB(ISYM)
      MAXBB = MAX(MAXBB,NSTRB(ISYM))
100   CONTINUE
      DO 110 ISYM=1,8
      NCI(ISYM) = 0
      DO 110 ISYMB=1,8
110   NCI(ISYM) = NCI(ISYM) + NSTRA(MULT(ISYMB,ISYM))*NSTRB(ISYMB)
C
C...  ORBITAL PAIR ADDRESSING, ETC..
      DO 120 I=1,256
120   ICF(I) = (I*(I-1))/2
      INTADT = ICORI(NACT*NACT)
      NINT2=0
      MAXPAR = 0
      DO 140 ISYM=1,8
      INTOFF(ISYM) = NINT2
      NPAIR(ISYM)=0
      NPB = 0
      DO 130 I=1,NACT
      DO 130 J=1,I
      IF (ISYM.EQ.MULT(ITYPEA(I),ITYPEA(J))) THEN
      NPAIR(ISYM) = NPAIR(ISYM)+1
      IQ(INTADT+(I-1)*NACT+(J-1)) = NPAIR(ISYM)
      IQ(INTADT+(I-1)+NACT*(J-1)) = NPAIR(ISYM)
      END IF
130   CONTINUE
      MAXPAR = MAX(NPAIR(ISYM),MAXPAR)
140   NINT2 = NINT2 + NPAIR(ISYM)**2
C...  UPPER BOUND ON NUMBER OF SINGLE REPLACEMENTS FROM ANY STRING
      MAXREP = MIN(MAXPAR,MAX(NA*(NACT+1-NA),NB*(NACT+1-NB)))
C
      WRITE (IOUT,141)
     > NACT,(NOCC(I),I=1,8)
     >,NA+NB,(NA-NB)*0.5,NPAIR,NSTRA,NSTRB,NCI
141   FORMAT(
     1' Active orbitals:',              T30,I3,'(',8I3,')'/
     2' Active electrons:',             T30,I3/
     3' Spin quantum number:',          T30,F5.1/
     4' Orbital pairs:',                T30,8I9/
     5' Strings:',                      T30,8I9/T30,8I9/
     6' Determinants:',                 T30,4I12/T30,4I12)
C      WRITE (IOUT,'('' ACTIVE ORBITAL SYMMETRIES:'')')
C      WRITE (IOUT,150) (ITYPEA(I),I=1,NACT)
C150   FORMAT(1X,40I2)
C
C...  LOAD INTEGRALS
      INTAA = ICORR(NINT2)
      INTHAA = ICORR(NACT)
      INTJAA = ICORR(NACT**2)
      INTKAA = ICORR(NACT**2)
      IBASE = ICORR(0)
      IONEA = ICORR(NPAIR(1))
      CALL INPINT (Q(IONEA),Q(INTAA))
C
C...  INTEGRALS FOR DIAGONAL ELEMENTS OF HAMILTONIAN
      DO 10 I=1,NACT
10    Q(INTHAA-1+I) = Q(IONEA-1+IQ(INTADT+(I-1)*(NACT+1)))
      DO 11 I=1,NACT
      DO 11 J=1,NACT
      IJ = IQ(INTADT+(I-1)*NACT+J-1)
      II = IQ(INTADT+(I-1)*NACT+I-1)
      JJ = IQ(INTADT+(J-1)*NACT+J-1)
      ISYMIJ = MULT(ITYPEA(I),ITYPEA(J))
      Q(INTJAA+(I-1)*NACT+J-1) = Q(INTAA+(II-1)*NPAIR(1)+JJ-1)
11    Q(INTKAA+(I-1)*NACT+J-1) =
     >     Q(INTAA+INTOFF(ISYMIJ)+(IJ-1)*(NPAIR(ISYMIJ)+1))
C
C...  MODIFY INTEGRALS SUCH THAT HAMILTONIAN BECOMES
C     SUM(IJKL)  E(IJ) E(KL) (IJ|KL)"
C     (NO 1 ELECTRON TERM REMAINS)
C
C...  (IJ|KL)' = 0.5 * (IJ|KL)
      DO 20 IJKL=1,NINT2
20    Q(INTAA-1+IJKL) = Q(INTAA-1+IJKL)*0.5D0
C
C...  H'(IJ) = H(IJ) - SUM(K) (IK|KJ)'
      IJ = IONEA-1
      DO 40 I=1,NACT
      DO 40 J=1,I
      IF (ITYPEA(I).NE.ITYPEA(J)) GOTO 40
      IJ = IJ+1
      DO 30 K=1,NACT
      ISYMIK = MULT(ITYPEA(I),ITYPEA(K))
      IK =IQ(INTADT+(I-1)*NACT+K-1)-1
      JK =IQ(INTADT+(J-1)*NACT+K-1)-1
30    Q(IJ) = Q(IJ) - Q(INTAA+INTOFF(ISYMIK)+IK*NPAIR(ISYMIK)+JK)
40    CONTINUE
C
C...  (IJ|KL)" = (IJ|KL)'+(1/2*NELEC)*(DEL(IJ)*H'(KL)+DEL(KL)*H'(IJ))
      FAC = 2*(NA+NB)
      FAC = 1D0/FAC
      DO 70 K=1,NACT
      KK = IQ(INTADT+(K-1)*(NACT+1))
      DO 50 IJ=1,NPAIR(1)
50    Q(INTAA+(KK-1)*NPAIR(1)+IJ-1) = Q(INTAA+(KK-1)*NPAIR(1)+IJ-1)
     > + FAC*Q(IONEA-1+IJ)
      DO 60 IJ=1,NPAIR(1)
60    Q(INTAA+(IJ-1)*NPAIR(1)+KK-1) = Q(INTAA+(IJ-1)*NPAIR(1)+KK-1)
     > + FAC*Q(IONEA-1+IJ)
70    CONTINUE
      CALL CORLSR(IBASE)
      WRITE (IOUT,'('' Core energy:'',T30,F20.12)') CORE
      RETURN
      END
      SUBROUTINE FSIGMA (C,S)
CSUBR RETURNS S = H * C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /CFULLR/ CORE,THR,THRRES
      COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
     >               ,INTOFF(8),ITYPEA(256)
     >               ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
     >               ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
     >               ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
     >               ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
      COMMON /TAPES / INP,IOUT
      COMMON /BIG/ Q(1)
      INTEGER IQ(1)
      EQUIVALENCE (Q(1),IQ(1))
      LOGICAL TIMING
      DIMENSION IOFFCI(8),C(*),S(*)
      COMMON /CHAMIL/ ISYMH,NNCI
      TIMING = IPRTIM.GE.0
       NOPMXM = 0
       NOPMX  = 0
       NOPOTH = 0
       NOPCCA = 0
       NOPCCB = 0
       NLOCCA = 0
       NLOCCB = 0
       NOPONE = 0
       T0 = SECOND()
       TMXM = 0
       TCCA = 0
       TCCB = 0
       TONE = 0
      IBASE = ICORR(0)
      NNCI = 0
      DO 10 ISYMB=1,8
      ISYMA = MULT(ISYMH,ISYMB)
      IOFFCI(ISYMB) = NNCI
      IF (IPRHAM.GT.0) WRITE (IOUT,*)'ISYMB,ISYMA,NNCI,NSTRA,NSTRB '
     1,ISYMB,ISYMA,NNCI,NSTRA(ISYMA),NSTRB(ISYMB)
10    NNCI = NNCI + NSTRA(ISYMA)*NSTRB(ISYMB)
      IF (IPRHAM.GT.0) CALL OUTVEC (C,NNCI,'c passed to sigma')
      IC = ICORI(NSTRAA)
      ITZ = ICORI(NSTRBB)
      CALL FCMIC (IQ(IC),IQ(ITZ),ISYMH)
      DO 170 ISYMA=1,8
      DO 170 ISYMB=1,8
      IBASES = ICORR(0)
      ISYMR = MULT(ISYMH,MULT(ISYMA,ISYMB))
      ISYMBE = MULT(ISYMB,ISYMR)
      NSBE = NSTRB(ISYMBE)
      NSA = NSTRA(ISYMA)
      NSB = NSTRB(ISYMB)
      NPAIRR = NPAIR(ISYMR)
      IF (NSA.EQ.0.OR.NSB.EQ.0.OR.NPAIRR.EQ.0) GOTO 170
      MAXREA = MIN(NA*(NACT-NA+1),NPAIRR)
      MAXREB = MAX(1,MIN(NB*(NACT-NB+1),NPAIRR))
      MAXR2A = MAXREA * 2
      MAXR3A = MAXREA * 3
      MAXR2B = MAXREB * 2
      MAXR3B = MAXREB * 3
      IBB = ICORR(0)
      IWA = ICORI(MAXR3A*NSA)
      NWA = ICORI(NSA)
      IWB = ICORI(MAXR3B*NSB)
      NWB = ICORI(NSB)
C...  CONSTRUCT ALL SINGLE REPLACEMENTS
      T1 = SECOND()
      INTER = ICORI(NA*NACT)
      CALL INITI (IQ(INTER),NA,NACT)
      CALL ONELAL (NA,ITYPEA,IQ(INTER),IQ(INTADT),ISYMA,ISYMR
     1 ,IQ(IWA),IQ(NWA),IQ(IC),MAXREA)
      CALL CORLSI(INTER)
      INTER = ICORI(NB*NACT)
      CALL INITI (IQ(INTER),NB,NACT)
      CALL ONELAL (NB,ITYPEA,IQ(INTER),IQ(INTADT),ISYMB,ISYMR
     1 ,IQ(IWB),IQ(NWB),IQ(ITZ),MAXREB)
      TONE = TONE+SECOND()-T1
      CALL CORLSI(INTER)
C...  DETERMINE BLOCKING
      NMAT = 2
      MWORD = ICORRM()/NMAT
      WORDS = MWORD
      NWORD = SQRT(WORDS/NPAIRR)
      NABLK = (MIN(NWORD,NSA)/2)*2+1
      DO 20 NBBLK=(NSB/2)*2+1,1,-2
      IF (NABLK*NBBLK*NPAIRR.LE.MWORD) GOTO 30
20    CONTINUE
30    IDA = ICORR(NABLK*NBBLK*NPAIRR)
      IEA = ICORR(NABLK*NBBLK*NPAIRR)
      IDB = IDA
      IEB = IEA
C...  LOOP OVER BLOCKS
      DO 160 IOFFA=0,NSA-1,NABLK
      NAA = MIN(NABLK,NSA-IOFFA)
      DO 160 IOFFB=0,NSB-1,NBBLK
      NBB = MIN(NBBLK,NSB-IOFFB)
      NAABB = NAA*NBB
      IF (TIMING) NOPOTH = NOPOTH+NAABB*NPAIRR
      CALL FZERO (Q(IDA),NAABB*NPAIRR)
C...  ALPHA REPLACEMENTS
      IF (TIMING) THEN
       T1 = SECOND()
       ICNT = 0
      END IF
      IATABB = IWA+IOFFA*MAXR3A
      IDAA = IDA-1-NAABB
      DO 61 IA=1,NAA
      IATAB = IATABB
      IF (TIMING) ICNT = ICNT+IQ(NWA+IOFFA-1+IA)
      DO 60 IWW=1,IQ(NWA+IOFFA-1+IA)
      ICCC = IOFFB + IQ(IATAB)
      IAAA = IDAA+IQ(IATAB+MAXR2A)*NAABB
      IF (IQ(IATAB+MAXREA).LT.0) THEN
C...  EQUATION (18)
      DO 40 IB=1,NBB
40    Q(IAAA+IB) = - C(ICCC+IB)
      ELSE
C...  EQUATION (18)
      DO 50 IB=1,NBB
50    Q(IAAA+IB) =   C(ICCC+IB)
      END IF
60    IATAB = IATAB+1
      IDAA = IDAA + NBB
61    IATABB = IATABB + MAXR3A
      IF (TIMING) THEN
       TCCA = TCCA + SECOND()-T1
       NOPCCA = NOPCCA + ICNT*NBB
       NLOCCA = NLOCCA + ICNT
      END IF
C...  BETA REPLACEMENTS
      IF (TIMING) THEN
       ICNT = 0
       T1 = SECOND()
      END IF
      IBTABB = IWB+IOFFB*MAXR3B
      ICCCC = IOFFCI(ISYMBE)+(IOFFA-1)*NSBE
      IDBB = IDB-NBB-NAABB
      DO 91 IB=1,NBB
      IBTAB = IBTABB
      IF (TIMING) ICNT = ICNT+IQ(NWB+IOFFB-1+IB)
      DO 90 IWW=1,IQ(NWB-1+IOFFB+IB)
      ICCC = ICCCC+IQ(IBTAB)
      IBBB = IDBB+IQ(IBTAB+MAXR2B)*NAABB
      IF (IQ(IBTAB+MAXREB).LT.0) THEN
C...  EQUATION (19)
      DO 70 IA=1,NAA
70    Q(IBBB+IA*NBB) = Q(IBBB+IA*NBB) - C(ICCC+IA*NSBE)
      ELSE
C...  EQUATION (19)
      DO 80 IA=1,NAA
80    Q(IBBB+IA*NBB) = Q(IBBB+IA*NBB) + C(ICCC+IA*NSBE)
      END IF
90    IBTAB = IBTAB + 1
      IDBB = IDBB + 1
91    IBTABB = IBTABB + MAXR3B
      IF (TIMING) THEN
       TCCB = TCCB + SECOND()-T1
       NOPCCB = NOPCCB + ICNT*NAA
       NLOCCB = NLOCCB + ICNT
      END IF
C
C...  MATRIX MULTIPLICATION WITH INTEGRALS
C...  EQUATION (20)
      IF (TIMING) T1 = SECOND()
      IF (IPRHAM.GT.3)
     1CALL OUTSQR (Q(IDA),NAABB,NAABB,NPAIRR,'DA')
      IF (IPRHAM.GT.3)
     1CALL OUTSQR (Q(INTAA+INTOFF(ISYMR)),NPAIRR,NPAIRR,NPAIRR,'INTAA')
      CALL FCMXM (Q(IDA),Q(INTAA+INTOFF(ISYMR)),Q(IEA),NAABB,NPAIRR,
     >NOPMXM)
      NOPMX = NOPMX+2*NAABB*NPAIRR**2
      IF (IPRHAM.GT.3)
     1CALL OUTSQR (Q(IEA),NAABB,NAABB,NPAIRR,'EA')
      IF (TIMING) TMXM = TMXM+SECOND()-T1
C
C...  ALPHA REPLACEMENTS
      IF (TIMING) THEN
       T1 = SECOND()
       ICNT = 0
      END IF
      IATABB = IWA+IOFFA*MAXR3A
      IEAA = IEA-1-NAABB
      DO 121 IA=1,NAA
      IATAB = IATABB
      IF (TIMING) ICNT = ICNT+IQ(NWA+IOFFA-1+IA)
      DO 120 IWW=1,IQ(NWA+IOFFA-1+IA)
      ICCC = IOFFB + IQ(IATAB)
      IAAA = IEAA+IQ(IATAB+MAXR2A)*NAABB
      IF (IQ(IATAB+MAXREA).LT.0) THEN
C...  EQUATION (21)
      DO 100 IB=1,NBB
100   S(ICCC+IB) = S(ICCC+IB) - Q(IAAA+IB)
      ELSE
C...  EQUATION (21)
      DO 110 IB=1,NBB
110   S(ICCC+IB) = S(ICCC+IB) + Q(IAAA+IB)
      END IF
120   IATAB = IATAB + 1
      IEAA = IEAA + NBB
121   IATABB = IATABB + MAXR3A
      IF (TIMING) THEN
       TCCA = TCCA + SECOND()-T1
       NOPCCA = NOPCCA + ICNT*NBB
       NLOCCA = NLOCCA + ICNT
      END IF
C...  BETA REPLACEMENTS
      IF (TIMING) THEN
       T1 = SECOND()
       ICNT = 0
      END IF
      IBTABB = IWB+IOFFB*MAXR3B
      ICCCC = IOFFCI(ISYMBE)+(IOFFA-1)*NSBE
      IEBB = IEB-NBB-NAABB
      DO 151 IB=1,NBB
      IBTAB = IBTABB
      IF (TIMING) ICNT = ICNT+IQ(NWB+IOFFB-1+IB)
      DO 150 IWW=1,IQ(NWB-1+IOFFB+IB)
      ICCC = ICCCC+IQ(IBTAB)
      IBBB = IEBB+IQ(IBTAB+MAXR2B)*NAABB
      IF (IQ(IBTAB+MAXREB).LT.0) THEN
C...  EQUATION (22)
      DO 130 IA=1,NAA
130   S(ICCC+IA*NSBE) = S(ICCC+IA*NSBE) - Q(IBBB+IA*NBB)
      ELSE
C...  EQUATION (22)
      DO 140 IA=1,NAA
140   S(ICCC+IA*NSBE) = S(ICCC+IA*NSBE) + Q(IBBB+IA*NBB)
      END IF
150   IBTAB = IBTAB + 1
      IEBB = IEBB + 1
151   IBTABB = IBTABB + MAXR3B
      IF (TIMING) THEN
       TCCB = TCCB + SECOND()-T1
       NOPCCB = NOPCCB + ICNT*NAA
       NLOCCB = NLOCCB + ICNT
      END IF
160   CONTINUE
      CALL CORLSR(IBASES)
170   CONTINUE
      CALL CORLSR(IBASE)
      IF (TIMING) THEN
      TTOT = SECOND()-T0
      TOTH = TTOT-TMXM-TCCA-TCCB-TONE
      NOPTOT = NOPOTH+NOPMXM+NOPCCA+NOPCCB+NOPONE
      CALL TPRINT ('String c. coeff',NOPONE,TONE,' ',0)
      CALL TPRINT ('Coupl coeffs(A)',NOPCCA,TCCA,'Av. v length'
     >,NOPCCA/NLOCCA)
      CALL TPRINT ('Coupl coeffs(B)',NOPCCB,TCCB,'Av. v length'
     >,NOPCCB/NLOCCB)
      IAVSP = NINT(1D2-(1D2*NOPMXM)/NOPMX)
      CALL TPRINT ('Matrix multiply',NOPMXM,TMXM,'Av. % sparsity',IAVSP)
      CALL TPRINT ('Other          ',NOPOTH,TOTH,' ',0)
      CALL TPRINT ('Total          ',NOPTOT,TTOT,' ',0)
      END IF
      IF (IPRHAM.GT.0)
     1CALL OUTVEC (S,NNCI,'S RETURNED BY SIGMA')
      RETURN
      END
      SUBROUTINE ONELAL (N,ITYPE,INTER,INTAD,ISYM,ISYMR,IW,NW,ITZ,MAXRE)
CSUBR CONSTRUCT AND STORE ALL ONE ELECTRON REPLACEMENTS
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...  MAXIMUM NUMBERS OF ORBITALS IN ANY SYMMETRY MUST BE LESS THANTHIS!
      PARAMETER (MXORB=64)
      COMMON /BIG/ Q(1)
      INTEGER IQ(1)
      EQUIVALENCE (Q(1),IQ(1))
      COMMON /CFULLR/ CORE,THR,THRRES
      COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
     >               ,INTOFF(8),ITYPEA(256)
     >               ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
     >               ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
     >               ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
     >               ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
      DIMENSION ITYPE(*),INTER(N,*),INTAD(*)
     1  ,IW(*),NW(*),ITZ(*)
C     1  ,IW(MAXRE,3,*),NW(*),ITZ(*)
      DIMENSION ICG(32),ICGI(33),JSTOR1(MXORB),JSTOR2(MXORB)
     >,JSTOR3(MXORB)
      nw(1) = 0
      if (n.le.0) return
      MAXRE2 = MAXRE*2
      MAXRE3 = MAXRE*3
C...  SET COUNTERS TO ZERO AND SET UP LOOK-UP FOR BRAS
      NINSYM = IBINOM(NACT,N)
      ITZI = ICORI(NINSYM)
      NI = 0
      ICGI(1) = 0
      DO 23 II=1,2999999
      IISYM = 0
      CALL FCSTRS (N,ICGI,IISYM,ITYPE)
      IF (ICGI(1).EQ.0) GOTO 24
      IF (IISYM.NE.ISYM) GOTO 23
      NI=NI+1
      IQ(ITZI-1+II) = NI
23    CONTINUE
24    CONTINUE
      DO 25 I=1,NI
25    NW(I) = 0
      IF (NI.EQ.0) GOTO 999
C...  LOOP OVER INTERMEDIATE N-1 ELECTRON STRINGS K
      ICG(1) = 0
      DO 50 ISK=1,2999999
      ISYMK = 0
      CALL FCSTRS (N-1,ICG,ISYMK,ITYPE)
      IF (ICG(1).EQ.0) GOTO 999
      ICG(N) = NACT+1
      ISYMIO = MULT(ISYMK,ISYM)
      ISYMJO = MULT(ISYMR,ISYMIO)
      ISYMJ = MULT(ISYM,ISYMR)
C...  CONSTRUCT CREATIONS TO KETS J (SYMMETRY ISYMJ)
C      DO 777 IREP=1,10
      IREPL = 1
      IADR = 1
      IPARIT = 1
      DO 10 J=1,N-1
10    IADR = IADR + INTER(J+1,ICG(J))
      NJ = 0
      DO 30 JO=1,NACT
      IF(JO.EQ.ICG(IREPL))THEN
      IADR = IADR+INTER(IREPL,JO)-INTER(IREPL+1,JO)
      IREPL = IREPL+1
      IPARIT = -IPARIT
      ELSE IF (ITYPE(JO).EQ.ISYMJO) THEN
C...  GET CONNECTED N ELECTRON DETERMINANT J
C...  STORE INFO
      NJ = NJ+1
      JSTOR1(NJ) = ITZ(IADR+INTER(IREPL,JO))
      JSTOR2(NJ) = IPARIT
      JSTOR3(NJ) = JO
      END IF
30    CONTINUE
C777   CONTINUE
      IF (NJ.EQ.0) GOTO 50
C...  CONSTRUCT CREATIONS TO BRAS I (SYMMETRY ISYM)
      IREPL = 1
      IADR = ITZI
      IPARIT = 1
      DO 110 J=1,N-1
110   IADR = IADR + INTER(J+1,ICG(J))
      DO 130 IO=1,NACT
      IF(IO.EQ.ICG(IREPL))THEN
      IADR = IADR+INTER(IREPL,IO)-INTER(IREPL+1,IO)
      IREPL = IREPL+1
      IPARIT = -IPARIT
      ELSE IF (ITYPE(IO).EQ.ISYMIO) THEN
C...  GET CONNECTED N ELECTRON DETERMINANT I
      I = IQ(IADR+INTER(IREPL,IO))
      IOFFS = (I-1)*MAXRE3+NW(I)
      IOO = (IO-1)*NACT
C...  STORE INFO FOR EACH J
      IF (IPARIT.GT.0) THEN
CDIR$ SHORTLOOP
CDIR$ IVDEP
C$DIR NO_RECURRENCE
*VOPTION INDEP VEC
      DO 45 J=1,NJ
      IW(IOFFS       +J) = JSTOR1(J)
      IW(IOFFS+MAXRE2+J) = JSTOR3(J)+IOO
C      IW(IOFFS+MAXRE2+J) = INTAD(JSTOR3(J)+IOO)
45    IW(IOFFS+MAXRE +J) = JSTOR2(J)
      ELSE
CDIR$ SHORTLOOP
CDIR$ IVDEP
C$DIR NO_RECURRENCE
*VOPTION INDEP VEC
      DO 46 J=1,NJ
      IW(IOFFS       +J) = JSTOR1(J)
      IW(IOFFS+MAXRE2+J) = JSTOR3(J)+IOO
C      IW(IOFFS+MAXRE2+J) = INTAD(JSTOR3(J)+IOO)
46    IW(IOFFS+MAXRE +J) = -JSTOR2(J)
      END IF
      NW(I) = NW(I) + NJ
      END IF
130   CONTINUE
50    CONTINUE
999   CALL CORLSI(ITZI)
C...  GATHER LOOP FOR PACKED INTEGRAL ADDRESSES
      IOFFS = MAXRE2
      DO 147 I=1,NI
      DO 148 J=1,NW(I)
148   IW(IOFFS+J) = INTAD(IW(IOFFS+J))
147   IOFFS = IOFFS + MAXRE3
      RETURN
      END
      SUBROUTINE FCDIAG (S)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /TAPES / INP,IOUT
      COMMON /CHAMIL/ ISYMH,NNCI
      COMMON /CFULLR/ CORE,THR,THRRES
      COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
     >               ,INTOFF(8),ITYPEA(256)
     >               ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
     >               ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
     >               ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
     >               ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
      COMMON /BIG/ Q(1)
      INTEGER IQ(1)
      EQUIVALENCE (Q(1),IQ(1))
      DIMENSION S(*)
      COMMON /CLOCAL/ ICGA(32),ICGB(32)
      IF (IPRTIM.GE.0) T0 = SECOND()
      NOP = 0
      IF (IPRDIA.GT.0) THEN
C      WRITE (IOUT,*) 'COMPUTE DIAGONAL ELEMENTS'
      CALL OUTVEC (Q(INTHAA),NACT,'HAA')
      CALL OUTSQR (Q(INTJAA),NACT,NACT,NACT,'INTJAA')
      CALL OUTSQR (Q(INTKAA),NACT,NACT,NACT,'INTKAA')
      END IF
C      IF (IPRDIA.GT.1) WRITE (IOUT,*) 'MAXAA,MAXBB ',MAXAA,MAXBB
      IZB = ICORR(NACT*MAXBB)
      IZA = ICORR(NACT*MAXAA)
      IZ  = ICORR(NACT*MAXBB)
      IF  = ICORR(MAXBB)
      ICCC=0
      DO 220 ISYMB=1,8
      ISYMA=MULT(ISYMB,ISYMH)
      NSA = NSTRA(ISYMA)
      NSB = NSTRB(ISYMB)
      IF (IPRDIA.GT.1) WRITE (IOUT,*) 'ISYMB,ISYMA,NSA,NSB ',ISYMB
     1,ISYMA,NSA,NSB
      IF (NSA*NSB.LE.0) GOTO 220
      CALL FZERO (Q(IZB),NACT*NSB)
C..   STORE BETA STRING OCCUPANCIES
      ICGB(1) = 0
      DO 20 IB=1,NSB
      CALL FCSTRS (NB,ICGB,ISYMB,ITYPEA)
      DO 10 II=1,NB
10    Q(IZB-1+IB+(ICGB(II)-1)*NSB)=1
20    CONTINUE
      IF (IPRDIA.GT.2) CALL OUTSQR (Q(IZB),NSB,NSB,NACT,'ZB')
C
C..   STORE ALPHA STRING OCCUPANCIES
      ICGA(1) = 0
      CALL FZERO (Q(IZA),NACT*NSA)
      DO 40 IA=1,NSA
      CALL FCSTRS (NA,ICGA,ISYMA,ITYPEA)
      DO 40 II=1,NA
40    Q(IZA-1+IA+(ICGA(II)-1)*NSA) = 1
      IF (IPRDIA.GT.2) CALL OUTSQR (Q(IZA),NSA,NSA,NACT,'ZA')
C
      DO 50 IAB=1,NSA*NSB
50    S(ICCC+IAB) = 0
C
C...   SPECIAL CODE FOR RHF ORBITALS
      DO 210 IA=1,NSA
C...  OCCUPATION NUMBERS OF DETERMINANTS
      DO 100 I=1,NACT
      DO 100 IB=1,NSB
100   Q(IZ-1+IB+(I-1)*NSB) = Q(IZB-1+IB+(I-1)*NSB)
     1   + Q(IZA-1+IA+(I-1)*NSA)
      IF (IPRDIA.GT.3) CALL OUTSQR (Q(IZ),NSB,NSB,NACT,'Z')
C
C...  CONSTANT EXCHANGE CONTRIBUTION & COULOMB
C...  1/2  SUM(IJ) N(I) N(J) ((II|JJ)-0.5(IJ|JI))
C...  EQUATIONS (23), (25)
      IJ = -1
      DO 130 I=1,NACT
      ZZ = Q(INTHAA-1+I)
      IF (IPRDIA.GT.3) WRITE (IOUT,*) '1567;I,HAA(I) ',I,ZZ
      DO 110 IB=1,NSB
110   S(ICCC+IB) = S(ICCC+IB) + ZZ*Q(IZ-1+IB+(I-1)*NSB)
      DO 130 J=1,I
      IJ = (I-1)*NACT-1+J
      ZZ = Q(INTJAA+IJ)-0.5D0*Q(INTKAA+IJ)
      IF (I.EQ.J) ZZ = ZZ*0.5D0
      IF (IPRDIA.GT.3) WRITE (IOUT,*) '1567; I,J,ZZ ',I,J,ZZ
      DO 120 IB=1,NSB
120    S(ICCC+IB) = S(ICCC+IB)
     1 + ZZ*Q(IZ-1+IB+(I-1)*NSB)*Q(IZ-1+IB+(J-1)*NSB)
130   CONTINUE
C
      IF (IPRDIA.GT.3) CALL OUTVEC (S(ICCC+1),NSB,'S AFTER CONSTS')
C...  VV = 4*S**2, SEE EQUATION (26)
      VV = (NA-NB)**2
C..   Z TO HOLD 1 IF SINGLY OCC, 0 OTHERWISE
C...  (CAPITAL N = N(2-N), EQUATION (25))
      CALL FZERO (Q(IF),NSB)
      DO 150 I=1,NACT
      DO 140 IB=1,NSB
140   Q(IZ-1+IB+(I-1)*NSB) = Q(IZ-1+IB+(I-1)*NSB)
     1 * (2.0D0-Q(IZ-1+IB+(I-1)*NSB))
      DO 150 IB=1,NSB
150   Q(IF-1+IB) = Q(IF-1+IB) + Q(IZ-1+IB+(I-1)*NSB)
      DO 160 IB=1,NSB
160   Q(IF-1+IB) = MAX(Q(IF-1+IB),1.1D0)
      DO 170 IB=1,NSB
170   Q(IF-1+IB) = (VV-Q(IF-1+IB))/( Q(IF-1+IB) * (Q(IF-1+IB)-1.0D0) )
C..   F NOW HOLDS F FACTOR, EQUATION (26)
C
      IF (IPRDIA.GT.3) CALL OUTVEC (Q(IF),NSB,'F')
      IJ = -1
      DO 200 I=1,NACT
      IF (I.EQ.1) GOTO 190
      DO 180 J=1,I-1
      IJ = (I-1)*NACT-1+J
      ZZ = -0.5D0*Q(INTKAA+IJ)
      DO 180 IB=1,NSB
180   S(ICCC+IB) = S(ICCC+IB)
     1 + ZZ*Q(IF-1+IB)*Q(IZ-1+IB+(I-1)*NSB)*Q(IZ-1+IB+(J-1)*NSB)
190   IJ = (I-1)*NACT-1+I
      ZZ = -0.25D0*Q(INTKAA+IJ)
      IF (IPRDIA.GT.3) WRITE (IOUT,*) '1875; I,ZZ ',I,ZZ
      DO 200 IB=1,NSB
200   S(ICCC+IB) = S(ICCC+IB) + ZZ*Q(IZ-1+IB+(I-1)*NSB)
210   ICCC = ICCC+NSB
      NOP = NOP+NSA*NSB*(NACT*7+3*NACT*(NACT+1)+5)
220   CONTINUE
      CALL CORLSR(IZB)
      IF (IPRDIA.GE.0) CALL OUTVEC (S,ICCC,'DIAGONAL ELEMENTS')
      IF (IPRTIM.GE.0) CALL TPRINT ('Hamiltonian diag. els.'
     >,NOP,SECOND()-T0,' ',0)
      RETURN
      END
      FUNCTION INICOR(MEMREQ)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER CORMON,CORLSR,CORLSI,ENDCOR
      COMMON /BIG/ Q(1)
      INTEGER IQ(1)
      EQUIVALENCE (Q(1),IQ(1))
      COMMON /TAPES / INP,IOUT
      COMMON /CORCTL/ INTREL,LWORD,LTOP,LMAX,LMIN,MREAL,MXBFF,IBASE,MEM
      SAVE ICOUNT
      DATA ICOUNT/0/
C
C...  EXECUTE ENDCOR CODE IF NOT FIRST CALL TO THIS ROUTINE
      IF (ICOUNT.NE.0) CALL FMAIN (Q(IBASE+1),MEM)
      ICOUNT = ICOUNT+1
C
      MEM = MEMREQ
      LENN = MEM+2
C
C
      CALL GMAINV (Q,IBASE,LENN)
      IBASE = IBASE-1
      IF (LENN.LT.MEM) THEN
       WRITE (IOUT,'('' Request for'',I9,'' words refused'',
     >  '', '',I9,'' available'')') MEM,LENN
       STOP 'Insufficient memory'
      END IF
      WRITE (IOUT,'('' Variable memory set to '',I10,'' words'')') MEM
      LWORD = MEM+IBASE
      LTOP = IBASE
      LMAX = IBASE
      LMIN = IBASE
      INICOR = LWORD-IBASE
      RETURN
C
      ENTRY ICORIM ()
      ICORIM = (LWORD - LTOP) * INTREL
      RETURN
C
      ENTRY ICORRM ()
      ICORRM = LWORD - LTOP
      RETURN
C
      ENTRY ICORI (NWOR)
      NWORD = IABS(NWOR)
      NW = (NWORD+INTREL-1)/INTREL
      MFR = LWORD - LTOP
      IF (NW .GT. MFR) GOTO 20
      ICORI = LTOP*INTREL+1
      LTOP=LTOP+NW
      IF (NWOR.LT.0) LMAX = MAX0(LMAX,LTOP)
      IF (NWOR.GT.0) LMIN = MAX0(LMIN,LTOP)
      RETURN
C
      ENTRY ICORR (NWOR)
      NWORD = IABS(NWOR)
      NW = NWORD
      MFR = LWORD - LTOP
      IF (NW .GT. MFR) GOTO 20
      ICORR = LTOP+1
      LTOP=LTOP+NW
      IF (NWOR.LT.0) LMAX = MAX0(LTOP,LMAX)
      IF (NWOR.GT.0) LMIN = MAX0(LTOP,LMIN)
      RETURN
C
20    CONTINUE
      WRITE (IOUT,*) 'insufficient memory available - require ',NW,
     >               ' have ',MFR
      IF (NW .NE. NWORD) THEN
         WRITE(IOUT,*) 'the request was for ',NWORD,' integer words'
      ELSE
         WRITE(IOUT,*) 'the request was for real words'
      END IF
       STOP 'Insufficient memory'
      END
      SUBROUTINE CORLSI (IAD)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /BIG/ Q(1)
      INTEGER IQ(1)
      EQUIVALENCE (Q(1),IQ(1))
      COMMON /TAPES / INP,IOUT
      COMMON /CORCTL/ INTREL,LWORD,LTOP,LMAX,LMIN,MREAL,MXBFF,IBASE,MEM
      LTOP = (IAD-2+INTREL)/INTREL
      RETURN
C
      ENTRY CORLSR (IAD)
      LTOP = IAD - 1
      RETURN
C
      ENTRY CORMON()
      WRITE (IOUT,40) LMIN-IBASE,LWORD-LMIN
40    FORMAT(/1X,'Minimal core high water =',I8,
     1' real numbers (spare core =',I8,')')
      RETURN
C
      ENTRY ENDCOR()
      CALL FMAIN (Q(IBASE+1),MEM)
      WRITE (IOUT,'('' Variable memory released'')')
      RETURN
      END
      FUNCTION FCISTR (IPAR,INTER,ICG,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION ICG(N),INTER(N,*)
      DIMENSION ICGG(32)
      DO 10 J=1,N
10    ICGG(J) = ICG(J)
      IF (N.LT.2) GOTO 40
      DO 30 JJ=2,N
      JJM = JJ-1
      I = JJM
      DO 20 II=1,JJM
      IF (ICGG(I).LT.ICGG(I+1)) GOTO 20
      J1=ICGG(I)
      ICGG(I)=ICGG(I+1)
      ICGG(I+1)=J1
      IPAR=-IPAR
20    I=I-1
30    CONTINUE
40    CONTINUE
      FCISTR=1
      DO 50 J=1,N
50    FCISTR = FCISTR+INTER(J,ICGG(J))
      RETURN
      END
      SUBROUTINE FCMIC (IC,ITZ,ISYM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /CFULLR/ CORE,THR,THRRES
      COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
     >               ,INTOFF(8),ITYPEA(256)
     >               ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
     >               ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
     >               ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
     >               ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
      DIMENSION ICG(32)
      DIMENSION IC(*),ITZ(*)
      DO 1 I=1,8
      NSTRA(I) = 0
1     NSTRB(I) = 0
      NSTRAT = IBINOM(NACT,NA)
      NSTRBT = IBINOM(NACT,NB)
C...  ALPHA STRING INFORMATION
      ICG(1) = 0
      DO 11 MT=1,NSTRAT
      JMT = 0
      CALL FCSTRS (NA,ICG,JMT,ITYPEA)
      NSTRA(JMT) = NSTRA(JMT)+1
C...  STORE SYMMETRY FOR LATER PROCESSING
11    IC(MT) = -JMT
C...  BETA STRING INFORMATION
      ICG(1) = 0
      DO 10 MT=1,NSTRBT
      JMT = 0
      CALL FCSTRS (NB,ICG,JMT,ITYPEA)
      NSTRB(JMT) = NSTRB(JMT)+1
10    ITZ(MT) = NSTRB(JMT)
C
C...  LOOP OVER SYMMETRIES BUILDING TOTAL OFFSETS
      NNCI = 0
      DO 110 ISYMB=1,8
      ISYMA = MULT(ISYM,ISYMB)
      DO 120 MT=1,NSTRAT
      IF (IC(MT).NE.-ISYMA) GOTO 120
      IC(MT) = NNCI
      NNCI = NNCI + NSTRB(ISYMB)
120   CONTINUE
110   CONTINUE
      RETURN
      END
      SUBROUTINE FCSPAD (C,IREF)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /CFULLR/ CORE,THR,THRRES
      COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
     >               ,INTOFF(8),ITYPEA(256)
     >               ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
     >               ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
     >               ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
     >               ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
      COMMON /BIG/ Q(1)
      INTEGER IQ(1)
      EQUIVALENCE (Q(1),IQ(1))
      DIMENSION C(*)
      COMMON /CHAMIL/ ISYMH,NNCI
      DIMENSION ICGA(32),ICGB(32),IDOCA(32)
     1         ,IDOCB(32),IALP(32),IBET(32)
      DO 10 I=1,NNCI
10    C(I)=0
      IBASE = ICORR(0)
      INTERA = ICORI(NA*NACT)
      INTERB = ICORI(NB*NACT)
      IC     = ICORI(NSTRAA)
      ITZ    = ICORI(NSTRBB)
      CALL INITI (IQ(INTERA),NA,NACT)
      CALL INITI (IQ(INTERB),NB,NACT)
      CALL FCMIC (IQ(IC),IQ(ITZ),ISYMH)
      CALL FCSTRG (IREF,ISYMH,ICGA,ICGB)
      NDOC=0
      NALP=0
      NBET=0
      DO 70 I=1,NACT
      ICASE=1
      DO 20 J=1,NA
      IF (ICGA(J).NE.I) GOTO 20
      ICASE=3
      JA=J
20    CONTINUE
      DO 30 J=1,NB
      IF (ICGB(J).NE.I) GOTO 30
      ICASE=ICASE+1
      JB=J
30    CONTINUE
      GOTO (70,40,50,60),ICASE
40    NBET=NBET+1
      IBET(NBET)=JB
      GOTO 70
50    NALP=NALP+1
      IALP(NALP)=JA
      GOTO 70
60    NDOC=NDOC+1
      IDOCA(NDOC)=JA
      IDOCB(NDOC)=JB
70    CONTINUE
      FAC = 1.0D0/SQRT(DBLE(2**NBET))
      IF (NBET.GT.5) STOP 5
      DO 120 I5=1,2
      DO 110 I4=1,2
      DO 100 I3=1,2
      DO 90 I2=1,2
      DO 80 I1=1,2
      IPAR = 1
      MTA = FCISTR (IPAR,IQ(INTERA),ICGA,NA)
      MTB = FCISTR (IPAR,IQ(INTERB),ICGB,NB)
      C(IQ(IC-1+MTA)+IQ(ITZ-1+MTB)) = DBLE(IPAR)*FAC
      IF (NBET.LT.1) GOTO 130
      J1=ICGB(IBET(1))
      ICGB(IBET(1))=ICGA(IALP(1))
      ICGA(IALP(1))=J1
80    CONTINUE
      IF (NBET.LT.2) GOTO 130
      J1=ICGB(IBET(2))
      ICGB(IBET(2))=ICGA(IALP(2))
      ICGA(IALP(2))=J1
90    CONTINUE
      IF(NBET.LT.3) GOTO 130
      J1=ICGB(IBET(3))
      ICGB(IBET(3))=ICGA(IALP(3))
      ICGA(IALP(3))=J1
100   CONTINUE
      IF (NBET.LT.4) GOTO 130
      J1=ICGB(IBET(4))
      ICGB(IBET(4))=ICGA(IALP(4))
      ICGA(IALP(4))=J1
110   CONTINUE
      IF (NBET.LT.5) GOTO 130
      J1=ICGB(IBET(5))
      ICGB(IBET(5))=ICGA(IALP(5))
      ICGA(IALP(5))=J1
120   CONTINUE
130   CONTINUE
      CALL CORLSR (IBASE)
      RETURN
      END
      SUBROUTINE FCSTRG (II,ISYM,ICGA,ICGB)
CSUBR EXTRACT STRINGS ICGA,ICGB FOR CONFIGURATION II IN SYMMETRY ISYM
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /CFULLR/ CORE,THR,THRRES
      COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
     >               ,INTOFF(8),ITYPEA(256)
     >               ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
     >               ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
     >               ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
     >               ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
      DIMENSION ICGA(*),ICGB(*)
      III = II
      DO 10 ISYMB=1,8
      ISYMA = MULT(ISYMB,ISYM)
      NSB = NSTRB(ISYMB)
      NSA = NSTRA(ISYMA)
      IF (III.LE.NSB*NSA) GOTO 20
10    III = III - NSB*NSA
20    ISA = (III-1)/NSB+1
      ISB = III - (ISA-1)*NSB
      ICGA(1) = 0
      DO 30 JSA=1,ISA
30    CALL FCSTRS (NA,ICGA,ISYMA,ITYPEA)
      ICGB(1) = 0
      DO 40 JSB=1,ISB
40    CALL FCSTRS (NB,ICGB,ISYMB,ITYPEA)
      RETURN
      END
      SUBROUTINE FCSTRS (N,ICG,ISYM,ITYPE)
C... GENERATE A NEW ALPHA OR BETA STRING WITH SYMMETRY ISYM
C...  ICG SHOULD BE PRESERVED BETWEEN CALLS
C... IF (ISYM.EQ.0) ON ENTRY, THEN JUST NEXT STRING GENERATED
C...  ITS SYMMETRY RETURNED IN ISYM
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /CFULLR/ CORE,THR,THRRES
      COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
     >               ,INTOFF(8),ITYPEA(256)
     >               ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
     >               ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
     >               ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
     >               ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
      DIMENSION ICG(*),ITYPE(*)
C
      IF (ICG(1).LE.0) THEN
       DO 10 I=1,N
10     ICG(I)=I
       IF (N.LE.0) THEN
        DO 111 JJ=1,NACT
        IF (ITYPE(JJ).EQ.1) GOTO 112
111     CONTINUE
112     ICG(1) = JJ
       END IF
       GOTO 40
      END IF
C
21    DO 20 IU=N,1,-1
      ICG(IU) = ICG(IU) + 1
      IF (ICG(IU).LE.NACT+IU-N) GOTO 30
20    CONTINUE
      GOTO 60
CDIR$ NEXTSCALAR
*VOPTION NOVEC
30    DO 31 JU=IU+1,N
31    ICG(JU) = ICG(IU)+JU-IU
C
40    JMT = ITYPE(ICG(1))
      DO 50 J=2,N
50    JMT = MULT(JMT,ITYPE(ICG(J)))
      IF (JMT.NE.ISYM.AND.ISYM.NE.0) GOTO 21
      IF (ISYM.EQ.0) ISYM = JMT
      RETURN
60    ICG(1) = 0
      RETURN
      END
      FUNCTION IBINOM (M,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      TOP = M
      NN = MIN0(M-N,N)
      IF (NN.LT.0) GOTO 20
      BOT = NN
      BINOM = 1.0D0
      DO 10 I=1,NN
      BINOM = BINOM * TOP/BOT
      TOP = TOP - 1.0D0
10    BOT = BOT - 1.0D0
      IBINOM = BINOM + 0.1D0
      IF (DABS(BINOM-DBLE(IBINOM)).GT.0.5D0)
     1 STOP 'ARITHMETIC PROBLEM IN COMPUTING BINOMIAL COEFFICIENT'
      RETURN
20    IBINOM = 0
      RETURN
      END
      SUBROUTINE INITI (INTER,NITEM,NBOX)
CSUBR SETS UP PARTIAL WEIGHT ARRAY FOR ADDRESSING BINOMIAL DISTRIBUTIONS
      DIMENSION INTER(NITEM,*)
      N1=NITEM+1
      DO 30 K=1,NITEM
      DO 10 L=1,NBOX
10    INTER(K,L)=0
      DO 20 L=K,NBOX-NITEM+K-1
      INTER(K,L+1) = IBINOM(NBOX-L,NITEM-K)+INTER(K,L)
20    CONTINUE
30    CONTINUE
      DO 40 K=1,NITEM-1
      DO 40 L=K,NBOX-NITEM+K
40    INTER(K,L) = INTER(K,L) - INTER(K+1,L+1)
      DO 50 L=NITEM,NBOX
50    INTER(NITEM,L) = L-NITEM
      RETURN
      END
      SUBROUTINE FZERO (A,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(*)
      DO 1 I=1,N
1     A(I) = 0D0
      RETURN
      END
      SUBROUTINE OUTSQR (Q,IDIM,IA,IB,TITLE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      CHARACTER*(*) TITLE
      DIMENSION Q(IDIM,1)
      CHARACTER*10 FORMAT
      DATA FORMAT /'(/1X,A000)'/
      WRITE (FORMAT(7:9),'(I3)') 55+LEN(TITLE)/2
      WRITE (6,FORMAT) TITLE
      M=1
      N=8
10    IF (IB.LT.M) RETURN
      N=MIN0(N,IB)
      WRITE (6,30) (I,I=M,N)
      WRITE (6,40)
      DO 20 J=1,IA
      WRITE (6,50) J,(Q    (J,I) ,I=M,N)
20    CONTINUE
      M=M+8
      N=N+8
      GOTO 10
30    FORMAT(//3X,8I14)
40    FORMAT(/)
50    FORMAT(7X,I3,8F14.7)
      END
      SUBROUTINE OUTVEC (P,IB,TITLE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /TAPES / INP,IOUT
      CHARACTER*(*) TITLE
      DIMENSION P(1)
      WRITE (IOUT,'(/1X,(A))') TITLE
      M=1
      N=10
10    IF (IB.LT.M) RETURN
      N=MIN0(N,IB)
      WRITE(IOUT,20) M,N,(P    (I) ,I=M,N)
      M=M+10
      N=N+10
      GOTO 10
20    FORMAT(I4,'-',I3,10F12.6)
      END
      SUBROUTINE DIAG2(M,N,D,X)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MAXDIM=200)
      DIMENSION   D(M), X(M,M)
      DIMENSION   E(MAXDIM)
      IF(M.GT.MAXDIM) THEN
        WRITE(6,10) M,MAXDIM
10      FORMAT(' DIMENSION TOO LARGE IN DIAG2:',2I8)
        STOP
      END IF
CSTART VAX UNIX-HP UNIX-CONVEX CRAY IBM UNIVAC
      EPS=2.5E-13
      TOL=2.5D-18
CEND
C
      IF(N.EQ.1) GO TO 400
C
C     HOUSEHOLDER'S REDUCTION
C     SIMULATION OF LOOP DO 150 I=N,2,(-1)
C
      DO 150 NI=2,N
      I=N+2-NI
      L=I-2
      H=0.0
      G=X(I,I-1)
      IF(L) 140,140,20
   20 DO 30 K=1,L
   30 H=H+X(I,K)**2
      S=H+G*G
      IF(S.GE.TOL) GO TO 50
   40 H=0.0
      GO TO 140
   50 IF(H) 140,140,60
   60 L=L+1
      F=G
      G=DSQRT(S)
      IF(F) 75,75,70
   70 G=-G
   75 H=S-F*G
      X(I,I-1)=F-G
      F=0.0
C
      DO 110 J=1,L
      X(J,I)=X(I,J)/H
      S=0.0
      DO 80 K=1,J
   80 S=S+X(J,K)*X(I,K)
      J1=J+1
      IF(J1.GT.L) GO TO 100
      DO 90 K=J1,L
   90 S=S+X(K,J)*X(I,K)
  100 E(J)=S/H
  110 F=F+S*X(J,I)
C
      F=F/(H+H)
C
      DO 120 J=1,L
  120 E(J)=E(J)-F*X(I,J)
C
      DO 130 J=1,L
      F=X(I,J)
      S=E(J)
      DO 130 K=1,J
  130 X(J,K)=X(J,K)-F*E(K)-X(I,K)*S
C
  140 D(I)=H
  150 E(I-1)=G
C
C     ACCUMULATION OF TRANSFORMATION MATRICES
C
  160 D(1)=X(1,1)
      X(1,1)=1.0
      DO 220 I=2,N
      L=I-1
      IF(D(I)) 200,200,170
  170 DO 190 J=1,L
      S=0.0
      DO 180 K=1,L
  180 S=S+X(I,K)*X(K,J)
      DO 190 K=1,L
  190 X(K,J)=X(K,J)-S*X(K,I)
  200 D(I)=X(I,I)
      X(I,I)=1.0
  210 DO 220 J=1,L
      X(I,J)=0.0
  220 X(J,I)=0.0
C
C     DIAGONALIZATION OF THE TRIDIAGONAL MATRIX
C
      B=0.0
      F=0.0
      E(N)=0.0
C
      DO 340 L=1,N
      H=EPS*(DABS(D(L))+DABS(E(L)))
      IF (H.GT.B) B=H
C
C     TEST FOR SPLITTING
C
      DO 240 J=L,N
      IF (DABS(E(J)).LE.B) GOTO 250
  240 CONTINUE
C
C     TEST FOR CONVERGENCE
C
  250 IF(J.EQ.L) GO TO 340
C
C     SHIFT FROM UPPER 2*2 MINOR
C
  260 P=(D(L+1)-D(L))*0.5/E(L)
      R=DSQRT(P*P+1.0)
      IF(P) 270,280,280
  270 P=P-R
      GO TO 290
  280 P=P+R
  290 H=D(L)-E(L)/P
      DO 300 I=L,N
  300 D(I)=D(I)-H
      F=F+H
C
C     QR TRANSFORMATION
C
      P=D(J)
      C=1.0
      S=0.0
C
C     SIMULATION OF LOOP DO 330 I=J-1,L,(-1)
C
      J1=J-1
      DO 330 NI=L,J1
      I=L+J1-NI
      G=C*E(I)
      H=C*P
C
C     PROTECTION AGAINST UNDERFLOW OF EXPONENTS
C
      IF (DABS(P).LT.DABS(E(I))) GOTO 310
      C=E(I)/P
      R=DSQRT(C*C+1.0)
      E(I+1)=S*P*R
      S=C/R
      C=1.0/R
      GO TO 320
  310 C=P/E(I)
      R=DSQRT(C*C+1.0)
      E(I+1)=S*E(I)*R
      S=1.0/R
      C=C/R
  320 P=C*D(I)-S*G
      D(I+1)=H+S*(C*G+S*D(I))
      DO 330 K=1,N
      H=X(K,I+1)
      X(K,I+1)=X(K,I)*S+H*C
  330 X(K,I)=X(K,I)*C-H*S
C
      E(L)=S*P
      D(L)=C*P
      IF (DABS(E(L)).GT.B) GO TO 260
C
C     CONVERGENCE
C
  340 D(L)=D(L)+F
C
C     ORDERING OF EIGENVALUES
C
      NI=N-1
  350 DO 380I=1,NI
      K=I
      P=D(I)
      J1=I+1
      DO 360J=J1,N
      IF(D(J).GE.P) GOTO 360
      K=J
      P=D(J)
  360 CONTINUE
      IF (K.EQ.I) GOTO 380
      D(K) =D(I)
      D(I)=P
      DO 370 J=1,N
      P=X(J,I)
      X(J,I)=X(J,K)
  370 X(J,K)=P
  380 CONTINUE
C
C     FIXING OF SIGN
C
      DO 385 I=1,N
      PM=0
      DO 386 J=1,N
      IF(PM.GT.DABS(X(J,I))) GOTO 386
      PM =DABS(X(J,I))
      K=J
  386 CONTINUE
      IF(X(K,I).GE.0) GOTO 385
      DO 387 J=1,N
  387 X(J,I)=-X(J,I)
  385 CONTINUE
  390 GO TO 410
C
C     SPECIAL TREATMENT OF CASE N = 1
C
  400 D(1)=X(1,1)
      X(1,1)=1.0
  410 RETURN
      END
      SUBROUTINE TPRINT (TITLE,NOP,T,STITLE,NS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /TAPES / INP,IOUT
      COMMON /CFULLR/ CORE,THR,THRRES
      COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
     >               ,INTOFF(8),ITYPEA(256)
     >               ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
     >               ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
     >               ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
     >               ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
      CHARACTER*(*) TITLE,STITLE
      FLOP = 0
      IF (T.NE.0D0) FLOP = 1D-6 * NOP/T
      IF (STITLE.EQ.' ') THEN
      WRITE (IOUT,666) TITLE,NOP,T,FLOP
      ELSE
      WRITE (IOUT,666) TITLE,NOP,T,FLOP,STITLE,NS
      END IF
666   FORMAT(' Timing: ',A,5X,'Operations',I12,5X,'Time',F10.2,
     >5X,'Mflops',F10.3,:,5X,A,I6)
      RETURN
      END