PROGRAM CAR
C     this program reads roa ouptput
      IMPLICIT INTEGER*4 (I-N) 
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (NT0=500,MX3=3*NT0,BOHR=0.529177d0)
      DIMENSION al(MX3,3,3),G(MX3,3,3),A(MX3,3,3,3),
     1alv(MX3,3,3),ROAAFP(3,3),ROAGNP(3,3),ROAAP(3,3,3),r(3,NT0),
     1ROAAFVP(3,3),al0(3,3),g0(3,3),alv0(3,3),alp(MX3,3,3),
     1a0(3,3,3)
      CHARACTER*1 ok
      CHARACTER*4 s4
      logical lex
c
      open(4,file='NAT')
      read(4,*)nat
c     steo in A
      read(4,*)stepA
      step=stepA/0.529177d0
      read(4,*)ic
      close(4)
c
      do 4 l=1,3*nat
      do 4 i=1,3
      do 4 j=1,3
      al(l,i,j)=0.0d0
      alv(l,i,j)=0.0d0
      G(l,i,j)=0.0d0
      do 4 k=1,3
4     A(l,i,j,k)=0.0d0
c
      ii=-1
      do 1 i2=1,ic
      do 1 iat=1,nat
      is=1
      if(i2.eq.1.and.iat.eq.1)is=0
      do 1 iax=is,3
      ico=3*(iat-1)+iax
      ii=ii+1
c
      open(4,file='scr')
      write(4,40)ii
40    format(i4)
      rewind 4
      read(4,41)s4(1:4)
41    format(a4)
      close(4)
      write(6,67)i2,iat,iax,ii,ico
67    format(5i4)
c
      if(ii.le.9)i1=1
      if(ii.gt.9.and.ii.le.99)i1=2
      if(ii.gt.99.and.ii.le.999)i1=3
      if(ii.gt.999)i1=4
c
c     polarizability
      OPEN(7,FILE=s4(5-i1:4)//'AL.TEN')
      read(7,*)
      read(7,*)
      read(7,*)
      read(7,*)((ROAAFP(I,J),J=1,3),I=1,3)
      close(7)
c
c     velocity polarizability
      OPEN(7,FILE=s4(5-i1:4)//'ALV.TEN')
      read(7,*)
      read(7,*)
      read(7,*)
      read(7,*)((ROAAFVP(I,J),J=1,3),I=1,3)
c     J-velocity index
      close(7)
c
c     GTEN
      OPEN(7,FILE=s4(5-i1:4)//'G.TEN')
      read(7,*)
      read(7,*)
      read(7,*)
      read(7,*)((ROAGNP(I,J),J=1,3),I=1,3)
c     J-magnetic index
      close(7)
c
c     A
      OPEN(7,FILE=s4(5-i1:4)//'A.TEN')
      read(7,*)
      read(7,*)
      read(7,*)
      do 3 I=1,3
      read(7,*)ROAAP(I,1,1),ROAAP(I,2,2),ROAAP(I,3,3),
     1         ROAAP(I,1,2),ROAAP(I,1,3),ROAAP(I,2,3)
      ROAAP(I,2,1)=ROAAP(I,1,2)
      ROAAP(I,3,1)=ROAAP(I,1,3)
      ROAAP(I,3,2)=ROAAP(I,2,3)
3     continue
      close(7)
c
      if(ii.eq.0)then
       do 5 i=1,3
       do 5 j=1,3
       al0(i,j)=ROAAFP(i,j)
       alv0(i,j)=ROAAFVP(i,j)
       g0(i,j)=ROAGNP(i,j)
       do 5 k=1,3
5      a0(i,j,k)=ROAAP(i,j,k)
      else
       sign=1.0d0
       if(i2.eq.2)sign=-1.0d0
       do 9 i=1,3
       do 9 j=1,3
       al(ico,i,j)=al(ico,i,j)+sign*(ROAAFP(i,j)-al0(i,j))
       alv(ico,i,j)=alv(ico,i,j)+sign*(ROAAFVP(i,j)-alv0(i,j))
       G(ico,i,j)=G(ico,i,j)+sign*(ROAGNP(i,j)-g0(i,j))
       do 9 k=1,3
9      A(ico,i,j,k)=A(ico,i,j,k)+sign*(ROAAP(i,j,k)-a0(i,j,k))
      endif
1     continue
c
      bs=real(ic)*step
      do 6 ico=1,3*nat
      do 6 i=1,3
      do 6 j=1,3
      al(ico,i,j)=al(ico,i,j)/bs
      alv(ico,i,j)=alv(ico,i,j)/bs
      G(ico,i,j)=G(ico,i,j)/bs
      do 6 k=1,3
6     A(ico,i,j,k)=A(ico,i,j,k)/bs
c
      write(6,*)'Origin transformation with FILE.X (y/n) ?'
      read(5,'(a)')ok
      if(ok.eq.'Y'.or.ok.eq.'y')then
       inquire(file='FILE.X',exist=lex)
       if(lex)then
        open(34,file='FILE.X')
        read(34,*)
        read(34,*)natx
        if(natx.ne.nat)then
         write(6,*)'nat <> natx !'
         stop
        endif
        do 7 i=1,nat
        read(34,*)idum,(r(j,i),j=1,3)
        do 7 j=1,3
7       r(j,i)=r(j,i)/BOHR
        close(34)
        call ttai(MX3,al,alv,A,G,nat,4,r)
       else
        write(6,*)'FILE.X does not exist, no trafo done !'
       endif 
      endif
c
      write(6,*)'Improvement with ALPHA.TTT (y/n) ?'
      read(5,'(a)')ok
      if(ok.eq.'Y'.or.ok.eq.'y')then
       inquire(file='FILE.X',exist=lex)
       if(lex)then
        open(34,file='FILE.X')
        read(34,*)
        read(34,*)natx
        if(natx.ne.nat)then
         write(6,*)'nat <> natx !'
         stop
        endif
        do 8 i=1,nat
        read(34,*)idum,(r(j,i),j=1,3)
        do 8 j=1,3
8       r(j,i)=r(j,i)/BOHR
        close(34)
        inquire(file='ALPHA.TTT',exist=lex)
        if(lex)then
         open(34,file='ALPHA.TTT')
         read(34,*)
         read(34,*)nata
         read(34,*)
         read(34,*)
         if(nat.ne.nata)then
          write(6,*)'nat different !'
          stop
         endif
         do 10 i=1,3
         read(34,*)
         do 10 l=1,nat
         do 10 ix=1,3
         IIND=3*(l-1)+ix
10       read(34,*)idum,idum,(alp(IIND,i,j),j=1,3)
         call ttai(MX3,alp,al,A,G,nat,5,r)
         call WRITETTT(MX3,NAT,alp,A,G,al)
         WRITE(*,*)' Tensors written into CAR.TTT'
         stop
        else
         write(6,*)'ALPHA.TTT does not exist, no trafo done !'
        endif
       else
        write(6,*)'FILE.X does not exist, no trafo done !'
       endif 
      endif
c
      call WRITETTT(MX3,NAT,AL,A,G,alv)
      WRITE(*,*)' Tensors written into CAR.TTT'
      call sumr(MX3,NAT,AL,A,G,alv)
      stop
999   write(6,*)'Unexpected end of file'
      close(7)
      stop
      END
c     ==============================================
      SUBROUTINE WRITETTT(N30,NAT,ALPHA,A,G,alv)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(N30,3,3),G(N30,3,3),A(N30,3,3,3)
      DIMENSION ALv(N30,3,3)
c
      OPEN(2,FILE='CAR.TTT')
      WRITE(2,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(2,2002)I
2002  FORMAT(' Alpha(',I1,',J):')
      DO 1 L=1,NAT
      DO 1 IX=1,3
      IIND=3*(L-1)+IX
1     WRITE(2,2001)L,IX,(ALPHA(IIND,I,J),J=1,3)
2001  FORMAT(I5,1H ,I1,3F13.7)
      WRITE(2,2004)
2004  FORMAT(' The electric dipole magnetic dipole polarizability:',/,
     1       ' Atom/x    jx           jy          jz (magn. index -->)')
      DO 2 I=1,3
      WRITE(2,2003)I
2003  FORMAT(' G(',I1,',J):')
      DO 2 L=1,NAT
      DO 2 IX=1,3
      IIND=3*(L-1)+IX
c     (last index magnetic)
2     WRITE(2,2001)L,IX,(G(IIND,I,J),J=1,3)
      WRITE(2,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(2,2006)I,J
2006  FORMAT(' A(',I1,',',I1,',K):')
      DO 3 L=1,NAT
      DO 3 IX=1,3
      IIND=3*(L-1)+IX
3     WRITE(2,2001)L,IX,(A(IIND,I,J,K),K=1,3)
      WRITE(2,2007)
2007  FORMAT(' The polarizability ("velocity"):',/,
     1       ' Atom/x   vx          vy         vz ')
      DO 4 I=1,3
      WRITE(2,2002)I
      DO 4 L=1,NAT
      DO 4 IX=1,3
      IIND=3*(L-1)+IX
4     WRITE(2,2001)L,IX,(ALv(IIND,I,J),J=1,3)
      CLOSE(2)
      RETURN
      END
c     ============================================
      SUBROUTINE sumr(N30,NAT,ALPHA,A,G,alv)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(N30,3,3),G(N30,3,3),A(N30,3,3,3),
     1tem(3,3)
      DIMENSION ALv(N30,3,3)
c
      OPEN(2,FILE='SUM.TXT')
      WRITE(2,*)'translational sum rules'
      WRITE(2,*)'alpha'
      DO 1 IX=1,3
      WRITE(2,*)IX,' coord'
      WRITE(2,*)'first atom:'
      do 11 i=1,3
11    WRITE(2,2001)(ALPHA(IX,I,J),J=1,3)
2001  FORMAT(3F13.7)
      do 12 i=1,3
      do 12 j=1,3
      tem(i,j)=0.0d0
      DO 12 L=2,NAT
      IIND=3*(L-1)+IX
12    tem(i,j)=tem(i,j)-ALPHA(IIND,I,J)
      WRITE(2,*)'sum rule:'
      do 14 i=1,3
14    WRITE(2,2001)(tem(I,J),J=1,3)
1     continue
c
      WRITE(2,*)'G'
      DO 2 IX=1,3
      WRITE(2,*)IX,' coord'
      WRITE(2,*)'first atom:'
      do 21 i=1,3
21    WRITE(2,2001)(G(IX,I,J),J=1,3)
      do 22 i=1,3
      do 22 j=1,3
      tem(i,j)=0.0d0
      DO 22 L=2,NAT
      IIND=3*(L-1)+IX
22    tem(i,j)=tem(i,j)-G(IIND,I,J)
      WRITE(2,*)'sum rule:'
      do 24 i=1,3
24    WRITE(2,2001)(tem(I,J),J=1,3)
2     continue
c
      WRITE(2,*)'A'
      DO 3 IX=1,3
      DO 3 K=1,3
      WRITE(2,*)IX,' coord; ',K,' ind'
      WRITE(2,*)'first atom:'
      do 31 i=1,3
31    WRITE(2,2001)(A(IX,K,I,J),J=1,3)
      do 32 i=1,3
      do 32 j=1,3
      tem(i,j)=0.0d0
      DO 32 L=2,NAT
      IIND=3*(L-1)+IX
32    tem(i,j)=tem(i,j)-A(IIND,K,I,J)
      WRITE(2,*)'sum rule:'
      do 34 i=1,3
34    WRITE(2,2001)(tem(I,J),J=1,3)
3     continue
c
      WRITE(2,*)'alphaV'
      DO 4 IX=1,3
      WRITE(2,*)IX,' coord'
      WRITE(2,*)'first atom:'
      do 41 i=1,3
41    WRITE(2,2001)(ALV(IX,I,J),J=1,3)
      do 42 i=1,3
      do 42 j=1,3
      tem(i,j)=0.0d0
      DO 42 L=2,NAT
      IIND=3*(L-1)+IX
42    tem(i,j)=tem(i,j)-ALV(IIND,I,J)
      WRITE(2,*)'sum rule:'
      do 44 i=1,3
44    WRITE(2,2001)(tem(I,J),J=1,3)
4     continue
      CLOSE(2)
      RETURN
      END
C                                                                      
      FUNCTION EPS(I,J,K)
      IMPLICIT INTEGER*4 (I-N) 
      REAL*8 EPS
      INTEGER*4 I,J,K
      EPS=0.0d0
      IF (I.EQ.1.AND.J.EQ.2.AND.K.EQ.3)EPS= 1.0d0
      IF (I.EQ.1.AND.J.EQ.3.AND.K.EQ.2)EPS=-1.0d0
      IF (I.EQ.2.AND.J.EQ.3.AND.K.EQ.1)EPS= 1.0d0
      IF (I.EQ.2.AND.J.EQ.1.AND.K.EQ.3)EPS=-1.0d0
      IF (I.EQ.3.AND.J.EQ.1.AND.K.EQ.2)EPS= 1.0d0
      IF (I.EQ.3.AND.J.EQ.2.AND.K.EQ.1)EPS=-1.0d0
      RETURN                                                          
      END                                                             
c
      SUBROUTINE ttai(N,ALPHA,AV,A,G,NAT,ICO,R)
C
C     ICO=1  ... transform A and G into local origins
C     ICO=2  ... transform A and G into common origin
C     ICO=3  ... zero out A and G 
C     ICO=4  ... distributed origin transformation for G
C     ICO=5  ... improvement with better alpha, in AV is the worse one
C
C     R supposed to be passed in atomic units, G=G/w
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),R(3,NAT),AV(N,3,3)
C
      IF(ICO.EQ.1)SG  =  1.0D0
      IF(ICO.EQ.1)SG2 =  0.0D0
      IF(ICO.EQ.2)SG  = -1.0D0
      IF(ICO.EQ.2)SG2 =  0.0D0
      IF(ICO.EQ.4)SG  =  1.0D0
      IF(ICO.EQ.4)SG2 =  1.0D0
      IF(ICO.EQ.5)SG  =  1.0D0
      IF(ICO.EQ.5)SG2 =  1.0D0
C
      DO 1 IA=1,3
      DO 1 IB=1,3
      DO 1 L=1,NAT
      DO 1 IE=1,3
      K=3*(L-1)+IE
      if(ICO.EQ.3)then
       G(K,IA,IB)=0.0d0
      else
       DO 2 IG=1,3
       DO 2 ID=1,3
2      G(K,IA,IB)=G(K,IA,IB)
     1       +SG *0.5D0*EPS(IB,IG,ID)*R(IG,L)*AV(   K,IA,ID)
     2       -SG2*0.5D0*EPS(IB,IG,ID)*R(IG,L)*ALPHA(K,IA,ID)
c      last index of G magnetic
      endif
1     continue
c
      IF(ICO.EQ.1)WRITE(6,*)' G transformed into atomic origins with AV'
      IF(ICO.EQ.2)WRITE(6,*)' G transformed into common origin with AV'
      IF(ICO.EQ.3)WRITE(6,*)' G set to zero'
      IF(ICO.EQ.4)then
       WRITE(6,*)' DOG for G'
       return
      ENDIF
      IF(ICO.EQ.5)WRITE(6,*)' G inproved with ALPHA.TTT'
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
      K=3*(L-1)+IE
      IF(ICO.EQ.3)then
       A(K,IA,IB,IC)=0.0D0
      else
       S1=0.0D0
       S2=0.0D0
       IF(IB.EQ.IC)THEN
        DO 4 ID=1,3
        S2=S2+R(ID,L)*   AV(K,IA,ID)
4       S1=S1+R(ID,L)*ALPHA(K,IA,ID)
       ENDIF
       if(ICO.EQ.5)then
        A(K,IA,IB,IC)=A(K,IA,IB,IC)
     1  -(S1-1.5D0*(R(IB,L)*ALPHA(K,IA,IC)+R(IC,L)*ALPHA(K,IA,IB)))
     1  +(S2-1.5D0*(R(IB,L)*AV(   K,IA,IC)+R(IC,L)*AV(   K,IA,IB)))
       else
        A(K,IA,IB,IC)=A(K,IA,IB,IC)
     1  +SG *(S1-1.5D0*(R(IB,L)*ALPHA(K,IA,IC)+R(IC,L)*ALPHA(K,IA,IB)))
       endif
      endif
3     CONTINUE
      IF(ICO.EQ.1)WRITE(6,*)' A transformed into atomic origins'
      IF(ICO.EQ.2)WRITE(6,*)' A transformed into common origin'
      IF(ICO.EQ.3)WRITE(6,*)' A set to zero'
      IF(ICO.EQ.5)WRITE(6,*)' A inproved with ALPHA.TTT'
c
      RETURN
      END