c----------------------------------------------------------------------
      subroutine get_bonds
c----------------------------------------------------------------------
c
c          This subroutine generates from the cartesian coordinates two
c     arrays used to specify the bonds in the molecule.  IBOND and
c     JBOND are respectively the I'th and J'th atoms in a bond, while
c     NUMBONDS is the total number of bonds found.
c
c          Two atoms are considered bonded if their distance is less or
c     equal to the sum of their radii.  Due to the sequence followed in
c     the algorithm, IBOND(k) is always lower than JBOND(k).
c
c
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
      logical changed_mol
      integer bndtyp, curr_mol
      dimension co_i(3)
c
      common/bonds/numbonds,ibond(nbond),jbond(nbond),nbonds(nbond),
     &             lastbond(natom),numpoint
      common/bond_type/bndtyp(nbond)
      common/mol_num/numat,num(natom),co(3,natom),atrad(natom)
      common/multi_mol/mol_numat(100),num_mol
      common/per_tbl_num/radii(105)
c
      numbonds=0
      curr_mol = 1
      last_atom = numat
      if (num_mol.gt.1) last_atom = mol_numat(1)
c
      do i=1,numat
         changed_mol = i.gt.last_atom
         rad_i=radii(num(i))
         co_i(1)=co(1,i)
         co_i(2)=co(2,i)
         co_i(3)=co(3,i)
         if (changed_mol) then
            curr_mol  = curr_mol+1
            last_atom = last_atom+mol_numat(curr_mol)
         endif
         do j=i+1,last_atom
            sumrad2 = (radii(num(j))+rad_i)**2
            dist2   = 0.0e0
            do k=1,3
               dist2 = dist2+(co(k,j)-co_i(k))**2
            enddo
            if (dist2.le.sumrad2) then
               numbonds         = numbonds+1
               ibond(numbonds)  = i
               jbond(numbonds)  = j
               bndtyp(numbonds) = 0
            endif
         enddo
      enddo
      return
      end
c----------------------------------------------------------------------
      subroutine bond_at_list
c----------------------------------------------------------------------
c
c          This subroutine generates a list of ALL atoms bonded to each 
c     other and stores it in the NBONDS array.  It utilizes a pointer 
c     array, LASTBOND to keep track of the origin atom. i. e. :
c
c        - atoms NBONDS(1,2,...,LASTBOND(1)) are bonded to 1,
c
c        - atoms NBONDS(LASTBOND(I-1),LASTBOND(I),...,LASTBOND(I)) are
c                bonded to atom I.
c
c          NUMPOINT is the total number of pointers generated.
c
c          This subroutine requires the previous generation of the IBOND
c     and JBOND arrays.  This information is generated in the GET_BONDS, 
c     ADD_BONDS, and DEL_BONDS subroutines.
c
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
      common/bonds/numbonds,ibond(nbond),jbond(nbond),nbonds(nbond),
     &             lastbond(natom),numpoint
      common/mol_num/numat,num(natom),co(3,natom),atrad(natom)
      common/per_tbl_num/radii(105)
c
c     Clear the arrays :
c
      do i=1,nbond
         nbonds(i)=0
      enddo
c
      do i=1,natom
         lastbond(i)=0
      enddo
      nnbonds=0
c
      do i=1,numat
         do j=1,numbonds
            if (ibond(j).eq.i .or. jbond(j).eq.i) then
               nnbonds=nnbonds+1
               nbonds(nnbonds)=j
            endif
         enddo
         lastbond(i)=nnbonds
      enddo                     
      numpoint=nnbonds
      return
      end
c-----------------------------------------------------------------------
      subroutine add_bond(atom1,atom2,full)
c-----------------------------------------------------------------------
c
c         This subroutine will add the specification corresponding to
c     a bond between ATOM1 and ATOM2 to the needed arrays (IBOND, JBOND)
c     and increment the numbonds counter by one.  It will
c     then call the  BOND_AT_LIST subroutine to update the NBONDS and
c     LASTBOND arrays.  Depending on the value of the FULL logical 
c     variable it will either create a full or a partial bond.
c
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
      integer atom1,atom2, bndtyp
      logical full
c
      common /bonds/       numbonds,ibond(nbond),jbond(nbond),
     &                     nbonds(nbond),lastbond(natom),numpoint
      common /bond_type/   bndtyp(nbond)
      common /other_bonds/ initial,mod_bonds(nbond)
      common /mol_num/     numat,num(natom),co(3,natom),atrad(natom)
c
c     first figure out which atom goes first (lower number) :
c
      if (atom1.lt.atom2) then
         i=atom1
         j=atom2
      else
         i=atom2
         j=atom1
      endif
c
c     add it to the modified bonds array :
c
      mod_bonds(initial)=i
      mod_bonds(initial+1)=j
      if (full) then
         mod_bonds(initial+2)=1
      else
         mod_bonds(initial+2)=2
      endif
      initial=initial+3
c
c     then add the atoms I & J to the IBOND and JBOND arrays :
c
      numbonds=numbonds+1
      ibond(numbonds)=i
      jbond(numbonds)=j
      if (full) then
           bndtyp(numbonds) = 0
      else
           bndtyp(numbonds) = 1
      endif
      call bond_at_list
      return
      end
c-----------------------------------------------------------------------
      subroutine del_bond(atom1,atom2)
c-----------------------------------------------------------------------
c
c          This subroutine deletes the bond between ATOM1 and ATOM2.  In
c     practice it locates the pointer BONDNUM that relates the two atoms
c     in the IBOND and JBOND array, then left-shifts all the following
c     specifications in the IBOND, JBOND, and BNDTYP arrays by one
c     position.  Afterwards it left-shifts the BONDS_CO array by two, and
c     finishes by calling the BOND_AT_LIST subroutine to update the 
c     NBONDS and LASTBOND arrays.
c
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
      integer atom1,atom2,bondnum
      integer bndtyp
      logical found_bond
c
      common/bonds/numbonds,ibond(nbond),jbond(nbond),nbonds(nbond),
     &             lastbond(natom),numpoint
      common/bond_type/bndtyp(nbond)
      common/other_bonds/initial,mod_bonds(nbond)
      common/mol_num/numat,num(natom),co(3,natom),atrad(natom)
c
c     first figure out which atom goes first (lower number) :
c
      if (atom1.lt.atom2) then
           i=atom1
           j=atom2
      else
           i=atom2
           j=atom1
      endif
c
c     add it to the modified bonds list
c
      mod_bonds(initial)=i
      mod_bonds(initial+1)=j
      mod_bonds(initial+2)=-1
      initial=initial+3
c
c     figure out the number of the bond that contains these atoms
c
      found_bond=.false.
      bondnum=1
      do while (.not.found_bond .and. bondnum.le.numbonds)
           found_bond=(ibond(bondnum).eq.i .and. jbond(bondnum).eq.j)
           bondnum=bondnum+1
      enddo
c
      if (.not. found_bond) then
           call notify('Bond was not found. ')
           return
      endif
      bondnum=bondnum-1
c
c     now left-shift all the bond-data containing arrays :
c
      if (bondnum.lt.numbonds) then
           do k=bondnum,(numbonds-1)
                ibond(k)  = ibond(k+1)
                jbond(k)  = jbond(k+1)
                bndtyp(k) = bndtyp(k+1)
           enddo
      endif
c
c     then add the atoms i & j to the ibond and jbond arrays :
c
      numbonds=numbonds-1
      call bond_at_list
      return
      end