program xsn
      character*80 s80,fn
      parameter (nat0=100000)
      dimension r(nat0,3),jo(nat0),jh(nat0)
      integer*4 as(nat0),ll(nat0),nll

c     write(6,600)
600   format(' Extracts solute atomic numbers from FILE.X')
c     nll ... number of atoms
c     ll... list of atoms

      ilone=0
      idiat=0
      is=0
      if(iargc().eq.0)then
       fn='FILE.X'
      else
       call getarg(1,fn)
      endif

      open(3,file=fn,status='old')
1     read(3,300,end=999,err=999)s80
300   format(a80)
      read(3,*,end=999,err=999)nat
      if(nat.gt.nat0)call report('too many atoms')
      is=is+1
c
c     read in all atoms:
      do 2 i=1,nat
2     read(3,*)as(i),(r(i,k),k=1,3)

      nll=0
      do 3 i=1,nat
      xi=r(i,1)
      yi=r(i,2)
      zi=r(i,3)
c
c     io, in, ih,it ... how many covalently bond Os,Ns,Hs and any atoms:
      io=0
      in=0
      ih=0
      it=0
      do 4 j=1,nat
      if(j.ne.i)then
       xj=r(j,1)
       yj=r(j,2)
       zj=r(j,3)
       dij=sqrt((xi-xj)**2+(yi-yj)**2+(zi-zj)**2)
       if(dij.lt.1.2.and.as(j).eq.8)then
        io=io+1
        jo(io)=j
       endif
       if(dij.lt.1.4.and.as(j).eq.7)in=in+1
       if(dij.lt.1.2.and.as(j).eq.1)then
        ih=ih+1
        jh(ih)=j
       endif
       if((dij.lt.1.6.and.(as(j).ne.1.and.as(i).ne.1))
     1 .or.(dij.lt.1.2.and.(as(j).eq.1.or.as(i).eq.1)))then
        it=it+1
        if(it.eq.1)ifirst=j
       endif
      endif
4     continue
c
c     water oxygen
      if(as(i).eq.8.and.ih.eq.2)goto 3
c
      if(as(i).eq.1.and.io.eq.1)then
c      check hydrogen bond to oxygen:
       xk=r(jo(1),1)
       yk=r(jo(1),2)
       zk=r(jo(1),3)
c      look for second hydrogen:
       do 41 jj=1,nat
       if(jj.ne.i.and.jj.ne.jo(1))then
        xl=r(jj,1)
        yl=r(jj,2)
        zl=r(jj,3)
        doh=sqrt((xk-xl)**2+(yk-yl)**2+(zk-zl)**2)
c       water hydrogen:
        if(doh.lt.1.2.and.as(jj).eq.1)goto 3
       endif
41     continue
      endif
c
c     check for diatomic
      if(it.eq.1)then
       idia=0
       xk=r(ifirst,1)
       yk=r(ifirst,2)
       zk=r(ifirst,3)
       do 51 jj=1,nat
       if(jj.ne.i.and.jj.ne.ifirst)then
        xl=r(jj,1)
        yl=r(jj,2)
        zl=r(jj,3)
        df=sqrt((xk-xl)**2+(yk-yl)**2+(zk-zl)**2)
        if((df.lt.1.6.and.as(jj).ne.1).or.(as(jj).eq.1.and.df.lt.1.2))
     1  idia=idia+1
       endif
51     continue
       if(idia.eq.0)then
        idiat=idiat+1
        goto 3
       endif
      endif
c
c     lone atoms:
      if(it.lt.1)then
       ilone=ilone+1
       goto 3
      endif

      nll=nll+1
      ll(nll)=i
3     continue
      write(6,6001)is
      write(6,*)
      write(6,6001)-nll
      write(6,*)
      write(6,6001)(ll(i),i=1,nll)
6001  format(i5,$)
      write(6,*)

c     go to next geometry in FILE.X:
      goto 1

999   close(3)

      write(6,6003)is,ilone,idiat
6003  format(/,i5,' structures',/,
     1         i5,' lone atoms',/,
     2         i5,' diatomic cases')
      stop
      end

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