program cibule
      implicit none
      integer*4 nat,nats,i,il,ih,ix,j,q,crit,qb,XCHARGE,ib,ibb,is
      real*8 xi,yi,zi,xj,yj,zj,tol,d,tol0
      integer*4,allocatable::ind(:),it(:),bt(:,:),nt(:),itb(:),btb(:,:),
     1ntb(:),indi(:)
      real*8,allocatable::r(:),rb(:)
      character*80 s80
      character*1 la
      character*1,allocatable:: si(:)
      write(6,600)
 600  format('  Cibule: Oniom for fragments',/,/,
     1       '  Input files:   FILE.X     big geometry',/,
     1       '                 FRAG.X     fragment geometry',/,
     1       '                 FRAG.OV    fragment overlap',/,
     1       '                 G.TXT      Gaussian header',/,/,
     1       '  Output     :   G.INP      Gaussian input',/)
      tol=1.6d0
      tol0=0.1d0
      crit=2
      open(7,file='FILE.X',status='old')
      read(7,*)
      read(7,*)nat
      allocate(ind(nat),rb(3*nat),itb(nat),si(nat),btb(nat,7),ntb(nat))
      q=0
      do 10 i=1,nat
      read(7,*)itb(i),(rb(ix+3*(i-1)),ix=1,3),(btb(i,ix),ix=1,7)
      q=q+itb(i)
      ntb(i)=0
      do 10 ix=1,7
10    if(btb(i,ix).ne.0)ntb(i)=ntb(i)+1
      write(6,601)nat,'FILE.X',q
601   format(i6,' atoms in ',a6,', sum of atomic charges:',i6)
      qb=XCHARGE(nat,itb,btb,ntb)
      close(7)

      open(7,file='FRAG.X',status='old')
      read(7,*)
      read(7,*)nats
      allocate(r(3*nats),it(nats),indi(nats),nt(nats),bt(nats,7))
      q=0
      do 11 i=1,nats
      read(7,*)it(i),(r(ix+3*(i-1)),ix=1,3)
11    q=q+it(i)
      write(6,601)nats,'FRAG.X',q
      close(7)

      open(7,file='FRAG.OV',status='old')
      read(7,*)ind
      read(7,*)ind
      close(7)
      do 4 i=1,nats
      indi(i)=0
      do 4 j=1,nat
 4    if(ind(j).eq.i)indi(i)=j

c     bond table of FRAG.X according to FILE.X:
      nt=0
      bt=0
      do 5 i=1,nats
      ib=indi(i)
      if(ib.ne.0)then
       do 51 ix=1,ntb(ib)
       ibb=btb(ib,ix)
       is=ind(ibb)
       if(is.ne.0)then
        nt(i)=nt(i)+1
        bt(i,nt(i))=is
       endif
 51    continue
      endif
 5    continue

      q=XCHARGE(nats,it,bt,nt)

      open(8,file='G.INP')

      open(7,file='G.TXT',status='old')
      il=0
1     read(7,80,end=88,err=88)s80
      il=il+1
      write(8,80)s80
80    format(a80)
      goto 1
88    close(7)
      write(6,602)il
602   format(i6,' lines in G.TXT')

      write(6,607)
607   format(' Charge and multiplicities all    high    low:')
      write(6,608)qb,1,q,1,qb-q,1
608   format(24x,3(2i3,2x))
      write(8,608)qb,1,q,1,qb-q,1

c     High or low level atom:
      ih=0
      do 3 i=1,nat
      if(crit.eq.1)then
c      criterium I ... based on overlap map
       if(ind(i).ne.0)then
        si(i)='H'
        ih=ih+1
       else
        si(i)='L'
       endif
      else
c      criterium II ... based on actual overlap
c      is this atom in the fragment?
       call cc(xi,yi,zi,i,rb)
       do 22 j=1,nats
       call cc(xj,yj,zj,j,r)
       d=dsqrt((xi-xj)**2+(yi-yj)**2+(zi-zj)**2)
       if(d.lt.tol0)then
        si(i)='H'
        ih=ih+1
        goto 222
       endif
 22    continue
       si(i)='L'
222    continue
      endif
3     continue
      write(6,603)ih
603   format(i6,' high level atoms',/)

      il=0
      do 2 i=1,nat
      la=' '
      call cc(xi,yi,zi,i,rb)
      do 21 j=1,nat
      call cc(xj,yj,zj,j,rb)
      d=dsqrt((xi-xj)**2+(yi-yj)**2+(zi-zj)**2)
      if(d.lt.tol.and.si(j).eq.'H'.and.si(i).eq.'L')then
       la='H'
       il=il+1
       write(6,605)il,i,j
605    format(i6,'. link between ',i5,'L and',i5,'H')
       goto 101
      endif
21    continue
101   write(8,81)itb(i),xi,yi,zi,si(i),la    
81    format(i4,3f12.6,2x,a1,1x,a1)
2     continue

      close(7)
      write(6,604)il
604   format(/,i6,' link atoms')

      write(8,*)
      close(8)
      end

      subroutine cc(x,y,z,i,r)
      real*8 x,y,z,r(*)
      integer*4 i
      x=r(1+3*(i-1))
      y=r(2+3*(i-1))
      z=r(3+3*(i-1))
      return
      end

      function XCHARGE(nat,ity,bt,nt)
c     find molecular charge in a molecule
      implicit none
      integer*4 nat,ia,q,ix,o1,o2,no,nh,n1,n2,n3,icoo,inp,iarg,ib,nn,
     1ihis,nnh,nnc,ic,c41,c42,nc4,n5,nn5,nh2,XCHARGE
      integer*4 ity(*),bt(nat,7),nt(*)
      q=0
      icoo=0
      iarg=0
      inp=0
      ihis=0
      do 3 ia=1,nat
      if(ity(ia).eq.6.and.nt(ia).eq.3)then
       o1=0
       o2=0
       no=0
       nn=0
       n1=0
       n2=0
       n3=0
       nh=0
       do 4 ix=1,nt(ia)
       ib=bt(ia,ix)
       if(ity(ib).eq.1)then
        nh=nh+1
       endif
       if(ity(ib).eq.8)then
        no=no+1
        if(o1.eq.0)then
         o1=ib
        else
         if(o2.eq.0)o2=ib
        endif
       endif
       if(ity(ib).eq.7)then
        nn=nn+1
        if(n1.eq.0)then
         n1=ib
        else
         if(n2.eq.0)then
          n2=ib
         else
          if(n3.eq.0)n3=ib
         endif
        endif
       endif
4      continue
C      Histidine+:
       if(n1.gt.0.and.n2.gt.0)then
       if(nn.eq.2.and.nh.eq.1.and.nt(n1).eq.3.and.nt(n2).eq.3)then
        nnh=0
        nnc=0
        ic=0
        do 8 ix=1,nt(n1)
        ib=bt(n1,ix)
        if(ity(ib).eq.1)nnh=nnh+1
        if(ity(ib).eq.6.and.ib.ne.ia)then
         nnc=nnc+1
         ic=ib
        endif
8       continue
        nh2=0
        do 11 ix=1,nt(n2)
11      if(ity(bt(n2,ix)).eq.1)nh2=nh2+1
        if(nnh.eq.1.and.nh2.eq.1.and.nnc.eq.1.and.nt(ic).eq.3)then
         nc4=0
         c41=0
         c42=0
         do 9 ix=1,nt(ic)
         ib=bt(ic,ix)
         if(ity(ib).eq.6)then
          nc4=nc4+1
          if(c41.eq.0)then
           c41=ib
          else
           c42=ib
          endif
         endif
9        continue
         if(nc4.eq.1.or.nc4.eq.2)then
          if(nt(c41).eq.3)then
           n5=0
           nn5=0
           do 10 ix=1,nt(c41)
           ib=bt(c41,ix)
           if(ity(ib).eq.7)then
            nn5=nn5+1
            n5=ib
           endif
10         continue
           if(nn5.eq.1.and.n5.eq.n2)then
            q=q+1
            ihis=ihis+1
           endif
          endif
          if(c42.ne.0)then
          if(c42.ne.0.and.nt(c42).eq.3)then
           n5=0
           nn5=0
           do 12 ix=1,nt(c42)
           ib=bt(c42,ix)
           if(ity(ib).eq.7)then
            nn5=nn5+1
            n5=ib
           endif
12         continue
           if(nn5.eq.1.and.n5.eq.n2)then
            q=q+1
            ihis=ihis+1
           endif
          endif
          endif
         endif
        endif
       endif
       endif
c      COO-:
       if(o1.ne.0.and.o2.ne.0)then
       if(no.eq.2.and.nt(o1).eq.1.and.nt(o2).eq.1)then
        q=q-1
        icoo=icoo+1
       endif
       endif
c      ARG+:
       if(n1.ne.0.and.n2.ne.0)then
       if(nn.eq.3.and.nt(ia).eq.3.and.nt(n1).eq.3.and.nt(n2).eq.3)then
        nh=0
        do 5 ix=1,nt(n1)
        ib=bt(n1,ix)
5       if(ity(ib).eq.1)nh=nh+1
        do 6 ix=1,nt(n2)
        ib=bt(n2,ix)
6       if(ity(ib).eq.1)nh=nh+1
        do 7 ix=1,nt(n3)
        ib=bt(n3,ix)
7       if(ity(ib).eq.1)nh=nh+1
        if(nh.eq.5)then
         q=q+1
         iarg=iarg+1
        endif
       endif
       endif
      endif

c     N+:
      if(ity(ia).eq.7.and.nt(ia).eq.4)then
       q=q+1
       inp=inp+1
      endif

3     continue

      XCHARGE=q

      write(6,600)q,icoo,inp,iarg,ihis
600   format(5i5,' (Q NCOO- NN+ NARG+ NHIS+)')
      return
      end