PROGRAM mcdv
c     vivrationally resolved MCD
      implicit none
      character*80 fn,s80
      logical lex,HW,ltab,LQ1,LDD
      integer*4 nat,NG,ip(10),ix,jx,IERR,iw,N,LEXCL,NExc,iex,
     1N7,NSL,ic,icen,NQ1,II,np0,i,kk,ik,nuk,
     1iwra,iwrb,iwr,nt,NRoot
      integer*4,allocatable::ty(:),si(:),sit(:)
      REAL*8 CM,GETDET,sp1,sp2,fac,sp,EJ,THR,E0,FC,pp,ut(3),u0(3),
     1proc,ds,qk,pm
      real*8,allocatable::r(:),J(:,:),K(:),A(:,:),B(:),C(:,:),D(:),
     1E(:,:),W(:),WP(:),G(:,:),GP(:,:),T(:,:),F(:,:),TP(:,:),
     1JT(:,:),EE(:,:),FI(:,:),GH(:,:),GHP(:,:),MYA(:,:),TV(:),TU(:),
     1MYB(:),dd(:),ddi(:)
      CM=219474.0d0
      E0=42256.0d0
c     gaussian output filename:
      fn='G.OUT'
c     maximal degree of excitations ("class"):
      LEXCL=3
c     use first dipole derivatives:
      LQ1=.true.
c     extensive output:
      HW=.true.
c     dimension of the temporary list:
      np0=10
c     optional writing option:
      iwr=1
c     write matrix A:
      iwra=0
c     write vector B:
      iwrb=0
c     write FC.TBA:
      ltab=.true.
c     threshold to consider Franck-Condon factors:
      THR=1.0d-9
c     read dipole derivatives:
      LDD=.true.
c     number of th eexcited state
      Nroot=1
      inquire(file='MCDV.OPT',exist=lex)
      if(lex)then
       open(9,file='MCDV.OPT')
1      read(9,900,err=99,end=99)s80
900    format(a80)
       if(s80(1:2).eq.'HW')read(9,*)HW
       if(s80(1:4).eq.'FILE')read(9,900)fn
       if(s80(1:4).eq.'LEXC')read(9,*)LEXCL
       if(s80(1:4).eq.'IWRA')read(9,*)iwra
       if(s80(1:4).eq.'IWRB')read(9,*)iwrb
       if(s80(1:4).eq.'IWRT')read(9,*)iwr
       if(s80(1:3).eq.'NP0' )read(9,*)np0
       if(s80(1:3).eq.'LQ1' )read(9,*)LQ1
       if(s80(1:3).eq.'LDD' )read(9,*)LDD
       if(s80(1:4).eq.'LTAB')read(9,*)ltab
       if(s80(1:4).eq.'THRE')read(9,*)THR
       if(s80(1:4).eq.'NROO')read(9,*)NRoot
       goto 1
99     close(9)
      endif

      call rdx(fn,nat,0,ty,r,NG,.true.)
      if(nat.eq.0)call report('Nat = 0')
      allocate(r(3*nat),ty(nat))
      call rdx(fn,nat,1,ty,r,NG,.true.)
      write(6,*)nat,' atoms, ',NG,' geomnetries'
      N=3*nat-6

      if(LDD)then
       allocate(dd(9*nat),ddi(9*nat))
       call rdd(fn,dd,nat,Nroot,u0,iwr)
       call trafo(dd,ddi,nat,N)
      endif

      allocate(J(N,N),K(N),A(N,N),B(N),C(N,N),D(N),E(N,N),W(N),WP(N),
     1JT(N,N))
      call zm(J,N)
      call zm(JT,N)
      call zm(A,N)
      call zm(C,N)
      call zm(E,N)
      call zv(K,N)
      call zv(B,N)
      call zv(D,N)
      call zv(W,N)
      call zv(WP,N)
      do 3 ix=1,10
3     ip(ix)=0
      iw=0
      open(9,file=fn,status='old')
2     read(9,900,end=88,err=88)s80
      if(s80(4:20).eq.'Duschinsky Matrix')call rdm(9,N,JT,ip(1))
      if(s80(4:15).eq.'Shift Vector')call rdv(9,N,K,ip(2))
      if(s80(4:11).eq.'A Matrix')call rdm(9,N,A,ip(3))
      if(s80(4:11).eq.'B Vector')call rdv(9,N,B,ip(4))
      if(s80(4:11).eq.'C Matrix')call rdm(9,N,C,ip(5))
      if(s80(4:11).eq.'D Vector')call rdv(9,N,D,ip(6))
      if(s80(4:11).eq.'E Matrix')call rdm(9,N,E,ip(7))
      if(s80(2:22).eq.'Harmonic frequencies ')call rw(9,N,W,WP,iw)
      goto 2
88    close(9)
      if(ip(1).eq.0)call report('Duschinsky Matrix not found')
      if(ip(2).eq.0)call report('Shift Vector not found')
      if(ip(3).eq.0)call report('A Matrix not found')
      if(ip(4).eq.0)call report('B Vector not found')
      if(ip(5).eq.0)call report('C Matrix not found')
      if(ip(6).eq.0)call report('D Vector not found')
      if(ip(7).eq.0)call report('E Matrix not found')
      if(iw.ne.2)call report('Frequencies not found')
      allocate(F(N,N),G(N,N),GP(N,N),T(N,N),TP(N,N),FI(N,N),EE(N,2*N))
      allocate(GH(N,N),GHP(N,N),MYA(N,N),TV(N),TU(N),MYB(N))
      call mz(G,N)
      call mz(GP,N)
      do 4 ix=1,N
      G(ix,ix)=W(ix)
4     GP(ix,ix)=WP(ix)
      call mt(J,JT,N)
      call mm(T,G,J,N)
      call mm(TP,JT,T,N)
      call ms(F,TP,GP,N)
      call INV(N,F,FI,N,EE,IERR)
      if(IERR.ne.0)call report('Inversion error')
      call mz(GH,N)
      call mz(GHP,N)
      do 5 ix=1,N
      GH(ix,ix) =dsqrt(G(ix,ix))
5     GHP(ix,ix)=dsqrt(GP(ix,ix))
      call mm(T,JT,GH,N)
      call mm(TP,FI,T,N)
      call mm(T,J,TP,N)
      call mm(TP,GH,T,N)
      do 6 ix=1,N
      do 6 jx=1,N
6     MYA(ix,jx)=2.0d0*TP(ix,jx)
      do 7 ix=1,N
7     MYA(ix,ix)=MYA(ix,ix)-1.0d0
      if(iwra.eq.1)call wdm(6,N,MYA)
      call mm(T,JT,G,N)
      call mm(TP,FI,T,N)
      call mm(T,J,TP,N)
      do 8 ix=1,N
8     T(ix,ix)=T(ix,ix)-1.0d0
      call mv(TV,T,K,N)
      call mv(TU,GH,TV,N)
      do 9 ix=1,N
9     MYB(ix)=-2.0d0*TU(ix)
      if(iwrb.eq.1)call wdv(6,N,MYB)

      if(ltab)then
       open(50,file='FC.TAB')
       write(50,500)
500    format(/,'Franck-Condon Factors',/,60(1h-))
      endif

      proc=0.0d0
      nt=0
      call mm(T,G,GP,N)
      fac=dsqrt(2.0d0**N)*dsqrt(dsqrt(dabs(GETDET(T,N))))
     1*dsqrt(dabs(GETDET(J,N)/GETDET(F,N)))
      call mv(TV,G,K,N)
      sp1=sp(K,TV,N)
      call mv(TU,G,K,N)
      call mv(TV,JT,TU,N)
      call mv(TU,FI,TV,N)
      call mv(TV,J,TU,N)
      call mv(TU,G,TV,N)
      sp2=sp(K,TU,N)
      fac=fac*exp(-0.50d0*(sp1-sp2))
      proc=fac*fac*100.0d0
      write(6,300)fac,proc
300   format('<0|0*> = ',f10.3,' (',f6.2,'%)')
      nt=nt+1
      NSL=N
      allocate(si(LEXCL+1),sit(LEXCL+1))
      do 11 kk=1,LEXCL+1
11    si(kk)=0

c     transition dipole <0|u|*>=u0<0|*>+du/dQk <0|Qk|*>
      do 12 ix=1,3
12    ut(ix)=u0(ix)*fac
      if(LQ1)then
       do 10 kk=1,N
       si(1)=kk
c      <0|Q*k|0*>=sqrt(h/(2*wk)) <0|1*k>
       pp=FC(fac,si,1,np0,LEXCL,iwr,C,D,N)/dsqrt(2.0d0*GP(kk,kk))
       if(LDD)then
        do 13 ix=1,3
13      ut(ix)=ut(ix)+ddi(ix+3*(kk-1))*pp
       endif
10     write(6,1600)kk,pp
1600   format('<0|Q_',i3,'|0*> = ',g13.4)
      endif

      ds=ut(1)**2+ut(2)**2+ut(3)**2
      if(ltab)write(50,550)nt,E0,ds,proc
550   format(i8,f10.2,g13.4,' ',f10.2,'%',10i3)

c     states Nexc excited
      N7=1
      NQ1=LEXCL
      do 30000 Nexc=1,LEXCL
c     distribute Nexc excitations upon centers N7..NSL:
c     initial indices - all to the first center:
      do 41 iex=1,Nexc
41    si(iex)=N7
50000 do 11000 I=1,Nexc
c     check that the state is allowed:
      icen=1
      do 11000 II=I+1,Nexc
      if(si(II).eq.si(I))icen=icen+1
11000 if(icen.gt.NQ1)goto 12000
c
      EJ=0.0d0
      DO 1002 II=1,Nexc
1002  EJ=EJ+GP(si(II),si(II))

      if(iwr.gt.1)write(6,604)ii,(si(iex),iex=1,NExc)
604   format(' state ',i2,':',10i3)

c     <0|*>
      pp=FC(fac,si,Nexc,np0,LEXCL,iwr,C,D,N)
      proc=proc+pp*pp*100.0d0

      if(HW)then
       write(6,599)
599    format('<0|',$)
       do 107 i=1,Nexc
107    write(6,601)si(i)
601    format(i3,$)
       write(6,598)
598    format('>: ',$)
       write(6,602)pp,proc,EJ*CM
602    format(g12.4,' (',f10.3,'%), E= ',F10.2)
      endif

c     <0|u|*>=u0<0|*>+du/dQk<0|Qk|*>
c     <0|Qk|*>=sqrt(h/(2wk))(sqrt(vk)<0|*-1k>+sqrt(vk+1)<0|*+1k>)
      do 14 ix=1,3
14    ut(ix)=u0(ix)*pp
      if(LQ1)then
       if(LDD)then
        do 15 kk=1,N
        nuk=0
        do 16 ic=1,Nexc
16      if(si(ic).eq.kk)nuk=nuk+1
        pm=0.0d0
        if(nuk.gt.0)then
         ip=0
         ik=0
         do 17 ic=1,Nexc
         if(si(ic).eq.kk.and.ik.eq.0)then
          ik=ik+1
         else
          ip=ip+1
          sit(ip)=si(ic)
         endif
17       continue
         pm=dsqrt(dble(nuk))*FC(fac,sit,Nexc-1,np0,LEXCL,iwr,C,D,N)
        endif
        do 18 ic=1,Nexc
18      sit(ic)=si(ic)
        sit(Nexc+1)=kk
        pp=dsqrt(dble(nuk+1))*FC(fac,sit,Nexc+1,np0,LEXCL,iwr,C,D,N)
        qk=(pm+pp)/dsqrt(2.0d0*GP(kk,kk))
        do 15 ix=1,3
15      ut(ix)=ut(ix)+ddi(ix+3*(kk-1))*qk
       endif
      endif

      ds=ut(1)**2+ut(2)**2+ut(3)**2
      if(ltab.and.dabs(ds).gt.THR)then
       nt=nt+1
       write(50,550)nt,E0+EJ*CM,ds,proc,(si(i),i=1,Nexc)
      endif
c
c     find index to be changed
12000 do 80000 ic=Nexc,1,-1
80000 if(si(ic).lt.NSL)goto 90000
      goto 30000
90000 do 10000 i=ic+1,Nexc
10000 si(i)=si(ic)+1
      si(ic)=si(ic)+1
c
      goto 50000
30000 continue

      if(ltab)then
       write(50,551)
551    format(60(1h-))
      endif
      end
      

      subroutine mz(A,N)
      implicit none
      integer*4 n,i,j
      real*8 A(N,N)
      do 3 i=1,N
      do 3 j=1,N
3     A(i,j)=0.0d0
      return
      end

      SUBROUTINE INV(nr,a,ai,n,e,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      dimension ai(nr,nr),a(nr,nr),e(nr,2*nr)
C
      TOL=1.0d-10
10000 IERR=0
      DO 1  ii=1,n
      DO 1  jj=1,n
      e(ii,jj)=a(ii,jj)
1     e(ii,jj+n)=0.0D0
      do 13 ii=1,n
13    e(ii,ii+n)=1.0D0
c
      DO 2 ii=1,n-1
      if (ABS(e(ii,ii)).LE.TOL) then
       DO 3 io=ii+1,n
3      if (ABS(e(io,ii)).GT.TOL) goto 11
       IERR=1
       write(6,*)ii
       write(6,*)'tol = ',tol
       tol=tol*0.50d0
       if(tol.gt.1.0d-20)goto 10000
       RETURN
c
11     CONTINUE
       DO 4 kk=1,2*n
       w=e(ii,kk)
       e(ii,kk)=e(io,kk)
4      e(io,kk)=w
      ENDIF
      eii=e(ii,ii)
      DO 5 jj=ii+1,n
      e1=e(jj,ii)/eii
      DO 6 kk=ii+1, 2*n
6     e(jj,kk)=e(jj,kk)-e(ii,kk)*e1
5     e(jj,ii)=0.0D0
2     CONTINUE
c
      DO 7  i2=n,2,-1
      eii=e(i2,i2)
      DO 7 j2=i2-1,1,-1
      e1=e(j2,i2)/eii
      DO 9 kk=1, n
9     e(j2,kk+n)=e(j2,kk+n)-e(i2,kk+n)*e1
7     e(j2,i2)=0.0d0
c
      DO 10 ii=1,n
      ei=1.0d0/e(ii,ii)
      DO 12  jj=1,n
12    ai(ii,jj)=e(ii,jj+n)*ei
10    CONTINUE
c
      RETURN
      END

      subroutine mm(A,B,C,N)
c     matrix multiplication A = B x C
      implicit none
      integer*4 N,i,j,k
      real*8 A(N,N),B(N,N),C(N,N)
      do 1 i=1,N
      do 1 j=1,N
      A(i,j)=0.0d0
      do 1 k=1,N
1     A(i,j)=A(i,j)+B(i,k)*C(k,j)
      return
      end

      subroutine mv(A,B,C,N)
c     multiplication of matrix B by vector C, A = B . C
      implicit none
      integer*4 N,i,k
      real*8 A(*),B(N,N),C(*)
      do 1 i=1,N
      A(i)=0.0d0
      do 1 k=1,N
1     A(i)=A(i)+B(i,k)*C(k)
      return
      end

      function sp(A,C,N)
c     scalar product of two vectors
      implicit none
      integer*4 N,k
      real*8 A(*),C(*),sp
      sp=0.0d0
      do 1 k=1,N
1     sp=sp+A(k)*C(k)
      return
      end


      subroutine ms(A,B,C,N)
c     matrix sum A = B + C
      implicit none
      integer*4 N,i,j
      real*8 A(N,N),B(N,N),C(N,N)
      do 1 i=1,N
      do 1 j=1,N
1     A(i,j)=B(i,j)+C(i,j)
      return
      end

      subroutine mt(AT,A,N)
c     matrix transposition AT = A^T
      implicit none
      integer*4 N,i,j
      real*8 A(N,N),AT(N,N)
      do 1 i=1,N
      do 1 j=1,N
1     AT(i,j)=A(j,i)
      return
      end

      subroutine rw(o,N,a,b,i)
      implicit none
      integer*4 o,N,i,ns,ne,j
      real*8 a(*),b(*),cmtoau
      character*80 s
      cmtoau=219474.63d0
      if(i.ge.2)return
      i=i+1
      ns=1
1     read(o,9)s
9     format(a80)
      if(s(2:15).eq.'Frequencies --')then
       ne=min(N,ns+2)
       if(i.eq.1)read(s(16:len(s)),*)(a(j),j=ns,ne)
       if(i.eq.2)read(s(16:len(s)),*)(b(j),j=ns,ne)
c      convert to Hartree:
       do 2 j=ns,ne
       if(i.eq.1)a(j)=a(j)/cmtoau
2      if(i.eq.2)b(j)=b(j)/cmtoau
       ns=min(ns+3,N)
       if(ne.eq.N)return
      endif
      goto 1
      end

      subroutine rdx(f,nat,ic,ty,r,NG,lzmat)
      implicit none
      integer*4 nat,ic,ty(*),NG,ig98,l,i
      real*8 r(*)
      character*(*) f
      character*50 FN
      logical lzmat
      NG=0
      nat=0
      OPEN(2,FILE=f)
1     READ(2,2000,END=1000)FN
2000  FORMAT(A50)
      IF((lzmat.and.(FN(19:39).EQ.'Z-Matrix orientation:'.OR.
     1               FN(26:46).EQ.'Z-Matrix orientation:'.OR.
     1               FN(20:37).EQ.'Input orientation:'.OR.
     1               FN(27:44).EQ.'Input orientation:'))
     1    .OR.
     1    ((.not.lzmat).and.
     2              (FN(20:40).EQ.'Standard orientation:'.OR.
     2               FN(26:46).EQ.'Standard orientation:')))THEN
       NG=NG+1
       ig98=0
       if(FN(26:46).EQ.'Z-Matrix orientation:'.OR.
     1    FN(27:44).EQ.'Input orientation:'.OR.
     1    FN(26:46).EQ.'Standard orientation:')ig98=1
       DO 4 i=1,4
4      READ(2,*)
       l=0
5      READ(2,2000)FN
       IF(FN(2:4).NE.'---')THEN
        l=l+1
        if(ic.eq.1)then
         BACKSPACE 2
         if(ig98.eq.0)then
          READ(2,*)ty(l),ty(l),(r(i+(l-1)*3),i=1,3)
         else
          READ(2,*)ty(l),ty(l),r(1+(l-1)*3),(r(i+(l-1)*3),i=1,3)
         endif
        endif
        GOTO 5
       ENDIF
      ENDIF
      GOTO 1
1000  CLOSE(2)
      nat=l
      return
      end

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

      SUBROUTINE rdm(io,N,M,ic)
      IMPLICIT none
      integer*4 N,N1,N3,io,LN,J,ic
      real*8 M(N,N)
      ic=ic+1
      read(io,*)
      N1=1
1     N3=min(N1+4,N)
      read(io,*)
      DO 130 LN=1,N
130   READ(io,*)M(LN,N1),(M(LN,J),J=N1,N3)
      N1=N1+5
      IF(N3.LT.N)GOTO 1
      return
      end

      SUBROUTINE wdv(io,N,V)
      IMPLICIT none
      integer*4 N,io,LN
      real*8 V(*)
      write(io,*)
      DO 130 LN=1,N
130   write(io,500)LN,V(LN)
500   format(i7,D15.6)
      return
      end

      SUBROUTINE wdm(io,N,M)
      IMPLICIT none
      integer*4 N,N1,N3,io,LN,J
      real*8 M(N,N)
      write(io,*)
      N1=1
1     N3=min(N1+4,N)
      write(io,400)(J,J=N1,N3)
400   format(3x,5i15)
      DO 130 LN=1,N
130   write(io,500)LN,(M(LN,J),J=N1,N3)
500   format(i7,5D15.6)
      N1=N1+5
      IF(N3.LT.N)GOTO 1
      return
      end

      SUBROUTINE rdv(io,N,V,ic)
      IMPLICIT none
      integer*4 N,io,LN,ic
      real*8 V(*)
      ic=ic+1
      read(io,*)
      read(io,*)
      DO 130 LN=1,N
130   READ(9,*)V(LN),V(LN)
      return
      end

      FUNCTION GETDET(A,N)
c     determinant of a matrix
      IMPLICIT none
      integer*4 N,I,J,K
      logical DETEXISTS
      real*8 A(N,N),M,TEMP,L,GETDET
      real*8, allocatable::ELEM(:,:)
      allocate(ELEM(N,N))
      DO 1 I=1,N 
      DO 1 J=1,N 
1     ELEM(I,J)=A(I,J)
      DETEXISTS=.TRUE.
      L=1.0d0
      DO 2 K=1,N-1
      IF(DABS(ELEM(K,K)).LE.1.0D-20)THEN
       DETEXISTS=.FALSE.
       DO 3 I=K+1,N
       IF(ELEM(I,K).NE.0.0d0)THEN
        DO 4 J=1,N
        TEMP=ELEM(I,J)
        ELEM(I,J)=ELEM(K,J)
4       ELEM(K,J)=TEMP
        DETEXISTS=.TRUE.
        L=-L
        EXIT
       ENDIF
3      CONTINUE
       IF (DETEXISTS.EQV..FALSE.)THEN 
        GETDET = 0.0d0
        RETURN
       ENDIF
      ENDIF
      DO 2 J=K+1,N
      M=ELEM(J,K)/ELEM(K,K)
      DO 2 I=K+1,N
2     ELEM(J,I)=ELEM(J,I)-M*ELEM(K,I)
      GETDET =L
      DO 5 I=1,N
5     GETDET=GETDET*ELEM(I,I)
      RETURN
      END

      function jeje(np,nx,it,NN,ee,np0,LEXCL)
      implicit none
      integer*4 np,nx,it(*),NN(*),np0,LEXCL,ee(np0,LEXCL),je,jeje,jj,ie
      je=0
      do 104 jj=1,np
      if(nx.ne.NN(jj))goto 104
      do 1041 ie=1,NN(jj)
1041  if(it(ie).ne.ee(jj,ie))goto 104
      je=jj
      goto 1042
104   continue
1042  jeje=je
      return
      end

      subroutine digest(np,ie,np0,LEXCL,p,NN,it,iexc,T)
      implicit none
      integer*4 je,np,np0,LEXCL,NN(*),it(*),ie(np0,LEXCL),jj,iexc,
     1jeje
      real*8 T,p(*)

c     does it exist already within 1 ... np?:
      je=jeje(np,iexc,it,NN,ie,np0,LEXCL)

      if(je.ne.0)then
c      the term already exist in the expansion as je^th - just add it
       p(je)=p(je)+T
      else
c      term not found - add as new term
       np=np+1
       if(np.gt.np0)call report('Too many terms')
       p(np)=T
       NN(np)=iexc
       do 1 jj=1,iexc
1      ie(np,jj)=it(jj)
      endif
      return
      end

      function FC(fac,si,Nexc,np0,LEXCL,iwr,C,D,N)
c     FC = Franc Condon factor for <0|si>
c     fac = <0|0*>
c     np0 ... dimension of working buffer
      implicit none
      integer*4 si(*),Nexc,np,iex,ii,nu,jj,ip,N,ic,jold,
     1nuj,jc,kk,iwr,LEXCL,np0,iq,jq
      real*8 FC,fac,D(*),C(N,N),pini
      real*8,allocatable::p(:)
      integer*4,allocatable::NN(:),ie(:,:),it(:)
      np=1
      allocate(p(np0),NN(np0),ie(np0,LEXCL),it(np0))
      p(np)=1.0d0
      NN(np)=NExc
c     if(NExc.eq.2.and.si(1).eq.1.and.si(2).eq.12)then
c      write(6,*)'D1',D(1),' D12',D(12),'C1 12',C(1,12)
c      iwr=2
c     else
c      iwr=1
c     endif
      if(iwr.gt.1)write(6,604)(si(iex),iex=1,NExc)
604   format(/,'<0|e>, e = ',10i3,' requested:')
      do 102 iex=1,NN(np)
102   ie(np,iex)=si(iex)
c     expand reccurent formula into a sum and shrink back to <0|0>
777   do 101 ii=1,np
c     term ii - reduce excitation one by one:
      pini=p(ii)
      if(iwr.gt.1)write(6,607)ii,p(ii),(ie(ii,iex),iex=1,NN(ii))
607   format(' term ',i2,':',f10.6,10i3)
      do 101 iex=1,NN(ii)
      iq=ie(ii,iex)
      if(iq.gt.0)then
c      <0| si_nu> reduce to <0| i-1>, <0|i-2>, <0|i-1,j-1>(j<>i)
       nu=0
       do 1032 jj=1,NN(ii)
1032   if(ie(ii,jj).eq.iq)nu=nu+1

c      write term <0|v'-1_iex> into  string it:
       ip=0
       ic=0
       do 1031 jj=1,NN(ii)
       if(iq.eq.ie(ii,jj).and.ic.lt.1)then
        ic=ic+1
       else
        ip=ip+1
        it(ip)=ie(ii,jj)
       endif
1031   continue
c       
       call digest(np,ie,np0,LEXCL,p,NN,it,NN(ii)-1,
     1 pini*D(iq)/dsqrt(dble(2*nu)))
       if(iwr.gt.1)write(6,606)(it(jj),jj=1,NN(ii)-1)
606    format('digest ',10i3)

       if(nu.gt.1)then
c       write term <0|v'-2_iex> into string it:
        ic=0
        ip=0
        do 1033 jj=1,NN(ii)
        if(ie(ii,jj).eq.iq.and.ic.lt.2)then
         ic=ic+1
        else
         ip=ip+1
         it(ip)=ie(ii,jj)
        endif
1033    continue
        call digest(np,ie,np0,LEXCL,p,NN,it,NN(ii)-2,
     1  pini*dsqrt(dble(nu-1)/dble(nu))*C(iq,iq))
        if(iwr.gt.1)write(6,606)(it(jj),jj=1,NN(ii)-2)
       endif

       jold=0
       do 106 jj=1,NN(ii)
       jq=ie(ii,jj)
       if(jq.ne.iq.and.jq.ne.jold)then
        nuj=0
        do 1034 kk=1,NN(ii)
1034    if(jq.eq.ie(ii,kk))nuj=nuj+1
c       write term <0|v'-1_iex-1_j, j<>iex> into it:
        ic=0
        jc=0
        ip=0
        do 1035 kk=1,NN(ii)
        if(ie(ii,kk).eq.iq.and.ic.lt.1)then
         ic=ic+1
        else
         if(ie(ii,kk).eq.jq.and.jc.lt.1)then
          jc=jc+1
         else
          ip=ip+1
          it(ip)=ie(ii,kk)
         endif
        endif
1035    continue
        call digest(np,ie,np0,LEXCL,p,NN,it,NN(ii)-2,
     1  pini*dsqrt(dble(nuj)/dble(nu))*C(iq,jq))
        if(iwr.gt.1)write(6,606)(it(kk),kk=1,NN(ii)-2)
       endif
106    jold=jq

       if(iwr.gt.1)then
        write(6,*)'np now ',np
        do 108 jj=1,np
108     write(6,605)jj,NN(jj),p(jj),(ie(jj,kk),kk=1,NN(jj))
       endif
       
c      eliminate the old term and start over
       do 105 jj=ii,np-1
       p(jj)=p(jj+1)
       NN(jj)=NN(jj+1)
       do 105 kk=1,NN(jj)
105    ie(jj,kk)=ie(jj+1,kk)
       np=np-1
       if(iwr.gt.1)then
        write(6,*)'np after elimination ',np
        do 107 jj=1,np
107     write(6,605)jj,NN(jj),p(jj),(ie(jj,kk),kk=1,NN(jj))
605     format(i3,i2,g9.3,10i3)
       endif

       goto 777
       
      endif
101   continue

      if(np.ne.1)call report('np <> 1')

      FC=p(1)*fac
c     if(Nexc.eq.2)then
c      write(6,*)p(1),fac
c      write(6,*)'input something'
c      read(5,*)jj
c     endif
      return
      end

      subroutine rdd(f,d,nat,Nroot,u0,iwr)
c     read transition dipole derivatives from excited state freq calc
      implicit none
      character*(*) f
      real*8 d(*),u(3),u0(*),step,sign,stepau
      integer*4 nat,ln,ia,xa,i,Nroot,xar,iar,ibr,ib,ic,ii,ix,iwr
      character*80 s80
      open(9,file=f,status='old')
      ln=0
c     where it starts:
1     read(9,900,end=88,err=88)s80
900   format(a80)
      ln=ln+1
      if(s80(1:2).eq.' #')then
       ic=0
       do 3 i=1,len(s80)-3
       if(s80(i:i+4).eq.' freq'.or.s80(i:i+4).eq.' freq')ic=ic+1
3      if(s80(i:i+2).eq.' td'  .or.s80(i:i+2).eq.' TD'  )ic=ic+1
       if(ic.eq.2)then
        write(6,*)'TD starts at line ',ln
        goto 2
       endif
      endif
      goto 1
88    close(88)
      call report('TD freq not found')

2     ia=1
      xa=1
      ib=0
      ii=0
      step=0.001d0
      do 7 ix=1,9*nat
7     d(ix)=0.0d0
4     read(9,900,end=78,err=78)s80
      if(s80(2:14).eq.'Nuclear step=')then
       write(6,*)s80(1:23)
       read(s80(15:23),*)step
      endif
      if(s80(2:24).eq.'Re-enter D2Numr: IAtom=')then
       read(s80(25:27),*)iar
       read(s80(34:34),*)xar
       read(s80(42:43),*)ibr
       if(ia.ne.iar)call report('ia <> iar')
       if(xa.ne.xar)call report('xa <> xar')
       if(ib.ne.ibr)call report('ib <> ibr')
      endif
      if(s80(2:65).eq.'Ground to excited state transition electric dipol
     1e moments (Au):')then
       if(ii.eq.0)then
        write(6,*)'Zero point'
        do 5 i=1,Nroot
5       read(9,*)
        read(9,*)u0(1),(u0(ix),ix=1,3)
       else
        ib=ib+1
        if(ib.gt.2)then
         ib=1
         xa=xa+1
         if(xa.gt.3)then
          xa=1
          ia=ia+1
         endif
        endif
        if(iwr.gt.1)write(6,*)' d ',ia,xa,ib,iar,xar,ibr
        do 6 i=1,Nroot
6       read(9,*)
        read(9,*)u(1),(u(ix),ix=1,3)
        if(ib.eq.1)then
         sign=1.0d0
        else
         sign=-1.0d0
        endif
        do 8 ix=1,3
8       d(ix+3*(xa-1)+9*(ia-1))=d(ix+3*(xa-1)+9*(ia-1))+sign*u(ix)
       endif
       ii=ii+1
       if(ia.eq.nat.and.ix.eq.3.and.ib.eq.2)goto 78
      endif
      goto 4
78    close(9)
      stepau=step/0.529177d0
      do 9 ix=1,9*nat
9     d(ix)=d(ix)/stepau
      open(40,file='DE.TEN')
      write(40,*)'atom coord dipx        dipy        dipz'
      do 10 ia=1,nat
      do 10 xa=1,3
 10   write(40,400)ia,xa,(d(ix+3*(xa-1)+9*(nat-1)),ix=1,3)
400   format(i5,i2,3g12.4)
      close(40)

      return
      end
c     ==============================================================
      SUBROUTINE READS(N3,S,E,NQ)
      IMPLICIT none
      integer*4 N3,NQ,N,nat3,nat,i,J,ix
      real*8 S(N3,N3),E(*),CM,SFAC
      CM=219470.0d0
      SFAC=0.0234280d0
C     1 a.u.= 2.1947E5 cm^-1
      open(4,file='F.INP')
      read(4,*)NQ,nat3,nat
      N=NAT3
      do 1 i=1,NAT
1     read(4,*)
      read(4,*)
      DO 2 I=1,NAT
      DO 2 J=1,NQ
2     read(4,*)(s(3*(i-1)+ix,J),ix=1,2),
     1(s(3*(i-1)+ix,J),ix=1,3)
      read(4,*)
      READ(4,4000)(E(I),I=1,NQ)
4000  FORMAT(6F11.6)
      close(4)
      DO 3 I=1,NQ
3     E(I)=E(I)/CM
      DO 5 I=1,N
      DO 5 J=1,N
5     S(I,J)=S(I,J)*SFAC
      RETURN
      end
c     ==============================================================
      subroutine trafo(dd,ddi,nat,N)
      implicit none
      integer*4 nat,NQ,N,ix,iq,ia,xa
      real*8 dd(*),ddi(*)
      real*8,allocatable::e(:),s(:,:)
      allocate(e(3*nat),s(3*nat,3*nat))
      call READS(3*nat,s,e,NQ)
      if(N.ne.NQ)call report('N <> NQ')
      do 1 ix=1,3
      do 1 iq=1,N
      ddi(ix+3*(iq-1))=0.0d0
      do 1 ia=1,nat
      do 1 xa=1,3
c     s-opposite order of normal modes:
1     ddi(ix+3*(iq-1))=ddi(ix+3*(iq-1))+dd(ix+3*(xa-1)+9*(ia-1))
     1*s(xa+3*(ia-1),N-iq+1)
      write(6,*)'dipoles transformed to normal modes'
      return
      end
c     ==============================================================
      subroutine zm(M,N)
      implicit none
      integer*4 N,i,j
      real*8 M(N,N)
      do 1 i=1,N
      do 1 j=1,N
1     M(i,j)=0.0d0
      return
      end
c     ==============================================================
      subroutine zv(J,N)
      implicit none
      real*8 J(*)
      integer*4 N,i
      do 1 i=1,N
1     J(i)=0.0d0
      return
      end