program prepmd2
      implicit none
      character*80 s80
      integer*4 NAT,NPOT,nx,nten,nttt,nfc,MA,NEXTRA,i,j,ii,jj,
     1IERR,M,N3,nats,ig,inow,ni0,skind
      parameter (ni0=1000)
      integer*4 iso(ni0)
      integer*4,allocatable::list(:),plist(:),ty(:),icoord(:,:),
     1is,ie,igeo,np,ISOT
      real*8,allocatable::a(:,:),ai(:,:),e(:,:),
     1fi(:,:),r(:),bv(:,:),tm(:,:,:),xv(:,:),xi(:)
      logical lex,loff,lzero,ltest
      real*8 wmin,wmax,riso(ni0),EXCA,temp,hw
      character*1 pt

      write(6,500)
500   format(' Prepares spectral-constrained MD computation',/,
     1       '    Input:  X.LST   list of .X files',/,
     1       '            *.X     *.X geometry files',/,
     1       '            FC.LST  list of *.FC files',/,
     1       '            *.FC    *.FC force field files',/,
     1       '            TEN.LST list of *.TEN files',/,
     1       '            *.TEN   *.TEN  tensors',/,
     1       '            TTT.LST list of *.TTT files',/,
     1       '            *.TTT   *.TTT  tensors',/)

      inquire(file='PREPMD2.OPT',exist=lex)
      if(.not.lex)call report('PREPMD2.OPT missing')

      igeo=0
      NAT=0
      NPOT=0
      lzero=.true.
      ltest=.true.
      skind=8
c     LZERO  .. sum of coeff = 0
c     NAT   
c     nat                  number of solute atoms
c     i1,i2, ....  ,inat   and their list
c     NPOT
c     npot                 number of atoms where potential is measured
c     i1,i2, ....  ,inpot  and their list
c     LOFF ... fit also off-diagonal FF
c     NEXTRA ... number of other conditions, saved for later used
c     IGEO ... generate spectrum for this geometry with the parameters:
c              WMIN
c              WMAX
c              NP
c              HW FWHM
c      EXCA  excitation freq in A
c      SKIND spectral kind: 1 absorption
c                          2 VCD
c                          3 Raman
c                          4 ROA
c      pt = L Lorentzian bands
c           G gaussian bands
c      temp temperature in K
c      ltest ... test the fit
c     
      temp=300.0d0
      hw=10.0d0
      pt='L'
      SKIND=8
      EXCA=4880.0d0
      NEXTRA=0
      loff=.false.
      wmin=100.d0
      wmax=4000.d0
      np=3901
      ISOT=0
      open(9,file='PREPMD2.OPT')
2     read(9,90,end=89,err=89)s80
90    format(a80)
      write(6,90)s80
      if(s80(1:4).eq.'LOFF')read(9,*)loff
      if(s80(1:4).eq.'LZER')read(9,*)lzero
      if(s80(1:4).eq.'LTES')read(9,*)ltest
      if(s80(1:4).eq.'WMIN')read(9,*)wmin
      if(s80(1:4).eq.'ISOT')then
       read(9,*)ISOT
       if(ISOT.gt.ni0)call report('Too many isotopic substitutions')
       read(9,*)(iso(i),i=1,ISOT)
       read(9,*)(riso(i),i=1,ISOT)
      endif
      if(s80(1:4).eq.'BAND')read(9,'(a)')pt
      if(s80(1:4).eq.'TEMP')read(9,*)temp
      if(s80(1:2).eq.'HW')read(9,*)hw
      if(s80(1:4).eq.'WMAX')read(9,*)wmax
      if(s80(1:4).eq.'EXCA')read(9,*)EXCA
      if(s80(1:4).eq.'NPOI')read(9,*)np
      if(s80(1:4).eq.'IGEO')read(9,*)igeo
      if(s80(1:4).eq.'SKIN')read(9,*)skind
      if(s80(1:3).eq.'NAT')then
       read(9,*)NAT
       write(6,*)NAT
       allocate(list(NAT))
       read(9,*)(list(i),i=1,NAT)
      endif
      if(s80(1:4).eq.'NPOT')then
       read(9,*)NPOT
       write(6,*)NPOT
       allocate(plist(NPOT))
       read(9,*)(plist(i),i=1,NPOT)
      endif
      goto 2
89    close(9)
      if(NAT .eq.0)call report('NAT  = 0')
      if(NPOT.eq.0)call report('NPOT = 0')
      
      call tf('X.LST',nx)
      if(nx.eq.0)call report('nothing to do')
      call tf('FC.LST',nfc)
      call tf('TEN.LST',nten)
      call tf('TTT.LST',nttt)
      if(nfc.ne.nx)call report('nfc <> nx')
      if(nten.ne.nx)call report('nten <> nx')
      if(nttt.ne.nx)call report('nttt <> nx')
c     
c     potential at each site and geometry
      allocate(fi(NPOT,nx))
      call tz(fi,NPOT,nx)
c
c     go through geometries
      open(9,file='X.LST')
      do 6 ig=1,nx
      read(9,90)s80

      call rdx(s80,nats,ty,r,0)
      allocate(r(3*nats),ty(nats))
      call rdx(s80,nats,ty,r,1)
c
c     potential at each site
      call setpot(NPOT,plist,r,nats,ty,fi,nx,ig,list,NAT)

6     deallocate(r,ty)
      close(9)
c     
c     initialize and calculate correlation matrix
c     MA dimension of the matrix:
c      NPOT   number of the potential sites
c      1      p0 ... constant part of the parameter
c      1      lambda ... for the boundary condition
c      NEXTRA saved for later use
      if(lzero)then
       MA=NPOT+NEXTRA+2
      else
       MA=NPOT+NEXTRA+1
      endif
      allocate(a(MA,MA),ai(MA,MA),e(MA,2*MA))
      do 3 i=1,MA
      do 3 j=1,MA
3     a(i,j)=0.0d0
      a(1,1)=dble(nx)
      if(lzero)then
       do 4 i=1,NPOT
       a(i+1,MA)=1.0d0
4      a(MA,i+1)=1.0d0
      endif
      do 10 i=2,NPOT+1
      do 101 ii=1,nx
101   a(1,i)=a(1,i)+fi(i-1,ii)
10    a(i,1)=a(1,i)
      do 11 i=2,NPOT+1
      do 11 j=i,NPOT+1
      do 111 ii=1,nx
111   a(i,j)=a(i,j)+fi(i-1,ii)*fi(j-1,ii)
11    a(j,i)=a(i,j)
      call wmtr(A,MA,MA,'correlation matrix',1)
      call INV(MA,a,ai,MA,e,IERR)
      if(IERR.eq.0)then
       call wmtr(ai,MA,MA,'correlation matrix inverted',1)
      else
       call report('inversion error')
      endif
c
c     define local coordinate systems for solute atoms
c     based on the first geometry
      open(9,file='X.LST')
      read(9,90)s80
      close(9)
      call rdx(s80,nats,ty,r,0)
      allocate(r(3*nats),ty(nats))
      call rdx(s80,nats,ty,r,1)
      allocate(icoord(NAT,2))
      call defcoord(icoord,NAT,r,list)
      deallocate(r,ty)
c
c     for each geometry, load in parameters and transfer to internal
c     from the internal parameters, make vectors B to solve A X = B
c     p_i = p_i0 + sum(j) c_ij f_j + Xk gk
c     p_i ... spectroscopic parameter, i = 1 .. M
c     f_j ... electrostatic potential on atom j, j=1...NPOT
c     Xik ... extra conditions to be defined later, k=1..NEXTRA
c     nf  ... dimension of the 
c         
      N3=3*NAT
      NEXTRA=0
c     M: number of parameters
c          FF         APT   AAT   alpha  G      A
      if(loff)then
       M = (N3*(N3+1))/2+N3*(3+3+9+9+27)
c      M = N3*(N3+1)/2  FC          1             ...    N3*(N3+1)/2
c           3*NAT * 3   APT         N3**2+1       ...    N3*(N3+1)/2+ 3*N3
c           3*NAT * 3   AAT         N3**2+3*N3 +1 ...    N3*(N3+1)/2+ 6*N3
c           3*NAT * 9   alpha_der   N3**2+6*N3 +1 ...    N3*(N3+1)/2+15*N3
c           3*NAT * 9   G_der       N3**2+15*N3+1 ...    N3*(N3+1)/2+24*N3
c           3*NAT *27   A_der       N3**2+24*N3+1 ...    N3*(N3+1)/2+51*N3
      else
       M = N3*( 1+3+3+9+9+27)
c      M =  3*NAT       FC          1             ...    N3   
c           3*NAT * 3   APT         N3   +1       ...    N3   + 3*N3
c           3*NAT * 3   AAT         N3   +3*N3 +1 ...    N3   + 6*N3
c           3*NAT * 9   alpha_der   N3   +6*N3 +1 ...    N3   +15*N3
c           3*NAT * 9   G_der       N3   +15*N3+1 ...    N3   +24*N3
c           3*NAT *27   A_der       N3   +24*N3+1 ...    N3   +51*N3
      endif
      allocate(bv(MA,M))
      call tz(bv,MA,M)
      write(6,*)
      allocate(xv(MA,M),xi(MA))
      inquire(file='XV.PRM',exist=lex)
      if(lex)then
       write(6,*)'XV.PRM found, skip its calculation'
       open(9,file='XV.PRM')
       do 20 i=1,M
20     read(9,*)XV(1,i),(XV(j,i),j=1,MA)
       close(9)
      else
       write(6,*)'reading tensors, making sums:'
       open(9,file='X.LST')
       open(10,file='FC.LST')
       open(11,file='TEN.LST')
       open(12,file='TTT.LST')
       do 15 ig=1,nx
       read(9,90)s80
       call iss(s80,is,ie)
       write(6,*)'Geometry ',ig,':'//s80(is:ie)
       call rdx(s80(is:ie),nats,ty,r,0)
       allocate(r(3*nats),ty(nats),tm(NAT,3,3))
       call rdx(s80(is:ie),nats,ty,r,1)
       write(6,*)nats,' atoms'

c      define transformation matrices:
       call deftm(NAT,list,icoord,r,tm)
       write(6,*)'transformation matrix done'

       inow=0
       read(10,90)s80
       call iss(s80,is,ie)
       write(6,*)s80(is:ie)
       call trff(s80(is:ie),
     1 nats,tm,list,NAT,bv,inow,loff,M,MA,fi,NPOT,nx,ig)
       write(6,*)'FF done'

       read(11,90)s80
       call iss(s80,is,ie)
       write(6,*)s80(is:ie)
       call trafoten(s80(is:ie),
     1 nats,list,tm,NAT,r,bv,inow,M,MA,fi,NPOT,nx,ig)
       write(6,*)'TEN done'

       read(12,90)s80
       call iss(s80,is,ie)
       write(6,*)s80(is:ie)
       call trafottt(s80(is:ie),
     1 nats,list,tm,NAT,r,bv,inow,M,MA,fi,NPOT,nx,ig)
       write(6,*)'TTT done'
       write(6,*)
       deallocate(r,ty,tm)
15     continue
       close(9)
       close(10)
       close(11)
       close(12)

c      solve AX = B, i.e. X = A-1 B
      
c      save X:
       do 1 i=1,M
       do 17 j=1,MA
       xi(j)=0.0d0
       do 17 jj=1,MA
17     xi(j)=xi(j)+ai(j,jj)*bv(jj,i)
       do 1 j=1,MA
1      xv(j,i)=xi(j)
      
       open(9,file='XV.PRM')
       do 18 i=1,M
       write(9,9001)i
9001   format(i9,$)
       do 181 j=1,MA
181    write(9,9000)XV(j,i)
9000   format(' ',g13.6,$)
18     write(9,*)
       close(9)
      endif
c     
c     test fit
      if(ltest)call tests(xv,MA,M,nx,loff,fi,NPOT,NAT,list)
c
c     generate spectrum for a given geometry
      if(igeo.ne.0)then
       open(9,file='X.LST')
       do 19 i=1,igeo
19     read(9,90)s80
       close(9)
       call iss(s80,is,ie)
       call rdx(s80,nats,ty,r,0)
       write(6,*)'Generating spectrum for '//s80(is:ie)
       allocate(r(3*nats),ty(nats))
       call rdx(s80,nats,ty,r,1)
       write(6,*)'Coordinates read'
       call gsp(nats,ty,r,NPOT,plist,list,NAT,icoord,ISOT,iso,riso,
     1 wmin,wmax,EXCA,skind,NP,pt,temp,hw,loff,xv,MA,M)
      endif

      end

      subroutine gsp(nats,ty,r,NPOT,plist,list,NAT,icoord,ISOT,iso,riso,
     1wmin,wmax,EXCA,skind,NP,pt,temp,hw,loff,xv,MA,M)
      implicit none
      integer*4 nats,ty(*),NPOT,NAT,plist(*),list(*),icoord(NAT,2),
     1ISOT,skind,iso(*),NP,MA,M,ii,jj,i,ix,j,jx,kx,lx,ia,id,ig,iia,
     1inow,ip,jja,ja,L,jxstart,N
      real*8 r(*),EXCA,riso(*),wmax,wmin,temp,hw
      real*8,allocatable::f(:,:),fi(:,:),p(:,:,:),a(:,:,:),al(:,:,:),
     1g(:,:,:),aa(:,:,:,:),rr(:,:),u(:,:),XM(:),FRQ(:),fu(:,:),
     1pu(:,:,:),au(:,:,:),alu(:,:,:),gu(:,:,:),aau(:,:,:,:),xmu(:),
     1st(:,:),SMAT(:,:),tm(:,:,:),ru(:)
      integer*4,allocatable::tyu(:)
      character*1 pt
      logical loff
      real*8 xv(MA,M),BOHR,SUM,EPS
c     
c     potential at this geometry
      allocate(fi(NPOT,1))
      call tz(fi,NPOT,1)
      call setpot(NPOT,plist,r,nats,ty,fi,1,1,list,NAT)
c     
c     FF, APT, AAT, al, Gp, A:
      allocate(f(3*nats,3*nats),p(nats,3,3),a(nats,3,3),al(3*nats,3,3),
     1g(3*nats,3,3),aa(3*nats,3,3,3),rr(3,nats),u(3,nats),XM(nats))
      N=3*NAT
      allocate(fu(N,N),pu(NAT,3,3),au(NAT,3,3),ru(N),
     1alu(N,3,3),gu(N,3,3),aau(N,3,3,3),xmu(NAT),tyu(NAT))
      inow=0
      if(loff)then
       do 2 i=1,NAT
       ia=list(i)
       do 2 ix=1,3
       ii=ix+3*(ia-1)
       do 2 j=i,NAT
       ja=list(j)
       if(j.eq.i)then
        jxstart=ix
       else
        jxstart=1
       endif
       do 2 jx=jxstart,3
       jj=jx+3*(ja-1)
       inow=inow+1
       f(ii,jj)=xv(1,inow)
       do 21 ip=1,NPOT
21     f(ii,jj)=f(ii,jj)+xv(ip+1,inow)*fi(ip,1)
2      f(jj,ii)=f(ii,jj)
      else
       do 3 i=1,NAT
       ia=list(i)
       do 3 ix=1,3
       ii=ix+3*(ia-1)
       inow=inow+1
       f(ii,ii)=xv(1,inow)
       do 3 ip=1,NPOT
 3     f(ii,ii)=f(ii,ii)+xv(ip+1,inow)*fi(ip,1)
      endif
      write(6,*)'ff'

      do 5 i=1,NAT
      ia=list(i)
      do 5 ix=1,3
      do 5 jx=1,3
      inow=inow+1
      p(ia,ix,jx)=xv(1,inow)
      do 5  ip=1,NPOT
5     p(ia,ix,jx)=p(ia,ix,jx)+xv(ip+1,inow)*fi(ip,1)
      write(6,*)'apt'

      do 1 i=1,NAT
      ia=list(i)
      do 1 ix=1,3
      do 1 jx=1,3
      inow=inow+1
      a(ia,ix,jx)=xv(1,inow)
      do 1  ip=1,NPOT
1     a(ia,ix,jx)=a(ia,ix,jx)+xv(ip+1,inow)*fi(ip,1)
      write(6,*)'aat'
      
      do 6 i=1,NAT
      ia=list(i)
      do 6 ix=1,3
      ii=ix+3*(ia-1)
      do 6 jx=1,3
      do 6 kx=1,3
      inow=inow+1
      al(ii,jx,kx)=xv(1,inow)
      do 6  ip=1,NPOT
6     al(ii,jx,kx)=al(ii,jx,kx)+xv(ip+1,inow)*fi(ip,1)
      write(6,*)'al'

      do 7 i=1,NAT
      ia=list(i)
      do 7 ix=1,3
      ii=ix+3*(ia-1)
      do 7 jx=1,3
      do 7 kx=1,3
      inow=inow+1
      g(ii,jx,kx)=xv(1,inow)
      do 7  ip=1,NPOT
7     g(ii,jx,kx)=g(ii,jx,kx)+xv(ip+1,inow)*fi(ip,1)
      write(6,*)'g'

      do 8 i=1,NAT
      ia=list(i)
      do 8 ix=1,3
      ii=ix+3*(ia-1)
      do 8 jx=1,3
      do 8 kx=1,3
      do 8 lx=1,3
      inow=inow+1
      aa(ii,jx,kx,lx)=xv(1,inow)
      do 8  ip=1,NPOT
8     aa(ii,jx,kx,lx)=aa(ii,jx,kx,lx)+xv(ip+1,inow)*fi(ip,1)
      write(6,*)'a'
c
c     tranformation to common coordinate system
      allocate(tm(NAT,3,3))
      call deftm(NAT,list,icoord,r,tm)
      call trfftoCOM(f,nats,tm,list,NAT)
      call trafote1(p ,nat,nats,list,tm,1)
      call trafote1(a ,nat,nats,list,tm,1)
      call trafott1(al,nat,nats,list,tm,1)
      call trafott1(g ,nat,nats,list,tm,1)
      call trafott2(aa,nat,nats,list,tm,1)
      write(6,*)'com'
c
c     correct origin dependence
      do 9 i=1,nats
      do 9 ix=1,3
9     rr(ix,i)=r(ix+3*(i-1))      
      CALL TTTT(3*nats,al,aa,g,2,rr)

      BOHR=0.52917705993d0
      DO 10 L=1,nats
      DO 10 I=1,3
10    u(I,L)=rr(I,L)/BOHR
      DO 2201 L=1,nats
      DO 2201 J=1,3
      DO 2201 I=1,3
      SUM=0.0d0
      DO 2202 ID=1,3
      DO 2202 IG=1,3
2202  SUM=SUM-0.25d0*EPS(I,IG,ID)*u(IG,L)*P(L,J,ID)
2201  A(L,J,I)=A(L,J,I)+SUM
      WRITE(6,*)' AAT shifted to common origin'
c
c     define masses
      call INPTRY(XM,nats,ty,iso,riso,ISOT)
      WRITE(6,*)' masses defined'

c     rewrite to solute-only fields:
      do 114 i=1,NAT
      ia=list(i)
      xmu(i)=XM(ia)
114   tyu(i)=ty(ia)

      do 115 i=1,NAT
      ia=list(i)
      do 115 ix=1,3
      ii =ix+3*(i -1)
      iia=ix+3*(ia-1)
      do 115 j=1,NAT
      ja=list(j)
      do 115 jx=1,3
      jj =jx+3*(j -1)
      jja=jx+3*(ja-1)
115   fu(ii,jj)=f(iia,jja)

      do 116 i=1,NAT
      ia=list(i)
      do 116 ix=1,3
      do 116 jx=1,3
      pu(i,ix,jx)=p(ia,ix,jx)
116   au(i,ix,jx)=a(ia,ix,jx)

      do 111 i=1,NAT
      ia=list(i)
      do 111 ix=1,3
      ii=ix+3*(i-1)
      iia=ix+3*(ia-1)
      ru(ii)=r(iia)
      do 111 jx=1,3
      do 111 kx=1,3
      alu(ii,jx,kx   )=al(iia,jx,kx   )
      gu( ii,jx,kx   )=g( iia,jx,kx   )
      do 111 lx=1,3
111   aau(ii,jx,kx,lx)=aa(iia,jx,kx,lx)

      open(20,file='FILE.FC')
      CALL WRITEFF(3*NAT,fu)
      close(20)
      WRITE(6,*)'FILE.FC'
      CALL WRITETTT(N,alu,aau,gu,'FILE.TTT')

      do 113 ii=1,N
      do 113 jj=1,N
113   fu(ii,jj)=fu(ii,jj)*4.359828d0/0.5291772d0/0.5291772d0
      
      open(9,file='FILE.X')
      write(9,*)'test geom'
      write(9,*)NAT
      do 112 i=1,NAT
112   write(9,900)tyu(i),(ru(ix+3*(i-1)),ix=1,3),(0,ix=1,7),0.0d0
900   format(I3,3F12.6,7(1x,i1),f4.1)
      close(9)
      write(6,*)'FILE.X'

c     diagonalize FF
      allocate(FRQ(N),SMAT(N,N))
      call new4(NAT,fu,xmu,SMAT,FRQ,1,tyu,ru)
      allocate(st(2,N))
      if(skind.ne.3.and.skind.ne.8)then
       write(6,*)'skind ',skind,' not implemented'
       stop
      else
c      get Raman/ROA intensities
       call new6(SMAT,wmin,wmax,alu,gu,aau,N,EXCA,FRQ,st,3)
c      make and save spectrum
       call mkprn(N,st,3,NP,wmin,wmax,pt,temp,hw)
       call new6(SMAT,wmin,wmax,alu,gu,aau,N,EXCA,FRQ,st,8)
c      make and save spectrum
       call mkprn(N,st,8,NP,wmin,wmax,pt,temp,hw)
      endif
      return
      end

      subroutine mkprn(N,st,ic,npoint,wsmin,wsmax,pt,temp,hw)
      implicit none
      integer*4 ic,npoint,i,i1,N
      real*8 wsmax,wsmin,dw,dd,hw,xx,const,wo,
     1sf,temp,st(2,N)
      real*8,allocatable::s(:),t(:),w(:)
      character*1 pt
      logical LG
c    1' Arguments: N number of frequencies
c                 st frequencies and intensities
c    1'            ic type - type of spectrum (column # in name)',/,
c    1'            npoint - number of spectral points',/,
c    1'            wsmin - wmin',/,
c    1'            wsmax -wmax',/,
c    1'            6 G - Gaussian peaks instead of Lorentzians',/,
c    1'            7 T - temperature',/,
c    1'            8 width - full width at half height')
      allocate(s(npoint),t(N),w(N))
      LG=pt.eq.'g'.or.pt.eq.'G'
      do 1 i1=1,N
      w(i1)=st(1,i1)
1     t(i1)=st(2,i1)
      
      do 3 i=1,npoint
3     s(i)=0.0d0
      dw=(wsmax-wsmin)/dble(npoint-1)
      if(LG)then
       dd=hw/2.0d0/dsqrt(log(2.0d0))
      else
       dd=hw/2.0d0
      endif

      do 4 i=1,npoint
      xx=wsmin+dw*dble(i-1)
      do 4 i1=1,N
      wo=w(i1)
      const=t(i1)
c     abs:
      if(ic.eq.1)const=108.7d0*t(i1)*wo
c     vcd:
      if(ic.eq.2)const=435.0d0*t(i1)*wo
c     mcd - B-terms:
      if(ic.eq.20)const=-5.98442d-3*t(i1)*wo
c     Raman:
      if(ic.eq.3)const=t(i1)/wo/(1.0d0-exp(-(1.44d0*wo/temp)))
c     roa:
      if(ic.eq.8)const=t(i1)/wo/(1.0d0-exp(-(1.44d0*wo/temp)))
c     ord:
      if(ic.eq.30)const=t(i1)*8.35d18
4     s(i)=s(i)+sf(dd,LG,xx,wo,ic)*const 
c
      if(ic.eq.1)call wrspec2('d.prn',npoint,wsmin,dw,s)
      if(ic.eq.2)call wrspec2('r.prn',npoint,wsmin,dw,s)
      if(ic.eq.3)call wrspec2('ram.prn',npoint,wsmin,dw,s)
      if(ic.eq.8)call wrspec2('roa.prn',npoint,wsmin,dw,s)
      end
c     ==================================================================
      subroutine wrspec2(sn,nw,wmin,dw,spec)
      implicit none
      character*(*) sn
      integer*4 nw,iw
      real*8 wmin,dw,spec(*),w
      open(8,file=sn)
      do 1 iw=1,nw
      w=wmin+dble(iw-1)*dw
1     write(8,805)w,spec(iw)
805   format(f12.2,' ',g25.12)
      close(8)
      write(6,*)sn,' written'
      return
      end
c     ==================================================================
      function sf(d,g,x,x0,ic)
      implicit none
      logical g
      integer*4 ic
      real*8 pi,spi,sf,d,x,x0,dd,ab
      pi=4.0d0*atan(1.0d0)
      spi=dsqrt(pi)
      dd=((x-x0)/d)**2
c     write(6,*)d,x,x0,g
      if(ic.eq.30)then
c      ORD, real parts, only Lorentz possible:
       ab=(x-x0)*(x+x0)*x**2*x0**2
       sf=ab/(ab**2+(d**2)**2)
      else
c      dispersion spectra, imaginary part
       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
      endif
      return
      end
c     
      subroutine new6(S,wmin,wmax,ALPHA,G,A,N,EXCA,FRQ,st,skind)
      IMPLICIT REAL*8 (A-H,O-Z)
      real*8 st(2,N),S(N,N),ALPHA(N,3,3),A(N,3,3,3),G(N,3,3),FRQ(*)
      integer*4 skind,N,i,ix,jx,kx
      real*8,allocatable::ali(:,:,:),ai(:,:,:,:),gi(:,:,:)
      allocate(ali(N,3,3),ai(N,3,3,3),gi(N,3,3))
      do 1 i=1,N
      do 1 ix=1,3
      do 1 jx=1,3
      ali(i,ix,jx)=ALPHA(i,ix,jx)
      gi( i,ix,jx)=G(    i,ix,jx)
      do 1 kx=1,3
1     ai(i,ix,jx,kx)=A(i,ix,jx,kx)
      write(6,*)'1x x x (x) ',ali(1,1,1),gi(1,1,1),ai(1,1,1,1)
      CALL TRTEN(N,S,ali,ai,gi)
      CALL DORAMAN(N,ali,gi,ai,EXCA,wmin,wmax,FRQ,st,skind,1)
      return 
      END

      SUBROUTINE DORAMAN(N,ALPHA,G,A,EXCA,wmin,wmax,E,st,skind,ipr)
c
c     supposedly G/W is supplied
c
      IMPLICIT REAL*8 (A-H,O-Z)
      implicit integer*4 (i-n)
      DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),E(*)
      real*8 st(2,N)
      integer*4 skind
c
      if(ipr.eq.1)then
       OPEN(4,FILE='ROA.TAB')
       WRITE(4,3304)
3304   FORMAT(
     1 'MODE   FREQ    Ramx (180)  Ramy    Ramtot   CID90',
     2 '     CID-X   CID180  DEP   ROA180 RamanDCPI',/,
     3 '       cm-1       A^4/AMU                    x10^4',
     4 '     x10^4     x10^4      A^5/AMUx1e4 au.10^4')
       write(4,3002)
3002   FORMAT(80(1H-))
      endif
c
      DO 1 IQ=1,N
      ECM=E(IQ)
      IF(ABS(ECM).GT.0.01d0)THEN
      if(wmin.eq.0.0d0.or.ECM.gt.wmin)then
      if(wmax.eq.0.0d0.or.ECM.lt.wmax)then
       SAL0=0.0d0
       SAL1=0.0d0
       SAG0=0.0d0
       SAG1=0.0d0
       SA1=0.0d0
       SPAL=0.0d0
       DO 3 I=1,3
       SPAL=SPAL+ALPHA(IQ,I,I)
       DO 3 J=1,3
       SAL0=SAL0+ALPHA(IQ,I,I)*ALPHA(IQ,J,J)
       SAL1=SAL1+ALPHA(IQ,I,J)*ALPHA(IQ,I,J)
       SAG0=SAG0+ALPHA(IQ,I,I)*    G(IQ,J,J)
       SAG1=SAG1+ALPHA(IQ,I,J)*    G(IQ,I,J)
       DO 3 K=1,3
       DO 3 L=1,3
3      SA1=SA1+EPS(I,K,L)*ALPHA(IQ,I,J)*A(IQ,K,L,J)/3.0d0
c
c      IL+IR, backscattering, DCP_I
c      SDCP180=3.0d0*SAL1-SAL0
c
c      IL+/-IR, backscattering
       S180=7.0d0*SAL1+SAL0
       D180=8.0d0*(3.0d0*SAG1-SAG0+SA1)
       D=D180
       IF(S180.gt.1.0d-9)D=D/S180
c
c      IL+/-IR, 90o, x
       S90X=7.0d0*SAL1+SAL0
       D90X=2.0d0*(7.0d0*SAG1+SAG0+SA1)
       DX=D90X
       IF(S90X.gt.1.0d-9)DX=DX/S90X
c
c      IL+/-IR, 90o, z
       S90Z=2.0d0*(3.0d0*SAL1-SAL0)
       D90Z=4.0d0*(3.0d0*SAG1-SAG0-SA1)
       DZ=D90Z
       IF(S90Z.gt.1.0d-9)DZ=DZ/S90Z
c
c      polarized backscattering components (YDY+YDX=S180):
       YDY=3.0d0*SAL1-SAL0
       YDX=4.0d0*SAL1+2.0d0*SAL0
       P1=YDY/YDX
c
c      Try to get Gaussian units:
       AMU=1822.0d0
       BOHR=0.529177d0
       SpA3=SPAL/3.0d0*dsqrt((AMU*BOHR**4))
       beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4)
       YDY=6.0d0*YDY*(AMU*BOHR**4)
       YDX=6.0d0*YDX*(AMU*BOHR**4)
c      gpisvejc=980.0
       gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/EXCA
c
c
c      beta(A)^2=Delta2 in Gaussian:
       betaa=SA1*3.0d0/2.0d0*gpisvejc
       betag=(3.0d0*SAG1-SAG0)/2.0d0*gpisvejc
c
c      DCP 180
c      ram3=24.0d0*beta2
       roa3=8.0d0*(12.0d0*betag+4.0d0*betaa)
c
c      unknown:
       S90X=S90X*(AMU*BOHR**4)
       D90X=D90x*gpisvejc
       DX=D90X
       IF(S90X.gt.1.0d-9)DX=DX/S90X
c
c      ICP(90)=90Z
       roa2=8.0d0*(3.0d0*betag-betaa)
       ram2=12.0d0*beta2
       CID2=0.0d0
       if(ram2.gt.1.0d-9)CID2=roa2/ram2
c
c      ICP(180):
c      8/c*(12*beta(G')^2 + 4*beta(A)^2):
       roa1=96.0d0*(betag+betaa/3.0d0)
       ram1=4.0d0*(45.0d0*SpA3**2+7.0d0*beta2)
       CID1=0.0d0
       if(ram1.gt.1.0d-9)CID1=roa1/ram1

       if(ipr.eq.1)WRITE(4,3300)IQ,ECM,YDX,YDY,
     1 ram1,CID2,DX,CID1,P1,roa1,roa3
3300   FORMAT(I4,f9.2,6f9.2,f6.3,2f11.2)
       st(1,IQ)=ECM
       if(skind.eq.3)then
        st(2,IQ)=ram1
       else
        if(skind.eq.8)then
         st(2,IQ)=roa1
        else
         call report('unknown spectral type required')
        endif
       endif
c
      endif
      endif
      ENDIF
1     CONTINUE
      if(ipr.eq.1)then
       WRITE(4,3002)
       CLOSE(4)
      endif
      RETURN
      END

      SUBROUTINE TRTEN(N,S,ALPHA,A,G)
      IMPLICIT REAL*8 (A-H,O-Z)
      implicit integer*4 (i-n)
      DIMENSION S(N,N),A(N,3,3,3),ALPHA(N,3,3),G(N,3,3)
      real*8,allocatable::V(:)
      allocate(V(N))
      SFAC=0.0234280d0
      DO 1 I=1,3
      DO 1 J=1,3
      DO 3 JQ=1,N
3     V(JQ)=ALPHA(JQ,I,J)
      DO 1 IQ=1,N
      ALPHAS=0.0d0
      DO 2 IX=1,N
2     ALPHAS=ALPHAS+S(IX,IQ)*V(IX)*SFAC
1     ALPHA(IQ,I,J)=ALPHAS
      DO 11 I=1,3
      DO 11 J=1,3
      DO 31 JQ=1,N
31    V(JQ)=G(JQ,I,J)
      DO 11 IQ=1,N
      GS=0.0d0
      DO 21 IX=1,N
21    GS=GS+S(IX,IQ)*V(IX)*SFAC
11    G(IQ,I,J)=GS
      DO 12 I=1,3
      DO 12 J=1,3
      DO 12 K=1,3
      DO 32 JQ=1,N
32    V(JQ)=A(JQ,I,J,K)
      DO 12 IQ=1,N
      AS=0.0D0
      DO 22 IX=1,N
22    AS=AS+S(IX,IQ)*V(IX)*SFAC
12    A(IQ,I,J,K)=AS
      RETURN
      END

      subroutine new4(NAT,FCAR,ZIM,SMAT,Q,iwr,ty,r)
      IMPLICIT none
      INTEGER*4 I,J,IM,IA,NAT,N,iwr,ty(*),IERR
      REAL*8 FCAR(3*NAT,3*NAT),ZIM(*),SMAT(3*NAT,3*NAT),
     1Q(*),r(*),SIGN
      real*8,allocatable::ZM(:),E(:),FNEW(:,:)
      N=3*NAT
      allocate(ZM(N),E(N),FNEW(N,N))
      DO 220 I=1,NAT
      DO 220 J=1,3
220   ZM(3*I-3+J)=1.0d0/SQRT(ZIM(I))
      CALL FAST(N,FCAR,FNEW,ZM)
      CALL TRED12(N,FNEW,SMAT,Q,2,IERR,E)
      IF(IERR.NE.0)call report('diagonalization error')
C     Make imaginary frequencies negative and in cm-1:
      DO 310 I=1,N
      SIGN=1.0d0
      IF(Q(I).LT.0.0d0)SIGN=-1.0d0
310   Q(I)=SIGN*SQRT(ABS(Q(I)))*1302.828d0
C     Calculate the true mass-weighted S-matrix:
      DO 210 IA=1,N
      DO 210 IM=1,N
210   SMAT(IA,IM)=SMAT(IA,IM)*ZM(IA)
      if(iwr.eq.1)then
       CALL SMATOUT(SMAT,Q,N,ty,r)
       WRITE(*,*)'F.INP'
      endif
      return
      END

      SUBROUTINE SMATOUT(S,W,N,ty,r)
C     dX  =  S dQ   S(nat3xnint)
      IMPLICIT none
      INTEGER*4 I,J,N,ty(*),ix,IU,K
      real*8 S(N,N),W(*),amx,r(*)
c     Phase correction:
      DO 101 J=1,N
      amx=S(1,J)
      DO 102 I=2,N
102   if(dabs(S(I,J)).gt.dabs(amx))amx=S(I,J)
      if(amx.lt.0.0d0)then
       do 103 I=1,N
103    S(I,J)=-S(I,J)
      endif
101   continue
c
      OPEN(34,FILE='F.INP')
      WRITE(34,10)N,N,N/3
10    FORMAT(3I7)
      DO 3 I=1,N/3
3     WRITE(34,11)ty(I),(r(ix+3*(I-1)),ix=1,3)
11    FORMAT(I7,3F12.6)
      WRITE(34,14)
14    FORMAT(' Atom Mode    X-disp.    Y-disp.      Z-disp.')
      DO 1 I=1,N/3
      IU=3*(I-1)
      DO 1 J=1,N
      K=N-J+1
1     WRITE(34,15)I,K,S(IU+1,J),S(IU+2,J),S(IU+3,J)
15    FORMAT(2I7,3F11.6)
      WRITE(34,11)N
      WRITE(34,13)(W(I),I=1,N)
13    format(6F11.3)
      CLOSE(34)

      RETURN
      END

      SUBROUTINE FAST(NA,FCAR,FNEW,ZM)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4(I-N)
      DIMENSION FCAR(NA,NA),FNEW(NA,NA),ZM(*)
      DO 214 I=1,NA
      UI=ZM(I)
      DO 214 J=I,NA
      FNEW(I,J)=FCAR(I,J)*ZM(J)*UI
214   FNEW(J,I)=FNEW(I,J)
      RETURN
      END
C
      SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR,E)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(N,N),Z(N,N),D(*),E(*)
      IF (IEIG .LT. 0) GO TO 110
      DO 100 I = 1, N
C
      DO 100 J = 1, I
      Z(I,J) = A(I,J)
  100 CONTINUE
C
  110 ICODE = IABS(IEIG)
      IF (N .EQ. 1) GO TO 320
C     :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- ::::::::::
      DO 300 II = 2, N
      I = N + 2 - II
      L = I - 1
      H = 0.0D0
      SCALE = 0.0D0
      IF (L .LT. 2) GO TO 130
C     :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) ::::::::::
      DO 120 K = 1, L
  120 SCALE = SCALE + DABS(Z(I,K))
C
      IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140
  130 E(I) = Z(I,L)
      GO TO 290
C
  140 DO 150 K = 1, L
      F = Z(I,K) / SCALE
      Z(I,K) = F
      H = H + F * F
  150 CONTINUE
C
      G = -DSIGN(DSQRT(H),F)
      E(I) = SCALE * G
      H = H - F * G
      Z(I,L) = F - G
      F = 0.0D0
C
      DO 240 J = 1, L
      IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H)
      G = 0.0D0
C     :::::::::: FORM ELEMENT OF A*U ::::::::::
      DO 180 K = 1, J
  180 G = G + Z(J,K) * Z(I,K)
C
      JP1 = J + 1
      IF (L .LT. JP1) GO TO 220
C
      DO 200 K = JP1, L
  200 G = G + Z(K,J) * Z(I,K)
C     :::::::::: FORM ELEMENT OF P ::::::::::
  220 E(J) = G / H
      F = F + E(J) * Z(I,J)
  240 CONTINUE
C
      HH = F / (H + H)
C     :::::::::: FORM REDUCED A ::::::::::
      DO 260 J = 1, L
      F = Z(I,J)
      G = E(J) - HH * F
      E(J) = G
C
      DO 260 K = 1, J
      Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K)
  260 CONTINUE
C
      DO 280 K = 1, L
  280 Z(I,K) = SCALE * Z(I,K)
C
  290 D(I) = H
  300 CONTINUE
C
  320 D(1) = 0.0D0
      E(1) = 0.0D0
C     :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES ::::::::::
      IF (ICODE .NE. 2) GO TO 600
      DO 500 I = 1, N
      L = I - 1
      IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380
C
      DO 360 J = 1, L
      G = 0.0D0
C
      DO 340 K = 1, L
  340 G = G + Z(I,K) * Z(K,J)
C
      DO 360 K = 1, L
      Z(K,J) = Z(K,J) - G * Z(K,I)
  360 CONTINUE
C
  380 D(I) = Z(I,I)
      Z(I,I) = 1.0D0
      IF (L .LT. 1) GO TO 500
C
      DO 400 J = 1, L
      Z(I,J) = 0.0D0
      Z(J,I) = 0.0D0
  400 CONTINUE
C
  500 CONTINUE
C
  620 CALL TQL12 (N,Z,D,IERR,ICODE,E)
      RETURN
C     :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY ::::::::::
  600 DO 610 I=1,N
  610 D(I) = Z(I,I)
      GO TO 620
      END
C
      SUBROUTINE TQL12(N,Z,D,IERR,ICODE,E)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION Z(N,N),D(*),E(*)
      REAL*8 MACHEP
      EPS = 1.0D0
   10 EPS = 0.50D0*EPS
      TOL1 = EPS + 1.0D0
      IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10
      IF(TOL1-EPS.EQ.1.0D0)  EPS = EPS + EPS
      MACHEP = EPS
C
C     MACHEP=16.0D0**(-13)
C
      IERR = 0
      IF (N .EQ. 1) GO TO 1001
C
      DO 100 I = 2, N
  100 E(I-1) = E(I)
C
      F = 0.0D0
      B = 0.0D0
      E(N) = 0.0D0
C
      DO 240 L = 1, N
      J = 0
      H = MACHEP * (DABS(D(L)) + DABS(E(L)))
      IF (B .LT. H) B = H
C     :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT ::::::::::
      DO 110 M = L, N
      IF (DABS(E(M)) .LE. B) GO TO 120
C     :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
C                THROUGH THE BOTTOM OF THE LOOP ::::::::::
  110 CONTINUE
C
  120 IF (M .EQ. L) GO TO 220
  130 IF (J .EQ. 30) GO TO 1000
      J = J + 1
C     :::::::::: FORM SHIFT ::::::::::
      L1 = L + 1
      G = D(L)
      P = (D(L1) - G) / (2.0D0 * E(L))
      R = DSQRT(P*P+1.0D0)
      D(L) = E(L) / (P + DSIGN(R,P))
      H = G - D(L)
C
      DO 140 I = L1, N
  140 D(I) = D(I) - H
C
      F = F + H
C     :::::::::: QL TRANSFORMATION ::::::::::
      P = D(M)
      C = 1.0D0
      S = 0.0D0
      MML = M - L
C     :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- ::::::::::
      DO 200 II = 1, MML
      I = M - II
      G = C * E(I)
      H = C * P
      IF (DABS(P) .LT. DABS(E(I))) GO TO 150
      C = E(I) / P
      R = DSQRT(C*C+1.0D0)
      E(I+1) = S * P * R
      S = C / R
      C = 1.0D0 / R
      GO TO 160
  150 C = P / E(I)
      R = DSQRT(C*C+1.0D0)
      E(I+1) = S * E(I) * R
      S = 1.0D0 / R
      C = C * S
  160 P = C * D(I) - S * G
      D(I+1) = H + S * (C * G + S * D(I))
C     :::::::::: FORM VECTOR ::::::::::
      IF (ICODE .NE. 2) GO TO 200
      DO 180 K = 1, N
      H = Z(K,I+1)
      Z(K,I+1) = S * Z(K,I) + C * H
      Z(K,I) = C * Z(K,I) - S * H
  180 CONTINUE
C
  200 CONTINUE
C
      E(L) = S * P
      D(L) = C * P
      IF (DABS(E(L)) .GT. B) GO TO 130
  220 D(L) = D(L) + F
  240 CONTINUE
C     :::::::::: ORDER EIGENVALUES AND EIGENVECTORS ::::::::::
      DO 300 II = 2, N
      I = II - 1
      K = I
      P = D(I)
C
      DO 260 J = II, N
      IF (D(J) .LE. P) GO TO 260
      K = J
      P = D(J)
  260 CONTINUE
C
      IF (K .EQ. I) GO TO 300
      D(K) = D(I)
      D(I) = P
C
      IF (ICODE .NE. 2) GO TO 300
      DO 280 J = 1, N
      P = Z(J,I)
      Z(J,I) = Z(J,K)
      Z(J,K) = P
  280 CONTINUE
C
  300 CONTINUE
C
      GO TO 1001
C     :::::::::: SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30 ITERATIONS ::::::::::
 1000 IERR = L
 1001 RETURN
      END

      subroutine INPTRY(XM,N,kt,iso,riso,ISOT)
      IMPLICIT none
      real*8 XM(*),riso(*)
      integer*4 N,kt(*),iso(*),ISOT,MENDELEV,i
      PARAMETER (MENDELEV=89)
      real*8 amas(MENDELEV)
      data amas/1.007825,4.003,
     2  6.941, 9.012,  10.810,12.000,14.003074,15.994915,18.9984,20.179,
     3 22.990,24.305,   26.981,28.086,30.9737,31.972,34.968852,39.948,
     4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847,
     4 58.933,58.700,63.546,65.380,
     4                  69.720,72.590,74.922,78.960,79.904,83.800,
     5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070,
     5 102.906,106.400,107.868,112.410,
     5                 114.82,118.69,121.75,127.600,126.905,131.300,
     6 132.905,137.330,138.906,
     6                 140.120,140.908,144.240,145.000,150.400,
     6 151.960,157.250,158.925,162.500,164.930,167.260,168.934,
     6 173.040,174.970,
     6 178.490,180.948,183.850,186.207,190.200,192.220,195.090,
     6 196.967,207.590,204.370,207.200,208.980,210.000,210.001,
     6 222.02,
     7 223.000,226.025,227.028/     
      do 1 i=1,N
1     XM(i)=amas(kt(i))
      do 2 i=1,ISOT
2     XM(iso(i))=riso(i)
      return
      end

      subroutine tests(xv,MA,M,nx,loff,fi,NPOT,NAT,list)
      implicit none
      integer*4 MA,M,nats,i,ia,ix,j,ja,jx,NPOT,NAT,ig,inow,ip,list(*),
     1is,ie,nx,kx,lx,ii,jj,jxstart
      real*8 xv(MA,M),pi,fi(NPOT,nx)
      character*80 s80
      real*8,allocatable::FCAR(:,:),r(:),p(:,:,:),a(:,:,:),rr(:,:),
     1al(:,:,:),aa(:,:,:,:),g(:,:,:)
      integer*4,allocatable::ty(:)
      logical loff
      write(6,*)
      write(6,*)'Test entered'
      open(9,file='X.LST')
      open(10,file='FC.LST')
      open(11,file='TEN.LST')
      open(12,file='TTT.LST')
      do 1 ig=1,nx
      read(9,90)s80
 90   format(a80)
      call rdx(s80,nats,ty,r,0)
      allocate(r(3*nats),ty(nats))
      call rdx(s80,nats,ty,r,1)

      allocate(FCAR(3*nats,3*nats))
      read(10,90)s80
      call iss(s80,is,ie)
      open(20,file=s80(is:ie)//'.i')
      call READFF(3*nats,FCAR)
      close(20)
      open(20,file=s80(is:ie)//'.i.test')
      inow=0
      if(loff)then
       do 2 i=1,NAT
       ia=list(i)
       do 2 ix=1,3
       ii=3*(ia-1)+ix
       do 2 j=i,NAT
       ja=list(j)
       if(j.eq.i)then
        jxstart=ix
       else
        jxstart=1
       endif
       do 2 jx=jxstart,3
       jj=3*(ja-1)+jx
       inow=inow+1
       pi=xv(1,inow)
       do 21 ip=1,NPOT
21     pi=pi+xv(ip+1,inow)*fi(ip,ig)
2      write(20,200)i,ix,j,jx,pi-xv(1,inow),
     1 FCAR(ii,jj)-xv(1,inow),pi,FCAR(ii,jj)
200    format(4i5,4f12.6)
      else
       do 3 i=1,NAT
       ia=list(i)
       do 3 ix=1,3
       ii=3*(ia-1)+ix
       inow=inow+1
       pi=xv(1,inow)
       do 31 ip=1,NPOT
31     pi=pi+xv(ip+1,inow)*fi(ip,ig)
3      write(20,200)i,ix,i,ix,pi-xv(1,inow),
     1 FCAR(ii,ii)-xv(1,inow),pi,FCAR(ii,ii)
      endif
      close(20)
      deallocate(FCAR)

      allocate(p(nats,3,3),a(nats,3,3),rr(3,nats))
      do 4 i=1,nats
      do 4 ix=1,3
4     rr(ix,i)=r(ix+3*(i-1))      
      read(11,90)s80
      call iss(s80,is,ie)
      open(15,file=s80(is:ie)//'.i')
      CALL READTEN(nats,p,a,rr,0)
      CLOSE(15)
      open(20,file=s80(is:ie)//'.i.p.test')
      do 5 i=1,NAT
      ia=list(i)
      do 5 ix=1,3
      do 5 jx=1,3
      inow=inow+1
      pi=xv(1,inow)
      do 51 ip=1,NPOT
51    pi=pi+xv(ip+1,inow)*fi(ip,ig)
5     write(20,200)i,ix,i,ix,pi-xv(1,inow),p(ia,ix,jx)-xv(1,inow)
      close(20)
      open(20,file=s80(is:ie)//'.i.a.test')
      do 6 i=1,NAT
      ia=list(i)
      do 6 ix=1,3
      do 6 jx=1,3
      inow=inow+1
      pi=xv(1,inow)
      do 61 ip=1,NPOT
61    pi=pi+xv(ip+1,inow)*fi(ip,ig)
6     write(20,200)i,ia,ix,jx,pi-xv(1,inow),a(ia,ix,jx)-xv(1,inow)
      close(20)
      deallocate(p,a,rr)

      allocate(al(3*nats,3,3),g(3*nats,3,3),aa(3*nats,3,3,3))
      read(12,90)s80
      call iss(s80,is,ie)
      open(15,file=s80(is:ie)//'.i')
      CALL READTTT(3*nats,al,aa,g)
      CLOSE(15)
      open(20,file=s80(is:ie)//'.i.al.test')
      do 7 i=1,NAT
      ia=list(i)
      do 7 ix=1,3
      ii=3*(ia-1)+ix
      do 7 jx=1,3
      do 7 kx=1,3
      inow=inow+1
      pi=xv(1,inow)
      if(ig.eq.3.and.i.eq.1.and.ix.eq.1.and.jx.eq.1.and.kx.eq.1)then
       write(6,*)'inow ',inow
       write(6,*)'P0',xv(1,inow)
       write(6,*)'P1',xv(2,inow)
       write(6,*)'pot',fi(1,ig)
      endif
      do 71 ip=1,NPOT
71    pi=pi+xv(ip+1,inow)*fi(ip,ig)
7     write(20,200)i,ix,jx,kx,pi-xv(1,inow),al(ii,jx,kx)-xv(1,inow)
      close(20)

      open(20,file=s80(is:ie)//'.i.g.test')
      do 8 i=1,NAT
      ia=list(i)
      do 8 ix=1,3
      ii=3*(ia-1)+ix
      do 8 jx=1,3
      do 8 kx=1,3
      inow=inow+1
      pi=xv(1,inow)
      do 81 ip=1,NPOT
81    pi=pi+xv(ip+1,inow)*fi(ip,ig)
8     write(20,200)i,ix,jx,kx,pi-xv(1,inow),g(ii,jx,kx)-xv(1,inow)
      close(20)
      open(20,file=s80(is:ie)//'.i.a.test')
      do 9 i=1,NAT
      ia=list(i)
      do 9 ix=1,3
      ia=list(i)
      do 9 jx=1,3
      do 9 kx=1,3
      do 9 lx=1,3
      inow=inow+1
      pi=xv(1,inow)
      do 91 ip=1,NPOT
91    pi=pi+xv(ip+1,inow)*fi(ip,ig)
9     write(20,202)i,ia,ix,jx,kx,lx,pi-xv(1,inow),
     1aa(ii,jx,kx,lx)-xv(1,inow)
202   format(6i5,2f12.6)
      close(20)
      deallocate(al,g,aa)

      deallocate(r,ty)
1     continue
      close(9)
      close(10)
      close(11)
      close(12)
      return
      end

      SUBROUTINE READTTT(N,ALPHA,A,G)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3)
      NAT=N/3
      READ(15,*)
      READ(15,*)NATC
      if(NAT.NE.NATC)call report('Wrong number of atoms in .TTT file')
      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

      subroutine trafottt(s80,nats,list,tm,NAT,r,
     1bv,inow,M,MA,fi,NPOT,nx,ig)
      implicit none
      character*(*) s80
      real*8,allocatable::alpha(:,:,:),a(:,:,:,:),g(:,:,:),rr(:,:)
      integer nats,i,ia,jx,ix,list(*),nat,ip,kx,lx,ii
      real*8 tm(NAT,3,3),r(*)
      integer*4 M,MA,inow,NPOT,nx,ig
      real*8 bv(MA,M),fi(NPOT,nx)
      allocate(alpha(3*nats,3,3),g(3*nats,3,3),a(3*nats,3,3,3),
     1rr(3,nats))

      do 1 i=1,nats
      do 1 ix=1,3
1     rr(ix,i)=r(ix+3*(i-1))

      OPEN(15,FILE=s80,STATUS='OLD')
      CALL READTTT(3*nats,ALPHA,A,G)
      CLOSE(15)
      CALL TTTT(3*nats,ALPHA,A,G,1,rr)
      call trafott1(alpha,nat,nats,list,tm,0)
      call trafott1(g    ,nat,nats,list,tm,0)
      call trafott2(a,nat,nats,list,tm,0)

      do 2 i=1,NAT
      ia=list(i)
      do 2 ix=1,3
      ii=3*(ia-1)+ix
      do 2 jx=1,3
      do 2 kx=1,3
      inow=inow+1
      bv(1   ,inow)=bv(1   ,inow)+alpha(ii,jx,kx)
      do 2 ip=1,NPOT
2     bv(1+ip,inow)=bv(1+ip,inow)+alpha(ii,jx,kx)*fi(ip,ig)

      do 3 i=1,NAT
      ia=list(i)
      do 3 ix=1,3
      ii=3*(ia-1)+ix
      do 3 jx=1,3
      do 3 kx=1,3
      inow=inow+1
      bv(1   ,inow)=bv(1   ,inow)+g(ii,jx,kx)
      do 3 ip=1,NPOT
3     bv(1+ip,inow)=bv(1+ip,inow)+g(ii,jx,kx)*fi(ip,ig)

      do 4 i=1,NAT
      ia=list(i)
      do 4 ix=1,3
      ii=3*(ia-1)+ix
      do 4 jx=1,3
      do 4 kx=1,3
      do 4 lx=1,3
      inow=inow+1
      bv(1   ,inow)=bv(1   ,inow)+a(ii,jx,kx,lx)
      do 4 ip=1,NPOT
4     bv(1+ip,inow)=bv(1+ip,inow)+a(ii,jx,kx,lx)*fi(ip,ig)

      CALL WRITETTT(3*nats,ALPHA,A,G,s80//'.i')
      return
      end
c     ============================================================ 
      SUBROUTINE TTTT(N,ALPHA,A,G,ICO,R)
c     Tensors derivatives
C     Transforms A and G to local origins (ICO=1) or 
c     common origin (ICO=2);
C     sets A and G to zero for ICO=3.
C     R supposed to be passed in in A
C
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),R(3,N/3)
      real*8,allocatable::X(:,:)
      NAT=N/3
      allocate(X(3,NAT))
C
C     transfer into atomic coordinates:
      DO 5 I=1,3
      DO 5 J=1,NAT
5     X(I,J)=R(I,J)/0.529177D0
C
      IF(ICO.EQ.1)SIGN=1.0d0
      IF(ICO.EQ.2)SIGN=-1.0d0
      IF(ICO.EQ.3)SIGN=0.0d0
      W=1.0d0
C     supposedly the static limit G/w is calculated
C
      DO 2 IA=1,3
      DO 2 IB=1,3
      DO 2 L=1,NAT
      DO 2 IE=1,3
      IIND=3*(L-1)+IE
      SUM=0.0d0
      DO 1 IG=1,3
      DO 1 ID=1,3
1     SUM=SUM+SIGN*(0.5d0*W)*EPS(IB,IG,ID)*X(IG,L)*ALPHA(IIND,IA,ID)
      G(IIND,IA,IB)=G(IIND,IA,IB)+SUM
      IF(ICO.EQ.3)G(IIND,IA,IB)=0.0d0
2     CONTINUE
      IF(ICO.EQ.1)WRITE(6,*)' G transformed into atomic origins'
      IF(ICO.EQ.2)WRITE(6,*)' G transformed into common origin'
      IF(ICO.EQ.3)WRITE(6,*)' G set to zero'
C
c     A IA index is electric
      DO 3 IA=1,3
      DO 3 IB=1,3
      DO 3 IC=1,3
      DO 3 L=1,NAT
      DO 3 IE=1,3
      IIND=3*(L-1)+IE
      SUM=0.0d0
      IF(IB.EQ.IC)THEN
       DO 4 ID=1,3
4      SUM=SUM+X(ID,L)*ALPHA(IIND,IA,ID)*SIGN
      ENDIF
      A(IIND,IA,IB,IC)=A(IIND,IA,IB,IC)+SUM
     1-1.5d0*SIGN*(X(IB,L)*ALPHA(IIND,IA,IC)+X(IC,L)*ALPHA(IIND,IA,IB))
      IF(ICO.EQ.3)A(IIND,IA,IB,IC)=0.0d0
3     CONTINUE
      IF(ICO.EQ.1)WRITE(6,*)' A transformed into atomic origins'
      IF(ICO.EQ.2)WRITE(6,*)' A transformed into common origin'
      IF(ICO.EQ.3)WRITE(6,*)' A set to zero'
      RETURN
      END
      
      subroutine trafoten(s80,nats,list,tm,NAT,r,
     1bv,inow,M,MA,fi,NPOT,nx,ig)
      implicit none
      character*(*) s80
      real*8,allocatable::p(:,:,:),a(:,:,:),rr(:,:)
      integer nats,i,ia,jx,ix,list(*),nat,ip
      real*8 tm(NAT,3,3),r(*)
      integer*4 M,MA,inow,NPOT,nx,ig
      real*8 bv(MA,M),fi(NPOT,nx)
      allocate(p(nats,3,3),a(nats,3,3),rr(3,nats))
      do 1 i=1,nats
      do 1 ix=1,3
1     rr(ix,i)=r(ix+3*(i-1))
      open(15,file=s80)
      CALL READTEN(nats,p,a,rr,1)
      CLOSE(15)
      call trafote1(p,nat,nats,list,tm,0)
      call trafote1(a,nat,nats,list,tm,0)

      do 2 i=1,NAT
      ia=list(i)
      do 2 ix=1,3
      do 2 jx=1,3
      inow=inow+1
      bv(1   ,inow)=bv(1   ,inow)+p(ia,ix,jx)
      do 2 ip=1,NPOT
2     bv(1+ip,inow)=bv(1+ip,inow)+p(ia,ix,jx)*fi(ip,ig)

      do 3 i=1,NAT
      ia=list(i)
      do 3 ix=1,3
      do 3 jx=1,3
      inow=inow+1
      bv(1   ,inow)=bv(1   ,inow)+a(ia,ix,jx)
      do 3 ip=1,NPOT
3     bv(1+ip,inow)=bv(1+ip,inow)+a(ia,ix,jx)*fi(ip,ig)

      
      open(15,file=s80//'.i')
      CALL WRITETEN(nats,p,a,rr,NATS,0)
      CLOSE(15)

      return
      end

      subroutine trafott1(g,NAT,nats,list,tm,ic)
c     ic = 0 ... to local
c          1 ... to common
      implicit none
      integer nats,i,ia,jx,ix,list(*),nat,kx,ii,ic
      real*8 tm(NAT,3,3),g(3*nats,3,3),m(3,3),t(3,3,3),a(3,3,3)

      do 1 i=1,NAT
      ia=list(i)

      if(ic.eq.0)then
       do 101 ix=1,3
       do 101 jx=1,3
101    m(ix,jx)=tm(i,ix,jx)
      else
       do 100 ix=1,3
       do 100 jx=1,3
100    m(jx,ix)=tm(i,ix,jx)
      endif

      do 200 ix=1,3
      ii=ix+3*(ia-1)
      do 200 jx=1,3
      do 200 kx=1,3
200   a(ix,jx,kx)=g(ii,jx,kx)

      do 22 ix=1,3
      do 22 jx=1,3
      do 22 kx=1,3
22    t(ix,jx,kx)=
     1a(1,jx,kx)*m(ix,1)+a(2,jx,kx)*m(ix,2)+a(3,jx,kx)*m(ix,3)

      do 23 ix=1,3
      do 23 jx=1,3
      do 23 kx=1,3
23    a(ix,jx,kx)=
     1t(ix,1,kx)*m(jx,1)+t(ix,2,kx)*m(jx,2)+t(ix,3,kx)*m(jx,3)

      do 1 ix=1,3
      ii=ix+3*(ia-1)
      do 1 jx=1,3
      do 1 kx=1,3
1     g(ii,jx,kx)=
     1a(ix,jx,1)*m(kx,1)+a(ix,jx,2)*m(kx,2)+a(ix,jx,3)*m(kx,3)

      return
      end

      subroutine trafott2(a,NAT,nats,list,tm,ic)
c     ic = 0  ... to local
c          1  ... to common origin
      implicit none
      integer nats,i,ia,jx,ix,ixx,jxx,list(*),nat,kx,kxx,lx,lxx,ii,ic
      real*8 tm(NAT,3,3),a(3*nats,3,3,3),t(3,3,3,3),m(3,3)
      do 21 i=1,NAT
      ia=list(i)

      if(ic.eq.0)then
       do 101 ix=1,3
       do 101 jx=1,3
101    m(ix,jx)=tm(i,ix,jx)
      else
       do 100 ix=1,3
       do 100 jx=1,3
100    m(jx,ix)=tm(i,ix,jx)
      endif

      do 22 ix=1,3
      do 22 jx=1,3
      do 22 kx=1,3
      do 22 lx=1,3
      t(ix,jx,kx,lx)=0.0d0
      do 22 ixx=1,3
22    t(ix,jx,kx,lx)=t(ix,jx,kx,lx)+a(ixx+3*(ia-1),jx,kx,lx)*m(ix,ixx)

      do 23 ix=1,3
      ii=ix+3*(ia-1)
      do 23 jx=1,3
      do 23 kx=1,3
      do 23 lx=1,3
      a(ii,jx,kx,lx)=0.0d0
      do 23 jxx=1,3
23    a(ii,jx,kx,lx)=a(ii,jx,kx,lx)+t(ix,jxx,kx,lx)*m(jx,jxx)

      do 24 ix=1,3
      ii=ix+3*(ia-1)
      do 24 jx=1,3
      do 24 kx=1,3
      do 24 lx=1,3
      t(ix,jx,kx,lx)=0.0d0
      do 24 kxx=1,3
24    t(ix,jx,kx,lx)=t(ix,jx,kx,lx)+a(ii,jx,kxx,lx)*m(kx,kxx)

      do 21 ix=1,3
      ii=ix+3*(ia-1)
      do 21 jx=1,3
      do 21 kx=1,3
      do 21 lx=1,3
      a(ii,jx,kx,lx)=0.0d0
      do 21 lxx=1,3
21    a(ii,jx,kx,lx)=a(ii,jx,kx,lx)+t(ix,jx,kx,lxx)*m(lx,lxx)

      return
      end
      
      subroutine trafote1(p,NAT,nats,list,tm,ic)
c     ic = 0 ... to local
c          1     to common origin
      implicit none
      integer nats,i,ia,jx,ix,list(*),nat,ic
      real*8 tm(NAT,3,3),p(nats,3,3),t(3,3),m(3,3)
      do 21 i=1,NAT
      ia=list(i)

      if(ic.eq.0)then
       do 101 ix=1,3
       do 101 jx=1,3
101    m(ix,jx)=tm(i,ix,jx)
      else
       do 100 ix=1,3
       do 100 jx=1,3
100    m(jx,ix)=tm(i,ix,jx)
      endif

      do 22 ix=1,3
      do 22 jx=1,3
22    t(ix,jx)=p(ia,1,jx)*m(ix,1)+p(ia,2,jx)*m(ix,2)+p(ia,3,jx)*m(ix,3)

      do 21 ix=1,3
      do 21 jx=1,3
21    p(ia,ix,jx)=t(ix,1)*m(jx,1)+t(ix,2)*m(jx,2)+t(ix,3)*m(jx,3)

      return
      end

      SUBROUTINE READTEN(N0,P,A,X,ic)
c     ic=1 transform to internal
c        0 just read
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION P(N0,3,3),A(N0,3,3),X(3,N0)
      real*8,allocatable::AI(:,:,:),AJ(:,:,:),R(:,:),PV(:,:,:)
      allocate(AI(N0,3,3),AJ(N0,3,3),R(3,N0),PV(N0,3,3))
      READ (15,*) NOAT
      if(NOAT.ne.N0)then
       write(6,*)N0,NOAT
       call report('NAT <> NOAT')
      endif
C
C     Read dipole derivatives:
      DO 10 L=1,NOAT
      DO 10 J=1,3                                                     
10    READ (15,*) (P(L,J,I),I=1,3)
C
      DO 220 L=1,NOAT                                                 
      DO 220 J=1,3                                                    
220   READ (15,*) (AI(L,J,I),I=1,3)
      DO 230 L=1,NOAT                                                 
      DO 230 J=1,3                                                    
230   READ (15,*) (AJ(L,J,I),I=1,3)                                   
      DO 100 L=1,NOAT                                                 
      DO 100 J=1,3                                                    
100   READ (15,*) (PV(L,J,I),I=1,3)                                   
C     Distributed-origin gauge:
c     CALL VCDD0(AI,PV,P,R,N0)   
C     
C     Total axial tensor = electronic + nuclear:
      DO 4 L=1,N0
      DO 4 I=1,3
      DO 4 J=1,3
4     A(L,I,J)=AI(L,I,J)+AJ(L,I,J)
C
      if(ic.eq.1)then
       BOHR=0.52917705993d0
       DO 3 L=1,N0
       DO 3 I=1,3
3      R(I,L)=X(I,L)/BOHR
C      Put axial tensor in the origin on the atom L:
       DO 2201 L=1,N0
       DO 2201 J=1,3
       DO 2201 I=1,3
       SUM=0.0d0
       DO 2202 ID=1,3
       DO 2202 IG=1,3
2202   SUM=SUM-0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID)
2201   A(L,J,I)=A(L,J,I)+SUM
       WRITE(6,*)' AAT shifted to atomic origin'
      endif

      RETURN
      END

      SUBROUTINE WRITETEN(N0,P,A,X,NAT,ic)
c     ic=1   transform back to laboratory system
c        0   just write
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION P(N0,3,3),A(N0,3,3),X(3,NAT)
      real*8,allocatable::R(:,:)
      allocate(R(3,NAT))
      if(ic.eq.1)then
       BOHR=0.52917705993d0
       DO 3 L=1,NAT
       DO 3 I=1,3
3      R(I,L)=X(I,L)/BOHR
       DO 2201 L=1,NAT
       DO 2201 J=1,3
       DO 2201 I=1,3
       SUM=0.0d0
       DO 2202 ID=1,3
       DO 2202 IG=1,3
2202   SUM=SUM+0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID)
2201   A(L,J,I)=A(L,J,I)+SUM
       WRITE(6,*)' AAT shifted back to laboratory system'
      endif
      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)
      DO 220 L=1,NAT                                                 
      DO 220 J=1,3                                                    
220   WRITE(15,1501) (A(L,J,I),I=1,3),L
      Z0=0.0d0
      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(6,*)' .TEN.i written ...'
      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 vp(A,X,Y)
      implicit none
      real*8 A(*),X(*),Y(*)
      A(1)=X(2)*Y(3)-X(3)*Y(2)
      A(2)=X(3)*Y(1)-X(1)*Y(3)
      A(3)=X(1)*Y(2)-X(2)*Y(1)
      return
      end

      subroutine assxx(x,y,z,r,ia)
      implicit none
      real*8 x,y,z,r(*)
      integer*4 ia
      x=r(1+3*(ia-1))
      y=r(2+3*(ia-1))
      z=r(3+3*(ia-1))
      return
      end

      SUBROUTINE INV(nr,a,ai,n,e,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      dimension ai(nr,nr),a(nr,nr),e(nr,2*nr)
C
      TOL=1.0d-9
10000 IERR=0
      DO 1  ii=1,n
      DO 1  jj=1,n
      e(ii,jj)=a(ii,jj)
1     e(ii,jj+n)=0.0D0
      do 13 ii=1,n
13    e(ii,ii+n)=1.0D0
c
      DO 2 ii=1,n-1
      if (ABS(e(ii,ii)).LE.TOL) then
       DO 3 io=ii+1,n
3      if (ABS(e(io,ii)).GT.TOL) goto 11
       IERR=1
       write(6,*)ii
       write(6,*)'tol = ',tol
       tol=tol*0.50d0
       if(tol.gt.1.0d-20)goto 10000
       RETURN
c
11     CONTINUE
       DO 4 kk=1,2*n
       w=e(ii,kk)
       e(ii,kk)=e(io,kk)
4      e(io,kk)=w
      ENDIF
      eii=e(ii,ii)
      DO 5 jj=ii+1,n
      e1=e(jj,ii)/eii
      DO 6 kk=ii+1, 2*n
6     e(jj,kk)=e(jj,kk)-e(ii,kk)*e1
5     e(jj,ii)=0.0D0
2     CONTINUE
c
      DO 7  i2=n,2,-1
      eii=e(i2,i2)
      DO 7 j2=i2-1,1,-1
      e1=e(j2,i2)/eii
      DO 9 kk=1, n
9     e(j2,kk+n)=e(j2,kk+n)-e(i2,kk+n)*e1
7     e(j2,i2)=0.0d0
c
      DO 10 ii=1,n
      ei=1.0d0/e(ii,ii)
      DO 12  jj=1,n
12    ai(ii,jj)=e(ii,jj+n)*ei
10    CONTINUE
c
      RETURN
      END

      subroutine wmtr(A,n0,n,st,ic)
c
c     ic=0 .. writes triangle of a supposedly symmetric matrix
c     ic=1 .. writes it all
c
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      dimension A(n0,n)
      character*(*) st
      write(6,*)st
      write(6,*)
      N1=1
1     N3=MIN(N1+4,N)
      WRITE(6,18)(I,I=N1,N3)
18    FORMAT(4X,5I14)
      lnst=n1
      if(ic.eq.1)lnst=1
      DO 130 LN=lnst,N
      IEND=MIN(LN,N3)
      if(ic.eq.1)iend=n3
130   WRITE(6,17)LN,(A(LN,J),J=N1,iend)
      N1=N1+5
      IF(N3.LT.N)GOTO 1
17    FORMAT(I4,5D14.6)
      write(6,*)
      return
      end

      subroutine tf(s,nx)
      integer*4 nx
      character*(*) s
      nx=0
      open(9,file=s,status='old')
3     read(9,*,end=79,err=79)
      nx=nx+1
      goto 3
79    close(9)
      write(6,*)nx,' lines in '//s
      return
      end

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

      SUBROUTINE READFF(N,FCAR)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(N,N)
      N1=1
1     N3=N1+4
      IF(N3.GT.N)N3=N
      DO 130 LN=N1,N
130   READ(20,17)(FCAR(LN,J),J=N1,MIN(LN,N3))
      N1=N1+5
      IF(N3.LT.N)GOTO 1
17    FORMAT(4X,5D14.6)
      DO 31 I=1,N
      DO 31 J=I+1,N
   31 FCAR(I,J)=FCAR(J,I)
      RETURN
      END

      SUBROUTINE WRITEFF(N,FCAR)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(N,N)
      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 tz(t,m,n)
      implicit none
      real*8 t(m,n)
      integer*4 m,n,i,j
      do 1 i=1,m
      do 1 j=1,n
1     t(i,j)=0.0d0
      return
      end

      subroutine trff(s80,nats,tm,list,NAT,bv,inow,loff,M,MA,fi,
     1NPOT,nx,ig)
      implicit none
      character*(*) s80
      real*8,allocatable::t(:,:),FCAR(:,:)
      integer nats,list(*),inow,NAT,i,ia,j,ja,ix,jx,ixx,jxx,M,MA,
     1NPOT,nx,ig,ip,ii,jj,jxstart
      real*8 tm(NAT,3,3),bv(MA,M),fi(NPOT,nx)
      logical loff
      allocate(t(3*nats,3*nats),FCAR(3*nats,3*nats))
      open(20,file=s80)
      CALL READFF(3*nats,FCAR)
      close(20)
      do 19 i=1,NAT
      ia=list(i)
      do 19 ix=1,3
      ii=ix+3*(ia-1)
      do 19 j=1,NAT
      ja=list(j)
      do 19 jx=1,3
      jj=jx+3*(ja-1)
      t(ii,jj)=0.0d0
      do 19 ixx=1,3
19    t(ii,jj)=t(ii,jj)+FCAR(ixx+3*(ia-1),jj)*tm(i,ix,ixx)

      do 20 i=1,NAT
      ia=list(i)
      do 20 ix=1,3
      ii=ix+3*(ia-1)
      do 20 j=1,NAT
      ja=list(j)
      do 20 jx=1,3
      jj=jx+3*(ja-1)
      FCAR(ii,jj)=0.0d0
      do 20 jxx=1,3
20    FCAR(ii,jj)=FCAR(ii,jj)+t(ii,jxx+3*(ja-1))*tm(j,jx,jxx)

      open(20,file=s80//'.i')
      CALL WRITEFF(3*nats,FCAR)
      close(20)
      if(loff)then
       do 1 i=1,NAT
       ia=list(i)
       do 1 ix=1,3
       ii=ix+3*(ia-1)
       do 1 j=i,NAT
       ja=list(j)
       if(j.eq.i)then
        jxstart=ix
       else
        jxstart=1
       endif
       do 1 jx=jxstart,3
       jj=jx+3*(ja-1)
       inow=inow+1
       bv(1,inow)=bv(1,inow)+FCAR(ii,jj)
       do 1 ip=1,NPOT
1      bv(1+ip,inow)=bv(1+ip,inow)+FCAR(ii,jj)*fi(ip,ig)
      else
       do 2 i=1,NAT
       ia=list(i)
       do 2 ix=1,3
       inow=inow+1
       bv(1,inow)=bv(1,inow)+FCAR(ix+3*(ia-1),ix+3*(ia-1))
       do 2 ip=1,NPOT
2      bv(1+ip,inow)=bv(1+ip,inow)+FCAR(ix+3*(ia-1),ix+3*(ia-1))
     1 *fi(ip,ig)
      endif
      return
      end

      subroutine trfftoCOM(FCAR,nats,tm,list,NAT)
      implicit none
      real*8,allocatable::t(:,:)
      integer nats,list(*),NAT,i,ia,j,ja,ix,jx,ixx,jxx,ii,jj
      real*8 tm(NAT,3,3),FCAR(3*nats,3*nats)
      allocate(t(3*nats,3*nats))
      do 19 i=1,NAT
      ia=list(i)
      do 19 ix=1,3
      ii=ix+3*(ia-1)
      do 19 j=1,NAT
      ja=list(j)
      do 19 jx=1,3
      jj=jx+3*(ja-1)
      t(ii,jj)=0.0d0
      do 19 ixx=1,3
19    t(ii,jj)=t(ii,jj)+FCAR(ixx+3*(ia-1),jj)*tm(i,ixx,ix)
      do 20 i=1,NAT
      ia=list(i)
      do 20 ix=1,3
      ii=ix+3*(ia-1)
      do 20 j=1,NAT
      ja=list(j)
      do 20 jx=1,3
      jj=jx+3*(ja-1)
      FCAR(ii,jj)=0.0d0
      do 20 jxx=1,3
20    FCAR(ii,jj)=FCAR(ii,jj)+t(ii,jxx+3*(ja-1))*tm(j,jxx,jx)
      return
      end

      subroutine setpot(NPOT,plist,r,nats,ty,fi,nx,ig,list,NAT)
      implicit none
      integer*4 NPOT,plist(*),j,ia,ja,jj,nats,list(*),ty(*),nx,ig,
     1NAT
      real*8 xa,ya,za,r(*),x,y,z,q,d,fi(NPOT,nx)
      do 8 j=1,NPOT
      ia=plist(j)
      call assxx(xa,ya,za,r,ia)
      do 8 ja=1,nats
      do 9 jj=1,NAT
9     if(ja.eq.list(jj))goto 8
      call assxx(x,y,z,r,ja)
      if(ty(ja).eq.1)then
       q=0.4d0
      else
       if(ty(ja).eq.8)then
        q=-0.8d0
       else
        write(6,*)ja,ty(ja)
        call report('unknown solvent atom')
       endif
      endif
      d=dsqrt((x-xa)**2+(y-ya)**2+(z-za)**2)
      fi(j,ig)=fi(j,ig)+q/d
8     continue
      return
      end

      subroutine defcoord(icoord,NAT,r,list)
      implicit none
      integer*4 NAT,icoord(NAT,2),list(*),i,ia,i1,ii,iia,i2
      real*8 x,y,z,r(*),d2min,d2,xa,ya,za
      open(9,file='COORD.PRM')
      write(9,600)
600   format('   Local coordinate systems:',/,
     1       '   list   atom     i1     i2',/,
     1       ' ---------------------------')
      do 12 i=1,NAT
      ia=list(i)
      call assxx(x,y,z,r,ia)
      i1=1
      d2min=1.0d9
      do 121 ii=1,NAT
      iia=list(ii)
      call assxx(xa,ya,za,r,iia)
      d2=(x-xa)**2+(y-ya)**2+(z-za)**2
      if(d2.lt.d2min.and.iia.ne.ia)then
       i1=iia
       d2min=d2
      endif
121   continue
      i2=1
      d2min=1.0d9
      do 122 ii=1,NAT
      iia=list(ii)
      call assxx(xa,ya,za,r,iia)
      d2=(x-xa)**2+(y-ya)**2+(z-za)**2
      if(d2.lt.d2min.and.iia.ne.ia.and.iia.ne.i1)then
       i2=iia
       d2min=d2
      endif
122   continue
      write(9,6001)i,ia,i1,i2
6001  format(4i7)
      icoord(i,1)=i1
12    icoord(i,2)=i2
      close(9)
      write(6,*)'Local coordinates to COORD.PRM'
      return
      end

      subroutine deftm(NAT,list,icoord,r,tm)
      implicit none
      integer*4 NAT,list(*),icoord(NAT,2),i,ia,i1,i2,ix
      real*8 r(*),a(3),b(3),c(3),tm(NAT,3,3)
      do 18 i=1,NAT
      ia=list(i)
      i1=icoord(i,1)
      i2=icoord(i,2)
      do 14 ix=1,3
      a(ix)=r(ix+3*(i1-1))-r(ix+3*(ia-1))
14    b(ix)=r(ix+3*(i2-1))-r(ix+3*(ia-1))
      call vp(c,a,b)
      call vp(b,c,a)
      call vnorm(a)
      call vnorm(b)
      call vnorm(c)
      do 17 ix=1,3
      tm(i,1,ix)=a(ix)
      tm(i,2,ix)=b(ix)
17    tm(i,3,ix)=c(ix)
18    continue
      return
      end

      subroutine vnorm(v)
      real*8 v(*),vn
      vn=dsqrt(v(1)**2+v(2)**2+v(3)**2)
      if(vn.ne.0.0d0)then
       v(1)=v(1)/vn
       v(2)=v(2)/vn
       v(3)=v(3)/vn
      endif
      return
      end

      subroutine rdx(s80,nats,ty,r,ic)
c     read geometry
c     ic=0 only nats
c        1 all
      implicit none
      character*(*) s80
      integer*4 nats,ty(*),ia,ix,ic
      real*8 r(*)
      open(91,file=s80)
      read(91,*)
      read(91,*)nats
      if(ic.ne.1)goto 99
      do 16 ia=1,nats
16    read(91,*)ty(ia),(r(ix+3*(ia-1)),ix=1,3)
99    close(91)
      return
      end

      SUBROUTINE WRITETTT(N,ALPHA,A,G,s)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3)
      character*(*)s
      NAT=N/3
      OPEN(2,FILE=s)
      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(IIND,I,J),J=1,3)
2001  FORMAT(I5,1H ,I1,3F18.7)
      WRITE(2,2004)
2004  FORMAT(' The electric dipole magnetic dipole polarizability:',/,
     1       ' Atom/x    jx           jy           jz')
      DO 2 I=1,3
      WRITE(2,2003)I
2003  FORMAT(' G(',I1,',J):')
      DO 2 L=1,NAT
      DO 2 IX=1,3
      IIND=3*(L-1)+IX
2     WRITE(2,2001)L,IX,(G(IIND,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
      IIND=3*(L-1)+IX
3     WRITE(2,2001)L,IX,(A(IIND,I,J,K),K=1,3)
      CLOSE(2)
      WRITE(6,*)s//' written '
      RETURN
      END

      subroutine iss(s80,is,ie)
      implicit none
      integer*4 is,ie,i
      character*(*) s80
      do 1 i=len(s80),1,-1
 1    if(s80(i:i).ne.' ')goto 2
 2    ie=i
      do 3 i=1,len(s80)
 3    if(s80(i:i).ne.' ')goto 4
 4    is=i
      return
      end