      SUBROUTINE DENROT
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      common /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
      common /GEOM  / GEO(3,NUMATM)
      common /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      common /ELEMTS/ ELEMNT(107)
      common /SCRACH/ B(MAXORB*MAXORB), BONDAB(MAXPAR**2-MAXORB*MAXORB)
************************************************************************
*
* DENROT PRINTS THE DENSITY MATRIX AS (S-SIGMA, P-SIGMA, P-PI) RATHER
*        THAN (S, PX, PY, PZ).
*
************************************************************************
      dimension AROT(9,9), C(3,5,5), PAB(9,9), VECT(9,9)
      dimension NATOM(MAXORB)
      dimension XYZ(3,NUMATM), IROT(5,35), ISP(9)
      character * 6 LINE(21)
      character ELEMNT*2,ATORBS(9)*7,ITEXT(MAXORB)*7,JTEXT(MAXORB)*2
      save ATORBS, IROT, ISP
      data ATORBS/'S-SIGMA','P-SIGMA','  P-PI ','  P-PI ','D-SIGMA',
     1            '  D-PI ','  D-PI ',' D-DELL',' D-DELL'/
***********************************************************************
* IROT IS A MAPPING LIST. FOR EACH ELEMENT OF AROT 5 NUMBERS ARE
* NEEDED. THESE ARE, IN ORDER, FIRST AND SECOND SUBSCRIPTS OF AROT,
* AND FIRST,SECOND, AND THIRD SUBSCRIPTS OF C, THUS THE FIRST
* LINE OF IROT DEFINES AROT(1,1)=C(1,3,3)
*
***********************************************************************
      data IROT/1,1,1,3,3, 2,2,2,4,3, 3,2,2,2,3, 4,2,2,3,3, 2,3,2,4,2,
     1          3,3,2,2,2, 4,3,2,3,2, 2,4,2,4,4, 3,4,2,2,4, 4,4,2,3,4,
     2          5,5,3,1,5, 6,5,3,4,3, 7,5,3,3,3, 8,5,3,2,3, 9,5,3,5,3,
     3          5,6,3,1,2, 6,6,3,4,2, 7,7,3,3,2, 8,6,3,2,2, 9,6,3,5,2,
     4          5,7,3,1,4, 6,7,3,4,4, 7,7,3,3,4, 8,7,3,2,4, 9,7,3,5,4,
     5          5,8,3,1,1, 6,8,3,4,1, 7,8,3,3,1, 8,8,3,2,1, 9,8,3,5,1,
     6          5,9,3,1,5, 6,9,3,4,5, 7,9,3,3,5, 8,9,3,2,5, 9,9,3,5,5/
      data ISP /1,2,3,3,4,5,5,6,6/
      call GMETRY(GEO,XYZ)
      IPRT=0
      do 120 I=1,NUMAT
         IF=NFIRST(I)
         IL=NLAST(I)
         IPQ=IL-IF-1
         II=IPQ+2
         if (II.EQ.0)go to 120
         do 10 I1=1,II
            J1=IPRT+ISP(I1)
            ITEXT(J1)=ATORBS(I1)
            JTEXT(J1)=ELEMNT(NAT(I))
            NATOM(J1)=I
   10    continue
         IPRT=J1
         if (IPQ.NE.2)IPQ=MIN(MAX(IPQ,1),3)
         do 110 J=1,I
            JF=NFIRST(J)
            JL=NLAST(J)
            JPQ=JL-JF-1
            JJ=JPQ+2
            if (JJ.EQ.0) go to 110
            if (JPQ.NE.2) JPQ=MIN(MAX(JPQ,1),3)
            do 20 I1=1,9
               do 20 J1=1,9
   20       PAB(I1,J1)=0.D0
            KK=0
            do 30 K=IF,IL
               KK=KK+1
               LL=0
               do 30 L=JF,JL
                  LL=LL+1
   30       PAB(KK,LL)=P(L+(K*(K-1))/2)
            call COE(XYZ(1,I),XYZ(2,I),XYZ(3,I),
     1                 XYZ(1,J),XYZ(2,J),XYZ(3,J),IPQ,JPQ,C,R)
            do 40 I1=1,9
               do 40 J1=1,9
   40       AROT(I1,J1)=0.D0
            do 50 I1=1,35
   50       AROT(IROT(1,I1),IROT(2,I1))=
     1            C(IROT(3,I1),IROT(4,I1),IROT(5,I1))
            L1=ISP(II)
            L2=ISP(JJ)
            do 60 I1=1,9
               do 60 J1=1,9
   60       VECT(I1,J1)=-1.D0
            do 70 I1=1,L1
               do 70 J1=1,L2
   70       VECT(I1,J1)=0.D0
            if (I.NE.J) then
               IJ=MAX(II,JJ)
               do 90 I1=1,II
                  do 90 J1=1,JJ
                     SUM=0.D0
                     do 80 L1=1,IJ
                        do 80 L2=1,IJ
   80                SUM=SUM+AROT(L1,I1)*PAB(L1,L2)*AROT(L2,J1)
   90          VECT(ISP(I1),ISP(J1))=
     1                        VECT(ISP(I1),ISP(J1))+SUM**2
            end if
            K=0
            do  100 I1=IF,IL
               K=K+1
               L=0
               do 100 J1=JF,JL
                  L=L+1
  100       if (J1.LE.I1) B(J1+(I1*(I1-1))/2)=VECT(K,L)
  110    continue
  120 continue
C
C NOW TO REMOVE ALL THE DEAD SPACE IN P, characterIZED BY -1.0
C
      LINEAR=(NORBS*(NORBS+1))/2
      L=0
      do 130 I=1,LINEAR
         if (B(I).GT.-0.1) then
            L=L+1
            B(L)=B(I)
         end if
  130 continue
C
C   PUT ATOMIC ORBITAL VALENCIES ONTO THE DIAGONAL
C
      do 160 I=1,IPRT
         SUM=0.D0
         II=(I*(I-1))/2
         do 140 J=1,I
  140    SUM=SUM+B(J+II)
         do 150 J=I+1,IPRT
  150    SUM=SUM+B((J*(J-1))/2+I)
  160 B((I*(I+1))/2)=SUM
      do 170 I=1,21
  170 LINE(I)='------'
      LIMIT=(IPRT*(IPRT+1))/2
      KK=8
      NA=1
  180 LL=0
      M=MIN0((IPRT+1-NA),6)
      MA=2*M+1
      M=NA+M-1
      write(MIOB,'(/16X,10(1X,A7,3X))')   (ITEXT(I),I=NA,M)
      write(MIOB,'(15X,10(2X,A2,I3,4X))') (JTEXT(I),NATOM(I),I=NA,M)
      write(MIOB,'(20A6)')                (LINE(K),K=1,MA)
      do 200 I=NA,IPRT
         LL=LL+1
         K=(I*(I-1))/2
         L=MIN0((K+M),(K+I))
         K=K+NA
         if ((KK+LL).LE.50) go to 190
         write(MIOB,'(''1'')')
         write(MIOB,'(/17X,10(1X,A7,3X))') (ITEXT(N),N=NA,M)
         write(MIOB,'( 17X,10(2X,A2,I3,4X))') (JTEXT(N),NATOM(N),N=NA,M)
         write(MIOB,'(20A6)') (LINE(N),N=1,MA)
         KK=4
         LL=0
  190    write(MIOB,'(1X,A7,1X,A2,I3,10F11.6)')
     %              ITEXT(I),JTEXT(I),NATOM(I),(B(N),N=K,L)
  200 continue
      if (L.GE.LIMIT) go to 210
      KK=KK+LL+4
      NA=M+1
      if ((KK+IPRT+1-NA).LE.50) go to 180
      KK=4
      write(MIOB,'(''1'')')
      go to 180
  210 return
      end
