      subroutine find_ints(natoms,nbonds,nangles,ndihedrals,nimpropers,
     &                     natomtype,i_bond,i_angle,i_dihedral,
     &                     i_improper,iwkmat,coord)

c  this subroutine written by alain st-amant of the
c  department of chemistry, university of ottawa, ottawa, canada.
c  all rights reserved.  this is part of the DeFT project.

      implicit real*8(a-h,o-z)

      dimension radii(0:54),radius(0:54)

      dimension natomtype(*),i_bond(2,*),i_angle(3,*),i_dihedral(4,*),
     &          i_improper(4,*),iwkmat(natoms,*)

      dimension coord(3,*)

      data zero,bohr,ang_cutoff,one80,fact/0.0d0,0.529177d0,
     &                                     170.0d0,180.0d0,1.3d0/
      data perturb/0.99999d0/

      data radii/1.00d0,0.32d0,0.60d0,1.20d0,1.05d0,0.81d0,0.77d0,
     &                  0.74d0,0.74d0,0.72d0,0.72d0,1.50d0,1.40d0,
     &                  1.30d0,1.17d0,1.10d0,1.04d0,0.99d0,0.99d0,
     &                  1.80d0,1.60d0,1.40d0,1.40d0,1.40d0,1.40d0,
     &                  1.40d0,1.40d0,1.40d0,1.40d0,1.40d0,1.40d0,
     &                  1.40d0,1.30d0,1.20d0,1.20d0,1.10d0,1.10d0,
     &                  2.20d0,2.00d0,1.60d0,1.60d0,1.60d0,1.60d0,
     &                  1.60d0,1.60d0,1.60d0,1.60d0,1.60d0,1.60d0,
     &                  1.55d0,1.45d0,1.45d0,1.40d0,1.40d0,1.40d0/

      do 1001 i=0,54
 1001 radius(i)=radii(i)/bohr

      nbonds=0
      nangles=0
      ndihedrals=0
      nimpropers=0

      do 1002 i=1,natoms
      do 1002 j=1,i

      if(i.eq.j) goto 1002

      r_ab=fact*(radius(natomtype(i))+radius(natomtype(j)))

      if(dist(coord(1,i),coord(1,j)).le.r_ab) then
                                                    nbonds=nbonds+1

                                                    i_bond(1,nbonds)=i
                                                    i_bond(2,nbonds)=j
                                              endif

 1002 continue

      do 1003 j=1,natoms
      do 1004 i=1,natoms
 1004 iwkmat(i,j)=0
 1003 iwkmat(j,j)=1

 1005 continue

      do 1006 i=1,nbonds

      ii=i_bond(1,i)
      jj=i_bond(2,i)

      iwkmat(ii,jj)=1
      iwkmat(jj,ii)=1

      do 1007 j=1,natoms
 1007 iwkmat(j,ii)=iwkmat(j,ii)+iwkmat(j,jj)

      do 1008 j=1,natoms
 1008 iwkmat(j,jj)=iwkmat(j,ii)

      do 1009 j=1,natoms
 1009 iwkmat(ii,j)=iwkmat(ii,j)+iwkmat(jj,j)

      do 1010 j=1,natoms
 1010 iwkmat(jj,j)=iwkmat(ii,j)

 1006 continue

      icheck=0

      do 1011 j=1,natoms
      do 1011 i=1,natoms
 1011 if(iwkmat(i,j).eq.0) icheck=1

      if(icheck.eq.0) goto 1012

      do 1013 i=1,natoms

      jcheck=0

      do 1014 j=1,natoms
 1014 if(iwkmat(j,i).eq.0) jcheck=j

      if(jcheck.eq.0) goto 1013

      dist_min=one80**10

      do 1015 j=1,natoms
      do 1015 k=1,natoms
      if((iwkmat(j,i).ne.0).and.
     &   (iwkmat(k,jcheck).ne.0).and.
     &   (iwkmat(k,i).eq.0).and.
     &   (iwkmat(j,jcheck).eq.0)) then

                           distance=dist(coord(1,j),coord(1,k))
                           if(distance.le.dist_min) dist_min=distance

                                  endif
 1015 continue

      dist_min=perturb*dist_min

      do 1016 j=1,natoms
      do 1016 k=1,j
      distance=dist(coord(1,j),coord(1,k))
      if((distance.ge.dist_min).and.(distance.le.fact*dist_min)) then

                     nbonds=nbonds+1

                     i_bond(1,nbonds)=j
                     i_bond(2,nbonds)=k

                     iwkmat(j,k)=1
                     iwkmat(k,j)=1

                     do 1017 l=1,natoms
 1017                iwkmat(l,j)=iwkmat(l,j)+iwkmat(l,k)

                     do 1018 l=1,natoms
 1018                iwkmat(l,k)=iwkmat(l,j)

                     do 1019 l=1,natoms
 1019                iwkmat(j,l)=iwkmat(j,l)+iwkmat(k,l)

                     do 1020 l=1,natoms
 1020                iwkmat(k,l)=iwkmat(j,l)

                                                                 endif
 1016 continue

      goto 1005

 1013 continue

 1012 continue

      nbonds_tmp=nbonds

      do 1021 i=1,nbonds
      do 1021 j=1,i

      if(i.eq.j) goto 1021

      ind_1i=i_bond(1,i)
      ind_2i=i_bond(2,i)
      ind_1j=i_bond(1,j)
      ind_2j=i_bond(2,j)

      if(ind_1i.eq.ind_1j) then
                                 a1=ang(coord(1,ind_2i),
     &                                  coord(1,ind_1i),
     &                                  coord(1,ind_2j))

                                 if(a1.gt.ang_cutoff) then

                                              nbonds=nbonds+1

                                              i_bond(1,nbonds)=ind_2i
                                              i_bond(2,nbonds)=ind_2j

                                              goto 1021

                                                      endif

                                 nangles=nangles+1

                                 i_angle(1,nangles)=ind_2i
                                 i_angle(2,nangles)=ind_1i
                                 i_angle(3,nangles)=ind_2j

                                 goto 1021
                           endif

      if(ind_1i.eq.ind_2j) then
                                 a1=ang(coord(1,ind_2i),
     &                                  coord(1,ind_1i),
     &                                  coord(1,ind_1j))

                                 if(a1.gt.ang_cutoff) then

                                              nbonds=nbonds+1

                                              i_bond(1,nbonds)=ind_2i
                                              i_bond(2,nbonds)=ind_1j

                                              goto 1021

                                                      endif

                                 nangles=nangles+1

                                 i_angle(1,nangles)=ind_2i
                                 i_angle(2,nangles)=ind_1i
                                 i_angle(3,nangles)=ind_1j

                                 goto 1021
                           endif

      if(ind_2i.eq.ind_1j) then
                                 a1=ang(coord(1,ind_1i),
     &                                  coord(1,ind_2i),
     &                                  coord(1,ind_2j))

                                 if(a1.gt.ang_cutoff) then

                                              nbonds=nbonds+1

                                              i_bond(1,nbonds)=ind_1i
                                              i_bond(2,nbonds)=ind_2j

                                              goto 1021

                                                      endif

                                 nangles=nangles+1

                                 i_angle(1,nangles)=ind_1i
                                 i_angle(2,nangles)=ind_2i
                                 i_angle(3,nangles)=ind_2j

                                 goto 1021
                           endif

      if(ind_2i.eq.ind_2j) then
                                 a1=ang(coord(1,ind_1i),
     &                                  coord(1,ind_2i),
     &                                  coord(1,ind_1j))

                                 if(a1.gt.ang_cutoff) then

                                              nbonds=nbonds+1

                                              i_bond(1,nbonds)=ind_1i
                                              i_bond(2,nbonds)=ind_1j

                                              goto 1021

                                                      endif

                                 nangles=nangles+1

                                 i_angle(1,nangles)=ind_1i
                                 i_angle(2,nangles)=ind_2i
                                 i_angle(3,nangles)=ind_1j

                                 goto 1021
                           endif

 1021 continue

      do 1022 i=1,nbonds
      do 1022 j=1,natoms
      do 1022 k=1,j

      if(k.eq.j) goto 1022

      ind_1i=i_bond(1,i)
      ind_2i=i_bond(2,i)

      if(ind_1i.eq.j) goto 1022
      if(ind_1i.eq.k) goto 1022
      if(ind_2i.eq.j) goto 1022
      if(ind_2i.eq.k) goto 1022

      dist_ab=dist(coord(1,ind_1i),coord(1,j))
      dist_cd=dist(coord(1,ind_2i),coord(1,k))

      r_ab=fact*(radius(natomtype(ind_1i))+radius(natomtype(j)))
      r_cd=fact*(radius(natomtype(ind_2i))+radius(natomtype(k)))

      if((dist_ab.le.r_ab).and.(dist_cd.le.r_cd)) then

                  if(i.gt.nbonds_tmp) then
                                            nbonds=nbonds+1

                                            i_bond(1,nbonds)=j
                                            i_bond(2,nbonds)=k

                                            goto 1022
                                      endif

                  ang1=ang(coord(1,j),coord(1,ind_1i),coord(1,ind_2i))
                  ang2=ang(coord(1,k),coord(1,ind_2i),coord(1,ind_1i))

                  angmax=max(ang1,ang2)
                  if(angmax.gt.ang_cutoff) goto 1022

                  angmin=min(ang1,ang2)
                  if(angmin.lt.(one80-ang_cutoff)) goto 1022

                  ndihedrals=ndihedrals+1

                  i_dihedral(1,ndihedrals)=j
                  i_dihedral(2,ndihedrals)=ind_1i
                  i_dihedral(3,ndihedrals)=ind_2i
                  i_dihedral(4,ndihedrals)=k

                  goto 1022

                                                  endif

      dist_ab=dist(coord(1,ind_1i),coord(1,k))
      dist_cd=dist(coord(1,ind_2i),coord(1,j))

      r_ab=fact*(radius(natomtype(ind_1i))+radius(natomtype(k)))
      r_cd=fact*(radius(natomtype(ind_2i))+radius(natomtype(j)))

      if((dist_ab.le.r_ab).and.(dist_cd.le.r_cd)) then

                  if(i.gt.nbonds_tmp) then
                                            nbonds=nbonds+1

                                            i_bond(1,nbonds)=j
                                            i_bond(2,nbonds)=k

                                            goto 1022
                                      endif

                  ang1=ang(coord(1,k),coord(1,ind_1i),coord(1,ind_2i))
                  ang2=ang(coord(1,j),coord(1,ind_2i),coord(1,ind_1i))

                  angmax=max(ang1,ang2)
                  if(angmax.gt.ang_cutoff) goto 1022

                  angmin=min(ang1,ang2)
                  if(angmin.lt.(one80-ang_cutoff)) goto 1022

                  ndihedrals=ndihedrals+1

                  i_dihedral(1,ndihedrals)=k
                  i_dihedral(2,ndihedrals)=ind_1i
                  i_dihedral(3,ndihedrals)=ind_2i
                  i_dihedral(4,ndihedrals)=j

                  goto 1022

                                                  endif

      if(i.gt.nbonds_tmp) goto 1022

      ang_max=max(ang(coord(1,ind_1i),coord(1,ind_2i),coord(1,j)),
     &            ang(coord(1,ind_1i),coord(1,ind_2i),coord(1,k)),
     &            ang(coord(1,ind_2i),coord(1,ind_1i),coord(1,j)),
     &            ang(coord(1,ind_2i),coord(1,ind_1i),coord(1,k)))

      if(ang_max.gt.ang_cutoff) goto 1022

      ang_min=min(ang(coord(1,ind_1i),coord(1,ind_2i),coord(1,j)),
     &            ang(coord(1,ind_1i),coord(1,ind_2i),coord(1,k)),
     &            ang(coord(1,ind_2i),coord(1,ind_1i),coord(1,j)),
     &            ang(coord(1,ind_2i),coord(1,ind_1i),coord(1,k)))

      if(ang_min.lt.(one80-ang_cutoff)) goto 1022

      dist_ab=dist(coord(1,ind_1i),coord(1,j))
      dist_ac=dist(coord(1,ind_1i),coord(1,k))

      r_ab=fact*(radius(natomtype(ind_1i))+radius(natomtype(j)))
      r_ac=fact*(radius(natomtype(ind_1i))+radius(natomtype(k)))

      if((ind_1i.lt.k).and.
     &   (dist_ab.le.r_ab).and.
     &   (dist_ac.le.r_ac)) then
                                  nimpropers=nimpropers+1

                                  i_improper(1,nimpropers)=ind_2i
                                  i_improper(2,nimpropers)=j
                                  i_improper(3,nimpropers)=k
                                  i_improper(4,nimpropers)=ind_1i

                                  goto 1022
                            endif

      dist_ab=dist(coord(1,ind_2i),coord(1,j))
      dist_ac=dist(coord(1,ind_2i),coord(1,k))

      r_ab=fact*(radius(natomtype(ind_2i))+radius(natomtype(j)))
      r_ac=fact*(radius(natomtype(ind_2i))+radius(natomtype(k)))

      if((ind_1i.lt.k).and.
     &   (dist_ab.le.r_ab).and.
     &   (dist_ac.le.r_ac)) then
                                  nimpropers=nimpropers+1

                                  i_improper(1,nimpropers)=ind_1i
                                  i_improper(2,nimpropers)=j
                                  i_improper(3,nimpropers)=k
                                  i_improper(4,nimpropers)=ind_2i

                                  goto 1022
                            endif

 1022 continue

      return
      end
