C 19 SEP 05 - MWS - ADD TRUE NUCLEAR CHARGE ARRAY TO INFOA COMMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 30 APR 05 - DGF - SYNCHRONISE SYMBLK COMMON BLOCK
C 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL
C 19 MAY 04 - DGF - CIPA: SUPPRESS SOME OUTPUT
C  9 DEC 03 - MWS - SPHSET: ADDED SAVE FOR SUN COMPILER
C  5 NOV 03 - DGF - ADD UNCONTRACTED BACK TRANSFORMATION, SPHBCKU
C 24 JAN 02 - DGF - PASS NECESSARY STORAGE INTO SPHBCK
C 25 JUN 01 - MWS - COMPLETELY DEACTIVATE BLKTYP OPTION
C 11 JUN 00 - MWS - SALCPR: PRINT SIZE OF VARIATION SPACE
C 21 DEC 99 - DGF - REHABILITATE F AND G TRANSFORMATION MATRICES
C 31 DEC 98 - MAF - SPHMAT: SUBROUTINE ADDED FROM HONDO
C 12 NOV 98 - GDF - REMOVE ALL BIT PACKING USAGE
C 27 OCT 98 - MAF - CIPA,SPHTRF: ALLOW FOR USE OF SPHERICAL HARMONICS;
C                   ACTIVATE SPHMAK; ADD MICHEL'S SPHPRT FOR DEBUGGING
C 14 AUG 97 - MD  - SPHSET,SPHMAK: CHANGES TO NORMALIZATIONS
C 13 JUN 96 - VAG - ADD VARIABLE FOR CI TYPE TO SCFOPT COMMON
C  6 MAR 96 - MWS - CHANGES FOR SHELL SYMMETRY PACKING
C 12 NOV 94 - MWS - REMOVE MOST FTNCHEK WARNINGS
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C  5 NOV 92 - MWS - CHANGE TO 16 BIT PACKING FOR SALC INDICES
C  3 NOV 92 - MWS - CIPA: INCLUDE MICHEL'S COMPLEX GROUP FIX
C 20 MAR 92 - MWS - NEW CODE FROM MICHEL
C
C*MODULE SYMSLC  *DECK CIPA
      SUBROUTINE CIPA(P,PR,PI,PNORM,ISGRP,CA,SAO,TWRK,UWRK,ISB,IPB,
     *                CCHAR,CC,LAB,SHELL,LOC1,LOC2,IWRK1,IWRK2,KLASS,
     *                IRPBAS,IDLAB,DIJPQ,LDIJPQ,GAM,CHARR,
     *                TRACEP,TRACED,TRACEF,TRACEG,DD,SS,VV,TT,PP,ROT,
     *                BUFF,PSPH,DSPH,FSPH,GSPH,
     *                NATM,L1,L3,NSH,MXCIPA,DBG,SOME)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL SKIP,DBG,SOME,ABEL
C
      PARAMETER (MXATM=2000, MXSH=5000, MXGTOT=20000)
C
      DIMENSION P(45,45),PR(45,45),PI(45,45),PNORM(45),
     *          ISGRP(NATM,6),CA(NATM,45),SAO(L3),TWRK(L3),UWRK(L1),
     *          ISB(NSH,14),IPB(NSH,14),CCHAR(5,48),CC(MXCIPA),
     *          LAB(MXCIPA,4),SHELL(NATM,6),LOC1(NAT,5),LOC2(NAT,5),
     *          IWRK1(NATM),IWRK2(NATM),KLASS(NATM),IRPBAS(*),
     *          DIJPQ(1000),LDIJPQ(1000),GAM(48,48),CHARR(48,14),
     *          TRACEP(48),TRACED(48),TRACEF(48),TRACEG(48),
     *          DD(35),ROT(15,15),BUFF(15,15),IDLAB(*),
     *          SS(35,35),VV(35,35),TT(35,35),PP(35,35),
     *          PSPH(3,3),DSPH(6,6),FSPH(10,10),GSPH(15,15)
      DIMENSION ISG(6),KBB(5),KCC(5),SKIP(5),MP(5),MF(5),MIRP(5,14)
C
      COMMON /BASSPH/ QMTTOL,ISPHER
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SYMBLK/ NIRRED,NSALC,NSALC2,NSALC3,NSAFMO
      COMMON /SYMMOL/ GROUP,COMPLX,IGROUP,NAXIS,ILABMO,ABEL
      COMMON /SYMREP/ IRPNAM(14),IPA(14),LAMBDA(14),LAMBD0(14),
     *                IADDR1(14),IADDR2(14),IADDR3(14)
      COMMON /SYMSPD/ PTR(3,144),DTR(6,288),FTR(10,480),GTR(15,720)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
C
      DATA CMPLX    /8HCOMPLEX /
      DATA MXSALC   /45/
      DATA MF       /0,1,4,10,20/
      DATA MP       /1,3,6,10,15/
      DATA TOL      /1.0D-10/
      DATA ZERO,HALF,ONE /0.0D+00,0.5D+00,1.0D+00/
      DATA BLOCKI   /8HBLOCKI  /, BLANK/8HBLANK   /
C
C     ----- CIPA GENERATES COEFFICIENTS OF -----
C        SYMMETRY ADAPTED LINEAR COMBINATIONS
C
C     -----------------------------------------------
C     ----- BEWARE THAT THE -CIPA- COEFFICIENTS -----
C     ----- REFER TO THE SPHERICAL HARMONICS    -----
C     ----- BASIS, NOT TO THE CARTESIAN BASIS.  -----
C     -----------------------------------------------
C
      DO 5 IRP=1,NIRRED
      DO 5 ISH=1,NSHELL
      ISB(ISH,IRP)=0
    5 IPB(ISH,IRP)=0
      DO 10 I=1,5
      DO 10 IAT=1,NAT
      LOC1(IAT,I)=-1
   10 LOC2(IAT,I)=-1
C
C     ----- READ -PTR, DTR, FTR, GTR- MATRICES -----
C
      CALL DAREAD(IDAF,IODA,PTR,  432, 7,0)
      CALL DAREAD(IDAF,IODA,DTR, 1728, 8,0)
      CALL DAREAD(IDAF,IODA,FTR, 4800, 9,0)
      CALL DAREAD(IDAF,IODA,GTR,10800,10,0)
C
      IF(ISPHER.GE.0) THEN
C
C     -----  TRANSFORM TO SPHERICAL HARMONICS BASIS.    -----
C
         CALL SPHSET(DBG,DD,SS,VV,TT,PP,ROT,BUFF,
     *               PSPH,DSPH,FSPH,GSPH)
         CALL SPHTRF(DBG)
         IF(DBG) CALL SPHPRT
C
C     ----- IF REQUESTED, ELIMINATE UNWANTED FUNCTIONS. -----
C
         IF(ISPHER.GT.0) THEN
            CALL SPHCUT(DBG)
            IF(SOME  .AND.  NFG.EQ.0) WRITE(IW,9980)
         END IF
      END IF
C
C     ----- GET TRACE OF -PTR, DTR, FTR, GTR- FOR EACH OPERATION -----
C
      DO 15 IT=1,NT
      N= 3*(IT-1)
      TRACEP(IT)=PTR(1,N+1)+PTR(2,N+2)+PTR(3,N+3)
      N= 6*(IT-1)
      TRACED(IT)=DTR(1,N+1)+DTR(2,N+2)+DTR(3,N+3)+
     1           DTR(4,N+4)+DTR(5,N+5)+DTR(6,N+6)
      N=10*(IT-1)
      TRACEF(IT)=FTR(1,N+1)+FTR(2,N+2)+FTR(3,N+3)+
     1           FTR(4,N+4)+FTR(5,N+5)+FTR(6,N+6)+
     1           FTR(7,N+7)+FTR(8,N+8)+FTR(9,N+9)+
     1           FTR(10,N+10)
      N=15*(IT-1)
      TRACEG(IT)=GTR(1,N+1)+GTR(2,N+2)+GTR(3,N+3)+
     1           GTR(4,N+4)+GTR(5,N+5)+GTR(6,N+6)+
     1           GTR(7,N+7)+GTR(8,N+8)+GTR(9,N+9)+
     1           GTR(10,N+10)+GTR(11,N+11)+GTR(12,N+12)+
     1           GTR(13,N+13)+GTR(14,N+14)+GTR(15,N+15)
   15 CONTINUE
      IF(DBG) THEN
         WRITE(IW,9019) (TRACEP(IT),IT=1,NT)
         WRITE(IW,9018) (TRACED(IT),IT=1,NT)
         WRITE(IW,9017) (TRACEF(IT),IT=1,NT)
         WRITE(IW,9016) (TRACEG(IT),IT=1,NT)
      END IF
C
C     ---- FOR EACH CENTER FIND SUBGROUP WHICH MAPS IT INTO ITSELF ----
C
      DO 66 IAT=1,NAT
C
      K1=1
      K2=1
      K3=1
      K4=1
      K5=1
      K6=1
      J1=0
      J2=0
      J3=0
      J4=0
      J5=0
      J6=0
      DO 64 IT=1,NT
C
C     ----- BUILD FIRST WORD -----
C
      IF(IT-24) 32,32,36
   32 IF(MAPCTR(IAT,IT).NE.IAT) GO TO 34
      J1=J1+K1
   34 K1=K1+K1
      GO TO 64
C
C     ----- BUILD SECOND WORD IF NEEDED -----
C
   36 IF(IT-48) 38,38,42
   38 IF(MAPCTR(IAT,IT).NE.IAT) GO TO 40
      J2=J2+K2
   40 K2=K2+K2
      GO TO 64
C
C     ----- BUILD THIRD WORD IF NEEDED -----
C
   42 IF(IT-72) 44,44,48
   44 IF(MAPCTR(IAT,IT).NE.IAT) GO TO 46
      J3=J3+K3
   46 K3=K3+K3
      GO TO 64
C
C     ----- BUILD FOURTH WORD IF NEEDED -----
C
   48 IF(IT-96) 50,50,54
   50 IF(MAPCTR(IAT,IT).NE.IAT) GO TO 52
      J4=J4+K4
   52 K4=K4+K4
      GO TO 64
C
C     ----- BUILD FIFTH WORD IF NEEDED -----
C
   54 IF(IT-120) 56,56,60
   56 IF(MAPCTR(IAT,IT).NE.IAT) GO TO 58
      J5=J5+K5
   58 K5=K5+K5
      GO TO 64
C
C     ----- BUILD SIXTH WORD IF NEEDED -----
C
   60 CONTINUE
      IF(MAPCTR(IAT,IT).NE.IAT) GO TO 62
      J6=J6+K6
   62 K6=K6+K6
C
   64 CONTINUE
C
C     ----- THESE SIX WORDS DESCRIBE THE ATOM CLASSES -----
C
      ISGRP(IAT,1)=J1
      ISGRP(IAT,2)=J2
      ISGRP(IAT,3)=J3
      ISGRP(IAT,4)=J4
      ISGRP(IAT,5)=J5
      ISGRP(IAT,6)=J6
   66 CONTINUE
C
      IC=0
      ICLASS=0
      DO 70 IAT=1,NAT
   70 KLASS(IAT)=0
      DO 80 IRP=1,NIRRED
   80 IPA(IRP)=0
C
C     ----- PICK CENTER IN -P1- LIST -----
C
      DO 1000 IAT=1,NAT
C
      DO 110 IT=1,NT
      IF(MAPCTR(IAT,IT).GT.IAT) GO TO 1000
  110 CONTINUE
      IF(KLASS(IAT).NE.0) GO TO 1000
      ICLASS=ICLASS+1
      KLASS(IAT)=ICLASS
      DO 120 NQ=1,6
  120 ISG(NQ)=ISGRP(IAT,NQ)
      IF(DBG) WRITE(IW,9985) IAT
C
C     ----- DETERMINE OTHER -P1- CENTERS IN SAME CLASS -----
C
      MIN=IAT+1
      IF(MIN.GT.NAT) GO TO 240
      DO 230 JAT=MIN,NAT
      DO 210 NQ=1,6
      IF(ISGRP(JAT,NQ).NE.ISG(NQ)) GO TO 230
  210 CONTINUE
      DO 220 IT=1,NT
      IF(MAPCTR(JAT,IT).GT.JAT) GO TO 230
  220 CONTINUE
      KLASS(JAT)=ICLASS
  230 CONTINUE
  240 CONTINUE
C
C     ----- CHECK WHICH ANGULAR MOMENTA ARE PRESENT ON THAT CENTER -----
C
      DO 250 KT=1,5
  250 SKIP(KT)=.TRUE.
      DO 260 ISH=1,NSHELL
      JAT=KATOM(ISH)
      IF(KLASS(JAT).NE.ICLASS) GO TO 260
      KT=KTYPE(ISH)
      SKIP(KT)=.FALSE.
      IF(KMIN(ISH).EQ.1) SKIP(1)=.FALSE.
  260 CONTINUE
C
C     ----- FIND CHARACTERS OF REDUCIBLE REPRESENTATION.     -----
C     ----- EQUIVALENT CENTER IATOM WILL CONTRIBUTE ONLY TO  -----
C           CHARACTER OF OPERATIONS WHICH MAP IT INTO ITSELF.
C
      DO 310 IT=1,NT
      DO 310 KT=1,5
  310 CCHAR(KT,IT)=ZERO
C
      DO 350 JAT=1,NAT
      DO 320 IT=1,NT
      IF(MAPCTR(IAT,IT).EQ.JAT) GO TO 330
  320 CONTINUE
      GO TO 350
  330 CONTINUE
      DO 340 IT=1,NT
      IF(MAPCTR(JAT,IT).NE.JAT) GO TO 340
      CCHAR(1,IT)=CCHAR(1,IT)+ONE
      CCHAR(2,IT)=CCHAR(2,IT)+TRACEP(IT)
      CCHAR(3,IT)=CCHAR(3,IT)+TRACED(IT)
      CCHAR(4,IT)=CCHAR(4,IT)+TRACEF(IT)
      CCHAR(5,IT)=CCHAR(5,IT)+TRACEG(IT)
  340 CONTINUE
  350 CONTINUE
C
      IF(DBG) THEN
         DO 360 IT=1,NT
            WRITE(IW,9997) IT,(CCHAR(KT,IT),KT=1,5)
  360    CONTINUE
      END IF
C
C     ----- DECOMPOSE REDUCIBLE REPRESENTATION -----
C           INTO IRREDUCIBLE REPRESENTATIONS.
C
C      -LOC1(ICLASS,ITYP)- = FIRST ADDRESS OF COEFFICIENTS FOR
C                            -ICLASS- AND -S,P,D,F,G- FUNCTIONS .
C      -LOC2(ICLASS,ITYP)- =  LAST ADDRESS OF COEFFICIENTS FOR
C                            -ICLASS- AND -S,P,D,F,G- FUNCTIONS .
C
C     ----- LOOP OVER ANGULAR MOMENTA -----
C
      DO 800 KT=1,5
      IF(SKIP(KT)) GO TO 800
      NF=MF(KT)
      NP=MP(KT)
      LOC1(ICLASS,KT)=IC
C
C     ----- LOOP OVER IRREPS. -----
C
      IRP0=0
      DO 750 IRP=1,NIRRED
      LL=LAMBDA(IRP)
      MIRP(KT,IRP)=0
C
      DUM=ZERO
      DO 570 IT=1,NT
  570 DUM=DUM+CCHAR(KT,IT)*CHARR(IT,IRP)
      DUM=DUM/NT
      IF(DUM.LT.TOL) GO TO 750
      NIRP=NINT(DUM)
      IF(COMPLX.EQ.CMPLX.AND.LL.EQ.2) NIRP=NIRP/2
      MIRP(KT,IRP)=NIRP
      IF(DBG) WRITE(IW,9981) KT,IRP,LL,NIRP,NINT(DUM),DUM
C
C     ----- IRREDUCIBLE REPRESENTATION IRP OCCURS -NIRP- TIMES. -----
C     ----- FORM PROJECTION MATRICES -P-                        -----
C
      DO 580 J=1,NP*LL
      DO 580 I=1,NP*LL
  580 P(I,J)=ZERO
C
      DO 610 IT=1,NT
      IF(MAPCTR(IAT,IT).NE.IAT) GO TO 610
      DO 600 LL1=1,LL
      DO 600 LL2=1,LL
      GM=(LL*GAM(IT,LL2+LL*(LL1-1)+IRP0))/NT
      DO 590 J=1,NP
      DO 590 I=1,NP
      IF(KT.EQ.1) DUM=GM
      IF(KT.EQ.2) DUM=GM*PTR(I,J+NP*(IT-1))
      IF(KT.EQ.3) DUM=GM*DTR(I,J+NP*(IT-1))
      IF(KT.EQ.4) DUM=GM*FTR(I,J+NP*(IT-1))
      IF(KT.EQ.5) DUM=GM*GTR(I,J+NP*(IT-1))
  590 P(I+NP*(LL1-1),J+NP*(LL2-1))=P(I+NP*(LL1-1),J+NP*(LL2-1))+DUM
  600 CONTINUE
  610 CONTINUE
C
C     ----- ORTHOGONALIZE COLUMNS OF -P- -----
C
      IF(COMPLX.EQ.CMPLX.AND.LL.EQ.2) THEN
         DO 620 J=1,NP
         DO 620 I=1,NP
            PR(I   ,J   )= P(I   ,J   )
            PI(I   ,J   )=-P(I+NP,J   )
            PR(I+NP,J   )= ZERO
            PI(I+NP,J   )= ZERO
            PR(I+NP,J+NP)= P(I+NP,J+NP)
            PI(I+NP,J+NP)=-P(I   ,J+NP)
            PR(I   ,J+NP)=ZERO
            PI(I   ,J+NP)=ZERO
  620       CONTINUE
         IF(DBG) CALL SYMPRP(PR,NP*LL,NP*LL,MXSALC)
         IF(DBG) CALL SYMPRP(PI,NP*LL,NP*LL,MXSALC)
                 CALL SYMSMC(PR,PI,PNORM,NP*LL,MXSALC)
         IF(DBG) CALL SYMPRP(PR,NP*LL,NP*LL,MXSALC)
         IF(DBG) CALL SYMPRP(PI,NP*LL,NP*LL,MXSALC)
         DO 630 J=1,NP
         DO 630 I=1,NP
            P(I   ,J   )= HALF * ( PR(I   ,J   ) + PR(I+NP,J+NP))
            P(I+NP,J   )= HALF * (-PI(I   ,J   ) + PI(I+NP,J+NP))
            P(I   ,J+NP)=ZERO
  630       P(I+NP,J+NP)=ZERO
      ELSE
         IF(DBG) CALL SYMPRP(P ,NP*LL,NP*LL,MXSALC)
                 CALL SYMSMD(P ,PNORM,NP*LL,MXSALC)
         IF(DBG) CALL SYMPRV(P ,PNORM,NP*LL,NP*LL,MXSALC)
      END IF
C
C     ----- EXTRACT -NIRP- SYMMETRY ADAPTED VECTOR COMPONENTS -----
C
      IPA0=0
      DO 690 J=1,NP*LL
      IF(PNORM(J).LT.TOL.OR.IPA0.GE.NIRP) GO TO 690
      IPA0    =IPA0    +1
      IPA(IRP)=IPA(IRP)+1
C
C     ----- LOOP OVER PARTNERS IN -IRREP- -----
C
      DO 680 L=1,LL
C
C     ----- PACK AND STORE COMPONENTS -----
C
      DO 670 I=1,NP
      IC=IC+1
      IF(IC.GT.MXCIPA) THEN
         IF(SOME) WRITE(IW,9986) MXCIPA
         CALL ABRT
         STOP
      END IF
C
C     -I1- SERIAL NUMBER OF FUNCTION WITHIN SHELL WHICH HAS THIS COEF.
C     -I2- SERIAL NUMBER OF -REP- WITHIN -IRREP- SET
C     -I3- -IRREP-
C     -I4- PARTNER IN -IRREP-
C     NOTE THAT -I2- CAN EXCEED 8 BITS WITH LARGE BASIS SETS
C
      I1=NF+I
      I2=IPA(IRP)
      I3=IRP
      I4=L
      LAB(IC,1)=I1
      LAB(IC,2)=I2
      LAB(IC,3)=I3
      LAB(IC,4)=I4
C
      DUM=P(I+NP*(L-1),J)
      IF( ABS(DUM).LT.TOL) DUM=ZERO
      CC(IC)=DUM
      IF(DBG) WRITE(IW,9992) IC,(LAB(IC,II),II=1,4),CC(IC)
  670 CONTINUE
  680 CONTINUE
C
  690 CONTINUE
C
C     ----- ACCOUNT FOR OTHER SHELLS WHICH HAVE -----
C           THE EXACT SAME CONTRIBUTIONS.
C
      NIPA= NIRP
      NIPB=-NIPA
      DO 740 ISH=1,NSHELL
      IF(KLASS(KATOM(ISH)).NE.KLASS(IAT)               ) GO TO 740
      IF(KTYPE(ISH).EQ.KT                              ) GO TO 700
      IF(KTYPE(ISH).EQ.2.AND.KMIN(ISH).EQ.1.AND.KT.EQ.1) GO TO 720
      GO TO 740
  700 NIPB=NIPB+NIPA
      IPB(ISH,IRP)=NIPB
      DO 710 IT=1,NT
      JSH=MAPSHL(ISH,IT)
  710 IPB(JSH,IRP)=NIPB
      GO TO 740
  720 NIPB=NIPB+NIPA
      ISB(ISH,IRP)=NIPB
      DO 730 IT=1,NT
      JSH=MAPSHL(ISH,IT)
  730 ISB(JSH,IRP)=NIPB
  740 CONTINUE
C
      IPA(IRP)=IPA(IRP)+NIPB
C
C     ----- END LOOP OVER IRREPS -----
C
  750 IRP0=IRP0+LL*LL
      LOC2(ICLASS,KT)=IC
C
C     ----- END OF LOOP OVER ANGULAR MOMENTUM TYPES -----
C
  800 CONTINUE
C
      DO 810 KT=1,5
      KBB(KT)=LOC1(ICLASS,KT)
  810 KCC(KT)=LOC2(ICLASS,KT)
      JCLASS=ICLASS
C
C     ----- ROTATE COEFFICIENTS ONTO EQUIVALENT CENTERS -----
C
      DO 950 IT=1,NT
C
C     ----- FIND EQUIVALENT CENTER -----
C
      JAT=MAPCTR(IAT,IT)
      IF(DBG) WRITE(IW,9995) JAT
      IF(KLASS(JAT).NE.0) GO TO 950
      ICLASS=ICLASS+1
      KLASS(JAT)=ICLASS
C
C     ----- FIND OTHER CENTERS IN SAME CLASS -----
C
      DO 820 I=1,NAT
      IF(KLASS(I).NE.JCLASS) GO TO 820
      J=MAPCTR(I,IT)
      KLASS(J)=ICLASS
  820 CONTINUE
C
      DO 900 KT=1,5
      IF(SKIP(KT)) GO TO 900
      NP=MP(KT)
      LOC1(ICLASS,KT)=IC
      KB=KBB(KT)
      KC=KCC(KT)
      IF(DBG) WRITE(IW,9994) KT,KB,KC,IC,NP
C
      IRP0=0
      DO 890 IRP=1,NIRRED
      LL=LAMBDA(IRP)
      NIRP=MIRP(KT,IRP)
      IF(NIRP.EQ.0) GO TO 890
      DO 880 J=1,NIRP
C
      DO 870 L=1,LL
      DO 860 I=1,NP
C
      IC=IC+1
      IF(IC.GT.MXCIPA) THEN
         WRITE(IW,9986) MXCIPA
         CALL ABRT
         STOP
      END IF
C
      DUM=ZERO
      D=ZERO
      DO 850 M=1,LL
      DO 850 K=1,NP
      IF(KT.EQ.1) D=ONE
      IF(KT.EQ.2) D=PTR(I,K+NP*(IT-1))
      IF(KT.EQ.3) D=DTR(I,K+NP*(IT-1))
      IF(KT.EQ.4) D=FTR(I,K+NP*(IT-1))
      IF(KT.EQ.5) D=GTR(I,K+NP*(IT-1))
      DUM=DUM+D*CC(K+NP*(M-1)+KB)*GAM(IT,M+LL*(L-1)+IRP0)
  850 CONTINUE
      IF( ABS(DUM).LT.TOL) DUM=ZERO
C
      CC(IC)=DUM
      LAB(IC,1)=LAB(I+NP*(L-1)+KB,1)
      LAB(IC,2)=LAB(I+NP*(L-1)+KB,2)
      LAB(IC,3)=LAB(I+NP*(L-1)+KB,3)
      LAB(IC,4)=LAB(I+NP*(L-1)+KB,4)
      IF(DBG) WRITE(IW,9992) IC,(LAB(IC,II),II=1,4),CC(IC)
  860 CONTINUE
  870 CONTINUE
C
      KB=KB+NP*LL
  880 CONTINUE
C
  890 IRP0=IRP0+LL*LL
C
      LOC2(ICLASS,KT)=IC
C
C     ----- END LOOP OVER ANGULAR MOMENTA -----
C
  900 CONTINUE
C
C     ----- END LOOP OVER EQUIVALENT CENTERS -----
C
  950 CONTINUE
C
C     ----- END OF LOOP OVER UNIQUE CENTERS -----
C
 1000 CONTINUE
      NCC=IC
C
      IF(DBG) THEN
         WRITE(IW,9993)
         DO 1010 I=1,NCC
            WRITE(IW,9992) I,(LAB(I,II),II=1,4),CC(I)
 1010    CONTINUE
         DO 1020 I=1,NSHELL
            WRITE(IW,9991) (ISB(I,J) ,J=1,NIRRED)
            WRITE(IW,9990) (IPB(I,J) ,J=1,NIRRED)
 1020    CONTINUE
      END IF
C
      DO 1040 I=1,NSHELL
      DO 1040 J=1,NIRRED
      ISB(I,J)=ISB(I,J)+1
 1040 IPB(I,J)=IPB(I,J)+1
C
      IF(DBG) THEN
         WRITE(IW,9983) (KLASS(IAT),
     1                  (LOC1(KLASS(IAT),KT),KT=1,5),
     2                  (LOC2(KLASS(IAT),KT),KT=1,5),IAT=1,NAT)
      END IF
C
C     ----- GENERATE INITIAL ADDRESS OF SYMMETRY BLOCKS -----
C          ( THIS MAY GET MODIFIED IN -SALCAO- )
C
      J=0
      K=0
      L=0
      DO 1050 I=1,NIRRED
      LL=LAMBDA(I)
      N =IPA(I)*LL
      IADDR1(I)=J
      IADDR2(I)=K
      IADDR3(I)=L
      J=J+N
      K=K+N*(N+1)/2
      L=L+N*N
 1050 CONTINUE
      NSALC =J
      NSALC2=K
      NSALC3=L
      NSAFMO=NSALC
      IF(DBG) THEN
         WRITE(IW,9989) (IPA(J)   ,J=1,NIRRED)
         WRITE(IW,9988) (IADDR1(J),J=1,NIRRED)
         WRITE(IW,9987) (IADDR2(J),J=1,NIRRED)
         WRITE(IW,9982) (IADDR3(J),J=1,NIRRED)
      END IF
C
C     ----- PRINT -SALC- BASIS FUNCTIONS -----
C
      CALL SALCPR(CA,ISB,IPB,CC,LAB,LOC1,LOC2,IWRK1,KLASS,
     *            NATM,NSH,NCC,MXCIPA,DBG,SOME)
C
C     ----- SET UP -SALC- MATRIX -----
C
      CALL SALCAO(CA,SAO,TWRK,UWRK,ISB,IPB,CC,LAB,LOC1,LOC2,IWRK1,IWRK2,
     *            KLASS,IRPBAS,NATM,L1,L3,NSH,NCC,MXCIPA,DBG)
C
C     ----- CREATE -SALC- CHARGE DISTRIBUTIONS -----
C     INTEGRAL BLOCKING IS A DEFUNCT OPTION IN GAMESS
C
      BLKTYP=BLANK
      IF(BLKTYP.EQ.BLOCKI) CALL DIJGEN(KLASS,IDLAB,
     *                                 CC,LAB,SHELL,LOC1,LOC2,
     *                                 IWRK1,IWRK2,DIJPQ,LDIJPQ,
     *                                 NCC,MXCIPA,NAT,SOME)
C
C     ----- READ -PTR, DTR, FTR, GTR- MATRICES -----
C
      IF(ISPHER.GE.0) THEN
         CALL DAREAD(IDAF,IODA,PTR,  432, 7,0)
         CALL DAREAD(IDAF,IODA,DTR, 1728, 8,0)
         CALL DAREAD(IDAF,IODA,FTR, 4800, 9,0)
         CALL DAREAD(IDAF,IODA,GTR,10800,10,0)
      END IF
C
      RETURN
 9019 FORMAT(/,16H TRACE OF -PTR- ,/,(12F10.5))
 9018 FORMAT(/,16H TRACE OF -DTR- ,/,(12F10.5))
 9017 FORMAT(/,16H TRACE OF -FTR- ,/,(12F10.5))
 9016 FORMAT(/,16H TRACE OF -GTR- ,/,(12F10.5))
 9997 FORMAT(' CHAR. OF REDUCIBLE REP. FOR -IT- AND -S,P,D,F,G- = ',
     1 I5,5F10.5)
 9995 FORMAT(21H EQUIVALENT CENTER = ,I5)
 9994 FORMAT(18H KT,KB,KC,IC,NP = ,5I5)
 9993 FORMAT(/,24H LABELS AND COEFFICIENTS,/,
     1 38H          NF IPA(IRP) IRP  L         C)
 9992 FORMAT(1X,I5,2X,4Z4,E20.12)
 9991 FORMAT(21H ISB(ISHELL,IRREP) = ,/,(14I5))
 9990 FORMAT(21H IPB(ISHELL,IRREP) = ,/,(14I5))
 9989 FORMAT(21H        IPA(IRREP) = ,/,(14I5))
 9988 FORMAT(21H     IADDR1(IRREP) = ,/,(14I5))
 9987 FORMAT(21H     IADDR2(IRREP) = ,/,(14I5))
 9986 FORMAT(23H TOO MANY CIPAS. MAX = ,I5)
 9985 FORMAT(//,14H ----- ATOM = ,I2)
 9983 FORMAT(35H KLASS(IAT) LOC1(KLASS(IAT),KT) AND,
     1                  20H LOC2(KLASS(IAT),KT),
     2  I7,5X,5I5,10X,5I5)
 9982 FORMAT(21H     IADDR3(IRREP) = ,/,(14I5))
 9981 FORMAT(' .... ANG.MOM., IRREP., DEG., NIRP = ',4I5,I5,F10.5)
 9980 FORMAT(/1X,'-- VARIATIONAL SPACE WILL BE RESTRICTED TO PURE',
     *           ' SPHERICAL HARMONICS ONLY --')
      END
C*MODULE SYMSLC  *DECK DIJGEN
      SUBROUTINE DIJGEN(KLASS,IDLAB,CC,LAB,SHELL,
     *                  LOC1,LOC2,ISG,JSG,DIJPQ,LDIJPQ,
     *                  NCC,MXCIPA,NAT,SOME)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL SHELL,SOME,ABEL
C
      PARAMETER (MXATM=2000, MXSH=5000, MXGTOT=20000)
C
      DIMENSION KLASS(NAT),IDLAB(*),
     *          CC(NCC),LAB(MXCIPA,4),SHELL(NAT,6),
     *          LOC1(NAT,5),LOC2(NAT,5),ISG(NAT),JSG(NAT),
     *          DIJPQ(1000),LDIJPQ(1000)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SYMMOL/ GROUP,COMPLX,IGROUP,NAXIS,ILABMO,ABEL
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                TDMY(432),INVT(48),NT
C
      DATA CMPLX  /8HCOMPLEX /
C
      IF(COMPLX.EQ.CMPLX) THEN
         IF(SOME) WRITE(IW,9999)
         CALL ABRT
         STOP
      END IF
C
      DO 5 J=1,6
      DO 5 I=1,NAT
    5 SHELL(I,J)=.TRUE.
      DO 6 II=1,NSHELL
      LIT=KTYPE(II)
      ICEN=KATOM(II)
      IF(LIT.EQ.2.AND.KMIN(II).EQ.1) LIT=6
      SHELL(KLASS(ICEN),LIT)=.FALSE.
    6 CONTINUE
C
C     ----- LOOP OVER -IATOM- FROM -P1- LIST -----
C
      DO 300 IATOM=1,NAT
      DO 10 IT=1,NT
      IF(MAPCTR(IATOM,IT).GT.IATOM) GO TO 300
   10 CONTINUE
      DO 20 I=1,NAT
   20 ISG(I)=0
      DO 30 IT=1,NT
      ICEN=MAPCTR(IATOM,IT)
   30 ISG(ICEN)=KLASS(ICEN)
C
C     ----- LOOP OVER -JATOM- FROM -P1- LIST -----
C
      DO 200 JATOM=1,NAT
      DO 110 IT=1,NT
      IF(MAPCTR(JATOM,IT).GT.JATOM) GO TO 200
  110 CONTINUE
      DO 120 I=1,NAT
  120 JSG(I)=0
      DO 130 IT=1,NT
      JCEN=MAPCTR(JATOM,IT)
  130 JSG(JCEN)=KLASS(JCEN)
C
C     ----- LOOP OVER TYPES OF SHELLS ON -IATOM- -----
C
      DO 180 KT1=1,6
      IF(SHELL(IATOM,KT1)) GO TO 180
C
C     ----- LOOP OVER TYPES OF SHELLS ON -JATOM- -----
C
      DO 170 KT2=1,6
      IF(SHELL(JATOM,KT2)) GO TO 170
      NUMPQA=0
C
C     ----- LOOP OVER ATOMS EQUIVALENT TO -IATOM- -----
C
      DO 160 I=1,NAT
      II=ISG(I)
      IF(II.EQ.0) GO TO 160
C
C     ----- LOOP OVER ATOMS EQUIVALENT TO -JATOM- -----
C
      DO 150 J=1,NAT
      JJ=JSG(J)
      IF(JJ.EQ.0) GO TO 150
      IF(JJ.GT.II) GO TO 150
      IF(JJ.EQ.II.AND.KT2.GT.KT1) GO TO 150
      CALL DIJPQA(CC,LAB,LOC1,LOC2,DIJPQ,LDIJPQ,IDLAB,
     *            NCC,MXCIPA,NAT,II,KT1,JJ,KT2,NUMPQA,SOME)
  150 CONTINUE
  160 CONTINUE
  170 CONTINUE
  180 CONTINUE
  200 CONTINUE
  300 CONTINUE
      RETURN
 9999 FORMAT(' -DIJGEN- CAN NOT HANDLE SYMMETRY GROUPS WITH',
     1       ' COMPLEX REPRESENTATIONS. STOP               ')
      END
C*MODULE SYMSLC  *DECK DIJPQA
      SUBROUTINE DIJPQA(CC,LAB,LOC1,LOC2,DIJPQ,LDIJPQ,IDLAB,
     *                  NCC,MXCIPA,NAT,ISG,LIT1,JSG,LJT1,
     *                  NUMPQA,SOME)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL SPI,SPJ,SOME
C
      DIMENSION CC(NCC),LAB(MXCIPA,4), IDLAB(4,225,2),
     *          LOC1(NAT,5),LOC2(NAT,5),DIJPQ(1000),LDIJPQ(1000)
      DIMENSION KMIN(6),KMAX(6),NF(6),LOCD(225),LOCMIN(225),LOCMAX(225)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /SYMREP/ IRPNAM(14),IPA(14),LAMBDA(14),LAMBD0(14),
     *                IADDR1(14),IADDR2(14),IADDR3(14)
C
      DATA KMIN /1,2, 5,11,21,1/
      DATA KMAX /1,4,10,20,35,4/
      DATA NF   /1,3, 6,10,15,4/
      DATA TOL  /1.0D-10/
      DATA MAXPQA,NDMAX /225,1000/
C
C           IDLAB IS USED BEFORE IT IS INITIALIZED
C           HONDO-8 HAS IDLAB IN STATIC COMMON /IJPQA/,
C           WHERE IT IS NOT INITIALIZED BY SYM.F
C
      CIPA2=0.0D+00
      CIPA3=0.0D+00
C
      LIT=LIT1
      LJT=LJT1
      I1S=6*(ISG-1)+LIT
      I2S=6*(JSG-1)+LJT
      IF(I1S.GE.I2S) GO TO 20
      I3=I1S
      I1S=I2S
      I2S=I3
   20 CONTINUE
      MINI=KMIN(LIT)
      MINJ=KMIN(LJT)
      MAXI=KMAX(LIT)
      MAXJ=KMAX(LJT)
      MINPI=MINI
      MINPJ=MINJ
      IF(LIT.EQ.6) LIT=2
      IF(LJT.EQ.6) LJT=2
      NF1=NF(LIT)
      NF2=NF(LJT)
C
C     ----- LOOP OVER -I- FUNCTIONS -----
C
      KOUNT=0
      IJ=0
      DO 1000 I=MINI,MAXI
      SPI=I.EQ.1.AND.LIT.EQ.2
      LT=LIT
      IF(SPI) LT=1
      LOCIA=LOC1(ISG,LT)
      LOCIB=LOC2(ISG,LT)
      NFI=NF1
      IF(SPI) NFI=1
      ID=MINI
      IF(LIT.EQ.2.AND.MINPI.EQ.1.AND.I.GT.1) ID=2
C
C     ----- LOOP OVER -J- FUNCTIONS -----
C
      DO 900 J=MINJ,MAXJ
      SPJ=J.EQ.1.AND.LJT.EQ.2
      LT=LJT
      IF(SPJ) LT=1
      LOCJA=LOC1(JSG,LT)
      LOCJB=LOC2(JSG,LT)
      NFJ=NF2
      IF(SPJ) NFJ=1
      JD=MINJ
      IF(LJT.EQ.2.AND.MINPJ.EQ.1.AND.J.GT.1) JD=2
      IJ=IJ+1
      LOCMIN(IJ)=KOUNT+1
      I1=LOCIA+I-ID+1
  100 IF(I1.GT.LOCIB) GO TO 800
C
C     -I4-     PARTNER IN -IRREP-
C     -I3-     -IRREP-
C     -I2-     SERIAL NUMBER AMONG -IRREP- SET
C     -I1-     COEFFICIENT LOCATION
C
      ILAB=I1
      I1=LAB(ILAB,1)
      I2=LAB(ILAB,2)
      I3=LAB(ILAB,3)
      I4=LAB(ILAB,4)
C
      IP=I2
      IA=I3
      IALPHA=I4
      LL=LAMBDA(IA)
      CIPA1=CC(I1)
      IF(LL.EQ.1) GO TO 120
      I1=I1+NFI
      CIPA2=CC(I1)
      IF(LL.EQ.2) GO TO 120
      I1=I1+NFI
      CIPA3=CC(I1)
  120 CONTINUE
      J1=LOCJA+J-JD+1
  300 IF(J1.GT.LOCJB) GO TO 750
C
      I1=LAB(J1,1)
      I2=LAB(J1,2)
      I3=LAB(J1,3)
      I4=LAB(J1,4)
C
      JA=I3
      IF(JA.NE.IA) GO TO 700
      IF(I4.NE.IALPHA) GO TO 650
      IQ=I2
      DUM=CIPA1*CC(J1)
      IF(LL.EQ.1) GO TO 340
      J1=J1+NFJ
      DUM=DUM+CIPA2*CC(J1)
      IF(LL.EQ.2) GO TO 320
      J1=J1+NFJ
      DUM=DUM+CIPA3*CC(J1)
  320 DUM=DUM/LL
  340 CONTINUE
      IF( ABS(DUM).LT.TOL) GO TO 650
      KOUNT=KOUNT+1
C
      IF(KOUNT.GT.NDMAX) THEN
         IF(SOME) WRITE(IW,9004) NDMAX
         CALL ABRT
         STOP
      END IF
C
      IF(NUMPQA.EQ.0) GO TO 380
      DO 360 M=1,NUMPQA
         N3=IDLAB(1,M,1)
         N4=IDLAB(1,M,2)
         N5=IDLAB(2,M,2)
         IF(N3.EQ.IP.AND.N4.EQ.IQ.AND.N5.EQ.IA) GO TO 400
  360 CONTINUE
  380 NUMPQA=NUMPQA+1
C
      IF(NUMPQA.GT.MAXPQA) THEN
         IF(SOME) WRITE(IW,9006)
         CALL ABRT
         STOP
      END IF
C
      IDUM=0
      IB=0
      IC=0
C
C     ----- THESE LABELS NEEDED FOR -CI- INTEGRALS -----
C     -IB- AND -IC-
C
      IDLAB(1,NUMPQA,1)=IDUM
      IDLAB(2,NUMPQA,1)=IC
      IDLAB(3,NUMPQA,1)=IB
      IDLAB(4,NUMPQA,1)=IA
C
      IDLAB(1,NUMPQA,2)=IQ
      IDLAB(2,NUMPQA,2)=IP
      IDLAB(3,NUMPQA,2)=J
      IDLAB(4,NUMPQA,2)=I
C
      LOCD(NUMPQA)=NUMPQA
      M=NUMPQA
  400 LDIJPQ(KOUNT)=M
      DIJPQ(KOUNT)=DUM
      IF(SOME) WRITE(IW,9001) ISG,LIT,JSG,LJT,I,J,IP,IQ,IA,DUM,M
  650 J1=J1+NFJ
      GO TO 300
  700 J1=J1+NFJ*LAMBDA(JA)
      GO TO 300
  750 I1=I1+NFI
      GO TO 100
  800 CONTINUE
      LOCMAX(IJ)=KOUNT
  900 CONTINUE
 1000 CONTINUE
      IF(SOME) THEN
         WRITE(IW,9002) (LOCD(I),I=1,NUMPQA)
         WRITE(IW,9002) (LOCMIN(I),I=1,IJ)
         WRITE(IW,9002) (LOCMAX(I),I=1,IJ)
         WRITE(IW,9003) ((IDLAB(J,I,1),J=1,4),I=1,NUMPQA)
         WRITE(IW,9003) ((IDLAB(J,I,2),J=1,4),I=1,NUMPQA)
      END IF
C
      IF(SOME) WRITE(IW,9000) I1S,I2S,KOUNT
      RETURN
C
 9000 FORMAT(4H I1=,I5,6H   I2=,I5,9H   KOUNT=,I5/1X,90(1H.)/)
 9001 FORMAT(1X,9I4,E25.14,I10)
 9002 FORMAT(36I3)
 9003 FORMAT(6(1X,4Z4))
 9004 FORMAT(47H TOO MANY CHARGE DENSITY COEFFICIENTS. NDMAX = ,I10)
 9006 FORMAT(' -NUMPQA- IS TOO LARGE. STOP')
      END
C*MODULE SYMSLC  *DECK SALCAO
      SUBROUTINE SALCAO(CA,SAO,TWRK,UWRK,ISB,IPB,CC,LAB,LOC1,LOC2,
     *                  MAPC,MAPS,KLASS,IRPBAS,NATM,L1,L3,NSH,NCC,
     *                  MXCIPA,DBG)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL DBG,ABEL
C
      PARAMETER (MXATM=2000, MXSH=5000, MXGTOT=20000)
C
      DIMENSION CA(NATM,45),SAO(L3),TWRK(L3),UWRK(L1),
     *          ISB(NSH,14),IPB(NSH,14),
     *          CC(MXCIPA),LAB(MXCIPA,4),LOC1(NATM,5),LOC2(NATM,5),
     *          MAPC(NATM),MAPS(NATM),KLASS(NATM),IRPBAS(L1,3)
      DIMENSION MP(5)
C
      COMMON /BASSPH/ QMTTOL,ISPHER
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SYMBLK/ NIRRED,NSALC,NSALC2,NSALC3,NSAFMO
      COMMON /SYMMOL/ GROUP,COMPLX,IGROUP,NAXIS,ILABMO,ABEL
      COMMON /SYMREP/ IRPNAM(14),IPP(14),LAMBDA(14),LAMBD0(14),
     *                IADDR1(14),IADDR2(14),IADDR3(14)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
C
      DATA CMPLX /8HCOMPLEX /
      DATA ZERO  /0.0D+00/
      DATA MP    /1,3,6,10,15/
C
C     ----- ROUTINE SETS UP MATRIX OF -SALC- ORBITALS -----
C
      IF(DBG) WRITE(IW,9999)
      DO 10 I=1,NUM*NUM
   10 SAO(I)=ZERO
      IPA0=0
C
C     ----- LOOP OVER -P1- CENTERS -----
C
      DO 210 IAT=1,NAT
      DO 20 IT=1,NT
      IF(MAPCTR(IAT,IT).GT.IAT) GO TO 210
   20 CONTINUE
C
      DO 30 JAT=1,NAT
   30 MAPC(JAT)=0
      DO 40 IT=1,NT
      JAT=MAPCTR(IAT,IT)
   40 MAPC(JAT)=1
      ICLASS=KLASS(IAT)
C
C     ----- LOOP OVER SHELL TYPES ON -P1- CENTER -----
C
      DO 200 KT=1,5
      NP=MP(KT)
C
      ISH=0
   50 ISH=ISH+1
      IF(ISH.GT.NSHELL) GO TO 200
      IF(KATOM(ISH).NE.IAT) GO TO 50
      IF(KTYPE(ISH).EQ.KT ) GO TO 60
      IF(KT.EQ.1.AND.KTYPE(ISH).EQ.2.AND.KMIN(ISH).EQ.1) GO TO 60
      GO TO 50
C
   60 CONTINUE
      DO 70 IT=1,NT
      JSH=MAPSHL(ISH,IT)
      JAT=KATOM(JSH)
   70 MAPS(JAT)=JSH
C
      KA=LOC1(ICLASS,KT)
      KAA=KA
      KB=LOC2(ICLASS,KT)
  100 CONTINUE
C
      I1=LAB(KA+1,1)
      I2=LAB(KA+1,2)
      I3=LAB(KA+1,3)
      I4=LAB(KA+1,4)
C
      IRP=I3
      IF(            KTYPE(ISH).EQ.KT                   )
     1 IPA0=IPB(ISH,IRP)
      IF(KT.EQ.1.AND.KTYPE(ISH).EQ. 2.AND.KMIN(ISH).EQ.1)
     1 IPA0=ISB(ISH,IRP)
      IPA=I2+IPA0-1
      LL=LAMBDA(IRP)
C
C     ----- GET -SALC- FOR SHELL ON -P1- CENTER -----
C
      MAX=NP*LL
      DO 110 I=1,MAX
  110 CA(IAT,I)=CC(KA+I)
      LOCI=KLOC(ISH)
      IF(KT.EQ.2.AND.KTYPE(ISH).EQ.2.AND.KMIN(ISH).EQ.1) LOCI=LOCI+1
      IC=0
      DO 115 L=1,LL
      DO 115 I=1,NP
      IC=IC+1
      IROW=(LOCI+I-1)
      LCOL=(IADDR1(IRP)+IPP(IRP)*(L-1)+IPA-1)
      IL=IROW+NUM*LCOL
      DUM=CC(KA+IC)
      IF(DBG) WRITE(IW,9991) IROW,LCOL,DUM
  115 SAO(IL)=DUM
C
C     ----- NOW GET -SALC- FOR EQUIVALENT SHELLS AND CENTERS -----
C
      DO 130 JAT=1,NAT
      IF(MAPC(JAT).EQ.0) GO TO 130
      IF(JAT.EQ.IAT) GO TO 130
      JCLASS=KLASS(JAT)
      KC=LOC1(JCLASS,KT)+KA-KAA
      MAX=NP*LL
      DO 120 I=1,MAX
  120 CA(JAT,I)=CC(KC+I)
      JSH=MAPS(JAT)
      LOCJ=KLOC(JSH)
      IF(KT.EQ.2.AND.KTYPE(JSH).EQ.2.AND.KMIN(JSH).EQ.1) LOCJ=LOCJ+1
      JC=0
      DO 125 L=1,LL
      DO 125 J=1,NP
      JC=JC+1
      JROW=(LOCJ+J-1)
      LCOL=(IADDR1(IRP)+IPP(IRP)*(L-1)+IPA-1)
      JL=JROW+NUM*LCOL
      DUM=CC(KC+JC)
      IF(DBG) WRITE(IW,9991) JROW,LCOL,DUM
  125 SAO(JL)=DUM
  130 CONTINUE
C
C     ----- CREATE SYMMETRY LABEL OF THE -SALC- BASIS FUNCTIONS -----
C
      DO 135 L=1,LL
         IRDEX=IADDR1(IRP)+IPP(IRP)*(L-1)+IPA
         IRPBAS(IRDEX,1)=IRP
         IRPBAS(IRDEX,2)=LL
         IRPBAS(IRDEX,3)=L
  135 CONTINUE
C
C     ----- NOW PRINT -SALC- COMBINATIONS -----
C
      IF(DBG) THEN
         IF(ISPHER.GE.0) THEN
            CALL SLCSPH(IPA,IRPNAM(IRP),MAPC,NP,LL,KT,CA,NATM)
            CALL SPHXYZ(IPA,IRPNAM(IRP),MAPC,NP,LL,KT,CA,NATM)
         ELSE
            CALL SLCXYZ(IPA,IRPNAM(IRP),MAPC,NP,LL,KT,CA,NATM)
         END IF
      END IF
C
      KA=KA+MAX
      IF(KA.LT.KB) GO TO 100
C
      GO TO 50
C
C     ----- END LOOP OVER SHELL TYPES ON -P1- CENTER -----
C     ----- END LOOP OVER -P1- CENTERS -----
C
  200 CONTINUE
  210 CONTINUE
C
      IF(DBG) THEN
         WRITE(IW,9996) (IRPNAM(I),IPP(I),I=1,NIRRED)
         DO 220 I=1,NSALC
            IRP = IRPBAS(I,1)
            LL  = IRPBAS(I,2)
            L   = IRPBAS(I,3)
            WRITE(IW,9995) I,IRP,IRPNAM(IRP),LL,L
  220    CONTINUE
      END IF
C
C     ----- AT THIS POINT THE -SALC- COEFFICIENTS MAY BE -----
C            EXPRESSED IN THE SPHERICAL HARMONICS BASIS.
C            BACKTRANSFORM TO CARTESIAN GAUSSIAN BASIS.
C
      IF(ISPHER.GE.0) CALL SPHBCK(SAO(1),TWRK,UWRK,NUM,DBG)
C     IF(NFG.NE.0) CALL FMOSALC(SAO,NSALC)
C
      IF(DBG) THEN
         WRITE(IW,*) 'SALC MATRIX IS'
         CALL SYMPRP(SAO,NUM,NUM,NUM)
      END IF
C
C     ----- SAVE -SALC- MATRIX ON DAF -----
C
      NDUM=NUM*NUM
      CALL DAWRIT(IDAF,IODA,SAO,NDUM,44,0)
C
C     ----- SAVE TRUE -IRREP- DEGENERACY -----
C           -LAMBDA- MAY BE CHANGED ( SEE BELOW )
C
      DO 225 IRP=1,NIRRED
  225 LAMBD0(IRP)=LAMBDA(IRP)
C
C     ----- SPECIAL CODE FOR COMPLEX REPRESENTATIONS:     -----
C           COMPLEX -E- (2D) OF DIMENSION - N- IS TREATED
C           AS REAL -E- (1D) OF DIMENSION -2N- .
C
      IF(COMPLX.NE.CMPLX) GO TO 300
C
      IF(DBG) WRITE(IW,9994)
      IF(DBG) WRITE(IW,9993)
      DO 230 IC=1,NCC
C
      I1=LAB(IC,1)
      I2=LAB(IC,2)
      I3=LAB(IC,3)
      I4=LAB(IC,4)
C
      IRP=I3
      IF(LAMBDA(IRP).NE.2.OR.I4.NE.2) GO TO 230
      I2=I2+IPP(IRP)
      I4=1
C
      LAB(IC,1)=I1
      LAB(IC,2)=I2
      LAB(IC,3)=I3
      LAB(IC,4)=I4
  230 IF(DBG) WRITE(IW,9992) IC,(LAB(IC,J),J=1,4),CC(IC)
C
      DO 240 IRP=1,NIRRED
         IF(LAMBDA(IRP).NE.2) GO TO 240
         LAMBDA(IRP)=1
         IPP(IRP)=IPP(IRP)*2
  240 CONTINUE
C
      IF(DBG) WRITE(IW,9990)
      IF(DBG) WRITE(IW,9996) (IRPNAM(I),IPP(I),I=1,NIRRED)
C
C     ----- DETERMINE INITIAL ADDRESS OF SYMMETRY BLOCKS -----
C
  300 CONTINUE
      J=0
      K=0
      L=0
      DO 400 I=1,NIRRED
         LL=LAMBDA(I)
         N =IPP(I)*LL
         IADDR1(I)=J
         IADDR2(I)=K
         IADDR3(I)=L
         J=J+N
         K=K+N*(N+1)/2
         L=L+N*N
  400 CONTINUE
C
      NSALC =J
      NSALC2=K
      NSALC3=L
      IF(DBG) THEN
         WRITE(IW,9989) (IPP(J)   ,J=1,NIRRED)
         WRITE(IW,9988) (IADDR1(J),J=1,NIRRED)
         WRITE(IW,9987) (IADDR2(J),J=1,NIRRED)
         WRITE(IW,9986) (IADDR3(J),J=1,NIRRED)
      END IF
      RETURN
C
 9999 FORMAT(/,' -- -SALC- MATRIX -- ',/)
 9996 FORMAT(/,35H DIMENSION OF SYMMETRY SUBSPACES = ,
     1       /,8(1X,A4,I4),/,6(1X,A4,I4))
 9995 FORMAT(' BF = ',I5,' LABEL = ',Z16,
     1       ' IRREP, DEG, PARTNER = ',1X,A4,I3,I3)
 9994 FORMAT(' ADJUSTMENT FOR COMPLEX REPRESENTATIONS ')
 9993 FORMAT(/1X,' PACKED LABELS AND COEFFICIENTS'/
     *        1X,'    NF       IPA(IRP) IRP L         C')
 9992 FORMAT(1X,I5,2X,4Z4,F20.10)
 9991 FORMAT(' SALC ELEMENT , I,J,C = ',2I6,F15.8)
 9990 FORMAT(' FOR GROUPS WITH COMPLEX REPRESENTATIONS',
     1       ' THE -E- REPRESENTATIONS ARE NOW TREATED',/,
     2       ' AS -1D- REPRESENTATIONS WITH DOUBLED DIMENSION')
 9989 FORMAT(21H        IPA(IRREP) = ,/,(14I5))
 9988 FORMAT(21H     IADDR1(IRREP) = ,/,(14I5))
 9987 FORMAT(21H     IADDR2(IRREP) = ,/,(14I5))
 9986 FORMAT(21H     IADDR3(IRREP) = ,/,(14I5))
      END
C*MODULE SYMSLC  *DECK SALCPR
      SUBROUTINE SALCPR(CA,ISB,IPB,CC,LAB,LOC1,LOC2,MAPC,KLASS,
     *                  NATM,NSH,NCC,MXCIPA,DBG,SOME)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL DBG,SOME
C
      PARAMETER (MXATM=2000, MXSH=5000, MXGTOT=20000)
C
      DIMENSION CA(NATM,45),ISB(NSH,14),IPB(NSH,14),CC(NCC),
     *          LAB(MXCIPA,4),LOC1(NATM,5),LOC2(NATM,5),MAPC(NATM),
     *          KLASS(NATM)
      DIMENSION MP(5)
C
      COMMON /BASSPH/ QMTTOL,ISPHER
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SYMBLK/ NIRRED,NSALC,NSALC2,NSALC3,NSAFMO
      COMMON /SYMREP/ IRPNAM(14),IPP(14),LAMBDA(14),LAMBD0(14),
     *                IADDR1(14),IADDR2(14),IADDR3(14)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
C
      DATA MP  /1,3,6,10,15/
C
C     ----- ROUTINE PRINTS SYMMETRY ADAPTED LINEAR COMBINATIONS -----
C
      IF(DBG) WRITE(IW,9999)
      IPA0=0
C
C     ----- LOOP OVER -P1- CENTERS -----
C
      DO 210 IAT=1,NAT
      DO 10 IT=1,NT
      IF(MAPCTR(IAT,IT).GT.IAT) GO TO 210
   10 CONTINUE
C
      DO 20 JAT=1,NAT
   20 MAPC(JAT)=0
      DO 30 IT=1,NT
      JAT=MAPCTR(IAT,IT)
   30 MAPC(JAT)=1
      ICLASS=KLASS(IAT)
C
C     ----- LOOP OVER SHELL TYPES ON -P1- CENTER -----
C
      DO 200 KT=1,5
      NP=MP(KT)
C
      ISH=0
   40 ISH=ISH+1
      IF(ISH.GT.NSHELL) GO TO 200
      IF(KATOM(ISH).NE.IAT) GO TO 40
      IF(KTYPE(ISH).EQ.KT ) GO TO 50
      IF(KT.EQ.1.AND.KTYPE(ISH).EQ.2.AND.KMIN(ISH).EQ.1) GO TO 50
      GO TO 40
C
   50 CONTINUE
C
      KA=LOC1(ICLASS,KT)
      KAA=KA
      KB=LOC2(ICLASS,KT)
  100 CONTINUE
C
      I2=LAB(KA+1,2)
      I3=LAB(KA+1,3)
C
      IRP=I3
      IF(            KTYPE(ISH).EQ.KT                  )
     1 IPA0=IPB(ISH,IRP)
      IF(KT.EQ.1.AND.KTYPE(ISH).EQ.2.AND.KMIN(ISH).EQ.1)
     1 IPA0=ISB(ISH,IRP)
      IPA=I2+IPA0-1
      LL=LAMBDA(IRP)
C
C     ----- GET -SALC- FOR SHELL ON -P1- CENTER -----
C
      MAX=NP*LL
      DO 110 I=1,MAX
  110 CA(IAT,I)=CC(KA+I)
C
C     ----- NOW GET -SALC- FOR EQUIVALENT SHELLS AND CENTERS -----
C
      DO 130 JAT=1,NAT
      IF(MAPC(JAT).EQ.0) GO TO 130
      IF(JAT.EQ.IAT) GO TO 130
      JCLASS=KLASS(JAT)
      KC=LOC1(JCLASS,KT)+KA-KAA
      MAX=NP*LL
      DO 120 I=1,MAX
  120 CA(JAT,I)=CC(KC+I)
  130 CONTINUE
C
C     ----- NOW PRINT -SALC- COMBINATION -----
C
      IF(DBG) THEN
         IF(ISPHER.GE.0) THEN
            CALL SLCSPH(IPA,IRPNAM(IRP),MAPC,NP,LL,KT,CA,NATM)
            CALL SPHXYZ(IPA,IRPNAM(IRP),MAPC,NP,LL,KT,CA,NATM)
         ELSE
            CALL SLCXYZ(IPA,IRPNAM(IRP),MAPC,NP,LL,KT,CA,NATM)
         END IF
      END IF
C
      KA=KA+MAX
      IF(KA.LT.KB) GO TO 100
C
C     ----- END LOOP OVER SHELL TYPES ON -P1- CENTER -----
C
  200 CONTINUE
C
C     ----- END LOOP OVER -P1- CENTERS -----
C
  210 CONTINUE
      IF(SOME) THEN
         IF(ISPHER.EQ.1   .AND.  NFG.EQ.0) WRITE(IW,9998) NSALC
         WRITE(IW,9997)
         WRITE(IW,9996) (IRPNAM(I),IPP(I),I=1,NIRRED)
      END IF
      RETURN
C
 9999 FORMAT(/43H -- SYMMETRY ADAPTED LINEAR COMBINATIONS --/)
 9998 FORMAT(1X,'AFTER EXCLUDING CONTAMINANT COMBINATIONS FROM THE',
     *          ' CARTESIAN GAUSSIAN BASIS'/
     *       1X,'SET, THE NUMBER OF SPHERICAL HARMONICS KEPT IN',
     *          ' THE VARIATION SPACE IS',I5)
 9997 FORMAT(/5X,'DIMENSIONS OF THE SYMMETRY SUBSPACES ARE')
 9996 FORMAT(5(1X,A4,'=',I5,4X))
      END
C*MODULE SYMSLC  *DECK SPHBCK
      SUBROUTINE SPHBCK(V,T,U,NBF,DBG)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION V(NBF,NBF),T(NBF,NBF),U(NBF)
C
      LOGICAL DBG
C
      PARAMETER (MXATM=2000, MXSH=5000, MXGTOT=20000)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SPHERI/ PSHELL(3,3),DSHELL(6,6),
     *                FSHELL(10,10),GSHELL(15,15),
     *                PIHELL(3,3),DIHELL(6,6),
     *                FIHELL(10,10),GIHELL(15,15)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
C     ----- TRANSFORM -SALC- IN SPHERICAL HARMONIC SPACE -----
C           BACK TO -SALC- IN CARTESIAN GUASSIAN SPACE.
C
C     ----- SET BACK-TRANSFORMATION MATRIX -----
C
      DO 20 J=1,NUM
      DO 10 I=1,NUM
   10 T(I,J)=ZERO
   20 T(J,J)=ONE
C
      DO 1000 ISH=1,NSHELL
      IT =KTYPE(ISH)
      MIN=KMIN(ISH)
      MAX=KMAX(ISH)
      LOC=KLOC(ISH)-MIN
      GO TO (100,200,300,400,500),IT
  100 CONTINUE
      GO TO 1000
  200 IF(MIN.EQ.1) THEN
      DO 210 I=MIN+1,MAX
      DO 210 J=MIN+1,MAX
      LOCI=LOC+I
      LOCJ=LOC+J
  210 T(LOCI,LOCJ)=PSHELL(I-MIN  ,J-MIN  )
      ELSE
      DO 220 I=MIN,MAX
      DO 220 J=MIN,MAX
      LOCI=LOC+I
      LOCJ=LOC+J
  220 T(LOCI,LOCJ)=PSHELL(I-MIN+1,J-MIN+1)
      END IF
      GO TO 1000
  300 DO 310 I=MIN,MAX
      DO 310 J=MIN,MAX
      LOCI=LOC+I
      LOCJ=LOC+J
  310 T(LOCI,LOCJ)=DSHELL(I-MIN+1,J-MIN+1)
      GO TO 1000
  400 DO 410 I=MIN,MAX
      DO 410 J=MIN,MAX
      LOCI=LOC+I
      LOCJ=LOC+J
  410 T(LOCI,LOCJ)=FSHELL(I-MIN+1,J-MIN+1)
      GO TO 1000
  500 DO 510 I=MIN,MAX
      DO 510 J=MIN,MAX
      LOCI=LOC+I
      LOCJ=LOC+J
  510 T(LOCI,LOCJ)=GSHELL(I-MIN+1,J-MIN+1)
 1000 CONTINUE
C
      IF(DBG) CALL SYMPRP(T,NUM,NUM,NUM)
      IF(DBG) CALL SYMPRP(V,NUM,NUM,NUM)
C
C     ----- BACK-TRANSFORM -----
C
      DO 1030 J=1,NUM
      DO 1020 I=1,NUM
      DUM=ZERO
      DO 1010 K=1,NUM
 1010 DUM=DUM+T(I,K)*V(K,J)
 1020 U(I)=DUM
      DO 1030 I=1,NUM
 1030 V(I,J)=U(I)
C
      IF(DBG) CALL SYMPRP(V,NUM,NUM,NUM)
C
      RETURN
      END
C*MODULE SYMSLC  *DECK SPHBCKU
      SUBROUTINE SPHBCKU(T,L1,NBF,NSPH)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION T(NBF,*)
C
      PARAMETER (MXSH=5000, MXGTOT=20000)
C
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SPHERI/ PSHELL(3,3),DSHELL(6,6),
     *                FSHELL(10,10),GSHELL(15,15),
     *                PIHELL(3,3),DIHELL(6,6),
     *                FIHELL(10,10),GIHELL(15,15)
C
C     ----- CREATE TRANSFORMATION MATRIX CARTESIAN->SPHERICAL
C           A CLONE OF SPHBCK.
C     T WILL CONTAIN THE TRANSFORMATION MATRIX, NSPH - THE NUMBER OF
C     SPHERICAL HARMONICS IN THE BASIS SET, T(L1,NSPH).
C     NOTE THAT T IS NOT UNITARY, IT IS DEVISED SO THAT T^T * S * T = 1
C     (FOR S - BLOCK OF OVERLAPS FOR A FIXED SHELL).
C
      NSPH=0
      NEL=0
      CALL VCLR(T,1,L1*L1)
C
      DO 1000 ISH=1,NSHELL
         IT =KTYPE(ISH)
         MIN=KMIN(ISH)
         MAX=KMAX(ISH)
         LOC=KLOC(ISH)-MIN
C
         NSPH=NSPH+(2*IT-1)
C                                     L-SHELLS
         IF(IT.EQ.2.AND.MIN.EQ.1) NSPH=NSPH+1
C
         GO TO (100,200,300,400,500),IT
C
  100    CONTINUE
         LOCI=LOC+1
         LOCJ=LOC+1-NEL
         T(LOCI,LOCJ)=1.0D+00
         GO TO 1000
C
  200    CONTINUE
         IF(MIN.EQ.1) THEN
            LOCI=LOC+1
            LOCJ=LOC+1-NEL
            T(LOCI,LOCJ)=1.0D+00
            DO I=MIN+1,MAX
               DO J=MIN+1,MAX
                  LOCI=LOC+I
                  LOCJ=LOC+J-NEL
                  T(LOCI,LOCJ)=PSHELL(I-MIN  ,J-MIN  )
               ENDDO
            ENDDO
         ELSE
            DO I=MIN,MAX
               DO J=MIN,MAX
                  LOCI=LOC+I
                  LOCJ=LOC+J-NEL
                  T(LOCI,LOCJ)=PSHELL(I-MIN+1,J-MIN+1)
               ENDDO
            ENDDO
         END IF
         GO TO 1000
C
  300    CONTINUE
         DO I=MIN,MAX
            DO J=MIN,MIN+4
               LOCI=LOC+I
               LOCJ=LOC+J-NEL
               T(LOCI,LOCJ)=DSHELL(I-MIN+1,J-MIN+1)
            ENDDO
         ENDDO
         NEL=NEL+1
         GO TO 1000
C
  400    CONTINUE
         DO I=MIN,MAX
            DO J=MIN,MIN+6
               LOCI=LOC+I
               LOCJ=LOC+J-NEL
               T(LOCI,LOCJ)=FSHELL(I-MIN+1,J-MIN+1)
            ENDDO
         ENDDO
         NEL=NEL+3
         GO TO 1000
C
  500    CONTINUE
         DO I=MIN,MAX
            DO J=MIN,MIN+8
               LOCI=LOC+I
               LOCJ=LOC+J-NEL
               T(LOCI,LOCJ)=GSHELL(I-MIN+1,J-MIN+1)
            ENDDO
         ENDDO
         NEL=NEL+6
 1000 CONTINUE
C
      RETURN
      END
C*MODULE SYMSLC  *DECK SPHCUT
      SUBROUTINE SPHCUT(DBG)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL DBG
C
      DIMENSION MF(5),NAM(35),INDSR(10)
C
      PARAMETER (MXSH=5000, MXATM=2000)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /SPHERI/ PSHELL(3,3),DSHELL(6,6),
     *                FSHELL(10,10),GSHELL(15,15),
     *                PIHELL(3,3),DIHELL(6,6),
     *                FIHELL(10,10),GIHELL(15,15)
      COMMON /SYMSPD/ PTR(3,144),DTR(6,288),FTR(10,480),GTR(15,720)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
C
      PARAMETER (ZERO=0.0D+00)
C
      DATA MF  /0,1,4,10,20/
      DATA NAM /4H  S ,4H  X ,4H  Y ,4H  Z ,4H  XX,4H  YY,4H  ZZ,
     1          4H  XY,4H  XZ,4H  YZ,4H XXX,4H YYY,4H ZZZ,4H XXY,
     2          4H XXZ,4H YYX,4H YYZ,4H ZZX,4H ZZY,4H XYZ,4HXXXX,
     3          4HYYYY,4HZZZZ,4HXXXY,4HXXXZ,4HYYYX,4HYYYZ,4HZZZX,
     4          4HZZZY,4HXXYY,4HXXZZ,4HYYZZ,4HXXYZ,4HYYXZ,4HZZXY/
      DATA INDSR /4HD S ,4HF P1,4HF P2,4HF P3,4HG S ,
     1            4HG D1,4HG D2,4HG D3,4HG D4,4HG D5/
C
C     ----- SET UNDESIRED BLOCKS IN THE ROTATION MATRIX TO ZERO -----
C
C     ---- LOOP OVER SYMMETRY OPERATIONS -----
C
      DO 60 IT=1,NT
      ND= 6*(IT-1)
      NF=10*(IT-1)
      NG=15*(IT-1)
C
C     ----- LOOP OVER S,P,D,F,G CLASSES -----
C
      DO 20 IPER=1,3
C
      N=6
      IF(IPER.EQ.2) N=10
      IF(IPER.EQ.3) N=15
C
      GO TO (21,22,23),IPER
C
 21   DO 31 I=1,N
      DTR(I,ND+6)=ZERO
      DTR(6,ND+I)=ZERO
 31   CONTINUE
      GO TO 20
C
 22   DO 32 I=8,10
      DO 32 J=1,N
      FTR(I,NF+J)=ZERO
      FTR(J,NF+I)=ZERO
 32   CONTINUE
      GO TO 20
C
 23   DO 33 I=10,15
      DO 33 J=1,N
      GTR(I,NG+J)=ZERO
      GTR(J,NG+I)=ZERO
 33   CONTINUE
 20   CONTINUE
 60   CONTINUE
C
C     ----- ELIMINATION DONE -----
C
      IF(DBG) THEN
C
C     ----  PRINTING  ----
C
      WRITE(IW,9996)
C
      DO 120 IPER=1,3
C
      N=6
      IF(IPER.EQ.2) N=10
      IF(IPER.EQ.3) N=15
      ID=MF(IPER+2)
      GO TO (121,122,123),IPER
C
C     ---- D S ----
C
 121  DO 131 J=6,6
      WRITE(IW,9998) J,INDSR(J-6+1)
      MAXNP=0
 141  MINNP=MAXNP+1
      MAXNP=MAXNP+6
      IF(MAXNP.GT.N) MAXNP=N
      WRITE(IW,9997) (DSHELL(I,6),NAM(ID+I),I=MINNP,MAXNP)
      IF(MAXNP.LT.N) GO TO 141
 131  CONTINUE
      GO TO 120
C
C     ---- F P ----
C
 122  DO 132 J=8,10
      WRITE(IW,9998) J,INDSR(J-8+2)
      MAXNP=0
 142  MINNP=MAXNP+1
      MAXNP=MAXNP+6
      IF(MAXNP.GT.N) MAXNP=N
      WRITE(IW,9997) (FSHELL(I,J),NAM(ID+I),I=MINNP,MAXNP)
      IF(MAXNP.LT.N) GO TO 142
 132  CONTINUE
      GO TO 120
C
C     ---- G S AND G D ----
C
 123  DO 133 J=10,15
      WRITE(IW,9998) J,INDSR(J-10+5)
      MAXNP=0
 143  MINNP=MAXNP+1
      MAXNP=MAXNP+6
      IF(MAXNP.GT.N) MAXNP=N
      WRITE(IW,9997) (GSHELL(I,J),NAM(ID+I),I=MINNP,MAXNP)
      IF(MAXNP.LT.N) GO TO 143
 133  CONTINUE
 120  CONTINUE
C
      END IF
C
      RETURN
 9998 FORMAT(/,5H CHI(,I4,2X,A4,2H)=)
 9997 FORMAT(12X,6(F8.4,1X,A4))
 9996 FORMAT(/,10X,30(1H-),/,10X,'ELIMINATED VECTORS IN -SPHCUT-',/,
     1         10X,30(1H-))
      END
C*MODULE SYMSLC  *DECK SPHMAK
      SUBROUTINE SPHMAK(DBG,DD,SS,VV,TT,PP,SSPH,PSPH,DSPH,FSPH,GSPH)
C
      IMPLICIT   DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL    NORM,DBG
C
      DIMENSION DD(35),SS(35,35),VV(35,35),TT(35,35),PP(35,35),
     *          SSPH(1,1),PSPH(3,3),DSPH(6,6),FSPH(10,10),GSPH(15,15)
      DIMENSION BC(5),NX(35),NY(35),NZ(35),
     *          ISSPH(1,1),IPSPH(3,3),IDSPH(6,6),IFSPH(10,10),
     *          IGSPH(15,15)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
C
      DATA ISSPH/
     1  1/
C
      DATA IPSPH/
     1  1,  0,  0,
     2  0,  1,  0,
     3  0,  0,  1/
C
      DATA IDSPH/
     1 -1, -1,  2,  0,  0,  0,
     2  1, -1,  0,  0,  0,  0,
     3  0,  0,  0,  1,  0,  0,
     4  0,  0,  0,  0,  1,  0,
     5  0,  0,  0,  0,  0,  1,
     6  1,  1,  1,  0,  0,  0/
C
      DATA IFSPH/
     1   1,  0,  0,  0,  0, -3,  0,  0,  0,  0,
     2   0,  0,  0,  0,  1,  0, -1,  0,  0,  0,
     3  -1,  0,  0,  0,  0, -1,  0,  4,  0,  0,
     4   0,  0,  2,  0, -3,  0, -3,  0,  0,  0,
     5   0, -1,  0, -1,  0,  0,  0,  0,  4,  0,
     6   0,  0,  0,  0,  0,  0,  0,  0,  0,  1,
     7   0, -1,  0,  3,  0,  0,  0,  0,  0,  0,
     8   1,  0,  0,  0,  0,  1,  0,  1,  0,  0,
     9   0,  1,  0,  1,  0,  0,  0,  0,  1,  0,
     1   0,  0,  1,  0,  1,  0,  1,  0,  0,  0/
C
      DATA IGSPH/
     1  1,  1,  0,  0,  0,  0,  0,  0,  0, -6,  0,  0,  0,  0,  0,
     2  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0, -3,  0,
     3 -1,  1,  0,  0,  0,  0,  0,  0,  0,  0,  6, -6,  0,  0,  0,
     4  0,  0,  0,  0, -3,  0,  0,  4,  0,  0,  0,  0,  0, -3,  0,
     5  3,  3,  8,  0,  0,  0,  0,  0,  0,  6,-24,-24,  0,  0,  0,
     6  0,  0,  0,  0,  0,  0, -3,  0,  4,  0,  0,  0, -3,  0,  0,
     7  0,  0,  0, -1,  0, -1,  0,  0,  0,  0,  0,  0,  0,  0,  6,
     8  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0, -3,  0,  0,
     9  0,  0,  0,  1,  0, -1,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     1  1,  1,  1,  0,  0,  0,  0,  0,  0,  2,  2,  2,  0,  0,  0,
     2 -1, -1,  2,  0,  0,  0,  0,  0,  0, -2,  1,  1,  0,  0,  0,
     3  1, -1,  0,  0,  0,  0,  0,  0,  0,  0,  1, -1,  0,  0,  0,
     4  0,  0,  0,  1,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  1,
     5  0,  0,  0,  0,  1,  0,  0,  1,  0,  0,  0,  0,  0,  1,  0,
     6  0,  0,  0,  0,  0,  0,  1,  0,  1,  0,  0,  0,  1,  0,  0/
C
      DATA  NX    / 1, 2, 1, 1, 3, 1, 1, 2, 2, 1,
     1              4, 1, 1, 3, 3, 2, 1, 2, 1, 2,
     2              5, 1, 1, 4, 4, 2, 1, 2, 1, 3,
     3              3, 1, 3, 2, 2/
      DATA  NY    / 1, 1, 2, 1, 1, 3, 1, 2, 1, 2,
     1              1, 4, 1, 2, 1, 3, 3, 1, 2, 2,
     2              1, 5, 1, 2, 1, 4, 4, 1, 2, 3,
     3              1, 3, 2, 3, 2/
      DATA  NZ    / 1, 1, 1, 2, 1, 1, 3, 1, 2, 2,
     1              1, 1, 4, 1, 2, 1, 2, 3, 3, 2,
     2              1, 1, 5, 1, 2, 1, 2, 4, 4, 1,
     3              3, 3, 2, 2, 3/
      DATA ZERO   /  0.0D+00/
      DATA BC     /5*0.0D+00/
C
      DBG=.FALSE.
C
      BC(1)=1
      BC(2)=1
      BC(3)=1*3
      BC(4)=1*3*5
      BC(5)=1*3*5*7
C
      NORM =NORMF.EQ.0.OR.NORMP.EQ.0
C
      NBF = 35
      DO J=1,NBF
         DO I=1,NBF
            IJX=NX(I)+NX(J)-2
            IJY=NY(I)+NY(J)-2
            IJZ=NZ(I)+NZ(J)-2
            IF(((IJX/2)*2.NE.IJX).OR.((IJY/2)*2.NE.IJY)
     1                           .OR.((IJZ/2)*2.NE.IJZ)
     2                                                 ) THEN
               SS(I,J)=ZERO
            ELSE
               MX=IJX/2
               MY=IJY/2
               MZ=IJZ/2
               IF(DBG) THEN
                  WRITE(IW,9997) I,J,MX,MY,MZ,
     *               INT(BC(MX+1)),INT(BC(MY+1)),INT(BC(MZ+1))
               END IF
               SIJ =(BC(MX+1)*BC(MY+1)*BC(MZ+1))/(2**MX*2**MY*2**MZ)
               SS(I,J)=SIJ
            END IF
         ENDDO
      ENDDO
C
C     ----- DEAL WITH NORMALIZATION OF -AO- BASIS -----
C
      IF(NORM) THEN
         DO I=1,NBF
            DD(I)=SQRT(SS(I,I))
         ENDDO
         DO J=1,NBF
            DO I=1,NBF
               SS(I,J)=SS(I,J)/(DD(I)*DD(J))
            ENDDO
         ENDDO
      END IF
C
C     ----- PRINT NORMALIZATION -----
C
      IF(DBG) THEN
         WRITE(IW,*) '-AO- NORMALIZATION'
         JMAX=0
   10    JMIN=JMAX+1
         JMAX=JMAX+10
         IF(JMAX.GT.NBF) JMAX=NBF
         WRITE(IW,9999)
            WRITE(IW,9998) (DD(J),J=JMIN,JMAX)
         IF(JMAX.LT.NBF) GO TO 10
         WRITE(IW,9999)
      END IF
C
C     ----- PRINT -----
C
      IF(DBG) THEN
         WRITE(IW,*) '- AO- BASIS OVERLAP'
         JMAX=0
   20    JMIN=JMAX+1
         JMAX=JMAX+10
         IF(JMAX.GT.NBF) JMAX=NBF
         WRITE(IW,9999)
         DO I=1,NBF
            WRITE(IW,9998) (SS(I,J),J=JMIN,JMAX)
         ENDDO
         IF(JMAX.LT.NBF) GO TO 20
         WRITE(IW,9999)
      END IF
C
C     ----- BUILD SPHERICAL HARMONICS -----
C
      DO J=1,35
         DO I=1,35
            VV(I,J)=ZERO
         ENDDO
      ENDDO
C
      DO J=1,1
         DO I=1,1
            VV(I   ,J   )=ISSPH(I,J)
         ENDDO
      ENDDO
      DO J=1,3
         DO I=1,3
            VV(I+ 1,J+ 1)=IPSPH(I,J)
         ENDDO
      ENDDO
      DO J=1,6
         DO I=1,6
            VV(I+ 4,J+ 4)=IDSPH(I,J)
         ENDDO
      ENDDO
      DO J=1,10
         DO I=1,10
            VV(I+10,J+10)=IFSPH(I,J)
         ENDDO
      ENDDO
      DO J=1,15
         DO I=1,15
            VV(I+20,J+20)=IGSPH(I,J)
         ENDDO
      ENDDO
C
C     ----- DEAL WITH NORMALIZATION OF -AO- BASIS -----
C
      IF(NORM) THEN
         DO J=1,NBF
            DO I=1,NBF
               VV(I,J)=VV(I,J)*DD(I)
            ENDDO
         ENDDO
      END IF
C
C     ----- PRINT -----
C
      IF(DBG) THEN
         WRITE(IW,*) '-SPH- DEFINITION'
         JMAX=0
   30    JMIN=JMAX+1
         JMAX=JMAX+10
         IF(JMAX.GT.NBF) JMAX=NBF
         WRITE(IW,9999)
         DO I=1,NBF
            WRITE(IW,9998) (VV(I,J),J=JMIN,JMAX)
         ENDDO
         IF(JMAX.LT.NBF) GO TO 30
         WRITE(IW,9999)
      END IF
C
C     ----- TRANSFORM OVERLAP MATRIX -----
C
      DO I=1,NBF
         DO J=1,NBF
            DUM=ZERO
            DO K=1,NBF
               DUM=DUM+SS(I,K)*VV(K,J)
            ENDDO
            TT(I,J)=DUM
         ENDDO
      ENDDO
      DO I=1,NBF
         DO J=1,NBF
            DUM=ZERO
            DO K=1,NBF
               DUM=DUM+VV(K,I)*TT(K,J)
            ENDDO
            PP(I,J)=DUM
         ENDDO
      ENDDO
C
C     ----- PRINT -----
C
      IF(DBG) THEN
         WRITE(IW,*) '-SPH- BASIS OVERLAP'
         JMAX=0
   40    JMIN=JMAX+1
         JMAX=JMAX+10
         IF(JMAX.GT.NBF) JMAX=NBF
         WRITE(IW,9999)
         DO I=1,NBF
            WRITE(IW,9998) (PP(I,J),J=JMIN,JMAX)
         ENDDO
         IF(JMAX.LT.NBF) GO TO 40
         WRITE(IW,9999)
      END IF
C
C     ----- NORMALIZE FURTHER -----
C
      DO J=1,NBF
         DUM=SQRT(PP(J,J))
         DO I=1,NBF
            VV(I,J)=VV(I,J)/DUM
         ENDDO
      ENDDO
C
C     ----- PRINT -----
C
      IF(DBG) THEN
         WRITE(IW,*) '-SPH- NORMALIZED'
         JMAX=0
   50    JMIN=JMAX+1
         JMAX=JMAX+10
         IF(JMAX.GT.NBF) JMAX=NBF
         WRITE(IW,9999)
         DO I=1,NBF
            WRITE(IW,9998) (VV(I,J),J=JMIN,JMAX)
         ENDDO
         IF(JMAX.LT.NBF) GO TO 50
         WRITE(IW,9999)
      END IF
C
C     ----- TRANSFORM OVERLAP MATRIX -----
C
      DO I=1,NBF
         DO J=1,NBF
            DUM=ZERO
            DO K=1,NBF
               DUM=DUM+SS(I,K)*VV(K,J)
            ENDDO
            TT(I,J)=DUM
         ENDDO
      ENDDO
      DO I=1,NBF
         DO J=1,NBF
            DUM=ZERO
            DO K=1,NBF
               DUM=DUM+VV(K,I)*TT(K,J)
            ENDDO
            PP(I,J)=DUM
         ENDDO
      ENDDO
C
C     ----- PRINT -----
C
      IF(DBG) THEN
         WRITE(IW,*) '-SPH- BASIS OVERLAP'
         JMAX=0
   60    JMIN=JMAX+1
         JMAX=JMAX+10
         IF(JMAX.GT.NBF) JMAX=NBF
         WRITE(IW,9999)
         DO I=1,NBF
            WRITE(IW,9998) (PP(I,J),J=JMIN,JMAX)
         ENDDO
         IF(JMAX.LT.NBF) GO TO 60
         WRITE(IW,9999)
      END IF
C
C     ----- BUILD -SSPH- -PSPH- -DSPH- -FSPH- -GSPH- -----
C
      DO J=1,1
         DO I=1,1
            SSPH(I,J)=VV(I   ,J   )
         ENDDO
      ENDDO
      DO J=1,3
         DO I=1,3
            PSPH(I,J)=VV(I+ 1,J+ 1)
         ENDDO
      ENDDO
      DO J=1,6
         DO I=1,6
            DSPH(I,J)=VV(I+ 4,J+ 4)
         ENDDO
      ENDDO
      DO J=1,10
         DO I=1,10
            FSPH(I,J)=VV(I+10,J+10)
         ENDDO
      ENDDO
      DO J=1,15
         DO I=1,15
            GSPH(I,J)=VV(I+20,J+20)
         ENDDO
      ENDDO
      RETURN
C
 9999 FORMAT(/)
 9998 FORMAT(10F12.6)
 9997 FORMAT(5I3,3X,3I5)
      END
C*MODULE SYMSLC  *DECK SPHMAT
      SUBROUTINE SPHMAT(SP,NDIM,IFINV,IFTSP)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXSH=5000, MXGTOT=20000)
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SPHERI/ PSHELL(3,3),DSHELL(6,6),
     *                FSHELL(10,10),GSHELL(15,15),
     *                PIHELL(3,3),DIHELL(6,6),
     *                FIHELL(10,10),GIHELL(15,15)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
      DIMENSION SP(NDIM,NDIM)
C
C
      DO J=1,NUM
         DO I=1,NUM
            SP(I,J)=ZERO
         ENDDO
         SP(J,J)=ONE
      ENDDO
C
C     ----- BRANCH ON -IFINV- AND -IFTSP- -----
C
      IF(IFINV.EQ.-1.AND.IFTSP.EQ.0) THEN
C
C     ----- INVERSE -----
C
         DO II=1,NSHELL
            LIT =KTYPE(II)
            MIN =KMIN(II)
            LOC =KLOC(II)
            IF(LIT.EQ.1) THEN
            ELSE IF(LIT.EQ.2) THEN
               IF(MIN.EQ.1) THEN
                  DO J=1,3
                     DO I=1,3
                        SP(I+LOC  ,J+LOC  )=PIHELL(I,J)
                     ENDDO
                  ENDDO
               ELSE
                  DO J=1,3
                     DO I=1,3
                        SP(I+LOC-1,J+LOC-1)=PIHELL(I,J)
                     ENDDO
                  ENDDO
               END IF
            ELSE IF(LIT.EQ.3) THEN
               DO J=1,6
                  DO I=1,6
                     SP(I+LOC-1,J+LOC-1)=DIHELL(I,J)
                  ENDDO
               ENDDO
            ELSE IF(LIT.EQ.4) THEN
               DO J=1,10
                  DO I=1,10
                     SP(I+LOC-1,J+LOC-1)=FIHELL(I,J)
                  ENDDO
               ENDDO
            ELSE IF(LIT.EQ.5) THEN
               DO J=1,15
                  DO I=1,15
                     SP(I+LOC-1,J+LOC-1)=GIHELL(I,J)
                  ENDDO
               ENDDO
            ELSE
               WRITE(IW,9999)
               CALL ABRT
            END IF
         ENDDO
C
      ELSE IF(IFINV.EQ.-1.AND.IFTSP.NE.0) THEN
C
C     ----- INVERSE TRANSPOSE -----
C
         DO II=1,NSHELL
            LIT =KTYPE(II)
            MIN =KMIN(II)
            LOC =KLOC(II)
            IF(LIT.EQ.1) THEN
            ELSE IF(LIT.EQ.2) THEN
               IF(MIN.EQ.1) THEN
                  DO J=1,3
                     DO I=1,3
                        SP(J+LOC  ,I+LOC  )=PIHELL(I,J)
                     ENDDO
                  ENDDO
               ELSE
                  DO J=1,3
                     DO I=1,3
                        SP(J+LOC-1,I+LOC-1)=PIHELL(I,J)
                     ENDDO
                  ENDDO
               END IF
            ELSE IF(LIT.EQ.3) THEN
               DO J=1,6
                  DO I=1,6
                     SP(J+LOC-1,I+LOC-1)=DIHELL(I,J)
                  ENDDO
               ENDDO
            ELSE IF(LIT.EQ.4) THEN
               DO J=1,10
                  DO I=1,10
                     SP(J+LOC-1,I+LOC-1)=FIHELL(I,J)
                  ENDDO
               ENDDO
            ELSE IF(LIT.EQ.5) THEN
               DO J=1,15
                  DO I=1,15
                     SP(J+LOC-1,I+LOC-1)=GIHELL(I,J)
                  ENDDO
               ENDDO
            ELSE
               WRITE(IW,9999)
               CALL ABRT
            END IF
         ENDDO
C
      ELSE IF(IFINV.NE.-1.AND.IFTSP.EQ.0) THEN
C
C     -----                   -----
C
         DO II=1,NSHELL
            LIT =KTYPE(II)
            MIN =KMIN(II)
            LOC =KLOC(II)
            IF(LIT.EQ.1) THEN
            ELSE IF(LIT.EQ.2) THEN
               IF(MIN.EQ.1) THEN
                  DO J=1,3
                     DO I=1,3
                        SP(I+LOC  ,J+LOC  )=PSHELL(I,J)
                     ENDDO
                  ENDDO
               ELSE
                  DO J=1,3
                     DO I=1,3
                        SP(I+LOC-1,J+LOC-1)=PSHELL(I,J)
                     ENDDO
                  ENDDO
               END IF
            ELSE IF(LIT.EQ.3) THEN
               DO J=1,6
                  DO I=1,6
                     SP(I+LOC-1,J+LOC-1)=DSHELL(I,J)
                  ENDDO
               ENDDO
            ELSE IF(LIT.EQ.4) THEN
               DO J=1,10
                  DO I=1,10
                     SP(I+LOC-1,J+LOC-1)=FSHELL(I,J)
                  ENDDO
               ENDDO
            ELSE IF(LIT.EQ.5) THEN
               DO J=1,15
                  DO I=1,15
                     SP(I+LOC-1,J+LOC-1)=GSHELL(I,J)
                  ENDDO
               ENDDO
            ELSE
               WRITE(IW,9999)
               CALL ABRT
            END IF
         ENDDO
C
      ELSE IF(IFINV.NE.-1.AND.IFTSP.NE.0) THEN
C
C     -----         TRANSPOSE -----
C
         DO II=1,NSHELL
            LIT =KTYPE(II)
            MIN =KMIN(II)
            LOC =KLOC(II)
            IF(LIT.EQ.1) THEN
            ELSE IF(LIT.EQ.2) THEN
               IF(MIN.EQ.1) THEN
                  DO J=1,3
                     DO I=1,3
                        SP(J+LOC  ,I+LOC  )=PSHELL(I,J)
                     ENDDO
                  ENDDO
               ELSE
                  DO J=1,3
                     DO I=1,3
                        SP(J+LOC-1,I+LOC-1)=PSHELL(I,J)
                     ENDDO
                  ENDDO
               END IF
            ELSE IF(LIT.EQ.3) THEN
               DO J=1,6
                  DO I=1,6
                     SP(J+LOC-1,I+LOC-1)=DSHELL(I,J)
                  ENDDO
               ENDDO
            ELSE IF(LIT.EQ.4) THEN
               DO J=1,10
                  DO I=1,10
                     SP(J+LOC-1,I+LOC-1)=FSHELL(I,J)
                  ENDDO
               ENDDO
            ELSE IF(LIT.EQ.5) THEN
               DO J=1,15
                  DO I=1,15
                     SP(J+LOC-1,I+LOC-1)=GSHELL(I,J)
                  ENDDO
               ENDDO
            ELSE
               WRITE(IW,9999)
               CALL ABRT
            END IF
         ENDDO
C
      END IF
C
      RETURN
 9999 FORMAT(' -SPHMAT- LIMITED TO -G- SHELLS . STOP .')
C 9998 FORMAT(/,1X,'SPHERICAL HARMONICS TRANSFORMATION MATRIX',
C     1            ' -IFINV- = ',I5,' -IFTSP- = ',I5,
C     2       /,1X,41(1H-))
      END
C*MODULE SYMSLC  *DECK SPHSET
      SUBROUTINE SPHSET(DBG,DD,SS,VV,TT,PP,ROT,BUFF,
     *                  PSPH,DSPH,FSPH,GSPH)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DOUBLE PRECISION DD(35),ROT(15,15),BUFF(15,15),
     *          SS(35,35),VV(35,35),TT(35,35),PP(35,35),
     *          PSPH(3,3),DSPH(6,6),FSPH(10,10),GSPH(15,15)
      DOUBLE PRECISION INDX(15),VWRK(15)
C
      LOGICAL DBG
C
      COMMON /FRAME / U1,U2,U3,V1,V2,V3,W1,W2,W3,X0,Y0,Z0
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /SPHERI/ PSHELL(3,3),DSHELL(6,6),
     *                FSHELL(10,10),GSHELL(15,15),
     *                PIHELL(3,3),DIHELL(6,6),
     *                FIHELL(10,10),GIHELL(15,15)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
C         THE NEXT IS FOR SUN 32 BIT F77 COMPILER VERSION 4.2
      SAVE INDX,VWRK,SSPH,A,D,I,J,K,IERR
C
C     -----  THIS ROUTINE DEFINES THE SPHERICAL HARMONICS -----
C     -----  AND ITS INVERSE.                             -----
C
      IF(DBG) WRITE(IW,9999)
C
C     ----- INITIALIZATIONS -----
C
      DO 1 I=1,15
       INDX(I)=0
 1    CONTINUE
      DO 10 I=1,3
      DO  5 J=1,3
        PSPH(I,J)=ZERO
      PSHELL(I,J)=ZERO
  5   PIHELL(I,J)=ZERO
 10   PIHELL(I,I)=ONE
      DO 20 I=1,6
      DO 15 J=1,6
        DSPH(I,J)=ZERO
      DSHELL(I,J)=ZERO
 15   DIHELL(I,J)=ZERO
 20   DIHELL(I,I)=ONE
      DO 30 I=1,10
      DO 25 J=1,10
        FSPH(I,J)=ZERO
      FSHELL(I,J)=ZERO
 25   FIHELL(I,J)=ZERO
 30   FIHELL(I,I)=ONE
      DO 40 I=1,15
      DO 35 J=1,15
        GSPH(I,J)=ZERO
      GSHELL(I,J)=ZERO
 35   GIHELL(I,J)=ZERO
 40   GIHELL(I,I)=ONE
C
C     ----- GET SPHERICAL HARMONICS DEFINITIONS -----
C
      CALL SPHMAK(DBG,DD,SS,VV,TT,PP,SSPH,PSPH,DSPH,FSPH,GSPH)
C
C     ----- P SPHERICAL HARMONICS -----
C
      IF(DBG) CALL SYMPRP(PSPH, 3, 3, 3)
C
C     ----- LOAD TRANSFORMATION OF FRAMES IN -ROT- -----
C
C     NEW X
      ROT(1,1)=U1
      ROT(2,1)=U2
      ROT(3,1)=U3
C     NEW Y
      ROT(1,2)=V1
      ROT(2,2)=V2
      ROT(3,2)=V3
C     NEW Z
      ROT(1,3)=W1
      ROT(2,3)=W2
      ROT(3,3)=W3
C
C     ----- EXPRESS SPHERICAL HARMONICS IN ABSOLUTE FRAME -----
C
      DO 80 I=1,3
      DO 80 J=1,3
      A=ZERO
      DO 70 K=1,3
      A=A+ROT(I,K)*PSPH(K,J)
  70  CONTINUE
        BUFF(I,J)=A
  80  PSHELL(I,J)=A
C
C     ----- DEFINE INVERSE TRANSFORMATION -----
C
      CALL LUDCMP(BUFF,15, 3,INDX,VWRK,D,IERR)
      DO I=1,3
         CALL LUBKSB(BUFF,15, 3,INDX,PIHELL(1,I))
      ENDDO
C
C     PRINT
C
      IF(DBG) CALL SYMPRP(PSHELL, 3, 3, 3)
      IF(DBG) CALL SYMPRP(PIHELL, 3, 3, 3)
C
C     ----- D SPHERICAL HARMONICS -----
C
      IF(DBG) CALL SYMPRP(DSPH, 6, 6, 6)
C
C     ----- LOAD TRANSFORMATION OF FRAMES IN -ROT- -----
C
C     NEW X2
      ROT(1,1)=U1*U1
      ROT(2,1)=U2*U2
      ROT(3,1)=U3*U3
      ROT(4,1)=U1*U2+U1*U2
      ROT(5,1)=U1*U3+U1*U3
      ROT(6,1)=U2*U3+U2*U3
C     NEW Y2
      ROT(1,2)=V1*V1
      ROT(2,2)=V2*V2
      ROT(3,2)=V3*V3
      ROT(4,2)=V1*V2+V1*V2
      ROT(5,2)=V1*V3+V1*V3
      ROT(6,2)=V2*V3+V2*V3
C     NEW Z2
      ROT(1,3)=W1*W1
      ROT(2,3)=W2*W2
      ROT(3,3)=W3*W3
      ROT(4,3)=W1*W2+W1*W2
      ROT(5,3)=W1*W3+W1*W3
      ROT(6,3)=W2*W3+W2*W3
C     NEW XY
      ROT(1,4)=U1*V1
      ROT(2,4)=U2*V2
      ROT(3,4)=U3*V3
      ROT(4,4)=U1*V2+U2*V1
      ROT(5,4)=U1*V3+U3*V1
      ROT(6,4)=U2*V3+U3*V2
C     NEW XZ
      ROT(1,5)=U1*W1
      ROT(2,5)=U2*W2
      ROT(3,5)=U3*W3
      ROT(4,5)=U1*W2+U2*W1
      ROT(5,5)=U1*W3+U3*W1
      ROT(6,5)=U2*W3+U3*W2
C     NEW YZ
      ROT(1,6)=V1*W1
      ROT(2,6)=V2*W2
      ROT(3,6)=V3*W3
      ROT(4,6)=V1*W2+V2*W1
      ROT(5,6)=V1*W3+V3*W1
      ROT(6,6)=V2*W3+V3*W2
C
C     ----- EXPRESS SPHERICAL HARMONICS IN ABSOLUTE FRAME -----
C
      DO 160 I=1,6
      DO 160 J=1,6
      A=ZERO
      DO 150 K=1,6
      A=A+ROT(I,K)*DSPH(K,J)
 150  CONTINUE
        BUFF(I,J)=A
 160  DSHELL(I,J)=A
C
C     ----- DEFINE INVERSE TRANSFORMATION -----
C
      CALL LUDCMP(BUFF,15, 6,INDX,VWRK,D,IERR)
      DO 180 I=1,6
      CALL LUBKSB(BUFF,15, 6,INDX,DIHELL(1,I))
 180  CONTINUE
C
C     PRINT
C
      IF(DBG) CALL SYMPRP(DSHELL,6,6,6)
      IF(DBG) CALL SYMPRP(DIHELL,6,6,6)
C
C     ----- F SPHERICAL HARMONICS -----
C
      IF(DBG) CALL SYMPRP(FSPH,10,10,10)
C
C     ----- LOAD TRANSFORMATION OF FRAMES IN -ROT- -----
C
C     NEW XXX
      ROT(1,1) =U1*U1*U1
      ROT(2,1) =U2*U2*U2
      ROT(3,1) =U3*U3*U3
      ROT(4,1) =U2*U1*U1+ U1*U2*U1+ U1*U1*U2
      ROT(5,1) =U3*U1*U1+ U1*U3*U1+ U1*U1*U3
      ROT(6,1) =U1*U2*U2+ U2*U1*U2+ U2*U2*U1
      ROT(7,1) =U3*U2*U2+ U2*U3*U2+ U2*U2*U3
      ROT(8,1) =U1*U3*U3+ U3*U1*U3+ U3*U3*U1
      ROT(9,1) =U2*U3*U3+ U3*U2*U3+ U3*U3*U2
      ROT(10,1)=U1*U2*U3+ U1*U3*U2+ U2*U1*U3+ U2*U3*U1+
     1            U3*U1*U2+ U3*U2*U1
C     NEW YYY
      ROT(1,2) =V1*V1*V1
      ROT(2,2) =V2*V2*V2
      ROT(3,2) =V3*V3*V3
      ROT(4,2) =V2*V1*V1+ V1*V2*V1+ V1*V1*V2
      ROT(5,2) =V3*V1*V1+ V1*V3*V1+ V1*V1*V3
      ROT(6,2) =V1*V2*V2+ V2*V1*V2+ V2*V2*V1
      ROT(7,2) =V3*V2*V2+ V2*V3*V2+ V2*V2*V3
      ROT(8,2) =V1*V3*V3+ V3*V1*V3+ V3*V3*V1
      ROT(9,2) =V2*V3*V3+ V3*V2*V3+ V3*V3*V2
      ROT(10,2)=V1*V2*V3+ V1*V3*V2+ V2*V1*V3+ V2*V3*V1+
     1            V3*V1*V2+ V3*V2*V1
C     NEW ZZZ
      ROT(1,3) =W1*W1*W1
      ROT(2,3) =W2*W2*W2
      ROT(3,3) =W3*W3*W3
      ROT(4,3) =W2*W1*W1+ W1*W2*W1+ W1*W1*W2
      ROT(5,3) =W3*W1*W1+ W1*W3*W1+ W1*W1*W3
      ROT(6,3) =W1*W2*W2+ W2*W1*W2+ W2*W2*W1
      ROT(7,3) =W3*W2*W2+ W2*W3*W2+ W2*W2*W3
      ROT(8,3) =W1*W3*W3+ W3*W1*W3+ W3*W3*W1
      ROT(9,3) =W2*W3*W3+ W3*W2*W3+ W3*W3*W2
      ROT(10,3)=W1*W2*W3+ W1*W3*W2+ W2*W1*W3+ W2*W3*W1+
     1            W3*W1*W2+ W3*W2*W1
C     NEW XXY
      ROT(1,4) =U1*U1*V1
      ROT(2,4) =U2*U2*V2
      ROT(3,4) =U3*U3*V3
      ROT(4,4) =U2*U1*V1+ U1*U2*V1+ U1*U1*V2
      ROT(5,4) =U3*U1*V1+ U1*U3*V1+ U1*U1*V3
      ROT(6,4) =U1*U2*V2+ U2*U1*V2+ U2*U2*V1
      ROT(7,4) =U3*U2*V2+ U2*U3*V2+ U2*U2*V3
      ROT(8,4) =U1*U3*V3+ U3*U1*V3+ U3*U3*V1
      ROT(9,4) =U2*U3*V3+ U3*U2*V3+ U3*U3*V2
      ROT(10,4)=U1*U2*V3+ U1*U3*V2+ U2*U1*V3+ U2*U3*V1+
     1            U3*U1*V2+ U3*U2*V1
C     NEW XXZ
      ROT(1,5) =U1*U1*W1
      ROT(2,5) =U2*U2*W2
      ROT(3,5) =U3*U3*W3
      ROT(4,5) =U2*U1*W1+ U1*U2*W1+ U1*U1*W2
      ROT(5,5) =U3*U1*W1+ U1*U3*W1+ U1*U1*W3
      ROT(6,5) =U1*U2*W2+ U2*U1*W2+ U2*U2*W1
      ROT(7,5) =U3*U2*W2+ U2*U3*W2+ U2*U2*W3
      ROT(8,5) =U1*U3*W3+ U3*U1*W3+ U3*U3*W1
      ROT(9,5) =U2*U3*W3+ U3*U2*W3+ U3*U3*W2
      ROT(10,5)=U1*U2*W3+ U1*U3*W2+ U2*U1*W3+ U2*U3*W1+
     1            U3*U1*W2+ U3*U2*W1
C     NEW YYX
      ROT(1,6) =V1*V1*U1
      ROT(2,6) =V2*V2*U2
      ROT(3,6) =V3*V3*U3
      ROT(4,6) =V2*V1*U1+ V1*V2*U1+ V1*V1*U2
      ROT(5,6) =V3*V1*U1+ V1*V3*U1+ V1*V1*U3
      ROT(6,6) =V1*V2*U2+ V2*V1*U2+ V2*V2*U1
      ROT(7,6) =V3*V2*U2+ V2*V3*U2+ V2*V2*U3
      ROT(8,6) =V1*V3*U3+ V3*V1*U3+ V3*V3*U1
      ROT(9,6) =V2*V3*U3+ V3*V2*U3+ V3*V3*U2
      ROT(10,6)=V1*V2*U3+ V1*V3*U2+ V2*V1*U3+ V2*V3*U1+
     1            V3*V1*U2+ V3*V2*U1
C     NEW YYZ
      ROT(1,7) =V1*V1*W1
      ROT(2,7) =V2*V2*W2
      ROT(3,7) =V3*V3*W3
      ROT(4,7) =V2*V1*W1+ V1*V2*W1+ V1*V1*W2
      ROT(5,7) =V3*V1*W1+ V1*V3*W1+ V1*V1*W3
      ROT(6,7) =V1*V2*W2+ V2*V1*W2+ V2*V2*W1
      ROT(7,7) =V3*V2*W2+ V2*V3*W2+ V2*V2*W3
      ROT(8,7) =V1*V3*W3+ V3*V1*W3+ V3*V3*W1
      ROT(9,7) =V2*V3*W3+ V3*V2*W3+ V3*V3*W2
      ROT(10,7)=V1*V2*W3+ V1*V3*W2+ V2*V1*W3+ V2*V3*W1+
     1            V3*V1*W2+ V3*V2*W1
C     NEW ZZX
      ROT(1,8) =W1*W1*U1
      ROT(2,8) =W2*W2*U2
      ROT(3,8) =W3*W3*U3
      ROT(4,8) =W2*W1*U1+ W1*W2*U1+ W1*W1*U2
      ROT(5,8) =W3*W1*U1+ W1*W3*U1+ W1*W1*U3
      ROT(6,8) =W1*W2*U2+ W2*W1*U2+ W2*W2*U1
      ROT(7,8) =W3*W2*U2+ W2*W3*U2+ W2*W2*U3
      ROT(8,8) =W1*W3*U3+ W3*W1*U3+ W3*W3*U1
      ROT(9,8) =W2*W3*U3+ W3*W2*U3+ W3*W3*U2
      ROT(10,8)=W1*W2*U3+ W1*W3*U2+ W2*W1*U3+ W2*W3*U1+
     1            W3*W1*U2+ W3*W2*U1
C     NEW ZZY
      ROT(1,9) =W1*W1*V1
      ROT(2,9) =W2*W2*V2
      ROT(3,9) =W3*W3*V3
      ROT(4,9) =W2*W1*V1+ W1*W2*V1+ W1*W1*V2
      ROT(5,9) =W3*W1*V1+ W1*W3*V1+ W1*W1*V3
      ROT(6,9) =W1*W2*V2+ W2*W1*V2+ W2*W2*V1
      ROT(7,9) =W3*W2*V2+ W2*W3*V2+ W2*W2*V3
      ROT(8,9) =W1*W3*V3+ W3*W1*V3+ W3*W3*V1
      ROT(9,9) =W2*W3*V3+ W3*W2*V3+ W3*W3*V2
      ROT(10,9)=W1*W2*V3+ W1*W3*V2+ W2*W1*V3+ W2*W3*V1+
     1            W3*W1*V2+ W3*W2*V1
C     NEW XYZ
      ROT(1,10) =U1*V1*W1
      ROT(2,10) =U2*V2*W2
      ROT(3,10) =U3*V3*W3
      ROT(4,10) =U2*V1*W1+ U1*V2*W1+ U1*V1*W2
      ROT(5,10) =U3*V1*W1+ U1*V3*W1+ U1*V1*W3
      ROT(6,10) =U1*V2*W2+ U2*V1*W2+ U2*V2*W1
      ROT(7,10) =U3*V2*W2+ U2*V3*W2+ U2*V2*W3
      ROT(8,10) =U1*V3*W3+ U3*V1*W3+ U3*V3*W1
      ROT(9,10) =U2*V3*W3+ U3*V2*W3+ U3*V3*W2
      ROT(10,10)=U1*V2*W3+ U1*V3*W2+ U2*V1*W3+ U2*V3*W1+
     1             U3*V1*W2+ U3*V2*W1
C
C     ----- EXPRESS SPHERICAL HARMONICS IN ABSOLUTE FRAME -----
C
      DO 260 I=1,10
      DO 260 J=1,10
      A=ZERO
      DO 250 K=1,10
      A=A+ROT(I,K)*FSPH(K,J)
 250  CONTINUE
        BUFF(I,J)=A
 260  FSHELL(I,J)=A
C
C     ----- DEFINE INVERSE TRANSFORMATION -----
C
      CALL LUDCMP(BUFF,15,10,INDX,VWRK,D,IERR)
      DO 280 I=1,10
      CALL LUBKSB(BUFF,15,10,INDX,FIHELL(1,I))
 280  CONTINUE
C
C     PRINT
C
      IF(DBG) CALL SYMPRP(FSHELL,10,10,10)
      IF(DBG) CALL SYMPRP(FIHELL,10,10,10)
C
C     ----- G SPHERICAL HARMONICS -----
C
      IF(DBG) CALL SYMPRP(GSPH,15,15,15)
C
C     ----- LOAD TRANSFORMATION OF FRAMES IN -ROT- -----
C
C     NEW XXXX
      ROT(1,1)  =U1*U1*U1*U1
      ROT(2,1)  =U2*U2*U2*U2
      ROT(3,1)  =U3*U3*U3*U3
      ROT(4,1)  =U2*U1*U1*U1+ U1*U2*U1*U1+ U1*U1*U2*U1+ U1*U1*U1*U2
      ROT(5,1)  =U3*U1*U1*U1+ U1*U3*U1*U1+ U1*U1*U3*U1+ U1*U1*U1*U3
      ROT(6,1)  =U1*U2*U2*U2+ U2*U1*U2*U2+ U2*U2*U1*U2+ U2*U2*U2*U1
      ROT(7,1)  =U3*U2*U2*U2+ U2*U3*U2*U2+ U2*U2*U3*U2+ U2*U2*U2*U3
      ROT(8,1)  =U1*U3*U3*U3+ U3*U1*U3*U3+ U3*U3*U1*U3+ U3*U3*U3*U1
      ROT(9,1)  =U2*U3*U3*U3+ U3*U2*U3*U3+ U3*U3*U2*U3+ U3*U3*U3*U2
      ROT(10,1) =U1*U1*U2*U2+ U2*U2*U1*U1+ U1*U2*U1*U2+ U1*U2*U2*U1
     1           + U2*U1*U1*U2+ U2*U1*U2*U1
      ROT(11,1) =U1*U1*U3*U3+ U3*U3*U1*U1+ U1*U3*U1*U3+ U1*U3*U3*U1
     1           + U3*U1*U1*U3+ U3*U1*U3*U1
      ROT(12,1) =U2*U2*U3*U3+ U3*U3*U2*U2+ U2*U3*U2*U3+ U2*U3*U3*U2
     1           + U3*U2*U2*U3+ U3*U2*U3*U2
      ROT(13,1) =U1*U1*U2*U3+ U1*U1*U3*U2+ U1*U2*U1*U3+ U1*U3*U1*U2
     1           + U1*U2*U3*U1+ U1*U3*U2*U1+ U2*U1*U1*U3+ U3*U1*U1*U2
     2           + U2*U1*U3*U1+ U3*U1*U2*U1+ U2*U3*U1*U1+ U3*U2*U1*U1
      ROT(14,1) =U2*U2*U1*U3+ U2*U2*U3*U1+ U2*U1*U2*U3+ U2*U3*U2*U1
     1           + U2*U1*U3*U2+ U2*U3*U1*U2+ U1*U2*U2*U3+ U3*U2*U2*U1
     2           + U1*U2*U3*U2+ U3*U2*U1*U2+ U1*U3*U2*U2+ U3*U1*U2*U2
      ROT(15,1) =U3*U3*U1*U2+ U3*U3*U2*U1+ U3*U1*U3*U2+ U3*U2*U3*U1
     1           + U3*U1*U2*U3+ U3*U2*U1*U3+ U1*U3*U3*U2+ U2*U3*U3*U1
     2           + U1*U3*U2*U3+ U2*U3*U1*U3+ U1*U2*U3*U3+ U2*U1*U3*U3
C     NEW YYYY
      ROT(1,2)  =V1*V1*V1*V1
      ROT(2,2)  =V2*V2*V2*V2
      ROT(3,2)  =V3*V3*V3*V3
      ROT(4,2)  =V2*V1*V1*V1+ V1*V2*V1*V1+ V1*V1*V2*V1+ V1*V1*V1*V2
      ROT(5,2)  =V3*V1*V1*V1+ V1*V3*V1*V1+ V1*V1*V3*V1+ V1*V1*V1*V3
      ROT(6,2)  =V1*V2*V2*V2+ V2*V1*V2*V2+ V2*V2*V1*V2+ V2*V2*V2*V1
      ROT(7,2)  =V3*V2*V2*V2+ V2*V3*V2*V2+ V2*V2*V3*V2+ V2*V2*V2*V3
      ROT(8,2)  =V1*V3*V3*V3+ V3*V1*V3*V3+ V3*V3*V1*V3+ V3*V3*V3*V1
      ROT(9,2)  =V2*V3*V3*V3+ V3*V2*V3*V3+ V3*V3*V2*V3+ V3*V3*V3*V2
      ROT(10,2) =V1*V1*V2*V2+ V2*V2*V1*V1+ V1*V2*V1*V2+ V1*V2*V2*V1
     1           + V2*V1*V1*V2+ V2*V1*V2*V1
      ROT(11,2) =V1*V1*V3*V3+ V3*V3*V1*V1+ V1*V3*V1*V3+ V1*V3*V3*V1
     1           + V3*V1*V1*V3+ V3*V1*V3*V1
      ROT(12,2) =V2*V2*V3*V3+ V3*V3*V2*V2+ V2*V3*V2*V3+ V2*V3*V3*V2
     1           + V3*V2*V2*V3+ V3*V2*V3*V2
      ROT(13,2) =V1*V1*V2*V3+ V1*V1*V3*V2+ V1*V2*V1*V3+ V1*V3*V1*V2
     1           + V1*V2*V3*V1+ V1*V3*V2*V1+ V2*V1*V1*V3+ V3*V1*V1*V2
     2           + V2*V1*V3*V1+ V3*V1*V2*V1+ V2*V3*V1*V1+ V3*V2*V1*V1
      ROT(14,2) =V2*V2*V1*V3+ V2*V2*V3*V1+ V2*V1*V2*V3+ V2*V3*V2*V1
     1           + V2*V1*V3*V2+ V2*V3*V1*V2+ V1*V2*V2*V3+ V3*V2*V2*V1
     2           + V1*V2*V3*V2+ V3*V2*V1*V2+ V1*V3*V2*V2+ V3*V1*V2*V2
      ROT(15,2) =V3*V3*V1*V2+ V3*V3*V2*V1+ V3*V1*V3*V2+ V3*V2*V3*V1
     1           + V3*V1*V2*V3+ V3*V2*V1*V3+ V1*V3*V3*V2+ V2*V3*V3*V1
     2           + V1*V3*V2*V3+ V2*V3*V1*V3+ V1*V2*V3*V3+ V2*V1*V3*V3
C     NEW ZZZZ
      ROT(1,3)  =W1*W1*W1*W1
      ROT(2,3)  =W2*W2*W2*W2
      ROT(3,3)  =W3*W3*W3*W3
      ROT(4,3)  =W2*W1*W1*W1+ W1*W2*W1*W1+ W1*W1*W2*W1+ W1*W1*W1*W2
      ROT(5,3)  =W3*W1*W1*W1+ W1*W3*W1*W1+ W1*W1*W3*W1+ W1*W1*W1*W3
      ROT(6,3)  =W1*W2*W2*W2+ W2*W1*W2*W2+ W2*W2*W1*W2+ W2*W2*W2*W1
      ROT(7,3)  =W3*W2*W2*W2+ W2*W3*W2*W2+ W2*W2*W3*W2+ W2*W2*W2*W3
      ROT(8,3)  =W1*W3*W3*W3+ W3*W1*W3*W3+ W3*W3*W1*W3+ W3*W3*W3*W1
      ROT(9,3)  =W2*W3*W3*W3+ W3*W2*W3*W3+ W3*W3*W2*W3+ W3*W3*W3*W2
      ROT(10,3) =W1*W1*W2*W2+ W2*W2*W1*W1+ W1*W2*W1*W2+ W1*W2*W2*W1
     1           + W2*W1*W1*W2+ W2*W1*W2*W1
      ROT(11,3) =W1*W1*W3*W3+ W3*W3*W1*W1+ W1*W3*W1*W3+ W1*W3*W3*W1
     1           + W3*W1*W1*W3+ W3*W1*W3*W1
      ROT(12,3) =W2*W2*W3*W3+ W3*W3*W2*W2+ W2*W3*W2*W3+ W2*W3*W3*W2
     1           + W3*W2*W2*W3+ W3*W2*W3*W2
      ROT(13,3) =W1*W1*W2*W3+ W1*W1*W3*W2+ W1*W2*W1*W3+ W1*W3*W1*W2
     1           + W1*W2*W3*W1+ W1*W3*W2*W1+ W2*W1*W1*W3+ W3*W1*W1*W2
     2           + W2*W1*W3*W1+ W3*W1*W2*W1+ W2*W3*W1*W1+ W3*W2*W1*W1
      ROT(14,3) =W2*W2*W1*W3+ W2*W2*W3*W1+ W2*W1*W2*W3+ W2*W3*W2*W1
     1           + W2*W1*W3*W2+ W2*W3*W1*W2+ W1*W2*W2*W3+ W3*W2*W2*W1
     2           + W1*W2*W3*W2+ W3*W2*W1*W2+ W1*W3*W2*W2+ W3*W1*W2*W2
      ROT(15,3) =W3*W3*W1*W2+ W3*W3*W2*W1+ W3*W1*W3*W2+ W3*W2*W3*W1
     1           + W3*W1*W2*W3+ W3*W2*W1*W3+ W1*W3*W3*W2+ W2*W3*W3*W1
     2           + W1*W3*W2*W3+ W2*W3*W1*W3+ W1*W2*W3*W3+ W2*W1*W3*W3
C     NEW XXXY
      ROT(1,4)  =U1*U1*U1*V1
      ROT(2,4)  =U2*U2*U2*V2
      ROT(3,4)  =U3*U3*U3*V3
      ROT(4,4)  =U2*U1*U1*V1+ U1*U2*U1*V1+ U1*U1*U2*V1+ U1*U1*U1*V2
      ROT(5,4)  =U3*U1*U1*V1+ U1*U3*U1*V1+ U1*U1*U3*V1+ U1*U1*U1*V3
      ROT(6,4)  =U1*U2*U2*V2+ U2*U1*U2*V2+ U2*U2*U1*V2+ U2*U2*U2*V1
      ROT(7,4)  =U3*U2*U2*V2+ U2*U3*U2*V2+ U2*U2*U3*V2+ U2*U2*U2*V3
      ROT(8,4)  =U1*U3*U3*V3+ U3*U1*U3*V3+ U3*U3*U1*V3+ U3*U3*U3*V1
      ROT(9,4)  =U2*U3*U3*V3+ U3*U2*U3*V3+ U3*U3*U2*V3+ U3*U3*U3*V2
      ROT(10,4) =U1*U1*U2*V2+ U2*U2*U1*V1+ U1*U2*U1*V2+ U1*U2*U2*V1
     1           + U2*U1*U1*V2+ U2*U1*U2*V1
      ROT(11,4) =U1*U1*U3*V3+ U3*U3*U1*V1+ U1*U3*U1*V3+ U1*U3*U3*V1
     1           + U3*U1*U1*V3+ U3*U1*U3*V1
      ROT(12,4) =U2*U2*U3*V3+ U3*U3*U2*V2+ U2*U3*U2*V3+ U2*U3*U3*V2
     1           + U3*U2*U2*V3+ U3*U2*U3*V2
      ROT(13,4) =U1*U1*U2*V3+ U1*U1*U3*V2+ U1*U2*U1*V3+ U1*U3*U1*V2
     1           + U1*U2*U3*V1+ U1*U3*U2*V1+ U2*U1*U1*V3+ U3*U1*U1*V2
     2           + U2*U1*U3*V1+ U3*U1*U2*V1+ U2*U3*U1*V1+ U3*U2*U1*V1
      ROT(14,4) =U2*U2*U1*V3+ U2*U2*U3*V1+ U2*U1*U2*V3+ U2*U3*U2*V1
     1           + U2*U1*U3*V2+ U2*U3*U1*V2+ U1*U2*U2*V3+ U3*U2*U2*V1
     2           + U1*U2*U3*V2+ U3*U2*U1*V2+ U1*U3*U2*V2+ U3*U1*U2*V2
      ROT(15,4) =U3*U3*U1*V2+ U3*U3*U2*V1+ U3*U1*U3*V2+ U3*U2*U3*V1
     1           + U3*U1*U2*V3+ U3*U2*U1*V3+ U1*U3*U3*V2+ U2*U3*U3*V1
     2           + U1*U3*U2*V3+ U2*U3*U1*V3+ U1*U2*U3*V3+ U2*U1*U3*V3
C     NEW XXXZ
      ROT(1,5)  =U1*U1*U1*W1
      ROT(2,5)  =U2*U2*U2*W2
      ROT(3,5)  =U3*U3*U3*W3
      ROT(4,5)  =U2*U1*U1*W1+ U1*U2*U1*W1+ U1*U1*U2*W1+ U1*U1*U1*W2
      ROT(5,5)  =U3*U1*U1*W1+ U1*U3*U1*W1+ U1*U1*U3*W1+ U1*U1*U1*W3
      ROT(6,5)  =U1*U2*U2*W2+ U2*U1*U2*W2+ U2*U2*U1*W2+ U2*U2*U2*W1
      ROT(7,5)  =U3*U2*U2*W2+ U2*U3*U2*W2+ U2*U2*U3*W2+ U2*U2*U2*W3
      ROT(8,5)  =U1*U3*U3*W3+ U3*U1*U3*W3+ U3*U3*U1*W3+ U3*U3*U3*W1
      ROT(9,5)  =U2*U3*U3*W3+ U3*U2*U3*W3+ U3*U3*U2*W3+ U3*U3*U3*W2
      ROT(10,5) =U1*U1*U2*W2+ U2*U2*U1*W1+ U1*U2*U1*W2+ U1*U2*U2*W1
     1           + U2*U1*U1*W2+ U2*U1*U2*W1
      ROT(11,5) =U1*U1*U3*W3+ U3*U3*U1*W1+ U1*U3*U1*W3+ U1*U3*U3*W1
     1           + U3*U1*U1*W3+ U3*U1*U3*W1
      ROT(12,5) =U2*U2*U3*W3+ U3*U3*U2*W2+ U2*U3*U2*W3+ U2*U3*U3*W2
     1           + U3*U2*U2*W3+ U3*U2*U3*W2
      ROT(13,5) =U1*U1*U2*W3+ U1*U1*U3*W2+ U1*U2*U1*W3+ U1*U3*U1*W2
     1           + U1*U2*U3*W1+ U1*U3*U2*W1+ U2*U1*U1*W3+ U3*U1*U1*W2
     2           + U2*U1*U3*W1+ U3*U1*U2*W1+ U2*U3*U1*W1+ U3*U2*U1*W1
      ROT(14,5) =U2*U2*U1*W3+ U2*U2*U3*W1+ U2*U1*U2*W3+ U2*U3*U2*W1
     1           + U2*U1*U3*W2+ U2*U3*U1*W2+ U1*U2*U2*W3+ U3*U2*U2*W1
     2           + U1*U2*U3*W2+ U3*U2*U1*W2+ U1*U3*U2*W2+ U3*U1*U2*W2
      ROT(15,5) =U3*U3*U1*W2+ U3*U3*U2*W1+ U3*U1*U3*W2+ U3*U2*U3*W1
     1           + U3*U1*U2*W3+ U3*U2*U1*W3+ U1*U3*U3*W2+ U2*U3*U3*W1
     2           + U1*U3*U2*W3+ U2*U3*U1*W3+ U1*U2*U3*W3+ U2*U1*U3*W3
C     NEW YYYX
      ROT(1,6)  =V1*V1*V1*U1
      ROT(2,6)  =V2*V2*V2*U2
      ROT(3,6)  =V3*V3*V3*U3
      ROT(4,6)  =V2*V1*V1*U1+ V1*V2*V1*U1+ V1*V1*V2*U1+ V1*V1*V1*U2
      ROT(5,6)  =V3*V1*V1*U1+ V1*V3*V1*U1+ V1*V1*V3*U1+ V1*V1*V1*U3
      ROT(6,6)  =V1*V2*V2*U2+ V2*V1*V2*U2+ V2*V2*V1*U2+ V2*V2*V2*U1
      ROT(7,6)  =V3*V2*V2*U2+ V2*V3*V2*U2+ V2*V2*V3*U2+ V2*V2*V2*U3
      ROT(8,6)  =V1*V3*V3*U3+ V3*V1*V3*U3+ V3*V3*V1*U3+ V3*V3*V3*U1
      ROT(9,6)  =V2*V3*V3*U3+ V3*V2*V3*U3+ V3*V3*V2*U3+ V3*V3*V3*U2
      ROT(10,6) =V1*V1*V2*U2+ V2*V2*V1*U1+ V1*V2*V1*U2+ V1*V2*V2*U1
     1           + V2*V1*V1*U2+ V2*V1*V2*U1
      ROT(11,6) =V1*V1*V3*U3+ V3*V3*V1*U1+ V1*V3*V1*U3+ V1*V3*V3*U1
     1           + V3*V1*V1*U3+ V3*V1*V3*U1
      ROT(12,6) =V2*V2*V3*U3+ V3*V3*V2*U2+ V2*V3*V2*U3+ V2*V3*V3*U2
     1           + V3*V2*V2*U3+ V3*V2*V3*U2
      ROT(13,6) =V1*V1*V2*U3+ V1*V1*V3*U2+ V1*V2*V1*U3+ V1*V3*V1*U2
     1           + V1*V2*V3*U1+ V1*V3*V2*U1+ V2*V1*V1*U3+ V3*V1*V1*U2
     2           + V2*V1*V3*U1+ V3*V1*V2*U1+ V2*V3*V1*U1+ V3*V2*V1*U1
      ROT(14,6) =V2*V2*V1*U3+ V2*V2*V3*U1+ V2*V1*V2*U3+ V2*V3*V2*U1
     1           + V2*V1*V3*U2+ V2*V3*V1*U1+ V1*V2*V2*U3+ V3*V2*V2*U1
     2           + V1*V2*V3*U2+ V3*V2*V1*U2+ V1*V3*V2*U2+ V3*V1*V2*U2
      ROT(15,6) =V3*V3*V1*U2+ V3*V3*V2*U1+ V3*V1*V3*U2+ V3*V2*V3*U1
     1           + V3*V1*V2*U3+ V3*V2*V1*U3+ V1*V3*V3*U2+ V2*V3*V3*U1
     2           + V1*V3*V2*U3+ V2*V3*V1*U3+ V1*V2*V3*U3+ V2*V1*V3*U3
C     NEW YYYZ
      ROT(1,7)  =V1*V1*V1*W1
      ROT(2,7)  =V2*V2*V2*W2
      ROT(3,7)  =V3*V3*V3*W3
      ROT(4,7)  =V2*V1*V1*W1+ V1*V2*V1*W1+ V1*V1*V2*W1+ V1*V1*V1*W2
      ROT(5,7)  =V3*V1*V1*W1+ V1*V3*V1*W1+ V1*V1*V3*W1+ V1*V1*V1*W3
      ROT(6,7)  =V1*V2*V2*W2+ V2*V1*V2*W2+ V2*V2*V1*W2+ V2*V2*V2*W1
      ROT(7,7)  =V3*V2*V2*W2+ V2*V3*V2*W2+ V2*V2*V3*W2+ V2*V2*V2*W3
      ROT(8,7)  =V1*V3*V3*W3+ V3*V1*V3*W3+ V3*V3*V1*W3+ V3*V3*V3*W1
      ROT(9,7)  =V2*V3*V3*W3+ V3*V2*V3*W3+ V3*V3*V2*W3+ V3*V3*V3*W2
      ROT(10,7) =V1*V1*V2*W2+ V2*V2*V1*W1+ V1*V2*V1*W2+ V1*V2*V2*W1
     1           + V2*V1*V1*W2+ V2*V1*V2*W1
      ROT(11,7) =V1*V1*V3*W3+ V3*V3*V1*W1+ V1*V3*V1*W3+ V1*V3*V3*W1
     1           + V3*V1*V1*W3+ V3*V1*V3*W1
      ROT(12,7) =V2*V2*V3*W3+ V3*V3*V2*W2+ V2*V3*V2*W3+ V2*V3*V3*W2
     1           + V3*V2*V2*W3+ V3*V2*V3*W2
      ROT(13,7) =V1*V1*V2*W3+ V1*V1*V3*W2+ V1*V2*V1*W3+ V1*V3*V1*W2
     1           + V1*V2*V3*W1+ V1*V3*V2*W1+ V2*V1*V1*W3+ V3*V1*V1*W2
     2           + V2*V1*V3*W1+ V3*V1*V2*W1+ V2*V3*V1*W1+ V3*V2*V1*W1
      ROT(14,7) =V2*V2*V1*W3+ V2*V2*V3*W1+ V2*V1*V2*W3+ V2*V3*V2*W1
     1           + V2*V1*V3*W2+ V2*V3*V1*W1+ V1*V2*V2*W3+ V3*V2*V2*W1
     2           + V1*V2*V3*W2+ V3*V2*V1*W2+ V1*V3*V2*W2+ V3*V1*V2*W2
      ROT(15,7) =V3*V3*V1*W2+ V3*V3*V2*W1+ V3*V1*V3*W2+ V3*V2*V3*W1
     1           + V3*V1*V2*W3+ V3*V2*V1*W3+ V1*V3*V3*W2+ V2*V3*V3*W1
     2           + V1*V3*V2*W3+ V2*V3*V1*W3+ V1*V2*V3*W3+ V2*V1*V3*W3
C     NEW ZZZX
      ROT(1,8)  =W1*W1*W1*U1
      ROT(2,8)  =W2*W2*W2*U2
      ROT(3,8)  =W3*W3*W3*U3
      ROT(4,8)  =W2*W1*W1*U1+ W1*W2*W1*U1+ W1*W1*W2*U1+ W1*W1*W1*U2
      ROT(5,8)  =W3*W1*W1*U1+ W1*W3*W1*U1+ W1*W1*W3*U1+ W1*W1*W1*U3
      ROT(6,8)  =W1*W2*W2*U2+ W2*W1*W2*U2+ W2*W2*W1*U2+ W2*W2*W2*U1
      ROT(7,8)  =W3*W2*W2*U2+ W2*W3*W2*U2+ W2*W2*W3*U2+ W2*W2*W2*U3
      ROT(8,8)  =W1*W3*W3*U3+ W3*W1*W3*U3+ W3*W3*W1*U3+ W3*W3*W3*U1
      ROT(9,8)  =W2*W3*W3*U3+ W3*W2*W3*U3+ W3*W3*W2*U3+ W3*W3*W3*U2
      ROT(10,8) =W1*W1*W2*U2+ W2*W2*W1*U1+ W1*W2*W1*U2+ W1*W2*W2*U1
     1           + W2*W1*W1*U2+ W2*W1*W2*U1
      ROT(11,8) =W1*W1*W3*U3+ W3*W3*W1*U1+ W1*W3*W1*U3+ W1*W3*W3*U1
     1           + W3*W1*W1*U3+ W3*W1*W3*U1
      ROT(12,8) =W2*W2*W3*U3+ W3*W3*W2*U2+ W2*W3*W2*U3+ W2*W3*W3*U2
     1           + W3*W2*W2*U3+ W3*W2*W3*U2
      ROT(13,8) =W1*W1*W2*U3+ W1*W1*W3*U2+ W1*W2*W1*U3+ W1*W3*W1*U2
     1           + W1*W2*W3*U1+ W1*W3*W2*U1+ W2*W1*W1*U3+ W3*W1*W1*U2
     2           + W2*W1*W3*U1+ W3*W1*W2*U1+ W2*W3*W1*U1+ W3*W2*W1*U1
      ROT(14,8) =W2*W2*W1*U3+ W2*W2*W3*U1+ W2*W1*W2*U3+ W2*W3*W2*U1
     1           + W2*W1*W3*U2+ W2*W3*W1*U2+ W1*W2*W2*U3+ W3*W2*W2*U1
     2           + W1*W2*W3*U2+ W3*W2*W1*U2+ W1*W3*W2*U2+ W3*W1*W2*U2
      ROT(15,8) =W3*W3*W1*U2+ W3*W3*W2*U1+ W3*W1*W3*U2+ W3*W2*W3*U1
     1           + W3*W1*W2*U3+ W3*W2*W1*U3+ W1*W3*W3*U2+ W2*W3*W3*U1
     2           + W1*W3*W2*U3+ W2*W3*W1*U3+ W1*W2*W3*U3+ W2*W1*W3*U3
C     NEW ZZZY
      ROT(1,9)  =W1*W1*W1*V1
      ROT(2,9)  =W2*W2*W2*V2
      ROT(3,9)  =W3*W3*W3*V3
      ROT(4,9)  =W2*W1*W1*V1+ W1*W2*W1*V1+ W1*W1*W2*V1+ W1*W1*W1*V2
      ROT(5,9)  =W3*W1*W1*V1+ W1*W3*W1*V1+ W1*W1*W3*V1+ W1*W1*W1*V3
      ROT(6,9)  =W1*W2*W2*V2+ W2*W1*W2*V2+ W2*W2*W1*V2+ W2*W2*W2*V1
      ROT(7,9)  =W3*W2*W2*V2+ W2*W3*W2*V2+ W2*W2*W3*V2+ W2*W2*W2*V3
      ROT(8,9)  =W1*W3*W3*V3+ W3*W1*W3*V3+ W3*W3*W1*V3+ W3*W3*W3*V1
      ROT(9,9)  =W2*W3*W3*V3+ W3*W2*W3*V3+ W3*W3*W2*V3+ W3*W3*W3*V2
      ROT(10,9) =W1*W1*W2*V2+ W2*W2*W1*V1+ W1*W2*W1*V2+ W1*W2*W2*V1
     1           + W2*W1*W1*V2+ W2*W1*W2*V1
      ROT(11,9) =W1*W1*W3*V3+ W3*W3*W1*V1+ W1*W3*W1*V3+ W1*W3*W3*V1
     1           + W3*W1*W1*V3+ W3*W1*W3*V1
      ROT(12,9) =W2*W2*W3*V3+ W3*W3*W2*V2+ W2*W3*W2*V3+ W2*W3*W3*V2
     1           + W3*W2*W2*V3+ W3*W2*W3*V2
      ROT(13,9) =W1*W1*W2*V3+ W1*W1*W3*V2+ W1*W2*W1*V3+ W1*W3*W1*V2
     1           + W1*W2*W3*V1+ W1*W3*W2*V1+ W2*W1*W1*V3+ W3*W1*W1*V2
     2           + W2*W1*W3*V1+ W3*W1*W2*V1+ W2*W3*W1*V1+ W3*W2*W1*V1
      ROT(14,9) =W2*W2*W1*V3+ W2*W2*W3*V1+ W2*W1*W2*V3+ W2*W3*W2*V1
     1           + W2*W1*W3*V2+ W2*W3*W1*V2+ W1*W2*W2*V3+ W3*W2*W2*V1
     2           + W1*W2*W3*V2+ W3*W2*W1*V2+ W1*W3*W2*V2+ W3*W1*W2*V2
      ROT(15,9) =W3*W3*W1*V2+ W3*W3*W2*V1+ W3*W1*W3*V2+ W3*W2*W3*V1
     1           + W3*W1*W2*V3+ W3*W2*W1*V3+ W1*W3*W3*V2+ W2*W3*W3*V1
     2           + W1*W3*W2*V3+ W2*W3*W1*V3+ W1*W2*W3*V3+ W2*W1*W3*V3
C     NEW XXYY
      ROT(1,10) =U1*U1*V1*V1
      ROT(2,10) =U2*U2*V2*V2
      ROT(3,10) =U3*U3*V3*V3
      ROT(4,10) =U2*U1*V1*V1+ U1*U2*V1*V1+ U1*U1*V2*V1+ U1*U1*V1*V2
      ROT(5,10) =U3*U1*V1*V1+ U1*U3*V1*V1+ U1*U1*V3*V1+ U1*U1*V1*V3
      ROT(6,10) =U1*U2*V2*V2+ U2*U1*V2*V2+ U2*U2*V1*V2+ U2*U2*V2*V1
      ROT(7,10) =U3*U2*V2*V2+ U2*U3*V2*V2+ U2*U2*V3*V2+ U2*U2*V2*V3
      ROT(8,10) =U1*U3*V3*V3+ U3*U1*V3*V3+ U3*U3*V1*V3+ U3*U3*V3*V1
      ROT(9,10) =U2*U3*V3*V3+ U3*U2*V3*V3+ U3*U3*V2*V3+ U3*U3*V3*V2
      ROT(10,10)=U1*U1*V2*V2+ U2*U2*V1*V1+ U1*U2*V1*V2+ U1*U2*V2*V1
     1           + U2*U1*V1*V2+ U2*U1*V2*V1
      ROT(11,10)=U1*U1*V3*V3+ U3*U3*V1*V1+ U1*U3*V1*V3+ U1*U3*V3*V1
     1           + U3*U1*V1*V3+ U3*U1*V3*V1
      ROT(12,10)=U2*U2*V3*V3+ U3*U3*V2*V2+ U2*U3*V2*V3+ U2*U3*V3*V2
     1           + U3*U2*V2*V3+ U3*U2*V3*V2
      ROT(13,10)=U1*U1*V2*V3+ U1*U1*V3*V2+ U1*U2*V1*V3+ U1*U3*V1*V2
     1           + U1*U2*V3*V1+ U1*U3*V2*V1+ U2*U1*V1*V3+ U3*U1*V1*V2
     2           + U2*U1*V3*V1+ U3*U1*V2*V1+ U2*U3*V1*V1+ U3*U2*V1*V1
      ROT(14,10)=U2*U2*V1*V3+ U2*U2*V3*V1+ U2*U1*V2*V3+ U2*U3*V2*V1
     1           + U2*U1*V3*V2+ U2*U3*V1*V2+ U1*U2*V2*V3+ U3*U2*V2*V1
     2           + U1*U2*V3*V2+ U3*U2*V1*V2+ U1*U3*V2*V2+ U3*U1*V2*V2
      ROT(15,10)=U3*U3*V1*V2+ U3*U3*V2*V1+ U3*U1*V3*V2+ U3*U2*V3*V1
     1           + U3*U1*V2*V3+ U3*U2*V1*V3+ U1*U3*V3*V2+ U2*U3*V3*V1
     2           + U1*U3*V2*V3+ U2*U3*V1*V3+ U1*U2*V3*V3+ U2*U1*V3*V3
C     NEW XXZZ
      ROT(1,11) =U1*U1*W1*W1
      ROT(2,11) =U2*U2*W2*W2
      ROT(3,11) =U3*U3*W3*W3
      ROT(4,11) =U2*U1*W1*W1+ U1*U2*W1*W1+ U1*U1*W2*W1+ U1*U1*W1*W2
      ROT(5,11) =U3*U1*W1*W1+ U1*U3*W1*W1+ U1*U1*W3*W1+ U1*U1*W1*W3
      ROT(6,11) =U1*U2*W2*W2+ U2*U1*W2*W2+ U2*U2*W1*W2+ U2*U2*W2*W1
      ROT(7,11) =U3*U2*W2*W2+ U2*U3*W2*W2+ U2*U2*W3*W2+ U2*U2*W2*W3
      ROT(8,11) =U1*U3*W3*W3+ U3*U1*W3*W3+ U3*U3*W1*W3+ U3*U3*W3*W1
      ROT(9,11) =U2*U3*W3*W3+ U3*U2*W3*W3+ U3*U3*W2*W3+ U3*U3*W3*W2
      ROT(10,11)=U1*U1*W2*W2+ U2*U2*W1*W1+ U1*U2*W1*W2+ U1*U2*W2*W1
     1           + U2*U1*W1*W2+ U2*U1*W2*W1
      ROT(11,11)=U1*U1*W3*W3+ U3*U3*W1*W1+ U1*U3*W1*W3+ U1*U3*W3*W1
     1           + U3*U1*W1*W3+ U3*U1*W3*W1
      ROT(12,11)=U2*U2*W3*W3+ U3*U3*W2*W2+ U2*U3*W2*W3+ U2*U3*W3*W2
     1           + U3*U2*W2*W3+ U3*U2*W3*W2
      ROT(13,11)=U1*U1*W2*W3+ U1*U1*W3*W2+ U1*U2*W1*W3+ U1*U3*W1*W2
     1           + U1*U2*W3*W1+ U1*U3*W2*W1+ U2*U1*W1*W3+ U3*U1*W1*W2
     2           + U2*U1*W3*W1+ U3*U1*W2*W1+ U2*U3*W1*W1+ U3*U2*W1*W1
      ROT(14,11)=U2*U2*W1*W3+ U2*U2*W3*W1+ U2*U1*W2*W3+ U2*U3*W2*W1
     1           + U2*U1*W3*W2+ U2*U3*W1*W2+ U1*U2*W2*W3+ U3*U2*W2*W1
     2           + U1*U2*W3*W2+ U3*U2*W1*W2+ U1*U3*W2*W2+ U3*U1*W2*W2
      ROT(15,11)=U3*U3*W1*W2+ U3*U3*W2*W1+ U3*U1*W3*W2+ U3*U2*W3*W1
     1           + U3*U1*W2*W3+ U3*U2*W1*W3+ U1*U3*W3*W2+ U2*U3*W3*W1
     2           + U1*U3*W2*W3+ U2*U3*W1*W3+ U1*U2*W3*W3+ U2*U1*W3*W3
C     NEW YYZZ
      ROT(1,12) =V1*V1*W1*W1
      ROT(2,12) =V2*V2*W2*W2
      ROT(3,12) =V3*V3*W3*W3
      ROT(4,12) =V2*V1*W1*W1+ V1*V2*W1*W1+ V1*V1*W2*W1+ V1*V1*W1*W2
      ROT(5,12) =V3*V1*W1*W1+ V1*V3*W1*W1+ V1*V1*W3*W1+ V1*V1*W1*W3
      ROT(6,12) =V1*V2*W2*W2+ V2*V1*W2*W2+ V2*V2*W1*W2+ V2*V2*W2*W1
      ROT(7,12) =V3*V2*W2*W2+ V2*V3*W2*W2+ V2*V2*W3*W2+ V2*V2*W2*W3
      ROT(8,12) =V1*V3*W3*W3+ V3*V1*W3*W3+ V3*V3*W1*W3+ V3*V3*W3*W1
      ROT(9,12) =V2*V3*W3*W3+ V3*V2*W3*W3+ V3*V3*W2*W3+ V3*V3*W3*W2
      ROT(10,12)=V1*V1*W2*W2+ V2*V2*W1*W1+ V1*V2*W1*W2+ V1*V2*W2*W1
     1           + V2*V1*W1*W2+ V2*V1*W2*W1
      ROT(11,12)=V1*V1*W3*W3+ V3*V3*W1*W1+ V1*V3*W1*W3+ V1*V3*W3*W1
     1           + V3*V1*W1*W3+ V3*V1*W3*W1
      ROT(12,12)=V2*V2*W3*W3+ V3*V3*W2*W2+ V2*V3*W2*W3+ V2*V3*W3*W2
     1           + V3*V2*W2*W3+ V3*V2*W3*W2
      ROT(13,12)=V1*V1*W2*W3+ V1*V1*W3*W2+ V1*V2*W1*W3+ V1*V3*W1*W2
     1           + V1*V2*W3*W1+ V1*V3*W2*W1+ V2*V1*W1*W3+ V3*V1*W1*W2
     2           + V2*V1*W3*W1+ V3*V1*W2*W1+ V2*V3*W1*W1+ V3*V2*W1*W1
      ROT(14,12)=V2*V2*W1*W3+ V2*V2*W3*W1+ V2*V1*W2*W3+ V2*V3*W2*W1
     1           + V2*V1*W3*W2+ V2*V3*W1*W2+ V1*V2*W2*W3+ V3*V2*W2*W1
     2           + V1*V2*W3*W2+ V3*V2*W1*W2+ V1*V3*W2*W2+ V3*V1*W2*W2
      ROT(15,12)=V3*V3*W1*W2+ V3*V3*W2*W1+ V3*V1*W3*W2+ V3*V2*W3*W1
     1           + V3*V1*W2*W3+ V3*V2*W1*W3+ V1*V3*W3*W2+ V2*V3*W3*W1
     2           + V1*V3*W2*W3+ V2*V3*W1*W3+ V1*V2*W3*W3+ V2*V1*W3*W3
C     NEW XXYZ
      ROT(1,13) =U1*U1*V1*W1
      ROT(2,13) =U2*U2*V2*W2
      ROT(3,13) =U3*U3*V3*W3
      ROT(4,13) =U2*U1*V1*W1+ U1*U2*V1*W1+ U1*U1*V2*W1+ U1*U1*V1*W2
      ROT(5,13) =U3*U1*V1*W1+ U1*U3*V1*W1+ U1*U1*V3*W1+ U1*U1*V1*W3
      ROT(6,13) =U1*U2*V2*W2+ U2*U1*V2*W2+ U2*U2*V1*W2+ U2*U2*V2*W1
      ROT(7,13) =U3*U2*V2*W2+ U2*U3*V2*W2+ U2*U2*V3*W2+ U2*U2*V2*W3
      ROT(8,13) =U1*U3*V3*W3+ U3*U1*V3*W3+ U3*U3*V1*W3+ U3*U3*V3*W1
      ROT(9,13) =U2*U3*V3*W3+ U3*U2*V3*W3+ U3*U3*V2*W3+ U3*U3*V3*W2
      ROT(10,13)=U1*U1*V2*W2+ U2*U2*V1*W1+ U1*U2*V1*W2+ U1*U2*V2*W1
     1           + U2*U1*V1*W2+ U2*U1*V2*W1
      ROT(11,13)=U1*U1*V3*W3+ U3*U3*V1*W1+ U1*U3*V1*W3+ U1*U3*V3*W1
     1           + U3*U1*V1*W3+ U3*U1*V3*W1
      ROT(12,13)=U2*U2*V3*W3+ U3*U3*V2*W2+ U2*U3*V2*W3+ U2*U3*V3*W2
     1           + U3*U2*V2*W3+ U3*U2*V3*W2
      ROT(13,13)=U1*U1*V2*W3+ U1*U1*V3*W2+ U1*U2*V1*W3+ U1*U3*V1*W2
     1           + U1*U2*V3*W1+ U1*U3*V2*W1+ U2*U1*V1*W3+ U3*U1*V1*W2
     2           + U2*U1*V3*W1+ U3*U1*V2*W1+ U2*U3*V1*W1+ U3*U2*V1*W1
      ROT(14,13)=U2*U2*V1*W3+ U2*U2*V3*W1+ U2*U1*V2*W3+ U2*U3*V2*W1
     1           + U2*U1*V3*W2+ U2*U3*V1*W2+ U1*U2*V2*W3+ U3*U2*V2*W1
     2           + U1*U2*V3*W2+ U3*U2*V1*W2+ U1*U3*V2*W2+ U3*U1*V2*W2
      ROT(15,13)=U3*U3*V1*W2+ U3*U3*V2*W1+ U3*U1*V3*W2+ U3*U2*V3*W1
     1           + U3*U1*V2*W3+ U3*U2*V1*W3+ U1*U3*V3*W2+ U2*U3*V3*W1
     2           + U1*U3*V2*W3+ U2*U3*V1*W3+ U1*U2*V3*W3+ U2*U1*V3*W3
C     NEW YYXZ
      ROT(1,14) =V1*V1*U1*W1
      ROT(2,14) =V2*V2*U2*W2
      ROT(3,14) =V3*V3*U3*W3
      ROT(4,14) =V2*V1*U1*W1+ V1*V2*U1*W1+ V1*V1*U2*W1+ V1*V1*U1*W2
      ROT(5,14) =V3*V1*U1*W1+ V1*V3*U1*W1+ V1*V1*U3*W1+ V1*V1*U1*W3
      ROT(6,14) =V1*V2*U2*W2+ V2*V1*U2*W2+ V2*V2*U1*W2+ V2*V2*U2*W1
      ROT(7,14) =V3*V2*U2*W2+ V2*V3*U2*W2+ V2*V2*U3*W2+ V2*V2*U2*W3
      ROT(8,14) =V1*V3*U3*W3+ V3*V1*U3*W3+ V3*V3*U1*W3+ V3*V3*U3*W1
      ROT(9,14) =V2*V3*U3*W3+ V3*V2*U3*W3+ V3*V3*U2*W3+ V3*V3*U3*W2
      ROT(10,14)=V1*V1*U2*W2+ V2*V2*U1*W1+ V1*V2*U1*W2+ V1*V2*U2*W1
     1           + V2*V1*U1*W2+ V2*V1*U2*W1
      ROT(11,14)=V1*V1*U3*W3+ V3*V3*U1*W1+ V1*V3*U1*W3+ V1*V3*U3*W1
     1           + V3*V1*U1*W3+ V3*V1*U3*W1
      ROT(12,14)=V2*V2*U3*W3+ V3*V3*U2*W2+ V2*V3*U2*W3+ V2*V3*U3*W2
     1           + V3*V2*U2*W3+ V3*V2*U3*W2
      ROT(13,14)=V1*V1*U2*W3+ V1*V1*U3*W2+ V1*V2*U1*W3+ V1*V3*U1*W2
     1           + V1*V2*U3*W1+ V1*V3*U2*W1+ V2*V1*U1*W3+ V3*V1*U1*W2
     2           + V2*V1*U3*W1+ V3*V1*U2*W1+ V2*V3*U1*W1+ V3*V2*U1*W1
      ROT(14,14)=V2*V2*U1*W3+ V2*V2*U3*W1+ V2*V1*U2*W3+ V2*V3*U2*W1
     1           + V2*V1*U3*W2+ V2*V3*U1*W2+ V1*V2*U2*W3+ V3*V2*U2*W1
     2           + V1*V2*U3*W2+ V3*V2*U1*W2+ V1*V3*U2*W2+ V3*V1*U2*W2
      ROT(15,14)=V3*V3*U1*W2+ V3*V3*U2*W1+ V3*V1*U3*W2+ V3*V2*U3*W1
     1           + V3*V1*U2*W3+ V3*V2*U1*W3+ V1*V3*U3*W2+ V2*V3*U3*W1
     2           + V1*V3*U2*W3+ V2*V3*U1*W3+ V1*V2*U3*W3+ V2*V1*U3*W3
C     NEW ZZXY
      ROT(1,15) =W1*W1*U1*V1
      ROT(2,15) =W2*W2*U2*V2
      ROT(3,15) =W3*W3*U3*V3
      ROT(4,15) =W2*W1*U1*V1+ W1*W2*U1*V1+ W1*W1*U2*V1+ W1*W1*U1*V2
      ROT(5,15) =W3*W1*U1*V1+ W1*W3*U1*V1+ W1*W1*U3*V1+ W1*W1*U1*V3
      ROT(6,15) =W1*W2*U2*V2+ W2*W1*U2*V2+ W2*W2*U1*V2+ W2*W2*U2*V1
      ROT(7,15) =W3*W2*U2*V2+ W2*W3*U2*V2+ W2*W2*U3*V2+ W2*W2*U2*V3
      ROT(8,15) =W1*W3*U3*V3+ W3*W1*U3*V3+ W3*W3*U1*V3+ W3*W3*U3*V1
      ROT(9,15) =W2*W3*U3*V3+ W3*W2*U3*V3+ W3*W3*U2*V3+ W3*W3*U3*V2
      ROT(10,15)=W1*W1*U2*V2+ W2*W2*U1*V1+ W1*W2*U1*V2+ W1*W2*U2*V1
     1           + W2*W1*U1*V2+ W2*W1*U2*V1
      ROT(11,15)=W1*W1*U3*V3+ W3*W3*U1*V1+ W1*W3*U1*V3+ W1*W3*U3*V1
     1           + W3*W1*U1*V3+ W3*W1*U3*V1
      ROT(12,15)=W2*W2*U3*V3+ W3*W3*U2*V2+ W2*W3*U2*V3+ W2*W3*U3*V2
     1           + W3*W2*U2*V3+ W3*W2*U3*V2
      ROT(13,15)=W1*W1*U2*V3+ W1*W1*U3*V2+ W1*W2*U1*V3+ W1*W3*U1*V2
     1           + W1*W2*U3*V1+ W1*W3*U2*V1+ W2*W1*U1*V3+ W3*W1*U1*V2
     2           + W2*W1*U3*V1+ W3*W1*U2*V1+ W2*W3*U1*V1+ W3*W2*U1*V1
      ROT(14,15)=W2*W2*U1*V3+ W2*W2*U3*V1+ W2*W1*U2*V3+ W2*W3*U2*V1
     1           + W2*W1*U3*V2+ W2*W3*U1*V2+ W1*W2*U2*V3+ W3*W2*U2*V1
     2           + W1*W2*U3*V2+ W3*W2*U1*V2+ W1*W3*U2*V2+ W3*W1*U2*V2
      ROT(15,15)=W3*W3*U1*V2+ W3*W3*U2*V1+ W3*W1*U3*V2+ W3*W2*U3*V1
     1           + W3*W1*U2*V3+ W3*W2*U1*V3+ W1*W3*U3*V2+ W2*W3*U3*V1
     2           + W1*W3*U2*V3+ W2*W3*U1*V3+ W1*W2*U3*V3+ W2*W1*U3*V3
C
C     ----- EXPRESS SPHERICAL HARMONICS IN ABSOLUTE FRAME -----
C
      DO 360 I=1,15
      DO 360 J=1,15
      A=ZERO
      DO 350 K=1,15
      A=A+ROT(I,K)*GSPH(K,J)
 350  CONTINUE
        BUFF(I,J)=A
 360  GSHELL(I,J)=A
C
C     ----- DEFINE INVERSE TRANSFORMATION -----
C
      CALL LUDCMP(BUFF,15,15,INDX,VWRK,D,IERR)
      DO 380 I=1,15
      CALL LUBKSB(BUFF,15,15,INDX,GIHELL(1,I))
  380 CONTINUE
C
C     PRINT
C
      IF(DBG) CALL SYMPRP(GSHELL,15,15,15)
      IF(DBG) CALL SYMPRP(GIHELL,15,15,15)
      RETURN
 9999 FORMAT(/,' DEFINITION OF SPHERICAL HARMONICS',/,
     1         ' ---------------------------------')
      END
C*MODULE SYMSLC  *DECK SPHTRF
      SUBROUTINE SPHTRF(DBG)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL DBG
C
      PARAMETER (MXSH=5000, MXATM=2000)
C
      DIMENSION V(15,15),T(15,15),TI(15,15),BUFF(15,15),
     *          PNAME(3),DNAME(6),FNAME(10),GNAME(15)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /SPHERI/ PSHELL(3,3),DSHELL(6,6),
     *                FSHELL(10,10),GSHELL(15,15),
     *                PIHELL(3,3),DIHELL(6,6),
     *                FIHELL(10,10),GIHELL(15,15)
      COMMON /SYMSPD/ PTR(3,144),DTR(6,288),FTR(10,480),GTR(15,720)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                TSYM(432),INVT(48),NT
C
      PARAMETER (ZERO=0.0D+00)
C
      DATA PNAME /4HP1  ,4HP2  ,4HP3  /
      DATA DNAME /4HD1  ,4HD2  ,4HD3  ,4HD4  ,4HD5  ,4HD6  /
      DATA FNAME /4HF1  ,4HF2  ,4HF3  ,4HF4  ,4HF5  ,4HF6  ,4HF7  ,
     1            4HF8  ,4HF9  ,4HF10 /
      DATA GNAME /4HG1  ,4HG2  ,4HG3  ,4HG4  ,4HG5  ,4HG6  ,4HG7  ,
     1            4HG8  ,4HG9  ,4HG10 ,4HG11 ,4HG12 ,4HG13 ,4HG14 ,
     2            4HG15 /
C
C     ---- THIS ROUTINE CALCULATES THE SYMMETRY TRANSFORMATION -----
C          MATRICES IN THE BASIS OF SPHERICAL HARMONICS.
C
      IF(DBG) WRITE(IW,9999)
C
C     ----- LOOP OVER SYMMETRY OPERATIONS -----
C
      DO 60 IT=1,NT
      NP= 3*(IT-1)
      ND= 6*(IT-1)
      NF=10*(IT-1)
      NG=15*(IT-1)
C
      DO 20 IPER=1,4
C
      N=3
      IF(IPER.EQ.2) N=6
      IF(IPER.EQ.3) N=10
      IF(IPER.EQ.4) N=15
C
      GO TO (21,22,23,24),IPER
C
 21   DO 31 I=1,N
      DO 31 J=1,N
      V(I,J)=PTR(I,NP+J)
      T(I,J)=PSHELL(I,J)
      TI(I,J)=PIHELL(I,J)
 31   CONTINUE
      GO TO 25
C
 22   DO 32 I=1,N
      DO 32 J=1,N
      V(I,J)=DTR(I,ND+J)
      T(I,J)=DSHELL(I,J)
      TI(I,J)=DIHELL(I,J)
 32   CONTINUE
      GO TO 25
C
 23   DO 33 I=1,N
      DO 33 J=1,N
      V(I,J)=FTR(I,NF+J)
      T(I,J)=FSHELL(I,J)
      TI(I,J)=FIHELL(I,J)
 33   CONTINUE
      GO TO 25
C
 24   DO 34 I=1,N
      DO 34 J=1,N
      V(I,J)=GTR(I,NG+J)
      T(I,J)=GSHELL(I,J)
      TI(I,J)=GIHELL(I,J)
 34   CONTINUE
      GO TO 25
C
C     ---- RIGHT-TRANSFORMATION ----
C
 25   DO 100 I=1,N
      DO 100 J=1,N
      A=ZERO
      DO 120 K=1,N
      A= A + V(I,K)*T(K,J)
 120  CONTINUE
 100  BUFF(I,J)=A
C
C     ---- LEFT-TRANSFORMATION ----
C
      DO 130 I=1,N
      DO 130 J=1,N
      A=ZERO
      DO 140 K=1,N
      A= A + TI(I,K)*BUFF(K,J)
 140  CONTINUE
 130  V(I,J)=A
C
C     ---- COPY BACK -----
C
      GO TO (41,42,43,44),IPER
C
 41   DO 51 I=1,N
      DO 51 J=1,N
      PTR(I,NP+J)=V(I,J)
 51   CONTINUE
      GO TO 20
C
 42   DO 52 I=1,N
      DO 52 J=1,N
      DTR(I,ND+J)=V(I,J)
 52   CONTINUE
      GO TO 20
C
 43   DO 53 I=1,N
      DO 53 J=1,N
      FTR(I,NF+J)=V(I,J)
 53   CONTINUE
      GO TO 20
C
 44   DO 54 I=1,N
      DO 54 J=1,N
      GTR(I,NG+J)=V(I,J)
 54   CONTINUE
      GO TO 20
C
C     DONE
C
 20   CONTINUE
 60   CONTINUE
C
C     ----- PRINT SECTION -----
C
      IF(.NOT.DBG) GO TO 300
C
      DO 255 IT=1,NT
      WRITE(IW,8860) IT
C
      NP=3*(IT-1)
      WRITE(IW,8862) (PNAME(J),J=1,3)
      WRITE(IW,8863)
      DO 250 I=1,3
  250 WRITE(IW,8861) PNAME(I),(PTR(I,NP+J),J=1,3)
      WRITE(IW,8863)
C
      ND=6*(IT-1)
      WRITE(IW,8862) (DNAME(J),J=1,6)
      WRITE(IW,8863)
      DO 251 I=1,6
  251 WRITE(IW,8861) DNAME(I),(DTR(I,ND+J),J=1,6)
      WRITE(IW,8863)
C
      NF=10*(IT-1)
      WRITE(IW,8862) (FNAME(J),J=1,10)
      WRITE(IW,8863)
      DO 252 I=1,10
  252 WRITE(IW,8861) FNAME(I),(FTR(I,NF+J),J=1,10)
      WRITE(IW,8863)
C
      NG=15*(IT-1)
      JMAX=0
  253 JMIN=JMAX+1
      JMAX=JMAX+10
      IF(JMAX.GT.15) JMAX=15
      WRITE(IW,8862) (GNAME(J),J=JMIN,JMAX)
      WRITE(IW,8863)
      DO 254 I=1,15
  254 WRITE(IW,8861) GNAME(I),(GTR(I,NG+J),J=JMIN,JMAX)
      WRITE(IW,8863)
      IF(JMAX.LT.15) GO TO 253
  255 CONTINUE
C
  300 CONTINUE
      RETURN
C
 9999 FORMAT(/,5X,' SYMMETRY TRANSFORMATION MATRICES IN SPHERICAL',
     1            ' HARMONICS BASIS',/,5X,1X,61(1H-))
 8863 FORMAT(/)
 8862 FORMAT(8X,10(3X,A4,3X))
 8861 FORMAT(2X,A4,2X,10F10.6)
 8860 FORMAT(/,21X,'TRANSFORMATION NUMBER',I4,/)
      END
C*MODULE SYMSLC  *DECK SPHPRT
      SUBROUTINE SPHPRT
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MXGTOT=20000)
      PARAMETER (MXSH=5000)
      CHARACTER*4 NAMXYZ
      CHARACTER*4 NAMSPH
      CHARACTER*4 NAMLAB
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SPHERI/ PSHELL(3,3),DSHELL(6,6),
     *                FSHELL(10,10),GSHELL(15,15),
     *                PIHELL(3,3),DIHELL(6,6),
     *                FIHELL(10,10),GIHELL(15,15)
      DIMENSION NF(5),NAMXYZ(35),NAMSPH(35),NAMLAB(35)
      DATA NF     /1,3,6,10,15/
      DATA NAMXYZ /'   S','   X','   Y','   Z','  XX','  YY','  ZZ',
     1             '  XY','  XZ','  YZ',' XXX',' YYY',' ZZZ',' XXY',
     2             ' XXZ',' YYX',' YYZ',' ZZX',' ZZY',' XYZ','XXXX',
     3             'YYYY','ZZZZ','XXXY','XXXZ','YYYX','YYYZ','ZZZX',
     4             'ZZZY','XXYY','XXZZ','YYZZ','XXYZ','YYXZ','ZZXY'/
      DATA NAMSPH /'  S ','P X ','P Y ','P Z ','D Z2','D X2','D XY',
     1             'D XZ','D YZ','D S ','  F1','  F2','  F3','  F4',
     2             '  F5','  F6','  F7','F P1','F P2','F P3','  G1',
     3             '  G2','  G3','  G4','  G5','  G6','  G7','  G8',
     4             '  G9','G S ','G D1','G D2','G D3','G D4','G D5'/
      DATA NAMLAB /'  S ','  P1','  P2','  P3','  D1','  D2','  D3',
     1             '  D4','  D5','  DS','  F1','  F2','  F3','  F4',
     2             '  F5','  F6','  F7','F P1','F P2','F P3','  G1',
     3             '  G2','  G3','  G4','  G5','  G6','  G7','  G8',
     4             '  G9','G S ','G D1','G D2','G D3','G D4','G D5'/
      DATA ONE    /1.0D+00/
C
C     ----  PRINT SPHERICAL HARMONICS DEFINITION -----
C
      WRITE(IW,9999)
C
      MAXTYP=0
      DO ISH=1,NSHELL
         IF(KTYPE(ISH).GT.MAXTYP) MAXTYP=KTYPE(ISH)
      ENDDO
C
      DO 100 ISPH=1,MAXTYP
         N=NF(ISPH)
         GO TO (10,20,30,40,50),ISPH
C
C     ---- S ----
C
  10     J=1
            WRITE(IW,9998) J,NAMSPH(J   ),J,NAMLAB(J   ),
     1                     ONE,NAMXYZ(1)
         GO TO 100
C
C     ---- P ----
C
  20     DO J=1,N
            MAX=0
  21        MIN=MAX+1
            MAX=MAX+6
            IF(MAX.GT.N) MAX=N
            IF(MIN.LT.6) THEN
               WRITE(IW,9998) J,NAMSPH(J+ 1),J,NAMLAB(J+ 1),
     1                        (PSHELL(I,J),NAMXYZ(I+ 1),I=MIN,MAX)
            ELSE
               WRITE(IW,9997) (PSHELL(I,J),NAMXYZ(I+ 1),I=MIN,MAX)
            END IF
            IF(MAX.LT.N) GO TO 21
         ENDDO
         GO TO 100
C
C     ---- D ----
C
  30     DO J=1,N
            MAX=0
  31        MIN=MAX+1
            MAX=MAX+6
            IF(MAX.GT.N) MAX=N
            IF(MIN.LT.6) THEN
               WRITE(IW,9998) J,NAMSPH(J+4),J,NAMLAB(J+4),
     1                        (DSHELL(I,J),NAMXYZ(I+ 4),I=MIN,MAX)
            ELSE
               WRITE(IW,9997) (DSHELL(I,J),NAMXYZ(I+ 4),I=MIN,MAX)
            END IF
            IF(MAX.LT.N) GO TO 31
         ENDDO
         GO TO 100
C
C     ---- F ----
C
  40     DO J=1,N
            MAX=0
  41        MIN=MAX+1
            MAX=MAX+6
            IF(MAX.GT.N) MAX=N
            IF(MIN.LT.6) THEN
               WRITE(IW,9998) J,NAMSPH(J+10),J,NAMLAB(J+10),
     1                        (FSHELL(I,J),NAMXYZ(I+10),I=MIN,MAX)
            ELSE
               WRITE(IW,9997) (FSHELL(I,J),NAMXYZ(I+10),I=MIN,MAX)
            END IF
            IF(MAX.LT.N) GO TO 41
         ENDDO
         GO TO 100
C
C     ---- G ----
C
  50     DO J=1,N
            MAX=0
  51        MIN=MAX+1
            MAX=MAX+6
            IF(MAX.GT.N) MAX=N
            IF(MIN.LT.6) THEN
               WRITE(IW,9998) J,NAMSPH(J+20),J,NAMLAB(J+20),
     1                        (GSHELL(I,J),NAMXYZ(I+20),I=MIN,MAX)
            ELSE
               WRITE(IW,9997) (GSHELL(I,J),NAMXYZ(I+20),I=MIN,MAX)
            END IF
            IF(MAX.LT.N) GO TO 51
         ENDDO
         GO TO 100
C
  100 CONTINUE
C
      RETURN
 9999 FORMAT(/,1X,'SPHERICAL HARMONICS',/,1X,19(1H-))
 9998 FORMAT(5H CHI(,I4,2X,A4,3H) =,1X,5H CHI(,I4,2X,A4,3H) =,
     1           1X,6(F7.3,3H * ,A4))
 9997 FORMAT(37X,1X,6(F7.3,3H * ,A4))
      END
