program trcd
      implicit none
      integer*4 nat,nq,n,i,j,k,iq,jq,n7,ii,jj,ndif1,ndif2
      real*8 t0
      real*8,allocatable::s(:,:),sn(:,:),w(:),wn(:),ss(:,:),
     1c(:,:,:),d(:,:,:),am(:),amn(:)
      integer*4,allocatable::z(:),zn(:),nml(:),nml1(:)
      real*8,allocatable::ALPHAA(:,:),GGA(:,:),AAA(:,:)
      real*8,allocatable::ALPHAAN(:,:),GGAN(:,:),AAAN(:,:)
      real*8,allocatable::alphaq(:,:,:),gq(:,:,:),aq(:,:,:,:)
      real*8,allocatable::alphaqn(:,:,:),gqn(:,:,:),aqn(:,:,:,:)
      real*8,allocatable::p0(:,:),g0(:,:),a0(:,:,:)
      logical lex
      write(6,900)
900   format(' Transforming anharmonic foce constants',/,/,
     1       '  Input:',/,
     1       '             F.INP  old S-matrix',/,
     1       '          FTRY.INP  old FTRY.INP (mass definition)',/,
     1       '          NEWF.INP  new S-matrix',/,
     1       '       NEWFTRY.INP  new FTRY.INP',/,
     1       '       CQQ.SCR.TXT  cubic force constants',/,
     1       '       DQQ.SCR.TXT  quartic force constants',/,
     1       '      TTTQ.TXT.SCR  polarizability derivatives',/,
     1       '     TTTAQ.TXT.SCR  second pol. derivatives',/,/,
     1       '  Output:',/,
     1       '    NEWCQQ.SCR.TXT  new cubic force constants',/,
     1       '    NEWDQQ.SCR.TXT  new quartic force constants',/,
     1       '   NEWTTTQ.TXT.SCR  polarizability derivatives',/,
     1       '  NEWTTTAQ.TXT.SCR  second pol. derivatives',/,/)
      open(9,file='F.INP',status='old')
      read(9,*)nq,n,nat
      close(9)
      if(nq.ne.n)call report('NQ <> 3NAT!')
      allocate(z(nat),w(n),s(n,n),sn(n,n),zn(nat),wn(n),ss(n,n),
     1c(n,n,n),d(n,n,n),am(n),amn(n))
      allocate(ALPHAA(n**2,9),GGA(n**2,9),AAA(n**2,27))
      allocate(ALPHAAN(n**2,9),GGAN(n**2,9),AAAN(n**2,27))
      allocate(alphaq(3,3,n),gq(3,3,n),aq(3,3,3,n))
      allocate(alphaqn(3,3,n),gqn(3,3,n),aqn(3,3,3,n))
      allocate(p0(3,3),g0(3,3),a0(3,3,3),nml(n),nml1(n))
      do 1 i=1,n
      nml(i)=i
1     nml1(i)=i
      call reads(n,s ,w ,   'F.INP',z )
      call reads(n,sn,wn,'NEWF.INP',zn)
      call readY(n,am ,   'FTRY.INP')
      call readY(n,amn,'NEWFTRY.INP')
      call ms(n,s,sn,ss,am,amn)
      inquire(file='NEWDQQ.SCR.TXT',exist=lex)
      if(lex)then
       write(6,*)' NEWDQQ.SCR.TXT exists, calculation skipped'
      else
       call load4(d,n,w)
       call tr4(d,n,ss)
       call save4(d,N,wn)
      endif
      call load3(c,n,w)
      call tc(c,n,ss)
      call save3(c,n,wn)
      call rdttqaq(N,N,ALPHAA,GGA,AAA,nml,ndif2)
      n7=7
      do 2 i=1,3
      do 2 j=1,3
      p0(i,j)=0.0d0
      g0(i,j)=0.0d0
      do 2 k=1,3
2     a0(i,j,k)=0.0d0

      t0=1.0d-6
      do iq=n7,n
       do jq=iq,n
        if(w(iq).gt.t0.and.w(jq).gt.t0)then
         do i=1,9
          ALPHAA(iq+N*(jq-1),i)=ALPHAA(iq+N*(jq-1),i)/dsqrt(w(iq)*w(jq))
          ALPHAA(jq+N*(iq-1),i)=ALPHAA(iq+N*(jq-1),i)
          GGA(iq+N*(jq-1),i)=GGA(iq+N*(jq-1),i)/dsqrt(w(iq)*w(jq))
          GGA(jq+N*(iq-1),i)=GGA(iq+N*(jq-1),i)
         enddo
         do i=1,27
          AAA(iq+N*(jq-1),i)=AAA(iq+N*(jq-1),i)/dsqrt(w(iq)*w(jq))
          AAA(jq+N*(iq-1),i)=AAA(iq+N*(jq-1),i)
         enddo
        endif
       enddo
      enddo

      do iq=n7,n
       do jq=iq,n
        if(w(iq).gt.t0.and.w(jq).gt.t0)then
         do i=1,9
          ALPHAAN(iq+N*(jq-1),i)=0.0d0
          GGAN(iq+N*(jq-1),i)=0.0d0
          do ii=1,n
           do jj=1,n
            ALPHAAN(iq+N*(jq-1),i)=ALPHAAN(iq+N*(jq-1),i)
     1      +ALPHAA(ii+N*(jj-1),i)*ss(iq,ii)*ss(jq,jj)
            GGAN(iq+N*(jq-1),i)=GGAN(iq+N*(jq-1),i)
     1      +GGA(ii+N*(jj-1),i)*ss(iq,ii)*ss(jq,jj)
           enddo
          enddo
         enddo
         do i=1,27
          AAAN(iq+N*(jq-1),i)=0.0d0
          do ii=1,n
           do jj=1,n
            AAAN(iq+N*(jq-1),i)=AAAN(iq+N*(jq-1),i)
     1      +AAA(ii+N*(jj-1),i)*ss(iq,ii)*ss(jq,jj)
           enddo
          enddo
         enddo
        endif
       enddo
      enddo

      do iq=n7,n
       do jq=iq,n
        do i=1,9
         ALPHAAN(iq+N*(jq-1),i)=
     1   ALPHAAN(iq+N*(jq-1),i)*dsqrt(wn(iq)*wn(jq))
         GGAN(iq+N*(jq-1),i)=GGAN(iq+N*(jq-1),i)*dsqrt(wn(iq)*wn(jq))
        enddo
        do i=1,27
         AAAN(iq+N*(jq-1),i)=AAAN(iq+N*(jq-1),i)*dsqrt(wn(iq)*wn(jq))
        enddo
       enddo
      enddo

      call wrttqaq(ndif2,nml,N,N,ALPHAAN,GGAN,AAAN)

      call readtttq(n,alphaq,gq,aq,N,N7,nml1,ndif1)

      do iq=n7,n
       if(w(iq).gt.t0)then
        do i=1,3
         do j=1,3
          alphaq(i,j,iq)=alphaq(i,j,iq)/dsqrt(w(iq))
          gq(i,j,iq)=gq(i,j,iq)/dsqrt(w(iq))
          do k=1,3
           aq(i,j,k,iq)=aq(i,j,k,iq)/dsqrt(w(iq))
          enddo
         enddo
        enddo
       endif
      enddo

      do iq=n7,n
       do i=1,3
        do j=1,3
         alphaqn(i,j,iq)=0.0d0
         gqn(i,j,iq)=0.0d0
         do ii=1,n
          alphaqn(i,j,iq)=alphaqn(i,j,iq)+alphaq(i,j,ii)*ss(iq,ii)
          gqn(i,j,iq)=gqn(i,j,iq)+gq(i,j,ii)*ss(iq,ii)
         enddo
         do k=1,3
          aqn(i,j,k,iq)=0.0d0
          do ii=1,n
           aqn(i,j,k,iq)=aqn(i,j,k,iq)+aq(i,j,k,ii)*ss(iq,ii)
          enddo
         enddo
        enddo
       enddo
      enddo

      do iq=n7,n
       do i=1,3
        do j=1,3
         alphaqn(i,j,iq)=alphaqn(i,j,iq)*dsqrt(w(iq))
         gqn(i,j,iq)=gqn(i,j,iq)*dsqrt(w(iq))
         do k=1,3
          aqn(i,j,k,iq)=aqn(i,j,k,iq)*dsqrt(w(iq))
         enddo
        enddo
       enddo
      enddo

      call wrttq(p0,ndif1,nml1,alphaqn,alphaqn,g0,a0,GQn,AQn,N)
      end

      subroutine save4(cx4,N,e)
      implicit none
      integer*4 N7,N,NQ,nqt,i,j,k
      real*8 cx4(N,N,N),e(*),cm,cqcn,t0
 
c     write sparse quartics
      cm=219470.0d0
      N7=7
      open(23,file='NEWDQQ.SCR.TXT')
      nQ=0
      nqt=0
      t0=1.0d-6
      do 621 i=n7,n
      do 621 j=n7,n
      do 621 k=j,n
      if(dabs(e(i)).gt.t0.and.dabs(e(j)).gt.t0.and.dabs(e(k)).gt.t0)then
       cqcn=cx4(i,j,k)*dsqrt(abs(e(i)**2*e(j)*e(k)))
       if(dabs(cqcn*cm).lt.0.1d0)goto 621
       write(23,5333)I,J,K,cqcn,cqcn*cm
5333   format(3I6,G16.8,' ',f9.1)
       nq=nq+1
      endif
621   nqt=nqt+1
      close(23)
      write(6,6001)nqt,nq,0.1d0
6001  format('number of normal quartics ',I10,/,
     1       'Selected                  ',I10,/,
     1       'Limit                     ',F10.7,/,
     1       'Writen in NEWDQQ.SCR.TXT')
      return
      end

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

      SUBROUTINE READS(N3,S,E,fn,z)
      IMPLICIT none
      integer*4 N3,z(*),i,ix,j
      real*8 S(N3,N3),E(*)
      character*(*)fn
      OPEN(4,FILE=fn,STATUS='OLD')
      read(4,*)
      do 1 i=1,N3/3
1     read(4,*)z(i)
      read(4,*)
      DO 2 I=1,N3/3
      DO 2 J=1,N3
2     read(4,*)(s(3*(i-1)+ix,N3-j+1),ix=1,2),
     1(s(3*(i-1)+ix,N3-j+1),ix=1,3)
      read(4,*)
      READ(4,4000)(E(N3-I+1),I=1,N3)
4000  FORMAT(6F11.6)
      close(4)
      write(6,*)fn//' read'
      RETURN
      end

      SUBROUTINE ms(n,s,sn,ss,am,amn)
      IMPLICIT none
      integer*4 n,i,j,k
      real*8 s(n,n),sn(n,n),ss(n,n),am(*),amn(*)
      do 1 i=1,n
      do 1 j=1,n
      ss(i,j)=0.0d0
      do 1 k=1,n
1     ss(i,j)=ss(i,j)+sn(k,i)*s(k,j)*dsqrt(am(k)*amn(k))
      return
      end

      subroutine load3(C,N,e)
      implicit none
      real*8 C(N,N,N),AC,e(*)
      integer*4 I,J,K,N,nqt
      call tz(c,N)
      nqt=0
      OPEN(4,file='CQQ.SCR.TXT')
8     READ(4,*,end=88,err=88)I,J,K,AC
      if(e(I).ne.0.0d0.and.e(J).ne.0.0d0.and.e(K).ne.0.0d0)
     1AC=AC/dsqrt(e(I)*e(J)*e(K))
      C(I,J,K)=AC
      C(I,K,J)=AC
      C(J,K,I)=AC
      C(J,I,K)=AC
      C(K,J,I)=AC
      C(K,I,J)=AC
      nqt=nqt+1
      goto 8
88    CLOSE(4)
      write(6,881)nqt
881   format(I10,' constants read from CQQ.SCR.TXT')
      return
      end

      subroutine load4(d,n,e)
      implicit none
      integer*4 n,I,J,K,IT
      real*8 d(n,n,n),AC,e(*)
      call tz(d,N)
      IT=0
      open(23,file='DQQ.SCR.TXT')
7     read(23,*,err=7000,end=7000)I,J,K,AC
      IT=IT+1
      if(e(I).ne.0.0d0.and.e(J).ne.0.0d0.and.e(K).ne.0.0d0)
     1AC=AC/(e(I)*dsqrt(e(J)*e(k)))
      d(I,J,K)=AC
      d(I,K,J)=AC
      goto 7
7000  close(23)
      WRITE(6,7004)IT
7004  FORMAT(I10,' quartics from DQQ.SCR.TXT')
      return
      end

      subroutine tz(d,N)
      implicit none
      integer*4 I,J,K,N
      real*8 d(N,N,N)
      do 8 I=1,N
      do 8 J=1,N
      do 8 K=1,N
8     d(I,J,K)=0.0d0
      return
      end

      subroutine save3(c,n,e)
      implicit none
      integer*4 i,j,k,n,nq,nqt,N7
      real*8 c(n,n,n),e(*),cqau,cqcn,cqcm,cm,tol3,t0
      tol3=0.1d0
      t0=1.0d-6
      cm=219470.0d0
      open(23,file='NEWCQQ.SCR.TXT')
      nq=0
      nqt=0
      N7=7
      do 62 i=N7,N
      do 62 j=i,N
      do 62 k=j,N
      if(dabs(e(i)).gt.t0.and.dabs(e(j)).gt.t0.and.dabs(e(k)).gt.t0)then
       cqau=C(I,J,K)
       cqcn=cqau
       cqcn=cqau*dsqrt(e(i)*e(j)*e(k))
       cqcm=cqcn*cm
       if(dabs(cqcm).gt.tol3)then
        write(23,2423)I,J,K,cqcn,cqcm
2423    format(3I5,G16.8,' ',F9.1)
        nq=nq+1
       endif
      endif
62    nqt=nqt+1
      close(23)
      write(6,6000)nqt,nq,tol3
6000  format('Number of normal cubics ',I10,/,
     1       'Selected                ',I10,/,
     1       'Limit                   ',F10.7,/,
     1       'Writen in NEWCQQ.SCR.TXT',/)
      return
      end

      subroutine tr4(CX4,N,S)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      real*8 CX4(N,N,N),S(N,N),TWO
      real*8,allocatable:: CQ4(:,:,:)
      allocate(CQ4(N,N,N))
      Z0=0.0d0
      TWO=2.0d0
      WRITE(6,*)' Transforming quartics, semidiagonals ...'
C     Calculate D(I,I,J,K)=S(A,I)*S(B,I)*S(C,J)*S(D,K)*DX(A,B,C,D):
      DO 46 I=1,N
      WRITE(6,*)I,' of ',N
      DO 46 J=1,N
      DO 46 K=J,N
      S1=Z0
      S2=Z0
      S3=Z0
      S4=Z0
      DO 2 IC=1,N
      SIA=S(I,IC)
      SJA=S(J,IC)
      SKA=S(K,IC)
C     iiii:
      S1=S1+CX4(IC,IC,IC)*SIA*SIA*SJA*SKA
      DO 3 JC=1,N
      IF(JC.EQ.IC)GOTO 3
      SIB=S(I,JC)
      SJB=S(J,JC)
      SKB=S(K,JC)
      IF(JC.GT.IC)THEN
C      iijj=jiij=ijij (and i>j):
       S2=S2+CX4(IC,JC,JC)*
     1  (SIA*(SIA*SJB*SKB+SIB*(SJA*SKB+SJB*SKA))+
     2   SIB*(SIB*SJA*SKA+SIA*(SJB*SKA+SJA*SKB)) )
C      iiij=iiji=ijii=jiii (and i>j):
       S3=S3+CX4(IC,IC,JC)*(SIA*(SIA*(SJA*SKB+SJB*SKA) +
     1                           SIB*(SJA*SKA+SJA*SKA)))+
     2       CX4(JC,JC,IC)*(SIB*(SIB*(SJB*SKA+SJA*SKB) +
     3                           SIA*(SJB*SKB+SJB*SKB)))
      ENDIF
      S01=(SIA*SIA*SJB+TWO*SIA*SJA*SIB)
      S02=(SIA*SIA*SKB+TWO*SIA*SKA*SIB)
      S03=TWO*(SIA*SJA*SKB+SIA*SKA*SJB+SIB*SJA*SKA)
C      iijk=ijik=jiik=jiki=ijki=jkii (and k<j):
      DO 4 KC=JC+1,N
      IF(KC.EQ.IC)GOTO 4
      S4=S4+CX4(IC,JC,KC)*(S(K,KC)*S01+S(J,KC)*S02+S(I,KC)*S03)
4     CONTINUE
3     CONTINUE
2     CONTINUE
46    CQ4(I,J,K)=S1+S2+S3+S4
C
      DO 54 I=1,N
      DO 54 J=1,N
      DO 54 K=J+1,N
54    CQ4(I,K,J)=CQ4(I,J,K)

      DO 55 I=1,N
      DO 55 J=1,N
      DO 55 K=1,N
55    CX4(I,K,J)=CQ4(I,J,K)

      return
      end

      subroutine tc(c,n,ss)
      implicit none
      integer*4 n,i,j,k,ii,jj,kk
      real*8 c(n,n,n),ss(n,n)
      real*8,allocatable::t(:,:,:)
      allocate(t(n,n,n))
      write(6,*)'first index-c'
      do 1 I=1,n
      do 1 j=1,n
      do 1 k=1,n
      t(I,j,k)=0.0d0
      do 1 ii=1,n
1     t(I,j,k)=t(I,j,k)+ss(I,ii)*c(ii,j,k)
      write(6,*)'second index-c'
      do 2 I=1,n
      do 2 J=1,n
      do 2 k=1,n
      c(I,J,k)=0.0d0
      do 2 jj=1,n
2     c(I,J,k)=c(I,J,k)+ss(J,jj)*t(I,jj,k)
      write(6,*)'third index-c'
      do 3 I=1,n
      do 3 J=1,n
      do 3 K=1,n
      t(I,J,K)=0.0d0
      do 3 kk=1,n
3     t(I,J,K)=t(I,J,K)+ss(K,kk)*c(I,J,kk)
      do 4 I=1,n
      do 4 J=1,n
      do 4 K=1,n
4     c(i,j,k)=t(I,J,K)
      return
      end

      subroutine readY(n,m,s)
      implicit none
      integer*4 n,i
      real*8 m(*)
      character*(*) s
      open(11,file=s,status='old')
      read(11,*)
      read(11,*)
      read(11,1000)(m(i),i=1,n,3)
1000  format(6f12.6)
      read(11,*)
      read(11,1000)(m(i),i=1,n,3)
      close(11)
      do 1 i=1,n,3
      m(i+1)=m(i)
1     m(i+2)=m(i)
      write(6,*)s//' read'
      return
      end

c     ============================================================
      subroutine readtttq(NT3,alphaq,gq,aq,N,N7,nml1,ndif1)
      implicit none
      integer*4 NT3,N,iq,i,j,k,l,N7,ndif1,ix,jx,kx,nml1(*)
      real*8 alphaq(3,3,NT3),gq(3,3,NT3),aq(3,3,3,NT3)

      do 4 iq=N7,N
      do 4 i=1,3
      do 4 j=1,3
      alphaq(i,j,iq)=0.0d0
      gq(i,j,iq)=0.0d0
      do 4 k=1,3
4     aq(i,j,k,iq)=0.0d0

      open(9,file='TTTQ.TXT.SCR')
      do 1 i=1,4
1     read(9,*)
      read(9,*)ndif1
      DO 19 i=1,ndif1
      do 19 j=1,3
      do 19 k=1,3
      read(9,*)iq,iq,ix,jx,alphaq(ix,jx,iq),alphaq(ix,jx,iq)
19    nml1(i)=iq
      do 2 i=1,4
2     read(9,*)
      DO 20 i=1,ndif1
      do 20 j=1,3
      do 20 k=1,3
20    read(9,*)iq,iq,ix,jx,gq(ix,jx,iq)
      do 3 i=1,10
3     read(9,*)
      DO 21 i=1,ndif1
      do 21 j=1,3
      do 21 k=1,3
      do 21 l=1,3
21    read(9,*)iq,iq,ix,jx,kx,aq(ix,jx,kx,iq)
      close(9)
      call w36('Polarizability derivatives from TTTQ.TXT.SCR')
      return
      end

      subroutine rdttqaq(N,NT3,ALPHAA,GGA,AAA,nml,ndif2)
      implicit none
      integer*4 ndif2,i,j,N,NT3,ix,ii,jj,jx,nml(*)
      real*8 ALPHAA(NT3**2,9),GGA(NT3**2,9),AAA(NT3**2,27)
      logical lex

      do 1 i=1,N**2
      do 1 ix=1,9
      ALPHAA(i,ix)=0.0d0
      GGA(i,ix)=0.0d0
      AAA(i,ix)=0.0d0
      AAA(i,ix+9)=0.0d0
1     AAA(i,ix+18)=0.0d0

      inquire(file='TTTAQ.TXT.SCR',exist=lex)
      if(lex)then
       open(21,file='TTTAQ.TXT.SCR')
       read(21,*)
       read(21,*)ndif2
       do 23 i=1,ndif2
       do 23 j=i,ndif2
       read(21,*)ii,jj,ii,jj
       nml(j)=jj
       read(21,*)(ALPHAA(ii+N*(jj-1),ix),ix=1,3)
       read(21,*)(ALPHAA(ii+N*(jj-1),ix),ix=4,6)
23     read(21,*)(ALPHAA(ii+N*(jj-1),ix),ix=7,9)
       do 24 i=1,ndif2
       do 24 j=i,ndif2
       read(21,*)ii,jj,ii,jj
       read(21,*)(GGA(ii+N*(jj-1),ix),ix=1,3)
       read(21,*)(GGA(ii+N*(jj-1),ix),ix=4,6)
24     read(21,*)(GGA(ii+N*(jj-1),ix),ix=7,9)
       do 25 i=1,ndif2
       do 25 j=i,ndif2
       read(21,*)ii,jj,ii,jj
       do 25 jx=0,8
25     read(21,*)(AAA(ii+N*(jj-1),ix+jx*3),ix=1,3)
       close(21)
       call w36('TTTAQ.TXT.SCR read')
      else
       call w36('TTTAQ.TXT.SCR not found')
      endif
      end

      subroutine w36(s)
      character*(*) s
      write(6,*)s
      return
      end

      subroutine wrttqaq(ndif2,nml,N,NT3,ALPHAA,GGA,AAA)
      implicit none
      integer*4 ndif2,i,j,nml(*),N,NT3,ix
      real*8 ALPHAA(NT3**2,9),GGA(NT3**2,9),AAA(NT3**2,27)
      open(21,file='NEWTTTAQ.TXT.SCR')
      write(21,*)' Normal mode polarizability second derivatives'
      write(21,*)ndif2
      do 23 i=1,ndif2
      do 23 j=i,ndif2
      write(21,2124)i,j,nml(i),nml(j)
2124  format(4i5,' modes')
23    write(21,2121)(ALPHAA(nml(I)+N*(nml(J)-1),ix),ix=1,9)
2121  format(3f15.6)
      do 24 i=1,ndif2
      do 24 j=i,ndif2
      write(21,2124)i,j,nml(i),nml(j)
24    write(21,2121)(GGA(nml(I)+N*(nml(J)-1),ix),ix=1,9)
      do 25 i=1,ndif2
      do 25 j=i,ndif2
      write(21,2124)i,j,nml(i),nml(j)
25    write(21,2121)(AAA(nml(I)+N*(nml(J)-1),ix),ix=1,27)
      close(21)
      call w36('NEWTTTAQ.TXT.SCR written')
      end

      subroutine wrttq(polar0,ndif1,nml1,ALPHA0Q,alphan,gten0,aten0,
     1G0Q,A0Q,NT3)
      implicit none
      integer*4 NT3,ndif1
      real*8 polar0(3,3),gten0(3,3),aten0(3,3,3),A0Q(3,3,3,NT3),
     1ALPHA0Q(3,3,NT3),G0Q(3,3,NT3),alphan(3,3,NT3)
      integer*4 i,j,k,nml1(*),l,iq
      open(9,file='NEWTTTQ.TXT.SCR')
      write(9,9004)((polar0(i,j),i=1,3),j=1,3)
9004  format(3(3f15.7,' au',/),
     1' Normal mode polarizability derivatives num anal')
      write(9,*)ndif1
      DO 19 i=1,ndif1
      iq=nml1(i)
      do 19 j=1,3
      do 19 k=1,3
19    write(9,9005)i,iq,j,k,alphan(j,k,iq),ALPHA0Q(j,k,iq)
9005  format(4i5,2f12.6)
      write(9,9006)((gten0(i,j),i=1,3),j=1,3)
9006  format(3(3f15.7,' G au',/),
     1' Normal mode G-tensor derivatives:')
      DO 20 i=1,ndif1
      iq=nml1(i)
      do 20 j=1,3
      do 20 k=1,3
20    write(9,9005)i,iq,j,k,G0Q(j,k,iq)
      write(9,9008)(((aten0(i,j,k),i=1,3),j=1,3),k=1,3)
9008  format(9(3f15.7,' A au',/),
     1' Normal mode A-tensor derivatives:')
      DO 21 i=1,ndif1
      iq=nml1(i)
      do 21 j=1,3
      do 21 k=1,3
      do 21 l=1,3
21    write(9,9007)i,iq,j,k,l,A0Q(j,k,l,iq)
9007  format(5i5,2f13.6)
      close(9)
      call w36('nm pol der in NEWTTTQ.TXT.SCR')
      return
      end