program dusch
c     Duschinsky transformation
      implicit none
      real*8 LOGDET
      real*8 bohr,u(3,3),proc,fac,sp,sp1,sp2,TOL,
     1lT,lJ,lF,lQ,overlap,wmin
      integer*4 ia,ix,N,nat,NE,NQ,IERR,maxit
      real*8,allocatable::rg(:),re(:),se(:,:),sg(:,:),JT(:,:),
     1eg(:),ep(:),k(:),m(:),F(:,:),G(:,:),GP(:,:),T(:,:),TP(:,:),
     1FI(:,:),E2(:,:),GH(:,:),GHP(:,:),A(:,:),TV(:),TU(:),
     1B(:),J(:,:),C(:,:),D(:),E(:,:)
      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

      bohr=0.529177d0
      TOL=1.0d-10

      call readopt(lit,maxit,lwr,wmin)
c
c     Read Geometry of the ground and excited states:
      call rdxx('ground/FILE.X',0,nat,z,rg)
      N=3*nat
      allocate(rg(N),z(nat),re(N),sg(N,N),se(N,N),eg(N),ep(N),m(nat))
      call rdxx('ground/FILE.X',1,nat,z,rg)
      call rdxx('excited/FILE.X',1,nat,z,re)
      do 6 ia=1,nat
6     m(ia)=dble(amas(z(ia)))
      
c     Find the transformation matrix, transform and shift re to rg:
      call xst(nat,rg,re,m,z,u)

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')
      write(6,*)'reading excited state S-matrix'
      call readsi(N,se,ep,NE,'excited/F.INP')
      if(NQ.ne.NE)call report('NE<>NQ')
      call fixmodes(N,NQ,eg,sg,ep,se,amas,z)

c     transform se to the ground system:
      call trse(N,se,NQ,u)

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))
      do ix=1,NQ
      do ia=1,2*NQ
      E2(ix,ia)=0.0d0
      enddo
      enddo

      call mz(A,NQ)
      call vz(B,NQ)
      call mz(C,NQ)
      call vz(D,NQ)
      call mz(E,NQ)
      call mz(F,NQ)
      call mz(FI,NQ)
      call mz(G,NQ)
      call mz(GH,NQ)
      call mz(GHP,NQ)
      call mz(GP,NQ)
      call mz(J,NQ)
      call mz(JT,NQ)
      call vz(K,NQ)
      call mz(T,NQ)
      call mz(TP,NQ)
      call vz(TV,NQ)
      call vz(TU,NQ)
      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)
      proc=fac*fac*100.0d0
      write(6,300)fac,proc
300   format('<0|0*> = ',g12.4,' (',f6.2,'%)')
      
      write(6,*)sp1,sp2
      write(6,*)'NQ',NQ
      write(6,*)'A',A(1,1),A(NQ,NQ)
      write(6,*)'B',B(1),B(NQ)
      write(6,*)'C',C(1,1),C(NQ,NQ)
      write(6,*)'D',D(1),D(NQ)
      write(6,*)'E',E(1,1),E(NQ,NQ)
      write(6,*)'F',F(1,1),F(NQ,NQ)
      write(6,*)'G',G(1,1),G(NQ,NQ)
      write(6,*)'GP',GP(1,1),GP(NQ,NQ)
      
      write(6,*)' Vertical approximation:'
      do 19 ix=1,NQ
19    K(ix)=0.0d0

      proc=0.0d0
      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)
      fac=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=fac*exp(-0.50d0*(sp1-sp2))
      proc=fac*fac*100.0d0
      write(6,300)fac,proc

      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)


     
      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     ==============================================================
      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 readopt(lit,maxit,lwr,wmin)
      IMPLICIT none
      logical lit,lwr,lex
      integer*4 maxit
      real*8 wmin
      character*3 key
      lit=.false.
      wmin=1.0d0
      maxit=100
      lwr=.true.
      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.'ITE')read(11,*)lit
       if(key.eq.'LWR')read(11,*)lwr
       if(key.eq.'MAX')read(11,*)maxit
       if(key.eq.'WMI')read(11,*)wmin
       goto 1
11     close(11)
      endif
      return
      end
c     ==============================================================
      subroutine fixmodes(N,NQ,g,sg,p,se,amas,z)
      IMPLICIT none
      integer*4 N,NQ,I,iz,ifx,j,ix,no,jmin,nblo,z(*)
      real*8 g(*),p(*),sg(N,N),se(N,N),CM,wlim,dij,dmin,t
      real amas(*)
      integer*4,allocatable::bl(:)
      character*3 key
      logical lex,lcor,lonly
      lcor=.false.
      lonly=.false.
      nblo=0
      wlim=600.0d0
      CM=219474.630d0
      no=NQ
      ifx=1
      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.'FIX')read(11,*)ifx
       if(key.eq.'WLI')read(11,*)wlim
       if(key.eq.'COR')read(11,*)lcor
       if(key.eq.'ONL')read(11,*)lonly
       if(key.eq.'BLO')then
        read(11,*)nblo
        allocate(bl(nblo))
        read(11,*)bl
       endif
       goto 1
11     close(11)
      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
      
      if(ifx.eq.1)then
       iz=0
       do 2 I=1,NQ
       if(g(I).lt.0.0d0)then
        g(I)=100.0d0/CM
        iz=iz+1
       endif
       if(p(I).lt.0.0d0)then
        p(I)=100.0d0/CM
        iz=iz+1
       endif
 2     continue
       if(iz.gt.0)write(6,*)iz,' neg. modes made 100 cm-1 modes'
      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
        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).lt.wlim.or.p(i).lt.wlim)then
        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

      if(nblo.ne.0)then
       do 19 I=1,NQ
       do 19 j=1,nblo
19     if(I.eq.bl(j))p(I)=-3.333d0
       iz=0
68     do 17 I=1,NQ
       if(p(i).eq.-3.333d0)then
        do 18 j=i,NQ-1
        g(j)=g(j+1)
        p(j)=p(j+1)
        do 18 ix=1,N
        sg(ix,j)=sg(ix,j+1)
18      se(ix,j)=se(ix,j+1)
        iz=iz+1
        NQ=NQ-1
        goto 68
       endif
 17    continue
       if(iz.gt.0)write(6,*)iz,' modes blocked (deleted)'
      endif

      if(NQ.ne.no)write(6,*)NQ,' modes now'
      return
      end

      SUBROUTINE readsi(N3,S,E,NQ,fn)
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
      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,*)
      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'

      iz=0
c     delete zero modes if exist:
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'

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

      DO 3 I=1,NQ
3     E(I)=E(I)/CM
      
      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,/))
      call trx(nat,rn,re,t)
      open(44,file='exc.as.ground.x')
      write(44,*)'Excited orientation to match the ground:'
      write(44,*)nat
      do 13 ia=1,nat
      do 131 i=1,3
131   re(i+3*(ia-1))=cmg(i)+rn(i+3*(ia-1))
13    write(44,602)z(ia),(re(i+3*(ia-1)),i=1,3)
602   format(i4,3f12.6)
      close(44)
      write(6,*)'exc.as.ground.x written'
      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)
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,iq
      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,
     4sp1,sp2,fac,lt,lF,lJ,lQ,LOGDET,sp,SFAC
      real*8,allocatable::ci(:),cj(:)
      call makeG(NQ,G,GP,GH,GHP,eg,ep)

      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)

      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)


      call mt(J,JT,NQ)
      write(6,*)'J',J(1,1),J(NQ,NQ)
      write(6,*)'JT',JT(1,1),JT(NQ,NQ)
      call mm(T,G,J,NQ)
      write(6,*)'T',T(1,1),T(NQ,NQ)
      call mm(TP,JT,T,NQ)
      write(6,*)'TP',TP(1,1),TP(NQ,NQ)
      call ms(F,TP,GP,NQ)
      write(6,*)'F',F(1,1),F(NQ,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,GH,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

      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)
      fac=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)
      overlap=fac*exp(-0.50d0*(sp1-sp2))

      if(ic.eq.1)then
       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')
       allocate(ci(NQ),cj(NQ))
       SFAC=0.0234280d0
       do 20 iq=1,NQ
       ci(iq)=0.0d0
       cj(iq)=0.0d0
       do 20 ix=1,N
       ia=(ix-1)/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,
     1 cj(iq)**2*ep(iq)/2.0d0
621    format(i5,2f12.6,2f8.1)
       close(10)
       write(6,*)' DUSCH.OUT written'
      endif

      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)
      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
      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)
      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)
      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)
      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