program mulliken
c
c     rewrites mullican atomic charges from gaussian output
c     into tinker ff
      character*80 fn,fna1,s80,fnaff
      character*2 sy
      parameter (nat0=3000)
      integer bt(nat0,7),ttype(nat0),ind(nat0)
      dimension nb(nat0),
     1iz(nat0),x(nat0),y(nat0),z(nat0),q(nat0)
      logical lex
c
      if(iargc().eq.0)then
       write(6,*)'Gaussian output filename:'
       read(5,'(a)')fn
       write(6,*)'Tinker geometry:'
       read(5,'(a)')fna1
       write(6,*)'Tinker force field:'
       read(5,'(a)')fnaff
      else
       call getarg(1,fn)
       call getarg(2,fNA1)
       call getarg(3,fNAff)
      endif

      ic=0
      open(8,file=fn)
1     read(8,800,end=888,err=888)s80
800   format(a80)
      if(s80(2:25).eq.'Mulliken atomic charges:')then
       read(8,*)
2      read(8,800,end=888,err=888)s80
       if(s80(1:4).ne.' Sum')then
        read(s80(11:21),*,end=888,err=888)qic
        ic=ic+1
        if(ic.gt.nat0)then
         write(6,*)'Too many charges'
         stop
        endif
        q(ic)=qic
        goto 2
       else
        goto 888
       endif
      endif
      goto 1
888   close(8)
      nc=ic
       
c     do ic=1,nc
c      write(6,*)ic,q(ic)
c     enddo
      write(6,*)nc,' charges read'

c     read in tinker file:
      open(4,file=fna1)
      read(4,*)nat
      if(nat.gt.nat0)then
       write(6,*)'Too many atoms'
       close(4)
       close(3)
       stop
      endif
      do 4 ia=1,nat
      read(4,400)sy,x(ia),y(ia),z(ia),ttype(ia),(bt(ia,ib),ib=1,7)
400   format(7x,a2,2x,3f12.6,i6,7i6)
      nb(ia)=0
      do 3 ib=1,7
3     if(bt(ia,ib).ne.0)nb(ia)=nb(ia)+1
      iz(ia)=0
      if(sy.eq.' N')iz(ia)=7
      if(sy.eq.' O')iz(ia)=8
      if(sy.eq.' C')iz(ia)=6
      if(sy.eq.' H')iz(ia)=1
      if(iz(ia).eq.0)then
       write(6,*)'Unknown atom ',sy,ia
       iz(ia)=1
      endif
4     continue
      close(4)
      write(6,*)nat,' atoms in tinker file'

      inquire(file='CCT.INP',exist=lex)
      if(.not.lex)then
       write(6,*)'CCT.INP not found'
       stop
      else
       open(9,file='CCT.INP')
       read(9,*)
       read(9,*)
       read(9,9000)(ind(i),i=1,nat)
       read(9,9000)(ind(i),i=1,nat)
9000   format(20i3)
       close(9)
       write(6,*)'CCT.INP read'
      endif
 
      do ia=1,nat
       write(6,600)ia,ttype(ia)
600    format(' atom ',i4,' type ',i4)
       if(ind(ia).ne.0)then
        write(6,*)'charge of this atom will be replaced in ff'
        isub=0
        open(9,file=fnaff)
6       read(9,800,end=99,err=99)s80
        if(s80(1:4).eq.'atom')read(s80(5:17),*)i1,i2
        if(i1.eq.ttype(ia))isub=i2
        goto 6
99      close(9)
        write(6,*)'subtype found: ',isub

        open(9,file=fnaff)
5       read(9,800,end=991,err=991)s80
        if(s80(1:6).eq.'charge')then
         read(s80(7:27),*)i1,qold
         if(i1.eq.ttype(ia))then
          write(6,6009)qold,q(ind(ia))
6009      format(' Old charge: ',f8.4,' New charge: ',f8.4)
          backspace 9
          write(9,9001)i1,q(ind(ia))
9001      format('charge',i8,f14.4)
         endif
        endif
        goto 5
991     close(9)
       endif
      enddo


      stop
      end