PROGRAM DODCGN
      IMPLICIT none
      integer*4 nt0,it,i,N,J,nc,ns,ix,jx,kx,k,ntem,ig,is,NG,nat,
     1ierr,ITER,NATB,ii,nm,ivel,imag,ntr,nd,imax,imin,nmol
      PARAMETER (nt0=3550)
c     Limits:
c     nt0 -# of transitions on one chromophore
c     npar0 -# of parameters to follow
c     n40 -# of potential sites
c     nwl0 -# of fitted transitions
      REAL*8 rc(3),mag(nt0,3),elel(nt0,3),elev(nt0,3),en(nt0),
     1en1(nt0,3),A,RTOL,um(3),ue(3),eau,ecm(nt0),eev(nt0),enm(nt0),
     1elelc(nt0,3),rcs(3),alr(3,3),gr(3,3),ar(3,3,3),anorm
      REAL*8,allocatable:: R(:),rn(:)
      integer*4,allocatable:: id(:)
      character*80 s80,output
      character*80,allocatable::listst(:)
      character*1 prog
      logical lvel,llen,lene,lmag,lgeo,lttt,ldistr
      COMMON/MATVAR/A(3,3),RTOL,NATB,ITER
      integer*4,allocatable::ntw(:),ics(:),ic(:)
      real*8,allocatable::alphad(:,:,:,:),Gpd(:,:,:,:),
     1Ad(:,:,:,:,:),tw(:,:),XS(:,:),XB(:,:),wc(:)
c
      write(6,6000)
6000  format(
     1' Program takes dipoles from gaussian/turbomole outputs,',/,
     1'          (F.INPs for vibrations),',/,
     1'      and spreads them on large molecule',/,/,
     1'     Input: DODCN.CON    options',/,
     1'            outputs listed in DODC.CON',/,
     1'            BIG.X  target molecule',/,
     1'            FILE.Q.TTT transition polarizabilities',/,
     1'            group list files listed in DODC.CON',/,/,
     1'     Output:DC.SC',/)
c
c     Read in transition energies and dipoles from Gaussian or Turbomole:
      lvel=.false.
      llen=.false.
      lene=.false.
      lmag=.false.
      lgeo=.false.
      lttt=.false.

c     read-in geometry
      OPEN(2,FILE='BIG.X',STATUS='OLD')
      read(2,*)
      read(2,*)N
      allocate (R(3*N))
      allocate (rn(3*N))
      do 1 i=1,N
1     READ(2,*)R(1+3*(i-1)),(R(J+3*(i-1)),J=1,3)
      close(2)
      write(6,*)N, ' atoms in BIG.X'
c
      open(4,file='DODCN.CON')
c     gaussian or turbomole:
      read(4,401)prog
401   format(a1)
c     velocity (0) magnetic moments (1):
      read(4,*)ivel,imag
      if(imag.eq.1)then
       write(6,*)' magnetic moments added'
      else
       write(6,*)' magnetic moments skipped'
      endif
c     number of chromophore types:
      read(4,*)nmol
      allocate(listst(nmol))

      NG=0
      do 1101 ii=1,nmol
      read(4,'(a)')listst(ii)
      open(11,file=listst(ii))
      read(11,*)i
      NG=NG+i
1101  close(11)

c     write DC.SC, reread DODC.CON:
      open(3,file='DC.SC')
      write(3,*)NG,' total groups'
      write(6,*)NG,' total groups'
      write(6,*)nmol,' source molecules'

      IG=0
      do 1000 ii=1,nmol
      read(4,'(a)')output
      write(6,'(a)')output,':'
      read(4,*)imin,imax

      if(prog.eq.'g'.or.prog.eq.'G')then
c      gaussian electric:
       open(30,file=output,status='old')
10     read(30,300,end=949,err=949)s80
300    format(a80)
       if(s80(12:44).eq.'excited state Transition electric')then
        call rd(nt0,it,elel,llen)
        write(6,*)it,' electric dipoles'
        goto 10
       endif
       if(s80(12:44).eq.'excited state transition electric')then
        call rd(nt0,it,elel,llen)
        write(6,*)it,' electric dipoles'
        goto 10
       endif
       if(s80(12:46).eq.'excited state transition velocity d')then
        call rd(nt0,it,elev,lvel)
        write(6,*)it,' electric velocity dipoles'
        goto 10
       endif
       if(s80(12:44).eq.'excited state transition magnetic')then
        call rd(nt0,it,mag,lmag)
        write(6,*)it,' magnetic dipoles'
        goto 10
       endif
       if(s80(2:15).eq.'Excited State '.and.s80(16:16).ne.'s')then
        lene=.true.
        do 203 i=1,79
        if(s80(i:i).eq.':')read(s80(15:i-1),*)nd
203     if(s80(i:i+1).eq.'nm')read(s80(i- 9:i-2),*)en(nd)
       endif
       if(s80(27:44).eq.'Input orientation:')call roi(rn,nat,lgeo)
       if(s80(26:46).eq.'Standard orientation:')call ro(rn,nat,lgeo)
       goto 10
949    close(30)
       write(6,*)nd,' energy transitions in ',output
       do 481 i=1,it
c      nm to hartrees:
       eau=1.0d7/(8065.54476345045d0*27.211384205943d0)/en(i)
       do 481 ix=1,3
c      correct also gaussian sign:
481    elev(i,ix)=-elev(i,ix)/eau
      else
       if(prog.eq.'t'.or.prog.eq.'T')then
c       turbomole electric:
        open(30,file=output,status='old')
100     read(30,300,end=9490,err=9490)s80
        if(s80(2:46).eq.'Transition dipole moments in length formalism')
     1  call rdt(nt0,it,elel,llen)
        if(s80(2:46).eq.'Transition dipole moments in velocity formali')
     1  call rdt(nt0,it,elev,lvel)
        if(s80(2:28).eq.'Transition magnetic moments')
     1  call rdt(nt0,it,mag,lmag)
        if(s80(2:24).eq.'Energies of Transitions')
     1  call rdt(nt0,it,en1,lene)
        if(s80(2:19).eq.'Atomic coordinates')call rot(rn,nat,lgeo)
        goto 100
9490    close(30)
        do 9432 i=1,it
9432    en(i)=en1(i,2)
        write(6,*)it,' transitions in moments'
       else
        if(prog.eq.'f'.or.prog.eq.'F')then
c        vibrational from F.INP:
         open(30,file=output,status='old')
         read(30,*)it,nat,nat
         allocate(id(nat))
         close(30)
         call rvd(output,nt0,nat,it,rn,id,mag,elel,en)
         inquire(file='FILE.Q.TTT',exist=lttt)
         if(lttt)then
          open(15,file='FILE.Q.TTT')
          READ(15,*)
          READ(15,*)nm
          close(15)
          allocate(alphad(1,nm,3,3),Gpd(1,nm,3,3),Ad(1,nm,3,3,3),
     1    tw(1,nm+1),ntw(1+1))
          ntw(1)=nm
c         read derivatives/transition polarizabilities/, molecular frame
          CALL READTTTQ(nm,1,alphad,Gpd,Ad,'FILE.Q.TTT',tw,1,ntw,nm)
          write(6,*)'Normal mode polarizability derivatives read'
         endif
        else
         call report(prog//' unknown option')
        endif
       endif
      endif

      if(prog.ne.'f'.and.prog.ne.'F')then
       if(.not.llen)call report(' Electric length dipoles not found')
       if(.not.lvel)call report(' Electric velocity dipoles not found')
       if(.not.lmag)call report(' Magnetic dipoles not found')
       if(.not.lene)call report(' Energies not found')
       if(.not.lgeo)call report(' Geometry not found')
       if(ivel.eq.1)then
        write(6,*)' velocity dipoles used'
        do 54 i=1,it
        do 54 ix=1,3
54      elel(i,ix)=elev(i,ix)
       else
        write(6,*)' length dipoles used'
       endif
c      transform aus to Bohr magnetons, debyes
       do 56 i=1,it
       do 56 ix=1,3
c      this is not sure, check the magnetic moment units!:
       mag(i,ix)=mag(i,ix)/2.0d0*2.541765d0
56     elel(i,ix)=elel(i,ix)*2.541765d0
       if(imag.eq.0)then
        do 55 i=1,it
        do 55 ix=1,3
55      mag(i,ix)=0.0d0
       endif
      else
      endif

      if(lttt)then
       open(33,file='DC.TT')
       write(33,*)NG,' total groups'
      endif
c
c     number of considered chromophore atoms, e.g.:   H N C O, and list:
      read(4,*)nc
c     distribute the dipoleis over all atoms:
      ldistr=nc.lt.0
      nc=abs(nc)
      backspace 4
      allocate(ics(nc),ic(nc),XS(nc,3),XB(nc,3),wc(nc))
      read(4,*)ics(1),(ics(i),i=1,nc)
      write(6,6009)(ics(i),i=1,nc)
6009  format('    atoms ',8i4)
c     atomic weights for placing of the dipole:
      read(4,*)(wc(i),i=1,nc)
      anorm=0.0d0
      do 59 i=1,nc
59    anorm=anorm+wc(i)
      do 60 i=1,nc
60    wc(i)=wc(i)/anorm
c     make center of the "small" part:
      call makecenter(rcs,rn,wc,nc,ics)
      open(8,file=listst(ii))
      read(8,*)ntem
      do 43 j=1,ntem
      read(8,*)ic(1),(ic(i),i=1,nc)
      write(6,6019)(ic(i),i=1,nc)
6019  format(' on atoms ',8i4)
c     make center of the "big" part:
      call makecenter(rc,R,wc,nc,ic)
c     write group reference atoms into XS and corresponding ones in XB:
      do 49 i=1,nc
      do 49 ix=1,3
      XB(i,ix)= R(ix+3*(ic( i)-1))-rc( ix)
49    XS(i,ix)=rn(ix+3*(ics(i)-1))-rcs(ix)
      NATB=nc
c     overlap XS,XB, make the rotation matrix, store in A:
      call domatrix(ierr,nc,XB,XS)
      if(ierr.ne.0)call report('rotation not successful')
      write(6,*)' Transformation matrix'
      do 67671 i=1,3
67671 write(6,6767)(A(i,k),k=1,3)
6767  format(3f10.5)

c     write transition energies and dipoles in temp variables:
      ntr=0
      do 472 i=1,it
      if(i.ge.imin.and.i.le.imax)then
       ntr=ntr+1
       if(prog.ne.'f'.and.prog.ne.'F')then
        enm(ntr)=en(i)
        ecm(ntr)=1.0d7/enm(ntr)
       else
        ecm(ntr)=en(i)
        enm(ntr)=1.0d7/ecm(ntr)
       endif
       eev(ntr)=ecm(ntr)/8065.54093d0
       do 473 ix=1,3
c      borrow elev variable for magnetic:
       if(imag.eq.0)mag(i,ix)=0.0d0
       elev(ntr,ix)=mag(i,ix)
473    elelc(ntr,ix)=elel(i,ix)
      endif
472   continue
      ns=ntr+1

      IG=IG+1
      if(ldistr)then
c      indicate distribution of the dipoles over the weighted atoms
       write(3,*)-IG,' . group:'
       write(3,1040)(rc(i)/10.0d0,i=1,3)
1040   format(3f15.6,' nm')
       write(3,*)nc,' sites (atoms):'
       do 58 i=1,nc
58     write(3,1041)(R(ix+3*(ic(i)-1))/10.0d0,ix=1,3),wc(i)
1041   format(4f15.6)
      else
       write(3,*)IG,' . group:'
       write(3,1040)(rc(i)/10.0d0,i=1,3)
      endif
      if(lttt)write(33,*)IG,' . group:'
      write(3,*)ns,' states'
      write(3,*)
c     ground state:
      write(3,*)'0 0.0'
      do 48 is=2,ns
48    write(3,1001)is,eev(is-1),ecm(is-1),enm(is-1)
1001  format(i5,f13.6,' eV =',f11.2,' cm-1 =',f11.2,' nm')
      write(3,*)
c     write rotated dipole moments:
      do 52 is=2,ns
      call rotv(A,ue,is-1,nt0,elelc)
      call rotv(A,um,is-1,nt0,elev )
52    write(3,1002)is,(ue(ix),ix=1,3),(um(ix),ix=1,3)
1002  format(i5,6f12.6)

      if(lttt)then
c      rotate and write transition polarizabilities:
       ntr=0
       do 57 i=1,it
       if(i.ge.imin.and.i.le.imax)then
        ntr=ntr+1
        call rot2(A,alr,i,1,1,nm,alphad)
        call rot2(A,Gr,i,1,1,nm,Gpd)
        call rot3(A,Ar,i,1,1,nm,Ad)
        write(33,1003)i,ecm(ntr),tw(1,i),
     1  ( (alr(ix,jx   ),ix=1,3),jx=1,3),
     1  ( (Gr( ix,jx   ),ix=1,3),jx=1,3),
     1  (((Ar( ix,jx,kx),ix=1,3),jx=1,3),kx=1,3)
1003    format(i5,2f12.2,/,/,
     1  3(3f12.6,/),/,3(3f12.6,/),/,9(3f12.6,/),/)
       endif
57     continue
      endif

43    continue
      close(8)
1000  continue
c
      close(3)
      if(lttt)close(33)
      close(4)
      call report('DC.SC written')
      end
c
      SUBROUTINE REPORT(RS)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER*(*) RS
      WRITE(6,3000)
3000  FORMAT(/,80(1H*))
      WRITE(6,*)RS
      WRITE(6,3001)
3001  FORMAT(80(1H*),/,/,'PROGRAM STOPPED')
      STOP
      END

      subroutine rd(nt0,it,e,l)
      implicit none
      integer*4 nt0,it,i
      character*80 s80
      real*8 e(nt0,3)
      logical nu,l
      l=.true.
      it=0
      read(30,*)
31    read(30,300)s80
300   format(a80)
      if(nu(s80(10:10)).or.s80(8:10).eq.'***')then
       it=it+1
       if(it.gt.nt0)call report(' too many transitions !')
       read(s80(11:len(s80)),*)(e(it,i),i=1,3)
       goto 31
      else
       backspace 30
       return
      endif
      end

      subroutine rdt(nt0,it,e,l)
      implicit none
      integer*4 nt0,it,i
      character*10 s10
      real*8 e(nt0,3)
      logical nu,l
      l=.true.
      it=0
      read(30,*)
31    read(30,300)s10
300   format(a10)
      if(nu(s10(4:4)))then
       it=it+1
       if(it.gt.nt0)then
        write(6,*)' too many transitions !'
        stop
       endif
       backspace 30
       read(30,*)(e(it,i),i=1,1),(e(it,i),i=1,3)
       goto 31
      endif
      backspace 30
      return
      end


      subroutine ro(rn,nat,l)
      implicit none
      integer*4 i,nat
      character*80 s80
      real*8 rn(*)
      logical nu,l

      l=.true.
      nat=0
      read(30,*)
      read(30,*)
      read(30,*)
      read(30,*)
31    read(30,300)s80
300   format(a80)
      if(nu(s80(7:7)))then
       nat=nat+1
       read(s80,*)(rn(i+3*(nat-1)),i=1,3),(rn(i+3*(nat-1)),i=1,3)
       goto 31
      endif
      write(6,*)nat,' atoms in chromophore'
      return
      end

      subroutine roi(rn,nat,l)
      implicit none
      integer*4 i,nat
      character*80 s80
      real*8 rn(*)
      logical nu,l

      l=.true.
      nat=0
      read(30,*)
      read(30,*)
      read(30,*)
      read(30,*)
31    read(30,300)s80
300   format(a80)
      if(nu(s80(7:7)))then
       nat=nat+1
       read(s80,*)(rn(i+3*(nat-1)),i=1,3),(rn(i+3*(nat-1)),i=1,3)
       goto 31
      endif
      write(6,*)nat,' atoms in chromophore'
      return
      end

      subroutine rot(rn,nat,l)
      implicit none
      integer*4 i,nat
      character*80 s80
      character*15 dum
      real*8 rn(*)
      logical l
      l=.true.
      nat=0
      read(30,*)
31    read(30,300)s80
300   format(a80)
      if(s80(1:15).ne.'               ')then
       nat=nat+1
       read(s80,*)dum,(rn(i+3*(nat-1)),i=1,3)
       goto 31
      endif
      write(6,*)nat,' atoms in chromophore'
      return
      end

      function nu(c)
      character*1 c
      logical nu
      nu=c.eq.'1'.or.c.eq.'2'.or.c.eq.'3'.or.c.eq.'4'.or.c.eq.'5'.
     1or.c.eq.'6'.or.c.eq.'7'.or.c.eq.'8'.or.c.eq.'9'.or.c.eq.'0'
      return
      end

      SUBROUTINE domatrix(ierr,nc,XB,XS)
C     This is a Fortran version of Petr Malon's subroutine Amoeba
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3)
      real*8 XS(nc,3),XB(nc,3)
      COMMON/MATVAR/A(3,3),RTOL,NAT,ITER
      IERR=0
C     STARTING VALUES FOR THE ITERATION:
      ANG(1)=0.2d0
      ANG(2)=0.0d0
      ANG(3)=0.0d0
      DO 1 I=1,3
1     P(1,I)=ANG(I)
      Y(1)=FU(ANG,nc,XS,XB)
      ANG(2)=ANG(1)
      ANG(1)=0.0d0
      DO 2 I=1,3
2     P(2,I)=ANG(I)
      Y(2)=FU(ANG,nc,XS,XB)
      ANG(3)=ANG(2)
      ANG(2)=0.0d0
      DO 3 I=1,3
3     P(3,I)=ANG(I)
      Y(3)=FU(ANG,nc,XS,XB)
      ANG(1)=ANG(3)
      ANG(2)=ANG(3)
      ANG(3)=0.0d0
      DO 4 I=1,3
4     P(4,I)=ANG(I)
      Y(4)=FU(ANG,nc,XS,XB)
      FTOL=0.0000001d0
      ITMAX=1500
      ITER=0
99999 ILO=1
      IF(Y(1).GT.Y(2))THEN
       IHI=1
       INHI=2
      ELSE
       IHI=2
       INHI=1
      ENDIF
      DO 5 I=1,4
      IF(Y(I).LT.Y(ILO))ILO=I
      IF(Y(I).GT.Y(IHI))THEN
       INHI=IHI
       IHI=I
      ELSE
       IF(Y(I).GT.Y(INHI))THEN
        IF(I.NE.IHI)INHI=I
       ENDIF
      ENDIF
5     CONTINUE
      IF(ABS(Y(IHI))+ABS(Y(ILO)).LT.0.00000000001d0)THEN
       RTOL=0.0d0
      ELSE
       RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0
      ENDIF
      IF(RTOL.LT.FTOL)RETURN
      IF(ITER.EQ.ITMAX)THEN
       WRITE(3,*)' Rotation has not converged !'
       IERR=1
       RETURN
      ENDIF
      ITER=ITER+1
      DO 55 I=1,3
55    PBAR(I)=0.0d0
      DO 6 I=1,4
      IF(I.NE.IHI)THEN
       DO 7 J=1,3
7      PBAR(J)=PBAR(J)+P(I,J)
      ENDIF
6     CONTINUE
      DO 8 J=1,3
      PBAR(J)=PBAR(J)/3.0d0
8     PR(J)=2.0d0*PBAR(J)-P(IHI,J)
      YPR=FU(PR,nc,XS,XB)
      IF(YPR.LE.Y(ILO))THEN
       DO 9 J=1,3
9      PRR(J)=2.0d0*PR(J)-PBAR(J)
       YPRR=FU(PRR,nc,XS,XB)
       IF(YPRR.LT.Y(ILO))THEN
        DO 10 J=1,3
10      P(IHI,J)=PRR(J)
        Y(IHI)=YPRR
       ELSE
        DO 11 J=1,3
11      P(IHI,J)=PR(J)
        Y(IHI)=YPR
       ENDIF
      ELSE
       IF(YPR.GE.Y(INHI))THEN
        IF(YPR.LT.Y(IHI))THEN
         DO 12 J=1,3
12       P(IHI,J)=PR(J)
         Y(IHI)=YPR
        ENDIF
        DO 13 J=1,3
13      PRR(J)=0.5d0*P(IHI,J)+0.5d0*PBAR(J)
        YPRR=FU(PRR,nc,XS,XB)
        IF(YPRR.LT.Y(IHI))THEN
         DO 14 J=1,3
14       P(IHI,J)=PRR(J)
         Y(IHI)=YPRR
        ELSE
         DO 15 I=1,4
         IF(I.NE.ILO)THEN
          DO 16 J=1,3
          PR(J)=0.5d0*(P(I,J)+P(ILO,J))
16        P(I,J)=PR(J)
          Y(I)=FU(PR,nc,XS,XB)
         ENDIF
15       CONTINUE
        ENDIF
       ELSE
        DO 17 J=1,3
17      P(IHI,J)=PR(J)
        Y(IHI)=YPR
       ENDIF
      ENDIF
      GOTO 99999
      END

      FUNCTION FU(ANG,nc,XS,XB)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      real*8 XS(nc,3),XB(nc,3)
      COMMON/MATVAR/A(3,3),RTOL,NAT,ITER
      DIMENSION ANG(3)
      S1=SIN(ANG(1))
      S2=SIN(ANG(2))
      S3=SIN(ANG(3))
      C1=COS(ANG(1))
      C2=COS(ANG(2))
      C3=COS(ANG(3))
C     ang1 = theta
C     ang2 = phi
C     ang3 = cappa
      A(1,1)=C1*C2*C3-S2*S3
      A(1,2)=C1*S2*C3+C2*S3
      A(1,3)=-S1*C3
      A(2,1)=-C1*C2*S3-S2*C3
      A(2,2)=-C1*S2*S3+C2*C3
      A(2,3)=S1*S3
      A(3,1)=S1*C2
      A(3,2)=S1*S2
      A(3,3)=C1
      R=0.0d0
      DO 1 I=1,NAT
      TX=A(1,1)*XS(I,1)+A(1,2)*XS(I,2)+A(1,3)*XS(I,3)
      TY=A(2,1)*XS(I,1)+A(2,2)*XS(I,2)+A(2,3)*XS(I,3)
      TZ=A(3,1)*XS(I,1)+A(3,2)*XS(I,2)+A(3,3)*XS(I,3)
1     R=R+(TX-XB(I,1))*(TX-XB(I,1))+
     1    (TY-XB(I,2))*(TY-XB(I,2))+(TZ-XB(I,3))*(TZ-XB(I,3))
      FU=R
      RETURN
      END
c
c

      subroutine rvd(fn,nt0,nat,nm,rn,id,m,u,en)
      implicit none
      integer nat,nm,ia,id(*),ix,it,nt0
      real*8 rn(*),u(nt0,3),m(nt0,3),en(*)
      character*(*) fn
      open(30,file=fn,status='old')
      read(30,*)nm,nat,nat
      if(nm.gt.nt0)call report('Too many transitions')
      do 1 ia=1,nat
1     read(30,*)id(ia),(rn(ix+(ia-1)*3),ix=1,3)
      do 2 ia=1,2+nat*nm
2     read(30,*)
      read(30,300)(en(it),it=1,nm)
300   format(6f11.3)
      do 3 it=1,nm
3     read(30,*)(u(it,ix),ix=1,3),(m(it,ix),ix=1,3)
      close(30)
      return
      end

      subroutine makecenter(c,r,w,n,ic)
      implicit none
      integer*4 n,i,ii,ic(*)
      real*8 c(*),r(*),w(*),sum
      sum=0.0d0
      do 45 ii=1,n
45    sum=sum+w(ii)
      do 46 i=1,3
      c(i)=0.0d0
      do 47 ii=1,n
47    c(i)=c(i)+r(i+3*(ic(ii)-1))*w(ii)
46    c(i)=c(i)/sum
      return
      end

      SUBROUTINE READTTTQ(nm,n,alphad,Gpd,ad,ft,tw,ip,ntw,ntmax)
c     read normal mode polarizability derivatives
      IMPLICIT none
      INTEGER*4 nm,nq,LL,I,J,K,n,ip,ntmax,ntw(*)
      real*8 alphad(n,nm,3,3),Gpd(n,nm,3,3),Ad(n,nm,3,3,3),tw(n,ntmax)
      character*(*) ft
      open(15,file=ft,status='old')
      READ(15,*)
      READ(15,*)nq

      if(nq.ne.ntw(ip))call report('nq<>ntw')

      READ(15,*)
      READ(15,*)

      DO 1 LL=1,nq
      READ(15,*)tw(ip,LL),tw(ip,LL)
      DO 1 I=1,3
1     READ(15,*)alphad(ip,LL,I,1),(alphad(ip,LL,I,J),J=1,3)
      tw(ip,nq+1)=0.0d0

      READ(15,*)
      READ(15,*)

      DO 2 LL=1,nq
      READ(15,*)
      DO 2 I=1,3
2     READ(15,*)GPD(ip,LL,I,1),(GPD(ip,LL,I,J),J=1,3)

      READ(15,*)
      READ(15,*)

      DO 3 LL=1,nq
      READ(15,*)
      DO 3 I=1,3
      DO 3 J=1,3
3     READ(15,*)(Ad(ip,LL,I,J,K),K=1,2),(Ad(ip,LL,I,J,K),K=1,3)
c     first index electric
      CLOSE(15)
      RETURN
      END

      subroutine rot2(A,G,s,t,n1,nm,e)
      implicit none
      integer*4 s,t,i,j,n1,nm
      real*8 A(3,3),G(3,3),e(n1,nm,3,3),m(3,3)
      do 1 i=1,3
      do 1 j=1,3
1     m(i,j)=A(i,1)*e(t,s,1,j)+A(i,2)*e(t,s,2,j)+A(i,3)*e(t,s,3,j)
      do 2 i=1,3
      do 2 j=1,3
2     G(i,j)=A(j,1)*m(i,1)+A(j,2)*m(i,2)+A(j,3)*m(i,3)
      return
      end

      subroutine rot3(A,G,s,t,n1,nm,e)
      implicit none
      integer*4 s,t,i,j,k,n1,nm
      real*8 A(3,3),G(3,3,3),e(n1,nm,3,3,3),m(3,3,3)
      do 1 i=1,3
      do 1 j=1,3
      do 1 k=1,3
1     G(i,j,k)=A(i,1)*e(t,s,1,j,k)
     1        +A(i,2)*e(t,s,2,j,k)
     2        +A(i,3)*e(t,s,3,j,k)
      do 2 i=1,3
      do 2 j=1,3
      do 2 k=1,3
2     m(i,j,k)=A(j,1)*G(i,1,k)+A(j,2)*G(i,2,k)+A(j,3)*G(i,3,k)
      do 3 i=1,3
      do 3 j=1,3
      do 3 k=1,3
3     G(i,j,k)=A(k,1)*m(i,j,1)+A(k,2)*m(i,j,2)+A(k,3)*m(i,j,3)
      return
      end

      subroutine rotv(A,u,s,nt0,e)
      implicit none
      integer*4 s,ix,nt0
      real*8 A(3,3),u(3),e(nt0,3)
      do 53 ix=1,3
53    u(ix)=A(ix,1)*e(s,1)+A(ix,2)*e(s,2)+A(ix,3)*e(s,3)
      return
      end