PROGRAM NEW6W
c     Tensors with multiple excitation frequencies in FILEW.TTT
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      real*8,allocatable::S(:,:),ALPHA(:,:,:,:),A(:,:,:,:,:),
     1G(:,:,:,:),ALPHAV(:,:,:,:),E(:),w(:),GV(:,:,:,:),
     1alphaw(:,:,:),GPw(:,:,:),Aw(:,:,:,:),
     1alphavw(:,:,:),GPvw(:,:,:),
     1ALPHA1(:,:,:),G1(:,:,:),A1(:,:,:,:),
     1ALPHA1V(:,:,:),G1V(:,:,:),
     1b2(:),alpha2(:),aG(:),bA(:),bG(:)
      LOGICAL RAMAN,LPOLAR,lpro,lgpro,lapro,
     1luseA,luseG,lcomplex
      integer*4 units
      character*10 s10
      CM=219470.0d0
C
      OPEN(3,FILE='ROA.OUT')
      WRITE(6,6666)
      WRITE(3,6666)
6666  FORMAT(/,/,
     A'          * * * * * * * * * * * * * * * * * * * * * *',/,
     1'        *                                         *  ',/,
     2'      *           RAMAN OPTICAL ACTIVITY        *    ',/,
     3'        *                                         *  ',/,
     4'          * * * * * * * * * * * * * * * * * * * * * *',/,
     8'                                                     ',/,
     8'                       Petr Bour 1997-2018           ',/,/,/,
     A'INPUT FILES: C5R.OUT, ROA.OPT, FILE.XYZ, FILE.TTT,F.INP',/,
     B'OUTPUT     : ROA.OUT, ROA.TAB ROA.TTT',/)

      CALL IOPAR(RAMAN,N1,N2,LPOLAR,
     1wmin,wmax,lpro,lgpro,lapro,N6,N7,IDIF,units,temp,ri,IKIND,
     1luseA,luseG,lcomplex)

      call getatw(NAT,np)
      N=3*NAT
      allocate(S(N,N),ALPHA(np,2*N,3,3),A(np,2*N,3,3,3),
     1G(np,2*N,3,3),ALPHAV(np,2*N,3,3),E(N),w(np),GV(np,2*N,3,3))
      CALL READTEN(np,N,NAT,ALPHA,ALPHAV,G,GV,A,w,lcomplex)
      CALL READS(N,S,E,M)
      write(6,*)N,' modes',np,' frequencies'
      if(iargc().gt.1)then
       call getarg(1,s10)
       read(s10,*)IQ
       call getarg(2,s10)
       read(s10,*)IR
      else
       write(6,605)
605    format(' Input mode and frequency of interest (M,N): ',$)
       read(5,*)IQ,IR
      endif
      write(6,609)IQ,E(IQ)*CM,IR,W(IR),w(IR)/8065.54093D0,1.0d7/w(IR)
609   format(i4,'. mode ',f10.2,' cm-1',/,
     1       i4,'. wexc ',f10.2,' cm-1',f10.2,' eV ',f10.2,' nm')
C
      CALL WRITETTT(N,np,IR,NAT,ALPHA,A,G,ALPHAV)
C
c     invariants for selected mode:
C     transform the tensors into IQ internal derivatives
      allocate(alphaw(2*np,3,3),GPw(2*np,3,3),Aw(2*np,3,3,3),
     1alphavw(2*np,3,3),GPvw(2*np,3,3))
      CALL TRTENii(IQ,np,S,N,ALPHA,ALPHAV,A,G,GV,
     1alphaw,alphavw,GPw,GPvw,Aw)
      allocate(b2(np),alpha2(np),aG(np),bA(np),bG(np))
      call vz(alpha2,np)
      call vz(b2,np)
      call vz(bA,np)
      call vz(aG,np)
      call vz(bG,np)
      call mkv(np,alpha2,aG,b2,bG,bA,alphaw,alphavw,GPw,GPvw,Aw)
      call wrw(IQ,alpha2,np,'alpha2',w)
      call wrw(IQ,aG    ,np,'aG'    ,w)
      call wrw(IQ,b2    ,np,'b2'    ,w)
      call wrw(IQ,bG    ,np,'bG'    ,w)
      call wrw(IQ,bA    ,np,'bA'    ,w)

c     spectrum for selected frequency:
C     transform the tensors into internal coordinate derivatives
      allocate(A1(2*N,3,3,3),ALPHA1(2*N,3,3),G1V(2*N,3,3),
     1ALPHA1V(2*N,3,3),G1(2*N,3,3))
      CALL TRTEN(IR,M,np,S,3*NAT,ALPHA,ALPHAV,A,G,GV,ALPHA1,
     1ALPHA1V,G1,G1V,A1)
      EXCNM=1.0d7/w(IR)
      EXCA=EXCNM*10.0d0
      CALL DORAMAN(M,ALPHA1,G1,A1,EXCA,wmin,wmax,E,
     1units,temp,ri,IKIND,luseA,luseG)
      END
C     ==========================================
      subroutine mkv(np,alpha2,aG,b2,bG,bA,alw,alv,GPw,GPvw,Aw)
      implicit none
      integer*4 np,iw,i,j,ig,id
      real*8 alpha2(*),aG(*),b2(*),bG(*),bA(*),spal,spG,
     1alw(2*np,3,3),GPw(2*np,3,3),Aw(2*np,3,3,3),alv(2*np,3,3),
     1GPvw(2*np,3,3)
      call vz(b2,np)
      call vz(bA,np)
      call vz(bG,np)
      do 1 iw=1,np
      spal=(alw(iw,1,1)+alw(iw,2,2)+alw(iw,3,3))/3.0d0
c     origin independent: need GPv
      spG =(GPvw(iw,1,1)+GPvw(iw,2,2)+GPvw(iw,3,3))/3.0d0
      do 103 i=1,3
      do 103 j=1,3
      b2(iw)=b2(iw)+1.5d0*alw(iw,i,j)*alw(iw,i,j)
c     origin independent: need alv
      bG(iw)=bG(iw)+1.5d0*alv(iw,i,j)*GPw(iw,i,j)
      ig=i+1
      if(ig.gt.3)ig=1
      id=ig+1
      if(id.gt.3)id=1
c     origin independent: need alv??
103   bA(iw)=bA(iw)+0.5d0*alv(iw,i,j)*(Aw(iw,ig,id,j)-Aw(iw,id,ig,j))
      alpha2(iw)=spal*spal
      aG(    iw)=spal*spG
      b2(iw)=b2(iw)-4.5d0*alpha2(iw)
1     bG(iw)=bG(iw)-4.5d0*aG(iw)
      return
      end
C     ==========================================
      SUBROUTINE WRITETTT(N,np,IR,NAT,ALPHA,A,G,ALPHAV)
      IMPLICIT none
      INTEGER*4 N,np,IR,NAT,I,IX,II,L,J,K
      real*8 ALPHA(np,2*N,3,3),G(np,2*N,3,3),A(np,2*N,3,3,3),
     1ALPHAV(np,2*N,3,3)
c
      open(22,file='ROA.TTT')
      WRITE(22,2000)NAT
2000  FORMAT(' ROA tensors, cartesian derivatives',/,I4,' atoms',/,
     1' The electric-dipolar electric-dipolar polarizability:',/,
     2' Atom/x    jx           jy           jz')
      DO 1 I=1,3
      WRITE(22,2002)I
2002  FORMAT(' Alpha(',I1,',J):')
      DO 1 L=1,NAT
      DO 1 IX=1,3
      II=3*(L-1)+IX
1     WRITE(22,2001)L,IX,(ALPHA(IR,II,I,J),J=1,3)
2001  FORMAT(I5,1H ,I1,3F13.7)
      WRITE(22,2004)
2004  FORMAT(' The electric dipole magnetic dipole polarizability:',/,
     1       ' Atom/x    jx(Bx)       jy(By)       jz(Bz)')
      DO 2 I=1,3
      WRITE(22,2003)I
2003  FORMAT(' G(',I1,',J):')
      DO 2 L=1,NAT
      DO 2 IX=1,3
      II=3*(L-1)+IX
c     last inner index magnetic
2     WRITE(22,2001)L,IX,(G(IR,II,I,J),J=1,3)
      WRITE(22,2005)
2005  FORMAT(' The electric dipole electric quadrupole polarizability:',
     2/,     ' Atom/x    kx           ky           kz')
      DO 3 I=1,3
      DO 3 J=1,3
      WRITE(22,2006)I,J
2006  FORMAT(' A(',I1,'(u),',I1,',K):')
      DO 3 L=1,NAT
      DO 3 IX=1,3
      II=3*(L-1)+IX
3     WRITE(22,2001)L,IX,(A(IR,II,I,J,K),K=1,3)
      WRITE(22,2009)
2009  FORMAT('The velocity form of alpha',/,
     1' Atom/x    vx           vy           vz')
      DO 11 I=1,3
      WRITE(22,2012)I
2012  FORMAT(' AlphaV(',I1,',J):')
      DO 11 L=1,NAT
      DO 11 IX=1,3
      II=3*(L-1)+IX
11    WRITE(22,2001)L,IX,(ALPHAV(IR,II,I,J),J=1,3)
      close(22)
      RETURN
      END

C     ================================================================
C     ==============================================================
C     ==============================================================
C     =================================================================
      FUNCTION EPS(I,J,K)
      IMPLICIT INTEGER*4 (I-N)
      REAL*8 EPS
      EPS=0.0D0
      IF(I.EQ.1)then
       if(J.EQ.2.AND.K.EQ.3)EPS= 1.0D0
       IF(J.EQ.3.AND.K.EQ.2)EPS=-1.0D0
      endif
      IF(I.EQ.2)then
       IF(J.EQ.3.AND.K.EQ.1)EPS= 1.0D0
       IF(J.EQ.1.AND.K.EQ.3)EPS=-1.0D0
      endif
      IF(I.EQ.3)then
       IF(J.EQ.1.AND.K.EQ.2)EPS= 1.0D0
       IF(J.EQ.2.AND.K.EQ.1)EPS=-1.0D0
      endif
      RETURN
      END
C     ==============================================================
      SUBROUTINE READS(N3,S,E,M)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION S(N3,N3),E(*)
      CM=219470.0d0
      SFAC=0.0234280d0
      OPEN(4,FILE='F.INP',STATUS='OLD')
      read(4,*)M,nat3,nat
      N=NAT3
      do 1 i=1,NAT+1
1     read(4,*)
      DO 2 I=1,NAT
      DO 2 J=1,M
2     read(4,*)(s(3*(i-1)+ix,j),ix=1,2),
     1(s(3*(i-1)+ix,j),ix=1,3)
      read(4,*)
      READ(4,4000)(E(I),I=1,M)
4000  FORMAT(6F11.6)
      DO 3 I=1,M
3     E(I)=E(I)/CM
      DO 5 I=1,N
      DO 5 J=1,M
5     S(I,J)=S(I,J)*SFAC
      NAT=N/3
      RETURN
      end
C
      SUBROUTINE READTEN(np,N,NAT,ALPHA,ALPHAV,G,GV,A,w,lcomplex)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION ALPHA(np,2*N,3,3),G(np,2*N,3,3),A(np,2*N,3,3,3),w(*)
      DIMENSION ALPHAv(np,2*N,3,3),GV(np,2*N,3,3)
      logical lcomplex
      ALPHA=0.0d0
      ALPHAV=0.0d0
      G=0.0d0
      A=0.0d0
      GV=0.0d0
      OPEN(2,FILE='FILEW.TTT',STATUS='OLD')
      do 4 iw=1,np
      READ(2,*)
      READ(2,200)NAT,w(iw)
200   format(i4,14x,f11.2)
      if(NAT.gt.N/3)call report('too many atoms')
      READ(2,*)
      READ(2,*)
      DO 1 I=1,3
      READ(2,*)
      DO 1 L=1,NAT
      DO 1 IX=1,3
      IND=3*(L-1)+IX
      if(lcomplex)then
       READ(2,*)(ALPHA(iw,IND,I,J),J=1,2),(ALPHA(iw,IND,I,J),J=1,3),
     1 (ALPHA(iw,N+IND,I,J),J=1,3)
      else
       READ(2,*)(ALPHA(iw,IND,I,J),J=1,2),(ALPHA(iw,IND,I,J),J=1,3)
      endif
1     continue
      READ(2,*)
      READ(2,*)
      DO 2 I=1,3
      READ(2,*)
      DO 2 L=1,NAT
      DO 2 IX=1,3
      IND=3*(L-1)+IX
c     last index magnetic, inner
      if(lcomplex)then
       READ(2,*)(G(iw,IND,I,J),J=1,2),(G(iw,IND,I,J),J=1,3),
     1 (G(iw,IND+N,I,J),J=1,3)
      else
       READ(2,*)(G(iw,IND,I,J),J=1,2),(G(iw,IND,I,J),J=1,3)
      endif
2     continue
      READ(2,*)
      READ(2,*)
      DO 3 I=1,3
      DO 3 J=1,3
      READ(2,*)
      DO 3 L=1,NAT
      DO 3 IX=1,3
      IND=3*(L-1)+IX
      if(lcomplex)then
       READ(2,*)(A(iw,IND,I,J,K),K=1,2),(A(iw,IND,I,J,K),K=1,3),
     1 (A(iw,N+IND,I,J,K),K=1,3)
      else
       READ(2,*)(A(iw,IND,I,J,K),K=1,2),(A(iw,IND,I,J,K),K=1,3)
      endif
3     continue
      READ(2,*)
      READ(2,*)
      DO 4 I=1,3
      READ(2,*)
      DO 4 L=1,NAT
      DO 4 IX=1,3
      IND=3*(L-1)+IX
      if(lcomplex)then
       READ(2,*)(ALPHAv(iw,IND,I,J),J=1,2),(ALPHAv(iw,IND,I,J),J=1,3),
     1 (ALPHAv(iw,N+IND,I,J),J=1,3)
      else
       READ(2,*)(ALPHAv(iw,IND,I,J),J=1,2),(ALPHAv(iw,IND,I,J),J=1,3)
      endif
4     continue
      CLOSE(2)
      WRITE(3,*)NAT,' atoms, tensors read in from FILE.TTT'
      RETURN
      END
C     ==============================================================
C
      SUBROUTINE IOPAR(RAMAN,N1,N2,LPOLAR,
     1wmin,wmax,lpro,lgpro,lapro,N6,N7,IDIF,units,temp,ri,IKIND,
     1luseA,luseG,lcomplex)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      LOGICAL RAMAN,LPOLAR,lpro,lex,lgpro,lapro,
     1luseA,luseG,lcomplex
      CHARACTER*15 STRING
      integer*4 units
c
      N7=7
      N1=7
      N2=0
      IDIF=-1
      wmin=0.0d0
      wmax=1.0d90
      lpro=.false.
      lcomplex=.true.
      lapro=.false.
      lgpro=.false.
      RAMAN=.FALSE.
      N6=3
      LPOLAR=.FALSE.
c     excitation light wavelength in nm:
c     units:
c     0  .... Gaussian
c     1  .... Kapitan
c     2  .... old units (Polavarapu)
      units=0
c     temperature (K):
      temp=293.0d0
c     refractive index
      ri=1.0d0
c     kind of ROA:
c     0 - 180 deg SCP (ordinary)
c     1 - 90 deg polarized
c     2 - 90 deg depolarized
      IKIND=0
c     use A-tensor derivatives:
      luseA=.true.
c     use G'-tensor derivatives:
      luseG=.true.

      INQUIRE(FILE='ROA.OPT',exist=lex)
      if(lex)then
       WRITE(6,*)' FILE ROA.OPT FOUND ...'
       OPEN(2,FILE='ROA.OPT')
      else
       INQUIRE(FILE='NEW6.OPT',exist=lex)
       if(lex)then
        OPEN(2,FILE='NEW6.OPT')
        WRITE(6,*)' FILE NEW6.OPT FOUND ...'
       else
        return
       endif
      endif
   
2001  READ(2,1111)STRING
1111  FORMAT(A15)
      WRITE(3,1111)STRING
      IF(STRING(1:4).EQ.'N1N2')THEN
       WRITE(3,*) ' Vibrations selected ...'
       READ(2,*)N1,N2
      ENDIF
      IF(STRING(1:4).EQ.'APRO')READ(2,*)lapro
      IF(STRING(1:4).EQ.'COMP')READ(2,*)lcomplex
      IF(STRING(1:4).EQ.'GPRO')READ(2,*)lgpro
      IF(STRING(1:4).EQ.'KIND')READ(2,*)IKIND
      IF(STRING(1:2).EQ.'N6')READ(2,*)N6
      IF(STRING(1:3).EQ.'PRO')READ(2,*)lpro
      IF(STRING(1:5).EQ.'POLAR')READ(2,*)LPOLAR
      IF(STRING(1:5).EQ.'RAMAN')READ(2,*)RAMAN
      IF(STRING(1:4).EQ.'REFR')READ(2,*)ri
      IF(STRING(1:4).EQ.'TEMP')READ(2,*)temp
      IF(STRING(1:4).EQ.'USEA')READ(2,*)luseA
      IF(STRING(1:4).EQ.'UNIT')READ(2,*)units
      IF(STRING(1:4).EQ.'USEG')READ(2,*)luseG
      IF(STRING(1:4).EQ.'WMIN')READ(2,*)wmin
      IF(STRING(1:4).EQ.'WMAX')READ(2,*)wmax
      IF(STRING(1:3).EQ.'END')GOTO 2002
      GOTO 2001
2002  CLOSE(2)
     
      RETURN
      END

      SUBROUTINE DORAMAN(N,ALPHA,G,A,EXCA,wmin,wmax,E,
     1units,temp,ri,IKIND,luseA,luseG)
c
c     supposedly G/W is supplied
c
      IMPLICIT REAL*8 (A-H,O-Z)
      implicit integer*4 (i-n)
      DIMENSION ALPHA(2*N,3,3),G(2*N,3,3),A(2*N,3,3,3),E(*)
      integer*4 units
      logical luseA,luseG
c
      CM=219470.0d0
      if(.not.luseA)then
       write(6,*)'A-tensor deleted'
       A=0.0d0
      endif

      if(.not.luseG)then
       write(6,*)'G-tensor deleted'
       G=0.0d0
      endif
c
      OPEN(4,FILE='ROA.TAB')
      if(units.eq.0)then
       write(6,*)'Gaussian units'
       WRITE(4,3304)
      elseif(units.eq.1)then
       write(6,*)'Absolute units'
       WRITE(4,3305)
      elseif(units.eq.2)then
       write(6,*)'Arbitrary units'
       WRITE(4,3004)
      else
       call report('Unknown unit option')
      endif

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')
3004  FORMAT(
     1'MODE   FREQ    IDX    IDY    RAMAN     ROA-Z',
     2'     ROA-X   ROA-180  DEP   raw180 RamanDCPI',/,
     3'       cm-1                            x10^4',
     4'     x10^4     x10^4       au.10^4')
3305   FORMAT(
     1'MODE   FREQ    Ramx (180)  Ramy    Ramtot   CID90',
     2'     CID-X   CID180  DEP   ROA180 RamanDCPI',/,
     3'       cm-1     Ram: 10^-12 cm^2/J  Roa: 10^-16 cm^2/J ',
     4'     x10^4     x10^4      A^5/AMUx1e4 au.10^4')

      write(4,3002)
3002  FORMAT(80(1H-))
c
      DO 1 IQ=1,N
      ECM=E(IQ)*CM
      IF(ABS(ECM).GT.0.01d0)THEN
      if(wmin.eq.0.0d0.or.ECM.gt.wmin)then
      if(wmax.eq.0.0d0.or.ECM.lt.wmax)then
       SAL0=0.0d0
       SAL1=0.0d0
       SAG0=0.0d0
       SAG1=0.0d0
       SA1=0.0d0
       DO 3 I=1,3
       DO 3 J=1,3
c      Re ALPHA Re ALPHA + Im ALPHA Im ALPHA, etc.:
       SAL0=SAL0+ALPHA(  IQ,I,I)*ALPHA(  IQ,J,J)
     1          +ALPHA(N+IQ,I,I)*ALPHA(N+IQ,J,J)
       SAL1=SAL1+ALPHA(  IQ,I,J)*ALPHA(  IQ,I,J)
     1          +ALPHA(N+IQ,I,J)*ALPHA(N+IQ,I,J)
       SAG0=SAG0+ALPHA(  IQ,I,I)*    G(  IQ,J,J)
     1          +ALPHA(N+IQ,I,I)*    G(N+IQ,J,J)
       SAG1=SAG1+ALPHA(  IQ,I,J)*    G(  IQ,I,J)
     1          +ALPHA(N+IQ,I,J)*    G(N+IQ,I,J)
       DO 3 K=1,3
       DO 3 L=1,3
c      if(EPS(I,K,L).ne.0.0d0)write(6,*)I,J,K,L,ALPHA(IQ,I,J),
c    1 A(IQ,K,L,J)/3.0d0
3      SA1=SA1+EPS(I,K,L)*(ALPHA(  IQ,I,J)*A(  IQ,K,L,J)
     1                    +ALPHA(N+IQ,I,J)*A(N+IQ,K,L,J))/3.0d0
c
c      IL+IR, backscattering, DCP_I
       SDCP180=3.0d0*SAL1-SAL0
c
c      IL+/-IR, backscattering
       S180=7.0d0*SAL1+SAL0
c      write(6,*)ECM,SAL0,SAL1,SAG0,SAG1,SA1
       D180=8.0d0*(3.0d0*SAG1-SAG0+SA1)
       D=D180
       IF(S180.gt.1.0d-9)D=D/S180
c
c      IL+/-IR, 90o, x
       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
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
c      polarized backscattering components (YDY+YDX=S180):
       YDY=3.0d0*SAL1-SAL0
       YDX=4.0d0*SAL1+2.0d0*SAL0
       if(YDX.gt.1.0d-9)then
        P1=YDY/YDX
       else
        P1=0.0d0
       endif
      
c      calculate degree of circularity
       down=        1.0d0*SAL0+ 7.0d0*SAL1
       doc=5.0d0*(  1.0d0*SAL0- 1.0d0*SAL1)
       if(down.ne.0.0d0)doc=doc/down

c      calculate depolarization
       down= 6.0d0*SAL0-12.0d0*SAL1
       ro  =-3.0d0*SAL0+ 9.0d0*SAL1
       if(down.ne.0.0d0)ro=ro/down
c
       if(units.eq.0.or.units.eq.1)then
        AMU=1822.0d0
        BOHR=0.529177d0
        SpA3=dsqrt(SAL0)/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)
c       gpisvejc=980.0
        gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/EXCA
c
c       Barron's tensor invariants:
c       alpha G = [(1/3)Sp alpha] [(1/3)Sp G'] etc
c       beta(A)^2=Delta2 in Gaussian:
        alphag=SAG0/9.0d0                *gpisvejc
        betaa=SA1*3.0d0/2.0d0            *gpisvejc
        betaalpha=0.5d0*(3.0d0*SAL1-SAL0)*gpisvejc
        betag=(3.0d0*SAG1-SAG0)/2.0d0    *gpisvejc
        alpha2=SAL0/9.0d0                *gpisvejc
c
c       DCP 180
c       ram3=24.0d0*beta2
        roa3=8.0d0*(12.0d0*betag+4.0d0*betaa)
c
c       polarized ROA,90deg:
        roa90x=45.0d0*alphag+7.0d0*betag+betaa
        ram90x=0.5d0*(45.0d0*alpha2+7.0d0*betaalpha)
c       depolarized ROA,90deg:
        roa90z=6.0d0*(betag-betaa/3.0d0)
        ram90z=0.5d0*6.0d0*betaalpha
c
c       unknown:
        S90X=S90X*(AMU*BOHR**4)
        D90X=D90x*gpisvejc
        DX=D90X
        IF(S90X.gt.1.0d-9)DX=DX/S90X
c
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
c       ICP(180):
c       8/c*(12*beta(G')^2 + 4*beta(A)^2):
        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
        doc=(3*P1-1)/(P1+1)
c       Gaussian units:
        if(IKIND.eq.1)then
c        replace 180SCP by 90 polarized:
         ram1=ram90x
         roa1=roa90x
        endif
        if(IKIND.eq.2)then
c        replace 180SCP by 90 depolarized:
         ram1=ram90z
         roa1=roa90z
        endif
        if(units.eq.1)then
c        Josef's units:
c        1.333d-21: cm^4 J^-1 Angstrom^-4 Amu
         pisvejc2=1.333d-21*(1.0d8/EXCA-ECM)**3/ECM
     1   /(1.0d0-exp(-(1.44d0*ECM/temp)))*
     1   ((ri**2+2.0d0)/3.0d0)**4
         YDX=YDX*pisvejc2*1.0d12
         YDY=YDY*pisvejc2*1.0d12
         ram1=ram1*pisvejc2*1.0d12
         roa1=roa1*pisvejc2*1.0d16
         roa3=roa3*pisvejc2*1.0d16
        endif
        WRITE(4,3000)IQ,ECM,YDX,YDY,ram1,CID2,DX,CID1,P1,roa1,doc,roa3
3000    FORMAT(I5,f9.2,6g12.3,f6.3,4g12.4)
       else
        PISVEJC=1500.0d0/7.0d0
        S180=PISVEJC*S180
        YDX =PISVEJC*YDX
        YDY =PISVEJC*YDY
c
        DZ   =1000.0d0*DZ
        DX   =1000.0d0*DX
        D    =1000.0d0*D 
        D180 =1000.0d0*D180
        SDCP180=SDCP180*1000.0d0
        doc=(3*P1-1)/(P1+1)
        WRITE(4,3000)IQ,ECM,YDX,YDY,S180,DZ,DX,D,P1,D180,doc,SDCP180
       endif
c
      endif
      endif
      ENDIF
1     CONTINUE
      WRITE(4,3002)
      CLOSE(4)
      RETURN
      END

      SUBROUTINE TRTEN(IR,M,np,S,N,ALPHA,ALPHAV,A,G,GV,ALPHA1,
     1ALPHA1V,G1,G1V,A1)
c     get all normal mode derivatives for frequency IR
      IMPLICIT REAL*8 (A-H,O-Z)
      implicit integer*4 (i-n)
      DIMENSION S(N,N),A(np,2*N,3,3,3),ALPHA(np,2*N,3,3),G(np,2*N,3,3)
      DIMENSION ALPHAV(np,2*N,3,3),GV(np,2*N,3,3)
      DIMENSION A1(2*N,3,3,3),ALPHA1(2*N,3,3),G1(2*N,3,3)
      DIMENSION ALPHA1V(2*N,3,3),G1V(2*N,3,3)
      
      A1=0.0d0
      ALPHA1=0.0d0
      ALPHA1V=0.0d0
      G1=0.0d0
      G1V=0.0d0
      DO 1 I=1,3
      DO 1 J=1,3
      DO 1 IQ=1,M
      DO 1 IX=1,N
      A1(IQ,I,J,1)=A1(IQ,I,J,1)+S(IX,IQ)*A(IR,IX,I,J,1)
      A1(IQ,I,J,2)=A1(IQ,I,J,2)+S(IX,IQ)*A(IR,IX,I,J,2)
      A1(IQ,I,J,3)=A1(IQ,I,J,3)+S(IX,IQ)*A(IR,IX,I,J,3)
      ALPHA1( IQ,I,J)=ALPHA1( IQ,I,J)+S(IX,IQ)*ALPHA( IR,IX,I,J)
      ALPHA1V(IQ,I,J)=ALPHA1V(IQ,I,J)+S(IX,IQ)*ALPHAV(IR,IX,I,J)
      G1(     IQ,I,J)=G1    ( IQ,I,J)+S(IX,IQ)*G(     IR,IX,I,J)
1     G1V(    IQ,I,J)=G1V   ( IQ,I,J)+S(IX,IQ)*GV(    IR,IX,I,J)
      DO 2 I=1,3
      DO 2 J=1,3
      DO 2 IQ=N+1,N+M
      DO 2 IX=N+1,2*N
      SI=S(IX-N,IQ-N)
      A1(IQ,I,J,1)=A1(IQ,I,J,1)+SI*A(IR,IX,I,J,1)
      A1(IQ,I,J,2)=A1(IQ,I,J,2)+SI*A(IR,IX,I,J,2)
      A1(IQ,I,J,3)=A1(IQ,I,J,3)+SI*A(IR,IX,I,J,3)
      ALPHA1( IQ,I,J)=ALPHA1( IQ,I,J)+SI*ALPHA( IR,IX,I,J)
      ALPHA1V(IQ,I,J)=ALPHA1V(IQ,I,J)+SI*ALPHAV(IR,IX,I,J)
      G1(     IQ,I,J)=G1    ( IQ,I,J)+SI*G(     IR,IX,I,J)
2     G1V(    IQ,I,J)=G1V   ( IQ,I,J)+SI*GV(    IR,IX,I,J)
      RETURN
      END
c     ============================================================
c     ============================================================
c     ====================
c     ====================

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

      subroutine  getatw(NAT,np)
      integer*4 NAT,np,i,j
      real*8 x
      OPEN(2,FILE='FILEW.TTT',STATUS='OLD')
      READ(2,*)
      READ(2,*)NAT
      rewind 2
      i=0
1     do 10 j=1,5
10    read(2,*,end=99,err=99)
      read(2,*,end=99,err=99)x,x,x,x,x
      do 11 j=1,9*NAT+3-2+2+9*NAT+3+2+27*NAT+9+2+9*NAT+3
11    read(2,*,end=99,err=99)
      i=i+1
      goto 1
99    np=i
      close(2)
      write(6,*)NAT,' atoms ',np,' frequencies'
      if(NAT.eq.0.or.np.eq.0)call report('nothing to do')
      return
      end

      SUBROUTINE  TRTENii(IQ,np,S,N,ALPHA,ALPHAV,A,G,GV,alphaw,alphavw,
     1GPw,GPvw,Aw)
c     for selected mode IQ transform all frequencies np
      IMPLICIT REAL*8 (A-H,O-Z)
      implicit integer*4 (i-n)
      DIMENSION S(N,N),A(np,2*N,3,3,3),ALPHA(np,2*N,3,3),G(np,2*N,3,3),
     1alphaw(2*np,3,3),GPw(2*np,3,3),Aw(2*np,3,3,3),GV(np,2*N,3,3),
     1alphavw(2*np,3,3),GPvw(2*np,3,3),ALPHAV(np,2*N,3,3)
      alphaw=0.0d0
      alphavw=0.0d0
      GPw=0.0d0
      GPvw=0.0d0
      Aw=0.0d0
      DO 1 I=1,3
      DO 1 J=1,3
      DO 1 IX=1,N
      DO 1 IW=1,np
      SI=S(IX,IQ)
      Aw(IW,I,J,1)=Aw(IW,I,J,1)+SI*A(IW,IX,I,J,1)
      Aw(IW,I,J,2)=Aw(IW,I,J,2)+SI*A(IW,IX,I,J,2)
      Aw(IW,I,J,3)=Aw(IW,I,J,3)+SI*A(IW,IX,I,J,3)
      ALPHAw( IW,I,J)=ALPHAw( IW,I,J)+SI*ALPHA( IW,IX,I,J)
      ALPHAvw(IW,I,J)=ALPHAvw(IW,I,J)+SI*ALPHAV(IW,IX,I,J)
      GPw(    IW,I,J)=GPw(    IW,I,J)+SI*G(     IW,IX,I,J)
1     GPvw(   IW,I,J)=GPvw(   IW,I,J)+SI*GV(    IW,IX,I,J)
      DO 2 I=1,3
      DO 2 J=1,3
      DO 2 IX=N+1,2*N
      DO 2 IW=np+1,2*np
      SI=S(IX-N,IQ)
      IO=IW-np
      Aw(IW,I,J,1)=Aw(IW,I,J,1)+SI*A(IO,IX,I,J,1)
      Aw(IW,I,J,2)=Aw(IW,I,J,2)+SI*A(IO,IX,I,J,2)
      Aw(IW,I,J,3)=Aw(IW,I,J,3)+SI*A(IO,IX,I,J,3)
      ALPHAw( IW,I,J)=ALPHAw( IW,I,J)+SI*ALPHA( IO,IX,I,J)
      ALPHAvw(IW,I,J)=ALPHAvw(IW,I,J)+SI*ALPHAV(IO,IX,I,J)
      GPw(    IW,I,J)=GPw(    IW,I,J)+SI*G(     IO,IX,I,J)
2     GPvw(   IW,I,J)=GPvw(   IW,I,J)+SI*GV(    IO,IX,I,J)
      RETURN
      END

      subroutine wrw(IQ,v,n,s,w)
      implicit none
      real*8 v(*),w(*)
      integer*4 n,i,IQ,is
      character*(*) s
      character*10 s10
      write(s10,10)IQ
10    format(i10)
      do 2 is=1,len(s10)
2     if(s10(is:is).ne.' ')goto 3
3     open(12,file=s//'.'//s10(is:len(s10))//'.prn')
      do 1 i=1,n
1     write(12,120)w(i)/8065.54093D0,v(i)
120   format(f12.6,' ',g12.4)
      close(12)
      write(6,*)s//'.'//s10(is:len(s10))//'.prn written'
      return
      end

      subroutine vz(x,n)
      implicit none
      integer*4 n,i
      real*8 x(*)
      do 1 i=1,n
1     x(i)=0.0d0
      return
      end