program polderdip2
      implicit none
      integer*4 na,nl,ne,nat,n,mm,np,natf,ie,iep
      real*8 gau,nmmin,w,dnm,g8
      real*8,allocatable::m(:,:),mi(:,:),we(:),wi(:),
     1al(:,:,:,:),g(:,:,:,:),a(:,:,:,:,:),
     1ald(:,:,:),gd(:,:,:),ad(:,:,:,:),aldv(:,:,:),
     1aldw(:,:,:,:),gpdw(:,:,:,:),adw(:,:,:,:,:),alvdw(:,:,:,:)
      logical lint,lcon
      write(6,600)
 600  format(' Resonance Raman optical activity -',
     1' tensor calculation',/,/,
     1'    Input: TTT.LST       list of FILE.TTT files',/,
     1'           TD.TAB        table of electronic energies',/,
     1'           POLDERDIP.OPT options',/,/,
     1'   Output: FILE.TTT, FILEW.TTT  polarizability derivatives')
      call readopt(nmmin,np,lint,lcon,w,gau,dnm,g8)
      na=nl('TTT.LST')
      if(na.eq.0)call report('No files in TTT.LST')
      ne=nl('TD.TAB')-4
      if(ne.le.0)call report('No frequencies in TD.TAB')
      write(6,*)ne,' energies in TD.TAB'
      allocate(we(ne))
      call rw(we,ne)
      nat=natf('TTT.LST')
      write(6,*)nat,' atoms'
      n=3*nat
      allocate(wi(na),al(na,n,3,3),g(na,n,3,3),a(na,n,3,3,3))
      call rp(na,n,al,g,a,wi)
c     correlation matrix:
      if(lcon)then
       mm=ne+1
      else
       mm=ne
      endif
      if(na.lt.mm)call report('Too few polarizabilities')
      allocate(m(mm,mm),mi(mm,mm))
      call dom(m,mm,ne,na,wi,we,lcon,g8)
      call INV(m,mi,mm,1.0d-10)
      write(6,*)'M-1:'
      do 3 ie=1,mm
 3    write(6,60)(mi(ie,iep),iep=1,mm)
   60 format(12e10.4)

c     do derivatives for one frequency:
      allocate(ald(n,3,3),aldv(n,3,3),gd(n,3,3),ad(n,3,3,3))
      call dt(na,ne,n,mm,ald,aldv,gd,ad,w,mi,we,lcon,gau,al,g,a,wi,g8)
      call WRITETTT('FILE.TTT',n,ald,aldv,gd,ad,w)

      if(lint)then
c      do derivatives for frequency interval:
       allocate(aldw(np,n,3,3),gpdw(np,n,3,3),adw(np,n,3,3,3),
     1 alvdw(np,n,3,3))
       call dw(ne,n,mm,np,na,ald,aldv,gd,ad,mi,we,lcon,gau,al,g,
     1 a,aldw,gpdw,adw,alvdw,nmmin,dnm,wi,g8)
       call WRITETTTW(np,n,aldw,alvdw,gpdw,adw,nmmin,dnm)
       call wrscale(np,nmmin,dnm)
      endif


      end

      subroutine wrscale(np,nmmin,dnm)
      implicit none
      real*8 nm,dnm,nmmin
      integer*4 np,i
      open(13,file='xnm.prn')
      nm=nmmin-dnm
      do 1 i=1,np
      nm=nm+dnm
1     write(13,121)nm
121   format(f12.2)
      close(13)
      write(6,*)'xnm.prn scale written'
      return
      end

      subroutine dw(ne,n,mm,np,na,ald,aldv,gd,ad,mi,we,lcon,gau,al,g,
     1a,aldw,gpdw,adw,alvdw,nmmin,dnm,wi,g8)
      implicit none
      integer*4 ne,n,mm,np,na,iw,ix,jx,l
      real*8 CM,ald(n,3,3),aldv(n,3,3),gd(n,3,3),ad(n,3,3,3),
     1wau,mi(mm,mm),we(*),gau,al(na,n,3,3),g(na,n,3,3),a(na,n,3,3,3),
     1aldw(np,n,3,3),gpdw(np,n,3,3),adw(np,n,3,3,3),alvdw(np,n,3,3),
     1nmmin,dnm,nm,wi(*),g8
      logical lcon
      CM=219474.0d0
      nm=nmmin-dnm
      do 1 iw=1,np
      nm=nm+dnm
      wau=1.0d7/nm/CM
      call dt(na,ne,n,mm,ald,aldv,gd,ad,wau,mi,we,lcon,gau,al,g,a,wi,g8)
      do 1 ix=1,3
      do 1 jx=1,3
      do 1 l=1,n
      aldw( iw,l,ix,jx)= ald(l,ix,jx)
      alvdw(iw,l,ix,jx)= ald(l,ix,jx)
      gpdw( iw,l,ix,jx)= gd( l,ix,jx)
      adw(  iw,l,ix,jx,1)=ad(l,ix,jx,1)
      adw(  iw,l,ix,jx,2)=ad(l,ix,jx,2)
 1    adw(  iw,l,ix,jx,3)=ad(l,ix,jx,3)
      return
      end

      subroutine dt(na,ne,n,mm,ald,aldv,gpd,ad,w,mi,we,lcon,gau,
     1al,g,aa,wi,g8)
      implicit none
      integer*4 ne,n,a,b,c,l,i,mm,na
      real*8 ald(n,3,3),gpd(n,3,3),aldv(n,3,3),ad(n,3,3,3),gau,wi(*),
     1mi(mm,mm),al(na,n,3,3),g(na,n,3,3),aa(na,n,3,3,3),w,we(*),g8
      real*8 ,allocatable::at(:,:,:,:),atd(:,:,:)
      logical lcon
      allocate(at(na,n,3,3),atd(n,3,3))
      call alg(na,ne,n,mm,mi,al,ald,we,wi,lcon,gau,w,g8)
      aldv=ald
      call alg(na,ne,n,mm,mi,g,gpd,we,wi,lcon,gau,w,g8)
      do 1 c=1,3
      do 2 a=1,3
      do 2 b=1,3
      do 2 l=1,n
      do 2 i=1,na
2     at(i,l,a,b)=aa(i,l,a,b,c)
      call alg(na,ne,n,mm,mi,at,atd,we,wi,lcon,gau,w,g8)
      do 1 a=1,3
      do 1 b=1,3
      do 1 l=1,n
1     ad(l,a,b,c)=atd(l,a,b)
      return
      end

      subroutine alg(na,ne,n,mm,mi,al,ald,we,wi,lcon,gau,w,g)
      implicit none
      integer*4 na,ne,n,mm,a,b,l,e,i,ep
      real*8 al(na,n,3,3),ald(n,3,3),we(ne),wi(na),w,mi(mm,mm),g,
     1gau
      real*8,allocatable:: v(:),AA(:)
      logical lcon
      allocate(v(mm),AA(mm))
      ald=0.0d0
      do 1 a=1,3
      do 1 b=1,3
      do 1 l=1,n
      v=0.0d0
      do 11 e=1,ne
      do 11 i=1,na
11    v(e)=v(e)+al(i,l,a,b)/((we(e)**2-wi(i)**2)**2+(g*we(e))**2)
      if(lcon)then
       do 111 i=1,na
111    v(mm)=v(mm)+al(i,l,a,b)
      endif
      do 12 e=1,mm
      AA(e)=0.0d0
      do 12 ep=1,mm
12    AA(e)=AA(e)+mi(e,ep)*v(ep)
c     fitted element at frequency w and bandwidth gau:
      if(lcon)ald(l,a,b)=ald(l,a,b)+AA(mm)
      do 1 e=1,ne
1     ald(l,a,b)=ald(l,a,b)+AA(e)/((we(e)**2-w**2)**2+(gau*we(e))**2)
      return
      end

      subroutine readopt(nmmin,np,lint,lcon,w,gau,dnm,g8)
      implicit none
      real*8 width,excnm,CM,w,gau,nmmin,nmmax,
     1dnm,widthnm,g8,g8nm
      integer*4 np
      logical lint,lex,lcon
      character*80 s80,st
      CM=219474.0d0
c     electronic band width in nm
      widthnm=25.0d0
c     excitation frequency in nm
      excnm=532.0d0
c     excitation wavelength limit for scan, in nm
      nmmin=500.0d0
c     excitation wavelength limit for scan, in nm
      nmmax=600.0d0
c     number of exc frequencies on the scan
      np=10
c     scan an interval of exc. frequencies:
      lint=.false.
c     add a constant term:
      lcon=.true.
c     limit width
      g8nm=5.0d0
      inquire(file='POLDERDIP.OPT',exist=lex)
      if(lex)then
       open(9,file='POLDERDIP.OPT')
1111   read(9,90,end=991,err=991)s80
90     format(a80)
       if(s80(1:5).eq.'WIDTH')read(9,*)widthnm
       if(s80(1:5).eq.'EXCNM')read(9,*)excnm
       if(s80(1:5).eq.'INTER')read(9,*)lint
       if(s80(1:5).eq.'NMMIN')read(9,*)nmmin
       if(s80(1:5).eq.'NMMAX')read(9,*)nmmax
       if(s80(1:5).eq.'NPOIN')read(9,*)np
       if(s80(1:5).eq.'LCONS')read(9,*)lcon
       if(s80(1:4).eq.'G8NM')read(9,*)g8nm
       if(s80(1:3).eq.'END')goto 991
       backspace 9
       read(9,90)st
       if(st.eq.s80)call report('unknown option '//st)
       goto 1111
991    close(9)
      endif

c     excitation in au:
      w=1.0d7/excnm/CM
c     width in cm-1:
      width=1.0d7/excnm**2*widthnm
c     width in au:
      gau=width/CM
      g8   =1.0d7/excnm**2*g8nm/CM
      if(np.gt.1)then
       dnm=(nmmax-nmmin)/dble(np-1)
      else
       dnm=0.0d0
      endif
      if(.not.lint)np=1

      write(6,60)widthnm,gau,excnm,w,nmmin,nmmax,lint,lcon,g8nm,np
 60   format(' Width     : ',f8.2,' nm   (',f14.6,' au)',/,
     1       ' Excitation: ',f8.2,' nm   (',f14.6,' au)',/,
     1       ' nmmin     : ',f8.2,' nmmax:',f14.6,' nm ',/,
     1       ' interval  : ',l8  ,' const:',l14  ,'    ',/,
     1       ' G8nm      : ',f8.2,' Np   :',i14  ,/)

      return
      end

      subroutine dom(m,mm,ne,na,wi,we,lcon,g)
      implicit none
      integer*4 na,ne,i,ie,mm,iep
      real*8 m(mm,mm),g,w1,w2,wi(na),we(ne)
      logical lcon
      m=0.0d0
      do 1 i=1,na
      do 1 ie=1,ne
      w1=    (we( ie)**2-wi(i)**2)**2+(g*we(ie ))**2
      do 1 iep=1,ne
      w2=w1*((we(iep)**2-wi(i)**2)**2+(g*we(iep))**2)
1     m(ie,iep)=m(ie,iep)+1.0d0/w2
      if(lcon)then
       do 2 ie=1,ne
       do 2 i=1,na
       w1=(we(ie)**2-wi(i)**2)**2+(g*we(ie))**2
       m(ie,mm)=m(ie,mm)+1.0d0/w1
2      m(mm,ie)=m(mm,ie)+1.0d0/w1
       m(mm,mm)=dble(na)
      endif
      write(6,*)'M:'
      do 3 ie=1,mm
 3    write(6,60)(m(ie,iep),iep=1,mm)
   60 format(12e10.4)
      return
      end

      subroutine rp(na,n,al,g,a,wi)
      implicit none
      integer*4 na,n,i,ii,ix,jx
      real*8 wi(na),al(na,n,3,3),g(na,n,3,3),a(na,n,3,3,3),w,CM
      character*80 s80
      real*8,allocatable::ALPHAS(:,:,:),GS(:,:,:),AS(:,:,:,:)
      CM=219474.0d0
      allocate(ALPHAS(n,3,3),GS(n,3,3),AS(n,3,3,3))
      OPEN(9,FILE='TTT.LST',STATUS='OLD')
      do 1 i=1,na
      read(9,90)s80
 90   format(a80)
      CALL READTTT(s80,n,ALPHAS,AS,GS,w)
      wi(i)=w
      do 1 ii=1,n
      do 1 ix=1,3
      do 1 jx=1,3
      al(i,ii,ix,jx)=ALPHAS(ii,ix,jx)
      g( i,ii,ix,jx)=GS(    ii,ix,jx)
      a( i,ii,ix,jx,1)=AS( ii,ix,jx,1)
      a( i,ii,ix,jx,2)=AS( ii,ix,jx,2)
 1    a( i,ii,ix,jx,3)=AS( ii,ix,jx,3)
      close(9)
      write(6,60)na
 60   format(i4,' polarizabilities at:')
      do 2 i=1,na
  2   write(6,61)i,1.0d7/(wi(i)*CM),wi(i)*CM,wi(i)
61    format(i3,f8.2,' nm,',f9.1,' cm-1,',f10.6,' au')
      return
      end

      SUBROUTINE READTTT(s,n,ALPHA,A,G,w)
      IMPLICIT none
      integer*4 n,NATC,NAT,I,L,IX,IIND,J,K
      real*8 ALPHA(n,3,3),G(n,3,3),A(n,3,3,3),w
      character*(*) s
      NAT=n/3
      OPEN(15,FILE=s,STATUS='OLD')
      READ(15,*)
      READ(15,1515)NATC,w
 1515 format(i4,14x,2x,f11.6)
      if(NAT.NE.NATC)call report('Wrong number of atoms in '//s)
      READ(15,*)
      READ(15,*)
      DO 1 I=1,3
      READ(15,*)
      DO 1 L=1,NAT
      DO 1 IX=1,3
      IIND=3*(L-1)+IX
1     READ(15,*)(ALPHA(IIND,I,J),J=1,2),(ALPHA(IIND,I,J),J=1,3)
      READ(15,*)
      READ(15,*)
      DO 2 I=1,3
      READ(15,*)
      DO 2 L=1,NAT
      DO 2 IX=1,3
      IIND=3*(L-1)+IX
2     READ(15,*)(G(IIND,I,J),J=1,2),(G(IIND,I,J),J=1,3)
      READ(15,*)
      READ(15,*)
      DO 3 I=1,3
      DO 3 J=1,3
      READ(15,*)
      DO 3 L=1,NAT
      DO 3 IX=1,3
      IIND=3*(L-1)+IX
3     READ(15,*)(A(IIND,I,J,K),K=1,2),(A(IIND,I,J,K),K=1,3)
      CLOSE(15)
      RETURN
      END

      function natf(s)
      implicit none
      character*(*) s
      character*80 f
      integer*4 natf,n
      open(9,file=s)
      read(9,90)f
 90   format(a80)
      close(9)
      open(9,file=f)
      read(9,*)
      read(9,*)n
      close(9)
      natf=n
      return
      end

      subroutine rw(we,ne)
      implicit none
      real*8 we(*),l,CM
      integer*4 ne,i
      CM=219474.0d0
      open(9,file='TD.TAB',status='old')
      read(9,*)
      read(9,*)
      read(9,*)
      do 1 i=1,ne
      read(9,*)l,l
      write(6,60)i,l,1.0d7/l,1.0d7/l/cm
60    format(i3,f8.2,' nm,',f9.1,' cm-1,',f10.6,' au')
1     we(i)=1.0d7/l/CM
      read(9,*)
      close(9)
      return
      end

      function nl(f)
c     number of lines in a file
      implicit none
      character*(*) f
      integer*4 nl,n
      open(9,file=f)
      n=0
1     read(9,*,end=99,err=99)
      n=n+1
      goto 1
99    close(9)
      nl=n
      return
      end
      
      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end

      SUBROUTINE INV(A,AI,N,TOL)                                     
      IMPLICIT none
      integer*4 N,oo,ii,jj,iw,io,kk,i2,j2
      REAL*8 TOL, A(N,N),AI(N,N),w                           
      real*8, allocatable ::E(:,:)

      allocate(E(N,2*N))

      DO 1  ii=1,N                                                 
      DO 1  jj=1,N                                                 
      e(ii,jj)=a(ii,jj)
      E(II,JJ+N)=0.0D0                                             
1     if (ii.EQ.jj)e(ii,jj+N)=1.0D0                                

      DO 2 ii=1,N-1                                                
      iw=ii
      if (DABS(e(iw,iw)).LT.TOL) then                                 
       DO 3 io=iw+1,N                                              
       oo=IO
       if (io.GT.N)oo=io-N                                    
3      if (DABS(e(oo,iw)).GE.TOL) goto 11                             
       call report('Inverse cannot be done')
11     CONTINUE                                                       
       DO 4 kk=1, 2*N                                              
       w=e(iw,kk)                                                     
       e(iw,kk)=e(oo,kk)                                              
4      e(oo,kk)=w                                                     
      ENDIF                                                           
      DO 2 jj=ii+1,N                                               
      DO 6 kk=ii+1, 2*N
6     e(jj,kk)=e(jj,kk)-e(ii,kk)*e(jj,ii)/e(ii,ii)                    
2     e(jj,ii)=0.0D0                                                  

      DO 7  ii=1, N-1
      i2=N-ii+1                                                    
      DO 7 jj=1,i2-1                                                  
      j2=i2-jj                                                        
      DO 9 kk=1, N                                                 
9     e(j2,kk+N)=e(j2,kk+N)-e(i2,kk+N)*e(j2,i2)/e(i2,i2)
7     e(j2,i2)=0.0D0
 
      DO 10 ii=1,N                                                 
      DO 10 jj=1,N                                                
10    AI(ii,jj)=e(ii,jj+N)/e(ii,ii)                                

      deallocate(E)
      RETURN                                                          
      END                     

      SUBROUTINE WRITETTTW(np,N3,ald,avd,gd,ad,nmmin,dnm)
      IMPLICIT none
      INTEGER*4 np,n3,L,IX,ii,J,I,K,NAT,iw
      real*8 ald(np,N3,3,3),gd(np,N3,3,3),ad(np,N3,3,3,3),fr,
     1ev,avd(np,N3,3,3),evcm,nmmin,dnm,nm

      evcm=8065.73d0
      NAT=N3/3
      OPEN(2,FILE='FILEW.TTT')
      nm=nmmin-dnm
      do 4 iw=1,np
      nm=nm+dnm
      fr=1.0d7/nm
      ev=fr/evcm
      WRITE(2,2000)NAT,fr,ev,nm
2000  FORMAT(' ROA tensors, cartesian derivatives',/,
     1I4,' atoms, freq. ',f11.2,' cm-1',f11.5,' ev,',f10.2,' nm'/,
     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
      ii=3*(L-1)+IX
1     WRITE(2 ,2001)L,IX,(ald(iw,ii,I,J),J=1,3)
2001  FORMAT(I5,1H ,I1,6g15.7)

      WRITE(2,2004)
2004  FORMAT(' The electric dipole magnetic dipole polarizability:',/,
     1       ' Atom/x    jx(Bx)       jy(By)       jz(Bz), v-form:')
      DO 2 I=1,3
      WRITE(2,2003)I
2003  FORMAT(' G(',I1,',J):')
      DO 2 L=1,NAT
      DO 2 IX=1,3
      ii=3*(L-1)+IX
2     WRITE(2 ,2001)L,IX,(GD(iw,ii,I,J),J=1,3),(GD(iw,ii,I,J),J=1,3)

      WRITE(2,2005)
2005  FORMAT(' The electric dipole electric quadrupole polarizability:',
     2/,     ' Atom/x    kx           ky           kz')
      DO 3 I=1,3
      DO 3 J=1,3
      WRITE(2,2006)I,J
2006  FORMAT(' A(',I1,',',I1,',K):')
      DO 3 L=1,NAT
      DO 3 IX=1,3
      ii=3*(L-1)+IX
3     WRITE(2,2007)L,IX,(ad(iw,ii,I,J,K),K=1,3),L,IX,I,J
2007  FORMAT(I5,1H ,I1,3g15.7,' ',4i3)

      write(2,*)
      write(2,*)' Alpha v:'
      DO 4 I=1,3
      WRITE(2,2002)I
      DO 4 L=1,NAT
      DO 4 IX=1,3
      ii=3*(L-1)+IX
4     WRITE(2 ,2001)L,IX,(avd(iw,ii,I,J),J=1,3)

      CLOSE(2)
      write(6,*)'FILEW.TTT written'
      
      RETURN
      END

      SUBROUTINE WRITETTT(fn,N3,ald,alv,GD,ad,fr)
      IMPLICIT none
      INTEGER*4 n3,L,IX,ii,J,I,K,NAT
      real*8 ald(N3,3,3),GD(N3,3,3),ad(N3,3,3,3),fr
      real*8 alv(N3,3,3)
      character*(*) fn

      NAT=N3/3
      OPEN(2,FILE=fn)
      WRITE(2,2000)NAT,fr
2000  FORMAT(' ROA tensors, cartesian derivatives',/,
     1I4,' atoms, freq. ',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
      ii=3*(L-1)+IX
1     WRITE(2 ,2001)L,IX,(ald(ii,I,J),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)',
     1       ' gradient alternative:')
      DO 2 I=1,3
      WRITE(2,2003)I
2003  FORMAT(' G(',I1,',J):')
      DO 2 L=1,NAT
      DO 2 IX=1,3
      ii=3*(L-1)+IX
2     WRITE(2 ,2001)L,IX,(GD(ii,I,J),J=1,3)

      WRITE(2,2005)
2005  FORMAT(' The electric dipole electric quadrupole polarizability:',
     2/,     ' Atom/x    kx           ky           kz')
      DO 3 I=1,3
      DO 3 J=1,3
      WRITE(2,2006)I,J
2006  FORMAT(' A(',I1,',',I1,',K):')
      DO 3 L=1,NAT
      DO 3 IX=1,3
      ii=3*(L-1)+IX
3     WRITE(2,2007)L,IX,(ad(ii,I,J,K),K=1,3),L,IX,I,J
2007  FORMAT(I5,1H ,I1,3g15.7,' ',4i3)

      write(2,*)
      write(2,*)'Alpha v:'
      DO 4 I=1,3
      WRITE(2,2002)I
      DO 4 L=1,NAT
      DO 4 IX=1,3
      ii=3*(L-1)+IX
4     WRITE(2 ,2001)L,IX,(alv(ii,I,J),J=1,3)

      write(6,*)fn//' written'

      CLOSE(2)
      
      RETURN
      END