program cubcut2
      implicit none
      integer*4 i,nd(3),nc,ix,iy,iz,N,nats,norb,iorb,jo,iu
      real*8 r0(3),dd,v(3),dr(3),s,vn,x(3),bohr,xx,yy,zz
      character*80 fn
      real*8,allocatable::rh(:),rhz(:)
      integer,allocatable::iol(:)
c
      write(6,6000)
6000  format(/,'     cut through Gaussian orbital .cub file',/)
 
      write(6,*)'Filename:'
      read(5,'(a)')fn
      open(9,file=fn)
      read(9,*)
      read(9,*)
      read(9,*)nats,(r0(i),i=1,3)
      if(nats.lt.0)then
       iu=0
       nats=-nats
      else
       iu=1
      endif
      read(9,*)nd(1),dr(1)
      read(9,*)nd(2),dr(2),dr(2)
      read(9,*)nd(3),dr(3),dr(3),dr(3)
      nc=nd(1)*nd(2)*nd(3)
      allocate(rh(nc))
      do 2 i=1,nats
2     read(9,*)
      read(9,*)norb
      allocate(iol(norb))
      backspace 9
      read(9,910)norb,(iol(i),i=1,norb)
910   format(10i5)
      write(6,911)norb
911   format(i5,' orbitals:  ')
      write(6,910)(iol(i),i=1,norb)
      write(6,912)
912   format(' Which one you want: ',$)
      read(5,*)iorb
      jo=0
      do 7 i=1,norb
7     if(iol(i).eq.iorb)jo=i
      allocate(rhz(nd(3)*norb))
      do 3 ix=1,nd(1)
      do 3 iy=1,nd(2)
      read(9,4000)((rhz(iz+nd(3)*(i-1)),i=1,norb),iz=1,nd(3))
      do 31 iz=1,nd(3)
31    rh(iz+(iy-1)*nd(3)+(ix-1)*nd(3)*nd(2))=rhz(iz+nd(3)*(jo-1))
3     continue
4000  format(6E13.5)
      close(9)
      if(iu.eq.0)then
       bohr=0.529177d0
       do 1 ix=1,3
       dr(ix)=dr(ix)*bohr
1      r0(ix)=r0(ix)*bohr
      endif
      write(6,*)'cube read in'
      write(6,600)r0(1),r0(1)+dble(nd(1))*dr(1),
     1r0(2),r0(2)+dble(nd(2))*dr(2),
     1r0(3),r0(3)+dble(nd(3))*dr(3),nats,nd(1),nd(2),nd(3),
     1dr(1),dr(2),dr(3)
600   format(' xmin ',f12.3,' xmax ',f12.4,/,
     1       ' ymin ',f12.3,' ymax ',f12.4,/,
     1       ' zmin ',f12.3,' zmax ',f12.4,/,
     1       i4,' atoms',/,
     1       i3,' x ',i3,' x ',i3,/,
     1       ' dx: ',f8.3,' dy: ',f8.3,' dz: ',f8.3,/)

      write(6,*)'Initial point x y z (A):'
      read(5,*)x(1),x(2),x(3)
      write(6,*)'Direction vector vx vy vz (A):'
      read(5,*)v(1),v(2),v(3)
      write(6,*)'Increment (A), number of points:'
      read(5,*)dd,N

      vn=dsqrt(v(1)**2+v(2)**2+v(3)**2)
      do 4 ix=1,3
4     v(ix)=v(ix)/vn

      open(9,file='X.PRN')
      do 6 i=1,N
      xx=x(1)+dd*v(1)*dble(i-1)
      yy=x(2)+dd*v(2)*dble(i-1)
      zz=x(3)+dd*v(3)*dble(i-1)
      ix=NINT((xx-r0(1))/dr(1))
      iy=NINT((yy-r0(2))/dr(2))
      iz=NINT((zz-r0(3))/dr(3))
      if(ix.ge.1.and.ix.le.nd(1).and.
     1   iy.ge.1.and.iy.le.nd(2).and.
     1   iz.ge.1.and.iz.le.nd(3))then
       s=rh(iz+(iy-1)*nd(3)+(ix-1)*nd(3)*nd(2))
      else
       s=0.0d0
      endif
6     write(9,900)dd*dble(i-1),s
900   format(f8.3,g15.5)
      close(9)
      write(6,*)'X.PRN written'

      end