program rroasum
      implicit none
      integer*4 nq,npx,ns0,i,ie,nn,ifr,nm,n0,ni0
      parameter (nn=10000)
      parameter (ns0=18,ni0=13)
      real*8 alr1(3,3),ali1(3,3),gr1(3,3),gi1(3,3),
     1ar1(3,3,3),ai1(3,3,3),gcr1(3,3),gci1(3,3),e(nn),
     1acr1(3,3,3),aci1(3,3,3),avr1(3,3),avi1(3,3),ff,inv(ni0),
     1THRC,wnm,wcm,cm,wau,wmin,wmax,fwhh,kelvin,ee(nn)
      logical lg,luseg,lusea,lstatic
      real*8,allocatable::alr(:,:,:),ali(:,:,:),gr(:,:,:),
     1gi(:,:,:),ar(:,:,:,:),ai(:,:,:,:),gcr(:,:,:),
     1gci(:,:,:),acr(:,:,:,:),aci(:,:,:,:),sr(:,:,:),
     1alrt(:,:,:),alit(:,:,:),grt(:,:,:),
     1git(:,:,:),art(:,:,:,:),ait(:,:,:,:),gcrt(:,:,:),
     1gcit(:,:,:),acrt(:,:,:,:),acit(:,:,:,:),
     1ls(:,:,:),gs(:,:,:),as(:,:,:,:),ws(:)

      character*80 s
      npx=3951
      THRC=1.0d-12
      wnm=532.0d0
      cm=219470.0d0
      wmin=50.0d0
      wmax=4000.0d0
      lg=.false.
      fwhh=10.0d0
      kelvin=300.0d0
      write(6,6001)
6001  format(' Reads polarizabilities from many files,',/,
     1       ' sums up and makes the spectra',/,/,
     1       '   Input: TTT.LST   list of files',/,
     1       '  Output: TTT.OUT   tensors',/,
     1       '          .prn      spectra ',/,
     1       '          INV.TXT   invariants ',/)
      call readopt('GRROA.OPT',luseg,lusea,lstatic,ff)
      ie=0
      open(9,file='TTT.LST',status='old')
1     read(9,800,end=99,err=99)s
800   format(a80)
      write(6,*)s
      ie=ie+1
      goto 1
 99   close(9)
      write(6,6003)ie
6003  format(i5,' files listed in TTT.LST')

      nm=0
      allocate(ls(1,3,3),gs(1,3,3),as(1,3,3,3),ws(1))
      if(lstatic)then
       call rst(1,nm,ls,gs,as,ws,0,ff)
       deallocate(ls,gs,as,ws)
       allocate(ls(nm,3,3),gs(nm,3,3),as(nm,3,3,3),ws(nm))
       call rst(nm,nm,ls,gs,as,ws,1,ff)
      endif


c     get maximal number of modes
      call  gee(nq,ie,e,ee)
c     maximal dimension in case some modes not assigned:
      n0=nq+nm

      allocate(alr(n0,3,3),ali(n0,3,3),gr(n0,3,3),gi(n0,3,3),
     1ar(n0,3,3,3),ai(n0,3,3,3),gcr(n0,3,3),gci(n0,3,3),
     1acr(n0,3,3,3),aci(n0,3,3,3))
      allocate(alrt(n0,3,3),alit(n0,3,3),grt(n0,3,3),
     1git(n0,3,3),art(n0,3,3,3),ait(n0,3,3,3),gcrt(n0,3,3),gcit(n0,3,3),
     1acrt(n0,3,3,3),acit(n0,3,3,3))
c     zero out total tensors:
      call zz1(n0,alrt,alit,grt,git,art,ait,gcrt,gcit,acrt,acit)

      ifr=1
      open(19,file='TTT.LST',status='old')
      do 2 i=1,ie
      read(19,800)s
      call zz1(n0,alr,ali,gr,gi,ar,ai,gcr,gci,acr,aci)
c     read this contribution:
      call ge(n0,nq,s,e,alr,ali,gr,gi,ar,ai,gcr,gci,acr,aci,wau,ifr)
c     if required, add the static part:
      if(lstatic)call adds(n0,nq,e,alr,gi,ar,gci,acr,nm,ls,gs,as,ws,ff,
     1wau)
c     add to total:
2     call aa1(n0,alrt,alit,grt,git,art,ait,gcrt,gcit,acrt,acit,
     1alr,ali,gr,gi,ar,ai,gcr,gci,acr,aci)

      wcm=ee(ifr)
      wau=wcm/cm
      wnm=1.0d7/wcm

      call initab
      allocate(sr(npx,2,ns0))
      sr=0.0d0
      open(78,file='INV.TXT')
      open(79,file='TTT.OUT')
      do 3 i=1,nq
      call bz2(n0,i,alrt,alit,grt,git,art,ait,gcrt,gcit,acrt,acit,
     1alr1,ali1,gr1,gi1,ar1,ai1,gcr1,gci1,acr1,aci1,avr1,avi1)
      call wrram3(wau,0.0d0,e(i),THRC,i,wnm,wmin,wmax,npx,lg,fwhh,
     1alr1,ali1,gr1,gi1,gcr1,gci1,ar1,ai1,acr1,aci1,avr1,avi1,sr,
     1luseg,lusea,inv)
3     call wzz(79,wau,0.0d0,e(i)/cm,alr1,ali1,gr1,gi1,gcr1,gci1,ar1,
     1ai1,acr1,aci1,avr1,avi1,inv)

      do 4 i=1,ns0
4     close(50+i)
      close(78)
      close(79)
      call c33(wmin,wmax,npx,sr,kelvin)
      end
c     ========================================
c     ==============================================================
      subroutine wzz(io,w,ei,ef,alr,ali,gr,gi,gcr,gci,ar,ai,
     1acr,aci,avr,avi,inv)
c     write ROA tensors into io
      implicit none
      integer*4 a,b,c,ni0,io
      parameter (ni0=13)
      real*8 w,ef,ei,CM,ECM,wcm,alr(3,3),ali(3,3),gr(3,3),gi(3,3),
     1gcr(3,3),gci(3,3),ar(3,3,3),ai(3,3,3),acr(3,3,3),aci(3,3,3),
     1avr(3,3),avi(3,3),inv(ni0)
      character*1 xyz(3)
      data xyz/'X','Y','Z'/
      character*9 si(ni0)
      data si/'a2       ','bs(a)2   ','ba(a)2   ','aG       ',
     1        'bs(G)2   ','ba(G)2   ','bs(A)2   ','ba(A)2   ',
     1        'aGc      ','bs(Gc)2  ','ba(Gc)2  ','bs(Ac)2  ',
     1        'ba(Ac)2  '/ 
      CM=219474.63d0
      ECM=(ef-ei)*CM
      WCM=w*CM
      write(io,701)ECM
701   format(/,' Energy =',f12.4,' cm^-1')
      write(io,702)WCM
702   format('   -> Omega =',f9.1,' cm^-1')
      write(io,706)'     ELECTRIC DIPOLE-ELECTRIC DIPOLE    '
706   format(1x,40(1h-),/,A40,/,1x,40(1h-),/,
     115x,'Real   < v -v >  Imag',11x,'Real   < r - r>  Imag')
      do 1 a=1,3
      do 1 b=1,3
1     write(io,707)xyz(b),'-',xyz(a),' ',avr(b,a),avi(b,a),
     1alr(b,a),ali(b,a)
707   format(2x,4a1,2E17.7,E15.7,E17.7)
      write(io,704)'     ELECTRIC DIPOLE-MAGNETIC DIPOLE    '
704   format(1x,40(1h-),/,A40,/,1x,40(1h-),/,15x,'Real',13x,'Imag')
      do 2 a=1,3
      do 2 b=1,3
2     write(io,705)xyz(b),'-',xyz(a),' ',gr(b,a),gi(b,a)
705   format(2x,4a1,2E17.7)
      write(io,704)'     MAGNETIC DIPOLE-ELECTRIC DIPOLE    '
      do 3 a=1,3
      do 3 b=1,3
3     write(io,705)xyz(b),'-',xyz(a),' ',gcr(b,a),gci(b,a)
      write(io,704)'   ELECTRIC DIPOLE-ELECTRIC QUADRUPOLE  '
      do 4 a=1,3
      do 4 b=1,3
      do 4 c=1,b
4     write(io,705)xyz(a),'-',xyz(c),xyz(b), ar(a,c,b), ai(a,c,b)
      write(io,704)'   ELECTRIC QUADRUPOLE-ELECTRIC DIPOLE  '
      do 5 a=1,3
      do 5 b=1,a
      do 5 c=1,3
5     write(io,705)xyz(b),xyz(a),'-',xyz(b),acr(c,b,a),aci(c,b,a)
      do 7 a=1,ni0
7     write(io,708)si(a),inv(a)
 708  format(3x,a9,' =',e14.6)
      return
      end
c     =====================================
      subroutine readopt(s,luseg,lusea,lstatic,ff)
      implicit none
      character*(*) s
      character*4 key
      logical lex,lusea,luseg,lstatic
      real*8 ff
      lusea=.true.
      luseg=.true.
      lstatic=.false.
      ff=1.0d0
      inquire(file=s,exist=lex)
      if(lex)then
       write(6,*)s
       open(9,file=s)
1      read(9,90,end=2,err=2)key
       write(6,90)key
90     format(a4)
       if(key.eq.'USEA')read(9,*)lusea
       if(key.eq.'USEG')read(9,*)luseg
       if(key.eq.'STAT')read(9,*)lstatic
       if(key.eq.'FACT')read(9,*)ff
       goto 1
2      close(9)
      endif
      return
      end
c     ========================================
      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end

      subroutine ap3(is,ns0,ECM,a,s,wrin,wrax,npx,fwhh,lglg)
      implicit none
      integer*4 ns0,npx,i,is
      real*8 ECM,wrin,wrax,fwhh,w,dw,sf,a(2),s(npx,2,ns0),t
      logical lglg
      dw=(wrax-wrin)/(npx-1)
      w=wrin-dw
      do 1 i=1,npx
      w=w+dw
      t=sf(fwhh,lglg,w,ECM)
      s(i,1,is)=s(i,1,is)+t*a(1)
1     s(i,2,is)=s(i,2,is)+t*a(2)
      return
      end
c     ========================================
      subroutine c33(wrin,wrax,npx,s,T)
      implicit none
      integer*4 npx,i,ns0,ii,ie
      parameter (ns0=18)
      real*8 wrin,wrax,s(npx,2,ns0),w,dw,tempcm,T
      character*11 st(ns0)
      data st/'ICP_0      ','ICPx_90    ','ICPz_90    ','ICPs_90    ',
     1        'ICPu_90    ','ICP_180    ','SCP_0      ','SCPx_90    ',
     1        'SCPz_90    ','SCPs_90    ','SCPu_90    ','SCP_180    ',
     1        'DCPI_0     ','DCPI_90    ','DCPI_180   ','DCPII_0    ',
     1        'DCPII_90   ','DCPII_180  '/
      tempcm=0.6950d0*T
      do 100 ii=1,ns0
      do 101 ie=len(st(ii)),1,-1
101   if(st(ii)(ie:ie).ne.' ')goto 102

102   open(45,file=st(ii)(1:ie)//'.ram.prn')
      dw=(wrax-wrin)/(npx-1)
      w=wrin-dw
      do 1 i=1,npx
      w=w+dw
1     write(45,4545)w,s(i,1,ii)/(1.0d0-exp(-w/tempcm))
 4545 format(f12.3,' ',g12.4)
      close(45)

      open(45,file=st(ii)(1:ie)//'.roa.prn')
      dw=(wrax-wrin)/(npx-1)
      w=wrin-dw
      do 2 i=1,npx
      w=w+dw
2     write(45,4545)w,s(i,2,ii)/(1.0d0-exp(-w/tempcm))

100   close(45)
      return
      end
c     ==============================================================
      subroutine wrram3(w,ei,ef,THRC,nt,EXCNM,wrin,wrax,npx,lglg,fwhh,
     1alr,ali,gtr,gti,gtcr,gtci,atr,ati,atcr,atci,avr,avi,sr,
     1luseg,lusea,inv)
c     EXCNM ... excitation frequency in nm
      implicit none
      integer*4 nt,npx,ns0,ni0,
     1ia,b,c,id,e,ii,is
      parameter (ns0=18,ni0=13)
c     ns0 : number of experimental setups
c     ni0 : number of invariants
      real*8 ef,ei,fwhh,THRC,AMU,BOHR,ECM,EXCNM,wrax,wrin,w,clight,CM,
     1YDY,YDX,sr(npx,2,ns0),gpisvejc,EXCA,roa1,ram1,tr,ti,a(2),eps,
     1inv(ni0),co(ni0+2,ns0),aaar,aaai,gtaar,gtaai,gtcaar,gtcaai,
     1alr(3,3),ali(3,3),gtr(3,3),gti(3,3),gtcr(3,3),gtci(3,3),
     1atr(3,3,3),ati(3,3,3),atcr(3,3,3),atci(3,3,3),
     1avr(3,3),avi(3,3),
     1alsr(3,3),alsi(3,3),gtsr(3,3),gtsi(3,3),gtcsr(3,3),gtcsi(3,3),
     1avsr(3,3),avsi(3,3),
     1alar(3,3),alai(3,3),gtar(3,3),gtai(3,3),gtcar(3,3),gtcai(3,3),
     1avar(3,3),avai(3,3),
     1ear(3,3),eapr(3,3),easr(3,3),eapsr(3,3),eaar(3,3),eapar(3,3),
     1eai(3,3),eapi(3,3),easi(3,3),eapsi(3,3),eaai(3,3),eapai(3,3),
     1eacr(3,3),eapcr(3,3),eacsr(3,3),eapcsr(3,3),eacar(3,3),
     1eapcar(3,3),eaci(3,3),eapci(3,3),eacsi(3,3),eapcsi(3,3),
     1eacai(3,3),eapcai(3,3)
      character*9 si(ni0)
      logical ldebug,lglg,luseg,lusea
      data si/'a2       ','bs(a)2   ','ba(a)2   ','aG       ',
     1        'bs(G)2   ','ba(G)2   ','bs(A)2   ','ba(A)2   ',
     1        'aGc      ','bs(Gc)2  ','ba(Gc)2  ','bs(Ac)2  ',
     1        'ba(Ac)2  '/ 
c             prefix RAM, prefix ROA
c                 a2  bs(a)2  ba(a)2      aG  bs(G)2 ba(G)2  bs(A)2
c             ba(A)2      aG  bs(G)2  ba(G)2  bs(A)2 ba(A)2
c     script:              ^     ^       ^       ^      ^
c
c    co: Raman prefactor, ROA prefactor,3 Raman invariants,
c        10 ROA invariants, all for ns0=18 spectral kinds
c    1 1 2 ICP(0o): Nafie
      data co/ 4.0d0,  8.0d0,
     1        45.0d0,  7.0d0,  5.0d0, 45.0d0,  7.0d0,  5.0d0, 1.0d0,
     1        -1.0d0,-45.0d0,  5.0d0, -5.0d0, -3.0d0,  1.0d0,
c    2 3 4 ICPx(90o): Nafie
     1         4.0d0,  2.0d0,
     1        45.0d0,  7.0d0,  5.0d0, 45.0d0,  7.0d0,  5.0d0,  1.0d0,
     1        -1.0d0,  0.0d0,  0.0d0,  0.0d0,  0.0d0,  0.0d0,
c    3 5 6 ICPz(90o): Nafie
     1         8.0d0,  4.0d0,
     1         0.0d0,  3.0d0,  5.0d0,  0.0d0,  3.0d0,  5.0d0, -1.0d0,
     1         1.0d0,  0.0d0,  0.0d0,  0.0d0,  0.0d0,  0.0d0,
c    4 7 8 ICP*(90o): (magic) Nafie
     1     13.3333d0,6.666d0,
     1         9.0d0,  2.0d0,  2.0d0,  9.0d0,  2.0d0,  2.0d0,  0.0d0,
     1         1.0d0,  0.0d0,  0.0d0,  0.0d0,  0.0d0,  0.0d0,
c    5 9 10 ICPu(90o): Nafie
     1         4.0d0,  4.0d0,
     1        45.0d0, 13.0d0, 15.0d0, 45.0d0, 13.0d0, 15.0d0, -1.0d0,
     1         1.0d0,  0.0d0,  0.0d0,  0.0d0,  0.0d0,  0.0d0,
c    6 11 12 ICP(180o): Nafie
     1         4.0d0,  8.0d0,
     1        45.0d0,  7.0d0,  5.0d0, 45.0d0,  7.0d0,  5.0d0,  1.0d0,
     1        -1.0d0, 45.0d0, -5.0d0,  5.0d0,  3.0d0, -1.0d0,
c    7 13 14 SCP(0o): nafie
     1         4.0d0,  8.0d0,
     1        45.0d0,  7.0d0,  5.0d0, 45.0d0, -5.0d0,  5.0d0, -3.0d0,
     1        -1.0d0,-45.0d0, -7.0d0, -5.0d0,  1.0d0,  1.0d0,
c    8 15 16 SCPx(90o): nafie
     1         2.0d0,  4.0d0, 
     1        45.0d0,  7.0d0,  5.0d0,  0.0d0,  0.0d0,  0.0d0,  0.0d0, 
     1         0.0d0,-45.0d0, -7.0d0, -5.0d0,  1.0d0,  1.0d0,
c    9 17 18 SCPz(90o): Nafie
     1         4.0d0,  8.0d0,
     1         0.0d0,  3.0d0,  5.0d0,  0.0d0,  0.0d0,  0.0d0,  0.0d0,
     1         0.0d0,  0.0d0, -3.0d0, -5.0d0, -1.0d0, -1.0d0,
c    10 19 20 SCP*(90o): Nafie
     1     13.3333d0,6.666d0,
     1         9.0d0,  2.0d0,  2.0d0,  0.0d0,  0.0d0,  0.0d0,  0.0d0,
     1         0.0d0, -9.0d0, -2.0d0, -2.0d0,  0.0d0,  0.0d0,
c    11 21 22 SCPu(90o): Nafie
     1         4.0d0,  4.0d0, 
     1        45.0d0, 13.0d0, 15.0d0,  0.0d0,  0.0d0,  0.0d0,  0.0d0,
     1         0.0d0,-45.0d0,-13.0d0,-15.0d0, -1.0d0, -1.0d0,
c    12 23 24 SCP(180o):Nafie
     1         4.0d0,  8.0d0, 
     1        45.0d0,  7.0d0,  5.0d0,-45.0d0,  5.0d0, -5.0d0,  3.0d0,
     1         1.0d0,-45.0d0, -7.0d0, -5.0d0,  1.0d0,  1.0d0,
c    13 25 26 DCPI(0o): Nafie
     1         4.0d0,  8.0d0,
     1        45.0d0,  1.0d0,  5.0d0, 45.0d0,  1.0d0,  5.0d0, -1.0d0,
     1        -1.0d0,-45.0d0, -1.0d0, -5.0d0, -1.0d0, +1.0d0,
c    14 27 28 DCPI(90o): Nafie
     1         2.0d0,  2.0d0, 
     1        45.0d0, 13.0d0, 15.0d0, 45.0d0, 13.0d0, 15.0d0, -1.0d0,
     1        +1.0d0,-45.0d0,-13.0d0,-15.0d0, -1.0d0, -1.0d0,
c    15 29 30 DCPI(180o): Nafie
     1        24.0d0, 16.0d0,
     1         0.0d0,  1.0d0,  0.0d0,  0.0d0,  3.0d0,  0.0d0,  1.0d0,
     1         0.0d0,  0.0d0, -3.0d0,  0.0d0,  1.0d0,  0.0d0,
c    16 31 32 DCPII(0o): Nafie
     1        24.0d0, 16.0d0,
     1         0.0d0,  1.0d0,  0.0d0,  0.0d0,  3.0d0,  0.0d0,  1.0d0,
     1         0.0d0,  0.0d0,  3.0d0,  0.0d0, -1.0d0,  0.0d0,
c    17 33 34 DCPII(90o): 
     1         2.0d0,  2.0d0,
     1        45.0d0, 13.0d0, 15.0d0, 45.0d0, 13.0d0, 15.0d0, -1.0d0,
     1         1.0d0, 45.0d0, 13.0d0, 15.0d0,  1.0d0, 1.0d0,
c    18 35 36 DCPII(180o): 
     1         4.0d0,  8.0d0,
     1        45.0d0,  1.0d0,  5.0d0, 45.0d0,  1.0d0,  5.0d0, -1.0d0,
     1        -1.0d0, 45.0d0,  1.0d0,  5.0d0,  1.0d0, -1.0d0/
c             prefix      a2  bs(a)2  ba(a)2      aG  bs(G)2 ba(G)2
c             bs(A)2  ba(A)2      aG  bs(G)2  ba(G)2  bs(A)2 ba(A)2
c     script:            ^          ^    ^       ^        ^     ^
c     invariants:
c     1 a2     =(1/9)Re (as_aa as*_bb)
c     2 bs(a)2 =(1/2)Re(3as_ab as*_ab-as_aa as*_bb)
c     3 ba(a)2 =(1/2)Re(3as_ab as*_ab-as_aa as*_bb)
c     4 aG     =(1/9)Im(as_aa Gs*_bb)
c     5 bs(G)2 =(1/2)Im(3as_ab Gs*_ab-as_aa Gs*_bb)
c     6 ba(G)2 =(3/2)Im(3aa_ab Ga*_ab)
c     7 bs(A)2 =(w/2)Im(i as_ab eas*_ab)
c                ea_ab=eps_adg A_d,gb
c     8 ba(A)2 =(w/2)Im(i aa_ab eaa*_ab)+i aa_ab eap_ab
c                eap_ab=eps_abg A_d,gd
c       script tensors:
c     9 aG     =(1/9)Im(as_aa Gs*_bb)
c    10 bs(G)2 =(1/2)Im(3as_ab Gs*_ab-as_aa Gs*_bb)
c    11 ba(G)2 =(3/2)Im(3aa_ab Ga*_ab)
c    12 bs(A)2 =(w/2)Im(i as_ab eas*_ab)
c    13 ba(A)2 =(w/2)Im(i aa_ab eaa*_ab)+i aa_ab eap_ab

      ldebug=.false.
      AMU=1822.0d0
      BOHR=0.529177d0
      EXCA=EXCNM*10.0d0
      ECM=ef-ei
      CM=219474.0d0
      clight=137.03599d0
      gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/EXCA

c     eap_ab = eps_abc A_d,cd
      call calcep(eapr ,atr )
      call calcep(eapi ,ati )
      call calcep(eapci,atci)
      call calcep(eapcr,atcr)
c      
c     ea_ab = eps_adc A_d,cb
      ear =0.0d0
      eai =0.0d0
      eacr=0.0d0
      eaci=0.0d0
      do 2 ia=1,3
      id=ia+1
      if(id.gt.3)id=1
      e=id+1
      if(e.gt.3)e=1
      do 2 b=1,3
      ear( ia,b)=atr (id,e,b)-atr (e,id,b)
      eai( ia,b)=ati (id,e,b)-ati (e,id,b)
      eacr(ia,b)=atcr(id,e,b)-atcr(e,id,b)
2     eaci(ia,b)=atci(id,e,b)-atci(e,id,b)

c     symmetric and antisymmetric tensor combinations
      call dsa(alr,ali,alsr,alsi,alar,alai)
      call dsa(avr,avi,avsr,avsi,avar,avai)
      call dsa(avr,avi,avsr,avsi,avar,avai)
      call dsa(gtr,gti,gtsr,gtsi,gtar,gtai)
      call dsa(gtcr,gtci,gtcsr,gtcsi,gtcar,gtcai)
      call dsa(ear  ,eai  ,easr  ,easi  ,eaar  ,eaai  )
      call dsa(eapr ,eapi ,eapsr ,eapsi ,eapar ,eapai )
      call dsa(eacr ,eaci ,eacsr ,eacsi ,eacar ,eacai )
      call dsa(eapcr,eapci,eapcsr,eapcsi,eapcar,eapcai)
      inv=0.0d0
c     a2: (1/9) Re (als_aa als*_bb)
      aaar=alsr(1,1)+alsr(2,2)+alsr(3,3)
      aaai=alsi(1,1)+alsi(2,2)+alsi(3,3)
      inv(1)=(aaar*aaar+aaai*aaai)/9.0d0
c     beta_s(alpha)2:(1/2)Re(3als_ab als*_ab-als_aa als*_bb)
      call abab(tr,ti,alsr,alsi,alsr,alsi)
      inv(2)=1.5d0*tr-4.5d0*inv(1)
c     beta_a(alpha)2:(3/2)Re(als_ab als*_ab)
      call abab(tr,ti,alar,alai,alar,alai)
      inv(3)=1.5d0*tr

      if(luseg)then
c     aG: (1/9) Im(als_aa gs*_bb)
      gtaar=gtsr(1,1)+gtsr(2,2)+gtsr(3,3)
      gtaai=gtsi(1,1)+gtsi(2,2)+gtsi(3,3)
      inv(4)=(aaai*gtaar-aaar*gtaai)/9.0d0
c     beta_s(G)2:(1/2)Im(3als_ab Gs*_ab-als_aa Gs*_bb)
      call abab(tr,ti,alsr,alsi,gtsr,gtsi)
      inv(5)=1.5d0*ti-4.5d0*inv(4)
c     beta_a(G)2: (3/2)Im(ala_ab Ga*_ab)
      call abab(tr,ti,alar,alai,gtar,gtai)
      inv(6)=1.5d0*ti
c     aGc:
      gtcaar=gtcsr(1,1)+gtcsr(2,2)+gtcsr(3,3)
      gtcaai=gtcsi(1,1)+gtcsi(2,2)+gtcsi(3,3)
      inv(9)=(aaai*gtcaar-aaar*gtcaai)/9.0d0
c     beta_s(Gc)2:
      call abab(tr,ti,alsr,alsi,gtcsr,gtcsi)
      inv(10)=1.5d0*ti-4.5d0*inv(9)
c     beta_a(Gc)2: (3/2)Im(ala_ab Gca*_ab)
      call abab(tr,ti,alar,alai,gtcar,gtcai)
      inv(11)=1.5d0*ti
      endif

      if(lusea)then
c     betasA2: (w/2)Im(i als_ab e_agd As*g,db)
      call abab(tr,ti,alsr,alsi,easr,easi)
      inv(7)=0.5d0*w*tr/clight
c     betaaA2: (w/2)Im(i[ala_ab[e_agd Aa*g,db+ e_abg Aa*d,gd)
      call abab(tr,ti,alar,alai,eaar,eaai)
      inv(8)=tr
      call abab(tr,ti,alar,alai,eapar,eapai)
      inv(8)=0.5d0*w*(inv(8)+tr)/clight
c     betasAc2: (w/2)Im(i als_ab e_agd Acs*g,db)
      call abab(tr,ti,alsr,alsi,eacsr,eacsi)
      inv(12)=0.5d0*w*tr/clight
c     betaaAc2:
      call abab(tr,ti,alar,alai,eacar,eacai)
      inv(13)=tr
      call abab(tr,ti,alar,alai,eapcar,eapcai)
      inv(13)=0.5d0*w*(inv(13)+tr)/clight
      endif


      if(ldebug.and.nt.eq.37)then
       write(6,*)'Omega:',w/clight

       call wre('Eacsr:',eacsr)
       eacr=0.0d0
       eaci=0.0d0
       do 61 ia=1,3
       do 61 b=1,3
       do 61 c=1,3
       do 61 id=1,3
       eacr(ia,b)=eacr(ia,b)+eps(ia,c,id)*atcr(c,id,b)
 61    eaci(ia,b)=eaci(ia,b)+eps(ia,c,id)*atci(c,id,b)
       call wre('Eacr:',eacr)
       call dsa(eacr ,eaci ,eacsr ,eacsi ,eacar ,eacai )
       call wre('Eacsr:',eacsr)

       do 6 ia=1,3
       do 6 is=1,3
6      write(6,*)ia,is,alr(ia,is),eacsr(ia,is),eacsi(ia,is)
      endif

      write(78,701)ECM
701   format(/,' Energy =',f12.4,' cm^-1')
      write(78,702)w*CM
702   format('   -> Omega =',f9.1,' cm^-1')
      do 7 ia=1,ni0
7     write(78,708)si(ia),inv(ia)
 708  format(3x,a9,' =',e14.6)

c     spectral intensities, a1 Raman, a2 ROA:      
      do 3 is=1,ns0
      a(1)=
     1co(1,is)*(co(3,is)*inv(1)+co(4,is)*inv(2)+co(5,is)*inv(3))
      a(2)=0.0d0
      do 4 ii=4,ni0
4     a(2)=a(2)+inv(ii)*co(2,is)*co(2+ii,is)
      if(dabs(a(1)).gt.THRC)then
       call ap3(is,ns0,ECM,a,sr,wrin,wrax,npx,fwhh,lglg)
       YDX=0.0d0
       YDY=0.0d0
       ram1=a(1)*gpisvejc
       roa1=a(2)*gpisvejc
       WRITE(50+is,3000)nt,ECM,YDX,YDY,ram1,0.0d0,0.0d0,0.0d0,0.0d0,
     1 roa1
3000   FORMAT(I7,f9.2,2f3.0,g12.4,4f3.0,g12.4,f3.0,2g12.4,$)
       write(50+is,*) 0.
      endif
3     continue

      return
      end

      subroutine abab(tr,ti,ar,ai,br,bi)
c     t = a_ab b*_ab
      implicit none
      integer*4 a
      real*8 tr,ti,ar(3,3),ai(3,3),br(3,3),bi(3,3)
      tr=0.0d0
      ti=0.0d0
      do 1 a=1,3
      tr=tr+ar(a,1)*br(a,1)+ai(a,1)*bi(a,1)
     1     +ar(a,2)*br(a,2)+ai(a,2)*bi(a,2)
     1     +ar(a,3)*br(a,3)+ai(a,3)*bi(a,3)
1     ti=ti-ar(a,1)*bi(a,1)+ai(a,1)*br(a,1)
     1     -ar(a,2)*bi(a,2)+ai(a,2)*br(a,2)
     1     -ar(a,3)*bi(a,3)+ai(a,3)*br(a,3)
      return
      end

      subroutine wre(s,e)
      character*(*) s
      real*8 e(3,3)
      write(6,'(a)')s
      write(6,609)e
609   format(3f12.6)
      return
      end

      subroutine dsa(tr,ti,tsr,tsi,tar,tai)
      implicit none
      real*8 tr(3,3),ti(3,3),tsr(3,3),tsi(3,3),tar(3,3),tai(3,3)
      integer*4 a,b
      do 1 a=1,3
      do 1 b=1,3
      tsr(a,b)=0.5d0*(tr(a,b)+tr(b,a))
      tsi(a,b)=0.5d0*(ti(a,b)+ti(b,a))
      tar(a,b)=0.5d0*(tr(a,b)-tr(b,a))
 1    tai(a,b)=0.5d0*(ti(a,b)-ti(b,a))
      return
      end

      subroutine bz2(nq,n,alr,ali,gr,gi,ar,ai,gcr,gci,acr,aci,
     1alr1,ali1,gr1,gi1,ar1,ai1,gcr1,gci1,acr1,aci1,avr1,avi1)
      implicit none
      integer*4 nq,n,i,j,k
      real*8 alr1(3,3),ali1(3,3),gr1(3,3),gi1(3,3),
     1ar1(3,3,3),ai1(3,3,3),gcr1(3,3),gci1(3,3),
     1acr1(3,3,3),aci1(3,3,3),avr1(3,3),avi1(3,3),
     1alr(nq,3,3),ali(nq,3,3),gr(nq,3,3),gi(nq,3,3),
     1ar(nq,3,3,3),ai(nq,3,3,3),gcr(nq,3,3),gci(nq,3,3),
     1acr(nq,3,3,3),aci(nq,3,3,3)
      do 1 i=1,3
      do 1 j=1,3
      alr1(i,j)=alr(n,i,j)
      ali1(i,j)=ali(n,i,j)
      gr1( i,j)=gr( n,i,j)
      gi1( i,j)=gi( n,i,j)
      gcr1(i,j)=gcr(n,i,j)
      gci1(i,j)=gci(n,i,j)
      avr1(i,j)=alr(n,i,j)
      avi1(i,j)=ali(n,i,j)
      do 1 k=1,3
      ar1( i,j,k)=ar( n,i,j,k)
      ai1( i,j,k)=ai( n,i,j,k)
      acr1(i,j,k)=acr(n,i,j,k)
1     aci1(i,j,k)=aci(n,i,j,k)
      return
      end


      subroutine initab
      implicit none
      integer*4 i,j,ns0
      parameter(ns0=18)
      character*11 st(ns0)
      data st/'ICP_0o     ','ICPx_90o   ','ICPz_90o   ','ICPs_90o   ',
     1        'ICPu_90o   ','ICP_180o   ','SCP_0o     ','SCPx_90o   ',
     1        'SCPz_90o   ','SCPs_90o   ','SCPu_90o   ','SCP_180o   ',
     1        'DCPI_0o    ','DCPI_90o   ','DCPI_180o  ','DCPII_0o   ',
     1        'DCPII_90o  ','DCPII_180o '/
      do 1 i=1,ns0
      do 2 j=len(st),1,-1
2     if(st(i)(j:j).ne.' ')goto 3
3     open(50+i,file=st(i)(1:j)//'.tab')
1     write(50+i,3000)st(i)
3000  format(20x,a11,/,'   n wavelength (cm)   Raman    ROA',/,80(1h-))
      return
      end

      subroutine gee(n,ie,em,ee)
c     get number of different transition energies in all files
      implicit none
      integer*4 n,nn,i,ie,nmax,is,n1,j,nex
      parameter (nn=10000)
      real*8 e,em(nn),e1(nn),ee(nn),wcm
      character*80 s,s80
c     nex ... number of excitation frequencies ee
c     n   ... number of transition frequencies em
      nex=0
      n=0
      em=0.0d0
      ee=0.0d0
      is=0
      e1=0.0d0
      nmax=0
      open(19,file='TTT.LST',status='old')
      do 4 i=1,ie
      read(19,800)s
 800  format(a80)
      open(9,file=s)
      n1=0
1     read(9,80,end=99,err=99)s80
80    format(a80)
      if(s80(1:10).eq.' Energy = ')then
       read(s80(10:21),*)e
       read(9,80,end=99,err=99)s80
       read(s80(14:22),*)wcm
       do 2 j=1,n
 2     if(dabs(e-em(j)).lt.1.0d-2)goto 3
       n=n+1
       if(n.gt.nn)call report('too many energies')
       em(n)=e
 3     continue 
       do 21 j=1,n1
 21    if(dabs(e-e1(j)).lt.1.0d-2)goto 31
       n1=n1+1
       if(n1.gt.nn)call report('too many energies')
       e1(n1)=e
 31    continue 
       do 22 j=1,nex
 22    if(dabs(wcm-ee(j)).lt.1.0d-2)goto 32
       nex=nex+1
       if(nex.gt.nn)call report('too many exc. energies')
       ee(nex)=wcm
 32    continue 
      endif
      goto 1
99    close(9)
      if(n1.gt.nmax)then
       nmax=n1
       is=ie
      endif
4     continue
      close(19)
      write(6,600)n,nmax,is,nex
 600  format(i4,' different transitions',/,
     1' maximum of ',i4,' in file number ',i4,/,
     1i4,' excitation frequencies',/)
      return
      end

      subroutine ge(n0,nq,s,e,alr,ali,gr,gi,ar,ai,gcr,gci,acr,aci,wau,
     1nwr)
      implicit none
      integer*4 n,nq,nwr,nw,nr,i,n0
      real*8 e(n0),alr(n0,3,3),ali(n0,3,3),gr(n0,3,3),gi(n0,3,3),
     1ar(n0,3,3,3),ai(n0,3,3,3),gcr(n0,3,3),gci(n0,3,3),eold,er,
     1acr(n0,3,3,3),aci(n0,3,3,3),wau,CM,wcm
      character*(*) s
      character*80 s80
      logical lf(5),lr
      CM=219470.0d0
      nw=0
      lf=.false.
      eold=0.0d0
      nr=0
      open(9,file=s)
1     read(9,80,end=99,err=99)s80
80    format(a80)
      if(s80(1:10).eq.' Energy = ')then
       read(s80(10:21),*)er
       read(9,80,end=99,err=99)s80
       read(s80(14:22),*)wcm
       if(nr.eq.0)eold=er
       if(nr.eq.0)nr=nr+1
       if(dabs(er-eold).lt.1.0d-2)then
        nw=nw+1
       else
        nw=1
        nr=nr+1
       endif
       lr=nw.eq.nwr
       if(lr)wau=wcm/CM
       n=0
       do 2 i=1,nq
2      if(dabs(er-e(i)).lt.1.0d-2)n=i
       if(n.eq.0)call report('energy not found')
       eold=er
      endif

c     if the required exc frequency is read:
      if(lr)then
      if(s80(6:36).eq.'ELECTRIC DIPOLE-ELECTRIC DIPOLE')then
       call r33(n0,n,alr,ali)
       lf(1)=.true.
      endif
      if(s80(6:36).eq.'ELECTRIC DIPOLE-MAGNETIC DIPOLE')then
       call r33(n0,n,gr,gi)
       lf(2)=.true.
      endif
      if(s80(6:36).eq.'MAGNETIC DIPOLE-ELECTRIC DIPOLE')then
       call r33(n0,n,gcr,gci)
       lf(3)=.true.
      endif
      if(s80(4:38).eq.'ELECTRIC DIPOLE-ELECTRIC QUADRUPOLE')then
       call r39(n0,n,ar,ai)
       lf(4)=.true.
      endif
      if(s80(4:38).eq.'ELECTRIC QUADRUPOLE-ELECTRIC DIPOLE')then
       call r93(n0,n,acr,aci)
       lf(5)=.true.
      endif
      endif
      goto 1
99    close(9)
c     a,ac:first index electric dipole
      if(     lf(1))write(6,*)' alpha'
      if(.not.lf(1))write(6,*)' alpha not found'
      if(     lf(2))write(6,*)' G '
      if(.not.lf(2))write(6,*)' G not found'
      if(     lf(3))write(6,*)' Gc '
      if(.not.lf(3))write(6,*)' Gc not found'
      if(     lf(4))write(6,*)' A '
      if(.not.lf(4))write(6,*)' A not found'
      if(     lf(5))write(6,*)' Ac '
      if(.not.lf(5))write(6,*)' Ac not found'
      if(lf(1).and.lf(2).and.lf(3).and.lf(4).and.lf(5))return
      end

      subroutine r93(nq,n,ar,ai)
c     quadrupole-dipole
      implicit none
      integer*4 n,nq,i,j,k
      real*8 ar(nq,3,3,3),ai(nq,3,3,3)
      character*80 s80
      read(9,*)
      read(9,*)
      do 1 j=1,3
      do 1 i=1,j
c     electric dipole, keep it in the first index:
      do 1 k=1,3
      read(9,80)s80
      read(s80(7:40),*)ar(n,k,i,j),ai(n,k,i,j)
      ar(n,k,j,i)=ar(n,k,i,j)
1     ai(n,k,j,i)=ai(n,k,i,j)
80    format(a80)
      return
      end

      subroutine r39(nq,n,ar,ai)
c     dipole-quadrupoe
      implicit none
      integer*4 n,nq,i,j,k
      real*8 ar(nq,3,3,3),ai(nq,3,3,3)
      character*80 s80
      read(9,*)
      read(9,*)
      do 1 i=1,3
      do 1 k=1,3
      do 1 j=1,k
      read(9,80)s80
      read(s80(7:40),*)ar(n,i,j,k),ai(n,i,j,k)
      ar(n,i,k,j)=ar(n,i,j,k)
1     ai(n,i,k,j)=ai(n,i,j,k)
80    format(a80)
      return
      end

      subroutine r33(nq,n,alr,ali)
      implicit none
      integer*4 j,i,n,nq
      real*8 alr(nq,3,3),ali(nq,3,3)
      character*80 s80
      read(9,*)
      read(9,*)
      do 1 j=1,3
      do 1 i=1,3
      read(9,80)s80
1     read(s80(7:40),*)alr(n,i,j),ali(n,i,j)
80    format(a80)
      return
      end

c     ==================================================================
      function sf(d,g,x,x0)
      implicit none
      logical g
      real*8 pi,spi,sf,d,x,x0,dd
      pi =3.14159265358979d0
      spi=1.77245385090552d0
      dd=((x-x0)/d)**2
      if(g)then
       if(dd.lt.20.0d0)then
        sf=exp(-dd)/d/spi
       else
        sf=0.0d0
       endif
      else
       if(dd.lt.1.0d10)then
        sf=1.0d0/d/(dd+1.0d0)/pi
       else
        sf=0.0d0
       endif
      endif
      return
      end
c     ==============================================================
      FUNCTION EPS(I,J,K)
      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

      subroutine calcep(e,a)
c     eap_ab = eps_abc A_d,cd
      implicit none
      real*8 e(3,3),a(3,3,3)
      e(1,1)=0.0d0
      e(2,2)=0.0d0
      e(3,3)=0.0d0
      e(1,2)= a(1,3,1)+a(2,3,2)+a(3,3,3)
      e(1,3)=-a(1,2,1)-a(2,2,2)-a(3,2,3)
      e(2,3)= a(1,1,1)+a(2,1,2)+a(3,1,3)
      e(2,1)=-e(1,2)
      e(3,1)=-e(1,3)
      e(3,2)=-e(2,3)
      return
      end

      subroutine zz1(nq,alrt,alit,grt,git,art,ait,gcrt,gcit,acrt,acit)
      integer*4 nq
      real*8 alrt(nq,3,3),alit(nq,3,3),grt(nq,3,3),
     1git(nq,3,3),art(nq,3,3,3),ait(nq,3,3,3),gcrt(nq,3,3),gcit(nq,3,3),
     1acrt(nq,3,3,3),acit(nq,3,3,3)
      alrt=0.0d0
      alit=0.0d0
      grt=0.0d0
      git=0.0d0
      art=0.0d0
      ait=0.0d0
      gcrt=0.0d0
      gcit=0.0d0
      acrt=0.0d0
      acit=0.0d0
      return
      end

      subroutine aa1(nq,alrt,alit,grt,git,art,ait,gcrt,gcit,acrt,acit,
     1alr,ali,gr,gi,ar,ai,gcr,gci,acr,aci)
      integer*4 nq
      real*8 alr(nq,3,3),ali(nq,3,3),gr(nq,3,3),gi(nq,3,3),
     1ar(nq,3,3,3),ai(nq,3,3,3),gcr(nq,3,3),gci(nq,3,3),
     1acr(nq,3,3,3),aci(nq,3,3,3),
     1alrt(nq,3,3),alit(nq,3,3),grt(nq,3,3),
     1git(nq,3,3),art(nq,3,3,3),ait(nq,3,3,3),gcrt(nq,3,3),gcit(nq,3,3),
     1acrt(nq,3,3,3),acit(nq,3,3,3)
      alrt=alrt+alr
      alit=alit+ali
      grt=grt+gr
      git=git+gi
      art=art+ar
      ait=ait+ai
      gcrt=gcrt+gcr
      gcit=gcit+gci
      acrt=acrt+acr
      acit=acit+aci
      return
      end

      subroutine rst(nm0,nm,ls,gs,as,ws,ic,ff)
      implicit none
      integer*4 nm0,nm,ic,i,ix,iy,iz
      real*8 ls(nm0,3,3),gs(nm0,3,3),as(nm0,3,3,3),ws(nm0),ff
      open(9,file='FILE.Q.TTT',status='old')
      read(9,*)
      read(9,*)nm
      if(ic.eq.0)then
       close(9)
       return
      endif
      read(9,*)
      read(9,*)
      do 9 i=1,nm
      read(9,*)ws(i),ws(i)
      do 9 ix=1,3
9     read(9,*)ls(i,ix,1),(ls(i,ix,iy),iy=1,3)
      read(9,*)
      read(9,*)
      do 91 i=1,nm
      read(9,*)
      do 91 ix=1,3
c     magnetic index iy, last
91    read(9,*)gs(i,ix,1),(gs(i,ix,iy),iy=1,3)
      read(9,*)
      read(9,*)
      do 92 i=1,nm
      read(9,*)
      do 92 ix=1,3
      do 92 iy=1,3
c     el index ix
92    read(9,*)(as(i,ix,iy,iz),iz=1,2),(as(i,ix,iy,iz),iz=1,3)
      close(9)
      write(6,600)ff
600   format('static polarizability read, factor used',e12.3)
      return
      end

      subroutine adds(nn,nq,e,lr,gi,ar,gci,acr,nm,ls,gs,as,ws,ff,
     1w)
      integer*4 nq,nm,i,j,ix,iy,iz,nn,ji,na
      real*8 ls(nm,3,3),gs(nm,3,3),as(nm,3,3,3),ws(nm),
     1e(nn),lr(nn,3,3),gi(nn,3,3),ar(nn,3,3,3),gci(nn,3,3),
     1acr(nn,3,3,3),t,ff,CM,f,d,w,clight,g
      integer*4 ,allocatable::asi(:)
      clight=137.03599d0
      CM=219474.0d0
      allocate(asi(nn))
      asi=0
c     na ... number of unassigned freuqencies
      na=0
      do 1 i=1,nm
      f=ff/dsqrt(2.0d0*ws(i)/CM)
      g=f*clight/w
      ji=0
      d=1.0d9
      do 2 j=1,nq
      t=dabs(ws(i)-e(j))
      if(t.lt.1.0d0.and.t.lt.d.and.asi(j).eq.0)then
       ji=j
       d=t
      endif
2     continue
      if(ji.eq.0)then
       write(6,600)i,ws(i)
600    format(' Fraquency',i6,f12.2,' cm-1 not assigned')
       na=na+1
       ji=nq+na
       e(ji)=ws(i)
       do 4 ix=1,3
       do 4 iy=1,3
       lr( ji,ix,iy   )=ls(i,ix,iy   )*f
c      gs is supposedly static limit, G'/w, so here it must be
c      multiplied by w:
       gi( ji,ix,iy   )=-gs(i,ix,iy   )*f*w
       gci(ji,ix,iy   )= gs(i,ix,iy   )*f*w
       do 4 iz=1,3
       ar( ji,ix,iy,iz)=as(i,ix,iy,iz)*g
4      acr(ji,ix,iy,iz)=as(i,ix,iy,iz)*g
      else
       asi(ji)=1
       do 3 ix=1,3
       do 3 iy=1,3
       lr( ji,ix,iy   )=lr( ji,ix,iy   )+ls(i,ix,iy   )*f
       gi( ji,ix,iy   )=gi( ji,ix,iy   )-gs(i,ix,iy   )*f*w
       gci(ji,ix,iy   )=gci(ji,ix,iy   )+gs(i,ix,iy   )*f*w
       do 3 iz=1,3
       ar( ji,ix,iy,iz)=ar( ji,ix,iy,iz)+as(i,ix,iy,iz)*g
3      acr(ji,ix,iy,iz)=acr(ji,ix,iy,iz)+as(i,ix,iy,iz)*g
      endif
1     continue
      nq=nq+na
      return
      end