PROGRAM NEW6C
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      PARAMETER (N30=4500)
      DIMENSION SI(N30,N30),SR(N30,N30),
     2ALPHA(N30,3,3),A(N30,3,3,3),G(N30,3,3),
     2ALPHAR(N30,3,3),AR(N30,3,3,3),GR(N30,3,3),
     2ALPHAI(N30,3,3),AI(N30,3,3,3),GI(N30,3,3),
     3R(3,N30/3),E(N30)
      integer*4 xstart,ystart,zstart
      CHARACTER*20 sx,sy,sz
      LOGICAL SCALED,RAMAN,TTT,XYZ,LPOLAR,LGAUSS,lpro,lgpro,lapro
C
      WRITE(6,6666)
6666  FORMAT(/,/,
     A'          * * * * * * * * * * * * * * * * * * * * * *',/,
     1'        *                                         *  ',/,
     2'      *       CRYSTAL RAMAN OPTICAL ACTIVITY    *    ',/,
     3'        *                                         *  ',/,
     4'          * * * * * * * * * * * * * * * * * * * * * *',/,
     5'                                                     ',/,
     6'                                                     ',/,
     9'                       PETR BOUR, Prague, 2012       ',/,/,
     A'INPUT FILES: ROA.OPT,Q.TXT,CELL.X,CELL.TTT,xxxF.INP',/,
     B'OUTPUT     : ROA.OUT, ROA.TAB ROA.TTT',/)
      N7=7
      N1=7
      N2=0
      CALL IOPAR(SCALED,TTT,RAMAN,XYZ,N1,N2,LPOLAR,LGAUSS,EXCA,
     1wmin,wmax,lpro,lgpro,lapro,N6)
      IF(TTT)CALL READTEN(N30,NAT,ALPHA,G,A)
      IF(XYZ)CALL READXYZ(R,NAT,N30)

      open(2,file='Q.TXT')
      read(2,*)qxmin,qxstep,nxstep
      read(2,*)qymin,qystep,nystep
      read(2,*)qzmin,qzstep,nzstep
      close(2)
      write(6,6005)qxmin,qxstep,nxstep,qymin,qystep,nystep,
     1qzmin,qzstep,nzstep,nxstep*nystep*nzstep
6005  format(/,' q-vectors:',/,
     1'        min/a     del q nsteps',/,
     1'X: ',f10.3,f10.3,i7,/,
     1'Y: ',f10.3,f10.3,i7,/,
     1'Z: ',f10.3,f10.3,i7,/,/,
     1i6,' in total',/)

      OPEN(44,FILE='ROAC.TAB')
      if(LGAUSS)then
       WRITE(44,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')
      else
       WRITE(44,3004)
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')
      endif
      write(44,3002)
3002  FORMAT(80(1H-))

      itr=0
      qx=qxmin-qxstep
      do 1001 iqx=1,nxstep
      qx=qx+qxstep
      qy=qymin-qystep
      do 1001 iqy=1,nystep
      qy=qy+qystep
      qz=qzmin-qzstep
      do 1001 iqz=1,nzstep
      qz=qz+qzstep
      write(6,6006)qx,qy,qz
6006  format('q:',3f10.3)
      call sset(sx,iqx,xstart)
      call sset(sy,iqy,ystart)
      call sset(sz,iqz,zstart)
c     read real and imaginary S-matrices:
      open(4,file=
     1sx(xstart:20)//'.'//sy(ystart:20)//'.'//sz(zstart:20)//'FR.INP')
      CALL READS(N30,SR,E,N7)
      close(4)
      open(4,file=
     1sx(xstart:20)//'.'//sy(ystart:20)//'.'//sz(zstart:20)//'FI.INP')
      CALL READS(N30,SI,E,N7)
      close(4)
      open(4,file=
     1sx(xstart:20)//'.'//sy(ystart:20)//'.'//sz(zstart:20)//'.TAB')
      WRITE(4,3304)
      
C     transform the tensors into internal coordinate derivatives
      CALL TRTEN(SR,SI,3*NAT,ALPHA,A,G,ALPHAR,AR,GR,ALPHAI,AI,GI)
c     calculate intensities:
      CALL DORAMAN(3*NAT,ALPHAR,GR,AR,LGAUSS,EXCA,wmin,wmax,E,itr)
      CALL DORAMAN(3*NAT,ALPHAI,GI,AI,LGAUSS,EXCA,wmin,wmax,E,itr)
      WRITE(4,3002)
1001  close(4)
      WRITE(44,3002)
      CLOSE(44)
      write(6,*)itr,' transitions'

      WRITE(6,*)' PROGRAM FINISHED OK'
      STOP
      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
      read(4,*)nint,nat3,nat
      N=NAT3
      N7=N-NINT+1
      do 1 i=1,NAT
1     read(4,*)
      read(4,*)
      DO 2 I=1,NAT
      DO 2 J=1,NINT
2     read(4,*)kdumm,kdumm,s(3*(i-1)+1,j+N7-1),
     1s(3*(i-1)+2,j+N7-1),s(3*(i-1)+3,j+N7-1)
      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
      NAT=N/3
      RETURN
      end
C     ==============================================================
      SUBROUTINE READTEN(N30,NAT,ALPHA,G,A)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION ALPHA(N30,3,3),G(N30,3,3),A(N30,3,3,3)
      OPEN(2,FILE='CELL.TTT',STATUS='OLD')
      READ(2,*)
      READ(2,*)NAT
      if(NAT.gt.N30/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
1     READ(2,*)LDUM,IXDUM,(ALPHA(IND,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
      IND=3*(L-1)+IX
c     last index magnetic, inner
2     READ(2,*)LDUM,IXDUM,(G(IND,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
      IND=3*(L-1)+IX
3     READ(2,*)LDUM,IXDUM,(A(IND,I,J,K),K=1,3)
      CLOSE(2)
      WRITE(6,*)NAT,' atoms, tensors read in from CELL.TTT'
      RETURN
      END
C     ==============================================================
C
      SUBROUTINE READXYZ(R,NAT,N30)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION R(3,N30/3)
      BOHR=0.529177d0
c
      OPEN(2,FILE='CELL.X',STATUS='OLD')
      READ(2,*)
      READ(2,*)nat
      do 1 i=1,nat
      READ(2,*)IDUMM,(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 CELL.X'
      RETURN
      END
C
      SUBROUTINE IOPAR(SCALED,TTT,RAMAN,XYZ,N1,N2,LPOLAR,LGAUSS,EXCA,
     1wmin,wmax,lpro,lgpro,lapro,N6)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      LOGICAL SCALED,TTT,RAMAN,XYZ,LPOLAR,LGAUSS,lpro,lex,lgpro,lapro
      CHARACTER*15 STRING
c
      wmin=50.0d0
      wmax=500000.0d0
      SCALED=.TRUE.
      lpro=.false.
      lapro=.false.
      lgpro=.false.
      RAMAN=.FALSE.
      N6=3
      LPOLAR=.FALSE.
      LGAUSS=.true.
      TTT=.TRUE.
      XYZ=.TRUE.
c     excitation light wavelength in nm:
      EXCNM=532.0d0

      INQUIRE(FILE='ROA.OPT',exist=lex)
      if(lex)then
       WRITE(6,*)' FILE ROA.OPT FOUND ...'
       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 ...'
        WRITE(6,*)' FILE NEW6.OPT FOUND ...'
       else
        goto 99
       endif
      endif
   
2001  READ(2,1111)STRING
1111  FORMAT(A15)
      WRITE(6,1111)STRING
      IF(STRING(1:4).EQ.'N1N2')THEN
       WRITE(6,*) ' Vibrations selected ...'
       READ(2,*)N1,N2
      ENDIF
      IF(STRING(1:6).EQ.'SCALED')READ(2,*)SCALED
      IF(STRING(1:5).EQ.'POLAR')READ(2,*)LPOLAR
      IF(STRING(1:5).EQ.'GAUSS')READ(2,*)LGAUSS
      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.'RAMAN')READ(2,*)RAMAN
      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:3).EQ.'XYZ')READ(2,*)XYZ
      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
c     ============================================================
      SUBROUTINE DORAMAN(N,ALPHA,G,A,LGAUSS,EXCA,wmin,wmax,E,itr)
c
c     supposedly G/W is supplied
c
      IMPLICIT REAL*8 (A-H,O-Z)
      implicit integer*4 (i-n)
      PARAMETER (N30=4500)
      DIMENSION ALPHA(N30,3,3),G(N30,3,3),A(N30,3,3,3),E(*)
      logical LGAUSS
c
      CM=219470.0d0
c
c
      DO 1 IQ=1,N
      ECM=E(IQ)*CM
      if(ECM.gt.wmin.and.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
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)/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.ne.0.0d0)then
        P1=YDY/YDX
       else
        P1=0.0d0
       endif
      

c
c      Try to get Gaussian units:
       if(LGAUSS)then
        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)
c       gpisvejc=980.0
        gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/EXCA
c
c       alphag=SAG0/9.0d0*gpisvejc
c
c       beta(A)^2=Delta2 in Gaussian:
        betaa=SA1*3.0d0/2.0d0*gpisvejc
        betag=(3.0d0*SAG1-SAG0)/2.0d0*gpisvejc
c
c       DCP 180
c       ram3=24.0d0*beta2
        roa3=8.0d0*(12.0d0*betag+4.0d0*betaa)
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)then
         CID1=roa1/ram1
         WRITE(44,3000)IQ,ECM,YDX,YDY,ram1,CID2,DX,CID1,P1,roa1,roa3
         WRITE(4,3000)IQ,ECM,YDX,YDY,ram1,CID2,DX,CID1,P1,roa1,roa3
3000     FORMAT(I5,f9.2,6g12.3,f6.3,2g12.4)
        endif
       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
        WRITE(44,3000)IQ,ECM,YDX,YDY,S180,DZ,DX,D,P1,D180,SDCP180
        WRITE(4,3000)IQ,ECM,YDX,YDY,S180,DZ,DX,D,P1,D180,SDCP180
       endif
       itr=itr+1
c
      endif
1     CONTINUE
      RETURN
      END
c     ============================================================
      SUBROUTINE TRTEN(SR,SI,N,ALPHA,A,G,ALPHAR,AR,GR,ALPHAI,AI,GI)
      IMPLICIT REAL*8 (A-H,O-Z)
      implicit integer*4 (i-n)
      PARAMETER (N30=4500)
      DIMENSION SR(N30,N30),SI(N30,N30),V(N30),
     1ALPHA(N30,3,3),G(N30,3,3),A(N30,3,3,3),
     1ALPHAR(N30,3,3),GR(N30,3,3),AR(N30,3,3,3),
     1ALPHAI(N30,3,3),GI(N30,3,3),AI(N30,3,3,3)
      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
      ALPHASR=0.0d0
      ALPHASI=0.0d0
      DO 2 IX=1,N
      ALPHASI=ALPHASI+SI(IX,IQ)*V(IX)
2     ALPHASR=ALPHASR+SR(IX,IQ)*V(IX)
      ALPHAI(IQ,I,J)=ALPHASI
1     ALPHAR(IQ,I,J)=ALPHASR

      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
      GSR=0.0d0
      GSI=0.0d0
      DO 21 IX=1,N
      GSI=GSI+SI(IX,IQ)*V(IX)
21    GSR=GSR+SR(IX,IQ)*V(IX)
      GR(IQ,I,J)=GSR
11    GI(IQ,I,J)=GSI

      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
      ASR=0.0D0
      ASI=0.0D0
      DO 22 IX=1,N
      ASI=ASI+SI(IX,IQ)*V(IX)
22    ASR=ASR+SR(IX,IQ)*V(IX)
      AR(IQ,I,J,K)=ASR
12    AI(IQ,I,J,K)=ASI

      RETURN
      END
c     ====================
      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end
c     ====================
      subroutine sset(s,i,i1)
      implicit none
      character*20 s
      integer*4 i,i1
      write(s,10)i
10    format(i20)
      do 1 i1=1,len(s)
1     if(s(i1:i1).ne.' ')return
      return
      end