program transro
      implicit none
      integer*4 cmax
      parameter (cmax=19)
      integer*4 i,ifa,nc(3),aij(cmax,2),daij(cmax,2),ni,nid,
     1nat,N,i2,iargc,j,Nexc,ie,ief
      real*8 c(cmax),d(cmax),r0(3),rc(3,3),cij,csum
      integer*4,allocatable::q(:)
      real*8,allocatable::r(:),o(:),oa(:),ob(:),zq(:)
      character*80 f,s80
      character*10 s10
      logical lopen
c 
      write(6,602)
602   format(/,' Transition charge density from Gaussian cub files',/,/,
     1         '        Input:    <Gaussian TD output> ',/,
     1         '                  .cub orbital cub files ',/,
     2         '        output:   trans.cub',/)
      if(iargc().ne.2)call report('Usage: transro <output> <Nstate>')
      call getarg(1,f)
      call getarg(2,s80)
      read(s80,*)Nexc
      call rdf(f,cmax,Nexc,c,d,aij,daij,ni,nid,lopen)

c     just look what is available:
      write(6,*)'Involved orbitals:'
      ifa=0
      do 1 i=1,ni
      do 1 i2=1,2
      write(6,630)aij(i,i2)
 630  format(i6,$)
      write(s10,100)aij(i,i2)
 100  format(i10)
      ie=ief(s10)
      if(lopen)then
       call inq(s10(ie:len(s10))//'A.cub',ifa,nat,r0,rc,nc)
      else
       call inq(s10(ie:len(s10))//'.cub',ifa,nat,r0,rc,nc)
      endif
 1    continue
      if(lopen)then
       do 2 i=1,nid
       do 2 i2=1,2
       write(s10,100)daij(i,i2)
       ie=ief(s10)
 2     call inq(s10(ie:len(s10))//'B.cub',ifa,nat,r0,rc,nc)
      endif

c     read cubes and make the sum:
      csum=0.0d0
      N=nc(1)*nc(2)*nc(3)
      allocate(o(N),oa(N),ob(N),zq(nat),q(nat),r(3*nat))
      o=0.0d0
      do 3 i=1,ni
      write(s10,100)aij(i,1)
      ie=ief(s10)
      if(lopen)then
       call rcu(s10(ie:len(s10))//'A.cub',N,nat,r0,rc,oa,r,q,zq,nc)
      else
       call rcu(s10(ie:len(s10))//'.cub',N,nat,r0,rc,oa,r,q,zq,nc)
      endif
      write(s10,100)aij(i,2)
      ie=ief(s10)
      if(lopen)then
       call rcu(s10(ie:len(s10))//'A.cub',N,nat,r0,rc,ob,r,q,zq,nc)
      else
       call rcu(s10(ie:len(s10))//'.cub',N,nat,r0,rc,ob,r,q,zq,nc)
      endif
      cij=c(i)
      csum=csum+cij**2
      write(6,600)i,aij(i,1),aij(i,2),cij
 600  format(i6,i6,' ->',i6,f12.5)
      do 3 j=1,N
 3    o(j)=o(j)+cij*oa(j)*ob(j)
      write(6,*)o(1),o(2)

      if(lopen)then
       do 4 i=1,nid
       write(s10,100)aij(i,1)
       ie=ief(s10)
       call rcu(s10(ie:len(s10))//'B.cub',N,nat,r0,rc,oa,r,q,zq,nc)
       write(s10,100)aij(i,2)
       ie=ief(s10)
       call rcu(s10(ie:len(s10))//'B.cub',N,nat,r0,rc,ob,r,q,zq,nc)
       cij=d(i)
       csum=csum+cij**2
       do 4 j=1,N
 4     o(j)=o(j)+cij*oa(j)*ob(j)
      endif

      write(6,601)csum
 601  format(' csum:',f12.5)

      call wh('r.cub',nat,r0,rc,o,q,zq,r,nc)
      end

      subroutine addc(ta,tb,tc,ni,cmax,aij,cij)
      implicit none
      integer*4 i,ta,tb,ni,cmax,aij(cmax,2)
      real*8 cij(*),tc
c     add to buffer, if many, take only the biggest ones
      if(ni+1.le.cmax)then
       ni=ni+1
       aij(ni,1)=ta
       aij(ni,2)=tb
       cij(ni)=tc
      else
       do 1 i=1,cmax
       if(dabs(cij(i)).lt.dabs(tc))then
        aij(i,1)=ta
        aij(i,2)=tb
        cij(i)=tc
        return
       endif
1      continue
      endif
      return
      end

      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end

      subroutine rdf(f,cmax,Nexc,cij,dij,aij,daij,ni,nid,lopen)
      implicit none
      integer*4 cmax,il,ie,Nexc,aij(cmax,2),daij(cmax,2),ni,nid,
     1ta,tb
      character*(*) f
      character*80 s80
      real*8 tc,cij(*),dij(*)
      logical lopen
      character*1 o
      lopen=.false.
      open(9,file=f,status='old')
      il=0
      ie=0
 1    read(9,80,end=99,err=99)s80
 80   format(a80)
      il=il+1
      if(s80(2:15).eq.'Excited State '.and.s80(16:16).ne.'s')then
       ie=ie+1
       write(6,80)s80
       if(ie.eq.Nexc)then
        write(6,*)'Required state found'
        ni=0
        nid=0
111     read(9,80,end=99,err=99)s80
        if(s80(10:11).eq.'->')then
         o=s80(8:8)
         if(o.eq.'A'.or.o.eq.'B')then
          lopen=.true.
c         open shell, alpha:
          read(s80( 1: 7),*)ta
          if(s80(15:15).eq.o)then
           read(s80(12:14),*)tb
          else
           if(s80(16:16).eq.o)then
            read(s80(12:15),*)tb
           else
            write(6,80)s80
            call report('Unknown format of excitations')
           endif
          endif
          read(s80(18:30),*)tc
          if(o.eq.'A')call addc(ta,tb,tc,ni ,cmax, aij,cij)
          if(o.eq.'B')call addc(ta,tb,tc,nid,cmax,daij,dij)
         else
c         closed shell:
          read(s80( 1: 8),*)ta     
          read(s80(12:15),*)tb  
          read(s80(18:30),*)tc
          call addc(ta,tb,tc,ni,cmax,aij,cij)
         endif
         goto 111
        endif
        goto 99
       endif
c      if(ie.eq.Nexc)then
      endif
      goto 1
 99   close(9)
      write(6,605)il,ni,nid
605   format(i6,' lines,',i3,' A orbitals,',i3,' B orbitals')
      if(ie.ne.Nexc)call report('Required state not found')
      return
      end


      subroutine rh(f,nat,r0,rc,nc)
c     read header of gaussian cub file:
      implicit none
      real*8 r0(3),rc(3,3)
      integer*4 nc(3),i,nat
      character*(*)f
      open(9,file=f,status='old')
      read(9,*)
      read(9,*)
      read(9,*)nat,r0(1),r0(2),r0(3)
      nat=abs(nat)
      read(9,*)nc(1),(rc(1,i),i=1,3)
      read(9,*)nc(2),(rc(2,i),i=1,3)
      read(9,*)nc(3),(rc(3,i),i=1,3)
      write(6,600)nc(1),nc(2),nc(3)
600   format(i6,' x',i6,' x',i6,' grid points')
      close(9)
      write(6,*)f//' header read'
      return
      end

      subroutine rcu(f,N,nat,r0,rc,o,r,q,zq,nc)
c     read gaussian cub file:
      implicit none
      integer*4 N,nc(3),i,i1,i2,i3,q(*),ix,nat
      real*8 r0(3),rc(3,3),o(N),r(*),zq(*)
      character*(*)f
      real*8,allocatable::fi(:)
      logical lex
      inquire(file=f,exist=lex)
      if(lex)then
       open(9,file=f,status='old')
       read(9,*)
       read(9,*)
       read(9,*)nat,r0(1),r0(2),r0(3)
       nat=abs(nat)
       read(9,*)nc(1),(rc(1,i),i=1,3)
       read(9,*)nc(2),(rc(2,i),i=1,3)
       read(9,*)nc(3),(rc(3,i),i=1,3)
       do 5 i=1,nat
5      read(9,*)q(i),zq(i),(r(ix+3*(i-1)),ix=1,3)
       read(9,*)
       allocate(fi(nc(3)))
       do 6 i1=1,nc(1)
       do 6 i2=1,nc(2)
       read(9,4000)(fi(i3),i3=1,nc(3))
4000   format(6E13.5)
       do 6 i3=1,nc(3)
6      o(i1+(i2-1)*nc(1)+(i3-1)*nc(1)*nc(2))=fi(i3)
       close(9)
      else
       o=0.0d0
       write(6,*)f//' skipped'
      endif
      return
      end

      subroutine wh(f,nat,r0,rc,o,q,zq,r,nc)
c     write gaussian cub file:
      implicit none
      real*8 r0(3),rc(3,3),o(*),r(*),zq(*)
      real*8,allocatable::fi(:)
      integer*4 i,i1,i2,i3,nc(*),ix,nat,q(*)
      character*(*)f
      open(9,file=f)
      write(9,91)
 91   format(' test MO=1')
      write(9,92)
 92   format(' MO coefficients')
      write(9,93)-nat,r0(1),r0(2),r0(3)
      write(9,93)nc(1),(rc(1,i),i=1,3)
      write(9,93)nc(2),(rc(2,i),i=1,3)
      write(9,93)nc(3),(rc(3,i),i=1,3)
 93   format(i5,4f12.6)
      do 5 i=1,nat
5     write(9,93)q(i),zq(i),(r(ix+3*(i-1)),ix=1,3)
      write(9,94)1,1
 94   format(2i5)
      allocate(fi(nc(3)))
      do 6 i1=1,nc(1)
      do 6 i2=1,nc(2)
      do 61 i3=1,nc(3)
61    fi(i3)=o(i1+(i2-1)*nc(1)+(i3-1)*nc(1)*nc(2))
6     write(9,4000)(fi(i3),i3=1,nc(3))
4000  format(6E13.5)
      close(9)
      write(6,*)f//' written'
      return
      end

      function ief(s10)
      implicit none
      character*(*) s10
      integer*4 ief,ie
      do 5 ie=1,len(s10)
5     if(s10(ie:ie).ne.' ')goto 6
6     ief=ie
      return
      end

      subroutine inq(s,ifa,nat,r0,rc,nc)
      implicit none
      integer*4 ifa,nat,nc(3)
      real*8 r0(3),rc(3,3)
      character*(*) s
      logical lex
      inquire(file=s,exist=lex)
      if(lex)then
       write(6,*)s//' found'
       ifa=ifa+1
       if(ifa.eq.1)call  rh(s,nat,r0,rc,nc)
      else
       write(6,*)s//' not found'
      endif
      return
      end