PROGRAM NEW6
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA0(3,3),A0(3,3,3),G0(3,3),ALPHA1(3,3),
     3A1(3,3,3),G1(3,3),alpha0v(3,3),alpha1v(3,3)
      real*8,allocatable::FCAR(:,:),S(:,:),ALPHA(:,:,:),A(:,:,:,:),
     1G(:,:,:),R(:,:),alphav(:,:,:),tem(:,:),E(:),
     1Gi(:,:,:),ALPHAi(:,:,:),Ai(:,:,:,:)
      CHARACTER*10 ATSTR
      CHARACTER*15 S15
      CHARACTER*50 STRING
      CHARACTER*16 S16
      EQUIVALENCE(S15,STRING)
      LOGICAL SCALED,RAMAN,TTT,XYZ,LPOLAR,lpro,lgpro,lapro,
     1FILEQ,luseA,luseG,lor,lintA,lintG,lcompl
      integer*4 units
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(SCALED,TTT,RAMAN,XYZ,N1,N2,LPOLAR,EXCA,
     1wmin,wmax,lpro,lgpro,lapro,N6,FILEQ,N7,IDIF,units,temp,ri,IKIND,
     1luseA,luseG,lor,lintA,lintG,lcompl)

c     read transformed tensors from FILE.Q.TTT:
      if(FILEQ)then
       NM=0
       CALL READTTTQ(NM,ALPHA,A,G,E,0)
       allocate(ALPHA(NM,3,3),A(NM,3,3,3),G(NM,3,3),E(NM))
       CALL READTTTQ(NM,ALPHA,A,G,E,1)
       goto 5005
      endif

      IF(TTT)then
       OPEN(2,FILE='FILE.TTT',STATUS='OLD')
       READ(2,*)
       READ(2,*)NAT
       close(2)
       N=3*NAT
       allocate(FCAR(N,N),S(N,N),ALPHA(N,3,3),A(N,3,3,3),G(N,3,3),
     1 R(3,NAT),alphav(N,3,3),tem(3,NAT),E(N))
       CALL READTEN(N,NAT,ALPHA,G,A,0)
       if(lcompl)then
        allocate(ALPHAi(N,3,3),Ai(N,3,3,3),Gi(N,3,3))
        CALL READTEN(N,NAT,ALPHAi,Gi,Ai,3)
        WRITE(6,*)'(Complex part)'
       endif
      else
       NAT=100
       N=3*NAT
       allocate(FCAR(N,N),S(N,N),ALPHA(N,3,3),A(N,3,3,3),G(N,3,3),
     1 R(3,NAT),alphav(N,3,3),tem(3,NAT),E(N))
      endif
      N30=N
      IF(XYZ)CALL READXYZ(R,NAT,N,1)
      IF(TTT.AND.XYZ.AND.SCALED) GOTO 5003
      
C
      IGEO=0
      OPEN(2,FILE='C5R.OUT')
  201 READ(2,1111,END=999,ERR=999)STRING
1111  FORMAT(A50)
C                         12345678901234
      IF (STRING(2:14).EQ.'Geometry, (in')THEN
       WRITE(3,*) ' Geometry found ...'
       WRITE(3,*)
       OPEN(4,FILE='FILE.XYZ')
       WRITE(4,441)
441    FORMAT(' MOLECULE STUDIED FOR ROA')
       DO 35 I=1,5
35     READ(2,*)
       IA=0
36     READ(2,1112)ATSTR
1112   FORMAT(A10)
       IF(ATSTR.NE.' ---------')THEN
        BACKSPACE 2
        READ(2,1113)S15,AX,X,Y,Z
        IA=IA+1
        R(1,IA)=X
        R(2,IA)=Y
        R(3,IA)=Z
        BOHR=0.529177d0
        I=-1
C                  123456789012345
        IF(S15.EQ.'    HYDROGEN   ')I=4
        IF(S15.EQ.'    CARBON     ')I=1
        IF(S15.EQ.'    OXYGEN     ')I=2
        IF(S15.EQ.'    NITROGEN   ')I=3
        IF(I.EQ.-1)WRITE(3,*)' UNKNOWN ATOM', S15
        WRITE(4,442)I,X*BOHR,Y*BOHR,Z*BOHR
442     FORMAT(I4,3F15.6,' 0 0 0 0 0.0')
1113   FORMAT(A15,F4.1,3F15.8,I5)
        WRITE(3,1113)S15,AX,X,Y,Z,IA
        GOTO 36
       ENDIF
       NAT=IA
       WRITE(4,442)9999,0.0,0.0,0.0
       CLOSE(4)
       IF(NAT.EQ.2)THEN
        N1=6
        N7=6
       ENDIF
       N=3*NAT
       IGEO=IGEO+1
      ENDIF
      IF (STRING(32:49).EQ.'Molecular Geometry')THEN
       WRITE(3,*) ' Geometry found ...'
       WRITE(3,*)
       OPEN(4,FILE='FILE.XYZ')
       WRITE(4,441)
       DO 351 I=1,2
351    READ(2,*)
       IA=0
361    READ(2,1112)ATSTR
       IF(ATSTR.NE.'          ')THEN
        BACKSPACE 2
        READ(2,11131)S16,X,Y,Z
        IA=IA+1
        R(1,IA)=X
        R(2,IA)=Y
        R(3,IA)=Z
        BOHR=0.529177d0
        I=-1
C                  1234567890123456
        IF(S16(7:12).EQ.'HYDROG')I=4
        IF(S16(7:12).EQ.'CARBON')I=1
        IF(S16(7:12).EQ.'OXYGEN')I=2
        IF(S16(7:12).EQ.'NITROG')I=3
        IF(I.EQ.-1)WRITE(3,*)' UNKNOWN ATOM', S15
        WRITE(4,442)I,X*BOHR,Y*BOHR,Z*BOHR
11131 FORMAT(A16,3F20.10)
        WRITE(3,11131)S16,X,Y,Z
        GOTO 361
       ENDIF
       NAT=IA
       WRITE(4,442)9999,0.0,0.0,0.0
       CLOSE(4)
       IF(NAT.EQ.2)THEN
        N1=6
        N7=6
       ENDIF
       N=3*NAT
       IGEO=IGEO+1
      ENDIF
      IF(IGEO.GT.0.AND.TTT.AND.SCALED)GOTO 5003
C      
c     Cadpac output:
      IF (STRING(2:15).EQ.'#RUNTYPE = CUB'.AND.(.NOT.TTT))THEN
       WRITE(3,*)
       WRITE(3,*) ' Cubic RUNTYPE found ...'
       DO 38 I=1,10
       IF(I.EQ.6)THEN
        READ(2,11141)STEP
11141   FORMAT(16X,F7.4)
        WRITE(3,11142)STEP
11142   FORMAT(1X,'  STEP OF THE DIFFERENTIATION : ',F7.4,' BOHR')
        GOTO 38
       ENDIF
       IF(I.EQ.5)THEN
        READ(2,11143)INUM
        WRITE(3,11143)INUM
11143 FORMAT(16X,I4)
        IF(INUM-2*(INUM/2).EQ.0)INUM=2
        IF(INUM.NE.2)INUM=1
        IF(INUM.EQ.1)WRITE(3,11144)
        IF(INUM.EQ.2)WRITE(3,21144)
11144 FORMAT(1X,'  ONE - POINT DIFFERENTIATION  ( X  X+H )   ')
21144 FORMAT(1X,'  TWO - POINT DIFFERENTIATION  ( X-H X X+H )')
        GOTO 38
       ENDIF
       READ(2,*)
38     CONTINUE
40     READ(2,1111)STRING
C                         123456789012345
       IF(STRING(2:15).EQ.'^^^^^^^^^^^^^^')THEN
        IDIF=IDIF+1
        READ(2,1115)IF
1115  FORMAT(I5)
        WRITE(3,1116)IF,IDIF
1116  FORMAT(I3,' (',I3,'). Point of the numerical differentiation.')
        ILAST=0
        IF(IDIF.EQ.3*NAT)ILAST=1
        IF(IDIF.EQ.2*NAT*3)ILAST=2
        IF(IDIF.EQ.0)THEN
         WRITE(3,*)' ZERO POINT TENSORS:'
         CALL READALPHA(ALPHA0)
         WRITE(3,*)' ALPHA'
         WRITE(3,5002)((ALPHA0(I,J),I=1,3),J=1,3)
5002  FORMAT(3(3F15.6,/),/)
         CALL READA(A0)
         WRITE(3,*)' A'
         WRITE(3,5002)(((A0(I,J,K),I=1,3),J=1,3),K=1,3)
         CALL READG(G0)
         WRITE(3,*)' G'
         WRITE(3,5002)((G0(I,J),J=1,3),I=1,3)
         CALL READALV(alpha0v)
         WRITE(3,*)' ALPHAV'
         WRITE(3,5002)((alpha0v(I,J),I=1,3),J=1,3)
         GOTO 40
        ENDIF
        CALL READALPHA(ALPHA1)
        CALL READA(A1)
        CALL READG(G1)
        CALL READALV(alpha1v)
        IC=IDIF
        fs=1.0d0/STEP
        IF(IDIF.LE.3*NAT)THEN
         DO 41 I=1,3
         DO 41 J=1,3
         ALPHA(IDIF,I,J)=(ALPHA1(I,J)-ALPHA0(I,J))*fs
         alphav(IDIF,I,J)=(alpha1v(I,J)-alpha0v(I,J))*fs
         G(IDIF,I,J)=(G1(I,J)-G0(I,J))*fs
         DO 41 K=1,3
41       A(IDIF,I,J,K)=(A1(I,J,K)-A0(I,J,K))*fs
        ELSE
         IC=IDIF-3*NAT
         DO 412 I=1,3
         DO 412 J=1,3
         ALPHA(IC,I,J)=0.5d0*(ALPHA(IC,I,J)-
     1                      (ALPHA1(I,J)-ALPHA0(I,J))*fs)
         alphav(IC,I,J)=0.5d0*(alphav(IC,I,J)-
     1                      (alpha1v(I,J)-alpha0v(I,J))*fs)
         G(IC,I,J)    =0.5d0*(    G(IC,I,J)-
     1                      (    G1(I,J)-      G0(I,J))*fs)
         DO 412 K=1,3
412      A(IC,I,J,K)=
     1        0.5d0*( A(IC,I,J,K) - (A1(I,J,K)-A0(I,J,K)) *fs )
        ENDIF
        IF(ILAST.GT.0)THEN
         IF(INUM.EQ.2.AND.ILAST.EQ.1)GOTO 40
         WRITE(3,*)' End of the numerical differentiation.'
         open(22,file='ROAR.TTT')
         WRITE(22,2000)NAT
         CALL WRITETTT(N,NAT,ALPHA,A,G,alphav)
         write(6,*)'ROAR.TTT written'
         close(22)
         IF(SCALED)GOTO 5003
         GOTO 201
        ENDIF
       ENDIF
       GOTO 40
      ENDIF
C                 1234567890123456789012345678901234567890
      IF((S15.EQ.' RUN TYPE = SEC'.OR.S15.EQ.
     1           ' RUN TYPE = FOR').AND.(.NOT.SCALED))THEN
       CALL READSEC(N,FCAR,3*NAT,1)
       CALL READSMAT(S,3*NAT,E)
       GOTO 5001
      ENDIF
      GOTO 201
C
5003  CALL READS(N,S,E,N7)
      GOTO 5004
5001  CALL CONTROL(S,FCAR,3*NAT,E,N7)
      IF(RAMAN)THEN
       WRITE(3,*)' ROA tensors zeroed out ...'
       DO 1 I=1,3
       DO 1 J=1,3
       G0(I,J)=0.0d0
       G1(I,J)=0.0d0
       DO 2 K=1,3
       A0(I,J,K)=0.0d0
2      A1(I,J,K)=0.0d0
       ALPHA0(I,J)=0.0d0
       ALPHA1(I,J)=0.0d0
       DO 1 K=1,N
       ALPHA(K,I,J)=0.0d0
       G(K,I,J)=0.0d0
       DO 1 L=1,3
1      A(K,I,J,L)=0.0d0
      ENDIF
5004  IF(.NOT.TTT)CALL READAAL(3*NAT,ALPHA)
      IF(.NOT.(TTT.AND.SCALED))CLOSE(2)
      IF(NAT*3.GT.N30)THEN
       WRITE(3,*)' Too many atoms, maximum ',N30
       STOP
      ENDIF
C
      IF(LPOLAR)THEN
c      set A, G to zero:
       CALL TTTT(N,ALPHA,A,G,NAT,3,R,tem)
c      make common A, G from alpha:
       CALL TTTT(N,ALPHA,A,G,NAT,2,R,tem)
       write(6,*)' A, G made from alpha !!!'
      ENDIF

      if(.not.lintA.or..not.lintG)then
c      set internal A, G to zero:
c      1) transform A, G to atomic origins:
       CALL TTTT(N,ALPHA,A,G,NAT,1,R,tem)
c      2) delete them
       if(.not.lintG)CALL TTTT(N,ALPHA,A,G,NAT,4,R,tem)
       if(.not.lintA)CALL TTTT(N,ALPHA,A,G,NAT,5,R,tem)
c      3) transform A, G back to common origin:
       CALL TTTT(N,ALPHA,A,G,NAT,2,R,tem)
      endif
c
      call rules(ALPHA,A,G,NAT)
c     CALL TTTT(N,ALPHA,A,G,NAT,1,R,tem)
      CALL PROJECT(ALPHA,A,G,NAT,R,lpro,lgpro,lapro,N6)
      call rules(ALPHA,A,G,NAT)
c     CALL TTTT(N,ALPHA,A,G,NAT,2,R,tem)
c
C     write used cartesian tensor derivatives into ROA.TTT
      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')
      CALL WRITETTT(N,NAT,ALPHA,A,G,alphav)
      close(22)
C
C     transform the tensors into internal coordinate derivatives
      CALL TRTEN(S,3*NAT,ALPHA,A,G)
      if(lcompl)then
       CALL TRTEN(S,3*NAT,ALPHAi,Ai,Gi)
       CALL WTQi(N,N7,3*NAT,ALPHA,A,G,ALPHAi,Ai,Gi,alphav,E,WMIN,WMAX)
      else
C      write used normal mode tensor derivatives into FILE.Q.TTT
       CALL WRITETTTQ(N,N7,3*NAT,ALPHA,A,G,alphav,E,WMIN,WMAX)
      endif
c

      NM=3*NAT

5005  continue

      if(lcompl)then
       CALL DORAMANi(NM,ALPHA,ALPHAi,G,Gi,A,Ai,EXCA,wmin,wmax,E,
     1 units,temp,ri,IKIND,luseA,luseG,lor)
      else
       CALL DORAMAN(NM,ALPHA,G,A,EXCA,wmin,wmax,E,
     1 units,temp,ri,IKIND,luseA,luseG,lor)
      endif
      CALL DORAYLEIGH
      WRITE(6,*)' PROGRAM FINISHED OK'
      CLOSE(3)
      STOP
999   close(2)
      stop
      END
C     ==========================================
      SUBROUTINE READTTTQ(N,ALPHA,A,G,E,ic)
      IMPLICIT none
      integer*4 N,ic,I,J,K,II
      real*8 ALPHA(N,3,3),A(N,3,3,3),G(N,3,3),CM,E(*)
c
      CM=219470.0d0
      open(22,file='FILE.Q.TTT',status='old')
      read(22,*)
      read(22,*)N
      if(ic.eq.1)then
       read(22,*)
       read(22,*)
       do 1 II=1,N
       read(22,*)E(II),E(II)
       E(II)=E(II)/CM
       DO 1  I=1,3
1      read(22,*)ALPHA(II,I,1),(ALPHA(II,I,J),J=1,3)

       read(22,*)
       read(22,*)
       DO 2 II=1,N
       read(22,*)
       DO 2 I=1,3
2      read(22,*)G(II,I,1),(G(II,I,J),J=1,3)
c      last inner index magnetic

       read(22,*)
       read(22,*)
       DO 3 II=1,N
       read(22,*)
       DO 3 I=1,3
       DO 3 J=1,3
3      read(22,*)(A(II,I,J,K),K=1,2),(A(II,I,J,K),K=1,3)

      write(6,*)' FILE.Q.TTT read, ',N,' modes'
      endif

      close(22)
      RETURN
      END
C     ==========================================
      SUBROUTINE WTQi(N,N7,NINT,ALPHA,A,G,ALPHAi,Ai,Gi,av,E,WMIN,WMAX)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),
     1ALPHAi(N,3,3),Gi(N,3,3),Ai(N,3,3,3),av(N,3,3),E(*)
c
      CM=219470.0d0
      im=0
      DO 4 II=N7,NINT
4     if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)im=im+1
      open(22,file='FILE.Q.TTT')
      WRITE(22,2003)im
2003  FORMAT(' ROA tensors, normal modes derivatives',/,I4,' modes',/,
     1' The electric-dipolar electric-dipolar polarizability:',/,
     2        ' mode e(cm-1)   jx           jy           jz')
      im=0
      DO 1 II=N7,NINT
      if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)then
       im=im+1
       write(22,221)im,E(II)*CM
221    format(i5,f12.2)
       DO 11 I=1,3
11     WRITE(22,2001)I,(ALPHA(II,I,J),J=1,3),(ALPHAi(II,I,J),J=1,3)
2001   FORMAT(I3,6F13.7)
      endif
1     continue

      WRITE(22,2004)
2004  FORMAT(' The electric dipole magnetic dipole polarizability:',/,
     1      ' mode  e(cm-1)  jx(Bx)       jy(By)       jz(Bz)')
      im=0
      DO 2 II=N7,NINT
      if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)then
       im=im+1
       write(22,221)im,E(II)*CM
       DO 21 I=1,3
21     WRITE(22,2001)I,(G(II,I,J),J=1,3),(Gi(II,I,J),J=1,3)
c      last inner index magnetic
      endif
2     continue

      WRITE(22,2005)
2005  FORMAT(' The electric dipole electric quadrupole polarizability:',
     2/,     ' mode  e(cm-1)  kx           ky           kz')
      im=0
      DO 3 II=N7,NINT
      if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)then
       im=im+1
       write(22,221)im,E(II)*CM
       DO 31 I=1,3
       DO 31 J=1,3
31     WRITE(22,2002)I,J,(A(II,I,J,K),K=1,3),(Ai(II,I,J,K),K=1,3)
2002   FORMAT(2I3,6F13.7)
      endif
3     continue

      WRITE(22,2009)
2009  FORMAT('The velocity form of alpha',/,
     1' mode e(cm-1)    vx           vy           vz')
      im=0
      DO 5 II=N7,NINT
      if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)then
       im=im+1
       write(22,221)im,E(II)*CM
       DO 51 I=1,3
51     WRITE(22,2001)I,(av(II,I,J),J=1,3),(av(II,I,J),J=1,3)
      endif
5     continue

      close(22)
      RETURN
      END
C     ==========================================
      SUBROUTINE WRITETTTQ(N,N7,NINT,ALPHA,A,G,alphav,E,WMIN,WMAX)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),
     1alphav(N,3,3),E(*)
c
      CM=219470.0d0
      im=0
      DO 4 II=N7,NINT
4     if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)im=im+1
      open(22,file='FILE.Q.TTT')
      WRITE(22,2003)im
2003  FORMAT(' ROA tensors, normal modes derivatives',/,I4,' modes',/,
     1' The electric-dipolar electric-dipolar polarizability:',/,
     2        ' mode e(cm-1)   jx           jy           jz')
      im=0
      DO 1 II=N7,NINT
      if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)then
       im=im+1
       write(22,221)im,E(II)*CM
221    format(i5,f12.2)
       DO 11 I=1,3
11     WRITE(22,2001)I,(ALPHA(II,I,J),J=1,3)
2001   FORMAT(I3,3F13.7)
      endif
1     continue

      WRITE(22,2004)
2004  FORMAT(' The electric dipole magnetic dipole polarizability:',/,
     1      ' mode  e(cm-1)  jx(Bx)       jy(By)       jz(Bz)')
      im=0
      DO 2 II=N7,NINT
      if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)then
       im=im+1
       write(22,221)im,E(II)*CM
       DO 21 I=1,3
21     WRITE(22,2001)I,(G(II,I,J),J=1,3)
c      last inner index magnetic
      endif
2     continue

      WRITE(22,2005)
2005  FORMAT(' The electric dipole electric quadrupole polarizability:',
     2/,     ' mode  e(cm-1)  kx           ky           kz')
      im=0
      DO 3 II=N7,NINT
      if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)then
       im=im+1
       write(22,221)im,E(II)*CM
       DO 31 I=1,3
       DO 31 J=1,3
31     WRITE(22,2002)I,J,(A(II,I,J,K),K=1,3)
2002   FORMAT(2I3,3F13.7)
      endif
3     continue

      WRITE(22,2009)
2009  FORMAT('The velocity form of alpha',/,
     1' mode e(cm-1)    vx           vy           vz')
      im=0
      DO 5 II=N7,NINT
      if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)then
       im=im+1
       write(22,221)im,E(II)*CM
       DO 51 I=1,3
51     WRITE(22,2001)I,(ALPHAV(II,I,J),J=1,3)
      endif
5     continue

      close(22)
      RETURN
      END
C     ==========================================
      SUBROUTINE WRITETTT(N,NAT,ALPHA,A,G,alphav)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),
     1alphav(N,3,3)
c
      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(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(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(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(II,I,J),J=1,3)
      RETURN
      END

C     ================================================================
      SUBROUTINE READAAL(N,ALPHA)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(N,3,3),ALPHAA(N,3,3)
      CHARACTER*25 STRING
201   READ(2,1111,END=1)STRING
1111  FORMAT(A25)
C                    1234567890123456789012345
      IF (STRING.EQ.' Cartesian Polarisability') THEN
       WRITE(3,*)' Cartesian Polarisability Derivatives found ... '
       READ(2,*)
       WRITE(3,3001)
3001  FORMAT(1X,'  POLARIZABILITY DERIVATIVES',/,
     1       1X,'            ANALYTICAL     NUMERICAL',/,
     2       1X,'-------------------------------------')
       DO 203 I=1,3
       DO 203 J=1,3
       DO 203 K=1,N
       READ(2,3002)ALPHAA(K,I,J)
3002  FORMAT(30X,F18.8)
       WRITE(3,3000)K,I,J,ALPHAA(K,I,J),ALPHA(K,I,J)
203    ALPHA(K,I,J)=ALPHAA(K,I,J)
3000  FORMAT(1X,3I3,2F15.8)
       WRITE(3,3003)
3003  FORMAT(1X,'THE ANALYTICAL FORM OF ALPHA WILL BE USED')
       RETURN
      ENDIF
      GOTO  201
1     WRITE(3,*)' Analytical alpha not found, numerical ',
     1'derivatives will be used ...'
      RETURN
      END
C     ==============================================================
      SUBROUTINE READALPHA(ALPHA)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(3,3)
      CHARACTER*20 STRING
      IC=0
201   READ(2,1111)STRING
      IC=IC+1
1111  FORMAT(A20)
C                    12345678901234567890
      IF (STRING.EQ.' Polarizability tens') THEN
       WRITE(3,*)' Alpha tensor found ... '
       DO 202 I=1,3
202    READ(2,*)
       DO 203 I=1,3
203    READ(2,*)(ALPHA(I,J),J=1,3)
       RETURN
      ENDIF
      IF(IC.GT.20000)THEN
       WRITE(2,*)' Alpha not found, program stopped ... '
       STOP
      ENDIF
      GOTO  201
      END
C     ==============================================================
      SUBROUTINE READA(A)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(3,3,3)
      CHARACTER*20 STRING
      IC=0
201   READ(2,1111)STRING
      IC=IC+1
1111  FORMAT(A20)
C                    12345678901234567890
      IF (STRING.EQ.' Dipole-Quadrupole P') THEN
       WRITE(3,*)' A tensor found ... '
       READ(2,*)
       DO 203 I=1,3
c      dipole slowest (outer) first index
203    READ(2,*)A(I,1,1),A(I,2,2),A(I,3,3),
     1          A(I,1,2),A(I,1,3),A(I,2,3)
       DO 204 I=1,3
       DO 204 J=1,3
       DO 204 K=J+1,3
204    A(I,K,J)=A(I,J,K)
       RETURN
      ENDIF
      IF(IC.GT.20000)THEN
       WRITE(2,*)' A not found, program stopped ... '
       STOP
      ENDIF
      GOTO  201
      END
C     ==============================================================
      SUBROUTINE READG(G)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION G(3,3)
      CHARACTER*20 STRING
      IC=0
201   READ(2,1111)STRING
      IC=IC+1
1111  FORMAT(A20)
C                    12345678901234567890
      IF (STRING.EQ.' LOW FREQUENCY LIMIT') THEN
       WRITE(3,3000)
3000  FORMAT(22H  G' tensor found ... )
       DO 202 I=1,3
202    READ(2,*)
       DO 203 I=1,3
       READ(2,*)
c      last inner index magnetic
 203   READ(2,3002)(G(I,J),J=1,3)
3002   FORMAT(12X,3F20.13)
       RETURN
      ENDIF
      IF(IC.GT.20000)THEN
       WRITE(2,3001)
3001  FORMAT(36H  G' not found, program stopped ... )
       STOP
      ENDIF
      GOTO  201
      END
C     ===============================================================
      SUBROUTINE READSEC(N0,FCAR,N,IP)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(N0,N)
      CHARACTER*15 STRING
      IC=0
201   READ(2,1111)STRING
      IC=IC+1
1111  FORMAT(A15)
      IF (STRING.EQ.' Total cartesia'.OR.
     1    STRING.EQ.' Cartesian Forc')THEN
       WRITE(3,*)' Second derivatives found ... '
       DO 202 I=1,7
202    READ(2,*)
       N0C=0
204    NC=6*N0C+1
       MC=NC+5
       IF (MC.GT.N)MC=N
       IF(IP.EQ.0)THEN
        DO 2030 IC=1,N
2030    READ(2,2700)(FCAR(IC,JC),JC=NC,MC)
       ELSE
        DO 203 IC=1,N
203     READ(2,270)(FCAR(IC,JC),JC=NC,MC)
       ENDIF
270    FORMAT(18X,6F10.5)
2700   FORMAT(18X,6F14.9)
       IF(MC.LT.N) THEN
        DO 205 I=1,7
205     READ(2,*)
        N0C=N0C+1
        GOTO 204
       ENDIF
       RETURN
      ENDIF
      IF(IC.GT.20000)THEN
       WRITE(2,*)' Seconds not found, program stopped ... '
       STOP
      ENDIF
      GOTO  201
      END
C     ==============================================================
      SUBROUTINE CONTROL(S,F,N,E,N7)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION S(N,N),F(N,N),E(*)
      real*8,allocatable::QL(:,:)
      allocate(QL(N,N))
      WRITE(3,*)' NORMALIZATION OF S-MATRIX'
C     SUMijk (SiqSjpFij) = OMp^2 p=q
      DO 1 I=1,N
      DO 1 J=1,N
      QL(I,J)=0.0d0
      DO 1 II=1,N
      DO 1 JJ=1,N
1     QL(I,J)=QL(I,J)+S(II,I)*S(JJ,J)*F(II,JJ)
      DO 3 J=1,N
3     WRITE(3,2)(QL(I,J),I=1,N)
2     FORMAT(11F7.4)
      SFAC=0.0d0
      DO 4 I=N7,N
4     SFAC=SFAC+E(I)*E(I)/QL(I,I)
      SFAC=SFAC/DBLE(N-N7+1)
      SFAC=SQRT(SFAC)
      DO 5 I=1,N
      DO 5 J=1,N
5     S(I,J)=S(I,J)*SFAC
      WRITE(3,6)SFAC
6     FORMAT(1X,'FACTOR FOR S-MATRIX: ',F20.10)
      NAT=N/3
      OPEN(4,FILE='S.MAT',FORM='UNFORMATTED')
      WRITE(4)NAT,(E(I),I=1,3*NAT),((S(I,J),J=1,3*NAT),I=1,3*NAT),SFAC
      CLOSE(4)
      WRITE(3,*)' S.MAT written on the disk'
      RETURN
      END
C     ==============================================================
      SUBROUTINE READSMAT(S,N,E)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION S(N,N),E(*)
      CHARACTER*15 STRING
      CM=219470.0d0
C     1 a.u.= 2.1947E5 cm^-1
      NAT=N/3
      IC=0
201   READ(2,1111)STRING
      IC=IC+1
1111  FORMAT(A15)
C                    1234567890123456789012345678901234567890
      IF (STRING.EQ.' Transformation')THEN
       WRITE(3,*)' S matrix found ... '
       DO 1 I=1,3
1      READ(2,*)
       imode=0
       i=0
800    CONTINUE
       NM=3*NAT-IMODE
       IF(NM.GT.6)NM=6
       READ(2,*)(E(IE),IE=IMODE+1,IMODE+NM)
       IMODE=IMODE+NM
       read(2,*)
       read(2,*)
       DO 2 J=1,NAT
       read(2,200)(S((J-1)*3+1,IMODE-NM+JJ),JJ=1,NM)
       read(2,200)(S((J-1)*3+2,IMODE-NM+JJ),JJ=1,NM)
       read(2,200)(S((J-1)*3+3,IMODE-NM+JJ),JJ=1,NM)
200    FORMAT(17x,6F10.5)
2      CONTINUE
       IF(Imode.LT.3*nat) then
        DO 3 i=1, 3
3       read(2,*)
        goto 800
       endIF
       DO 4 I=1,3*NAT
4      E(I)=E(I)/CM
       RETURN
      ENDIF
      IF (IC.GT.20000)THEN
       WRITE(3,*)' S - matrix not found ...'
       STOP
      ENDIF
      GOTO 201
      end
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,N7)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION S(N3,N3),E(*)
      CM=219470.0d0
      SFAC=0.0234280d0
C     1 a.u.= 2.1947E5 cm^-1
      WRITE(3,*)' S matrix read from the file F.INP:'
      IS=0
      OPEN(4,FILE='F.INP',STATUS='OLD',IOSTAT=IS)
      IF(IS.NE.0)THEN
       IF(IS.EQ.6)WRITE(3,*)' FILE "F.INP" NOT FOUND'
       WRITE(3,*)' PROGRAM ERROR STOP'
       STOP
      ENDIF
      read(4,*)nint,nat3,nat
      N=NAT3
      N7=N-NINT+1
      WRITE(3,3002)
3002  FORMAT(1X,'THE GEOMETRY (CHECK IF THE SAME AS THE AB INITIO):')
      do 1 i=1,NAT
      read(4,*)ki,xi1,xi2,xi3
      WRITE(3,3001)i,ki,xi1,xi2,xi3
3001  FORMAT(2I3,3F15.6)
1     CONTINUE
      read(4,*)
      DO 2 I=1,NAT
      DO 2 J=1,NINT
2     read(4,*)(s(3*(i-1)+ix,j+N7-1),ix=1,2),
     1(s(3*(i-1)+ix,j+N7-1),ix=1,3)
      read(4,*)nint
      READ(4,4000)(E(I),I=N7,N)
4000  FORMAT(6F11.6)
      DO 3 I=N7,N
3     E(I)=E(I)/CM
      DO 7 I=1,N7-1
      E(I)=0
      DO 7 J=1,N
7     S(J,I)=0.0d0
      DO 5 I=1,N
      DO 5 J=1,N
5     S(I,J)=S(I,J)*SFAC
      WRITE(3,6)SFAC
6     FORMAT(1X,'FACTOR FOR S-MATRIX: ',F20.10)
      NAT=N/3
      OPEN(5,FILE='S.MAT',FORM='UNFORMATTED')
      WRITE(5)NAT,(E(I),I=1,3*NAT),((S(I,J),J=1,3*NAT),I=1,3*NAT),SFAC
      CLOSE(5)
      WRITE(3,*)' S.MAT written on the disk'
      RETURN
      end
C
      SUBROUTINE READTEN(N,NAT,P,G,A,ic)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION P(N,3,3),G(N,3,3),A(N,3,3,3)
      OPEN(2,FILE='FILE.TTT',STATUS='OLD')
      READ(2,*)
      READ(2,*)NAT
      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
      M=3*(L-1)+IX
1     READ(2,*)(P(M,I,J),J=1,2),(P(M,I,J),J=1,ic),(P(M,I,J),J=1,3)
      READ(2,*)
      READ(2,*)
      DO 2 I=1,3
      READ(2,*)
      DO 2 L=1,NAT
      DO 2 IX=1,3
      M=3*(L-1)+IX
c     last index magnetic, inner
2     READ(2,*)(G(M,I,J),J=1,2),(G(M,I,J),J=1,ic),(G(M,I,J),J=1,3)
      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
      M=3*(L-1)+IX
3     READ(2,*)(A(M,I,J,K),K=1,2),(A(M,I,J,K),K=1,ic),
     1(A(M,I,J,K),K=1,3)
      CLOSE(2)
      WRITE(6,*)NAT,' atoms, tensors read in from FILE.TTT'
      RETURN
      END
C     ==============================================================
      SUBROUTINE READXYZ(R,NAT,N,ic)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION R(3,N/3)
      BOHR=0.529177d0
c
      OPEN(2,FILE='FILE.X',STATUS='OLD')
      READ(2,*)
      READ(2,*)nat
      if(ic.eq.0)then
       WRITE(6,*)' Number of atoms ',NAT,' read from FILE.X'
      else
       do 1 I=1,nat
       READ(2,*)R(1,I),(R(J,I),J=1,3)
       DO 1 J=1,3
1      R(J,I)=R(J,I)/BOHR
       CLOSE(2)
       WRITE(6,*)NAT,' atoms, geometry read in from FILE.X'
      endif
      RETURN
      END
C
      SUBROUTINE IOPAR(SCALED,TTT,RAMAN,XYZ,N1,N2,LPOLAR,EXCA,
     1wmin,wmax,lpro,lgpro,lapro,N6,FILEQ,N7,IDIF,units,temp,ri,IKIND,
     1luseA,luseG,lor,lintA,lintG,lcompl)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      LOGICAL SCALED,TTT,RAMAN,XYZ,LPOLAR,lpro,lex,lgpro,lapro,
     1FILEQ,luseA,luseG,lor,lintA,lintG,lcompl
      CHARACTER*15 STRING
      integer*4 units
c
      N7=7
      N1=7
      N2=0
      IDIF=-1
      wmin=0.0d0
      wmax=1.0d90
      SCALED=.TRUE.
      lpro=.false.
      lapro=.false.
      lcompl=.false.
      lgpro=.false.
      RAMAN=.FALSE.
      N6=3
      LPOLAR=.FALSE.
      TTT=.TRUE.
      XYZ=.TRUE.
c     read FILE.TTT.Q instead of FILE.TTT:
      FILEQ=.false.
c     excitation light wavelength in nm:
      EXCNM=532.0d0
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.
c     use intrinsic part of A-tensor derivatives:
      lintA=.true.
c     use intrinsic part of G'-tensor derivatives:
      lintG=.true.
c     generate also spectra for oriented sampels:
      lor=.false.

      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
        goto 99
       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.'KIND')READ(2,*)IKIND
      IF(STRING(1:6).EQ.'SCALED')READ(2,*)SCALED
      IF(STRING(1:5).EQ.'POLAR')READ(2,*)LPOLAR
      IF(STRING(1:4).EQ.'WMIN')READ(2,*)wmin
      IF(STRING(1:4).EQ.'WMAX')READ(2,*)wmax
      IF(STRING(1:5).EQ.'EXCNM')READ(2,*)EXCNM
      IF(STRING(1:5).EQ.'COMPL')READ(2,*)lcompl
      IF(STRING(1:5).EQ.'RAMAN')READ(2,*)RAMAN
      IF(STRING(1:5).EQ.'FILEQ')READ(2,*)FILEQ
      IF(STRING(1:3).EQ.'TTT')READ(2,*)
      IF(STRING(1:3).EQ.'PRO')READ(2,*)lpro
      IF(STRING(1:2).EQ.'N6')READ(2,*)N6
      IF(STRING(1:4).EQ.'GPRO')READ(2,*)lgpro
      IF(STRING(1:4).EQ.'APRO')READ(2,*)lapro
      IF(STRING(1:4).EQ.'UNIT')READ(2,*)units
      IF(STRING(1:4).EQ.'TEMP')READ(2,*)temp
      IF(STRING(1:4).EQ.'REFR')READ(2,*)ri
      IF(STRING(1:4).EQ.'INTA')READ(2,*)lintA
      IF(STRING(1:4).EQ.'INTG')READ(2,*)lintG
      IF(STRING(1:4).EQ.'USEA')READ(2,*)luseA
      IF(STRING(1:4).EQ.'USEG')READ(2,*)luseG
      IF(STRING(1:3).EQ.'XYZ')READ(2,*)XYZ
      IF(STRING(1:5).EQ.'ORIEN')READ(2,*)lor
      IF(STRING(1:3).EQ.'END')GOTO 2002
      GOTO 2001
2002  CLOSE(2)
     
99    write(6,6009)EXCNM
 6009 format(' Excitation ',f5.1,' nm')
      EXCA=EXCNM*10.0D0
      RETURN
      END

      SUBROUTINE DORAMAN(N,ALPHA,G,A,EXCA,wmin,wmax,E,units,temp,ri,
     1IKIND,luseA,luseG,lor)
c
c     supposedly G/W is supplied
c
      IMPLICIT none
      integer*4 IQ,I,J,K,L,N,IKIND,units
      real*8 ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),E(*),EXCA,wmin,wmax,temp,
     1ri,CM,ECM,SAL0,SAL1,SAG0,SAG1,SA1,SPAL,SDCP180,SDCPII180,
     1S180,D180,D,S90X,D90X,DX,S90Z,D90Z,DZ,YDY,YDX,P1,down,doc,
     1ro,AMU,BOHR,SpA3,beta2,gpisvejc,alphag,betaa,
     1betag,alpha2,roa3,roa90x,ram90x,roa90z,ram90z,roa2,ram2,CID2,roa1,
     1ram1,pisvejc2,EPS,CID1,PISVEJC,IRR180,IRL180
      logical luseA,luseG,lor
      character*8 OF(3)
      data OF/'ROAX.TAB','ROAY.TAB','ROAZ.TAB'/
c
      CM=219470.0d0
      if(.not.luseA)then
       write(6,*)'A-tensor deleted'
       do 2 IQ=1,N
       do 2 I=1,3
       do 2 J=1,3
       do 2 K=1,3
2      A(IQ,I,J,K)=0.0d0
      endif

      if(.not.luseG)then
       write(6,*)'G-tensor deleted'
       do 4 IQ=1,N
       do 4 I=1,3
       do 4 J=1,3
4      G(IQ,I,J)=0.0d0
      endif
c
      if(lor)then
       do 5 i=1,3
       OPEN(9+i,FILE=OF(i))
       WRITE(9+i,3304)
5      write(9+i,3002)
      endif
      OPEN(4,FILE='ROA.TAB')
      OPEN(44,FILE='INV.TXT')
      WRITE(44,4304)
4304  FORMAT(4x,'IQ',9x,'W',5x,'Alpha2',6x,'Beta2',5x,'AlphaG',
     1'     Gamma2     Delta2')
      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',6x,'Ramtot',11x,'CID90',
     27x,'CID-X',6x,'CID180',3x,'DEP',5x,'ROA180',5x,'DOC',
     15x,'RamDCPI',5x,'RoaDCPI',3x,'RamnDCPII',4x,'IRR(180)',
     14x,'IRL(180)',/,10x,'cm-1')
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(149(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
       SPAL=0.0d0
       DO 3 I=1,3
       SPAL=SPAL+ALPHA(IQ,I,I)
       DO 3 J=1,3
       SAL0=SAL0+ALPHA(IQ,I,I)*ALPHA(IQ,J,J)
       SAL1=SAL1+ALPHA(IQ,I,J)*ALPHA(IQ,I,J)
       SAG0=SAG0+ALPHA(IQ,I,I)*    G(IQ,J,J)
       SAG1=SAG1+ALPHA(IQ,I,J)*    G(IQ,I,J)
       DO 3 K=1,3
       DO 3 L=1,3
3      SA1=SA1+EPS(I,K,L)*ALPHA(IQ,I,J)*A(IQ,K,L,J)/3.0d0
c
c      IL+IR, backscattering, DCP_I,DCP_II
       SDCP180  =6.0d0*(6.0d0*SAL1-2.0d0*SAL0)
       SDCPII180=6.0d0*(      SAL1+3.0d0*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=        7.0d0*SAL1+ 1.0d0*SAL0
       doc=5.0d0*(  1.0d0*SAL1- 1.0d0*SAL0)
       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
        Alpha2=SAL0/9.0d0*(AMU*BOHR**4)
        SpA3=SPAL/3.0d0*dsqrt((AMU*BOHR**4))
        beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4)
        SDCP180  =SDCP180  *(AMU*BOHR**4)
        SDCPII180=SDCPII180*(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
        betag=(3.0d0*SAG1-SAG0)/2.0d0    *gpisvejc
        WRITE(44,4305)IQ,ECM,Alpha2,beta2,alphag,betag,betaa
4305    FORMAT(i6,f10.2,5f11.3)
c
c       DCPI 180
c       ram3=24.0d0*beta2 = 12 (3al_ab al_ab -al_aa al_bb)
        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*beta2)
c       depolarized ROA,90deg:
        roa90z=6.0d0*(betag-betaa/3.0d0)
        ram90z=0.5d0*6.0d0*beta2
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/SCP(180):
c       8/c*(12*beta(G')^2 + 4*beta(A)^2):
        roa1=96.0d0*(betag+betaa/3.0d0)
c       ram1 = 6(7 al_ab al_ab + al_aa al_bb)
        ram1=4.0d0*(45.0d0*SpA3**2+7.0d0*beta2)
        CID1=0.0d0
        if(ram1.gt.1.0d-9)CID1=roa1/ram1

        IRR180 = 24.0d0*beta2
        IRL180 =180.0d0*SpA3**2+4.0d0*beta2

c       NON SCP modes:
        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

c       Josef's units:
        if(units.eq.1)then
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,
     1  roa1,doc,SDCP180,roa3,SDCPII180,IRR180,IRL180
3000    FORMAT(I5,f9.2,6e12.4,f6.3,e12.4,f7.3,7e12.4)
        if(lor)then
         do 6 i=1,3
         call setjk(i,j,k)
c        S0(x)=(1/2)(azz azz + ayy ayy + azy azy + azy azy)
c        S0(y)=(1/2)(axx axx + azz azz + axz axz + axz axz)
c        S0(z)=(1/2)(axx axx + ayy ayy + axy axy + axy axy)
c                     jj  jj    kk  kk    jk  jk    jk  jk
         ram1=180.0d0*0.5d0*AMU*BOHR**4*(
     1   ALPHA(IQ,j,j)*ALPHA(IQ,j,j)+ALPHA(IQ,k,k)*ALPHA(IQ,k,k)+
     1   ALPHA(IQ,j,k)*ALPHA(IQ,j,k)+ALPHA(IQ,j,k)*ALPHA(IQ,j,k))
c       -2ayx(Gxy+Gyx)+(ayy-axx)(Gxx-Gyy)
c          kj  jk  kj    kk  jj   jj  kk
c        +(2w/3)(ayyAx zy-axxAy zx+ayxAx zx-axyAy zy)
c                 kk j ik  jj k ij  kj j ij  jk k ik
         roa1=-180.0d0*gpisvejc*(
     1   -2.0d0*ALPHA(IQ,j,k)*(G(IQ,j,k)+G(IQ,k,j))
     1   +(ALPHA(IQ,k,k)-ALPHA(IQ,j,j))*(G(IQ,j,j)-G(IQ,k,k))
     1   +2.0d0/3.0d0*
     1   (ALPHA(IQ,k,k)*A(IQ,j,i,k)-ALPHA(IQ,j,j)*A(IQ,k,i,j)
     1   +ALPHA(IQ,k,j)*A(IQ,j,i,j)-ALPHA(IQ,j,k)*A(IQ,k,i,k)))
6        WRITE(9+i,3000)IQ,ECM,YDX,YDY,ram1,CID2,DX,CID1,P1,roa1
        endif
       else
        PISVEJC=1500.0d0/7.0d0
        S180=PISVEJC*S180
        SDCP180=PISVEJC*SDCP180
        SDCPII180=PISVEJC*SDCPII180
        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
        SDCPII180=SDCPII180*1000.0d0
        doc=(3*P1-1)/(P1+1)
        WRITE(4,3000)IQ,ECM,YDX,YDY,S180,DZ,DX,D,P1,D180,doc,SDCP180,
     1  SDCPII180
       endif
c
      endif
      endif
      ENDIF
1     CONTINUE
      WRITE(4,3002)
      CLOSE(4)
      CLOSE(44)
      if(lor)then
       do 7 i=1,3
       WRITE(9+i,3002)
7      CLOSE(9+i)
      endif
      RETURN
      END

      SUBROUTINE DORAMANi(N,ALr,ALi,Gr,Gi,Ar,Ai,EXCA,wmin,wmax,
     1E,units,temp,ri,IKIND,luseA,luseG,lor)
c
c     supposedly G/W is supplied
c
      IMPLICIT none
      integer*4 IQ,I,J,K,L,N,IKIND,units
      real*8 ALr(N,3,3),Gr(N,3,3),Ar(N,3,3,3),E(*),EXCA,wmin,wmax,
     1ALi(N,3,3),Gi(N,3,3),Ai(N,3,3,3),temp,
     1ri,CM,ECM,SAL0,SAL1,SAG0,SAG1,SA1,SDCP180,SDCPII180,
     1S180,D180,D,S90X,D90X,DX,S90Z,D90Z,DZ,YDY,YDX,P1,down,doc,
     1ro,AMU,BOHR,beta2,gpisvejc,alphag,betaa,
     1betag,alpha2,roa3,roa90x,ram90x,roa90z,ram90z,roa2,ram2,CID2,roa1,
     1ram1,pisvejc2,EPS,CID1,PISVEJC
      logical luseA,luseG,lor
      character*8 OF(3)
      data OF/'ROAX.TAB','ROAY.TAB','ROAZ.TAB'/
c
      CM=219470.0d0
      if(.not.luseA)then
       write(6,*)'A-tensor deleted'
       Ar=0.0d0
       Ai=0.0d0
      endif

      if(.not.luseG)then
       write(6,*)'G-tensor deleted'
       Gr=0.0d0
       Gi=0.0d0
      endif
c
      if(lor)then
       do 5 i=1,3
       OPEN(9+i,FILE=OF(i))
       WRITE(9+i,3304)
5      write(9+i,3002)
      endif
      OPEN(44,FILE='INV.TXT')
      WRITE(44,4304)
4304  FORMAT(4x,'IQ',9x,'W',5x,'Alpha2',6x,'Beta2',5x,'AlphaG',
     1'     Gamma2     Delta2')
      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',6x,'Ramtot',11x,'CID90',
     27x,'CID-X',6x,'CID180',3x,'DEP',5x,'ROA180',5x,'DOC',
     15x,'RamDCPI',5x,'RoaDCPI   RamnDCPII',/,
     3'          cm-1')
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(149(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
       SAL0=SAL0+ALr(IQ,I,I)*ALr(IQ,J,J)+ALi(IQ,I,I)*ALi(IQ,J,J)
       SAL1=SAL1+ALr(IQ,I,J)*ALr(IQ,I,J)+ALi(IQ,I,J)*ALi(IQ,I,J)

c      - Im(alpha G*)
       SAG0=SAG0+ALr(IQ,I,I)* Gi(IQ,J,J)-ALi(IQ,I,I)* Gr(IQ,J,J)
       SAG1=SAG1+ALr(IQ,I,J)* Gi(IQ,I,J)-ALi(IQ,I,J)* Gr(IQ,I,J)

       DO 3 K=1,3
       DO 3 L=1,3
3      SA1=SA1+EPS(I,K,L)*( ALr(IQ,I,J)*Ar(IQ,K,L,J)
     1                     +ALi(IQ,I,J)*Ai(IQ,K,L,J))/3.0d0
c
c      IL+IR, backscattering, DCP_I,DCP_II
       SDCP180  =6.0d0*(6.0d0*SAL1-2.0d0*SAL0)
       SDCPII180=6.0d0*(      SAL1+3.0d0*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=        7.0d0*SAL1+ 1.0d0*SAL0
       doc=5.0d0*(  1.0d0*SAL1- 1.0d0*SAL0)
       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
        alpha2=SAL0/9.0d0*(AMU*BOHR**4)
        beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4)
        SDCP180  =SDCP180  *(AMU*BOHR**4)
        SDCPII180=SDCPII180*(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
        betag=(3.0d0*SAG1-SAG0)/2.0d0    *gpisvejc
        betaa=SA1*3.0d0/2.0d0            *gpisvejc

        WRITE(44,4305)IQ,ECM,alpha2,beta2,alphag,betag,betaa
4305    FORMAT(i6,f10.2,5f11.3)
c
c       DCPI 180
c       ram3=24.0d0*beta2 = 12 (3al_ab al_ab -al_aa al_bb)
        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*beta2)
c       depolarized ROA,90deg:
        roa90z=6.0d0*(betag-betaa/3.0d0)
        ram90z=0.5d0*6.0d0*beta2
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/SCP(180):
c       8/c*(12*beta(G')^2 + 4*beta(A)^2):
        roa1=96.0d0*(betag+betaa/3.0d0)
c       ram1 = 6(7 al_ab al_ab + al_aa al_bb)
        ram1=4.0d0*(45.0d0*alpha2+7.0d0*beta2)
        CID1=0.0d0
        if(ram1.gt.1.0d-9)CID1=roa1/ram1

c       NON SCP modes:
        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

c       Josef's units:
        if(units.eq.1)then
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,
     1  roa1,doc,SDCP180,roa3,SDCPII180
3000    FORMAT(I5,f9.2,6e12.4,f6.3,e12.4,f7.3,5e12.4)
        if(lor)then
         do 6 i=1,3
         call setjk(i,j,k)
c        S0(x)=(1/2)(azz azz + ayy ayy + azy azy + azy azy)
c        S0(y)=(1/2)(axx axx + azz azz + axz axz + axz axz)
c        S0(z)=(1/2)(axx axx + ayy ayy + axy axy + axy axy)
c                     jj  jj    kk  kk    jk  jk    jk  jk
         ram1=180.0d0*0.5d0*AMU*BOHR**4*(
     1   ALr(IQ,j,j)*ALr(IQ,j,j)+ALr(IQ,k,k)*ALr(IQ,k,k)+
     1   ALr(IQ,j,k)*ALr(IQ,j,k)+ALr(IQ,j,k)*ALr(IQ,j,k)+
     1   ALi(IQ,j,j)*ALi(IQ,j,j)+ALi(IQ,k,k)*ALi(IQ,k,k)+
     1   ALi(IQ,j,k)*ALi(IQ,j,k)+ALi(IQ,j,k)*ALi(IQ,j,k))
c       -2ayx(Gxy+Gyx)+(ayy-axx)(Gxx-Gyy)
c          kj  jk  kj    kk  jj   jj  kk
c        +(2w/3)(ayyAx zy-axxAy zx+ayxAx zx-axyAy zy)
c                 kk j ik  jj k ij  kj j ij  jk k ik
         roa1=-180.0d0*gpisvejc*(
     1   -2.0d0*ALr(IQ,j,k)*(Gr(IQ,j,k)+Gr(IQ,k,j))
     1   +(ALr(IQ,k,k)-ALr(IQ,j,j))*(Gr(IQ,j,j)-Gr(IQ,k,k))
     1   -2.0d0*ALi(IQ,j,k)*(Gi(IQ,j,k)+Gi(IQ,k,j))
     1   +(ALi(IQ,k,k)-ALi(IQ,j,j))*(Gi(IQ,j,j)-Gi(IQ,k,k))
     1   +2.0d0/3.0d0*
     1   (ALr(IQ,k,k)*Ar(IQ,j,i,k)-ALr(IQ,j,j)*Ar(IQ,k,i,j)
     1   +ALr(IQ,k,j)*Ar(IQ,j,i,j)-ALr(IQ,j,k)*Ar(IQ,k,i,k))
     1   +2.0d0/3.0d0*
     1   (ALi(IQ,k,k)*Ai(IQ,j,i,k)-ALi(IQ,j,j)*Ai(IQ,k,i,j)
     1   +ALi(IQ,k,j)*Ai(IQ,j,i,j)-ALi(IQ,j,k)*Ai(IQ,k,i,k)))
6        WRITE(9+i,3000)IQ,ECM,YDX,YDY,ram1,CID2,DX,CID1,P1,roa1
        endif
       else
        PISVEJC=1500.0d0/7.0d0
        S180=PISVEJC*S180
        SDCP180=PISVEJC*SDCP180
        SDCPII180=PISVEJC*SDCPII180
        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
        SDCPII180=SDCPII180*1000.0d0
        doc=(3*P1-1)/(P1+1)
        WRITE(4,3000)IQ,ECM,YDX,YDY,S180,DZ,DX,D,P1,D180,doc,SDCP180,
     1  SDCPII180
       endif
c
      endif
      endif
      ENDIF
1     CONTINUE
      WRITE(4,3002)
      CLOSE(4)
      CLOSE(44)
      if(lor)then
       do 7 i=1,3
       WRITE(9+i,3002)
7      CLOSE(9+i)
      endif
      RETURN
      END

      SUBROUTINE TRTEN(S,N,ALPHA,A,G)
      IMPLICIT REAL*8 (A-H,O-Z)
      implicit integer*4 (i-n)
      DIMENSION S(N,N),A(N,3,3,3),ALPHA(N,3,3),G(N,3,3)
      real*8,allocatable::V(:)
      allocate(V(N))
      DO 1 I=1,3
      DO 1 J=1,3
      DO 3 JQ=1,N
3     V(JQ)=ALPHA(JQ,I,J)
      DO 1 IQ=1,N
      ALPHAS=0.0d0
      DO 2 IX=1,N
2     ALPHAS=ALPHAS+S(IX,IQ)*V(IX)
1     ALPHA(IQ,I,J)=ALPHAS
      DO 11 I=1,3
      DO 11 J=1,3
      DO 31 JQ=1,N
31    V(JQ)=G(JQ,I,J)
      DO 11 IQ=1,N
      GS=0.0d0
      DO 21 IX=1,N
21    GS=GS+S(IX,IQ)*V(IX)
11    G(IQ,I,J)=GS
      DO 12 I=1,3
      DO 12 J=1,3
      DO 12 K=1,3
      DO 32 JQ=1,N
32    V(JQ)=A(JQ,I,J,K)
      DO 12 IQ=1,N
      AS=0.0D0
      DO 22 IX=1,N
22    AS=AS+S(IX,IQ)*V(IX)
12    A(IQ,I,J,K)=AS
      WRITE(3,3001)
3001  FORMAT(1X,'TENSORS HAVE BEEN TRANSFORMED TO Q COORDINATES',/,
     11X,'ALPHA DERIVATIVES LISTED:')
      DO 4 IQ=1,N
      SP=ALPHA(IQ,1,1)+ALPHA(IQ,2,2)+ALPHA(IQ,3,3)
4     WRITE(3,3000)IQ,SP,((ALPHA(IQ,I,J),I=1,3),J=1,3)
3000  FORMAT(1X,I4,' MODE; Trace: ',f15.8,/,3(3F15.8,/),/)
      WRITE(3,3002)
3002  FORMAT(1X, 11X,'G (magn index -->)')
      DO 5 IQ=1,N
      SP=G(IQ,1,1)+G(IQ,2,2)+G(IQ,3,3)
5     WRITE(3,3000)IQ,SP,((G(IQ,J,I),I=1,3),J=1,3)
      WRITE(3,3003)
3003  FORMAT(1X, 11X,'A DERIVATIVES LISTED:')
      DO 6 IQ=1,N
6     WRITE(3,3004)IQ,(((A(IQ,I,J,K),I=1,3),J=1,3),K=1,3)
3004  FORMAT(1X,I4,' MODE',/,3(9F8.5,/),/)
      RETURN
      END
c     ============================================================
      SUBROUTINE readpol2(a,fn)
      character*(*)fn
      real*8 a(3,3)
      integer*4 i,j
      logical lex
      inquire(file=fn,exist=lex)
      if(lex)then
       open(90,file=fn)
       read(90,*)
       read(90,*)
       read(90,*)((a(i,j),i=1,j),j=1,3)
       close(90)
       do 1 j=1,3
       do 1 i=1,j
1      a(j,i)=a(i,j)
       close(90)
      else
       a=0.0d0
       write(6,*)fn//' not found'
      endif
      return
      end
c     ============================================================
      SUBROUTINE readpolg(a,fn)
      character*(*)fn
      real*8 a(3,3)
      integer*4 i,j
      logical lex
      inquire(file=fn,exist=lex)
      if(lex)then
       open(90,file=fn)
       read(90,*)
       read(90,*)
       read(90,*)((a(j,i),i=1,3),j=1,3)
       close(90)
      else
       a=0.0d0
       write(6,*)fn//' not found'
      endif
      return
      end
c     ============================================================
      SUBROUTINE readpola(a,fn)
      character*(*)fn
      real*8 a(3,3,3)
      integer*4 i,j,k
      logical lex
      inquire(file=fn,exist=lex)
      if(lex)then
       open(90,file=fn)
       read(90,*)
       read(90,*)
       do 1 k=1,3
       read(90,*)((a(i,j,k),i=1,j),j=1,3)
       do 1 j=1,3
       do 1 i=1,j
1      a(j,i,k)=a(i,j,k)
       close(90)
      else
       a=0.0d0
       write(6,*)fn//' not found'
      endif
      return
      end
c     ============================================================
      SUBROUTINE DORAYLEIGH
c
c     supposedly G/W is supplied
c     see DORAMAN for the meaning of the symbols
c
      IMPLICIT none
      integer*4 N,npol,i,ix,iy,iz,IQ,i1,J,K,L
      real*8,allocatable::ALPHA(:,:,:),G(:,:,:),A(:,:,:,:),E(:)
      real*8 al(3,3),gp(3,3),ap(3,3,3),CM,ECM,ENM,SAL0,SAL1,SAG0,SAG1,
     1SA1,SPAL,S180,D180,D,YDY,YDX,P1,down,doc,ro,AMU,BOHR,
     1SpA3,beta2,gpisvejc,betaa,betag,roa1,ram1,cid1,EPS,EXCA
      character*10 s10
      character*80 f1,f2,f3
      logical lex
c
      N=npol('POL.TTT')
      if(N.eq.0)goto 999
      allocate(ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),E(N))

      inquire(file='EP.LST',exist=lex)
      if(.not.lex)goto 999
      open(7,file='EP.LST')
      read(7,*,err=999,end=999)
      read(7,701)(E(i),i=1,N)
701   format(5f12.6)
      close(7)

      do 2 i=1,N
      write(s10,'(i10)')i
      if(i.eq.1)then
       f1='POL.TTT'
       f2='GP.TTT'
       f3='A.TTT'
      else
       i1=10
       if(i.gt.9)i1=9
       if(i.gt.99)i1=8
       if(i.gt.999)i1=7
       f1='POL.TTT.f'//s10(i1:10)
       f2='GP.TTT.f'//s10(i1:10)
       f3='A.TTT.f'//s10(i1:10)
      endif       
      call readpol2(al,f1)
      call readpolg(gp,f2)
      call readpola(ap,f3)
      do 2 ix=1,3
      do 2 iy=1,3
      ALPHA(i,ix,iy)=al(ix,iy)
      G(i,ix,iy)=gp(ix,iy)
      do 2 iz=1,3
2     A(i,ix,iy,iz)=ap(ix,iy,iz)

      CM=219470.0d0
c
      OPEN(4,FILE='RAYLEIGH.TAB')
      WRITE(4,3304)
3304  FORMAT(
     1'number   FREQ    Rx Ry R (180)  ... dR (180)     CID',/,
     3'         nm      A^4/AMU    x10^4')
      write(4,3002)
3002  FORMAT(80(1H-))
c
      DO 1 IQ=1,N
      ECM=E(IQ)*CM
      ENM=1.0d7/ECM
      EXCA=ENM*10.0d0
      SAL0=0.0d0
      SAL1=0.0d0
      SAG0=0.0d0
      SAG1=0.0d0
      SA1=0.0d0
      SPAL=0.0d0
      DO 3 I=1,3
      SPAL=SPAL+ALPHA(IQ,I,I)
      DO 3 J=1,3
      SAL0=SAL0+ALPHA(IQ,I,I)*ALPHA(IQ,J,J)
      SAL1=SAL1+ALPHA(IQ,I,J)*ALPHA(IQ,I,J)
      SAG0=SAG0+ALPHA(IQ,I,I)*    G(IQ,J,J)
      SAG1=SAG1+ALPHA(IQ,I,J)*    G(IQ,I,J)
      DO 3 K=1,3
      DO 3 L=1,3
3     SA1=SA1+EPS(I,K,L)*ALPHA(IQ,I,J)*A(IQ,K,L,J)/3.0d0
      S180=7.0d0*SAL1+SAL0
      D180=8.0d0*(3.0d0*SAG1-SAG0+SA1)
      D=D180
      IF(S180.gt.1.0d-9)D=D/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
      down=        1.0d0*SAL0+ 7.0d0*SAL1
      doc=5.0d0*(  1.0d0*SAL0- 1.0d0*SAL1)
      if(down.ne.0.0d0)doc=doc/down
      down= 6.0d0*SAL0-12.0d0*SAL1
      ro  =-3.0d0*SAL0+ 9.0d0*SAL1
      if(down.ne.0.0d0)ro=ro/down
      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)/EXCA
      betaa=SA1*3.0d0/2.0d0*gpisvejc
      betag=(3.0d0*SAG1-SAG0)/2.0d0*gpisvejc
      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)
1     WRITE(4,3000)IQ,ENM,YDX,YDY,ram1,ram1,ram1,ram1,P1,roa1,doc,CID1
3000  FORMAT(I5,f9.2,6g12.3,f6.3,4g12.4)
      WRITE(4,3002)
      CLOSE(4)
999   RETURN
      END


      SUBROUTINE READALV(ALPHA)
      IMPLICIT REAL*8 (A-H,O-Z)
      implicit integer*4 (i-n)
      DIMENSION ALPHA(3,3)
      CHARACTER*50 ST
      IC=0
201   READ(2,1111)ST
      IC=IC+1
1111  FORMAT(A50)
C                    12345678901234567890
      IF(ST.EQ.' POLARISABILITY IN LENGTH-VELOCITY REPRESENTATION ')THEN
       WRITE(3,*)ST
       READ(2,*)
       READ(2,*)
       DO 203 I=1,3
       READ(2,*)
203    READ(2,2000)(ALPHA(I,J),J=1,3)
2000  FORMAT(12X,3F20.13)
       RETURN
      ENDIF
      IF(IC.GT.20000)THEN
       WRITE(2,*)' AlphaV not found, program stopped ... '
       STOP
      ENDIF
      GOTO  201
      END
c     ====================
      subroutine readpol(a)
      real*8 a(3,3)
      open(90,file='POL.TTT')
      read(90,*)
      read(90,*)
      read(90,*)a(1,1),a(1,2),a(2,2),a(1,3),a(2,3),a(3,3)
      close(90)
      a(2,1)=a(1,2)
      a(3,1)=a(1,3)
      a(3,2)=a(2,3)
      return
      end
c     ====================
      function npol(f)
      implicit none
      integer*4 npol,i
      character*(*) f
      character*80 s80
      character*10 s10
      logical lex
      i=1
1     write(s10,'(i10)')i
      if(i.eq.1)then
       s80=f
      else
       s80=f//'.f'//s10(10:10)
       if(i.gt.9)s80=f//'.f'//s10(9:10)
       if(i.gt.99)s80=f//'.f'//s10(8:10)
       if(i.gt.999)s80=f//'.f'//s10(7:10)
      endif
      inquire(file=s80,exist=lex)
      if(lex)then
       i=i+1
       goto 1
      endif
      i=i-1
      write(6,*)i,' frequency-dependent polarizabilities'
      npol=i
      return
      end
c     ====================
      SUBROUTINE rules(ALPHA,A,G,NAT)
C     project translations/rotations
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(3*NAT,3,3),G(3*NAT,3,3),A(3*NAT,3,3,3),
     1ap(3,3),aap(3,3,3),APOL(3,3)
      character*2 xyz(3)
      logical lpol
      data xyz/' X',' Y',' Z'/
c
      write(6,6077)' Sum rules-alpha:'
6077  format(a17)
      do 11 ii=1,3
      do 21 ix=1,3
      do 21 jx=1,3
      ap(ix,jx)=0.0d0
      do 21 ia=2,NAT
21    ap(ix,jx)=ap(ix,jx)+ALPHA((ia-1)*3+ii,ix,jx)
      write(6,6017)xyz(ii)
6017  format(a2)
      write(6,600)((  -ap(ix,jx),jx=1,3),ix=1,3)
600   format(9F8.2)
11    write(6,600)((ALPHA(ii,ix,jx),jx=1,3),ix=1,3)

      inquire(file='POL.TTT',exist=lpol)
      if(lpol)then
       call readpol(APOL)
       write(6,*)'Static polarizability read in'
      else
       do 61 ix=1,3
       do 61 jx=1,3
61     APOL(ix,jx)=0.0d0
      endif

      write(6,6077)' Sum rules-G    :'
      do 1 ii=1,3
      do 2 ix=1,3
      do 2 jx=1,3
      ap(ix,jx)=0.0d0
      do 2 ia=2,NAT
2     ap(ix,jx)=ap(ix,jx)+    G((ia-1)*3+ii,ix,jx)
      write(6,6017)xyz(ii)
      write(6,600)((-ap(ix,jx),jx=1,3),ix=1,3)
      do 3 ix=1,3
      do 3 jx=1,3
      ap(ix,jx)=0.0d0
      do 3 kx=1,3
3     ap(ix,jx)=ap(ix,jx)+0.5d0*EPS(jx,ii,kx)*APOL(ix,kx)
1     write(6,600)((G(ii,ix,jx)+ap(ix,jx),jx=1,3),ix=1,3)
c
      write(6,6077)' Sum rules-A    :'
      do 14 ii=1,3
      do 12 ix=1,3
      do 12 jx=1,3
      do 12 kx=1,3
      aap(ix,jx,kx)=0.0d0
      do 12 ia=2,NAT
12    aap(ix,jx,kx)=aap(ix,jx,kx)+A(3*(ia-1)+ii,ix,jx,kx)
      do 14 kx=1,3
      write(6,6018)xyz(ii),xyz(kx)
6018  format(2a2)
      write(6,600)((-aap(ix,jx,kx),jx=1,3),ix=1,3)
      do 41 ix=1,3
      do 41 jx=1,3
      ap(ix,jx)=0.0d0
      if(jx.eq.kx)ap(ix,jx)=ap(ix,jx)      +APOL(ix,ii)
      if(ii.eq.jx)ap(ix,jx)=ap(ix,jx)-1.5d0*APOL(ix,kx)
41    if(ii.eq.kx)ap(ix,jx)=ap(ix,jx)-1.5d0*APOL(ix,jx)
14    write(6,600)((A(ii,ix,jx,kx)+ap(ix,jx),jx=1,3),ix=1,3)

      return
      end
c     ====================
      SUBROUTINE PROJECT(ALPHA,A,G,NAT,R,lpro,lgpro,lapro,N6)
C     project translations/rotations
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(3*NAT,3,3),G(3*NAT,3,3),A(3*NAT,3,3,3),R(3,NAT)
      real*8,allocatable::aproj(:,:),t(:)
      logical lpro,lapro,lgpro
      allocate(aproj(3*NAT,3*NAT),t(3*NAT))

      call DOMA(aproj,R,NAT,N6)

      if(lpro)then
       do 1 ix=1,3
       do 1 jx=1,3
       do 2 i=1,3*NAT
       t(i)=0.0d0
       do 2 j=1,3*NAT
2      t(i)=t(i)+aproj(i,j)*ALPHA(j,ix,jx)
       do 1 i=1,3*NAT
1      ALPHA(i,ix,jx)=t(i)
       write(6,*)' Alpha derivatives projected'
      else
       write(6,*)' Alpha projection skipped'
      endif

      if(lgpro)then
       do 3 ix=1,3
       do 3 jx=1,3
       do 4 i=1,3*NAT
       t(i)=0.0d0
       do 4 j=1,3*NAT
4      t(i)=t(i)+aproj(i,j)*G(j,ix,jx)
       do 3 i=1,3*nat
3      G(i,ix,jx)=t(i)
       write(6,*)' G derivatives projected'
      else
       write(6,*)' G projection skipped'
      endif

      if(lapro)then
       do 5 ix=1,3
       do 5 jx=1,3
       do 5 kx=1,3
       do 6 i=1,3*NAT
       t(i)=0.0d0
       do 6 j=1,3*NAT
6      t(i)=t(i)+aproj(i,j)*A(j,ix,jx,kx)
       do 5 i=1,3*nat
5      A(i,ix,jx,kx)=t(i)
       write(6,*)' A derivatives projected'
      else
       write(6,*)' A projection skipped'
      endif

      return

      end
c     ====================
      SUBROUTINE TTTT(N,ALPHA,A,G,NAT,ICO,R,X)
c
C     ICO = 1: Transforms A and G to local origins
C     ICO = 2: Transforms A and G to common origin
C     ICO = 3: Sets A and G to zero 
c     ICO = 4: Delete G
c     ICO = 5: Delete A
C     R supposed to be passed in AU
C
      IMPLICIT none
      INTEGER*4 N,NAT,ICO,I,J,L,IE,II,K,M,IA,IB,IC,ID
      real*8 ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),R(3,N/2),X(3,N/3),
     1S,EPS,SUM
C
      DO 6 I=1,3
      DO 6 J=1,NAT
6     X(I,J)=R(I,J)
C
      IF(ICO.EQ.1)then
       S= 1.0d0
       WRITE(6,*)' A,G transformed into atomic origins'
      elseif(ICO.EQ.2)then
       S=-1.0d0
       WRITE(6,*)' A,G transformed into common origin'
      elseif(ICO.ge.3.and.ICO.LE.5)then
       IF(ICO.eq.3.or.ICO.EQ.4)then
        do 5 II=1,N
        do 5 I=1,3
        do 5 J=1,3
5       G(II,I,J)=0.0d0
        WRITE(6,*)' G set to zero'
       endif
       IF(ICO.eq.3.or.ICO.EQ.5)then
        do 7 II=1,N
        do 7 I=1,3
        do 7 J=1,3
        do 7 K=1,3
7       A(II,I,J,K)=0.0d0
        WRITE(6,*)' A set to zero'
       endif
       return
      else
       write(6,*)ICO
       call report('TTTT unknown instruction')
      endif
C
      DO 1 I=1,3
      DO 1 J=1,3
      DO 1 L=1,NAT
      DO 1 IE=1,3
      II=3*(L-1)+IE
      DO 1 K=1,3
      DO 1 M=1,3
c     G-second index-magnetic
1     G(II,I,J)=G(II,I,J)+S*0.5d0*EPS(J,K,M)*X(K,L)*ALPHA(II,I,M)
C
      DO 3 IA=1,3
      DO 3 IB=1,3
      DO 3 IC=1,3
      DO 3 L=1,NAT
      DO 3 IE=1,3
      II=3*(L-1)+IE
      SUM=0.0d0
      IF(IB.EQ.IC)THEN
       DO 4 ID=1,3
4      SUM=SUM+X(ID,L)*ALPHA(II,IA,ID)
      ENDIF
      A(II,IA,IB,IC)=A(II,IA,IB,IC)+SUM*S
     1-1.5d0*S*(X(IB,L)*ALPHA(II,IA,IC)+X(IC,L)*ALPHA(II,IA,IB))
3     CONTINUE
c
      RETURN
      END
c     ====================================
      SUBROUTINE DOMA(A,C,NAT,N6)
c     N6 - number of coordinates to be projected
c     1..3 translations, 4..6 rotations
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION C(3,NAT),A(3*NAT,3*NAT)
      real*8,allocatable::TEM(:,:)
      allocate(TEM(3*NAT,3*NAT))
      N=NAT*3
      AMACH=1.0d-12
      DO 12 I=1,N
      DO 12 J=1,N
12    A(I,J)=0.0d0
      DO 1 I=1,NAT
      A(3*I-2,1)= 1.00d0
      A(3*I-1,2)= 1.00d0
      A(3*I  ,3)= 1.00d0
      A(3*I-2,4)=-C(2,I)
      A(3*I-1,4)= C(1,I)
      A(3*I-1,5)=-C(3,I)
      A(3*I  ,5)= C(2,I)
      A(3*I-2,6)= C(3,I)
1     A(3*I  ,6)=-C(1,I)
      CALL ORT(A,N,N6,NAT)
      ICOL=1
5     IF(ICOL.GT.N6)GOTO 8
      S=0.0d0
      DO 6 J=1,N
6     S=S+A(J,ICOL)**2
      IF(S.LT.1000.0d0*AMACH)THEN
       DO 7 J=1,N
7      A(J,ICOL)=A(J,N6)
       N6=N6-1
       GOTO 5
      ENDIF
      ICOL=ICOL+1
      GOTO 5
8     CONTINUE
      WRITE(*,*)N6,' coordinates projected'
      DO 2 I=1,N
      DO 2 J=1,N
      TEM(I,J)=0.0d0
      DO 2 II=1,N6
2     TEM(I,J)=TEM(I,J)+A(I,II)*A(J,II)
      DO 3 I=1,N
      DO 4 J=1,N
4     A(I,J)=-TEM(I,J)
3     A(I,I)=A(I,I)+1.0d0
      RETURN
      END
c     ===============================
      SUBROUTINE ORT(A,NDIM,N6,NAT)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(NDIM,NDIM)
      AMACH=1.0d-13
      NAT3=3*NAT
      C=0.0d0
      DO 1 I=1,NAT3
1     C=C+A(I,1)**2
      IF(C.GT.AMACH)C=1.0d0/SQRT(C)
      DO 2 I=1,NAT3
2     A(I,1)=A(I,1)*C
      DO 3 IC=2,N6
C     Column IC, orthonormalize to previous columns:
      DO 3 ICP=1,IC-1
      SP=0.0d0
      DO 5 J=1,NAT3
5     SP=SP+A(J,IC)*A(J,ICP)
C     Subtract projection to ICP:
      DO 6 J=1,NAT3
6     A(J,IC)=A(J,IC)-SP*A(J,ICP)
C     Normalize :
      SP=0.0d0
      DO 7 J=1,NAT3
7     SP=SP+A(J,IC)**2
      IF(SP.GT.AMACH)SP=1.0d0/SQRT(SP)
      DO 3 J=1,NAT3
3     A(J,IC)=A(J,IC)*SP
      RETURN
      END

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

      subroutine setjk(i,j,k)
      integer*4 i,j,k
      j=i+1
      if(j.gt.3)j=1
      k=j+1
      if(k.gt.3)k=1
      return
      end