program ch3n
c     two ch3 group rotation, molecular spectra
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      parameter (N0=20,nf0=1500,nw0=10,nh0=1200,nd0=2)
      DIMENSION wc(nd0,nw0),ws(nd0,nw0),
     2APTu(nw0**nd0,N0*3,3),AATu(nw0**nd0,N0*3,3),
     3ALu(nw0**nd0,N0*3,3,3),Gu(nw0**nd0,N0*3,3,3),
     1Au(nw0**nd0,N0*3,3,3,3),wu(nw0**nd0,3*N0),
     4W(3*N0),hr(nh0,nh0),qr(nh0),sr(nh0,nh0),
     5EM(nh0),imr(nh0),mr(nh0),anr(nh0),hri(3*N0,nh0,nh0),
     5inr(nh0),nrs(nh0),
     6sri(3*N0,nh0,nh0),qri(3*N0,nh0),ht(nh0,nh0),st(nh0,nh0),
     7qt(nh0),fis(nd0),fie(nd0),dfi(nd0),ncos(nd0),nsin(nd0)
      character*9 s9
      integer*4 state(3*N0)
      logical ldiag,ldia2,lbath,lscale
      parameter (nt0=100000)
      integer*4 s1(3*N0,nt0),s2(3*N0,nt0),i44(8,nt0)
      dimension et(nt0),rt(nt0)
      common/trb/s1,s2,et,rt,itb,i44,lbath

      Z0=0.0d0
      pi=4.0d0*atan(1.0d0)

      open(9,file='SSDN.OPT')
      read(9,*)ndf
      do 100 ir=1,ndf
      read(9,*)ncos(ir)
      if(ncos(ir).gt.nw0)call report('too many cosinus freqs')
      read(9,*)(wc(ir,i),i=1,ncos(ir))
      read(9,*)nsin(ir)
      if(nsin(ir)+ncos(ir).gt.nw0)call report('too many c and s freqs')
      read(9,*)(ws(ir,i),i=1,nsin(ir))
      read(9,*)fis(ir),fie(ir),dfi(ir)
      fis(ir)=fis(ir)*pi/180.0d0
      fie(ir)=fie(ir)*pi/180.0d0
100   dfi(ir)=dfi(ir)*pi/180.0d0
      read(9,*)nf
      if(nf.gt.nf0)call report('too many points')
      close(9)

c     load in polarizability normal mode Fourier coefficients
      do 103 in=1,ncos(1)+nsin(1)
      do 103 jn=1,ncos(2)+nsin(2)
      injn=in+(jn-1)*(ncos(1)+nsin(1))
      write(s9,499)injn
 499  format(i9)
      open(15,file=s9//'.u.ten')
      read(15,*) NAT
      if(NAT.gt.N0)call report('too many atoms')
      DO 106 L=1,NAT*3
106   read(15,*) (APTu(injn,L,II),II=1,3)
      DO 2244 L=1,NAT*3
2244  read(15,*) (AATu(injn,L,II),II=1,3)
      CLOSE(15)

      open(2,file=s9//'.u.ttt')
      read(2,*)
      read(2,*)
      read(2,*)
      read(2,*)
      DO 15 II=1,3
      READ(2,*)
      DO 15 IND=1,NAT*3
15    READ(2,*)LDUM,LDUM,(ALu(injn,IND,II,J),J=1,3)
      READ(2,*)
      READ(2,*)
      DO 13 II=1,3
      READ(2,*)
      DO 13 IND=1,NAT*3
13    READ(2,*)LDUM,LDUM,(Gu(injn,IND,II,J),J=1,3)
      READ(2,*)
      READ(2,*)
      DO 14 II=1,3
      DO 14 J=1,3
      READ(2,*)
      DO 14 IND=1,NAT*3
14    READ(2,*)LDUM,LDUM,(Au(injn,IND,II,J,K),K=1,3)
      close(2)
      an1=0.0d0
      an2=0.0d0
      do 151 IND=1,NAT*3
      do 151 II=1,3
      an1=an1+APTu(injn,IND,II)**2
      DO 151 J=1,3
151   an2=an2+ALu(injn,IND,II,J)**2
      if(in.le.ncos(1))then
       if(jn.le.ncos(2))then
        write(6,6009)in,jn,'cos',wc(1,in),'cos',wc(2,jn),an1,an2
6009    format(2i3,' ',a3,'(',f2.0,' x f) x ',a3,'(',f2.0,' x p), 
     1  norm1=',g9.4,' norm2=',g9.4)
       else
        write(6,6009)in,jn,'cos',wc(1,in),'sin',ws(2,jn-ncos(2)),an1,an2
       endif
      else
       if(jn.le.ncos(2))then
        write(6,6009)in,jn,'sin',wc(1,in-ncos(1)),'cos',wc(2,jn),an1,an2
       else
        write(6,6009)in,jn,'sin',wc(1,in-ncos(1)),'sin',ws(2,jn-ncos(2))
     1  ,an1,an2
       endif
      endif
103   continue

      rc1=5.537d0
      rc2=5.537d0
c     V3C:
      v31=73.0d0
      v32=79.0d0
      mmax=5
      T=298.0d0
      Etol=140.0d0
      plim=1.0d-3
      nmin=0
      nmax=0
c     ldiag - diagonalize rotational Hamiltonian,
c     otherwise use free rotor
      ldiag=.true.
c     ldia2 - diagonalize each vibr state i - incorporate
c     dependence of harmonic frequency on rotation
      ldia2=.true.
      lbath=.false.
      lscale=.false.
      open(9,file='CH3N.OPT')
 111  read(9,900,end=599,err=599)s9
 900  format(a9)
      if(s9(1:4).eq.'ITE1')read(9,*)rc1
      if(s9(1:4).eq.'ITE2')read(9,*)rc2
      if(s9(1:4).eq.'MMAX')read(9,*)mmax
      if(s9(1:4).eq.'NMAX')read(9,*)nmax
      if(s9(1:4).eq.'NMIN')read(9,*)nmin
      if(s9(1:4).eq.'TEMP')read(9,*)T
      if(s9(1:4).eq.'ETOL')read(9,*)Etol
      if(s9(1:4).eq.'DIAG')read(9,*)ldiag
      if(s9(1:3).eq.'V31')read(9,*)v31
      if(s9(1:3).eq.'V32')read(9,*)v32
      if(s9(1:4).eq.'DIA2')read(9,*)ldia2
      if(s9(1:4).eq.'BATH')read(9,*)lbath
      if(s9(1:4).eq.'PLIM')read(9,*)plim
      if(s9(1:4).eq.'SCAL')read(9,*)lscale
      goto 111
 599  close(9)
c
      if(lbath)write(6,*)'Last six freqs at F.INP used as a bath.'
c
c     generate basis of rotational states A x cos/sin(m f1) x cos/sin(m f2)
      write(6,*)(2*mmax+1)**2, ' rotational states'
      i=0
      do 101 i1=1,2*mmax+1
      if(i1.eq.1)then
       m1=0
       im1=1
       a1=1.0d0/dsqrt(2.0d0*pi)
      else
       a1=1.0d0/dsqrt(pi)
       if(i1.le.1+mmax)then
        im1=1
        m1=i1-1
       else
        im1=2
        m1=i1-1-mmax
       endif
      endif
      do 101 i2=1,2*mmax+1
      if(i2.eq.1)then
       m2=0
       im2=1
       a2=1.0d0/dsqrt(2.0d0*pi)
      else
       a2=1.0d0/dsqrt(pi)
       if(i2.le.1+mmax)then
        im2=1
        m2=i2-1
       else
        im2=2
        m2=i2-1-mmax
       endif
      endif
      i=i+1
      if(i.gt.nh0)call report('Hr dimension too big')
      mr(i)=m1
      nrs(i)=m2
      imr(i)=im1
      inr(i)=im2
101   anr(i)=a1*a2
c
      nr=i
      write(6,*)'Rotational Hamiltonian dimension:',nr

      do 1042 i=1,nr
      do 1042 j=1,nr
1042  hr(i,j)=Z0

      if(ldia2)then
       do 1043 i=1,nr
       do 1043 j=1,nr
       do 1043 iq=1,3*NAT
1043   hri(iq,i,j)=Z0
       open(22,file='u.finp')
       read(22,*)
       do  109 iq=1,3*NAT
       read(22,*)
       do  109 jn=1,ncos(2)+nsin(2)
109    read(22,*)(wu(in+(jn-1)*(ncos(1)+nsin(1)),iq),
     1 in=1,ncos(1)+nsin(1))
       write(6,*)'u.finp read'
       close(22)
      endif
C
c     get harmonic frequencies
      call READSW(W)
      do 1091 iq=1,3*NAT
      write(6,1099)iq,wu(1,iq),w(iq)
1099  format(i4,2f10.2)
1091  if(lscale)wu(1,iq)=w(iq)
      if(lscale)write(6,*)'Average replaced by frequencies from F.INP'

c
c     calculate rotational Hamiltonian in a basis of free rotor functions:
      do  104 i=1,nr
      an1=anr(i)
      m1=mr(i)
      n1=nrs(i)
      im1=imr(i)
      in1=inr(i)
      do  104 j=1,nr
      an2=an1*anr(j)
      m2=mr(j)
      n2=nrs(j)
      im2=imr(j)
      in2=inr(j)
c     H = V31/2(1-cos(3f1))+RC1.n^2+V32/2(1-cos(3f2))+RC2 m^2
c     sin, cos are not true eigenfunction of the rotation,
c     but we neglect it now
      hr(i,j)=hr(i,j)+an2*(
     1 0.5d0*v31*(ame(im1,1,im2,m1,0,m2)-ame(im1,1,im2,m1,3,m2))
     2*ame(in1,1,in2,n1,0,n2)
     1+0.5d0*v32*(ame(in1,1,in2,n1,0,n2)-ame(in1,1,in2,n1,3,n2))
     2*ame(im1,1,im2,m1,0,m2)
     3+(rc1* dble(m2)**2+rc2* dble(n2)**2)
     4*ame(im1,1,im2,m1,0,m2)*ame(in1,1,in2,n1,0,n2))
      if(ldia2)then
       do 1044 iq=1,3*NAT
       hri(iq,i,j)=hr(i,j)
       do 1044 in=1,ncos(1)+nsin(1)
       if(in.le.ncos(1))then
        ip=nint(wc(1,in))
        iu=1
       else
        ip=nint(ws(1,in-ncos(1)))
        iu=2
       endif
       do 1044 jn=1,ncos(2)+nsin(2)
       injn=in+(jn-1)*(ncos(1)+nsin(1))
       if(jn.le.ncos(2))then
        jp=nint(wc(2,jn))
        ju=1
       else
        jp=nint(ws(2,jn-ncos(2)))
        ju=2
       endif
1044   if(in.ne.1.or.jn.ne.1)hri(iq,i,j)=hri(iq,i,j)+an2*
     1 wu(injn,iq)*ame(im1,iu,im2,m1,ip,m2)*ame(in1,ju,in2,n1,jp,n2)
      endif
104   continue
      do 108 i=1,nr
108   write(6,610)(hr(i,j),j=1,nr)
610   format(10f8.2)

      if(ldiag)then
       CALL TRED12(nr,hr,sr,qr,2,IERR,EM,nh0)
       if(IERR.ne.0)call report('Diagonalization error')
       write(6,*)'Rotational Energies:'
       write(6,610)(qr(i),i=1,nr)
       if(ldia2)then 
        do 1047 iq=1,3*NAT
	write(6,*)' Mode ',iq
	do 1046 i=1,nr
	do 1046 j=1,nr
1046    ht(i,j)=hri(iq,i,j)
        CALL TRED12(nr,ht,st,qt,2,IERR,EM,nh0)
        if(IERR.ne.0)call report('Diagonalization error')
        write(6,*)'Rotational Energies:'
        write(6,610)(qt(i),i=1,nr)
	do 1047 i=1,nr
	qri(iq,i)=qt(i)
	do 1047 j=1,nr
1047    sri(iq,i,j)=st(i,j)
       endif
      else
       do 107 i=1,nr
  107  qr(i)=hr(i,i)
      endif

      open(20,file='DOG.TAB')
      open(30,file='ROA.TAB')
      WRITE (20,501)
  501 FORMAT (/,'          FREQ.     ',
     1'  D [D^2]     ','  R [D^-2]DO         R - CO      I[km/mol]')
      write(20,886)
      WRITE(30,3304)
3304  FORMAT(
     1'MODE   FREQ    Ramx (180)  Ramy    Ramtot   CID90',
     2'     CID-X   CID180  DEP   ROA180 RamanDCPI',/,
     3'       cm-1       A^4/AMU                    x10^4',
     4'     x10^4     x10^4      A^5/AMUx1e4 au.10^4')
      write(30,886)
C
c
c     calculate partition function:
      EkT=T*8.31441d0/4184.0d0
      cmtokcal=627.5d0/219470.0d0

      Qo=0.0d0
      do 26 ir=1,nr
26    Qo=Qo+exp(-qr(ir)*cmtokcal/EkT)

      Qv=0.0d0
c     vibrational ground state:
      E=0.0d0
      Qv=Qv+1.0d0
c     first vibrational excitation
      do 21 i1=1,3*NAT
      if(w(i1).gt.Etol.or.(lbath.and.i1.gt.3*NAT-6))then
      
      E=E+w(i1)
      Qv=Qv+exp(-E*cmtokcal/EkT)
c     second vibrational excitation
      do 22 i2=i1,3*NAT
      if(w(i2).gt.Etol.or.(lbath.and.i2.gt.3*NAT-6))then
      E=E+w(i2)
      Qv=Qv+exp(-E*cmtokcal/EkT)
c     third vibrational excitation
      do 23 i3=i2,3*NAT
      if(w(i3).gt.Etol.or.(lbath.and.i3.gt.3*NAT-6))then
      E=E+w(i3)
      Qv=Qv+exp(-E*cmtokcal/EkT)
c     fourth vibrational excitation
      do 24 i4=i3,3*NAT
      if(w(i4).gt.Etol.or.(lbath.and.i4.gt.3*NAT-6))then
      E=E+w(i4)
      Qv=Qv+exp(-E*cmtokcal/EkT)
c     fifth vibrational excitation
      do 25 i5=i4,3*NAT
      if(w(i5).gt.Etol.or.(lbath.and.i5.gt.3*NAT-6))then
      E=E+w(i5)
      Qv=Qv+exp(-E*cmtokcal/EkT)
      E=E-w(i5)
      endif
25    continue
      E=E-w(i4)
      endif
24    continue
      E=E-w(i3)
      endif
23    continue
      E=E-w(i2)
      endif
22    continue
      E=E-w(i1)
      endif
21    continue

      Q=Qo*Qv
      write(6,6500)' approx. 1',Q,Qo,Qv
6500  format(' Partition functions tot rot vib, ',A10,': ',3f14.3)

      Qar=(pi*EkT/rc1/rc2/cmtokcal)*exp(-qr(i)*cmtokcal/EkT)
      Qav=1.0d0
      do 3 iq=1,3*NAT
      if(w(iq).gt.Etol.or.(lbath.and.iq.gt.3*NAT-6))then
       Ekcal=w(iq)*cmtokcal
       Qav=Qav/(1.0d0-exp(-Ekcal/EkT))
      endif
3     continue
      write(6,6500)' approx. 2',Qar*Qav,Qar,Qav
c
c     generate all initial states, as above:
c     and calculate spectra
      iskip=0
      itran=0
      itb=0
      E=0.0d0
      Qc=0.0d0
      do 4 iq=1,3*NAT
4     state(iq)=0
      do 30 ir=1,nr
      m=mr(ir)
      n=nrs(ir)
      im=imr(ir) 
      in=inr(ir) 
      an=anr(ir)
      call addqs(E,W,EkT,state,Q,Qc,m,im,n,in,APTu,AATu,wc,ws,plim,
     1ncos,nsin,Gu,ALu,Au,NAT,nw0,iskip,itran,Etol,ldia2,qri,sri,
     2nmin,nmax,0 ,mr,imr,nrs,inr,qr,sr,nh0,anr,an,ldiag,ir,nr,ic)
      if(ic.ne.0)goto 999
      do 31 i1=1,3*NAT
      if(w(i1).gt.Etol.or.(lbath.and.i1.gt.3*NAT-6))then
      call addqs(E,W,EkT,state,Q,Qc,m,im,n,in,APTu,AATu,wc,ws,plim,
     1ncos,nsin,Gu,ALu,Au,NAT,nw0,iskip,itran,Etol,ldia2,qri,sri,
     2nmin,nmax,i1,mr,imr,nrs,inr,qr,sr,nh0,anr,an,ldiag,ir,nr,ic)
      if(ic.ne.0)goto 999
      do 32 i2=i1,3*NAT
      if(w(i2).gt.Etol.or.(lbath.and.i2.gt.3*NAT-6))then
      call addqs(E,W,EkT,state,Q,Qc,m,im,n,in,APTu,AATu,wc,ws,plim,
     1ncos,nsin,Gu,ALu,Au,NAT,nw0,iskip,itran,Etol,ldia2,qri,sri,
     2nmin,nmax,i2,mr,imr,nrs,inr,qr,sr,nh0,anr,an,ldiag,ir,nr,ic)
      if(ic.ne.0)goto 999
      do 33 i3=i2,3*NAT
      if(w(i3).gt.Etol.or.(lbath.and.i3.gt.3*NAT-6))then
      call addqs(E,W,EkT,state,Q,Qc,m,im,n,in,APTu,AATu,wc,ws,plim,
     1ncos,nsin,Gu,ALu,Au,NAT,nw0,iskip,itran,Etol,ldia2,qri,sri,
     2nmin,nmax,i3,mr,imr,nrs,inr,qr,sr,nh0,anr,an,ldiag,ir,nr,ic)
      if(ic.ne.0)goto 999
      do 34 i4=i3,3*NAT
      if(w(i4).gt.Etol.or.(lbath.and.i4.gt.3*NAT-6))then
      call addqs(E,W,EkT,state,Q,Qc,m,im,n,in,APTu,AATu,wc,ws,plim,
     1ncos,nsin,Gu,ALu,Au,NAT,nw0,iskip,itran,Etol,ldia2,qri,sri,
     2nmin,nmax,i4,mr,imr,nrs,inr,qr,sr,nh0,anr,an,ldiag,ir,nr,ic)
      if(ic.ne.0)goto 999
      do 35 i5=i4,3*NAT
      if(w(i5).gt.Etol.or.(lbath.and.i5.gt.3*NAT-6))then
      call addqs(E,W,EkT,state,Q,Qc,m,im,n,in,APTu,AATu,wc,ws,plim,
     1ncos,nsin,Gu,ALu,Au,NAT,nw0,iskip,itran,Etol,ldia2,qri,sri,
     2nmin,nmax,i5,mr,imr,nrs,inr,qr,sr,nh0,anr,an,ldiag,ir,nr,ic)
      if(ic.ne.0)goto 999
      state(i5)=state(i5)-1
      E=E-w(i5)
      endif
35    continue
      state(i4)=state(i4)-1
      E=E-w(i4)
      endif
34    continue
      state(i3)=state(i3)-1
      E=E-w(i3)
      endif
33    continue
      state(i2)=state(i2)-1
      E=E-w(i2)
      endif
32    continue
      state(i1)=state(i1)-1
      E=E-w(i1)
      endif
31    continue
30    continue
999   write(6,6500)'     check',Qc

      write(20,886)
      write(30,886)
886   format(80(1h-))
      close(20)
      close(30)
 
      write(6,889)iskip,itran,itb
889   format(' Transitions from ',i10,' high states skipped',/,
     1I10,' transitions',I10,' recorded')

      write(6,*)'ordering'
      ntb=itb
      do 36 itb=1,ntb
      do 36 ip=itb+1,ntb
      if(et(ip).gt.et(itb))then
       t=et(ip)
       et(ip)=et(itb)
       et(itb)=t
       t=rt(ip)
       rt(ip)=rt(itb)
       rt(itb)=t
       do 37 iq=1,3*nat
       it=s1(iq,ip)
       s1(iq,ip)=s1(iq,itb)
       s1(iq,itb)=it
       it=s2(iq,ip)
       s2(iq,ip)=s2(iq,itb)
37     s2(iq,itb)=it
       do 39 ii=1,8
       it=i44(ii,ip)
       i44(ii,ip)=i44(ii,itb)
39     i44(ii,itb)=it
      endif
36    continue

      open(9,file='T.LST')
      do 38 itb=1,ntb
      write(9,603)et(itb),rt(itb)
603   format(f8.2,f10.3,$)
      do 41 i=1,3*nat
41    if(s1(i,itb).gt.0)write(9,601)i,s1(i,itb)
      write(9,602)
602   format('->',$)
      do 40 i=1,3*nat
40    if((w(i).gt.Etol.or.(lbath.and.i.gt.3*NAT-6))
     1.and.s2(i,itb).gt.0)write(9,601)i,s2(i,itb)
601   format(i3,'_',i1,' ',$)
38    write(9,604)(i44(ii,itb),ii=1,8)
604   format(8i4)
      close(9)
      write(6,*)'Transition list T.LST written.'

      stop
      end


      SUBROUTINE REPORT(RS)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER*(*) RS
      WRITE(6,3000)
3000  FORMAT(/,80(1H*))
      WRITE(6,*)RS
      WRITE(6,3001)
3001  FORMAT(80(1H*),/,/,'PROGRAM STOPPED')
      STOP
      END

      SUBROUTINE READSW(W)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION W(*)
      OPEN(4,FILE='F.INP')
      read(4,*)nint,N
      NAT=N/3
      do 1 i=1,NAT
1     read(4,*)
      read(4,*)
      DO 2 I=1,NAT
      DO 2 J=1,NINT
2     read(4,*)
      read(4,*)nint
      READ(4,4000)(W(I),I=1,NINT)
4000  FORMAT(6F11.6)
      DO 7 I=NINT+1,N
7     W(I)=0.0d0
      CLOSE(5)
      WRITE(6,*)' F.INP read'
      RETURN
      end

      subroutine addqs(E,W,EkT,state,Q,Qc,m,im,n,nn,APTu,AATu,wc,ws,
     1plim,
     1ncos,nsin,Gu,ALu,Au,NAT,nw0,iskip,itran,Etol,ldia2,qri,sri,
     2nmin,nmax,ia,mr,imr,nrs,inr,qr,sr,nh0,anr,an1,ldiag,ir,nr,ic)
c
c     add excitation and get all allowed transitions
c
      implicit none
      integer*4 m,im,i,ii,imp,in,inp,ix,jj,jx,k,kx,l,mp,nat,
     1ncos(*),nsin(*),state(*),nw0,iskip,itran,ip,nmin,nmax,ia,nh0,
     1mr(*),imr(*),ir,iy,jy,irp,nt0,N0,ic,nrs(*),inr(*),nn,n,
     1injn,jn,jnp,np,jp,nr,nd0
      parameter (nt0=100000,N0=20,nd0=2)
      real*8 E,W(*),EkT,Q,P,Ekcal,APTt(3),AATt(3),wc(nd0,nw0),
     1ws(nd0,nw0),Qc,Stol,
     3ALt(3,3),Gt(3,3),At(3,3,3),APTu(nw0**nd0,N0*3,3),
     1AATu(nw0**nd0,N0*3,3),
     3ALu(nw0**nd0,N0*3,3,3),Gu(nw0**nd0,N0*3,3,3),
     1Au(nw0**nd0,N0*3,3,3,3),ame,amu,
     1beta2,betaa,betag,bohr,cid1,cid2,d90x,d90z,dx,dz,
     1ep,eps,fc,fd,ff,fr,gpisvejc,PI,ram1,ram2,roa1,
     1roa2,roa3,s90x,s90z,sa1,sag0,sag1,sal1,spa3,spal,sumd,wr,
     1ydx,ydy,rdo,sal0,p1,Etol,Ep0,an1,qr(*),sr(nh0,nh0),anr(*),
     1an2,ci,cj,Eg,sri(3*N0,nh0,nh0),qri(3*N0,nh0),plim
      logical ldiag,ldia2
      integer*4 s1(3*N0,nt0),s2(3*N0,nt0),itb,iq,i44(8,nt0)
      real*8 et(nt0),rt(nt0)
      logical lbath
      common/trb/s1,s2,et,rt,itb,i44,lbath
c
      ic=0
      if(ia.gt.0)then
       state(ia)=state(ia)+1
       E=E+W(ia)
      endif

      Eg=E+qr(ir)

c
      pi=4.0d0*atan(1.0d0)
      Ekcal=Eg/219470.0d0*627.5d0
      FD=dsqrt(388.9219d0)
      FR=2586.285D-8
      FC=131.14D08
      P=exp(-Ekcal/EkT)/Q
      Qc=Qc+exp(-Ekcal/EkT)
c
      if((P.lt.plim.and..not.lbath).or.P.lt.1.0d-4)then
       iskip=iskip+1
       return
      endif
c
      write(6,600)
600   format(1h<,$)
      do 2 i=1,3*nat
2     if(w(i).gt.Etol.or.(lbath.and.i.gt.3*NAT-6))write(6,601)state(i)
601   format(i1,$)
      write(6,602)m,im,n,nn,Eg,P
602   format(1h|,1h<,4i1,1h|,f8.2,' cm-1, p = ',f6.3)

c
c     look at excited states |state> --> |state*>
c     harmonic transitions only (n->n+1)
      do 1 i=1,3*nat
      if(w(i).gt.Etol.or.(lbath.and.i.gt.3*NAT-6))then
      if(i.ge.nmin.or.nmin.eq.0)then
      if(i.le.nmax.or.nmax.eq.0)then
      state(i)=state(i)+1
      Ep0=E+w(i)
      
      do 10 irp=1,nr
      imp=imr(irp)
      mp=mr(irp)
      inp=inr(irp)
      np=nrs(irp)
      an2=anr(irp)
      Ep=Ep0+qr(irp)
      if(ldia2)Ep=Ep0+qri(i,irp)
c
c     transition tensors
      do 4 ix=1,3
      APTt(ix)=0.0d0
      AATt(ix)=0.0d0
      do 4 jx=1,3
      ALt(ix,jx)=0.0d0
      Gt(ix,jx)=0.0d0
      do 4 kx=1,3
4     At(ix,jx,kx)=0.0d0

      if(ldiag)then
       do 61 iy=1,nr
       ci=sr(iy,ir)
       if(ldia2)ci=sri(i,iy,ir)
       an1=anr(iy)
       do 61 jy=1,nr
       cj=ci*sr(jy,irp)
       if(ldia2)cj=ci*sri(i,jy,irp)
       an2=anr(jy)
       if(dabs(cj).gt.1.0d-6)then
        do 62 in=1,ncos(1)+nsin(1)
        if(in.le.ncos(1))then
         ip=nint(wc(1,in))
         inp=1
        else
         ip=nint(ws(1,in-ncos(1)))
         inp=2
        endif
        do 62 jn=1,ncos(2)+nsin(2)
        injn=in+(jn-1)*(ncos(1)+nsin(1))
        if(jn.le.ncos(2))then
         jp=nint(wc(2,jn))
         jnp=1
        else
         jp=nint(ws(2,jn-ncos(2)))
         jnp=2
        endif
        ff=ame(imr(iy),inp,imr(jy),mr(iy),ip,mr(jy))*cj*an1*an2
     1  *ame(inr(iy),jnp,inr(jy),nrs(iy),jp,nrs(jy))
     2  *dsqrt(dble(state(i)))
        do 62 ix=1,3
        APTt(ix)=APTt(ix)+APTu(injn,i,ix)*ff
        AATt(ix)=AATt(ix)+AATu(injn,i,ix)*ff
        do 62 jx=1,3
        ALt(ix,jx)=ALt(ix,jx)+ALu(injn,i,ix,jx)*ff
        Gt(ix,jx)=Gt(ix,jx)+Gu(injn,i,ix,jx)*ff
        do 62 kx=1,3
62      At(ix,jx,kx)=At(ix,jx,kx)+Au(injn,i,ix,jx,kx)*ff
       endif
61     continue
      else
       do 6 in=1,ncos(1)+nsin(1)
       if(in.le.ncos(1))then
        ip=nint(wc(1,in))
        inp=1
       else
        ip=nint(ws(1,in-ncos(1)))
        inp=2
       endif
       do 6 jn=1,ncos(2)+nsin(2)
       injn=in+(jn-1)*(ncos(1)+nsin(1))
       if(jn.le.ncos(2))then
        jp=nint(wc(2,jn))
        jnp=1
       else
        jp=nint(ws(2,jn-ncos(2)))
        jnp=2
       endif
       ff=ame(im,inp,imp,m,ip,mp)*an1*an2*ame(nn,jnp,inp,n,jp,np)
     1 *dsqrt(dble(state(i)))
       do 6 ix=1,3
       APTt(ix)=APTt(ix)+APTu(injn,i,ix)*ff
       AATt(ix)=AATt(ix)+AATu(injn,i,ix)*ff
       do 6 jx=1,3
       ALt(ix,jx)=ALt(ix,jx)+ALu(injn,i,ix,jx)*ff
       Gt(ix,jx)=Gt(ix,jx)+Gu(injn,i,ix,jx)*ff
       do 6 kx=1,3
6      At(ix,jx,kx)=At(ix,jx,kx)+Au(injn,i,ix,jx,kx)*ff
      endif

      fr=w(i)
      WR=dsqrt(dabs(fr))
      do 7 ix=1,3
      if(WR.gt.0.1d0)APTt(ix)=APTt(ix)/WR*FD
7     AATt(ix)=AATt(ix)*FC*WR
      SUMD=0.0d0
      RDO=0.0d0
      do 8 ix=1,3
      SUMD=SUMD+APTt(ix)**2
8     RDO=RDO+APTt(ix)*AATt(ix)

      SAL0=0.0d0
      SAL1=0.0d0
      SAG0=0.0d0
      SAG1=0.0d0
      SA1=0.0d0
      SPAL=0.0d0
      DO 9 II=1,3
      SPAL=SPAL+ALt(II,II)
      DO 9 JJ=1,3
      SAL0=SAL0+ALt(II,II)*ALt(JJ,JJ)
      SAL1=SAL1+ALt(II,JJ)*ALt(II,JJ)
      SAG0=SAG0+ALt(II,II)*    Gt(JJ,JJ)
      SAG1=SAG1+ALt(II,JJ)*    Gt(II,JJ)
      DO 9 K=1,3
      DO 9 L=1,3
9     SA1=SA1+EPS(II,K,L)*ALt(II,JJ)*At(K,L,JJ)/3.0d0
      S90X=7.0d0*SAL1+SAL0
      D90X=2.0d0*(7.0d0*SAG1+SAG0+SA1)
      DX=D90X
      IF(S90X.gt.1.0d-9)DX=DX/S90X
c     IL+/-IR, 90o, z
      S90Z=2.0d0*(3.0d0*SAL1-SAL0)
      D90Z=4.0d0*(3.0d0*SAG1-SAG0-SA1)
      DZ=D90Z
      IF(S90Z.gt.1.0d-9)DZ=DZ/S90Z
c     polarized backscattering components (YDY+YDX=S180):
      YDY=3.0d0*SAL1-SAL0
      YDX=4.0d0*SAL1+2.0d0*SAL0
      P1=YDY/YDX
      AMU=1822.0d0
      BOHR=0.529177d0
      SpA3=SPAL/3.0d0*dsqrt((AMU*BOHR**4))
      beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4)
      YDY=6.0d0*YDY*(AMU*BOHR**4)
      YDX=6.0d0*YDX*(AMU*BOHR**4)
      gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/4880.0d0
c     alphag=SAG0/9.0d0*gpisvejc
      betaa=SA1*3.0d0/2.0d0*gpisvejc
      betag=(3.0d0*SAG1-SAG0)/2.0d0*gpisvejc
c     DCP 180
c     ram3=24.0d0*beta2
      roa3=8.0d0*(12.0d0*betag+4.0d0*beta2)
      S90X=S90X*(AMU*BOHR**4)
      D90X=D90x*gpisvejc
      DX=D90X
      IF(S90X.gt.1.0d-9)DX=DX/S90X
c     ICP(90)=90Z
      roa2=8.0d0*(3.0d0*betag-betaa)
      ram2=12.0d0*beta2
      CID2=0.0d0
      if(ram2.gt.1.0d-9)CID2=roa2/ram2
c     ICP(180):
      roa1=96.0d0*(betag+betaa/3.0d0)
      ram1=4.0d0*(45.0d0*SpA3**2+7.0d0*beta2)
      CID1=0.0d0
      if(ram1.gt.1.0d-9)CID1=roa1/ram1
      ram1=ram1*P
      roa1=roa1*P
      roa3=roa3*P
      YDX=YDX*P
      YDY=YDY*P
      Stol=0.0d-6
      if(ram1.gt.Stol.or.SUMD.gt.Stol)then
       itran=itran+1
       if(ram1.gt.Stol)
     1 WRITE(30,3300)itran,Ep-Eg,YDX,YDY,ram1,CID2,DX,CID1,P1,roa1,roa3
3300   FORMAT(I8,f9.2,6g13.4,f6.3,2g11.4)
       SUMD=SUMD*P
       RDO=RDO*P
       if(SUMD.gt.Stol)
     1 write(20,523)itran,Ep-Eg,SUMD,RDO,RDO,dsqrt(2.0d0*pi)*fr*SUMD
523    format(i8,5G15.6)
       itb=itb+1
       if(itb.gt.nt0)then
        write(6,*)'Too many transitions'
	ic=1
	return
       endif
       do 11 iq=1,3*nat
       s1(iq,itb)=state(iq)
11     s2(iq,itb)=state(iq)
       s1(i,itb)=s1(i,itb)-1
       et(itb)=Ep-Eg
       rt(itb)=ram1
       i44(1,itb)=m
       i44(2,itb)=im
       i44(3,itb)=mp
       i44(4,itb)=imp
       i44(5,itb)=n
       i44(6,itb)=nn
       i44(7,itb)=np
       i44(8,itb)=inp
      endif

10    continue

      state(i)=state(i)-1
      endif
      endif
      endif
1     continue
      return
      end

      function ame(ii,jj,kk,i1,j1,k1)
c     int(cos/sin(if)*cos/sin(jf)*cos/sin(kf),f=0..2pi) 
      implicit none
      real*8 ame,pi,sum
      integer*4 ii,jj,kk,i1,j1,k1
      sum=0.0d0
      pi=4.0d0*atan(1.0d0)
      if(ii.eq.1.and.jj.eq.1.and.kk.eq.2)goto 99
      if(ii.eq.1.and.jj.eq.2.and.kk.eq.1)goto 99
      if(ii.eq.2.and.jj.eq.1.and.kk.eq.1)goto 99
      if(ii.eq.2.and.jj.eq.2.and.kk.eq.2)goto 99
      if(ii.eq.1.and.jj.eq.1.and.kk.eq.1)then
       if(i1+j1+k1.eq.0)sum=sum+pi/2.0d0
       if(i1+j1-k1.eq.0)sum=sum+pi/2.0d0
       if(i1-j1+k1.eq.0)sum=sum+pi/2.0d0
       if(i1-j1-k1.eq.0)sum=sum+pi/2.0d0
       goto 99
      endif
      if(ii.eq.1.and.jj.eq.2.and.kk.eq.2)then
       if(i1+j1+k1.eq.0)sum=sum-pi/2.0d0
       if(i1+j1-k1.eq.0)sum=sum+pi/2.0d0
       if(i1-j1+k1.eq.0)sum=sum+pi/2.0d0
       if(i1-j1-k1.eq.0)sum=sum-pi/2.0d0
       goto 99
      endif
      if(ii.eq.2.and.jj.eq.1.and.kk.eq.2)then
       if(i1+j1+k1.eq.0)sum=sum-pi/2.0d0
       if(i1+j1-k1.eq.0)sum=sum+pi/2.0d0
       if(i1-j1+k1.eq.0)sum=sum-pi/2.0d0
       if(i1-j1-k1.eq.0)sum=sum+pi/2.0d0
       goto 99
      endif
      if(ii.eq.2.and.jj.eq.2.and.kk.eq.1)then
       if(i1+j1+k1.eq.0)sum=sum-pi/2.0d0
       if(i1+j1-k1.eq.0)sum=sum-pi/2.0d0
       if(i1-j1+k1.eq.0)sum=sum+pi/2.0d0
       if(i1-j1-k1.eq.0)sum=sum+pi/2.0d0
       goto 99
      endif
      write(6,*)ii,jj,kk,i1,j1,k1
      call report('ame-unknown combination')
 99   ame=sum
      return
      end

      FUNCTION EPS(I,J,K)
      REAL*8 EPS
      EPS=0.0D0
      IF (I.EQ.1.AND.J.EQ.2.AND.K.EQ.3)EPS= 1.0D0
      IF (I.EQ.1.AND.J.EQ.3.AND.K.EQ.2)EPS=-1.0D0
      IF (I.EQ.2.AND.J.EQ.3.AND.K.EQ.1)EPS= 1.0D0
      IF (I.EQ.2.AND.J.EQ.1.AND.K.EQ.3)EPS=-1.0D0
      IF (I.EQ.3.AND.J.EQ.1.AND.K.EQ.2)EPS= 1.0D0
      IF (I.EQ.3.AND.J.EQ.2.AND.K.EQ.1)EPS=-1.0D0
      RETURN
      END

      SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR,E,MX3)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(MX3,N),Z(MX3,N),D(MX3),E(MX3)
      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,MX3)
      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,MX3)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION Z(MX3,N),D(MX3),E(MX3)
      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) .GE. 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