PROGRAM CONTACT
      IMPLICIT none
      logical lex
      character*4 key
      real*4 bonding(88),r,d1,d2,xi,xj,yi,yj,zi,zj,rt,rb
      integer*4 ti,tj,isol,ik,iout,i,j,ic,nat,natu,natc,nc
      allocatable ::isol(:),ik(:),r(:,:),iout(:)
      data bonding/
     10.50,0.32,0.98,1.28,0.95,0.87,0.82,0.80,0.78,0.77,
     10.76,1.59,1.41,1.23,1.16,1.11,1.20,1.04,1.03,2.08,1.79,1.49,1.37,
     11.27,1.23,1.22,1.22,1.21,1.20,1.22,1.30,1.31,1.27,1.25,1.21,1.19,
     11.17,2.21,1.96,1.67,1.50,1.39,1.35,1.32,1.30,1.30,1.33,1.39,1.53,
     11.49,1.46,1.45,1.41,1.38,1.36,2.40,2.03,1.74,1.70,1.70,1.69,1.68,
     11.67,1.90,1.66,1.64,1.64,1.63,1.62,1.61,1.75,1.61,1.49,1.39,1.35,
     11.33,1.31,1.32,1.35,1.39,1.54,1.53,1.52,1.51,1.51,1.50,1.49,0.10/
c
      WRITE(6,3000)
3000  FORMAT('  Takes geometry FILE.X, solute defined in SOLU.LST',/,
     1       '    and selects solute and molecules in contact ',/
     1       '      Input: FILE.X(.FC,,.TEN,.TTT)',/,
     3       '     Output: CONT.X(.FC,.TTT,.TEN,ACRYSTAL.FC)',/,/)

      d1=4.0

      inquire(file='CONTACT.OPT',exist=lex)
      if(lex)then
       open(2,file='CONTACT.OPT')
100    read(2,2000,end=2200,err=2200)key
2000   format(a4)
       if(key.eq.'DIST')read(2,*)d1
       goto 100
2200   close(2)
      endif

      d2=d1**2
      write(6,606)d1
606   format(' Limit: ',F10.1,' A')

      open(2,file='SOLU.LST')
      read(2,*)natu
      allocate(isol(natu))
      read(2,*)(isol(j),j=1,natu)
      close(2)

      open(2,file='FILE.X')
      read(2,*)
      read(2,*)nat
      allocate(ik(nat),r(nat,3),iout(nat))
      do 3 i=1,nat
      iout(i)=0
3     read(2,*)ik(i),(r(i,j),j=1,3)
      close(2)

      ic=0
      do 1 i=1,nat
c     keep solute atoms:
      do 1 j=1,natu
      if(isol(j).eq.i)then
       iout(i)=1
       ic=ic+1
       goto 1
      endif
1     continue
      write(6,*)ic,' solute atoms'

c     keep atoms close to solute
      ic=0
      do 12 i=1,nat
      if(iout(i).eq.0)then
       xi=r(i,1)
       yi=r(i,2)
       zi=r(i,3)
       do 121 j=1,natu
       xj=r(isol(j),1)
       yj=r(isol(j),2)
       zj=r(isol(j),3)
       if((xi-xj)**2+(yi-yj)**2+(zi-zj)**2.lt.d2)then
        iout(i)=2
        ic=ic+1
        goto 12
       endif
121    continue
      endif
12    continue
      write(6,*)ic,' atoms close to solute'

c     keep atoms bound in a molecule close to solute (iteration)
      nc=0
99    ic=0
      do 13 i=1,nat
      if(iout(i).eq.0)then
       xi=r(i,1)
       yi=r(i,2)
       zi=r(i,3)
       ti=ik(i)+1
       if(ti.lt.1.or.ti.gt.88)call report('type out of range')
       rb=bonding(ti)
       do 2 j=1,nat
       if(iout(j).ne.0.and.j.ne.i)then
        tj=ik(j)+1
        rt=(rb+bonding(tj))**2
        xj=r(j,1)
        yj=r(j,2)
        zj=r(j,3)
        if((xi-xj)**2+(yi-yj)**2+(zi-zj)**2.lt.rt)then
c        atom i is bound to j, which is close=take it
         iout(i)=3
         ic=ic+1
         nc=nc+1
         goto 13
        endif
       endif
2      continue
      endif
13    continue
      if(ic.ne.0)goto 99
      write(6,*)nc,' bound atoms'
       
      natc=0
      do 4 i=1,nat
4     if(iout(i).gt.0)natc=natc+1
      write(6,*)natc,' selected atoms from ',nat

      open(2,file='CONT.X')
      write(2,*)'contact molecules'
      write(2,*)natc
c     solute first
      do 5 i=1,nat
5     if(iout(i).eq.1)write(2,2001)ik(i),(r(i,j),j=1,3)
2001  format(i4,3f12.6,' 0 0 0 0 0 0 0 0.0')
      do 6 i=1,nat
6     if(iout(i).gt.1)write(2,2001)ik(i),(r(i,j),j=1,3)
      close(2)

      stop
      end
      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end