program ecda
      implicit none
      integer*4 nes,ne,na,nb,i,nas,nat,j,iargc
      real*8 debye,n(3),r,rp,ds,c,D(3),p(3),l(3)
      real*8,allocatable::m(:,:),u(:,:),e(:),v(:,:),x(:,:),z(:)
      character*80 fn
      debye=2.541765d0
      c=137.5d0
      if(iargc().ne.1)then
       write(6,600)
600    format(' Usage: ecda <Gaussian output>')
       stop
      endif
      call getarg(1,fn)
      ne=nes(fn)
      nat=nas(fn)
      allocate(e(ne),m(ne,3),u(ne,3),v(ne,3),x(nat,3),z(nat))
      call rtm(fn,ne,e,m,u,v,D,na,nb,x,z,nat)
      open(9,file='ECDA.TAB')
      write(9,90)
90    format('  transition l / nm   D          R           R',/,
     1       '                                LORG      GENERAL',/,
     160(1h-))
c     nuclear electric dipole moment:
      p=0.0d0
      do 2 i=1,3
      do 2 j=1,nat
 2    p(i)=p(i)+z(j)*x(j,i)
c     electronic dipole moment:
      l=D-p
      do 1 i=1,ne
      ds=u(i,1)*u(i,1)+u(i,2)*u(i,2)+u(i,3)*u(i,3)
c     (1/2) r x grad  . r:
      r=u(i,1)*m(i,1)+u(i,2)*m(i,2)+u(i,3)*m(i,3)
c     (1/2) r_nn x grad /Ne:
      n(1)=(l(2)*v(i,3)-l(3)*v(i,2))/dble(na+nb)/2.0d0
      n(2)=(l(3)*v(i,1)-l(1)*v(i,3))/dble(na+nb)/2.0d0
      n(3)=(l(1)*v(i,2)-l(2)*v(i,1))/dble(na+nb)/2.0d0
c     (1/2) [r x grad - r_nn x grad /Ne]  .r:
      rp=r-(u(i,1)*n(1)+u(i,2)*n(2)+u(i,3)*n(3))
      r =r *debye**2/c
      rp=rp*debye**2/c
      write(6,900)i,e(i),ds*debye**2,rp,r
1     write(9,900)i,e(i),ds*debye**2,rp,r
900   format(i6,f12.2,3e12.4)
      write(9,91)
91    format(60(1h-))
      close(9)
      end

      subroutine rtm(fn,ne,e,m,u,v,D,na,nb,x,z,nat)
      implicit none
      integer*4 ne,i,j,ir,na,nb,nat
      real*8 m(ne,3),u(ne,3),e(ne),v(ne,3),D(3),debye,x(nat,3),z(nat),
     1bohr
      character*(*) fn
      character*80 s80
      debye=2.541765d0
      bohr=0.529177d0
      u=0.0d0
      v=0.0d0
      m=0.0d0
      ir=0
      na=0
      nb=0
      open(8,file=fn)
1     read(8,80,end=99,err=99)s80
80    format(a80)
      if(s80(8:22).eq.'alpha electrons')then
       read(s80(1:6),*)na
       read(s80(25:31),*)nb
       write(6,601)'a'
      endif
      if(s80(26:65).eq.'transition electric dipole moments (Au):')then
       read(8,*)
       do 2 i=1,ne
2      read(8,*)u(i,1),(u(i,j),j=1,3)
       write(6,601)'u'
601    format(a1,$)
      endif
      if(s80(26:65).eq.'transition velocity dipole moments (Au):')then
       read(8,*)
       do 3 i=1,ne
3      read(8,*)v(i,1),(v(i,j),j=1,3)
       write(6,601)'v'
      endif
      if(s80(26:65).eq.'transition magnetic dipole moments (Au):')then
       read(8,*)
       do 4 i=1,ne
4      read(8,*)m(i,1),(m(i,j),j=1,3)
       write(6,601)'m'
      endif
      do 5 i=1,len(s80)-3
      if(s80(2:3).eq.'Ex'.and.s80(i:i+3).eq.' nm ')then
       ir=ir+1
       read(s80(i-8:i-1),*)e(ir)
       if(ir.eq.ne)write(6,601)'e'
      endif
5     continue
      if(s80(2:42).eq.'Dipole moment (field-independent basis, D')then
       read(8,80,end=99,err=99)s80
       read(s80( 7:26),*)D(1)
       read(s80(33:52),*)D(2)
       read(s80(59:78),*)D(3)
       if(ir.eq.ne)write(6,601)'D'
       D=D/debye
      endif
      if(s80(27:44).eq.'Input orientation:')then
       read(8,*)
       read(8,*)
       read(8,*)
       read(8,*)
       do 7 i=1,nat
7      read(8,*)z(i),z(i),(x(i,j),j=1,3)
       write(6,601)'X'
       x=x/bohr
      endif
      goto 1
99    close(8)
      do 6 i=1,ne
      do 6 j=1,3
      v(i,j)=-v(i,j)
 6    m(i,j)=-0.5d0*m(i,j)
      write(6,*)
      return
      end

      function nes(fn)
      implicit none
      integer*4 nes,ne
      character*(*) fn
      character*80 s80
      ne=0
      open(8,file=fn)
1     read(8,80,end=99,err=99)s80
80    format(a80)
      if(s80(26:65).eq.'transition electric dipole moments (Au):')then
2      read(8,80,end=99,err=99)s80
       if(s80(2:35).eq.'Ground to excited state transition')goto 99
       ne=ne+1
       goto 2
      endif
      goto 1
99    close(8)
      nes=ne-1
      write(6,*)ne-1,' transitions'
      return
      end

      function nas(fn)
      implicit none
      integer*4 nas,ne
      character*(*) fn
      character*80 s80
      ne=0
      open(8,file=fn)
1     read(8,80,end=99,err=99)s80
80    format(a80)
      if(s80(27:44).eq.'Input orientation:')then
       read(8,80,end=99,err=99)
       read(8,80,end=99,err=99)
       read(8,80,end=99,err=99)
       read(8,80,end=99,err=99)
2      read(8,80,end=99,err=99)s80
       if(s80(2:3).eq.'--')goto 99
       ne=ne+1
       goto 2
      endif
      goto 1
99    close(8)
      nas=ne
      write(6,*)ne,' atoms'
      return
      end