program rwbc
      parameter (nat0=2000)
      character*80 filename,ts
      character*2 sy(nat0),sys(nat0)
      dimension r(3*nat0)
      LOGICAL*4 notemptyline
c
      qh= 0.42
      qo=-0.84
      write(6,*)' Replace waters in gaussian input by charges:'
      write(6,*)'Filename'
      write(6,*)'N5.INP'
c     read(5,1000)filename
1000  format(a80)
      filename='N5.INP'
      open(10,file=filename)
      open(11,file='INP.NEW')
c
c     rewrite header
2     read(10,1000,end=899,err=899)ts
      if(notemptyline(ts))then
       write(11,1000)ts
       goto 2
      else
       write(11,*)' nosymm charge '
       write(11,1000)ts
      endif 
c
c     name
21    read(10,1000,end=899,err=899)ts
      write(11,1000)ts
      if(notemptyline(ts))goto 21
c
c     charge and multiplicity
      read(10,1000,end=899,err=899)ts
      write(11,1000)ts
c
c     Cartesian coordinates 
      ia=0 
22    read(10,1000,end=899,err=899)ts
      if(notemptyline(ts))then
       ia=ia+1
       if(ia.gt.nat0)then
        write(6,*)'Too many atoms'
        stop
       endif
       do 3 is=1,len(ts)
3      if(ts(is:is).ne.' ')goto 4
4      sy(ia)(1:1)=ts(is:is)
       sy(ia)(2:2)=ts(is+1:is+1)
       if(sy(ia).eq.'1 ')sy(ia)='h '
       if(sy(ia).eq.'8 ')sy(ia)='o '
       if(sy(ia).eq.'H ')sy(ia)='h '
       if(sy(ia).eq.'O ')sy(ia)='o '
       open(9,file='scr')
       write(9,1000)ts(is+2:len(ts))
       rewind 9
       read(9,*)(r(3*(ia-1)+ix),ix=1,3)
       close(9)
       goto 22
      endif
899   close(10)
      write(6,*)ia,' atoms read'
      nat=ia
      do 61 i=1,ia
61    sys(i)=sy(i)
      
      do 5 ia=1,nat
      if(sy(ia).eq.'o ')then
       x=r(3*(ia-1)+1)
       y=r(3*(ia-1)+2)
       z=r(3*(ia-1)+3)
       nh=0
       do 7 ja=1,nat
       if(sy(ja).eq.'h ')then
        xj=r(3*(ja-1)+1)
        yj=r(3*(ja-1)+2)
        zj=r(3*(ja-1)+3)
        d=(x-xj)**2+(y-yj)**2+(z-zj)**2
        if(d.lt.1.3**2)then
         nh=nh+1
         if(nh.eq.1)nh1=ja
         if(nh.eq.2)nh2=ja
        endif
       endif
7      continue
       if(nh.eq.2)then
c       we struck water, mark the atoms
        sys(ia)='OQ'
        sys(nh1)='HQ'
        sys(nh2)='HQ'
       endif
      endif
5     continue
c
c     write solute atoms
      ias=0
      do 9 ia=1,nat
      if(sys(ia).ne.'OQ'.and.sys(ia).ne.'HQ')then
       ias=ias+1
       write(11,1010)sy(ia),(r(3*(ia-1)+ix),ix=1,3)
1010   format(a2,2x,3f15.8)
      endif
9     continue
       
      write(11,*)
c
c     write solvent charges
      iac=0
      do 91 ia=1,nat
      if(sys(ia).eq.'OQ')then
       iac=iac+1
       write(11,1011)(r(3*(ia-1)+ix),ix=1,3),qo
1011   format(3f15.8,f10.4)
      endif
      if(sys(ia).eq.'HQ')then
       iac=iac+1
       write(11,1011)(r(3*(ia-1)+ix),ix=1,3),qh
      endif
91    continue
      write(11,*)
      close(11)
      write(6,*)nat, ' atoms'
      write(6,*)iac, ' charges'
      write(6,*)ias, ' solute atoms'
      write(6,*)' INP.NEW written'
      stop
      end

      function notemptyline(ts)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      LOGICAL*4 notemptyline
      CHARACTER*(*) ts
c
      l=len(ts)
      notemptyline=.false.
      do 1 i=1,l
1     if(ts(i:i).ne.' ')notemptyline=.true.
      return
      end