PROGRAM GAR
C     this program reads the archive of the Gaussian 94 output
C     version adapted for f77 on RS6000, from the garsgi.f that
C     ran in Calgary, on Silicon Graphics Irix; 2-12-1996; 
C     Petr Bour, UOCHB, Praha
      IMPLICIT INTEGER*4 (I-N) 
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (NT0=500,MX3=3*NT0)
      DIMENSION F(MX3,MX3),P(NT0,3,3),X(3,NT0),ALPHA(3,3,MX3),
     1AAT(NT0,3,3),GTENS(3,3,MX3),ATENS(3,3,3,MX3),iaa(6),jaa(6)
      CHARACTER*80 filename
      CHARACTER*1 BUFF(140),NUMBER(80),L,LS,LE
      CHARACTER*2 AT(NT0)
      CHARACTER*5 SSTR
      CHARACTER*2 LG,L2
c jk-start
      parameter (MFR=9) !max number of frequencies
      dimension xg09(6,3*MX3),alpha09(MFR,3,3,MX3),gtens09(MFR,3,3,MX3),
     1 atens09(MFR,3,3,3,MX3),APOL(9),GPOL(9),AAPOL(27),ORT(9)
      character*90 radek
      character*1 fchar
c jk-end
      COMMON/LINE/I
      LOGICAL LVCD,LGEO,LARCH,auto,lgpol,lapol
      character*2 atsy(89)
      data atsy/'H ','He','Li','Be','B ','C ','N ','O ','F ','Ne',
     3'Na','Mg','Al','Si','P ','S ','Cl','Ar',
     4'K ','Ca','Sc','Ti','V ','Cr','Mn','Fe','Co','Ni','Cu','Zn',
     4          'Ga','Ge','As','Se','Br','Kr',
     5'Rb','Sr','Y ','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd',
     5          'In','Sn','Sb','Te','I ','Xe',
     6'Cs','Ba','La',
     6               'Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho',
     6               'Er','Tm','Yb','Lu',
     6'Hf','Ta','W ','Re','Os','Ir','Pt','Au','Hg',
     6          'Tl','Pb','Bi','Po','At','Rn',
     7'Fr','Ra','Ac'/
      data iaa/1,2,3,1,1,2/
      data jaa/1,2,3,2,3,3/
c
      lapol=.false.
      lgpol=.false.
      OPEN(8,FILE='scr',FORM='FORMATTED')
c
      inquire(file='AUTO',exist=auto)
      if(auto)then
       open(7,file='FILE.X')
       read(7,*) 
       read(7,*)NAT
       close(7)
      else
       WRITE(*,*)' How many atoms ?'
       READ(*,*)NAT
      endif
      N=3*NAT
c jk-start
      do ifr=1,MFR
        do i=1,N
          do ix=1,3
            do iy=1,3
              alpha09(ifr,ix,iy,i)=0.d0
              gtens09(ifr,ix,iy,i)=0.d0
              do iz=1,3
                atens09(ifr,ix,iy,iz,i)=0.d0
              enddo
            enddo
          enddo
        enddo          
      enddo
c jk-end
      do 101 i=1,N
      do 101 ix=1,3
      do 101 iy=1,3
      ALPHA(ix,iy,i)=0.0d0
      AAT(i,ix,iy)=0.0d0
      P(i,ix,iy)=0.0d0
      GTENS(ix,iy,i)=0.0d0
      do 101 iz=1,3
101   ATENS(ix,iy,iz,i)=0.0d0
      NB=70
      LVCD=.FALSE.
      LGEO=.FALSE.
      LS='\\'
      LE='='
      LG='%%'
      L2='\\\\'
      if(auto)then
       FILENAME='FRE.OUT'
      else
       WRITE(*,*)
     1 ' Geometry looked for after string ',LG,' (only for G94)'
       WRITE(*,*)' Filename ? '
       READ(*,'(A)')filename
      endif
      OPEN(7,FILE=FILENAME,STATUS='OLD',FORM='FORMATTED')
C
c jk-start
8000  read(7,'(a)',end=8999,err=8999)radek
8001  format(8x,5e14.6)
      if(radek(23:45).eq.'Alpha(-w,w) derivatives')then
        read(radek(59:59),'(i1)')ifr
        read(radek(62:71),'(f10.6)')fr
        write(fchar,'(i1)')ifr
c       read tensor
        write(*,*)'G09 alpha-tensor derivatives found'
        write(*,'(a,i1,a,f10.6,a)')'  frequency ',ifr,' = ',fr,' a.u.'
        ii1=0
8002    read(7,*)radek
        do 9 ix=1,3
9       read(7,8001)(xg09(ix,ii2),ii2=ii1+1,min(ii1+5,9*nat))
        ii1=ii1+5
        if(ii1.lt.9*nat) goto 8002
        do 8 ix=1,3
        do 8 iy=1,3
        do 8 iat=1,n
8       alpha09(ifr,iy,ix,iat)=xg09(ix,3*(iat-1)+iy)
c       write FILE.TTT
        CALL WRITETTT09(NAT,ALPHA09,GTENS09,ATENS09,MX3,MFR,ifr,fr)
        WRITE(*,*)' ALPHA written'
        goto 8000
      endif
      if(radek(23:60).eq.'FD Optical Rotation Tensor derivatives')then
        read(radek(74:74),'(i1)')ifr
        read(radek(77:86),'(f10.6)')fr
        write(fchar,'(i1)')ifr
c       read tensor
        write(*,*)'G09 G-tensor derivatives found'
        write(*,'(a,i1,a,f10.6,a)')'  frequency ',ifr,' = ',fr,' a.u.'
        ii1=0
8003    read(7,*)radek
        do 7 ix=1,3
7       read(7,8001)(xg09(ix,ii2),ii2=ii1+1,min(ii1+5,9*nat))
        ii1=ii1+5
        if(ii1.lt.9*nat) goto 8003
        do 6 ii1=1,3
        do 6 ii2=1,3*n
6       gtens09(ifr,mod(ii2-1,3)+1,ii1,(ii2+2)/3)=xg09(ii1,ii2)
c       write FILE.TTT
        CALL WRITETTT09(NAT,ALPHA09,GTENS09,ATENS09,MX3,MFR,ifr,fr)
        WRITE(*,*)' GTENS written'
        goto 8000
      endif
      if(radek(1:44).eq.
     1' Property number 4 -- D-Q polarizability fre')then
       read(radek(51:53),*)ifr
       read(7,*)
       do 13 ii1=1,6
       ix=iaa(ii1)
       iy=jaa(ii1)
       read(7,*)ii2,(AAPOL(ix+3*(iy-1)+9*(iz-1)),iz=1,3)
       do 13 iz=1,3
13     AAPOL(iy+3*(ix-1)+9*(iz-1))=AAPOL(ix+3*(iy-1)+9*(iz-1))
       call WRITEPOL(AAPOL,ifr,'A.TTT')
       if(lapol.and.lgpol)then
c       optical rotation tensor
        do 14 i=1,3
        do 14 j=1,3
        sum=gpol(i+3*(j-1))+gpol(j+3*(i-1))
        do 12 ig=1,3
        do 12 id=1,3
12      sum=sum-(EPS(i,ig,id)*AAPOL(id+3*(j-1)+9*(ig-1))
     1          +EPS(j,ig,id)*AAPOL(id+3*(i-1)+9*(ig-1)))/3.0d0
14      ORT(i+3*(j-1))=sum/2.0d0
        call WRITEPOL(ORT,ifr,'ORT.TTT')
       endif
      endif
      if(radek(1:52).eq.
     1' Property number 2 -- FD Optical Rotation Tensor fre')then
       lgpol=.true.
       read(radek(59:61),*)ifr
       read(7,*)
       do 11 ix=1,3
11     read(7,*)ii1,(GPOL(ix+3*(iy-1)),iy=1,3)
       call WRITEPOL(GPOL,ifr,'GP.TTT')
      endif
      if(radek(1:37).eq.' Property number 1 -- Alpha(-w,w) fre')then
       lapol=.true.
       read(radek(44:46),*)ifr
       read(7,*)
       do 2 ix=1,3
2      read(7,*)ii1,(APOL(ix+3*(iy-1)),iy=1,3)
       call WRITEPOL(APOL,ifr,'POL.TTT')
      endif
      if(radek(23:52).eq.'D-Q polarizability derivatives')then
        read(radek(66:66),'(i1)')ifr
        read(radek(69:78),'(f10.6)')fr
        write(fchar,'(i1)')ifr
c       read tensor
        write(*,*)'G09 A-tensor derivatives found'
        write(*,'(a,i1,a,f10.6,a)')'  frequency ',ifr,' = ',fr,' a.u.'
        ii1=0
8004    read(7,*)radek
        do 15 ix=1,6
15      read(7,8001)(xg09(ix,ii2),ii2=ii1+1,min(ii1+5,9*nat))
        ii1=ii1+5
        if(ii1.lt.9*nat) goto 8004
        do 16 ii1=1,6
        do 16 ii2=1,3*n
        ix=iaa(ii1)
        iy=jaa(ii1)
        iz=mod(ii2-1,3)+1
        iat=(ii2+2)/3
        atens09(ifr,ix,iy,iz,iat)=xg09(ii1,ii2)
16      atens09(ifr,iy,ix,iz,iat)=xg09(ii1,ii2)
c       write FILE.TTT
        CALL WRITETTT09(NAT,ALPHA09,GTENS09,ATENS09,MX3,MFR,ifr,fr)
        WRITE(*,*)' ATENS written'
        write(*,*)
        goto 8000
      endif
      goto 8000
8999  rewind(7)
c jk-end
      LARCH=.FALSE.
5444  READ(7,7000,END=8888)(BUFF(I),I=1,NB)
      READ(7,7000,END=8888)(BUFF(I),I=NB+1,2*NB)
7000  FORMAT(1X,70A1)
      BACKSPACE 7
      DO 51 I=1,2*NB-5
      DO 52 J=1,5
52    SSTR(J:J)=BUFF(I+J-1) 
      IF(SSTR(1:2).EQ.L2)LARCH=.TRUE.
c     IF(SSTR.EQ.'leDer')THEN
      IF(SSTR.EQ.'leDer'.and.BUFF(I-3).eq.'i')THEN
       IF(.NOT.LARCH)GOTO 5444
       WRITE(*,*)' Atomic polar tensor found'
       GOTO 6555
      ENDIF
      IF(SSTR(3:5).EQ.'AAT')THEN
       IF(.NOT.LARCH)GOTO 5444
       WRITE(*,*)' Atomic axial tensor found'
       LVCD=.TRUE.
       GOTO 9555
      ENDIF
      IF(SSTR.EQ.'arDer')THEN
       IF(.NOT.LARCH)GOTO 5444
       WRITE(*,*)' Polarization Derivatives found'
       GOTO 8555
      ENDIF
      IF(SSTR.EQ.'olDer')THEN
       IF(.NOT.LARCH)GOTO 5444
       WRITE(*,*)' A-tensor Derivatives found'
       GOTO 11555
      ENDIF
      IF(SSTR.EQ.'otDer'.or.SSTR.EQ.'nsDer')THEN
       IF(.NOT.LARCH)GOTO 5444
       WRITE(*,*)' G-tensor Derivatives found'
       GOTO 10555
      ENDIF
C     NImag for g94, NIMAG for g92
      IF(SSTR.EQ.'NImag'.OR.SSTR.EQ.'NIMAG')THEN
       IF(.NOT.LARCH)GOTO 5444
       WRITE(*,*)' Second derivatives found'
       GOTO 5555
      ENDIF
      IF((.NOT.LGEO).AND.SSTR(1:2).EQ.LG)THEN
       IF(.NOT.LARCH)GOTO 5444
       WRITE(*,*)' Geometry found'
       LGEO=.TRUE.
       GOTO 7555
      ENDIF
51    CONTINUE
      GOTO 5444
C
7555  I=I+1
C     I .. index of the letter to be gotten
      IF(I.GT.NB)THEN
       DO 74 K=1,NB
74     BUFF(K)=BUFF(K+NB)
       I=I-NB
       READ(7,*)
      ENDIF
      DO 73 IA=1,NAT
      AT(IA)(1:1)=' '
      AT(IA)(2:2)=' '
C
C     Get element:
      I=I+1
      CALL GETLET(L,BUFF,NB)
      AT(IA)(1:1)=L
      I=I+1
      CALL GETLET(L,BUFF,NB)
      if(L.ne.',')then
       AT(IA)(2:2)=L
       I=I+1
       CALL GETLET(L,BUFF,NB)
      endif
C
      DO 773 IX=1,3
C     READ X(IX,IA):
      DO 75 J=1,20
75    NUMBER(J)=' '
      IN=0
7558  IN=IN+1
      I=I+1
      CALL GETLET(L,BUFF,NB)
      NUMBER(IN)=L
      IF(L.NE.','.AND.L.NE.LS)GOTO 7558
      NUMBER(IN)=' '
      REWIND 8
      WRITE(8,*)(NUMBER(J),J=1,20)
      REWIND 8
      READ(8,*)FU
      X(IX,IA)=FU
773   CONTINUE
73    CONTINUE
c
      OPEN(15,FILE='FILE.XYZ')
      WRITE(15,*)' FILE.XYZ from Gaussian output'
      DO 77 IA=1,NAT
      K=0
      IF(AT(IA).EQ.'C ')K=1
      IF(AT(IA).EQ.'O ')K=2
      IF(AT(IA).EQ.'N ')K=3
      IF(AT(IA).EQ.'H ')K=4
      IF(AT(IA).EQ.'S ')K=5
      IF(AT(IA).EQ.'P ')K=9
      IF(AT(IA).EQ.'Fe')K=26
      IF(K.EQ.0)WRITE(*,*)' Unknown atom ',AT(IA)
77    WRITE(15,1500)K,(X(IX,IA),IX=1,3),0,0,0,0,0.0
1500  FORMAT(I4,3F15.8,4I2,F4.1)
      WRITE(15,1500)9999,0.0,0.0,0.0,0,0,0,0,0.0
      CLOSE(15)
      WRITE(*,*)' X written into FILE.XYZ'
c
      OPEN(15,FILE='FILE.X')
      WRITE(15,*)' FILE.X from Gaussian output'
      WRITE(15,*)NAT
      DO 777 IA=1,NAT
      K=0
      do 1 ii=1,89
1     IF(AT(IA).EQ.atsy(ii))K=ii
      IF(K.EQ.0)then
       WRITE(*,*)' Unknown atom ',AT(IA)
       stop
      endif
777   WRITE(15,1501)K,(X(IX,IA),IX=1,3)
1501  FORMAT(I4,3F15.8,' 0 0 0 0 0 0 0 0.0')
      CLOSE(15)
      WRITE(*,*)' X written into FILE.X'
c
      BACKSPACE 7
      GOTO 5444
C
C
C     G' derivatives
10555 I=I+1
C     I .. index of the letter to be gotten
      IF(I.GT.NB)THEN
       DO 786 K=1,NB
786    BUFF(K)=BUFF(K+NB)
       I=I-NB
       READ(7,*)
      ENDIF
85561 I=I+1
      CALL GETLET(L,BUFF,NB)
      IF(L.NE.LE)GOTO 85561
      DO 801 IA=1,N
      DO 801 IX=1,3
      DO 801 IY=1,3
C
C     READ G'(IX,IY,IA)
      DO 791 J=1,20
791   NUMBER(J)=' '
      IN=0
85581 IN=IN+1
      I=I+1
      CALL GETLET(L,BUFF,NB)
      NUMBER(IN)=L
      IF(L.NE.','.AND.L.NE.'\\')GOTO 85581
      NUMBER(IN)=' '
      REWIND 8
      WRITE(8,*)(NUMBER(J),J=1,20)
      REWIND 8
      READ(8,*)FU
801   GTENS(IX,IY,IA)=FU
      CALL WRITETTT(NAT,ALPHA,GTENS,ATENS,MX3)
      WRITE(*,*)' GTENS written into FILE.TTT'
      BACKSPACE 7
      GOTO 5444
C
C
C     A (dipole-quadrupole) derivatives
11555 I=I+1
C     I .. index of the letter to be gotten
      IF(I.GT.NB)THEN
       DO 787 K=1,NB
787    BUFF(K)=BUFF(K+NB)
       I=I-NB
       READ(7,*)
      ENDIF
85562 I=I+1
      CALL GETLET(L,BUFF,NB)
      IF(L.NE.LE)GOTO 85562
      DO 802 IA=1,N
      DO 802 IZ=1,3
      DO 802 I6=1,6
      IX=iaa(I6)
      IY=jaa(I6)
C
C     READ A(IX,IY,IZ,IA)
      DO 792 J=1,20
792   NUMBER(J)=' '
      IN=0
85582 IN=IN+1
      I=I+1
      CALL GETLET(L,BUFF,NB)
      NUMBER(IN)=L
      IF(L.NE.','.AND.L.NE.'\\')GOTO 85582
      NUMBER(IN)=' '
      REWIND 8
      WRITE(8,*)(NUMBER(J),J=1,20)
      REWIND 8
      READ(8,*)FU
      ATENS(IY,IX,IZ,IA)=FU
802   ATENS(IX,IY,IZ,IA)=FU
      CALL WRITETTT(NAT,ALPHA,GTENS,ATENS,MX3)
      WRITE(*,*)' ATENS written into FILE.TTT'
      BACKSPACE 7
      GOTO 5444
C
C
C     Polarization derivatives
8555  I=I+1
C     I .. index of the letter to be gotten
      IF(I.GT.NB)THEN
       DO 78 K=1,NB
78     BUFF(K)=BUFF(K+NB)
       I=I-NB
       READ(7,*)
      ENDIF
8556  I=I+1
      CALL GETLET(L,BUFF,NB)
      IF(L.NE.LE)GOTO 8556
      DO 80 IA=1,N
      DO 80 IX=1,3
      DO 80 IY=1,IX
C
C     READ ALPHA(IX,IY,IA)
      DO 79 J=1,20
79    NUMBER(J)=' '
      IN=0
8558  IN=IN+1
      I=I+1
      CALL GETLET(L,BUFF,NB)
      NUMBER(IN)=L
      IF(L.NE.','.AND.L.NE.'\\')GOTO 8558
      NUMBER(IN)=' '
      REWIND 8
      WRITE(8,*)(NUMBER(J),J=1,20)
      REWIND 8
      READ(8,*)FU
      ALPHA(IY,IX,IA)=FU
80    ALPHA(IX,IY,IA)=FU
      CALL WRITETTT(NAT,ALPHA,GTENS,ATENS,MX3)
      WRITE(*,*)' ALPHA written into FILE.TTT'
      BACKSPACE 7
      GOTO 5444
C
C     Dipole derivatives     
6555  I=I+1
C     I .. index of the letter to be gotten
      IF(I.GT.NB)THEN
       DO 54 K=1,NB
54     BUFF(K)=BUFF(K+NB)
       I=I-NB
       READ(7,*)
      ENDIF
6556  I=I+1
      CALL GETLET(L,BUFF,NB)
      IF(L.NE.LE)GOTO 6556
      DO 53 IA=1,NAT
      DO 53 IX=1,3
      DO 53 IY=1,3
C
C     READ P(IA,IX,IY):
      DO 55 J=1,20
55    NUMBER(J)=' '
      IN=0
6558  IN=IN+1
      I=I+1
      CALL GETLET(L,BUFF,NB)
      NUMBER(IN)=L
      IF(L.NE.','.AND.L.NE.'\\')GOTO 6558
      NUMBER(IN)=' '
      REWIND 8
      WRITE(8,*)(NUMBER(J),J=1,20)
      REWIND 8
      READ(8,*)FU
53    P(IA,IX,IY)=FU
      CALL WRITETEN(NT0,P,AAT,LVCD,NAT)
      BACKSPACE 7
      GOTO 5444
C
C
C     Magnetic dipole derivatives     
9555  I=I+1
C     I .. index of the letter to be gotten
      IF(I.GT.NB)THEN
       DO 56 K=1,NB
56     BUFF(K)=BUFF(K+NB)
       I=I-NB
       READ(7,*)
      ENDIF
9556  I=I+1
      CALL GETLET(L,BUFF,NB)
      IF(L.NE.LE)GOTO 9556
      DO 58 IA=1,NAT
      DO 58 IX=1,3
      DO 58 IY=1,3
C
C     READ AAT(IA,IX,IY):
      DO 59 J=1,20
59    NUMBER(J)=' '
      IN=0
9558  IN=IN+1
      I=I+1
      CALL GETLET(L,BUFF,NB)
      NUMBER(IN)=L
      IF(L.NE.','.AND.L.NE.'\\')GOTO 9558
      NUMBER(IN)=' '
      REWIND 8
      WRITE(8,*)(NUMBER(J),J=1,20)
      REWIND 8
      READ(8,*)FU
58    AAT(IA,IX,IY)=FU
      CALL WRITETEN(NT0,P,AAT,LVCD,NAT)
      WRITE(*,*)' AAT written into FILE.TEN'
      BACKSPACE 7
      GOTO 5444
C
C     Force Field:
5555  I=I+1
C     I .. index of the letter to be gotten
      IF(I.GT.70)THEN
       DO 4 K=1,70
4      BUFF(K)=BUFF(K+70)
       I=I-70
       READ(7,*)
      ENDIF
C     NB .. line length
5556  I=I+1
      CALL GETLET(L,BUFF,NB)
      IF(L.NE.LS)GOTO 5556
5557  I=I+1
      CALL GETLET(L,BUFF,NB)
      IF(L.NE.LS)GOTO 5557
c     write(6,*)'numbers',I
      DO 3 IX=1,N
      DO 3 IY=1,IX
C
C     READ F(IX,IY):
      DO 5 J=1,20
5     NUMBER(J)=' '
      IN=0
5558  IN=IN+1
      I=I+1
      CALL GETLET(L,BUFF,NB)
      NUMBER(IN)=L
      IF(L.NE.','.AND.L.NE.'\\')GOTO 5558
      NUMBER(IN)=' '
      REWIND 8
      WRITE(8,*)(NUMBER(J),J=1,20)
      REWIND 8
      READ(8,*)FU
C     WRITE(9,*)IX,IY,FU
      F(IX,IY)=FU
3     F(IY,IX)=FU
      OPEN(20,FILE='FILE.FC')
      CALL WRITEFF(MX3,N,F)
      CLOSE(20)
      WRITE(*,*)' FF written into FILE.FC'
      BACKSPACE 7
      GOTO 5444
C
C
8888  CLOSE(7)
      STOP
      END
 
      SUBROUTINE WRITEFF(MX3,N,FCAR)
      IMPLICIT INTEGER*4 (I-N) 
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION FCAR(MX3,N)
C     CONST=4.359828/0.5291772**2
      CONST=1.0d0
      DO 6 I=1,N
      DO 6 J=1,N
6     FCAR(I,J)=FCAR(I,J)/CONST
      N1=1
1     N3=N1+4
      IF(N3.GT.N)N3=N
      DO 130 LN=N1,N
130   WRITE(20,17)LN,(FCAR(LN,J),J=N1,MIN(LN,N3))
      N1=N1+5
      IF(N3.LT.N)GOTO 1
17    FORMAT(I4,5D14.6)
      RETURN
      END

      SUBROUTINE GETLET(L,BUFF,NB)
      IMPLICIT INTEGER*4 (I-N) 
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER*1 BUFF(140)
      CHARACTER*1 L
      COMMON/LINE/I
      IF(I.GT.NB)THEN
       IF(I.GT.NB+1)THEN
        WRITE(*,*)' Cannot read letter ',I
c       STOP
       ENDIF
       I=1
       READ(7,7000)(BUFF(J),J=1,70)
7000   FORMAT(1X,70A1)
      ENDIF
      L=BUFF(I)
      RETURN
      END
        
      SUBROUTINE WRITETEN(N0,P,A,LVCD,NAT)
      IMPLICIT INTEGER*4 (I-N) 
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION P(N0,3,3),A(N0,3,3)
      LOGICAL LVCD
c     BOHR=0.52917705993d0
      Z0=0.0d0
      OPEN(15,FILE='FILE.TEN')
      WRITE(15,1500) NAT,NAT-6,0
1500  FORMAT(3I5)
      DO 10 L=1,NAT
      DO 10 J=1,3
10    WRITE(15,1501) (P(L,J,I),I=1,3),L
1501  FORMAT(3F14.8,I5)
      IF(.NOT.LVCD)THEN
       DO 220 L=1,NAT                                                 
       DO 220 J=1,3                                                    
c      write(6,*)L,J,nat
220    WRITE(15,1501) (Z0,I=1,3),L
      ELSE
       DO 221 L=1,NAT                                                 
       DO 221 J=1,3                                                    
c      write(6,*)L,J
221    WRITE(15,1501) (A(L,J,I),I=1,3),L
      ENDIF
      DO 230 L=1,NAT                                                 
      DO 230 J=1,3                                                    
230   WRITE(15,1501) (Z0,I=1,3),L
      DO 100 L=1,NAT
      DO 100 J=1,3
100   WRITE(15,1501) (P(L,J,I),I=1,3),L
      WRITE(*,*)' Dipole derivatives written into FILE.TEN'
      CLOSE(15)
      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 
c     SUBROUTINE VCDD0(VCD,VEL,DIP,C,NAT)
C     STOLEN FROM CADPAC, INDICES ARE DIFFERENT !!, VEL WITH NUCLEI
c     IMPLICIT REAL*8(A-H,O-Z)
c     IMPLICIT INTEGER*4 (I-N)
c     INTEGER*4 ALPHA,BETA,GAMMA,DELTA,LAMBDA                           
c     PARAMETER (MX=100)
c     DIMENSION VCD(MX,3,3),VEL(MX,3,3),DIP(MX,3,3),                  
c    1C(3,NAT)
c     DO 90 BETA =1,3                                                 
c     DO 90 GAMMA=1,3                                                 
c     DO 90 DELTA=1,3                                                 
c     SKEW=0.25D0*EPS(BETA,GAMMA,DELTA)                               
c     IF (DABS(SKEW).GT.1.0D-10)THEN
c     DO 80 ALPHA=1,3                                                 
c     DO 80 LAMBDA=1,NAT                                              
c80    VCD(LAMBDA,ALPHA,BETA)=VCD(LAMBDA,ALPHA,BETA)                   
c    1 -SKEW*C(GAMMA,LAMBDA)                                          
c    2 *(VEL(LAMBDA,ALPHA,DELTA)-DIP(LAMBDA,ALPHA,DELTA))             
c     ENDIF                                                           
c90    CONTINUE
c     RETURN
c     END
C     =======================================
      SUBROUTINE WRITETTT(NAT,ALPHA,GTENS,ATENS,MX3)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(3,3,MX3),GTENS(3,3,MX3),ATENS(3,3,3,MX3)
      OPEN(2,FILE='FILE.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(I,J,IIND),J=1,3)
2001  FORMAT(I5,1H ,I1,3g15.7)
      WRITE(2,2004)
2004  FORMAT(' The electric dipole magnetic dipole polarizability:',/,
     1       ' Atom/x    jx(Bx)       jy(By)       jz(Bz)')
      DO 2 I=1,3
      WRITE(2,2003)I
2003  FORMAT(' G(',I1,',J):')
      DO 2 L=1,NAT
      DO 2 K=1,3
2     WRITE(2,2001)L,K,(GTENS(I,J,3*(L-1)+K),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 M=1,3
3     WRITE(2,2007)L,M,(ATENS(K,J,I,3*(L-1)+M)*3.d0/2.d0,K=1,3),L,M,I,J
2007  FORMAT(I5,1H ,I1,3F15.7,' ',4i3)
      write(2,*)
      write(2,*)'dummy alpha v:'
      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,(ALPHA(I,J,IIND),J=1,3)
      CLOSE(2)
      RETURN
      END
c jk-start
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      SUBROUTINE WRITETTT09(NAT,ALPHA09,GTENS09,ATENS09,MX3,MFR,ifr,fr)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA09(MFR,3,3,MX3),GTENS09(MFR,3,3,MX3),
     1 ATENS09(MFR,3,3,3,MX3)
      character*80 filename
      character*1 fchar
      write(fchar,'(i1)')ifr
      if(ifr.eq.1)then
       filename='FILE.TTT'
      else
       filename='FILE.TTT.f'//fchar
      endif
      write(*,*)
      write(*,*)filename
      OPEN(2,FILE=filename)
      WRITE(2,2000)NAT,ifr,fr
2000  FORMAT(' ROA tensors, cartesian derivatives',/,
     1I4,' atoms, freq. ',i2,f11.6/,
     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,(ALPHA09(ifr,I,J,IIND),J=1,3)
2001  FORMAT(I5,1H ,I1,3g15.7)
      WRITE(2,2004)
2004  FORMAT(' The electric dipole magnetic dipole polarizability:',/,
     1       ' Atom/x    jx(Bx)       jy(By)       jz(Bz)')
      DO 2 I=1,3
      WRITE(2,2003)I
2003  FORMAT(' G(',I1,',J):')
      DO 2 L=1,NAT
      DO 2 K=1,3
2     WRITE(2,2001)L,K,(GTENS09(ifr,I,J,3*(L-1)+K),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 M=1,3
3     WRITE(2,2007)L,M,(ATENS09(ifr,K,J,I,3*(L-1)+M)*3.d0/2.d0,K=1,3),L,
     1 M,I,J
2007  FORMAT(I5,1H ,I1,3g15.7,' ',4i3)
      write(2,*)
      write(2,*)'dummy alpha v:'
      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,(ALPHA09(ifr,I,J,IIND),J=1,3)
      CLOSE(2)
      RETURN
      END

c     ===============================      
      SUBROUTINE WRITEPOL(APOL,ifr,ty)
      IMPLICIT INTEGER*4 (I-N) 
      IMPLICIT REAL*8(A-H,O-Z)
      real*8 APOL(*)
      character*80 filename
      character*1 fchar
      character*(*) ty
      write(fchar,'(i1)')ifr
      if(ifr.eq.1)then
       filename=ty
      else
       filename=ty//'.f'//fchar
      endif
      OPEN(90,FILE=filename)
      if(ty.eq.'A.TTT')then
       write(90,*)'A, dipole-quadrupole polarizability'
       write(90,902)
       do 1 IZ=1,3
1      write(90,900)((3.0d0/2.0d0*
     1 APOL(IY+(IX-1)*3+(IZ-1)*9),IY=1,IX),IX=1,3)
      else
       if(ty.eq.'POL.TTT')then
        write(90,*)'Polarizability'
        write(90,902)
        write(90,900)((APOL(IY+(IX-1)*3),IY=1,IX),IX=1,3)
       else
        if(ty.eq.'ORT.TTT')then
         write(90,*)'Optical rotation tensor'
         write(90,902)
         write(90,900)((APOL(IY+(IX-1)*3),IY=1,IX),IX=1,3)
        else 
         if(ty.eq.'GP.TTT')then
          write(90,*)'G tensor'
          write(90,901)
          write(90,900)((APOL(IY+(IX-1)*3),IY=1,3),IX=1,3)
         else
          write(90,*)'Unspecified polarizability'
          write(90,901)
          write(90,900)((APOL(IY+(IX-1)*3),IY=1,3),IX=1,3)
         endif
        endif
       endif
      endif
901   format(
     1'            XX            XY            XZ            YX',
     2'            YY            YZ            ZX            ZY',
     3'            ZZ')
902   format(
     1'            XX            XY            YY            XZ',
     2'            YZ            ZZ')
900   format(9F14.6)
      close(90)
      WRITE(*,*)filename//' written'
      return
      end

c jk-end
C
C\Version=SGI-G94RevB.3\HF=-966.6745649\RMSD=0.000e+00\RMSF=1.196e-01\D
Cipole=0.,0.,0.\PG=C01 [X(C38H52N16O24P4)]\NImag=84\\0.76752455,0.00053
C432,0.03734306,0.23733561,0.00042079,0.13291623,-0.70920516,0.00065060
C,-0.21702343,1.33224819,0.00170226,-0.03103517,0.00051274,0.02526703,1
C.09059550,-0.21867837,0.00010553,-0.09804719,0.00890088,-0.02513206,1.