C     program redaddxyz
C
C     reduced reference potentials to 1/2 or 1/3 for esp optimization
C     add bond middle xyz coordinates
C
C      2018.1.09 
C
C     ibon=999 : Delete Atoms 
C     read atom number  Capped Hydrogen Removal
C
C      delete atoms      999        delete atoms     999
C      one CHR atom      1 0 0      two CHR atoms    1 5 0
C
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION sx(2050),sy(2050),sz(2050),pot(2050)
      DIMENSION ssx(2050),ssy(2050),ssz(2050),ppot(2050)
      DIMENSION x(250),y(200),z(200),c(200),i1(200)
      dimension ib(200,10),lin(3)
      character*25 cha
      bohr=0.529177D0
C
      WRITE(6,*) '% start Reduced reference potentials to 1/2 %'
C
      IU1=1
      IU2=2
      IU3=3
      OPEN(IU1,FORM='FORMATTED',ACCESS='SEQUENTIAL')
      OPEN(IU3,FORM='FORMATTED',
     &    ACCESS='SEQUENTIAL')
C
      lin(1)=0
      lin(2)=0
      lin(3)=0
C
      kk=0
      nl=0
      write(6,*)  'reduded 1/2 or 1/3  add xyz cood'
      read(IU2,*) ibon
      if(ibon.eq.999) then
        ibon=0
        read(IU2,*) (lin(n),n=1,3)
        do i=1,3
          if(lin(i).gt.0) nl=nl+1
        enddo
      endif
C
      do i=1,ibon
        read(IU2,'(10i4)') (ib(i,m),m=1,10)
        write(6,*) '#ib',(ib(i,m),m=1,10)
      enddo

C     
C     do j=1,2050
C##        READ(IU1,'(20x,i5,a25)',end=999) NAT,cha
        READ(IU1,*,end=999) NAT,cha
C       write(IU3,'(8x,i5,a24)') NAT,cha
        kk=kk+1
        write(6,*) '#na',nat
        DO I=1,NAT
          READ(IU1,'(i2,4(f16.10,i2))') 
     &    i1(i),X(i),i2,Y(i),i3,Z(i),i4,C(i),i5
          write(6,'(i2,4(f16.10,i2))') 
     &    i1(i),X(i),i2,Y(i),i3,Z(i),i4,C(i),i5
        ENDDO 
        ki=nat+1
        write(6,*) '#ki',ki
        cc=0.1d0
        do i=1,ibon
          do m=2,10
            if(ib(i,m).ne.0)then
              mb1=ib(i,1)
              mb2=ib(i,m)
              x(ki)=(x(mb1)+x(mb2))/2.0d0
              y(ki)=(y(mb1)+y(mb2))/2.0d0
              z(ki)=(z(mb1)+z(mb2))/2.0d0
              i1(ki)=0
              c(ki)=cc
              write(6,*) mb1,mb2,ki,x(ki),y(ki),z(ki)
              cc=cc+0.01d0
              ki=ki+1
            endif
          enddo
        enddo
        write(IU3,'(8x,i5,2x,a25)') ki-1-nl,cha
        do i=1,ki-1
          if(i.eq.lin(1))then
            write(6,*) '#delete',i,lin(1)
          else if(i.eq.lin(2))then
            write(6,*) '#delete',i,lin(2)
          else if(i.eq.lin(3))then
            write(6,*) '#delete',i,lin(3)
          else
            write(IU3,'(i2,4(f16.10,i2))') 
     &      i1(i),X(i),1,Y(i),1,Z(i),1,C(i),0
          endif
        enddo

C
        READ(IU1,'(20x,i5)') NSF
        write(6,*) '% nsf',NSF
        do i=1,nsf
          READ(IU1,'(4f15.9)') SX(i),SY(i),SZ(i),POT(i)
        enddo
C
C     Reduced to 1/2
C
        if(nsf.lt.340)then
          nn=1
        else if((nsf/2).lt.340)then
          nn=2
        else
          nn=3
        endif
        write(6,*) '% nn ',nn
        k=0
        do i=1,nsf,nn
C         write(6,'(i5,4f15.9)') i,SX(i),SY(i),SZ(i),pot(i)
          k=k+1
        enddo
C
        nc=0
        nk=0
        do i=1,nsf,nn
          dis1=100.0d0
          dis2=100.0d0
          dis3=100.0d0
          if(lin(1).ne.0)then
          dis1=dsqrt(((x(lin(1))-SX(i))**2+(Y(lin(1))-SY(i))**2+
     &                (Z(lin(1))-SZ(i))**2))*bohr
          endif
          if(lin(2).ne.0)then
          dis2=dsqrt(((x(lin(2))-SX(i))**2+(Y(lin(2))-SY(i))**2+
     &                (Z(lin(2))-SZ(i))**2))*bohr
          endif
          if(lin(3).ne.0)then
          dis3=dsqrt(((x(lin(3))-SX(i))**2+(Y(lin(3))-SY(i))**2+
     &                (Z(lin(3))-SZ(i))**2))*bohr
          endif
          if(dis1.lt.2.20d0) then
          write(6,*) i,'# close 1 < 2.20'
          nc=nc+1
          else if(dis2.lt.2.20d0) then
          write(6,*) i,'# close 2 < 2.20'
          nc=nc+1
          else if(dis3.lt.2.20d0) then
          write(6,*) i,'# close 3 < 2.20'
          nc=nc+1
          else
          nk=nk+1
          ssx(nk)=sx(i)
          ssy(nk)=sy(i)
          ssz(nk)=sz(i)
          ppot(nk)=pot(i)
          endif
        enddo
        write(IU3,'(20x,i5)') k-nc
        do i=1,nk
          write(IU3,'(4f15.9)') ssx(i),ssy(i),ssz(i),ppot(i)
        enddo
C
C     enddo
C
  999 continue
      write(6,*) '% kk & k - nc & nk =',kk,k,nc,nk
      CLOSE(UNIT=IU1)
      CLOSE(UNIT=IU2)
      CLOSE(UNIT=IU3)
      WRITE(6,*) '% stop Reduced reference potentials  add xyz %'
      STOP
      END