program dusch
c     Duschinsky transformation
      implicit none
      real*8 bohr,u(3,3),fac,TOL,facv,wlim,wfix,
     1overlap,wmin,TOLN,del,e00,wsmin,wsmax,kfac
      integer*4 ia,N,nat,NE,NQ,IERR,maxit,NB,NP,nps,ncut,LEXCL,isu,
     1ifx,NF
      real*8,allocatable::rg(:),re(:),se(:,:),sg(:,:),JT(:,:),sgf(:,:),
     1eg(:),ep(:),K(:),m(:),F(:,:),G(:,:),GP(:,:),T(:,:),TP(:,:),
     1FI(:,:),E2(:,:),GH(:,:),GHP(:,:),A(:,:),TV(:),TU(:),
     1B(:),J(:,:),C(:,:),D(:),E(:,:),ref(:),rgf(:),K0(:)
      integer*4,allocatable::z(:)
      real amas(89)
      data amas/1.007825,4.002603,
     2  6.941, 9.012,   10.810,12.000,14.003074,15.9949,18.998,20.179,
     3 22.990,24.305,   26.981,28.086,30.974,32.060,35.453,39.948,
     4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847,
     4 58.933,58.700,63.546,65.380,
     4                  69.720,72.590,74.922,78.960,79.904,83.800,
     5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070,
     5 102.906,106.400,107.868,112.410,
     5                 114.82,118.69,121.75,127.600,126.905,131.300,
     6 132.905,137.330,138.906,
     6                 140.120,140.908,144.240,145.000,150.400,
     6 151.960,157.250,158.925,162.500,164.930,167.260,168.934,
     6 173.040,174.970,
     6 178.490,180.948,183.850,186.207,190.200,192.220,195.090,
     6 196.967,207.590,204.370,207.200,208.980,210.000,210.001,
     6 222.02,
     7 223.000,226.025,227.028/     
      logical lit,lwr,lfix,lmscan,lspec,lshifte,lshifti,lshiftq,
     1llin,lcor,lonly,lsfix,lcre,lcom,ldz

      bohr=0.529177d0
      TOL=1.0d-10

      call readopt(lit,maxit,lwr,wmin,lfix,lmscan,NB,NP,TOLN,
     1lspec,del,e00,wsmin,wsmax,nps,lshifte,lshifti,lshiftq,
     1ncut,LEXCL,isu,llin,ifx,wlim,lcor,lonly,lsfix,lcre,kfac,
     1lcom,wfix,ldz)
c
c     Read Geometry of the ground and excited states:
      allocate(z(1),rg(1))
      if(lfix)then
       call rdxx('ground/fixed.x',0,nat,z,rg)
      else
       call rdxx('ground/FILE.X',0,nat,z,rg)
      endif
      deallocate(z,rg)
      N=3*nat
      allocate(rg(N),z(nat),re(N),sg(N,N),se(N,N),eg(N),ep(N),m(nat),
     1ref(N),rgf(N),sgf(N,N))
      if(lfix)then
       call rdxx('ground/fixed.x',1,nat,z,rg)
       call rdxx('excited/fixed.x',1,nat,z,re)
      else
       call rdxx('ground/FILE.X',1,nat,z,rg)
       call rdxx('excited/FILE.X',1,nat,z,re)
      endif
c     assign masses, eventually do isotopic substitution:
      do 6 ia=1,nat
6     m(ia)=dble(amas(z(ia)))
      call readm(m)

c     read the ground (sg) and excited (se) S-matrices:
      write(6,*)'reading ground state S-matrix'
      call readsi(N,sg,eg,NQ,'ground/F.INP',rgf,ldz)
      if(llin)then
       write(6,*)'excited state S-matrix same as ground'
       call readsi(N,se,ep,NE,'ground/F.INP',ref,ldz)
      else
       write(6,*)'reading excited state S-matrix'
       call readsi(N,se,ep,NE,'excited/F.INP',ref,ldz)
      endif

      call ortog(N,NE,se,m)
      call ortog(N,NQ,sg,m)
      write(6,*)'S-matrices corrected'

      call fixmodes(N,NQ,NE,eg,sg,ep,se,amas,z,ifx,wlim,lcor,lonly,
     1sgf,NF,wfix)
c     fix se for deleted modes:
      if(lsfix.and.NF.gt.0)call seproj(N,NF,NQ,se,sgf,m)

      if(lfix)then
c      Find the transformation matrix from rgf to rg
       call xst(nat,rg,rgf,m,z,u)
c      transform ground S-matrix to the rg system:
       call trse(N,sg,NQ,u)
c      Find the transformation matrix from ref to re
       call xst(nat,re,ref,m,z,u)
c      transform excited S-matrix to the re system:
       call trse(N,se,NQ,u)
      else
c      Find the transformation matrix, transform and shift re to rg:
       call xst(nat,rg,re,m,z,u)
c      transform se to the ground system:
       call trse(N,se,NQ,u)
      endif

c     experimental option-correct excited state so that it can be 
c     reached with available modes
      if(lcre)then
       call cre(N,NQ,m,rg,re,se,'e')
       call wrx(nat,re,z,'test2.x','test')
       call cre(N,NQ,m,rg,re,sg,'g')
       call wrx(nat,re,z,'test3.x','test')
      endif

c     if required, replace excited S-vectors by ground:
      if(lcom)call rpl(N,NQ,nat,m,sg,se)
c     Make Duschinsk's and derived matrices and vectors:
      allocate(JT(NQ,NQ),K(NQ),G(NQ,NQ),GP(NQ,NQ),GH(NQ,NQ),GHP(NQ,NQ),
     1J(NQ,NQ),T(NQ,NQ),TP(NQ,NQ),F(NQ,NQ),FI(NQ,NQ),E2(NQ,2*NQ),
     1A(NQ,NQ),B(NQ),TV(NQ),TU(NQ),C(NQ,NQ),D(NQ),E(NQ,NQ),K0(NQ))
      call zzz(NQ,JT,K,G,GP,GH,GHP,J,T,TP,F,FI,E2,A,B,TV,TU,C,D,E,K0)

      IERR=0
      fac=overlap(1,NQ,N,nat,G,GP,GH,GHP,eg,ep,JT,J,sg,se,
     1m,k,rg,re,bohr,T,TP,F,FI,E2,TOL,IERR,A,TV,TU,B,C,D,E,kfac)
      write(6,300)fac,fac*fac*100.0d0
300   format('<0|0*> = ',g12.4,' (',f6.2,'%)')
      
      write(6,*)' Vertical approximation:'
      call dofac(NQ,facv,J,JT,F,FI,G,GP,K0)
      write(6,300)facv,facv*facv*100.0d0

      if(lit)call iterate(maxit,NQ,lwr,wmin,
     1N,nat,G,GP,GH,GHP,eg,ep,JT,J,sg,se,
     1m,k,rg,re,bohr,T,TP,F,FI,E2,TOL,IERR,A,TV,TU,B,C,D,E,kfac)

c     calculate normal mode gradient:
      call gradq(nat,N,m,ldz)
c     normal mode gradient in excited state but ground state geometry:
      call grade(nat,N,ldz)
c     calculate normal mode shift corresponding to geometry change:
      if(lshiftq)call shiftq(nat,N,NQ,se,ep,sg,eg)
c     different route:
      if(lshifti)call shifti(nat,N,NQ,se,ep)
c     different route, excited n's:
      if(lshifte)call shifte(nat,N,NQ,sg,se,eg,ep,m,lmscan,fac,
     1A,B,C,D,E,F,FI,G,GP,GH,GHP,J,JT,K,NB,NP,TOLN,
     2lspec,del,e00,wsmin,wsmax,nps,ncut,LEXCL,
     1isu)
     
      end
c     ==============================================================
      subroutine ortog(N,NQ,s,m)
c     orthogonalize s-vectors, to correct numerical inaccuracies:
      implicit none
      integer*4 N,NQ,i,j,ii,a
      real*8 s(N,N),m(N/3),sij,sn,dm
      real*8,allocatable::st(:,:)
      allocate(st(N,N))
      do 1 a=1,N/3
      dm=dsqrt(m(a))
      do 1 i=1,3
      ii=i+3*(a-1)
      do 1 j=1,NQ
 1    st(ii,j)=dm*s(ii,j)

      do 15 i=1,NQ
      do 16 j=1,i-1
      sij=0.0d0
      do 17 ii=1,N
17    sij=sij+st(ii,i)*st(ii,j)
      do 16 ii=1,N
16    st(ii,i)=st(ii,i)-sij*st(ii,j)
      sn=0.0d0
      do 18 ii=1,N
18    sn=sn+st(ii,i)**2
      sn=1.0d0/dsqrt(sn)
      do 15 ii=1,N
15    st(ii,i)=st(ii,i)*sn

      do 2 a=1,N/3
      dm=dsqrt(m(a))
      do 2 i=1,3
      ii=i+3*(a-1)
      do 2 j=1,NQ
 2    s(ii,j)=st(ii,j)/dm
      return
      end
c     ==============================================================
      subroutine cre(N,NQ,m,rg,re,s,y)
      implicit none
      integer*4 N,NQ,i,j,ia,ii
      real*8 rg(N),re(N),s(N,N),xn,sp,t(3,3),m(N/3),dm
      real*8, allocatable::x(:),xv(:),v(:),st(:,:)
      character*1 y
      allocate(x(N),xv(NQ),v(N),st(N,N))
      if(y.eq.'e')then
       open(44,file='TMX.TXT')
       do 9 i=1,3
9      read(44,*)(t(j,i),j=1,3)
       close(44)
c      excited state S-matrix in ground-like orientation:
       do 14 ia=1,N/3
       dm=dsqrt(m(ia))
       do 14 j=1,NQ
       do 14 i=1,3
       ii=3*(ia-1)
14     st(i+ii,j)=(s(1+ii,j)*t(1,i)+s(2+ii,j)*t(2,i)+s(3+ii,j)*t(3,i))
     1 *dm
      else
       do 141 ia=1,N/3
       dm=dsqrt(m(ia))
       do 141 j=1,NQ
       do 141 i=1,3
141    st(i+3*(ia-1),j)=s(i+3*(ia-1),j)*dm
      endif

      x=re-rg
      xn=sp(x,x,N)
      write(6,600)xn
600   format(' Correcting excited geometry, old shift',f12.3)
      do 1 i=1,NQ
      do 2 j=1,N
2     v(j)=st(j,i)
1     xv(i)=sp(x,v,N)
      x=0.0d0
      do 3 i=1,NQ
      do 4 j=1,N
4     v(j)=st(j,i)
      do 3 j=1,N
3     x(j)=x(j)+xv(i)*v(j)
      xn=sp(x,x,N)
      write(6,601)xn
601   format('                              new shift',f12.3)
      re=rg+x
      return
      end
c     ==============================================================
      subroutine wdxx(s,n,z,r)
      implicit none
      character*(*) s
      integer*4 n,z(*),ia,ix
      real*8 r(*)
      open(9,file=s)
      write(9,*)'geometry'
      write(9,*)n
      do 1 ia=1,n
1     write(9,90)z(ia),(r(3*(ia-1)+ix),ix=1,3)
90    format(i3,3f12.6)
      close(9)
      return
      end
c     ==============================================================
      subroutine rdxx(s,ic,n,z,r)
      implicit none
      character*(*) s
      integer*4 ic,n,z(*),ia,ix
      real*8 r(*)
      open(9,file=s)
      read(9,*)
      read(9,*)n
      if(ic.ne.0)then
       do 1 ia=1,n
1      read(9,*)z(ia),(r(3*(ia-1)+ix),ix=1,3)
      endif
      close(9)
      return
      end
c     ==============================================================
      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
c     ==============================================================
      subroutine norm(v,n)
      implicit none
      integer*4 n
      real*8 v(n),sp,vn
      vn=dsqrt(sp(v,v,n))
      v=v/vn
      return
      end
c     ==============================================================
      FUNCTION LOGDET(A,N,NQ)
c     logarithm of a determinant of a matrix
      IMPLICIT none
      integer*4 N,I,J,K,NQ
      logical DETEXISTS
      real*8 A(N,N),M,TEMP,LOGDET
      real*8, allocatable::ELEM(:,:)
      allocate(ELEM(NQ,NQ))
      DO 1 I=1,NQ 
      DO 1 J=1,NQ 
1     ELEM(I,J)=A(I,J)
      DETEXISTS=.TRUE.
c     L=1.0d0
      DO 2 K=1,NQ-1
      IF(ABS(ELEM(K,K)).LE.1.0d-20)THEN
       DETEXISTS=.FALSE.
       DO 3 I=K+1,NQ
       IF(ELEM(I,K).NE.0.0d0)THEN
        DO 4 J=1,NQ
        TEMP=ELEM(I,J)
        ELEM(I,J)=ELEM(K,J)
4       ELEM(K,J)=TEMP
        DETEXISTS=.TRUE.
c       L=-L
        EXIT
       ENDIF
3      CONTINUE
       IF (DETEXISTS.EQV..FALSE.)THEN 
        LOGDET = 0.0d0
        RETURN
       ENDIF
      ENDIF
      DO 2 J=K+1,NQ
      M=ELEM(J,K)/ELEM(K,K)
      DO 2 I=K+1,NQ
2     ELEM(J,I)=ELEM(J,I)-M*ELEM(K,I)
      LOGDET =0.0d0
      DO 5 I=1,NQ
      if(ELEM(I,I).lt.0.0d0)then
c      L=-L
       LOGDET=LOGDET+log(-ELEM(I,I))
      else
       LOGDET=LOGDET+log( ELEM(I,I))
      endif
5     continue
      RETURN
      END
c     ==============================================================
      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
c     ==============================================================
      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
c     ==============================================================
      SUBROUTINE wdv(io,N,V,s)
      IMPLICIT none
      integer*4 N,io,LN
      real*8 V(*)
      character*(*) s
      call wlabel(io,s)
      write(io,501)1
501   format(i14)
      DO 130 LN=1,N
130   write(io,500)LN,V(LN)
500   format(i7,1x,D14.6)
      return
      end
c     ==============================================================
      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
c     ==============================================================
      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
c     ==============================================================
      subroutine vz(A,N)
      implicit none
      integer*4 n,i
      real*8 A(*)
      do 3 i=1,N
3     A(i)=0.0d0
      return
      end
c     ==============================================================
      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
c     ==============================================================
      subroutine trse(N,s,NQ,t)
      implicit none
      integer*4 N,NQ,i,a,ia,ix
      real*8 s(N,N),t(3,3),v(3)
      do 1 i=1,NQ
      do 1 ia=1,N/3
      a=3*(ia-1)
      do 2 ix=1,3
2     v(ix)=s(ix+a,i)
      do 1 ix=1,3
1     s(ix+a,i)=t(1,ix)*v(1)+t(2,ix)*v(2)+t(3,ix)*v(3)
      return
      end
c     ==============================================================
      subroutine wlabel(io,s)
      implicit none
      character*(*) s
      integer*4 J,io
      write(io,*)
      write(io,101)
      do 1 J=1,len(s)
1     write(io,103)s(J:J)
      write(io,*)
      write(io,101)
101   format(' ',$)
      do 2 J=1,len(s)
2     write(io,103)'-'
103   format(a1,$)
      write(io,*)
      return
      end
c     ==============================================================
      SUBROUTINE wwm(io,N,e,s)
      IMPLICIT none
      integer*4 io,N,i
      real*8 e(*),CM
      character*(*) s
      CM=219474.630d0
      call wlabel(io,s)
      write(io,100)(e(i)*CM,i=1,N)
100   format(6f11.2)
      return
      end
c     ==============================================================
      SUBROUTINE wds(io,N0,N,M,s)
      IMPLICIT none
      integer*4 N,N1,N3,io,LN,J,N0
      real*8 M(N0,N0)
      character*(*) s
      call wlabel(io,s)
      N1=1
1     N3=min(N1+4,N)
      write(io,100)(J,J=N1,N3)
100   format(10x,5i14)
      DO 130 LN=1,N0
130   write(io,200)LN,(M(LN,J),J=N1,N3)
200   format(i7,1x,5e14.6)
      N1=N1+5
      IF(N3.LT.N)GOTO 1
      return
      end
c     ==============================================================
      SUBROUTINE wdm(io,N,M,s)
      IMPLICIT none
      integer*4 N,N1,N3,io,LN,J
      real*8 M(N,N)
      character*(*) s
      call wlabel(io,s)
      N1=1
1     N3=min(N1+4,N)
      write(io,100)(J,J=N1,N3)
100   format(10x,5i14)
      DO 130 LN=1,N
130   write(io,200)LN,(M(LN,J),J=N1,N3)
200   format(i7,1x,5e14.6)
      N1=N1+5
      IF(N3.LT.N)GOTO 1
      return
      end
c     ==============================================================
      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end
c     ==============================================================
      subroutine readm(m)
      IMPLICIT none
      real*8 m(*)
      logical lex
      character*3 key
      integer*4 i
      inquire(file='DUSCH.OPT',exist=lex)
      if(lex)then
       open(11,file='DUSCH.OPT')
1      read(11,900,end=11,err=11)key
 900   format(a3)
       if(key.eq.'MAS')read(11,*)i,m(i)
       goto 1
11     close(11)
      endif
      return
      end
c     ==============================================================
      subroutine readopt(lit,maxit,lwr,wmin,lfix,lmscan,NB,NP,TOL,
     1lspec,del,e00,wsmin,wsmax,nps,lshifte,lshifti,lshiftq,ncut,
     1LEXCL,isu,llin,ifx,wlim,lcor,lonly,lsfix,lcre,kfac,lcom,wfix,
     1ldz)
      IMPLICIT none
      integer*4 maxit,NB,NP,nps,ncut,LEXCL,isu,ifx,nblo
      real*8 wmin,TOL,e00,wsmin,wsmax,del,wlim,kfac,wfix
      logical lit,lwr,lex,lfix,lmscan,lspec,lshifte,lshifti,lshiftq,
     1llin,lcor,lonly,lsfix,lcre,lcom,ldz
      character*13 key
c     delete zero modes:
      ldz=.true.
c     replace negative modes by this value if required (in cm-1):
      wfix=100.0d0
c     ensure that Q Q' and K are in the same space
      lcom=.false.
c     factor for K shift
      kfac=1.0d0

c     fix final geometry according to available modes:
      lcre=.false.
c
c     fix se matrix for deleted modes:
      lsfix=.false.
c     
c     ilinear shift = excited modes same as ground
      llin=.false.
c
c     type of FC sum
      isu=2
c     maximal excitation for statesum1:
      LEXCL=3
c     if ncut>0, modify excited q definition according to the shift
      ncut=0
c     various paths to get the normal modeshift:
      lshifte=.true.
      lshifti=.true.
      lshiftq=.true.
c     make a spectrum:
      lspec=.false.
c     spectral parameters:e00,wmin,wmax (all cm-1) and number of points:
      e00=0.0d0
      wsmin=-1000.0d0
      wsmax= 1000.0d0
      nps=4001
c     spectral bandwidth:
      del=10.0d0
c     reorder excited modes according to ground:
      lcor=.false.
c     to CORREL, only show best match, no action:
      lonly=.false.
c     use fixed geometries, do not alling:
      lfix=.false.
c     scan modes -NP...NP:
      NP=10
c     quantum number biffer size:
      NB=10
c     take quantum number cut off:
      TOL=0.1d0
      lit=.false.
      wmin=1.0d0
      maxit=100
      lwr=.true.
      lmscan=.false.
c     what to do with negative modes:
c     ifx = 1 make them +100 cm-1
c           2 delete
c           3 borrow ground or excited state frequency
c           4 delete if w < wlim
      ifx=1
c     limit for mode deletion cm-1:
      wlim=600.0d0
      inquire(file='DUSCH.OPT',exist=lex)
      if(lex)then
       open(11,file='DUSCH.OPT')
1      read(11,900,end=11,err=11)key
 900   format(a13)
       if(key(1:3).eq.'BLO')then
        write(6,*)' Bloccked modes will be read again in fixmodes'
        read(11,*)nblo
        read(11,*)
        if(nblo.lt.0)read(11,*)
       endif
       if(key(1:3).eq.'LDZ')read(11,*)ldz
       if(key(1:3).eq.'COR')read(11,*)lcor
       if(key(1:3).eq.'DEL')read(11,*)del
       if(key(1:3).eq.'E00')read(11,*)e00
       if(key(1:3).eq.'FIX')read(11,*)ifx
       if(key(1:3).eq.'ISU')read(11,*)isu
       if(key(1:3).eq.'ITE')read(11,*)lit
       if(key(1:3).eq.'MAX')read(11,*)maxit
       if(key(1:2).eq.'NB')read(11,*)NB
       if(key(1:4).eq.'NCUT')read(11,*)ncut
       if(key(1:2).eq.'NP')read(11,*)NP
       if(key(1:3).eq.'NSP')read(11,*)nps
       if(key(1:5).eq.'LEXCL')read(11,*)LEXCL
       if(key(1:4).eq.'KFAC')read(11,*)kfac
       if(key(1:4).eq.'LCOM')read(11,*)lcom
       if(key(1:4).eq.'WFIX')read(11,*)wfix
       if(key(1:3).eq.'LFI')read(11,*)lfix
       if(key(1:5).eq.'LINEA')read(11,*)llin
       if(key(1:3).eq.'LWR')read(11,*)lwr
       if(key(1:3).eq.'ONL')read(11,*)lonly
       if(key(1:3).eq.'SCA')read(11,*)lmscan
       if(key(1:5).eq.'LSFIX')read(11,*)lsfix
       if(key(1:4).eq.'LCRE')read(11,*)lcre
       if(key(1:6).eq.'SHIFTQ')read(11,*)lshiftq
       if(key(1:6).eq.'SHIFTI')read(11,*)lshifti
       if(key(1:6).eq.'SHIFTE')read(11,*)lshifte
       if(key(1:5).eq.'SPECT')read(11,*)lspec
       if(key(1:3).eq.'TOL')read(11,*)tol
       if(key(1:3).eq.'WLI')read(11,*)wlim
       if(key(1:3).eq.'WMI')read(11,*)wmin
       if(key(1:5).eq.'WSMIN')read(11,*)wsmin
       if(key(1:5).eq.'WSMAX')read(11,*)wsmax
       goto 1
11     close(11)
      endif
      return
      end
c     ==============================================================
      subroutine fixmodes(N,NQ,NE,g,sg,p,se,amas,z,ifx,wlim,lcor,lonly,
     1sgf,NF,wfix)
      IMPLICIT none
      integer*4 N,NE,NQ,I,iz,ifx,j,ix,no,jmin,nblo,z(*),NU,NF
      real*8 g(*),p(*),sg(N,N),se(N,N),CM,wlim,dij,dmin,t,sgf(N,N),
     1wfix
      real amas(*)
      integer*4,allocatable::bl(:),be(:)
      character*3 key
      logical lex,lcor,lonly
      CM=219474.630d0
      no=NQ

      write(6,6009)ifx,NQ
6009  format(/,' Fixmodes',/,' IFIX = ',I4,/,' NQ = ',I4)

      if(NQ.ne.NE)then
       NU=min(NE,NQ)
       write(6,60091)NU
60091  format(' NE<>NQ, take minimum ',i4,' modes')
       NE=NU
       NQ=NU
      endif

      if(lcor)then
       write(6,*)'Correlating ground and excited state modes'
       if(lonly)then
        write(6,*)'Correlation only, no action'
        do 121 i=NQ,1,-1
        call findbest(i,NQ,jmin,dmin,dij,sg,se,N,amas,z)
121     write(6,609)i,jmin,dmin
609     format(' Ground mode',i4,' corresponds to',i4,' s =',f10.4)
        write(6,*)
        do 122 i=NQ,1,-1
        call findbest(i,NQ,jmin,dmin,dij,se,sg,N,amas,z)
122     write(6,6091)i,jmin,dmin
6091    format(' Exc. mode',i4,' corresponds to',i4,' s =',f10.4)
       else
c       loop starting from higher modes that migh be more reliable:
        do 12 i=NQ,1,-1
c       find best match for i within i...1:
        call findbest(i,i,jmin,dmin,dij,se,sg,N,amas,z)
        if(i.ne.jmin)then
         t=g(i)
         g(i)=g(jmin)
         g(jmin)=t
         do 15 ix=1,N
         t=sg(ix,i)
         sg(ix,i)=sg(ix,jmin)
15       sg(ix,jmin)=t
         write(6,619)i,jmin,dmin
619      format(' Exc. mode',i4,' corresponds to',i4,' s =',f10.4,
     1   ' modes switched')
        else
         write(6,6091)i,jmin,dmin
        endif
12      continue
       endif
      endif
      
c     ifx = 0   do nothing
c           1   change negative frequencies to wfix
c           2   delete negatives
c           3   borrow exc freqs from ground and vice versa
c           4   delete smaller than a limit
c
c     how many modes deleted:
      NF=0

      if(ifx.eq.1)then
       iz=0
       do 2 I=1,NQ
       if(g(I).lt.0.0d0)then
        g(I)=wfix/CM
        iz=iz+1
       endif
       if(p(I).lt.0.0d0)then
        p(I)=wfix/CM
        iz=iz+1
       endif
 2     continue
       if(iz.gt.0)write(6,46)iz,wfix
46     format(i6,' neg. modes made',f12.2,' cm-1 ')
      endif

      if(ifx.eq.2)then
       iz=0
66     do 3 I=1,NQ
       if(g(I).lt.0.0d0.or.p(i).lt.0.0d0)then
c       save deleted sg vector:
        NF=NF+1
        do 41 ix=1,N
41      sgf(ix,NF)=sg(ix,i)
        do 4 j=i,NQ-1
        g(j)=g(j+1)
        p(j)=p(j+1)
        do 4 ix=1,N
        sg(ix,j)=sg(ix,j+1)
4       se(ix,j)=se(ix,j+1)
        iz=iz+1
        NQ=NQ-1
        goto 66
       endif
 3     continue
       if(iz.gt.0)write(6,*)iz,' neg. modes deleted'
      endif

      if(ifx.eq.3)then
       iz=0
       do 5 I=1,NQ
       if(p(i).lt.0.0d0)then
        if(g(I).gt.0.0d0)then
         p(i)=g(i)
         iz=iz+1
        else
         write(6,*)I
         call report(' mode cannot be fixed')
        endif
       endif
       if(g(i).lt.0.0d0)then
        if(p(I).gt.0.0d0)then
         g(i)=p(i)
         iz=iz+1
        else
         write(6,*)I
         call report(' mode cannot be fixed')
        endif
       endif
5      continue
       if(iz.gt.0)write(6,*)iz,' neg. modes borrowed'
      endif

      if(ifx.eq.4)then
       iz=0
67     do 6 i=1,NQ
       if(g(i)*CM.lt.wlim.or.p(i)*CM.lt.wlim)then
c       save deleted sg vector:
        NF=NF+1
        do 71 ix=1,N
71      sgf(ix,NF)=sg(ix,i)
        do 7 j=i,NQ-1
        g(j)=g(j+1)
        p(j)=p(j+1)
        do 7 ix=1,N
        sg(ix,j)=sg(ix,j+1)
7       se(ix,j)=se(ix,j+1)
        iz=iz+1
        NQ=NQ-1
        goto 67
       endif
 6     continue
       if(iz.gt.0)write(6,*)iz,' small modes skipped'
      endif

      nblo=0
      inquire(file='DUSCH.OPT',exist=lex)
      if(lex)then
       open(11,file='DUSCH.OPT')
1      read(11,900,end=11,err=11)key
 900   format(a3)
       if(key.eq.'BLO')then
        read(11,*)nblo
        allocate(bl(abs(nblo)),be(abs(nblo)))
        read(11,*)bl
        if(nblo.lt.0)read(11,*)be
       endif
       goto 1
11     close(11)
      endif

      if(nblo.ne.0)then
       call dm(N,NQ,abs(nblo),bl,sg,g)
       if(nblo.gt.0)then
c       blocked same modes in excited:
        call dm(N,NQ,abs(nblo),bl,se,p)
       else
c       special list for excited:
        call dm(N,NQ,abs(nblo),be,se,p)
       endif
       NQ=NQ-abs(nblo)
       if(NQ.ne.no)write(6,*)NQ,' modes now'
      endif

      return
      end

      subroutine dm(N,NQO,nblo,b,s,e)
      implicit none
      integer*4 nblo,N,NQ,i,j,ix,b(nblo),NQO
      real*8 s(N,N),e(N)
      NQ=NQO
      do 1 i=1,NQ
      do 1 j=1,nblo
1     if(i.eq.b(j))e(i)=-3.333d0
68    do 2 i=1,NQ
      if(e(i).eq.-3.333d0)then
       do 3 j=i,NQ-1
       e(j)=e(j+1)
       do 3 ix=1,N
3      s(ix,j)=s(ix,j+1)
       NQ=NQ-1
       goto 68
      endif
2     continue
      write(6,*)nblo,' modes blocked (deleted)'
      return
      end

      SUBROUTINE readsi(N3,S,E,NQ,fn,r,ldz)
c     switch order of modes to match Gaussian convention:
      IMPLICIT none
      integer*4 N3,NQ,nat,i,J,ix,iz
      real*8 S(N3,N3),E(*),CM,r(*)
      logical ldz
      character*(*) fn
      CM=219474.630d0
C     1 a.u.= 2.1947E5 cm^-1
      open(4,file=fn)
      read(4,*)NQ,nat,nat
      do 1 i=1,NAT
1     read(4,*)r(3*(i-1)+1),(r(3*(i-1)+ix),ix=1,3)
      read(4,*)
      DO 2 I=1,NAT
      DO 2 J=1,NQ
2     read(4,*)(s(3*(i-1)+ix,NQ-J+1),ix=1,2),
     1(s(3*(i-1)+ix,NQ-J+1),ix=1,3)
      read(4,*)
      READ(4,4000)(E(NQ-J+1),J=1,NQ)
4000  FORMAT(6F11.6)
      close(4)
c
      write(6,*)NQ,' modes found'

c     delete zero modes if exist:
      if(ldz)then
       iz=0
66     do 6 i=1,NQ
       if(dabs(E(i)).lt.0.1d0)then
        do 7 j=i,NQ-1
        E(j)=E(j+1)
        do 7 ix=1,N3
7       s(ix,j)=s(ix,j+1)
        iz=iz+1
        NQ=NQ-1
        goto 66
       endif
6      continue
       if(iz.gt.0)write(6,*)iz,' zero modes: deleted'
      endif

      write(6,*)NQ,' vibrational modes considered'

      DO 3 I=1,NQ
3     E(I)=E(I)/CM
      
      RETURN
      end
c     ==============================================================
      SUBROUTINE writesi(NAT,N3,S,E,NQ,fn,x,z)
c     switch order of modes to match Gaussian convention:
      IMPLICIT none
      integer*4 N3,NQ,nat,i,J,ix,z(*)
      real*8 S(N3,N3),E(*),CM,x(*)
      character*(*) fn
      CM=219474.630d0
C     1 a.u.= 2.1947E5 cm^-1
      open(4,file=fn)
      write(4,*)NQ,nat,nat
      do 1 i=1,NAT
1     write(4,401)z(i),(x(ix+3*(i-1)),ix=1,3)
401   format(i6,3f12.6)
      write(4,*)'Atom mode x y z'
      DO 2 I=1,NAT
      DO 2 J=1,NQ
2     write(4,402)I,J,(s(3*(i-1)+ix,NQ-J+1),ix=1,3)
402   format(2i6,3f12.6)
      write(4,402)NQ
      write(4,4000)(E(NQ-J+1)*CM,J=1,NQ)
4000  FORMAT(6F11.3)
      close(4)
      RETURN
      end
c     ==============================================================
      subroutine xst(nat,rgo,re,m,z,t)
c     standard orientation of re with respect to rg
c     J. Phys. Chem. A 2001, 105, 5326-5333
      implicit none
      integer*4 nat,i,j,ia,IERR,ix,iy,iz,ii,ib,ic,z(*)
      real*8 rgo(*),re(*),c(3,3),m(*),cc(3,3),r(3,3),l(3),ci(3,3),
     1TOL,e(3,6),lrc(3,3),bl(3),rc(3,3),err,t0(3,3),det,t(3,3),em,
     1cmg(3),cme(3),mass
      real*8,allocatable::rn(:),rg(:)
      TOL=1.0d-10
      allocate(rg(3*nat))
c     calculate mass centers:
      do 14 i=1,3
      cmg(i)=0.0d0
      cme(i)=0.0d0
      mass=0.0d0
      do 16 ia=1,nat
      mass=mass+m(ia)
      cmg(i)=cmg(i)+m(ia)*rgo(i+3*(ia-1))
16    cme(i)=cme(i)+m(ia)* re(i+3*(ia-1))
      cmg(i)=cmg(i)/mass
14    cme(i)=cme(i)/mass
      write(6,603)cmg,cme
603   format(' Mass centers: g: ',3f12.4,/,
     1       '               e: ',3f12.4)
c     mass-center coordinates:
      do 15 ia=1,nat
      do 15 i=1,3
      rg(i+3*(ia-1))=rgo(i+3*(ia-1))-cmg(i)
15    re(i+3*(ia-1))= re(i+3*(ia-1))-cme(i)
c     form C =X M X':
      do 1 i=1,3
      do 1 j=1,3
      c(i,j)=0.0d0
      do 1 ia=1,nat
1     c(i,j)=c(i,j)+rg(i+3*(ia-1))*m(ia)*re(j+3*(ia-1))
c     correction for planar molecules:
      do 2 i=1,3
2     if(dabs(c(i,i)).lt.1.0d-9)c(i,i)=1.0d0
c     C^TC:
      call mm3(c,c,cc)
      CALL TRED12(3,cc,r,l,2,IERR)
      IF(IERR.NE.0)call report(' cannot diagonalize C^TC')
      call INV(3,c,ci,3,TOL,e,IERR)
      IF(IERR.NE.0)call report(' cannot invert C')
c     R^TC^-1:
      call mm3(r,ci,rc)
c     temporary coordinate string:
      allocate(rn(3*nat))
c     exlore all eight axis switching possibilities:
      ib=0
      ic=0
      do 4 ix=-1,1,2
      do 4 iy=-1,1,2
      do 4 iz=-1,1,2
      ic=ic+1
      bl(1)=dble(ix)
      bl(2)=dble(iy)
      bl(3)=dble(iz)
c     L l^(1/2)R^TC^-1:
      do 6 i=1,3
      do 6 j=1,3
6     lrc(i,j)=bl(i)*dsqrt(l(i))*rc(i,j)
      call md(r,lrc,t0)
      call trx(nat,rn,re,t0)
      err=0.0d0
      do 8 ii=1,3*nat
8     err=err+(rn(ii)-rg(ii))**2
      if(ic.eq.1)then
       em=err
       do 9 i=1,3
       do 9 j=1,3
9      t(i,j)=t0(i,j)
       ib=ic
      endif
      if(det(t0).gt.0.0d0.and.err.lt.em)then
       em=err
       do 11 i=1,3
       do 11 j=1,3
11     t(i,j)=t0(i,j)
       ib=ic
      endif
4     write(6,600)ix,iy,iz,err,det(t0)
600   format(3i2,'err:',f10.3,' det:',f10.3)
      write(6,601)ib
601   format(' Best trial number ',i2)
      write(6,604)t
604   format(' Transformation matrix:',/,3(3f10.3,/))
      open(44,file='TMX.TXT')
      do 17 i=1,3
17    write(44,605)(t(j,i),j=1,3)
605   format(3f10.3)
      close(44)
      call trx(nat,rn,re,t)
      do 131 ia=1,nat
      do 131 i=1,3
131   re(i+3*(ia-1))=cmg(i)+rn(i+3*(ia-1))
      call wrx(nat,re,z,'exc.as.ground.x',
     1'Excited orientation to match the ground:')
      write(6,*)'exc.as.ground.x written'
      return
      end
c     ==============================================================
      subroutine wrx(nat,r,z,s,sp)
      implicit none
      integer*4 nat,z(nat),a,i
      real*8 r(3*nat)
      character*(*) s,sp
      open(44,file=s)
      write(44,'(a)')sp
      write(44,*)nat
      do 13 a=1,nat
13    write(44,602)z(a),(r(i+3*(a-1)),i=1,3)
602   format(i4,3f12.6)
      close(44)
      return
      end
c     ==============================================================
      subroutine trx(nat,n,e,t)
      implicit none
      integer*4 nat,ia,i,a
      real*8 n(*),e(*),t(3,3)
      do 7 ia=1,nat
      a=3*(ia-1)
      do 7 i=1,3
7     n(i+a)=t(1,i)*e(1+a)+t(2,i)*e(2+a)+t(3,i)*e(3+a)
      return
      end
c     ==============================================================
      function det(t)
      implicit none
      real*8 det,t(3,3)
      det=t(1,1)*(t(2,2)*t(3,3)-t(2,3)*t(3,2))
     1   -t(1,2)*(t(2,1)*t(3,3)-t(2,3)*t(3,1))
     1   +t(1,3)*(t(2,1)*t(3,2)-t(2,2)*t(3,1))
      return
      end
c     ==============================================================
      subroutine md(a,b,c)
      implicit none
c     c=a.b
      real*8 a(3,3),b(3,3),c(3,3)
      integer*4 i,j,ii
      do 3 i=1,3
      do 3 j=1,3
      c(i,j)=0.0d0
      do 3 ii=1,3
3     c(i,j)=c(i,j)+a(i,ii)*b(ii,j)
      return
      end
c     ==============================================================
      subroutine mm3(a,b,c)
      implicit none
c     c=aT.b
      real*8 a(3,3),b(3,3),c(3,3)
      integer*4 i,j,ii
      do 3 i=1,3
      do 3 j=1,3
      c(i,j)=0.0d0
      do 3 ii=1,3
3     c(i,j)=c(i,j)+a(ii,i)*b(ii,j)
      return
      end
c     ==============================================================
      SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(N,N),Z(N,N),D(*)
      real*8,allocatable::E(:)
      allocate(E(N))
      IF (IEIG .LT. 0) GO TO 110
      DO 100 I = 1, N
C
      DO 100 J = 1, I
      Z(I,J) = A(I,J)
  100 CONTINUE
C
  110 ICODE = IABS(IEIG)
      IF (N .EQ. 1) GO TO 320
C     :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- ::::::::::
      DO 300 II = 2, N
      I = N + 2 - II
      L = I - 1
      H = 0.0D0
      SCALE = 0.0D0
      IF (L .LT. 2) GO TO 130
C     :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) ::::::::::
      DO 120 K = 1, L
  120 SCALE = SCALE + DABS(Z(I,K))
C
      IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140
  130 E(I) = Z(I,L)
      GO TO 290
C
  140 DO 150 K = 1, L
      F = Z(I,K) / SCALE
      Z(I,K) = F
      H = H + F * F
  150 CONTINUE
C
      G = -DSIGN(DSQRT(H),F)
      E(I) = SCALE * G
      H = H - F * G
      Z(I,L) = F - G
      F = 0.0D0
C
      DO 240 J = 1, L
      IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H)
      G = 0.0D0
C     :::::::::: FORM ELEMENT OF A*U ::::::::::
      DO 180 K = 1, J
  180 G = G + Z(J,K) * Z(I,K)
C
      JP1 = J + 1
      IF (L .LT. JP1) GO TO 220
C
      DO 200 K = JP1, L
  200 G = G + Z(K,J) * Z(I,K)
C     :::::::::: FORM ELEMENT OF P ::::::::::
  220 E(J) = G / H
      F = F + E(J) * Z(I,J)
  240 CONTINUE
C
      HH = F / (H + H)
C     :::::::::: FORM REDUCED A ::::::::::
      DO 260 J = 1, L
      F = Z(I,J)
      G = E(J) - HH * F
      E(J) = G
C
      DO 260 K = 1, J
      Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K)
  260 CONTINUE
C
      DO 280 K = 1, L
  280 Z(I,K) = SCALE * Z(I,K)
C
  290 D(I) = H
  300 CONTINUE
C
  320 D(1) = 0.0D0
      E(1) = 0.0D0
C     :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES ::::::::::
      IF (ICODE .NE. 2) GO TO 600
      DO 500 I = 1, N
      L = I - 1
      IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380
C
      DO 360 J = 1, L
      G = 0.0D0
C
      DO 340 K = 1, L
  340 G = G + Z(I,K) * Z(K,J)
C
      DO 360 K = 1, L
      Z(K,J) = Z(K,J) - G * Z(K,I)
  360 CONTINUE
C
  380 D(I) = Z(I,I)
      Z(I,I) = 1.0D0
      IF (L .LT. 1) GO TO 500
C
      DO 400 J = 1, L
      Z(I,J) = 0.0D0
      Z(J,I) = 0.0D0
  400 CONTINUE
C
  500 CONTINUE
C
  620 CALL TQL12 (N,Z,D,IERR,ICODE,E)
      RETURN
C     :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY ::::::::::
  600 DO 610 I=1,N
  610 D(I) = Z(I,I)
      GO TO 620
      END
C
      SUBROUTINE TQL12(N,Z,D,IERR,ICODE,E)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION Z(N,N),D(*),E(*)
      REAL*8 MACHEP
      EPS = 1.0D0
   10 EPS = 0.50D0*EPS
      TOL1 = EPS + 1.0D0
      IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10
      IF(TOL1-EPS.EQ.1.0D0)  EPS = EPS + EPS
      MACHEP = EPS
C
C     MACHEP=16.0D0**(-13)
C
      IERR = 0
      IF (N .EQ. 1) GO TO 1001
C
      DO 100 I = 2, N
  100 E(I-1) = E(I)
C
      F = 0.0D0
      B = 0.0D0
      E(N) = 0.0D0
C
      DO 240 L = 1, N
      J = 0
      H = MACHEP * (DABS(D(L)) + DABS(E(L)))
      IF (B .LT. H) B = H
C     :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT ::::::::::
      DO 110 M = L, N
      IF (DABS(E(M)) .LE. B) GO TO 120
C     :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
C                THROUGH THE BOTTOM OF THE LOOP ::::::::::
  110 CONTINUE
C
  120 IF (M .EQ. L) GO TO 220
  130 IF (J .EQ. 30) GO TO 1000
      J = J + 1
C     :::::::::: FORM SHIFT ::::::::::
      L1 = L + 1
      G = D(L)
      P = (D(L1) - G) / (2.0D0 * E(L))
      R = DSQRT(P*P+1.0D0)
      D(L) = E(L) / (P + DSIGN(R,P))
      H = G - D(L)
C
      DO 140 I = L1, N
  140 D(I) = D(I) - H
C
      F = F + H
C     :::::::::: QL TRANSFORMATION ::::::::::
      P = D(M)
      C = 1.0D0
      S = 0.0D0
      MML = M - L
C     :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- ::::::::::
      DO 200 II = 1, MML
      I = M - II
      G = C * E(I)
      H = C * P
      IF (DABS(P) .LT. DABS(E(I))) GO TO 150
      C = E(I) / P
      R = DSQRT(C*C+1.0D0)
      E(I+1) = S * P * R
      S = C / R
      C = 1.0D0 / R
      GO TO 160
  150 C = P / E(I)
      R = DSQRT(C*C+1.0D0)
      E(I+1) = S * E(I) * R
      S = 1.0D0 / R
      C = C * S
  160 P = C * D(I) - S * G
      D(I+1) = H + S * (C * G + S * D(I))
C     :::::::::: FORM VECTOR ::::::::::
      IF (ICODE .NE. 2) GO TO 200
      DO 180 K = 1, N
      H = Z(K,I+1)
      Z(K,I+1) = S * Z(K,I) + C * H
      Z(K,I) = C * Z(K,I) - S * H
  180 CONTINUE
C
  200 CONTINUE
C
      E(L) = S * P
      D(L) = C * P
      IF (DABS(E(L)) .GT. B) GO TO 130
  220 D(L) = D(L) + F
  240 CONTINUE
C     :::::::::: ORDER EIGENVALUES AND EIGENVECTORS ::::::::::
      DO 300 II = 2, N
      I = II - 1
      K = I
      P = D(I)
C
      DO 260 J = II, N
      IF (D(J) .LE. P) GO TO 260
      K = J
      P = D(J)
  260 CONTINUE
C
      IF (K .EQ. I) GO TO 300
      D(K) = D(I)
      D(I) = P
C
      IF (ICODE .NE. 2) GO TO 300
      DO 280 J = 1, N
      P = Z(J,I)
      Z(J,I) = Z(J,K)
      Z(J,K) = P
  280 CONTINUE
C
  300 CONTINUE
C
      GO TO 1001
C     :::::::::: SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30 ITERATIONS ::::::::::
 1000 IERR = L
 1001 RETURN
      END
C
      SUBROUTINE INV(nr,a,ai,n,TOL,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
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
c     ============================================================
      subroutine makeG(NQ,G,GP,GH,GHP,eg,ep)
      implicit none
      real*8 eg(*),ep(*),G(NQ,NQ),GP(NQ,NQ),GH(NQ,NQ),GHP(NQ,NQ)
      integer*4 ix,NQ
      call mz(G,NQ)
      call mz(GP,NQ)
      call mz(GH,NQ)
      call mz(GHP,NQ)
      do 10 ix=1,NQ
      G(ix,ix)=eg(ix)
      GP(ix,ix)=ep(ix)
      GH(ix,ix) =dsqrt(eg(ix))
10    GHP(ix,ix)=dsqrt(ep(ix))
      return
      end
c     ============================================================
      function overlap(ic,NQ,N,nat,G,GP,GH,GHP,eg,ep,JT,J,sg,se,
     1m,k,rg,re,bohr,T,TP,F,FI,E2,TOL,IERR,A,TV,TU,B,C,D,E,kfac)
c     calculates <0|0'>
c     ic  ... 0   only the overlap
c             1   also write DUSCH.OUT
      implicit none
      integer*4 ic,NQ,N,ig,ie,ia,nat,ix,jx,ii,IERR
      real*8 eg(*),ep(*),G(NQ,NQ),GP(NQ,NQ),GH(NQ,NQ),GHP(NQ,NQ),
     1J(NQ,NQ),JT(NQ,NQ),sg(N,N),se(N,N),m(*),k(*),rg(*),re(*),
     2bohr,T(NQ,NQ),TP(NQ,NQ),F(NQ,NQ),FI(NQ,NQ),E2(NQ,2*NQ),TOL,
     3A(NQ,NQ),TV(*),TU(*),B(*),C(NQ,NQ),D(*),E(NQ,NQ),overlap,
     4fac,kfac

      call makeG(NQ,G,GP,GH,GHP,eg,ep)

c     K = s'^-1 sqrt(M) dX      
      do 4 ig=1,NQ
      k(ig)=0
      do 4 ia=1,nat
      do 4 ix=1,3
      ii=ix+3*(ia-1)
4     k(ig)=k(ig)+sg(ii,ig)*(re(ii)-rg(ii))/bohr
     1*m(ia)*dsqrt(1822.89d0)

c     Duschinsky transformation:
c     Q' = JQ'' + K    (Q' initial, Q'' final)

c     J = s'^-1 s'', JT = J^t
      do 3 ig=1,NQ
      do 3 ie=1,NQ
      JT(ie,ig)=0.0d0
      do 3 ia=1,nat
      do 3 ix=1,3
      ii=ix+3*(ia-1)
3     JT(ie,ig)=JT(ie,ig)+sg(ii,ig)*m(ia)*se(ii,ie)

c     scale the shift if kfac <> 1:
      do 41 ig=1,NQ
41    k(ig)=k(ig)*kfac

      call mt(J,JT,NQ)
      call mm(T,G,J,NQ)
      call mm(TP,JT,T,NQ)
      call ms(F,TP,GP,NQ)
c     F = Jt G J + GP
      call INV(NQ,F,FI,NQ,TOL,E2,IERR)
      if(IERR.ne.0)call report('Inversion error')
      call mm(T,JT,GH,NQ)
      call mm(TP,FI,T,NQ)
      call mm(T,J,TP,NQ)
      call mm(TP,GH,T,NQ)
      do 11 ix=1,NQ
      do 11 jx=1,NQ
11    A(ix,jx)=2.0d0*TP(ix,jx)
      do 12 ix=1,NQ
12    A(ix,ix)=A(ix,ix)-1.0d0
c     A=2 G^1/2 J Fi Jt G^1/2 - I

      call mm(T,JT,G,NQ)
      call mm(TP,FI,T,NQ)
      call mm(T,J,TP,NQ)
      do 13 ix=1,NQ
13    T(ix,ix)=T(ix,ix)-1.0d0
      call mv(TV,T,K,NQ)
      call mv(TU,GH,TV,NQ)
      do 14 ix=1,NQ
14    B(ix)=-2.0d0*TU(ix)
c     B= - 2 G^1/2 ( J Fi Jt G - I) K

      call mm(T,FI,GHP,NQ)
      call mm(TP,GHP,T,NQ)
      do 15 ix=1,NQ
      do 15 jx=1,NQ
15    C(ix,jx)=2.0d0*TP(ix,jx)
      do 16 ix=1,NQ
16    C(ix,ix)=C(ix,ix)-1.0d0
c     C= 2 GP^1/2 Fi GP^1/2 - I

      call mv(TV,G,K,NQ)
      call mv(TU,JT,TV,NQ)
      call mv(TV,FI,TU,NQ)
      call mv(TU,GHP,TV,NQ)
      do 17 ix=1,NQ
17    D(ix)=-2.0d0*TU(ix)
c     D= - 2 GP^1/2 Fi Jt G K

      call mm(T,JT,GH,NQ)
      call mm(TP,FI,T,NQ)
      call mm(T,GHP,TP,NQ)
      do 18 ix=1,NQ
      do 18 jx=1,NQ
18    E(ix,jx)=4.0d0*T(ix,jx)
c     E= 4 GP^1/2 Fi Jt G^1/2

      write(6,6555)'A ',1,1,NQ,NQ,A(1,1),A(NQ,NQ)
      write(6,6555)'B ',1,1,NQ,NQ,B(1),B(NQ)
      write(6,6555)'C ',1,1,NQ,NQ,C(1,1),C(NQ,NQ)
      write(6,6555)'D ',1,1,NQ,NQ,D(1),D(NQ)
      write(6,6555)'E ',1,1,NQ,NQ,E(1,1),E(NQ,NQ)
      write(6,6555)'F ',1,1,NQ,NQ,F(1,1),F(NQ,NQ)
      write(6,6555)'Fi',1,1,NQ,NQ,FI(1,1),FI(NQ,NQ)
      write(6,6555)'G ',1,1,NQ,NQ,G(1,1),G(NQ,NQ)
      write(6,6555)'GP',1,1,NQ,NQ,GP(1,1),GP(NQ,NQ)
      write(6,6555)'J ',1,1,NQ,NQ,J(1,1),J(NQ,NQ)
      write(6,6555)'JT',1,1,NQ,NQ,JT(1,1),JT(NQ,NQ)
      write(6,6555)'K ',1,1,NQ,NQ,K(1),K(NQ)
6555  format(1x,a2,2i3,1x,2i3,2f12.6)

      call mm(T,G,GP,NQ)

      call dofac(NQ,fac,J,JT,F,FI,G,GP,K)
      overlap=fac

      if(ic.eq.1)
     1call wdusch(NQ,N,A,B,C,D,E,JT,K,eg,ep,sg,se,rg,re,m)
      
      return
      end
c     ============================================================
      subroutine iterate(maxit,NQ,lwr,wmin,
     1N,nat,G,GP,GH,GHP,eg,ep,JT,J,sg,se,
     1m,k,rg,re,bohr,T,TP,F,FI,E2,TOL,IERR,A,TV,TU,B,C,D,E,kfac)
      implicit none
      integer*4 iter,maxit,iw,NQ,ichange,IERR,N,nat
      real*8 dd,wp,wm,o,om,op,w,CM,wmin,delw,oold,kfac
      real*8 eg(*),ep(*),G(NQ,NQ),GP(NQ,NQ),GH(NQ,NQ),GHP(NQ,NQ),
     1J(NQ,NQ),JT(NQ,NQ),sg(N,N),se(N,N),m(*),k(*),rg(*),re(*),
     2bohr,T(NQ,NQ),TP(NQ,NQ),F(NQ,NQ),FI(NQ,NQ),E2(NQ,2*NQ),TOL,
     3A(NQ,NQ),TV(*),TU(*),B(*),C(NQ,NQ),D(*),E(NQ,NQ),overlap
      logical lwr
      CM=219474.630d0
      dd=0.5d0
      oold=1.0d0
      do 1 iter=1,maxit
      ichange=0
      if(lwr)write(6,*)'Iteration ',iter,' of ',maxit
      if(lwr)write(6,601)
601   format(3x,'iw',11x,'w',10x,'wp',10x,'wm',11x,'o',10x,'op',
     110x,'om')
      do 2 iw=1,NQ
      w=ep(iw)
      delw=max(wmin/CM,dd*w)
      o=overlap(0,NQ,N,nat,G,GP,GH,GHP,eg,ep,JT,J,sg,se,
     1m,k,rg,re,bohr,T,TP,F,FI,E2,TOL,IERR,A,TV,TU,B,C,D,E,kfac)
      wp=w+delw
      ep(iw)=wp
      op=overlap(0,NQ,N,nat,G,GP,GH,GHP,eg,ep,JT,J,sg,se,
     1m,k,rg,re,bohr,T,TP,F,FI,E2,TOL,IERR,A,TV,TU,B,C,D,E,kfac)
      wm=max(wmin/CM,w-delw)
      ep(iw)=wm
      om=overlap(0,NQ,N,nat,G,GP,GH,GHP,eg,ep,JT,J,sg,se,
     1m,k,rg,re,bohr,T,TP,F,FI,E2,TOL,IERR,A,TV,TU,B,C,D,E,kfac)
      if(lwr)write(6,600)iw,w*CM,wp*CM,wm*CM,o,op,om
600   format(i5,3f12.3,3g12.3)
      if(om.gt.o.or.op.gt.o)then
       if(om.gt.op)then
        ep(iw)=wm
       else
        ep(iw)=wp
       endif
       ichange=ichange+1
      else
       ep(iw)=w
      endif
2     continue
      if(ichange.eq.0)dd=dd/2.0d0
      if(oold.eq.o)return
1     oold=o
      return
      end
c     ============================================================
      function spsv(sg,se,i,j,amas,z,N)
      implicit none
      integer*4 N,z(*),i,j,ix
      real*8 spsv,sg(N,N),se(N,N),dij
      real amas(*)
      dij=0.0d0
      do 1 ix=1,N
1     dij=dij+sg(ix,i)*se(ix,j)*dble(amas(z((ix-1)/3+1)))
      spsv=dij
      return
      end
c     ============================================================
      subroutine findbest(i,jstart,jmin,dmin,dij,sg,se,N,amas,z)
      implicit none
      integer*4 i,jstart,jmin,N,z(*),j
      real*8 dmin,dij,sg(N,N),se(N,N),spsv
      real amas(*)
c     find best match for i within jstart...1
      jmin=0
      dmin=0.0d0
      dij=0.0d0
      do 14 j=jstart,1,-1
      dij=spsv(sg,se,i,j,amas,z,N)
      if(dabs(dij).gt.dmin.or.j.eq.jstart)then
       dmin=dabs(dij)
       jmin=j
      endif
14    continue
      return
      end
c     ============================================================
      subroutine gradq(nat,N,m,ldz)
      implicit none
      integer*4 nat,ia,ix,N,NQ,iq,i,j
      real*8 SFAC,CM,wep,a,ap,x0,ef,o,m(*)
      logical lexg,lexe,ldz
      real*8,allocatable::g(:),q(:),u(:,:),qe(:),qt(:),r(:),
     1sg(:,:),se(:,:),eg(:),ep(:)
      allocate(sg(N,N),se(N,N),eg(N),ep(N))
      SFAC=0.0234280d0
      CM=219474.0d0
      inquire(file='ground/FILE.GR',exist=lexg)
      if(lexg)then
       write(6,*)
       write(6,*)' Ground state gradient found'
       write(6,*)' GG.TXT opened'
       open(66,file='GG.TXT')
       allocate(g(N),q(N),r(N))
       open(8,file='ground/FILE.GR')
       do 1 ia=1,nat
1      read(8,*)(g(ix+3*(ia-1)),ix=1,3)
       close(8)
       call readsi(N,sg,eg,NQ,'ground/F.INP',r,ldz)
       do 2 iq=1,NQ
       q(iq)=0.0d0
       do 2 ia=1,nat
       do 2 ix=1,3
2      q(iq)=q(iq)+sg(ix+3*(ia-1),iq)*g(ix+3*(ia-1))*SFAC
       write(66,*)' ground state:'
       write(66,*)' mode        wg gradient (au)'
       do 3 iq=1,NQ
3      write(66,600)iq,eg(iq)*CM,q(iq)
600    format(i4,f12.2,f12.7)
       inquire(file='excited/FILE.GR',exist=lexe)
       if(lexe)then
        write(6,*)' Excited state gradient found, too'
        allocate(qe(N),qt(N))
        open(8,file='excited/FILE.GR')
        do 4 ia=1,nat
4       read(8,*)(g(ix+3*(ia-1)),ix=1,3)
        close(8)
c       reload because se might be transformed:
        call readsi(N,se,ep,NQ,'excited/F.INP',r,ldz)

        do 5 iq=1,NQ
        qe(iq)=0.0d0
        do 5 ia=1,nat
        do 5 ix=1,3
5       qe(iq)=qe(iq)+se(ix+3*(ia-1),iq)*g(ix+3*(ia-1))*SFAC
        write(66,*)' excited state:'
        write(66,*)' mode        we gradient (au)'
        do 6 iq=1,NQ
6       write(66,600)iq,ep(iq)*CM,qe(iq)
        write(6,*)' .... gradient written'
        write(66,*)' Independent mode approximation:'
        write(66,*)' Excited state in ground state modes:'
c       u=se.s transformation matrix between ground and excited modes:
        allocate(u(NQ,NQ))
        do 7 i=1,NQ
        do 7 j=1,NQ
        u(i,j)=0.0d0
        do 7 ia=1,nat
        do 7 ix=1,3
7       u(i,j)=u(i,j)+sg(ix+3*(ia-1),i)*m(ia)*se(ix+3*(ia-1),j)
        write(66,602)
602     format(' mode         wg         wep gradient (au)'
     1  //' q0    <f00|F0e>')
        do 8 iq=1,NQ
        qt(iq)=0.0d0
c       diagonal part of Wij:
        wep=0.0d0
        do 9 i=1,NQ
        qt(iq)=qt(iq)+u(iq,i)*qe(i)
9       wep   =wep   +u(iq,i)*ep(i)*u(iq,i)
        a =eg(iq)/2.0d0
        ap=wep   /2.0d0
        x0=-qe(iq)/wep**2
        if(wep.lt.0.0d0)x0=999999999.99d0
        if(a.gt.0.0d0.and.ap.gt.0.0d0)then
         ef=x0**2*a*ap/(a+ap)
         if(ef.gt.20.0d0)then
          o=0.0d0
         else
          o=dsqrt(2.0d0*dsqrt(a*ap)/(a+a))*exp(-ef)
         endif
        else
         o=-999.99d0
        endif
8       write(66,601)iq,eg(iq)*CM,wep*CM,qe(iq),x0,o
601     format(i4,2f12.2,3g12.4)
        write(6,*)' .... transformed into ground modes'
       endif
       close(66)
       write(6,*)' GG.TXT closed'
       write(6,*)
      endif
      return
      end
c     ============================================================
      subroutine grade(nat,N,ldz)
      implicit none
      integer*4 nat,ia,ix,N,NQ,iq,ii
      real*8 SFAC,CM,y,ani,w,bohr
      logical lex,ldz
      real*8,allocatable::g(:),q(:),rg(:),x(:),qt(:),ee(:),r(:),
     1sg(:,:),eg(:)
      integer*4 ,allocatable::z(:)
      allocate(sg(N,N),eg(N))
      SFAC=0.0234280d0
      CM=219474.0d0
      y=50.0d0/CM
      bohr=0.529177d0
      inquire(file='groundtd/FILE.GR',exist=lex)
      if(lex)then
       write(6,*)' Excited state gradient for ground state geometry'
       allocate(g(N),q(N),qt(N),ee(N),r(N))
c      cartesian gradient:
       open(8,file='groundtd/FILE.GR')
       do 1 ia=1,nat
1      read(8,*)(g(ix+3*(ia-1)),ix=1,3)
       close(8)
       call readsi(N,sg,eg,NQ,'groundtd/F.INP',r,ldz)
c      normal mode gradient:
       do 2 iq=1,NQ
       q(iq)=0.0d0
       do 2 ii=1,N
2      q(iq)=q(iq)+sg(ii,iq)*g(ii)*SFAC
       open(9,file='PGRAD')
       write(6,*)' mode        wg          wgp gradient (au)  Qt   n'
       write(9,*)' mode        wg          wgp gradient (au)  Qt   n'
       do 3 iq=1,NQ
       if(eg(iq).lt.y)then
        w=(q(iq)**2/2.0d0)**(1.0d0/3.0d0)
       else
        w=eg(iq)
       endif
       ee(iq)=w
       qt(iq)=q(iq)/w**2
       ani=w*qt(iq)**2/2.0d0-0.50d0
       write(9,600)iq,eg(iq)*CM,w*cm,q(iq),qt(iq),nint(ani+0.50d0)
3      write(6,600)iq,eg(iq)*CM,w*cm,q(iq),qt(iq),nint(ani+0.50d0)
600    format(i4,2f12.2,2g12.4,i6)
       close(9)
       write(6,*)'PGRAD written'

       allocate(z(nat),rg(N),x(N))
       call rdxx('ground/FILE.X',1,  nat,z,rg)
       do 4 ii=1,N
       x(ii)=rg(ii)
       do 4 iq=1,NQ
4      x(ii)=x(ii)+qt(iq)*sg(ii,iq)*SFAC*bohr
       call system('mkdir virtual')
       call wdxx('virtual/FILE.X',nat,z,x)
       call writesi(nat,N,sg,ee,NQ,'virtual/F.INP',x,z)
       write(6,*)'virtual geometry made'
      endif
      return
      end
c     ============================================================
      subroutine shiftq(nat,N,NQ,se,ep,sg,eg)
      implicit none
      integer*4 nat,N,NQ,i,j,ii,IERR,nt,it,jc,ia,jp
      real*8 se(N,N),ep(*),CM,ani,d,bohr,SFACR,TOL,t(3,3),pt,
     1sp,sg(N,N),eg(*),dc,pj,sgn,sn
      real*8,allocatable::sms(:,:),qt(:),xt(:),re(:),rg(:),x(:),
     1v(:),si(:,:),e(:,:),pc(:),st(:,:),u(:)
      integer*4,allocatable::z(:),iind(:)
      TOL=1.0d-10
      CM=219474.0d0
      bohr=0.529177d0
      SFACR=0.02342179d0
      allocate(sms(NQ,NQ),qt(NQ),xt(N),z(nat),re(N),rg(N),x(N),
     1v(NQ),si(NQ,NQ),e(NQ,2*NQ),pc(NQ),u(N),iind(NQ),st(N,NQ))

      call rdxx('ground/FILE.X',1,  nat,z,rg)
      call rdxx('exc.as.ground.x',1,nat,z,re)
      d=0.0d0
      do 5 i=1,N
      x(i)=(rg(i)-re(i))/bohr
5     d=d+x(i)**2
      write(6,607)d
607   format(' d = ',g12.4)

      do 1  i=1,NQ
      v(i)=0.0d0
      do 11 ii=1,N
11    v(i)=v(i)+x(ii)*se(ii,i)*SFACR
      do 1 j=1,NQ
      sms(i,j)=0.0d0
      do 1 ii=1,N
1     sms(i,j)=sms(i,j)+se(ii,i)*SFACR**2*se(ii,j)

      call INV(NQ,sms,si,NQ,TOL,e,IERR)
      if(IERR.ne.0)call report('inverson matrix cannot be found')

      do 2 i=1,NQ
      qt(i)=0.0d0
      do 2 j=1,NQ
2     qt(i)=qt(i)+si(i,j)*v(j) 

      d=0.0d0
      do 7 ii=1,N
      xt(ii)=re(ii)/bohr
      do 71 j=1,NQ
71    xt(ii)=xt(ii)+se(ii,j)*SFACR*qt(j)
      d=d+(rg(ii)/bohr-xt(ii))**2
7     xt(ii)=xt(ii)*bohr
      write(6,607)d
      open(44,file='test1.x')
      write(44,*)'Excited orientation to match the ground:'
      write(44,*)nat
      do 13 i=1,nat
13    write(44,608)z(i),(xt(ii+3*(i-1)),ii=1,3)
608   format(i4,3f12.6)
      close(44)
      
      nt=0
      open(9,file='QSTATE')
      write(9,600)
600   format(/,' mode       freq      Qt      n  nround:')
      do 3 i=1,NQ
      ani=ep(i)*qt(i)**2/2.0d0-0.50d0
      nt=nt+nint(ani+0.50d0)
3     write(9,601)i,ep(i)*CM,qt(i),ani,nint(ani+0.50d0)
 601  format(i4,f12.2,2G12.4,i6)
      close(9)
      write(6,*)nt,' - excited'
      write(6,*)'QSTATE written'

      write(6,610)'Excited state vibrations'
610   format(a24,'Projection of geometry difference to normal modes')
      open(44,file='TMX.TXT')
      do 9 i=1,3
9     read(44,*)(t(j,i),j=1,3)
      close(44)
c     excited state S-matrix in ground-like orientation:
      do 14 j=1,NQ
      do 14 ia=1,N/3
      do 14 i=1,3
      ii=3*(ia-1)
14    st(i+ii,j)=se(1+ii,j)*t(1,i)+se(2+ii,j)*t(2,i)+se(3+ii,j)*t(3,i)
c     as_ground_ai=tj,i exc_aj
      do 4 j=1,NQ
      iind(j)=j
c     unit vector along the mode:
      do 17 i=1,N
17    u(i)=st(i,j)
      call norm(u,N)
c     projection to the geometry change:
      pc(j)=dabs(sp(u,x,N))
c     in units of normal mode amplitude (~ 1/sqrt(w)):      
4     pc(j)=pc(j)*dsqrt(ep(j))
c     normalize:
      call norm(pc,NQ)
c     order
      do 10 j=1,NQ-1
      do 10 jp=j+1,NQ
      if(pc(jp).gt.pc(j))then
       pt=pc(jp)
       pc(jp)=pc(j)
       pc(j)=pt
       it=iind(jp)
       iind(jp)=iind(j)
       iind(j)=it
      endif
10    continue
      write(6,609)
609   format(' Mode exc   contrinution (%)  closest ground' )
      do 12 ii=1,NQ
      j=iind(ii)
c     most similar ground state mode:
      sgn=0.0d0
      sn=0.0d0
      jc=1
      dc=0.0d0
      do 15 i=1,N
      sgn=sgn+sg(i,1)**2
      sn =sn +st(i,j)**2
15    dc=dc+sg(i,1)*st(i,j)
      dc=dabs(dc)/dsqrt(sn*sgn)
      do 19 jp=2,NQ
      pj=0.0d0
      sgn=0.0d0
      sn=0.0d0
      do 16 i=1,N
      sgn=sgn+sg(i,jp)**2
      sn =sn +st(i,j)**2
16    pj=pj+sg(i,jp)*st(i,j)
      pj=dabs(pj)/dsqrt(sn*sgn)
      if(pj.gt.dc)then
       dc=pj
       jc=jp
      endif
19    continue
12    write(6,612)ii,j,ep(j)*CM,pc(ii)**2*100.0d0,jc,eg(jc)*CM
612   format(2i5,'(',f12.3,' cm-1)',f10.2,i6,'(',f12.2,' cm-1)')

      write(6,610)'Ground state vibrations '
      do 20 j=1,NQ
      iind(j)=j
      do 21 i=1,N
21    u(i)=sg(i,j)
      call norm(u,N)
      pc(j)=dabs(sp(u,x,N))
20    pc(j)=pc(j)*dsqrt(eg(j))
      call norm(pc,NQ)
      do 24 j=1,NQ-1
      do 24 jp=j+1,NQ
      if(pc(jp).gt.pc(j))then
       pt=pc(jp)
       pc(jp)=pc(j)
       pc(j)=pt
       it=iind(jp)
       iind(jp)=iind(j)
       iind(j)=it
      endif
24    continue
      write(6,619)
619   format(' Mode gr   contrinution (%)  closest exc' )
      do 25 ii=1,NQ
      j=iind(ii)
      sgn=0.0d0
      sn=0.0d0
      jc=1
      dc=0.0d0
      do 26 i=1,N
      sgn=sgn+st(i,1)**2
      sn =sn +sg(i,j)**2
26    dc=dc+st(i,1)*sg(i,j)
      dc=dabs(dc)/dsqrt(sn*sgn)
      do 27 jp=2,NQ
      pj=0.0d0
      sgn=0.0d0
      sn=0.0d0
      do 28 i=1,N
      sgn=sgn+st(i,jp)**2
      sn =sn +sg(i,j)**2
28    pj=pj+st(i,jp)*sg(i,j)
      pj=dabs(pj)/dsqrt(sn*sgn)
      if(pj.gt.dc)then
       dc=pj
       jc=jp
      endif
27    continue
25    write(6,612)ii,j,eg(j)*CM,pc(ii)**2*100.0d0,jc,ep(jc)*CM

      return
      end
c     ============================================================
      subroutine shifti(nat,N,NQ,se,ep)
c     iterative way
      implicit none
      integer*4 nat,ia,N,NQ,i,j,ii
      real*8 se(N,N),ep(*),CM,ani,sx,bohr,s2,SFACR,d
      real*8,allocatable::qt(:),rg(:),re(:),x(:),xt(:)
      integer*4,allocatable::z(:)
      CM=219474.0d0
      bohr=0.529177d0
      SFACR=0.02342179d0

      allocate(z(nat),rg(N),re(N),x(N),xt(N))
      call rdxx('ground/FILE.X',1  ,nat,z,rg)
      call rdxx('exc.as.ground.x',1,nat,z,re)
      d=0.0d0
      do 5 i=1,N
      x(i)=(rg(i)-re(i))/bohr
5     d=d+x(i)**2
      write(6,607)d
607   format(' d = ',g12.4)

c     delta(i,j)=se(ii,j)*se(ii,i)*m(ia)
c     m is in g/mol

      allocate(qt(NQ))
      call vz(qt,NQ)
      do 1 i=1,NQ

      sx=0.0d0
      do 2 ii=1,N
      sx=sx+x(ii)*se(ii,i)*SFACR
      do 2 j=1,i-1
2     sx=sx-se(ii,i)*SFACR**2*se(ii,j)*qt(j)

      s2=0.0d0
      do 6 ii=1,N
6     s2=s2+se(ii,i)**2*SFACR**2
1     qt(i)=sx/s2

      d=0.0d0
      do 7 ii=1,N
      xt(ii)=re(ii)/bohr
      do 71 j=1,NQ
71    xt(ii)=xt(ii)+se(ii,j)*SFACR*qt(j)
      d=d+(rg(ii)/bohr-xt(ii))**2
7     xt(ii)=xt(ii)*bohr
      write(6,607)d
      open(44,file='test.x')
      write(44,*)'Excited orientation to match the ground:'
      write(44,*)nat
      do 13 ia=1,nat
13    write(44,608)z(ia),(xt(ii+3*(ia-1)),ii=1,3)
608   format(i4,3f12.6)
      close(44)

      open(9,file='ISTATE')
      write(9,600)
600   format(/,' mode       freq      Qt      n  nround:')
      do 3 i=1,NQ
      ani=ep(i)*qt(i)**2/2.0d0-0.50d0
3     write(9,601)i,ep(i)*CM,qt(i),ani,nint(ani+0.50d0)
 601  format(i4,f12.2,2G12.4,i6)
      close(9)
      write(6,*)'ISTATE written'
      return
      end
c     ============================================================
      subroutine shifte(nat,N,NQ,sg,se,eg,ep,m,lmscan,fac,
     1A,B,C,D,E,F,FI,G,GP,GH,GHP,J,JT,K,NB,NP,TOL,
     2lspec,del,e00,wsmin,wsmax,nps,ncut,LEXCL,
     1isu)
      implicit none
      integer*4 nat,N,NQ,i,ii,nt,ia,ix,NB,NP,nps,ncut,LEXCL,isu
      real*8 se(N,N),ep(*),CM,bohr,m(*),fac,TOL,EH,
     1del,e00,wsmin,wsmax,sg(N,N),eg(*)
      real*8 A(NQ,NQ),B(*),C(NQ,NQ),D(*),E(NQ,NQ),F(NQ,NQ),FI(NQ,NQ),
     1G(NQ,NQ),GP(NQ,NQ),GH(NQ,NQ),GHP(NQ,NQ),J(NQ,NQ),JT(NQ,NQ),K(*)
      real*8,allocatable::qt(:),re(:),rg(:),x(:),KT(:),EK(:),ani(:)
      integer*4,allocatable::z(:),nopt(:),noptnew(:)
      logical lmscan,lspec
      CM=219474.0d0
      bohr=0.529177d0
      allocate(qt(NQ),z(nat),re(N),rg(N),x(N),EK(NQ))
      call vz(EK,NQ)

c     x = xe - xg:
      call rdxx('ground/FILE.X',1,  nat,z,rg)
      call rdxx('exc.as.ground.x',1,nat,z,re)
      do 5 ii=1,N
5     x(ii)=(re(ii)-rg(ii))/bohr

c     qt(e) = -se^(-1) x:
      do 2 i=1,NQ
      qt(i)=0.0d0
      do 2 ia=1,nat
      do 2 ix=1,3
      ii=ix+3*(ia-1)
2     qt(i)=qt(i)-se(ii,i)*x(ii)*m(ia)*dsqrt(1822.89d0)

c     Energy = 0.5 sum(i) wi^2 qti^2
      EH=0.d0
      do 4 i=1,NQ
4     EH=EH+0.50d0*qt(i)**2*ep(i)**2
      write(6,6007)EH,EH*219474.0d0
6007  format('EH = ',f12.6,'au = ',f10.2,' cm-1')

      nt=0
      allocate(nopt(NQ),noptnew(NQ),ani(NQ))
      do 3 i=1,NQ
      ani(i)=0.5d0*(ep(i)*qt(i)**2-1.0d0)
      nopt(i)=max(0,nint(ani(i)+0.50d0))
      noptnew(i)=nopt(i)
3     nt=nt+nopt(i)
      write(6,*)nt,' - excited'

      if(ncut.gt.0)then
       allocate(KT(NQ))
       do 12 i=1,NQ
12     KT(i)=K(i)
       do 1 i=1,NQ
       if(nopt(i).gt.ncut)then
c       redefinition of excited normal mode Qk -> Qk
        write(6,6008)i,nopt(i)
6008    format(i4,':',i4,' -> 0')
        EK(i)=dble(nopt(i))*ep(i)
        noptnew(i)=0
        do 11 ii=1,NQ
11      KT(ii)=KT(ii)+J(ii,i)*qt(i)
       endif
1      continue
       call dofac(NQ,fac,J,JT,F,FI,G,GP,KT)
       write(6,300)fac,fac**2*100.0d0
300    format('<0|0*p> = ',g12.4,' (',f6.2,'%)')
       call newbd(NQ,G,KT,J,JT,FI,GH,GHP,B,D)
       call wdusch(NQ,N,A,B,C,D,E,JT,KT,eg,ep,sg,se,rg,re,m)
      endif

      open(9,file='PSTATE')
      write(9,600)
600   format(/,' mode       freq      Qt      ne  nround nnew  EK:')
      do 31 i=1,NQ
31    write(9,601)i,ep(i)*CM,qt(i),ani(i),nopt(i),noptnew(i),EK(i)
 601  format(i4,f12.2,2G12.4,2i6,g14.6)
      write(9,601)NB
      close(9)
      write(6,*)'PSTATE written'

      if(lmscan)call mscan(noptnew,NQ,fac,C,D,NB,NP,TOL,
     1lspec,del,e00,wsmin,wsmax,nps,ep,EK,LEXCL,isu)

      return
      end
c     ============================================================
      subroutine mscan(so,NQ,fac,C,D,NB,NP,TOL,
     1lspec,del,e00,wsmin,wsmax,nps,ep,EK,LEXCL,isu)
      implicit none
      real*8 fac,FC,TOL,r,ojmax,del,e00,wsmin,wsmax,ep(*),
     1C(NQ,NQ),D(NQ),of,EK(*)
      integer*4 Ni,Nj,so(*),NQ,i,np0,nstart,ii,j,les,
     1kk,nnmax,jj,nend,NB,NP,nps,Nm,LEXCL,isu
      logical lspec
      integer*4,allocatable::si(:),sj(:),st(:),nbi(:),nn(:,:)
      real*8,allocatable::ob(:),oj(:)

      np0=1000000
      write(6,*)'scan called'

      write(6,*)'mother state:'
      do 3 kk=1,NQ
3     if(so(kk).ne.0)write(6,607)kk,so(kk)
607   format(i4,':',i2,$)
      write(6,*)
      Nm=les(NQ,so)
      write(6,*)Nm,' excited'
      allocate(sj(Nm+1))
      call viz(sj,Nm+1)
      call puts(NQ,so,sj)
      of=FC(fac,sj,Nm,np0,C,D,NQ)
      deallocate(sj)
      write(6,6093)of,tol
6093  format(' < 0 | mo> = ',e12.4,' tol: ',g12.4)
      write(6,*)


      Ni=0
      allocate(si(Ni+1),st(NQ),nbi(NQ),nn(NQ,NB),ob(NB),oj(2*NP+1))
      call viz(si,Ni+1)

      open(9,file='PSTATE.TAB')
      open(91,file='PROFILES.TXT')
      do 1 i=1,NQ
c     reference state to st:
      nstart=max(0,so(i)-NP)
      do 11 ii=-NP,NP
c     <0|j>:
      oj(ii+NP+1)=0.0d0
      if(ii+so(i).ge.nstart)then
       do 10 j=1,NQ
10     st(j)=so(j)
       st(i)=so(i)+ii
c      put to short record:
       Nj=les(NQ,st)
       allocate(sj(Nj+1))
       call viz(sj,Nj+1)
       call puts(NQ,st,sj)
       oj(ii+NP+1)=FC(fac,sj,Nj,np0,C,D,NQ)
       deallocate(sj)
      endif
11    continue
c
c     select maximum N:
      nnmax=1
      ojmax=0.0d0
      do 12 ii=-NP,NP
      if(dabs(oj(ii+NP+1)).gt.ojmax)then
       ojmax=dabs(oj(ii+NP+1))
       nnmax=so(i)+ii
      endif
12    continue

c     buffer with quantum number giving the best overlaps:
      nbi(i)=0
      do 13 ii=-NP,NP
      if(ojmax.gt.0.0d0)then
       r=dabs(oj(ii+NP+1))/ojmax
       if(r.gt.tol)then
        if(nbi(i).eq.0)then
         nbi(i)=nbi(i)+1
         nn(i,nbi(i))=so(i)+ii
         ob(nbi(i))=r
        else
         do 14 jj=1,nbi(i)
14       if(r.gt.ob(jj))goto 141
         if(nbi(i).lt.NB)then
          nbi(i)=nbi(i)+1
          nn(i,nbi(i))=so(i)+ii
          ob(nbi(i))=r
         endif
         goto 151
141      nend=min(NB,nbi(i)+1)
         do 15 kk=nend,jj+1,-1
         ob(kk)=ob(kk-1)
15       nn(i,kk)=nn(i,kk-1)
         ob(jj)=r
         nn(i,jj)=so(i)+ii
         nbi(i)=nend   
151      continue
        endif
       endif
      endif
13    continue

      write(9 ,601)i,nnmax,nbi(i),(nn(i,kk),kk=1,nbi(i))
601   format(100i4)
1     write(91,691)(oj(ii+NP+1),ii=-NP,NP)
691   format(100E10.2)
      close(9)
      write(6,*)'PSTATE.TAB written'

      
      if(isu.eq.1)call statesum1(fac,NQ,C,D,
     1lspec,del,e00,wsmin,wsmax,nps,ep,EK,LEXCL)

      if(isu.eq.2)call statesum2(fac,NQ,np0,C,D,NB,nbi,nn,
     1lspec,del,e00,wsmin,wsmax,nps,ep,EK)
      return
      end

c     ==============================================================
      subroutine viz(v,n)
      integer*4 v(*),i,n
      do 1  i=1,n
1     v(i)=0
      return 
      end
c     ==============================================================
      subroutine  puts(NQ,sd,sj)
c     rewrite vib state from long to compact notation:
      integer*4 NQ,sd(*),sj(*),kk,i,ii
      kk=0
      do 21 i=1,NQ
      do 21 ii=1,sd(i)
      kk=kk+1
21    sj(kk)=i
      return
      end
c     ==============================================================
      function les(N,s)
      integer*4 les,N,s(*),u,i
      u=0
      do 1 i=1,N
1     u=u+s(i)
      les=u
      return
      end
c     ==============================================================
      function FC(fac,si,Nexc,np0,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,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,Nexc),it(np0))
      p(np)=1.0d0
      NN(np)=Nexc
      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)
      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,Nexc,p,NN,it,NN(ii)-1,
     1 pini*D(iq)/dsqrt(dble(2*nu)))

       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,Nexc,p,NN,it,NN(ii)-2,
     1  pini*dsqrt(dble(nu-1)/dble(nu))*C(iq,iq))
       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,Nexc,p,NN,it,NN(ii)-2,
     1  pini*dsqrt(dble(nuj)/dble(nu))*C(iq,jq))
       endif
106    jold=jq

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

       goto 777
       
      endif
101   continue

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

      FC=p(1)*fac
      return
      end
c     ==============================================================
      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
c     ==============================================================
      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
c     ==============================================================
      subroutine statesum2(fac,NQ,np0,C,D,NB,nbi,nn,
     1lspec,del,e00,wsmin,wsmax,nps,ep,EK)
      implicit none
      integer*4 NQ,ii,NB,nbi(*),LEXP,nn(NQ,NB),les,np0,nps
      real*8 fac,proc,FC,procold,del,e00,wsmin,wsmax,CM,dx,EJ,
     1C(NQ,NQ),D(NQ),y,ep(*),ECM,emin,emax,EK(*),estart,p8
      logical lspec
      integer*4,allocatable::sj(:),st(:),is(:)
      integer*8 nt,i8
      real*8,allocatable::s(:)
      allocate(st(NQ),is(NQ))
      call viz(is,NQ)

      estart=e00
      do 100 ii=1,NQ
100   estart=estart+EK(ii)

      proc=0.0d0
      CM=219474.630d0
      procold=0.0d0

      if(lspec)then
       emin=1.0d99
       emax=0.0d0
       allocate(s(nps))
       call vz(s,nps)
       dx=(wsmax-wsmin)/dble(nps-1)
      endif

      nt=1
      do 1 ii=1,NQ
1     nt=int(nt,8)*int(nbi(ii),8)
      write(6,*)nt,' terms'

      do 2 i8=1,nt
      if(proc.gt.procold+0.1d0)then
       p8=100.0d0*dble(i8)/dble(nt)
       write(6,6000)i8,p8,proc
6000   format(i20,f9.4,'%',f9.4,'%')
       procold=proc
      endif
      call getis(i8,is,st,NQ,nn,NB,nbi)
      
c     rewrite final state st into short form in sj:
      LEXP=les(NQ,st)
      allocate(sj(LEXP+1))
      call puts(NQ,st,sj)
c

c     <0|*>
      y=FC(fac,sj,LEXP,np0,C,D,NQ)**2
      proc=proc+100.0d0*y

      if(lspec)then
c      energy:
       EJ=estart
       DO 1002 ii=1,LEXP
1002   EJ=EJ+ep(sj(ii))
       ECM=EJ*CM
       if(ECM.lt.emin)emin=ECM
       if(ECM.gt.emax)emax=ECM
       call spread(nps,y,dx,wsmin,ECM,del,s)
      endif

2     deallocate(sj)

      if(lspec)then
       call ws(s,wsmin,dx,nps)
       write(6,604)emin,emax
604    format(' emin: ',f12.2,' emax: ',f12.2,' cm-1')
      endif


      return
      end
c     ==============================================================
      subroutine spread(nps,y,dx,wmin,E,del,s)
      implicit none
      real*8 y,dx,wmin,E,del,s(*),x
      integer*4 nps,i,nc,n3,nstart,nend
c     closest point, approx:
      nc=nint((E-wmin)/dx)
c     points over three delta:
      n3=nint(3.0d0*del/dx)
      nstart=max(nc-n3,1)
      nend  =min(nc+n3,nps)
      x=wmin+dx*dble(nstart-2)
      do 1 i=nstart,nend
      x=x+dx
1     s(i)=s(i)+y*exp(-((E-x)/del)**2)
      return
      end 
c     ==============================================================
      subroutine ws(s,wmin,dx,nps)
      implicit none
      real*8 s(*),wmin,dx,x,y
      integer*4 nps,i
      open(9,file='S.PRN')
      x=wmin-dx
      do 1 i=1,nps
      x=x+dx
      y=s(i)
      if(dabs(y).lt.1.0d-10)y=0.0d0
 1    write(9,90)x,y
 90   format(f10.2,e12.4)
      close(9)
      return
      end
c     ==============================================================
      subroutine dofac(NQ,fac,J,JT,F,FI,G,GP,K)
      integer*4 NQ
      real*8 fac,G(NQ,NQ),GP(NQ,NQ),K(*),lT,lF,LOGDET,lQ,sp1,sp2,
     1F(NQ,NQ),J(NQ,NQ),JT(NQ,NQ),FI(NQ,NQ),lJ,f1,sp
      real*8,allocatable:: T(:,:),TV(:),TU(:)
      allocate(T(NQ,NQ),TV(NQ),TU(NQ))
      call mm(T,G,GP,NQ)
      lT=LOGDET(T,NQ,NQ)
      lF=LOGDET(F,NQ,NQ)
      lJ=LOGDET(J,NQ,NQ)
      lQ=dble(NQ)*log(2.0d0)
      f1=exp((lQ+lt/2.0d0+lJ-lF)/2.0d0)
      call mv(TV,G,K,NQ)
      sp1=sp(K,TV,NQ)
      call mv(TU,G,K,NQ)
      call mv(TV,JT,TU,NQ)
      call mv(TU,FI,TV,NQ)
      call mv(TV,J,TU,NQ)
      call mv(TU,G,TV,NQ)
      sp2=sp(K,TU,NQ)
      fac=f1*exp(-0.50d0*(sp1-sp2))
      return
      end
c     ==============================================================
      subroutine newbd(NQ,G,K,J,JT,FI,GH,GHP,B,D)
      implicit none
      integer*4 NQ,ix
      real*8 G(NQ,NQ),K(*),JT(NQ,NQ),FI(NQ,NQ),GHP(NQ,NQ),D(*),
     1B(*),GH(NQ,NQ),J(NQ,NQ)
      real*8,allocatable::TU(:),TV(:),T(:,:),TP(:,:)
      allocate(TV(NQ),TU(NQ),T(NQ,NQ),TP(NQ,NQ))
      call mv(TV,G,K,NQ)
      call mv(TU,JT,TV,NQ)
      call mv(TV,FI,TU,NQ)
      call mv(TU,GHP,TV,NQ)
      do 17 ix=1,NQ
17    D(ix)=-2.0d0*TU(ix)
c     D= - 2 GP^1/2 Fi Jt G K

      call mm(T,JT,G,NQ)
      call mm(TP,FI,T,NQ)
      call mm(T,J,TP,NQ)
      do 13 ix=1,NQ
13    T(ix,ix)=T(ix,ix)-1.0d0
      call mv(TV,T,K,NQ)
      call mv(TU,GH,TV,NQ)
      do 14 ix=1,NQ
14    B(ix)=-2.0d0*TU(ix)
c     B= - 2 G^1/2 ( J Fi Jt G - I) K
      return
      end
c     ==============================================================
      subroutine statesum1(fac,NQ,C,D,
     1lspec,del,e00,wsmin,wsmax,nps,ep,EK,LEXCL)
      implicit none
      integer*4 NQ,ii,NB,nm1,nm2,i,ic,iex,iq,jq,Nexc,
     1nps,LEXCL,nne,nu,nuj
      real*8 fac,proc,procold,del,e00,wsmin,wsmax,CM,dx,EJ,
     1C(NQ,NQ),D(NQ),y,ep(*),ECM,emin,emax,EK(*),estart,fci
      logical lspec
      real*8,allocatable::s(:),fcs(:),fcm1(:),fcm2(:)
      integer*4,allocatable::si(:),sls(:,:),slm1(:,:),slm2(:,:),sl(:)

      estart=e00
      do 100 ii=1,NQ
100   estart=estart+EK(ii)

      proc=0.0d0
      CM=219474.630d0
      procold=0.0d0

      call statereport(NQ,LEXCL)

      if(lspec)then
       emin=1.0d99
       emax=0.0d0
       allocate(s(nps))
       call vz(s,nps)
       dx=(wsmax-wsmin)/dble(nps-1)
      endif

c     LEXCL: maximal number of excitations
      nb=NQ**LEXCL
      allocate(si(LEXCL+1),fcm1(nb),fcm2(nb),slm1(nb,NQ),
     1slm2(nb,NQ),fcs(nb),sls(nb,NQ),sl(NQ))

      si(1)=0
c     lc:current number of excitations
c     <0|0*>:
      y=fac**2
      proc=100.0d0*y

      if(lspec)then
c      energy:
       EJ=estart
       ECM=EJ*CM
       if(ECM.lt.emin)emin=ECM
       if(ECM.gt.emax)emax=ECM
       call spread(nps,y,dx,wsmin,ECM,del,s)
      endif
      nm2=0
      nm1=0
      nne=1
      do 77 ii=1,NQ
77    sls(1,ii)=0
      fcs(1)=fac
      fcm1(1)=0
      slm1(1,1)=0

      do 30000 Nexc=1,LEXCL
      write(6,6000)NExc,proc
6000  format(i2,f9.4,'%')
c     rewrite strings:
      do 71 ii=1,nm1
      fcm2(ii)=fcm1(ii)
      do 71 iq=1,NQ
71    slm2(ii,iq)=slm1(ii,iq)
      nm2=nm1
      do 81 ii=1,nne
      fcm1(ii)=fcs(ii)
      do 81 iq=1,NQ
81    slm1(ii,iq)=sls(ii,iq)
      nm1=nne
c     number of states Nexc excited
      nne=0
c     distribute Nexc excitations upon centers 1..NQ:
c     initial indices - all to the first center:
      do 41 iex=1,Nexc
41    si(iex)=1
50000 nne=nne+1
c     transcript to long:
      call trl(si,Nexc,sl,NQ)
      fci=0.0d0
c     reduce first excitation on center iq
      iq=si(1)
      nu=sl(iq)
c     from <0|i-1>, find it in previous list:
      do 44 ii=1,nm1
      do 441 jq=1,NQ
      if(jq.ne.iq)then
       if(sl(jq).ne.slm1(ii,jq))goto 44
      else
       if(nu.ne.slm1(ii,jq)+1)goto 44
      endif
441   continue
      fci=fci+D(iq)*fcm1(ii)/dsqrt(dble(2*nu))
      goto 442
44    continue
442   continue
c     from <0|i-2>, find it in pre-previous list:
      if(nu.gt.1)then
       do 54 ii=1,nm2
       do 541 jq=1,NQ
       if(jq.ne.iq)then
        if(sl(jq).ne.slm2(ii,jq))goto 54
       else
        if(nu.ne.slm2(ii,jq)+2)goto 54
       endif
541    continue
       fci=fci+C(iq,iq)*fcm2(ii)*dsqrt(dble(nu-1)/dble(nu))
       goto 542
54     continue
542    continue
      endif
c     from <0|i-1 j-1>
      do 64 ii=1,nm2
      nuj=0
      do 641 jq=1,NQ
      if(jq.ne.iq)then
       if(sl(jq).ne.slm2(ii,jq))then
        if(nuj.eq.0.and.sl(jq).eq.slm2(ii,jq)+1)then
         nuj=sl(jq)
         goto 641
        else
         goto 64
        endif
       else
        goto 64
       endif
      else
       if(nu.ne.slm2(ii,jq)+1)goto 64
      endif
641   continue
      fci=fci+C(iq,jq)*fcm2(ii)*dsqrt(dble(nuj)/dble(nu))
      goto 642
64    continue
642   continue
c     save for later:
      do 42 ii=1,NQ
42    sls(nne,ii)=sl(ii)
      fcs(nne)=fci
c
c     <0|*>
      y=fci**2
      proc=proc+100.0d0*y
      if(proc.gt.procold+1.0d0)then
       write(6,6000)Nexc,proc
       procold=proc
      endif
      if(lspec)then
c      energy:
       EJ=estart
       DO 1002 ii=1,Nexc
1002   EJ=EJ+ep(si(ii))
       ECM=EJ*CM
       if(ECM.lt.emin)emin=ECM
       if(ECM.gt.emax)emax=ECM
       call spread(nps,y,dx,wsmin,ECM,del,s)
      endif
c
c     find index to be changed
      do 80000 ic=Nexc,1,-1
80000 if(si(ic).lt.NQ)goto 90000
      goto 30000
90000 do 10000 i=ic+1,Nexc
10000 si(i)=si(ic)+1
      si(ic)=si(ic)+1
      goto 50000
30000 continue

      if(lspec)then
       call ws(s,wsmin,dx,nps)
       write(6,604)emin,emax
604    format(' emin: ',f12.2,' emax: ',f12.2,' cm-1')
      endif

      return
      end
c     ==============================================================
      subroutine trl(si,Nexc,sl,NQ)
      integer*4 si(*),Nexc,sl(*),NQ,i
      do 1 i=1,NQ
1     sl(i)=0
      do 2 i=1,Nexc
2     sl(si(i))=sl(si(i))+1
      return
      end
c     ==============================================================
      subroutine statereport(N,LEXCL)
      implicit none
      integer*4 N,LEXCL
      integer*8 i,sum
      integer*8 ,allocatable::NS(:)
      allocate(NS(LEXCL+1))
      NS(1)=1
      do 1 i=1,LEXCL
1     NS(i+1)=(NS(i)*(N+i-1))/i
      WRITE(6,4012)0,NS(1),NS(1)
4012  FORMAT(' Number of states excited',i3,'-times: ',2i20)
      sum=1
      do 2 i=1,LEXCL
      sum=sum+NS(i+1)
2     WRITE(6,3012)i,NS(i+1),sum
3012  FORMAT('                         ',i3,'-times: ',2i20)
      WRITE(6,5012)sum
5012  FORMAT('                           -------- ',/,
     1       '                             total: ',i20)
      return
      end
c     ==============================================================
      subroutine wdusch(NQ,N,A,B,C,D,E,JT,K,eg,ep,sg,se,rg,re,m)
      implicit none
      integer*4 NQ,N,iq,ix,ia
      real*8 A(NQ,NQ),B(*),C(NQ,NQ),D(*),E(NQ,NQ),JT(NQ,NQ),m(*),
     1eg(*),ep(*),K(*),sg(N,N),se(N,N),rg(*),re(*),SFAC,bohr
      real*8,allocatable::ci(:),cj(:)
      SFAC=0.0234280d0
      bohr=0.529177d0
      allocate(ci(NQ),cj(NQ))
      open(10,FILE='DUSCH.OUT')
      write(10,*)NQ
      call wdm(10,NQ,JT,'Duschinsky Matrix')
      call wdv(10,NQ,K,'Shift Vector')
      call wdm(10,NQ,A,'A Matrix')
      call wdv(10,NQ,B,'B Vector')
      call wdm(10,NQ,C,'C Matrix')
      call wdv(10,NQ,D,'D Vector')
      call wdm(10,NQ,E,'E Matrix')
      call wwm(10,NQ,eg,'Ground w')
      call wwm(10,NQ,ep,'Excited w')
      call wds(10,N,NQ,sg,'Ground S-Matrix')
      call wds(10,N,NQ,se,'Excited S-Matrix')
      do 20 iq=1,NQ
      ci(iq)=0.0d0
      cj(iq)=0.0d0
      do 20 ix=1,N
      ia=(ix+2)/3
      cj(iq)=cj(iq)+(rg(ix)-re(ix))/bohr
     1*se(ix,iq)*SFAC*m(ia)*1822.89d0
20    ci(iq)=ci(iq)+(rg(ix)-re(ix))/bohr
     1*sg(ix,iq)*SFAC*m(ia)*1822.89d0
      write(10,*)
      write(10,*)' mode  SG*dX   SE*dX  ng    ne'
      do 21 iq=1,NQ
21    write(10,621)iq,ci(iq),cj(iq),ci(iq)**2*eg(iq)/2.0d0,
     1cj(iq)**2*ep(iq)/2.0d0
621   format(i5,2f12.6,2f8.1)
      close(10)
      write(6,*)' DUSCH.OUT written'
      return
      end
c     ==============================================================
      subroutine getis(i8,is,st,NQ,nn,NB,nbi)
      implicit none
      integer*4 is(*),st(*),NQ,NB,nn(NQ,NB),nbi(*),ii,k,kk
      integer*8 i8,isum,ifa,ifa2
      do 21 ii=NQ,1,-1
      isum=i8-int(1,8)
      do 221 k=1,NQ-ii
      ifa=int(nbi(1),8)
      do 223 kk=2,NQ-k
223   ifa=ifa*int(nbi(kk),8)
221   isum=isum-ifa*int(is(NQ-k+1)-1,8)
      ifa2=1
      do 222 k=1,ii-1
222   ifa2=ifa2*int(nbi(k),8)
21    is(ii)=int(isum/ifa2+int(1,8),4)
      do 224 k=1,NQ
224   st(k)=nn(k,is(k))
      return
      end
c     ==============================================================
      subroutine seproj(N,NF,NQ,se,sgf,m)
      implicit none
      integer*4 N,NF,NQ,i,j,ix
      real*8 se(N,N),sgf(N,N),norm,m(*),spm,sij
c     each mode in se:
      do 1 i=1,NQ
c     each deleted mode:
      do 2 j=1,NF
      sij=spm(sgf,se,j,i,N,m)
      do 2 ix=1,N
2     se(ix,i)=se(ix,i)-sij*sgf(ix,j)
c     previous modes
      do 3 j=1,i-1
      sij=spm(se,se,j,i,N,m)
      do 3 ix=1,N
3     se(ix,i)=se(ix,i)-sij*se(ix,j)
c     renormalize:
      norm=1.0d0/dsqrt(spm(se,se,i,i,N,m))
      do 1 ix=1,N
1     se(ix,i)=se(ix,i)*norm
      write(6,*)NF,' modes projected from se'
      return
      end
c     ==============================================================
      function spm(s1,s2,i,j,N,m)
      implicit none
      integer*4 i,j,N,ix,ia
      real*8 s1(N,N),s2(N,N),a,m(*),spm
      a=0.0d0
      do 1 ia=1,N/3
      do 1 ix=1,3
1     a=a+s1(ix+3*(ia-1),i)*s2(ix+3*(ia-1),j)*m(ia)
      spm=a
      return
      end
c     ==============================================================
      subroutine zzz(NQ,JT,K,G,GP,GH,GHP,J,T,TP,F,FI,E2,A,B,TV,TU,C,D,
     1E,K0)
      implicit none
      integer*4 NQ
      real*8 JT(NQ,NQ),K(NQ),G(NQ,NQ),GP(NQ,NQ),GH(NQ,NQ),GHP(NQ,NQ),
     1J(NQ,NQ),T(NQ,NQ),TP(NQ,NQ),F(NQ,NQ),FI(NQ,NQ),E2(NQ,2*NQ),
     1A(NQ,NQ),B(NQ),TV(NQ),TU(NQ),C(NQ,NQ),D(NQ),E(NQ,NQ),K0(NQ)
      A=0.0d0
      B=0.0d0
      C=0.0d0
      D=0.0d0
      E=0.0d0
      E2=0.0d0
      F=0.0d0
      FI=0.0d0
      G=0.0d0
      GH=0.0d0
      GHP=0.0d0
      GP=0.0d0
      J=0.0d0
      JT=0.0d0
      K=0.0d0
      K0=0.0d0
      T=0.0d0
      TP=0.0d0
      TV=0.0d0
      TU=0.0d0
      return
      end

      subroutine rpl(N,NQ,nat,m,sg,se)
      implicit none
      integer*4  N,NQ,nat,ie,ig,ix,ii,ia
      real*8 m(nat),sg(N,N),se(N,N),cl,sp
      integer*4 ,allocatable::ige(:),it(:)
      allocate(ige(NQ),it(NQ))
      write(6,600)
600   format(' excited mode / closest ground product')
      it=0
      ie=1
111   ige(ie)=0    
      cl=0.0d0
      do 2 ig=1,NQ
      if(it(ig).eq.0)then
       sp=0.0d0
       do 21 ia=1,nat
       do 21 ix=1,3
       ii=ix+3*(ia-1)
21     sp=sp+sg(ii,ig)*m(ia)*se(ii,ie)
       if(dabs(sp).gt.cl)then
        ige(ie)=ig
        cl=dabs(sp)
       endif
      endif
2     continue
      if(ige(ie).eq.0)call report('mode cannot be assigned')
      it(ige(ie))=ie
      write(6,601)ie,ige(ie),cl
601   format(i5,' /',i5,f6.3,$)
      if(mod(ie,3).eq.0)write(6,*)
      do 1 ia=1,nat
      do 1 ix=1,3
      ii=ix+3*(ia-1)
1     se(ii,ie)=sg(ii,ige(ie))
      if(ie.lt.NQ)then
       ie=ie+1
       goto 111
      endif
      write(6,602)NQ
602   format(i6,' S-vectors replaced')
      return
      end