PROGRAM G94BBB
c
c     rewrites gaussian cube file into BBB.SCR
c
      parameter (nat0=200,ngr0=200)  
      IMPLICIT REAL*8 (A-H,O-Z)
      implicit integer*4 (i-n)
      character*80 filename,fpot,ppp
      character*40 n1,n2
      character*1 ok1,ok
      dimension q(nat0),r(3,nat0),x0(3),x1(3),x2(3),x3(3),rt(3,nat0)
      real*4  rofo(ngr0,ngr0,2)
c
      write(*,*)'Job type (3..density, 4..potential):'
      read(*,*)nd  
      write(*,*)'Formatted input (Y/N)?'
      read(*,'(a)')ok
      write(*,*)'Density file:'
      read(*,'(a)')filename
      if(nd.eq.4)write(*,*)'Potential file:'
      if(nd.eq.4)read(*,'(a)')fpot
      write(*,*)'Number of atoms:'
      read(*,*)nat
      open(33,file='BBB.SCR',form='unformatted')
      write(*,*)'Units in the input file (B/A):'
      read(*,'(a)')ok1
      bohr=0.529177d0 
      if(ok1.eq.'b'.or.ok1.eq.'B')bohr=1.0d0
c
c     formatted gaussian output:
      if(ok.eq.'y'.or.ok.eq.'Y')then
       open(2,file=filename,status='old')
       read(2,'(a40)')n1
       write(*,*)n1
       read(2,'(a40)')n2
       read(2,*)ifl,(x0(i),i=1,3)
       read(2,*)nx,(x1(i),i=1,3)
       read(2,*)ny,(x2(i),i=1,3)
       read(2,*)nz,(x3(i),i=1,3)
       do 101 i=1,3
       x0(i)=x0(i)/bohr
       x1(i)=x1(i)/bohr
       x2(i)=x2(i)/bohr
101    x3(i)=x3(i)/bohr
       do 4 i=1,nat
4      read(2,*)idum,q(i),(r(j,i),j=1,3)
c
c      when individual orbitals, extra line skipped:
       read(2,'(a)')ppp
       do 2223 i=1,80
       if(ppp(i:i).eq.'.')then
        write(6,*)'no line skipped after atoms'
        backspace 2
        goto 2224
       endif
2223   continue
       write(6,*)'one line skipped after atoms'
2224   continue
c
       if(nx.ne.ny.or.nx.ne.nz)write(6,*)'nx,ny,nz are not equal'
       if(nx.gt.ngr0.or.ny.gt.ngr0.or.nz.gt.ngr0)then
        write(*,*)'Too many points'
        goto 1
       endif
       dx=sqrt(x1(1)**2+x1(2)**2+x1(3)**2)
       dy=sqrt(x2(1)**2+x2(2)**2+x2(3)**2)
       dz=sqrt(x3(1)**2+x3(2)**2+x3(3)**2)
       if(dx.ne.dy.or.dx.ne.dz)then
        write(6,*)'dx,dy,dz are not equal - fatal error'
        goto 1
       endif 
       write(33)n1,n2,nd,nx-1,ny-1,nz-1,(x0(i),i=1,3),dx,dy,dz
       if(nd.eq.4)then
        open(4,file=fpot,status='old')
        do 20 i=1,nat+6
20      read(4,*)
       endif
c
       do 2 k=1,nx
       write(*,*)k
       do 3 j=1,ny
       if(nd.eq.4)read(4,'(6E13.5)')(rofo(i,j,2),i=1,nz)
3      read(2,'(6E13.5)')(rofo(i,j,1),i=1,nz)
       if(nd.eq.4)then
        write(33)(((rofo(i,j,l),l=1,2),i=1,nz),j=1,ny)
       else
        write(33)((rofo(i,j,1),i=1,nz),j=1,ny)
       endif
2      continue
c
c     unformatted gaussian output:
      else
       open(2,file=filename,status='old',form='unformatted')
       read(2)n1
       read(2)n2
       read(2)ifl,(x0(i),i=1,3)
       read(2)nx,(x1(i),i=1,3)
       read(2)ny,(x2(i),i=1,3)
       read(2)nz,(x3(i),i=1,3)
       do 100 i=1,3
       x0(i)=x0(i)/bohr
       x1(i)=x1(i)/bohr
       x2(i)=x2(i)/bohr
100    x3(i)=x3(i)/bohr
       do 6 i=1,nat
6      read(2)idum,q(i),(r(j,i),j=1,3)
       if(nx.ne.ny.or.nx.ne.nz)write(6,*)'nx,ny,nz are not equal'
       if(nx.gt.ngr0.or.ny.gt.ngr0.or.nz.gt.ngr0)then
        write(*,*)'Too many points'
        goto 1
       endif
       dx=sqrt(x1(1)**2+x1(2)**2+x1(3)**2)
       dy=sqrt(x2(1)**2+x2(2)**2+x2(3)**2)
       dz=sqrt(x3(1)**2+x3(2)**2+x3(3)**2)
       if(dx.ne.dy.or.dx.ne.dz)then
        write(6,*)'dx,dy,dz are not equal - fatal error'
        goto 1
       endif 
       write(33)n1,n2,nd,nx-1,ny-1,nz-1,(x0(i),i=1,3),dx,dy,dz
       if(nd.eq.4)then
        open(4,file=fpot,status='old')
        do 21 i=1,nat+6
21      read(4,*)
       endif
       do 7 k=1,nx
       do 8 j=1,ny
       if(nd.eq.4)read(4)(rofo(i,j,2),i=1,nz)
8      read(2)(rofo(i,j,1),i=1,nz)
       if(nd.eq.4)then
        write(33)(((rofo(i,j,l),l=1,2),i=1,nz),j=1,ny)
       else
        write(33)((rofo(i,j,1),i=1,nz),j=1,ny)
       endif
7      continue
      endif
c
      do 10 i=1,nat
      do 10 j=1,3
10    r(j,i)=r(j,i)/bohr
c
      do 11 i=1,nat
      do 12 j=1,3
12    rt(j,i)=r(j,i)-x0(j)
      x=rt(1,i)*x1(1)+rt(2,i)*x1(2)+rt(3,i)*x1(3)
      y=rt(1,i)*x2(1)+rt(2,i)*x2(2)+rt(3,i)*x2(3)
      z=rt(1,i)*x3(1)+rt(2,i)*x3(2)+rt(3,i)*x3(3)
      rt(1,i)=x/dx
      rt(2,i)=y/dy
11    rt(3,i)=z/dz
c
      write(33)nat,((rt(i,l),i=1,3),l=1,nat),(q(i),i=1,nat)
c
      write(*,*)'Outputfile BBB.SCR written.'
1     close(2)
      if(nd.eq.4)close(4)
      close(33)
      stop
      end