c-----------------------------------------------------------------------
      subroutine hbpotlist
c-----------------------------------------------------------------------
c
c          This subroutine will generate arrays containing the atoms
c     that have the potential to be involved in hydrogen bonds.  The
c     only requirements considered so far are that the hydrogen in
c     question {HBND_H(i)} must be bonded to a heteroatom donor {HBND_D(i)}.
c     The acceptor atom {HBND_A(i)} must also be a heteroatom.
c
c          This subroutine is to be run on initialization, after the
c     data file has been read, and after the subroutines GET_BONDS and
c     BOND_AT_LIST, since some of the data needed is generated by them.
c
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
      integer hbnd_d,hbnd_h,hbnd_a,atom,at_type,parent,currbond
      integer initial,final
c
      common /hyd_bonds/ hbndlength,hbangle,hbnd_d(natom/3),
     &			 hbnd_h(natom/3),hbnd_a(natom/3),
     &			 nhbdon,nhbhyd,nhbacc
      common/mol_num/numat,num(natom),co(3,natom),atrad(natom)
      common/bonds/numbonds,ibond(nbond),jbond(nbond),nbonds(nbond),
     &             lastbond(natom),numpoint
c
      nhbacc  = 0
      nhbdon  = 0
      nhbhyd  = 0
      initial = 1
      do atom=1,numat
	 at_type = num(atom)
         final = lastbond(atom)
	 if (at_type.eq.1) then
            do currbond=initial,final
               parent = jbond(nbonds(currbond))
               if (parent.eq.atom) parent = ibond(nbonds(currbond))
               if (num(parent).ne.1 .and. num(parent).ne.6) then
                  nhbhyd = nhbhyd+1
                  hbnd_h(nhbhyd) = atom
                  hbnd_d(nhbhyd) = parent
               endif
            enddo
         else if (at_type.ne.6) then
            nhbacc = nhbacc+1
            hbnd_a(nhbacc) = atom
         endif
         initial = final+1
      enddo
      nhbdon = nhbhyd
      return
      end
c-----------------------------------------------------------------------
      subroutine calchbonds
c-----------------------------------------------------------------------
c
c          This subroutine will do the real calculations of hydrogen
c     bonds.  The lists of potential donors, hydrogen, and acceptors
c     generated in HBPOTLIST are followed, and a hydrogen bond is
c     considered formed if the H - - - A distance is less than the
c     cutoff HBLENGTH (default 2.5 A), and the D-H - - - A angle is
c     .le.120 degrees.  The angle itself is not used directly, but
c     used in conjunction w/the law of cosines to evaluate a
c     distance.
c
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
      integer bndtyp
c     integer atom,at_type,parent,currbond
      integer hbnd_d,hbnd_h,hbnd_a
      integer initial,final,curr_mol
      integer hydrogen,donor,acceptor
      logical wireframe,ball_and_stick,need_sort,label_num,
     &        label_cha,need_lim,perspective,round_bonds,
     &        autoscl,showhb,havehb
      logical changed_mol
c     logical first_draw
c
      common/bonds/numbonds,ibond(nbond),jbond(nbond),nbonds(nbond),
     &             lastbond(natom),numpoint
      common/bond_type/bndtyp(nbond)
      common /hyd_bonds/ hbndlength,hbangle,hbnd_d(natom/3),
     &			 hbnd_h(natom/3),hbnd_a(natom/3),
     &			 nhbdon,nhbhyd,nhbacc
      common/mol_num/numat,num(natom),co(3,natom),atrad(natom)
      common/multi_mol/mol_numat(100),num_mol
      common/draw_log/wireframe,ball_and_stick,need_sort,label_num,
     &                label_cha,need_lim,perspective,round_bonds,
     &                autoscl,showhb,havehb
c
      dimension hy_co(3),do_co(3)
c
      numhbonds = 0
      curr_mol  = 1
      final     = numat
      initial   = 1
      hbln2     = hbndlength*hbndlength
      if (num_mol.gt.1) final = mol_numat(1)
      do numhyd=1,nhbhyd
         hydrogen = hbnd_h(numhyd)
         donor    = hbnd_d(numhyd)
         ddh2     = 0.0e0
         do i=1,3
            hy_co(i) = co(i,hydrogen)
            do_co(i) = co(i,donor)
            ddh2     = ddh2 + (hy_co(i)-do_co(i))*(hy_co(i)-do_co(i))
         enddo
         changed_mol = hydrogen.gt.final
         if (changed_mol) then
            initial   = final+1
            curr_mol  = curr_mol+1
            final     = final+mol_numat(curr_mol)
         endif
         do numacc=1,nhbacc
            acceptor = hbnd_a(numacc)
            if (acceptor.ge.initial .and. acceptor.le.final .and.
     &                                    acceptor.ne.donor) then
               dha2 = 0.0e0
               dda2 = 0.0e0
               do i=1,3
                  dha2 = dha2 + (hy_co(i)-co(i,acceptor))**2
                  dda2 = dda2 + (do_co(i)-co(i,acceptor))**2
               enddo
               if (dha2.le.hbln2) then
                  dalim2 = dha2 + ddh2 - (2.0e0*sqrt(dha2*ddh2)*
     &                                  cosd(hbangle))
                  if (dda2.ge.dalim2) then
                     numbonds         = numbonds+1
                     ibond(numbonds)  = hydrogen
                     jbond(numbonds)  = acceptor
                     bndtyp(numbonds) = 2
                  endif
               endif
            endif
         enddo
      enddo
      havehb = .true.
      call bond_at_list
      return
      end
c-----------------------------------------------------------------------
      subroutine clearhbonds
c-----------------------------------------------------------------------
c
c          This subroutine will eliminate ALL hydrogen bonds from
c     the IBOND and JBOND array, then call the BOND_AT_LIST subroutine
c     to regenerate the NBONDS array.
c
c
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
      integer bndtyp
      integer hbnd_d,hbnd_h,hbnd_a
      logical wireframe,ball_and_stick,need_sort,label_num,
     &        label_cha,need_lim,perspective,round_bonds,
     &        autoscl,showhb,havehb
c
c     logical first_draw
      common /bonds/      numbonds,ibond(nbond),jbond(nbond),
     &                    nbonds(nbond),lastbond(natom),numpoint
      common /bond_type/  bndtyp(nbond)
      common /hyd_bonds/  hbndlength,hbangle,hbnd_d(natom/3),
     &			  hbnd_h(natom/3),hbnd_a(natom/3),
     &			  nhbdon,nhbhyd,nhbacc
      common /draw_log/   wireframe,ball_and_stick,need_sort,label_num,
     &                    label_cha,need_lim,perspective,round_bonds,
     &                    autoscl,showhb,havehb
c
      i = 1
      do while (i.le.numbonds)
         if (bndtyp(i).eq.2) then
            do j=i,numbonds-1
               ibond(j)  = ibond(j+1)
               jbond(j)  = jbond(j+1)
               bndtyp(j) = bndtyp(j+1)
            enddo
            numbonds = numbonds-1
         else
            i = i+1
         endif
      enddo
      havehb = .false.
      showhb = .false.
      call bond_at_list
      return
      end