      program top_sym
      implicit none
      integer natoms,attract_label
      integer, dimension (:), allocatable :: nat
      double precision, dimension (:), allocatable :: xat,yat,zat
      double precision xmin,xmax,ymin,ymax,zmin,zmax,vol
      double precision xminold,xmaxold,yminold,ymaxold,zminold,zmaxold
      double precision trans(3,3),dx,dy,dz
      double precision, dimension (:), allocatable :: volume 
      double precision xc, yc, zc, dist
      double precision v(64,3)
      character*16, dimension (:), allocatable :: basin_name
      character*4, dimension (:), allocatable :: atom_name
      character*40 file_in
      character*1 molecule(40),blank
      character*3 file_type
      integer, dimension (:,:), allocatable :: index_basin,aux_label
      integer indexi,indexk
      integer npx,npy,npz,nop,sym(3,3,3)
      integer npxold,npyold,npzold
      integer i, j, k, l, m, n, maxv, iop, np,npts
      integer ix1,ix2,iy1,iy2,iz1,iz2,ix,iy,iz,ijk
      logical if_file
      integer*2, dimension (:,:,:), allocatable :: ig
      integer*2, dimension (:), allocatable :: iaux
      real, dimension (:,:,:), allocatable :: g
      real, dimension (:), allocatable :: aux
      real r(8,3)
      double precision, dimension (:), allocatable :: x,y,z
      character*3 name_op(3)
      common /mat/xc,yc,zc,trans,sym,name_op
      equivalence (molecule(1),file_in)
      data blank/' '/
      write(*,*) ' input sbf file '
      read(*,'(a)') file_in
      inquire(file=file_in,exist=if_file)
      if(.not.if_file) then
        write(*,'(2x,"requested file does not exist:",a40)') file_in
        stop
      endif
      i=1
      do while (molecule(i).ne.blank)
         i=i+1
      enddo
      file_type=molecule(i-7)//molecule(i-6)//molecule(i-5)
      open(unit=1,file=file_in,status='old',form='unformatted')
      read(1) npxold,npyold,npzold
      read(1) xminold,xmaxold,yminold,ymaxold,zminold,zmaxold
      write(*,'(2x,"old limits ",6f10.3)')
     .                  xminold,xmaxold,yminold,ymaxold,zminold,zmaxold
      write(*,*) ' number of symmetry operations '
      read(*,*) nop
      do i=1,nop
         write(*,*) ' symmetry operations symbol '
         read(*,'(a3)') name_op(i)
         write(*,*) ' symmetry operations coordinates '
         read(*,*) xc,yc,zc
         call symmat(i)
      enddo
      do i=1,nop
         write(*,'(//)')
         do j=1,3
            write(*,'(3i5,5x,f12.6)')(sym(i,j,k),k=1,3),trans(i,j)
         enddo
      enddo
      dx=(xmaxold-xminold)/(npxold-1)
      dy=(ymaxold-yminold)/(npyold-1)
      dz=(zmaxold-zminold)/(npzold-1)
      vol=dx*dy*dz
      v(1,1)=xminold
      v(1,2)=yminold
      v(1,3)=zminold
      v(2,1)=xminold
      v(2,2)=yminold
      v(2,3)=zmaxold
      v(3,1)=xminold
      v(3,2)=ymaxold
      v(3,3)=zminold
      v(4,1)=xminold
      v(4,2)=ymaxold
      v(4,3)=zmaxold
      v(5,1)=xmaxold
      v(5,2)=yminold
      v(5,3)=zminold
      v(6,1)=xmaxold
      v(6,2)=yminold
      v(6,3)=zmaxold
      v(7,1)=xmaxold
      v(7,2)=ymaxold
      v(7,3)=zminold
      v(8,1)=xmaxold
      v(8,2)=ymaxold
      v(8,3)=zmaxold
      do i=1,nop
         maxv=2**(i+2)
         do j=1,maxv
            l=maxv+j
            do k=1,3
               v(l,k)=trans(i,k)
               do m=1,3
                  v(l,k)=v(l,k)+sym(i,k,m)*v(j,m)
               enddo
            enddo
         enddo
      enddo
      maxv=2*maxv
      xmin=xminold
      ymin=yminold
      zmin=zminold
      xmax=xmaxold
      ymax=ymaxold
      zmax=zmaxold
      do i=1,maxv
         if(v(i,1).lt.xmin) xmin=v(i,1)
         if(v(i,1).gt.xmax) xmax=v(i,1)
         if(v(i,2).lt.ymin) ymin=v(i,2)
         if(v(i,2).gt.ymax) ymax=v(i,2)
         if(v(i,3).lt.zmin) zmin=v(i,3)
         if(v(i,3).gt.zmax) zmax=v(i,3)
      enddo
      npx=nint((xmax-xmin)/dx+1)
      npy=nint((ymax-ymin)/dy+1)
      npz=nint((zmax-zmin)/dz+1)
      npts=npx*npy*npz
      if(file_type.eq.'elf') then
         allocate(aux(npts))
         npts=npxold*npyold*npzold
         read(1)(aux(i),i=1,npts)
      else if(file_type.eq.'bas') then
         read(1) natoms,attract_label
         allocate(xat(natoms),yat(natoms),zat(natoms),nat(natoms))
         allocate(atom_name(natoms))
         allocate(volume(attract_label),basin_name(attract_label))
         allocate(index_basin(attract_label,5),
     .            aux_label(attract_label,64))
         allocate(x(attract_label),y(attract_label),z(attract_label))
         do i=1,natoms
            read(1) nat(i),atom_name(i),xat(i),yat(i),zat(i)
         enddo
         do i=1,attract_label
            read(1) index_basin(i,1),basin_name(i),volume(i),
     .              x(i),y(i),z(i)
         enddo
         allocate(iaux(npts))
         npts=npxold*npyold*npzold
         read(1)(iaux(i),i=1,npts)
      else if(file_type.eq.'syn') then
         allocate(iaux(npts))
         npts=npxold*npyold*npzold
         read(1)(iaux(i),i=1,npts)
      else if(file_type.eq.'rho') then
         allocate(aux(npts))
         npts=npxold*npyold*npzold
         read(1)(aux(i),i=1,npts)
      else if(file_type.eq.'lap') then
         allocate(aux(npts))
         npts=npxold*npyold*npzold
         read(1)(aux(i),i=1,npts)
      endif
      rewind 1
      close(unit=1)
      write(*,'(2x,"new limits ",6f10.3)') xmin,xmax,ymin,ymax,zmin,
     .                                     zmax
      write(*,'(3i5)') npxold,npyold,npzold
      write(*,'(3i5)') npx,npy,npz
      ix1=nint((xminold-xmin)/dx)
      ix2=nint((xmaxold-xmin)/dx)
      iy1=nint((yminold-ymin)/dy)
      iy2=nint((ymaxold-ymin)/dy)
      iz1=nint((zminold-zmin)/dz)
      iz2=nint((zmaxold-zmin)/dz)
      ijk=0
      if((file_type.eq.'elf').or.(file_type.eq.'rho')
     .                       .or.(file_type.eq.'lap')) then
         allocate(g(0:npx,0:npy,0:npz))
         do i=ix1,ix2
            do j=iy1,iy2
               do k=iz1,iz2
                  ijk=ijk+1
                   g(i,j,k)=aux(ijk)
               enddo
            enddo
         enddo
      else
         allocate(ig(0:npx,0:npy,0:npz))
         do i=ix1,ix2
            do j=iy1,iy2
               do k=iz1,iz2
                  ijk=ijk+1
                  ig(i,j,k)=iaux(ijk)
               enddo
            enddo
         enddo
      endif
      if((file_type.eq.'elf').or.(file_type.eq.'rho')
     .                       .or.(file_type.eq.'lap')) then
         do i=1,npts
            aux(i)=0.0
         enddo
      else
         do i=1,npts
            iaux(i)=0
         enddo
      endif
      if((file_type.eq.'elf').or.(file_type.eq.'rho')
     .                       .or.(file_type.eq.'lap')) then
      do i=ix1,ix2
         r(1,1)=xmin+i*dx
         do j=iy1,iy2
            r(1,2)=ymin+j*dy
            do k=iz1,iz2
               r(1,3)=zmin+k*dz
               maxv=1
               do iop=1,nop
                  do np=1,maxv
                     l=maxv+np
                     do m=1,3
                        r(l,m)=trans(iop,m)
                        do n=1,3
                           r(l,m)=r(l,m)+sym(iop,m,n)*r(np,n)
                        enddo
                     enddo
                  enddo
                  maxv=maxv+maxv
               enddo
               do np=1,maxv
                  ix=nint((r(np,1)-xmin)/dx)
                  iy=nint((r(np,2)-ymin)/dy)
                  iz=nint((r(np,3)-zmin)/dz)
                  g(ix,iy,iz)=g(i,j,k)
               enddo
            enddo
         enddo
      enddo
      else if(file_type.eq.'bas') then
c        allocate(ig(0:npx,0:npy,0:npz))
         do i=1,attract_label
            r(1,1)=x(i)
            r(1,2)=y(i)
            r(1,3)=z(i)
            maxv=1
            aux_label(i,1)=i
            indexi=index_basin(i,1)
            do iop=1,nop
               do np=1,maxv
                  l=maxv+np
                  do m=1,3
                     r(l,m)=trans(iop,m)
                     do n=1,3
                        r(l,m)=r(l,m)+sym(iop,m,n)*r(np,n)
                     enddo
                  enddo
                  aux_label(i,l)=i
                  do k=1,attract_label
                     indexk=index_basin(k,1)
                     dist=(x(k)-r(l,1))**2+(y(k)-r(l,2))**2+
     .                    (z(k)-r(l,3))**2
                     if((dist.lt.1.0d-1).and.(indexi.eq.indexk)) then
                        aux_label(i,l)=k
                     endif
                  enddo
                  k=aux_label(i,l)
               enddo
               maxv=maxv+maxv
            enddo
         enddo
         do i=ix1,ix2
            r(1,1)=xmin+i*dx
            do j=iy1,iy2
               r(1,2)=ymin+j*dy
               do k=iz1,iz2
                  r(1,3)=zmin+k*dz
                  maxv=1
                  do iop=1,nop
                     do np=1,maxv
                        l=maxv+np
                        do m=1,3
                           r(l,m)=trans(iop,m)
                           do n=1,3
                              r(l,m)=r(l,m)+sym(iop,m,n)*r(np,n)
                           enddo
                        enddo
                     enddo
                     maxv=maxv+maxv
                  enddo
                  do np=1,maxv
                     ix=nint((r(np,1)-xmin)/dx)
                     iy=nint((r(np,2)-ymin)/dy)
                     iz=nint((r(np,3)-zmin)/dz)
                     l=ig(i,j,k)
                     if(l.gt.0) then
                        ig(ix,iy,iz)=aux_label(l,np)
                     else
                        ig(ix,iy,iz)=0
                     endif
                  enddo
               enddo
            enddo
         enddo
      else
c     allocate(ig(0:npx,0:npy,0:npz))
      do i=ix1,ix2
         r(1,1)=xmin+i*dx
         do j=iy1,iy2
            r(1,2)=ymin+j*dy
            do k=iz1,iz2
               r(1,3)=zmin+k*dz
               maxv=1
               do iop=1,nop
                  do np=1,maxv
                     l=maxv+np
                     do m=1,3
                        r(l,m)=trans(iop,m)
                        do n=1,3
                           r(l,m)=r(l,m)+sym(iop,m,n)*r(np,n)
                        enddo
                     enddo
                  enddo
                  maxv=maxv+maxv
               enddo
               do np=1,maxv
                  ix=nint((r(np,1)-xmin)/dx)
                  iy=nint((r(np,2)-ymin)/dy)
                  iz=nint((r(np,3)-zmin)/dz)
                  ig(ix,iy,iz)=ig(i,j,k)
               enddo
            enddo
         enddo
      enddo
      endif
      open(unit=2,file=file_in,status='old',form='unformatted')
      rewind 2
      write(2) npx,npy,npz
      write(2) xmin,xmax,ymin,ymax,zmin,zmax
      npts=npx*npy*npz
      ijk=0
      if((file_type.eq.'elf').or.(file_type.eq.'rho')
     .                       .or.(file_type.eq.'lap')) then
      do i=0,npx-1
         do j=0,npy-1
            do k=0,npz-1
               ijk=ijk+1
               aux(ijk)=g(i,j,k)
            enddo
         enddo
      enddo
      else
      do i=0,npx-1
         do j=0,npy-1
            do k=0,npz-1
               ijk=ijk+1
               iaux(ijk)=ig(i,j,k)
            enddo
         enddo
      enddo
      do i=1,attract_label
         volume(i)=0.0d0
      enddo
      do i=1,npts
         k=iaux(i)
         if(k.gt.0) then
             volume(k)=volume(k)+vol
         endif
      enddo
      endif
      if(file_type.eq.'bas') then
         write(2) natoms,attract_label
         do i=1,natoms
            write(2) nat(i),atom_name(i),xat(i),yat(i),zat(i)
         enddo
         do i=1,attract_label
            write(2) index_basin(i,1),basin_name(i),volume(i),
     .               x(i),y(i),z(i)
         enddo
      endif
      if((file_type.eq.'elf').or.(file_type.eq.'rho')
     .                       .or.(file_type.eq.'lap')) then
         write(2)(aux(i),i=1,npts)
      else
         write(2)(iaux(i),i=1,npts)
      endif
      close(unit=2)
      end
