PROGRAM NEW8
      IMPLICIT none
      integer*4 N,NAT,NQ
      real*8,allocatable::S(:,:),E(:),P(:,:,:),A(:,:,:),
     1PI(:,:),AI(:,:),ALPHA(:,:,:),G(:,:,:),AT(:,:,:,:)
C
      WRITE(6,6666)
6666  FORMAT(/,/,
     A'         * * * * * * * * * * * * * * * * * * * * * * *',/,
     1'       *                                          *   ',/,
     2'     *  MAGNETIC VIBRATIONAL CIRCULAR DICRHOISM *     ',/,
     3'       *                                          *   ',/,
     4'         * * * * * * * * * * * * * * * * * * * * * * *',/,
     5'                                                      ',/,
     6'                                                      ',/,
     7'                 Petr Bour 2017                       ',/,
     8'                                                      ',/,
     8'                Input:  F.INP                         ',/,
     8'                        FILE.TEN                      ',/,
     8'                        FILE.TTT                      ',/,
     8'                Output: MVCD.TAB                      ',/,
     8'                                                      ',/,
     8'                                                      ',/)
     
c     

      N=0
      call  READS(S,E,N,NAT,NQ,0)
      allocate(S(N,N),E(N),P(NAT,3,3),A(NAT,3,3),PI(NQ,3),AI(NQ,3))
      call  READS(S,E,N,NAT,NQ,1)
      call rdten(NAT,P,A)
      call tentr(NQ,N,NAT,P,A,S,PI,AI)
      call dogtab(NQ,PI,AI,E)
      allocate(ALPHA(N,3,3),G(N,3,3),AT(N,3,3,3))
      call READTTT(N,NAT,ALPHA,G,AT)
      CALL TRTEN(S,N,NQ,ALPHA,AT,G)
      CALL DORAMAN(N,NQ,ALPHA,G,AT,5320.0d0,50.0d0,1.0d6,E)
      call mvcdtab(N,NQ,PI,G,E)
      end

      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end
      
      SUBROUTINE READS(S,E,N,NAT,NQ,ic)
c     ic=0  .... dimensions only
c        1  .... whole matrix
      IMPLICIT none
      integer*4 NAT,NQ,ic,N,i,j,ix
      real*8 S(N,N),E(*),SFAC
      SFAC=0.0234280d0
C     1 a.u.= 2.1947E5 cm^-1
      OPEN(4,FILE='F.INP',STATUS='OLD')
      read(4,*)NQ,N,NAT
      if(ic.eq.0)then
       close(4)
       return
      endif
      do 1 i=1,NAT
1     read(4,*)
      read(4,*)
      DO 2 I=1,NAT
      DO 2 J=1,NQ
2     read(4,*)(s(3*(i-1)+ix,j),ix=1,2),(s(3*(i-1)+ix,j),ix=1,3)
      read(4,*)
      READ(4,4000)(E(I),I=1,NQ)
4000  FORMAT(6F11.6)
      DO 5 I=1,N
      DO 5 J=1,N
5     S(I,J)=S(I,J)*SFAC
      close(4)
      WRITE(*,*)' S-matrix read in'
      RETURN
      end

      subroutine rdten(NAT,P,A)
      implicit none
      integer*4 NAT,I,J,L,NOAT
      real*8 P(NAT,3,3),A(NAT,3,3)
      OPEN(15,FILE='FILE.TEN',STATUS='OLD')
      READ (15,*)NOAT
      if(NAT.NE.NOAT)call report('NAT<>NOAT')
C     READ DIPOLE DERIVATIVES = TENSOR P(A,B):
      DO 10 L=1,NOAT
      DO 10 J=1,3
10    READ (15,*) (P(L,J,I),I=1,3)
C     READ AAT:
      DO 220 L=1,NOAT
      DO 220 J=1,3
 220  READ (15,*) (A(L,J,I),I=1,3)
      CLOSE(15)
      WRITE(*,*)' Dipole derivatives read in'
      return
      end

      subroutine tentr(NQ,N,NAT,P,A,S,PI,AI)
      implicit none
      integer*4 N,NAT,I,L,LJ,J,B,NQ
      real*8 P(NAT,3,3),A(NAT,3,3),S(N,N),PI(NQ,3),AI(NQ,3)
      DO 1 I=1,NQ
      DO 1 B=1,3
      PI(I,B)=0.0d0
      AI(I,B)=0.0d0
      DO 1 L=1,NAT
      DO 1 J=1,3
      LJ=3*(L-1)+J
      PI(I,B)=PI(I,B)+P(L,J,B)*S(LJ,I)
1     AI(I,B)=AI(I,B)+A(L,J,B)*S(LJ,I)
      write(6,*)'APT,AAT transformed into normal modes'
      return
      end

      subroutine dogtab(NQ,P,A,E)
      implicit none
      integer*4 NQ,I
      real*8 P(NQ,3),A(NQ,3),E(*),D,R,WLIM,FD,W,ABSW,WR,FC,pi,po,
     1SFAC
      OPEN(18,FILE='DOG.TAB')
      WRITE (18,501)
  501 FORMAT (/,'          FREQ.     ',
     1'  D [D^2]     ','  R [D^-2]DO         R - CO      I[km/mol]',
     2/,'-----------------------------------------------------------')
      WLIM=1.0d-6
      FD =SQRT(388.9219d0)
      FC=131.14D-8
      pi=4.0d0*atan(1.0d0)
      po=dsqrt(2.0d0*pi)
      SFAC=0.0234280d0
      DO 1 I=1,NQ
      D=P(I,1)*P(I,1)+P(I,2)*P(I,2)+P(I,3)*P(I,3)
      R=P(I,1)*A(I,1)+P(I,2)*A(I,2)+P(I,3)*A(I,3)
      W=E(I)
      ABSW=DABS(W)
      WR=DSQRT(ABSW)
      IF(ABSW.GT.WLIM)D=D/WR**2*FD**2/SFAC**2
      R=R*FD*FC/SFAC**2
 1    WRITE (18,523) I,W,D,R,R,po*W*D
523   FORMAT(I5,f16.7,4G15.6)
      WRITE(18,1819)
1819  FORMAT('-------------------------------------------------')
      CLOSE(18)
      write(6,*)'DOG.TAB written'
      return
      end

      SUBROUTINE READTTT(N,NAT,ALPHA,G,A)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION ALPHA(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
      IND=3*(L-1)+IX
1     READ(2,*)(ALPHA(IND,I,J),J=1,2),(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,*)(G(IND,I,J),J=1,2),(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,*)(A(IND,I,J,K),K=1,2),(A(IND,I,J,K),K=1,3)
      CLOSE(2)
      WRITE(6,*)NAT,' atoms, tensors read in from FILE.TTT'
      RETURN
      END
c     ============================================================
      SUBROUTINE TRTEN(S,N,NQ,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,NQ
      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,NQ
      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,NQ
      AS=0.0D0
      DO 22 IX=1,N
22    AS=AS+S(IX,IQ)*V(IX)
12    A(IQ,I,J,K)=AS
      RETURN
      END
c     ============================================================
      SUBROUTINE DORAMAN(N,NQ,ALPHA,G,A,EXCA,wmin,wmax,E)
c
c     supposedly G/W is supplied
c
      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),E(*)
      CM=219470.0d0
      OPEN(4,FILE='ROA.TAB')
      WRITE(4,3304)
3304  FORMAT(
     1'MODE   FREQ    Ramx (180)  Ramy    Ramtot   CID90',
     2'     CID-X   CID180  DEP   ROA180 RamanDCPI',/,
     3'       cm-1       A^4/AMU                    x10^4',
     4'     x10^4     x10^4      A^5/AMUx1e4 au.10^4')
      write(4,3002)
3002  FORMAT(80(1H-))
c
      DO 1 IQ=1,NQ
      IF(E(IQ).GT.wmin.and.E(IQ).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
       S180=7.0d0*SAL1+SAL0
       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
       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      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      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)
       WRITE(4,3000)IQ,E(IQ),YDX,YDY,ram1,CID2,DX,CID1,P1,roa1,doc,roa3
3000   FORMAT(I5,f9.2,6g12.3,f6.3,4g12.4)
c
      endif
1     CONTINUE
      WRITE(4,3002)
      CLOSE(4)
      WRITE(6,*)'ROA.TAB written'
      RETURN
      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 mvcdtab(N,NQ,P,G,E)
      implicit none
      integer*4 N,NQ,I,A,B,C
      real*8 P(NQ,3),G(N,3,3),E(*),D,R,WLIM,FD,W,ABSW,WR,FC,pi,po,
     1SFAC,EPS
      OPEN(18,FILE='MVCD.TAB')
      WRITE (18,501)
  501 FORMAT (/,'          FREQ.     ',
     1'  D [D^2]     ','  R [D^-2]DO         R - CO      I[km/mol]',
     2/,'-----------------------------------------------------------')
      WLIM=1.0d-6
      FD =SQRT(388.9219d0)
      FC=131.14D-8
      pi=4.0d0*atan(1.0d0)
      po=dsqrt(2.0d0*pi)
      SFAC=0.0234280d0
      DO 1 I=1,NQ
      D=P(I,1)*P(I,1)+P(I,2)*P(I,2)+P(I,3)*P(I,3)
      R=0.0d0
      do 2 A=1,3
      do 2 B=1,3
      do 2 C=1,3
2     R=R+EPS(A,B,C)*P(I,A)*G(I,B,C)
      W=E(I)
      ABSW=DABS(W)
      WR=DSQRT(ABSW)
      IF(ABSW.GT.WLIM)then
       D=D/WR**2*FD**2/SFAC**2
       R=R*FD*FC/SFAC**2/W
      endif
 1    WRITE (18,523) I,W,D,R,R,po*W*D
523   FORMAT(I5,f16.7,4G15.6)
      WRITE(18,1819)
1819  FORMAT('-------------------------------------------------')
      CLOSE(18)
      write(6,*)'MVCD.TAB written'
      return
      end