      SUBROUTINE BONDS(P)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      PARAMETER (NATMS2=MAXPAR*MAXPAR-MAXORB*MAXORB)
      dimension P(*)
      common /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      common /SCRACH/ B(MAXORB,MAXORB), BONDAB(NATMS2)
C***********************************************************************
C
C   CALCULATES, AND PRINTS, THE BOND INDICES AND VALENCIES OF ATOMS
C
C  FOR REFERENCE, SEE "BOND INDICES AND VALENCY", J. C. S. DALTON,
C  ARMSTRONG, D.R., PERKINS, P.G., AND STEWART, J.J.P., 838 (1973)
C
C   ON INPUT
C            P = DENSITY MATRIX, LOWER HALF TRIANGLE, PACKED.
C            P   IS NOT ALTERED BY BONDS.
C
C***********************************************************************
      write(MIOB,10)
   10 format(//20X,'BOND ORDERS AND VALENCIES',//)
      K=0
      do 20 I=1,NORBS
         do 20 J=1,I
            K=K+1
            B(I,J)=P(K)
   20 B(J,I)=P(K)
      IJ = 0
      do 60 I=1,NUMAT
         L=NFIRST(I)
         LL=NLAST(I)
         do 40 J=1,I
            IJ = IJ + 1
            K=NFIRST(J)
            KK=NLAST(J)
            X=0.0
            do 30 IL=L,LL
               do 30 IH=K,KK
   30       X=X+B(IL,IH)*B(IL,IH)
   40    BONDAB(IJ)=X
         X=-BONDAB(IJ)
         do 50 J=L,LL
   50    X=X+2.D0*B(J,J)
         BONDAB(IJ)=X
   60 continue
      call VECPRT( BONDAB, NUMAT)
C
C THIS SUBROUTINE call ADDED BY VIC L. TO WRITE OUT THE BOND INDICES
C  TO THE SYBYL MOPAC OUTPUT FILE *.SYB. (29JULY86)
C  SUBROUTINE MPC_BONDS IS IN MODULE WRITE.FOR.
C
      if (ISYBYL.EQ.1) call MPCBDS(BONDAB,NUMAT,NATMS2)
      return
      end
