      PROGRAM STFTOCUBE
*
*     interface for Molekel http://www.cscs.ch/molekel
      IMPLICIT NONE

      INTEGER maxatoms,maxbasin,maxcell
      INTEGER i,j,k,natoms,truelen,katoms,xatt,yatt,zatt
      integer attract_label
      integer npx,npy,npz,npxyz,natt
      integer, dimension (:), allocatable :: nat
      integer, dimension (:,:), allocatable :: index_basin
      double precision, dimension (:), allocatable :: x,y,z,volume
      double precision, dimension (:), allocatable :: xat,yat,zat
      double precision xmin,xmax,ymin,ymax,zmin,zmax
      double precision stepx,stepy,stepz,stepnull
      integer ix,iy,iz,isx,isy,isz,indy,ij
      integer*2, dimension (:), allocatable :: attract_code
      real, dimension (:), allocatable :: elf
      real rnat
      character*3 file_type
      character*16, dimension (:), allocatable :: basin_name
      character*4, dimension (:), allocatable :: atom_name
      character*4 atom_namet
      CHARACTER*40 filein,filein2,fileout,dunname 
      character*1 molecule(40),blank
      logical if_file
      equivalence (molecule(1),dunname)
      data blank/' '/ 
      WRITE(*,'(/,a,/)')'SBF_TO_CUBE 1.0 - Laurent Joubert'

      WRITE(*,'(a)')'WFN INPUT FILE ?'
      READ(*,'(a40)') filein
      inquire(file=filein,exist=if_file)
      if(.not.if_file) then
        write(*,'(2x,"requested file does not exist:",a40)') filein
        stop
      endif
      OPEN(1,file=filein,status='old')

      READ(1,'(/,61x,i3)')natoms
      allocate(xat(natoms),yat(natoms),zat(natoms),nat(natoms))
      allocate(atom_name(natoms))
      DO i=1,natoms
         READ(1,'(a4,20x,3f12.8,11x,f4.1)') atom_name(i),xat(i),
     .        yat(i),zat(i),rnat
         nat(i)=int(rnat)
      ENDDO
      close(1)

      WRITE(*,'(/,a)')'SBF FILE ?'
      READ(*,'(a40)') filein2
      inquire(file=filein2,exist=if_file)
      if(.not.if_file) then
        write(*,'(2x,"requested file does not exist:",a40)') filein2
        stop
      endif
      dunname=filein2
      i=1
      do while (molecule(i).ne.blank)
         i=i+1
      enddo
      file_type=molecule(i-7)//molecule(i-6)//molecule(i-5)
*     WRITE(*,'(a)')'TYPE OF FILE : elf(rho)/syn/bas ?'
*     READ(*,'(a3)') file_type

      OPEN(2,file=filein2,status='old',form='unformatted')
      READ(2) npx,npy,npz
      READ(2) xmin,xmax,ymin,ymax,zmin,zmax
      npxyz=npx*npy*npz
      if((file_type.eq.'elf').or.(file_type.eq.'rho').or.
     .                           (file_type.eq.'lap')) then
         allocate(elf(npxyz))
         read(2)(elf(i),i=1,npxyz)
         write(*,*) npxyz, ' values read'
      else if(file_type.eq.'bas') then
         read(2) katoms,attract_label
         allocate(index_basin(attract_label,5),
     .            basin_name(attract_label))
         allocate(volume(attract_label),x(attract_label),
     .             y(attract_label))
         allocate(z(attract_label))
         do i=1,katoms
            read(2) natt,atom_namet,xatt,yatt,zatt
         enddo
         do i=1,attract_label
            read(2) index_basin(i,1),basin_name(i),volume(i),x(i),y(i),
     .              z(i)
         enddo
         allocate(attract_code(npxyz))
         read(2)(attract_code(i),i=1,npxyz)
      else if(file_type.eq.'syn') then
         allocate(attract_code(npxyz))
         read(2)(attract_code(i),i=1,npxyz)
      endif
      close(2)
***
*     read steps for x, y, z
***
      write(*,*) ' steps for x, y z (defaults 1 1 1)'
       isx = 1
       isy = 1
       isz = 1
      read(*,*) isx,isy,isz
      if((file_type.eq.'elf').or.(file_type.eq.'rho').or.
     .   (file_type.eq.'lap')) then
         ij=0
         ix=0
         do i=1,npx,isx
            ix=ix+1
            iy=0
            do j=1,npy,isy
               iy=iy+1
               iz=0
               do k=1,npz,isz
                  iz=iz+1
                  ij=ij+1
                  indy=k+(j-1)*npz+(i-1)*npz*npy
                  elf(ij)=elf(indy)
               enddo
            enddo
         enddo
         npx=ix
         npy=iy
         npz=iz
         npxyz=npx*npy*npz
      endif
      if((file_type.eq.'bas').or.(file_type.eq.'syn')) then
        ij=0
         ix=0
         do i=1,npx,isx
            ix=ix+1
            iy=0
            do j=1,npy,isy
               iy=iy+1
               iz=0
               do k=1,npz,isz
                  iz=iz+1
                  ij=ij+1
                  indy=k+(j-1)*npz+(i-1)*npz*npy
                  attract_code(ij)=attract_code(indy)
               enddo
            enddo
         enddo
         npx=ix
         npy=iy
         npz=iz
         npxyz=npx*npy*npz
      endif
      
      fileout=filein2(1:truelen(filein2))//'.cube'
      OPEN(3,file=fileout,status='unknown')  
      WRITE(3,*)' Cube File'
      WRITE(3,*)' ELF grid'
      WRITE(3,'(i5,3f12.6)')natoms,xmin,ymin,zmin
      stepx=(xmax-xmin)/(npx-1)
      stepy=(ymax-ymin)/(npy-1)
      stepz=(zmax-zmin)/(npz-1)
      stepnull=0d0
      WRITE(3,'(i5,3f12.6)')npx,stepx,stepnull,stepnull
      WRITE(3,'(i5,3f12.6)')npy,stepnull,stepy,stepnull
      WRITE(3,'(i5,3f12.6)')npz,stepnull,stepnull,stepz
      DO i=1,natoms
         WRITE(3,'(i5,4f12.6)') nat(i),float(nat(i)),
     .         xat(i),yat(i),zat(i)
      ENDDO
C HSR added this to get the sign convention for Lap the same as dgrid
      IF (file_type.eq.'lap') THEN
      DO i=1,npxyz
      elf(i) = -elf(i)
      ENDDO
      ENDIF
      IF ((file_type.eq.'elf').or.(file_type.eq.'rho').or.
     .    (file_type.eq.'lap')) THEN
         WRITE(3,'(6(1pe13.5))')(elf(i),i=1,npxyz)
      ELSE IF (file_type.eq.'bas') THEN
         WRITE(3,'(20i4)') (attract_code(i),i=1,npxyz)
      ELSE IF (file_type.eq.'syn') THEN 
         WRITE(3,'(40i2)') (attract_code(i),i=1,npxyz)      
      ENDIF
      WRITE(*,'(/,a)')'OUTPUT WRITTEN : '//fileout 
 
      CLOSE(1)
      CLOSE(2)

      END
      
*     ********************************************************

      FUNCTION TRUELEN(word)

      CHARACTER*30 word
      INTEGER TRUELEN

      TRUELEN=index(word,'.')-1
      IF (truelen.lt.0) truelen=len(word)

      END

*     *********************************************************
