PROGRAM VSCALE
      IMPLICIT none
      integer*4 IRANGE,JRANGE,NIT,NOAT,NA,np,ni,ity,NNU,NOSC
      real*8 DEL,hw,wmin,wmax
      integer*4,allocatable::NM(:),li(:,:),ZQ(:)
      real*8,allocatable::FCAR(:,:),SCFAC(:),ZM(:),
     1ZNU(:),ZIM(:),CHR(:),C(:,:),b(:,:),fi(:),s(:),sx(:)
      character*80 spec
      logical lproject,LG,lconj
C
      WRITE(*,60000)
60000 FORMAT('                 FORCE FIELD REFINEMENT ',/,/,
     1       '  Run before: new1 a, define coordinates -> FILE.UMA',/,
     1       '              bumatd -> INTY.FC and B.MAT',/,
     1       '              goin, def scale factors -> FTRY.INP',/,/,
     1       '       Input: FILE.X      FILE.FC',/,
     1       '              INTY.FC     VSCALE.OPT',/,
     1       '              FILE.UMA    FTRY.INP',/)
C
      NOAT=1
      allocate(ZQ(1),C(3,1))
      call rdx(NOAT,0,ZQ,C)
      deallocate(ZQ,C)
      NA=3*NOAT
      allocate(FCAR(NA,NA),SCFAC(3*NOAT),ZM(NA),ZNU(NA),
     1NM(3*NOAT),ZIM(NOAT),CHR(NOAT),C(3,NOAT),ZQ(NOAT))
      call rdx(NOAT,1,ZQ,C)
C
C     Read-in cartesian FF from FILE.FC:
      CALL READFF(NA,FCAR)
C
C     Read-in scalefactors etc from FTRY.INP:
      call readinp(NOSC,SCFAC,NM,CHR,NOAT,NA,NNU,ZIM,ZM,ZNU)
C
      call readopt(lproject,IRANGE,JRANGE,DEL,NIT,spec,ity,hw,LG,
     1wmin,wmax,lconj)

c     read spectrum
      np=1
      allocate(sx(np),s(np))
      call rsp(0,np,sx,s,spec)
      deallocate(sx,s)
      allocate(sx(np),s(np))
      call rsp(1,np,sx,s,spec)

c     read FILE.UMA, ni .. number of internal coordinates
      ni=1
      allocate(li(3,ni))
      call ru(NOAT,0,ni,li)
      deallocate(li)
      allocate(li(3,ni))
      call ru(NOAT,1,ni,li)

c     read B-matrix, stretching only
      allocate(b(ni,6))
      call rb(ni,li,b)

c     read internal force constants,diagonal only
      allocate(fi(ni))
      call rFF(ni,fi)

      if(lconj)then
       CALL IMPROVEG(NIT,NOSC,DEL,FCAR,NA,ity,np,sx,s,
     1 SCFAC,NM,ZNU,ZQ,C,ni,li,fi,b,lproject,ZM,hw,LG,
     1 IRANGE,JRANGE,wmin,wmax)
      else
       CALL IMPROVE(NIT,NOSC,DEL,FCAR,NA,ity,np,sx,s,
     1 SCFAC,NM,ZNU,ZQ,C,ni,li,fi,b,lproject,ZM,hw,LG,
     1 IRANGE,JRANGE,wmin,wmax)
      endif

c     write new scale factors:
      call wrfi(ni,NOSC,SCFAC,NM,NOAT,ZIM,CHR)

      END
C
      SUBROUTINE IMPROVE(NIT,NOSC,DEL,FCAR,NA,ity,np,sx,s,
     1SCFAC,NM,ZNU,ZQ,C,ni,li,fi,b,lproject,ZM,hw,LG,i1,i2,wmin,wmax)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      integer*4 ZQ(*),NIT,NOSC,NA,ity,np,ni,li(3,ni),NM(*),i1,i2
      real*8 ZM(NA),FCAR(NA,NA),SCFAC(*),C(3,NA/3),ZNU(*),sx(np),
     1s(np),fi(ni),b(ni,6),hw,wmin,wmax
      logical lproject,LG
      IIT=0
 1    IIT=IIT+1
      ICHANGE=0
      ie=0
      DO 2 I=1,NOSC
      X0=SCFAC(I)
      S0=SUMk(NOSC,0,ie,NA,FCAR,ZNU,SCFAC,NM,0,ZQ,
     1C,ni,li,fi,b,lproject,ZM,ity,s,sx,hw,LG,np,i1,i2,wmin,wmax)
      SCFAC(I)=X0-DEL
      SM=SUMk(NOSC,0,ie,NA,FCAR,ZNU,SCFAC,NM,0,ZQ,
     1C,ni,li,fi,b,lproject,ZM,ity,s,sx,hw,LG,np,i1,i2,wmin,wmax)
      SCFAC(I)=X0+DEL
      SP=SUMk(NOSC,0,ie,NA,FCAR,ZNU,SCFAC,NM,0,ZQ,
     1C,ni,li,fi,b,lproject,ZM,ity,s,sx,hw,LG,np,i1,i2,wmin,wmax)
      SCFAC(I)=X0
      write(6,607)I,X0,X0-DEL,X0+DEL,S0,SM,SP
 607  format(i3,3f7.4,3e12.4)
      IF (SP.LT.S0)THEN
       S0=SP
       SCFAC(I)=X0+DEL
       ICHANGE=1
      ENDIF
      IF (SM.LT.S0)THEN
       S0=SM
       SCFAC(I)=X0-DEL
       ICHANGE=1
      ENDIF
2     CONTINUE
      IF (ICHANGE.EQ.0)DEL=DEL/2.0D0
      WRITE(*,6000)IIT,DEL,S0
6000  FORMAT(I4,' del = ',F15.9,',  error = ',F15.4)
      IF (DEL.GT.0.00001d0.and.IIT.le.NIT)GOTO 1
      SP=SUMk(NOSC,1,ie,NA,FCAR,ZNU,SCFAC,NM,1,ZQ,
     1C,ni,li,fi,b,lproject,ZM,ity,s,sx,hw,LG,np,i1,i2,wmin,wmax)
      RETURN
      END

      subroutine setg(m,ie,NA,FCAR,ZNU,c,NM,ZQ,g,gnorm,
     1X,ni,li,fi,b,lproject,ZM,ity,s,sx,hw,LG,np,i1,i2,wmin,wmax)
      implicit none
      integer*4 m,ie,NA,NM(*),ZQ(*),ni,li(3,ni),ity,np,i1,i2,i
      real*8 FCAR(NA,NA),ZNU(*),c(m),g(m),X(3,NA/3),fi(ni),
     1b(ni,6),ZM(*),s(*),sx(*),hw,wmin,wmax,dd,ep,em,SUMk,gnorm
      real*8,allocatable::sft(:)
      logical lproject,LG
      allocate(sft(m))
      dd=0.01d0
      sft=c
      do 1 i=1,m
      sft(i)=c(i)+dd
      ep=SUMk(m,0,ie,NA,FCAR,ZNU,sft,NM,0,ZQ,
     1X,ni,li,fi,b,lproject,ZM,ity,s,sx,hw,LG,np,i1,i2,wmin,wmax)
      sft(i)=c(i)-dd
      em=SUMk(m,0,ie,NA,FCAR,ZNU,sft,NM,0,ZQ,
     1X,ni,li,fi,b,lproject,ZM,ity,s,sx,hw,LG,np,i1,i2,wmin,wmax)
      sft(i)=c(i)
 1    g(i)=(ep-em)/dd/2.0d0
      gnorm=0.0d0
      do 2 i=1,m
  2   gnorm=gnorm+g(i)**2
      gnorm=dsqrt(gnorm)
      return
      end

      SUBROUTINE IMPROVEG(maxit,NOSC,DEL,FCAR,NA,ity,np,sx,s,
     1SCFAC,NM,ZNU,ZQ,C,ni,li,fi,b,lproject,ZM,hw,LG,i1,i2,wmin,wmax)
      IMPLICIT none
      integer*4 ZQ(*),maxit,NOSC,NA,ity,np,ni,li(3,ni),NM(*),i1,i2,
     1nstep,m,lstep,i,ie
      real*8 ZM(NA),FCAR(NA,NA),SCFAC(NOSC),C(3,NA/3),ZNU(*),sx(np),
     1s(np),fi(ni),b(ni,6),hw,wmin,wmax,maxstep,e0,e00,SUMk,DEL,
     1dh,secder,eold,h,hh0,etot,gnorm0,enew,e1,hh2,hh1,e2,gradh,
     1gnorm1,glim
      logical lproject,LG
      real*8,allocatable::g0(:),h0(:),g(:),c0(:),ck(:)
c     conjugate gradient
      maxstep=DEL/5.0d0
      glim=1.0d-5
      ie=0
      nstep=0
c     number of parameters to optimize:
      m=NOSC
      allocate(g0(m),h0(m),g(m),c0(NOSC),ck(m))
c     initial values 
      c0=SCFAC
c     "energy" to be minimized"
      etot=SUMk(m,0,ie,NA,FCAR,ZNU,c0,NM,0,ZQ,
     1C,ni,li,fi,b,lproject,ZM,ity,s,sx,hw,LG,np,i1,i2,wmin,wmax)
c     gradient:
      call setg(m,ie,NA,FCAR,ZNU,c0,NM,ZQ,g0,gnorm0,
     1C,ni,li,fi,b,lproject,ZM,ity,s,sx,hw,LG,np,i1,i2,wmin,wmax)
      e0=etot
      h0=g0

222   e00=etot
      dh=maxstep/gnorm0
      secder=-1.0d0
      nstep=nstep+1
      write(6,600)nstep,etot
 600  format(' Linear search at point ',i6,', E = ',e12.4)
      eold=etot
      e0=etot
      h=0.0d0
      hh0=0.0d0
      lstep=0

111   lstep=lstep+1
      h=h+dh
      do 1 i=1,m
1     ck(i)=c0(i)-h*h0(i)
      etot=SUMk(m,0,ie,NA,FCAR,ZNU,ck,NM,0,ZQ,
     1C,ni,li,fi,b,lproject,ZM,ity,s,sx,hw,LG,np,i1,i2,wmin,wmax)
      enew=etot
c     possibly end of linear search?:
      if(lstep.gt.100)then
       write(6,*)'Too many linear steps'
       goto 444
      endif
      if(dabs(enew-eold)<1.0d-5*eold)goto 444
      if(abs((enew-eold)/dh)<0.5*glim*gnorm0) goto 444
      if( lstep.eq.1) then
       e1=etot
       hh1=h
      else
       if(lstep.eq.2) then
        e2=etot
        hh2=h
       else 
        e0=e1
        hh0=hh1
        e1=e2
        hh1=hh2
        e2=etot
        hh2=h
       endif
       if(hh1.eq.hh0.or.hh2.eq.hh0.or.hh2.eq.hh1)then
        secder=-1.0d0
       else
        secder=-2.0d0*((hh2-hh1)*(e1-e0)-(e2-e1)*(hh1-hh0))
     1  /((hh2-hh0)*(hh1-hh0)*(hh2-hh1))
        gradh=(e0-e1-secder/2.0d0*(hh0**2-hh1**2))/(hh0-hh1)
       endif
      endif
      if(secder.gt.0.0d0)then
       dh=-gradh/secder-h
       if(dabs(dh).gt.maxstep/gnorm0)then
        if(dh.gt.0.0d0)then
         dh=maxstep/gnorm0
        else
         dh=-maxstep/gnorm0
        endif
       endif
       eold=enew
       goto 111
      endif
      if(enew.gt.eold)then
       eold=enew
       dh=-dh/2
       goto 111
      endif
      eold=enew
      goto 111

444   write(6,601)lstep,etot
 601  format(i6,' linear steps ',e12.4)
      etot=SUMk(m,0,ie,NA,FCAR,ZNU,ck,NM,0,ZQ,
     1C,ni,li,fi,b,lproject,ZM,ity,s,sx,hw,LG,np,i1,i2,wmin,wmax)
      call setg(m,ie,NA,FCAR,ZNU,ck,NM,ZQ,g,gnorm1,
     1C,ni,li,fi,b,lproject,ZM,ity,s,sx,hw,LG,np,i1,i2,wmin,wmax)
      c0=ck
      do 2 i=1,m
      h0(i)=g(i)+(gnorm1/gnorm0)**2*h0(i)
   2  g0(i)=g(i)
      write(6,602)gnorm1,etot,etot-e00
 602  format(' Grad:',e12.4,' E:',e12.4,' dE:',e12.4)
      if(dabs(e00-etot).lt.1.0d-5*e00)then
       write(6,*)'Optimization converged'
       goto 555
      endif
      if(nstep.gt.maxit)then
       write(6,*)'too many steps'
       goto 555
      endif
      e0=etot
      gnorm0=gnorm1
      goto 222
555   write(6,605)etot
 605  format(' Error:',e12.4)
      SCFAC=c0
      etot=SUMk(NOSC,1,ie,NA,FCAR,ZNU,SCFAC,NM,1,ZQ,
     1C,ni,li,fi,b,lproject,ZM,ity,s,sx,hw,LG,np,i1,i2,wmin,wmax)
      RETURN
      END

      FUNCTION SUMk(NOSC,iwr,iee,NA,FCAR,ZNU,SCFAC,NM,is,ZQ,
     1C,ni,li,fi,b,lproject,ZM,ity,s,sx,hw,LG,np,i1,i2,wmin,wmax)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      integer*4 ni,li(3,ni),ity,ZQ(*),NM(*),is,ir,np
      real*8 FCAR(NA,NA),SCFAC(*),hw,ZNU(*),C(3,NA/3),saup(3,3),
     1sadown(3,3),fi(ni),b(ni,6),ZM(*),s(*),sx(*),wmin,wmax
      logical lproject,LG
      real*8,allocatable::FNEW(:,:),F(:,:),SMAT(:,:),Q(:),AS(:,:),
     1se(:)
c
      allocate(SMAT(NA,NA),Q(NA),AS(3,NA))
      wtol=0.01d0
      c1=1302.828d0
      NOAT=NA/3

c     FCAR:original FF
c     F:scaled force field:
      allocate(F(NA,NA))
      F=FCAR
      do 1 ia=1,NOAT
      saup=0.0d0
      sadown=0.0d0
      ib=0
c     bonds containing atom ia:
      do 2 ii=1,ni
      if(li(1,ii).eq.1.and.(li(2,ii).eq.ia.or.li(3,ii).eq.ia))then
       ib=ib+1
c      scale factor for this bond
       if(NM(1).ge.ii)then
        ir=1
       endif
       if(NM(1).ge.ii)goto 200
       do 3 ir=2,NOSC
3      if(NM(ir-1)+1.le.ii.and.NM(ir).ge.ii)goto 200
       write(6,*)' atom ',ia,' bond ',ii
       call report('scale factor not found')
200    sfb=SCFAC(ir)
       if(ia.eq.li(2,ii))then
        io=0
       else
        io=3
       endif
       do 4 ix=1,3
       do 4 jx=1,3
       saup(  ix,jx)=saup(  ix,jx)+sfb*fi(ii)*b(ii,ix+io)*b(ii,jx+io)
4      sadown(ix,jx)=sadown(ix,jx)+    fi(ii)*b(ii,ix+io)*b(ii,jx+io)
      endif
2     continue
      if(ib.gt.0)then
       do 5 ix=1,3
       do 5 jx=1,3
5      F(ix+3*(ia-1),jx+3*(ia-1))=F(ix+3*(ia-1),jx+3*(ia-1))
     1 *saup(ix,jx)/sadown(ix,jx)
      endif
1     continue

c     project translations, rotations:
      IF(lproject)CALL PROJECT(F,NA,C)
c     FNEW: mass-weighted FF:
      allocate(FNEW(NA,NA))
      CALL FAST(NA,F,FNEW,ZM)

      CALL TRED12(NA,FNEW,SMAT,Q,2,IERR)
      IF(IERR.NE.0)call report('diagonalization not ok')
C     True mass-weighted S-matrix:
      DO 210 IA=1,NA
      DO 210 IM=1,NA
210   SMAT(IA,IM)=SMAT(IA,IM)*ZM(IA)
c      
      RMSW=0.0D0
      ie=0
      DO 310 I=1,NA
      T=Q(I)
      SIGN=1.0d0
C     Make imaginary frequencies negative:
      IF(T.LT.0.0d0)SIGN=-1.0d0
      T=SIGN*SQRT(ABS(T))*c1
      Q(I)=T
      TP=ZNU(I)
      DEL=T-TP
      if(ABS(TP).le.wtol)then
       DEL=0.0d0
      else
       ie=ie+1
      endif
      RMSW=RMSW+DEL**2
      AS(1,I)=T
      AS(2,I)=TP
310   AS(3,I)=T-TP
      iee=ie
      if(ie.gt.0)then
       RMSW=DSQRT(RMSW/DFLOAT(ie))
      else
       if(NA.gt.0)RMSW=DSQRT(RMSW/DFLOAT(NA))
      endif
      SO1=0.0D0
      SO2=0.0D0
      DO 5151 J=1,NA
      SO1=SO1+AS(1,J)*AS(2,J)
5151  SO2=SO2+AS(1,J)*AS(1,J)
      AK=SO1/SO2

      allocate(se(np))
      if(ity.eq.3.or.ity.eq.8)then
       call mksp(NA,ity,np,sx,se,AS,SMAT,C,hw,LG)
      else
       write(6,*)ity
       call report('unknown spectral type')
      endif
      an1=0.0d0
      an2=0.0d0
      do 504 i=1,np
      if(sx(i).ge.wmin.and.sx(i).le.wmax)then
       an1=an1+dabs(s(i))
       an2=an2+dabs(se(i))
      endif
504   continue
      ser=0.0d0
      do 505 i=1,np
505   if(sx(i).ge.wmin.and.sx(i).le.wmax)
     1ser=ser+dabs(s(i)/an1-se(i)/an2)

      call wrsp(np,sx,se,an1,an2)

      SUMk=ser
C
      if(iwr.gt.0)then
       WRITE(*,5555)
5555   FORMAT(1X,'************  NEW FREQUENCIES ******************')
       DO 311 J=NA,1,-1
       a3=0.0d0
       if(AS(2,NA-J+1).gt.wtol)a3=AS(3,NA-J+1)
311    WRITE(*,6666)J,NA-J+1,(AS(I,NA-J+1),I=1,2),a3
6666   FORMAT(1X,
     1 I5,' (',I4,') ...... ',F8.2,'   -',F8.2,'/',F8.2)
       WRITE(*,5551)
5551   FORMAT(1X,'************************************************')
       WRITE(6,6667)ie,RMSW,AK
6667   FORMAT(i6,'experimental frequencies found',/,
     1  ' >>>>><< >> ',F8.2,' <<<<<<  ',F8.6)
      endif

      if(is.gt.0)then
       CALL SMATOUT(SMAT,AS,NA,i1,i2,ZQ,C)
       OPEN(20,FILE='NEW.FC')
       CALL WRITEFF(NA,F)
       CLOSE(20)
      endif

c
      RETURN
      END

      SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(N,N),Z(N,N),D(*)
      real*8,allocatable::E(:)
      allocate(E(N))
      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
C
      SUBROUTINE READFF(N,FCAR)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(N,N)
      CONST=4.359828d0/0.5291772d0**2
      OPEN(20,FILE='FILE.FC',STATUS='OLD')
      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)
C
      DO 3 I=1,N
      DO 3 J=I,N
3     FCAR(J,I)=FCAR(J,I)*CONST
      DO 31 I=1,N
      DO 31 J=I+1,N
   31 FCAR(I,J)=FCAR(J,I)
      CLOSE(20)
      WRITE(*,*)' Cartesian FF read in ... '
      RETURN
      END

      SUBROUTINE WRITEFF(N,FCAR)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(N,N)
      CONST=4.359828d0/0.5291772d0**2
      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)/CONST,J=N1,MIN(LN,N3))
      N1=N1+5
      IF(N3.LT.N)GOTO 1
17    FORMAT(I4,5D14.6)
      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 PROJECT(F,N,C)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION F(N,N),C(3,N/3)
      real*8,allocatable::TEM(:,:),A(:,:)
C
      allocate(TEM(N,N),A(N,N))
      CALL DOMA(A,N,C)
      DO 11 I=1,N
      DO 11 J=1,N
11    TEM(I,J)=0.0d0

      DO 1 I=1,N
      DO 1 II=1,N
      if(dabs(A(I,II)).gt.1.0d-10)then
       DO 12 J=1,N
12     TEM(I,J)=TEM(I,J)+A(I,II)*F(II,J)
      endif
1     continue

      DO 22 I=1,N
      DO 22 J=1,N
22    F(I,J)=0.0d0

      DO 2 J=1,N
      DO 2 II=1,N
      if(dabs(A(J,II)).gt.1.0d-10)then
       DO 23 I=1,N
23     F(I,J)=F(I,J)+TEM(I,II)*A(J,II)
      endif
2     continue
      RETURN
      END

      SUBROUTINE DOMA(A,N,C)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(N,N),C(3,N/3)
      real*8,allocatable::TEM(:,:)
C
      N=N
      allocate(TEM(N,N))
      AMACH=0.00000000000001d0
      NAT=N/3
      DO 12 I=1,N
      DO 12 J=1,6
12    A(I,J)=0.0d0
      DO 1 I=1,NAT
      A(3*I-2,1)= 1.00d0
      A(3*I-1,2)= 1.00d0
      A(3*I  ,3)= 1.00d0
      A(3*I-2,4)=-C(2,I)
      A(3*I-1,4)= C(1,I)
      A(3*I-1,5)=-C(3,I)
      A(3*I  ,5)= C(2,I)
      A(3*I-2,6)= C(3,I)
1     A(3*I  ,6)=-C(1,I)
      N6=6
      CALL ORT(A,N6,NAT)
      ICOL=1
5     IF(ICOL.GT.N6)GOTO 8
      S=0.0d0
      DO 6 J=1,N
6     S=S+A(J,ICOL)**2
      IF(S.LT.1000.0d0*AMACH)THEN
       DO 7 J=1,N
7      A(J,ICOL)=A(J,N6)
       N6=N6-1
       GOTO 5
      ENDIF
      ICOL=ICOL+1
      GOTO 5
8     CONTINUE
      DO 2 I=1,N
      DO 2 J=1,N
      TEM(I,J)=0.0d0
      DO 2 II=1,N6
2     TEM(I,J)=TEM(I,J)+A(I,II)*A(J,II)
      DO 3 I=1,N
      DO 4 J=1,N
4     A(I,J)=-TEM(I,J)
3     A(I,I)=A(I,I)+1.0d0
      RETURN
      END
C     ============================================================ 
      SUBROUTINE ORT(A,N6,NAT)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(3*NAT,N6)
      AMACH=0.00000000000001d0
      NAT3=3*NAT
      C=0.0d0
      DO 1 I=1,NAT3
1     C=C+A(I,1)**2
      IF(C.GT.AMACH)C=1.0d0/SQRT(C)
      DO 2 I=1,NAT3
2     A(I,1)=A(I,1)*C
      DO 3 IC=2,N6
C     Column IC, orthonormalize to previous columns:
      DO 3 ICP=1,IC-1
      SP=0.0d0
      DO 5 J=1,NAT3
5     SP=SP+A(J,IC)*A(J,ICP)
C     Subtract projection to ICP:
      DO 6 J=1,NAT3
6     A(J,IC)=A(J,IC)-SP*A(J,ICP)
C     Normalize :
      SP=0.0d0
      DO 7 J=1,NAT3
7     SP=SP+A(J,IC)**2
      IF(SP.GT.AMACH)SP=1.0d0/SQRT(SP)
      DO 3 J=1,NAT3
3     A(J,IC)=A(J,IC)*SP
      RETURN
      END

C     ============================================================ 
      subroutine readinp(NOSC,SCFAC,NM,CHR,NOAT,NA,NNU,ZIM,
     1ZM,ZNU)
      implicit none
      integer*4 NOSC,NM(*),NOAT,I,J,NMOL,NNU,NA
      real*8 SCFAC(*),CHR(*),CHATOT,ZIM(*),BM,ZM(*),ZNU(*)
      OPEN(15,FILE='FTRY.INP')
      READ(15,*)
      READ(15,*)NOSC
      READ(15,*) (SCFAC(I),I=1,NOSC)
      READ(15,*) (NM(I),I=1,NOSC)
      WRITE(*,*)NOSC,' scale factors'
C
C     Read-in atomic charges for FPC:
c     READ(15,1040) (CHR(I),I=1,NOAT)
      READ(15,*) (CHR(I),I=1,NOAT)
      CHATOT=0.0d0
      DO 77 I=1,NOAT
77    CHATOT=CHATOT+CHR(I)
      WRITE(*,777)CHATOT
777   FORMAT(' TOTAL CHARGE',F12.2)
C
C     Read-in atomic masses:
      READ(15,*)NMOL
      IF(NMOL.NE.1)THEN
       CLOSE(16)
       CLOSE(15)
       call report(' Only one molecule allowed here !')
      ENDIF
      READ(15,1040)(ZIM(I),I=1,NOAT)
1040  FORMAT(6G12.6)

      BM=0.0d0
      DO 102 I=1,NOAT
102   BM=BM+ZIM(I)
      WRITE(*,778)BM
778   FORMAT(' TOTAL MASS  ',F20.6)
      DO 220 I=1,NOAT
      DO 220 J=1,3
220   ZM(3*I-3+J)=1/SQRT(ZIM(I))
C
C     Read-in experimental frequencies:
      NNU=0
      READ(15,*,end=99,err=99)NNU
      IF(NNU.EQ.0)NNU=NA-6
      READ(15,1040)(ZNU(I),I=1,NNU)
      DO 211 I=NNU+1,NA
211   ZNU(I)=0.0d0
99    CLOSE(15)
      WRITE(*,*)' FTRY.INP READ IN'
      return
      end

      subroutine rdx(NOAT,ic,ZQ,C)
      implicit none
      integer*4 NOAT,i,ix,ic,ZQ(*)
      real*8 C(3,NOAT)
      OPEN(17,FILE='FILE.X')
      READ(17,*)
      READ(17,*)NOAT
      if(ic.eq.0)goto 99
      do 17171 i=1,NOAT
17171 read(17,*)ZQ(i),(C(ix,i),ix=1,3)
99    CLOSE(17)
      if(ic.eq.1)write(6,*)NOAT,' atoms in FILE.X'
      return
      end

      subroutine readopt(lproject,IRANGE,JRANGE,DEL,NIT,spec,ity,hw,LG,
     1wmin,wmax,lconj)
      implicit none
      integer*4 IRANGE,JRANGE,NIT,ity
      real*8 DEL,hw,wmin,wmax
      logical lproject,LG,lconj
      character*80 spec
      character*4 key,kp
      spec='sp.prn'
      IRANGE=0
      JRANGE=0
      lproject=.true.
      LG=.false.
      lconj=.true.
      hw=15.0d0
      DEL=0.05d0
      NIT=10
      ity=3
      wmin=200
      wmax=1900.0d0
      open(9,file='VSCALE.OPT',status='old')
1     read(9,90,end=99,err=99)key
90    format(a4)
      if(key.eq.'PROJ')read(9,*)lproject
      if(key.eq.'RANG')read(9,*)IRANGE,JRANGE
      if(key.eq.'IRAN')read(9,*)IRANGE
      if(key.eq.'JRAN')read(9,*)JRANGE
      if(key.eq.'DELT')read(9,*)DEL
      if(key.eq.'MAXI')read(9,*)NIT
      if(key.eq.'SPEC')read(9,'(A)')spec
      if(key.eq.'TYPE')read(9,*)ity
      if(key.eq.'GAUS')read(9,*)LG
      if(key.eq.'WIDT')read(9,*)hw
      if(key.eq.'CONJ')read(9,*)lconj
      if(key.eq.'WMIN')read(9,*)wmin
      if(key.eq.'WMAX')read(9,*)wmax
      backspace 9
      read(9,90)kp
      if(kp.eq.key)call report(' Unkown option: '//key)
      goto 1
99    close(9)
      return
      end

      SUBROUTINE SMATOUT(S,AS,NAT3,i1,i2,ZQ,C)
C     dX  =  S dQ   S(nat3xnint)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      integer*4 ZQ(*)
      DIMENSION S(NAT3,NAT3),AS(3,NAT3),C(3,NAT3/3)
      IRANGE=1
      if(i1.ne.0)IRANGE=i1
      JRANGE=NAT3
      if(i2.ne.0)JRANGE=i2
      NI=JRANGE-IRANGE+1
c
c     look at the phase and make the biggest element always positive:
      DO 101 J=IRANGE,JRANGE
      amx=S(1,J)
      DO 102 I=2,NAT3
102   if(abs(S(I,J)).gt.abs(amx))amx=S(I,J)
      if(amx.lt.0)then
       do 103 I=1,NAT3
103    S(I,J)=-S(I,J)
      endif
101   continue
c
      OPEN(34,FILE='F.INP')
      WRITE(34,10)NI,NAT3,NAT3/3
10    FORMAT(3I7)
      DO 3 I=1,NAT3/3
3     WRITE(34,11)ZQ(I),C(1,I),C(2,I),C(3,I)
11    FORMAT(I7,3F12.6)
      WRITE(34,14)
14    FORMAT(' Atom Mode    X-disp.    Y-disp.      Z-disp.')
      DO 1 I=1,NAT3/3
      DO 1 J=IRANGE,JRANGE
      IU=3*(I-1)
      K=NAT3-J+1
      WRITE(34,15)I,K,S(IU+1,J),S(IU+2,J),S(IU+3,J)
15    FORMAT(2I7,3F11.6)
1     continue
      WRITE(34,11)NI
      do 4 I=IRANGE,JRANGE
      WRITE(34,13)AS(1,I)
13    format(F11.3,$)
      if(mod(I,6).eq.0)write(34,*)
4     continue
      write(34,*)
      CLOSE(34)
      WRITE(*,*)' S-matrix written out ...'
      RETURN
      END

      subroutine rsp(ic,np,sx,s,spec)
      implicit none
      integer*4 ic,np
      real*8 sx(np),s(np),x,y
      character*80 spec
      open(9,file=spec)
      np=0
1     read(9,*,end=99,err=99)x,y
      np=np+1
      if(ic.eq.1)then
       sx(np)=x
       s(np)=y
      endif
      goto 1
99    close(9)
      return
      end

      subroutine mksp(N,ity,np,sx,se,AS,SMAT,C,hw,LG)
      implicit none
      integer*4 ity,np,N
      real*8 sx(np),se(np),AS(3,N),SMAT(N,N),C(3,N/3),hw
      real*8,allocatable::iint(:)
      logical LG
      allocate(iint(N))
      if(ity.eq.3.or.ity.eq.8)then
       call new66(N,SMAT,AS,C,iint,ity)
      else
       write(6,*)ity
       call report('unknown type')
      endif
      call tabprnff(N,iint,AS,LG,hw,np,se,sx,ity)
      return
      end

      subroutine tabprnff(N,t,AS,LG,hw,np,s,sx,ic)
      implicit none
      integer*4 ic,i,i1,N,np
      real*8 x,dd,hw,const,wo,sf,temp,AS(3,N),t(N),s(np),sx(np)
      logical LG
      temp=300.0d0
      s=0.0d0
      if(LG)then
       dd=hw/2.0d0/dsqrt(log(2.0d0))
      else
       dd=hw/2.0d0
      endif
      do 4 i=1,np
      x=sx(i)
      do 4 i1=1,N
      wo=AS(1,i1)
      const=0.0d0
      if(ic.eq.1)const=108.7d0*t(i1)*wo
      if(ic.eq.2)const=435.0d0*t(i1)*wo
      if(ic.eq.3)const=t(i1)/wo/(1.0d0-exp(-(1.44d0*wo/temp)))
      if(ic.eq.8)const=t(i1)/wo/(1.0d0-exp(-(1.44d0*wo/temp)))
4     s(i)=s(i)+const*sf(dd,LG,x,wo)
      return
      end
c     ==================================================================
      function sf(d,g,x,x0)
      implicit none
      logical g
      real*8 pi,spi,sf,d,x,x0,dd
      pi=4.0d0*atan(1.0d0)
      spi=dsqrt(pi)
      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

      subroutine new66(N,SI,AS,R,iint,ity)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      integer*4 ity,N
      real*8 SI(N,N),AS(3,N),R(3,N/3)
      real*8 iint(*)
      real*8,allocatable::S(:,:),ALPHA(:,:,:),A(:,:,:,:),
     1G(:,:,:),tem(:,:),E(:)
      LOGICAL RAMAN,LPOLAR,lpro,lgpro,lapro,
     1luseA,luseG,lintA,lintG
C
      CM=219470.0d0
      SFAC=0.0234280d0
      CALL IOPAR(RAMAN,LPOLAR,EXCA,
     1wmin,wmax,lpro,lgpro,lapro,N6,N7,IDIF,ri,
     1luseA,luseG,lintA,lintG)
      NAT=N/3
      NM=3*NAT
      allocate(S(N,N),ALPHA(N,3,3),A(N,3,3,3),G(N,3,3),
     1tem(3,NAT),E(N))
      CALL READTEN(N,NAT,ALPHA,G,A)
      do 501 I=1,N
501   E(I)=AS(1,I)/CM
      DO 502 I=1,N
      DO 502 J=1,N
502   S(I,J)=SI(I,J)*SFAC
      IF(LPOLAR)THEN
c      set A, G to zero:
       CALL TTTT(N,ALPHA,A,G,NAT,3,R,tem)
c      make common A, G from alpha:
       CALL TTTT(N,ALPHA,A,G,NAT,2,R,tem)
       write(6,*)' A, G made from alpha !!!'
      ENDIF
      if(.not.lintA.or..not.lintG)then
c      set internal A, G to zero:
c      1) transform A, G to atomic origins:
       CALL TTTT(N,ALPHA,A,G,NAT,1,R,tem)
c      2) delete them
       if(.not.lintG)CALL TTTT(N,ALPHA,A,G,NAT,4,R,tem)
       if(.not.lintA)CALL TTTT(N,ALPHA,A,G,NAT,5,R,tem)
c      3) transform A, G back to common origin:
       CALL TTTT(N,ALPHA,A,G,NAT,2,R,tem)
      endif
C     transform the tensors into internal coordinate derivatives
      CALL TRTEN(S,3*NAT,ALPHA,A,G)
      CALL DORAMAN(NM,ALPHA,G,A,EXCA,wmin,wmax,E,
     1luseA,luseG,iint,ity)
      END
C     
      FUNCTION EPS(I,J,K)
      IMPLICIT INTEGER*4 (I-N)
      REAL*8 EPS
      EPS=0.0D0
      IF(I.EQ.1)then
       if(J.EQ.2.AND.K.EQ.3)EPS= 1.0D0
       IF(J.EQ.3.AND.K.EQ.2)EPS=-1.0D0
      endif
      IF(I.EQ.2)then
       IF(J.EQ.3.AND.K.EQ.1)EPS= 1.0D0
       IF(J.EQ.1.AND.K.EQ.3)EPS=-1.0D0
      endif
      IF(I.EQ.3)then
       IF(J.EQ.1.AND.K.EQ.2)EPS= 1.0D0
       IF(J.EQ.2.AND.K.EQ.1)EPS=-1.0D0
      endif
      RETURN
      END
C
      SUBROUTINE READTEN(N,NAT,ALPHA,G,A)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3)
      OPEN(2,FILE='FILE.TTT',STATUS='OLD')
      READ(2,*)
      READ(2,*)NAT
      if(NAT.gt.N/3)call report('too many atoms')
      READ(2,*)
      READ(2,*)
      DO 1 I=1,3
      READ(2,*)
      DO 1 L=1,NAT
      DO 1 IX=1,3
      IND=3*(L-1)+IX
1     READ(2,*)(ALPHA(IND,I,J),J=1,2),(ALPHA(IND,I,J),J=1,3)
      READ(2,*)
      READ(2,*)
      DO 2 I=1,3
      READ(2,*)
      DO 2 L=1,NAT
      DO 2 IX=1,3
      IND=3*(L-1)+IX
c     last index magnetic, inner
2     READ(2,*)(G(IND,I,J),J=1,2),(G(IND,I,J),J=1,3)
      READ(2,*)
      READ(2,*)
      DO 3 I=1,3
      DO 3 J=1,3
      READ(2,*)
      DO 3 L=1,NAT
      DO 3 IX=1,3
      IND=3*(L-1)+IX
3     READ(2,*)(A(IND,I,J,K),K=1,2),(A(IND,I,J,K),K=1,3)
      CLOSE(2)
      RETURN
      END
C
      SUBROUTINE IOPAR(RAMAN,LPOLAR,EXCA,
     1wmin,wmax,lpro,lgpro,lapro,N6,N7,IDIF,ri,
     1luseA,luseG,lintA,lintG)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      LOGICAL RAMAN,LPOLAR,lpro,lex,lgpro,lapro,
     1luseA,luseG,lintA,lintG
      CHARACTER*15 STRING
c
      N7=7
      IDIF=-1
      wmin=20.0d0
      wmax=1.0d90
      lpro=.false.
      lapro=.false.
      lgpro=.false.
      RAMAN=.FALSE.
      N6=3
      LPOLAR=.FALSE.
c     excitation light wavelength in nm:
      EXCNM=532.0d0
c     refractive index
      ri=1.0d0
c     use A-tensor derivatives:
      luseA=.true.
c     use G'-tensor derivatives:
      luseG=.true.
c     use intrinsic part of A-tensor derivatives:
      lintA=.true.
c     use intrinsic part of G'-tensor derivatives:
      lintG=.true.

      INQUIRE(FILE='ROA.OPT',exist=lex)
      if(lex)then
       WRITE(6,*)' FILE ROA.OPT FOUND ...'
       OPEN(2,FILE='ROA.OPT')
      else
       INQUIRE(FILE='NEW6.OPT',exist=lex)
       if(lex)then
        OPEN(2,FILE='NEW6.OPT')
        WRITE(6,*)' FILE NEW6.OPT FOUND ...'
       else
        goto 99
       endif
      endif
   
2001  READ(2,1111)STRING
1111  FORMAT(A15)
      IF(STRING(1:5).EQ.'POLAR')READ(2,*)LPOLAR
      IF(STRING(1:4).EQ.'WMIN')READ(2,*)wmin
      IF(STRING(1:4).EQ.'WMAX')READ(2,*)wmax
      IF(STRING(1:5).EQ.'EXCNM')READ(2,*)EXCNM
      IF(STRING(1:5).EQ.'RAMAN')READ(2,*)RAMAN
      IF(STRING(1:3).EQ.'PRO')READ(2,*)lpro
      IF(STRING(1:2).EQ.'N6')READ(2,*)N6
      IF(STRING(1:4).EQ.'GPRO')READ(2,*)lgpro
      IF(STRING(1:4).EQ.'APRO')READ(2,*)lapro
      IF(STRING(1:4).EQ.'REFR')READ(2,*)ri
      IF(STRING(1:4).EQ.'INTA')READ(2,*)lintA
      IF(STRING(1:4).EQ.'INTG')READ(2,*)lintG
      IF(STRING(1:4).EQ.'USEA')READ(2,*)luseA
      IF(STRING(1:4).EQ.'USEG')READ(2,*)luseG
      IF(STRING(1:3).EQ.'END')GOTO 2002
      GOTO 2001
2002  CLOSE(2)
     
99    EXCA=EXCNM*10.0D0
      RETURN
      END


      SUBROUTINE DORAMAN(N,ALPHA,G,A,EXCA,wmin,wmax,E,
     1luseA,luseG,iint,ity)
c
c     supposedly G/W is supplied
c
      IMPLICIT none
      integer*4 IQ,I,J,K,L,N,ity
      real*8 ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),E(*),EXCA,wmin,wmax,
     1CM,ECM,SAL0,SAL1,SAG0,SAG1,SA1,SPAL,SDCP180,SDCPII180,
     1S180,D180,D,S90X,D90X,DX,S90Z,D90Z,DZ,YDY,YDX,down,doc,
     1ro,AMU,BOHR,SpA3,beta2,gpisvejc,betaa,iint(N),
     1betag,roa1,
     1ram1,EPS
      logical luseA,luseG
c
      CM=219470.0d0
      AMU=1822.0d0
      BOHR=0.529177d0
      if(.not.luseA)A=0.0d0
      if(.not.luseG)G=0.0d0
c
      DO 1 IQ=1,N
      iint(IQ)=0.0d0
      ECM=E(IQ)*CM
      IF(ABS(ECM).GT.0.01d0.and.(wmin.eq.0.0d0.or.ECM.gt.wmin).and.
     1(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,DCP_II
       SDCP180  =6.0d0*(6.0d0*SAL1-2.0d0*SAL0)
       SDCPII180=6.0d0*(      SAL1+3.0d0*SAL0)
c
c      IL+/-IR, backscattering
       S180=7.0d0*SAL1+SAL0
c      write(6,*)ECM,SAL0,SAL1,SAG0,SAG1,SA1
       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
      
c      calculate degree of circularity
       down=        7.0d0*SAL1+ 1.0d0*SAL0
       doc=5.0d0*(  1.0d0*SAL1- 1.0d0*SAL0)
       if(down.ne.0.0d0)doc=doc/down

c      calculate depolarization
       down= 6.0d0*SAL0-12.0d0*SAL1
       ro  =-3.0d0*SAL0+ 9.0d0*SAL1
       if(down.ne.0.0d0)ro=ro/down
c
       SpA3=SPAL/3.0d0*dsqrt((AMU*BOHR**4))
       beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4)
       SDCP180  =SDCP180  *(AMU*BOHR**4)
       SDCPII180=SDCPII180*(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      Barron's tensor invariants:
c      alpha G = [(1/3)Sp alpha] [(1/3)Sp G'] etc
c      beta(A)^2=Delta2 in Gaussian:
       betaa=SA1*3.0d0/2.0d0            *gpisvejc
       betag=(3.0d0*SAG1-SAG0)/2.0d0    *gpisvejc
c
c      ICP/SCP(180):
c      8/c*(12*beta(G')^2 + 4*beta(A)^2):
       roa1=96.0d0*(betag+betaa/3.0d0)
c      ram1 = 6(7 al_ab al_ab + al_aa al_bb)
       ram1=4.0d0*(45.0d0*SpA3**2+7.0d0*beta2)
       if(ity.eq.3)then
        iint(IQ)=ram1
       else
        iint(IQ)=roa1
       endif
      ENDIF
1     CONTINUE
      RETURN
      END

      SUBROUTINE TRTEN(S,N,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))
      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)
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)
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)
12    A(IQ,I,J,K)=AS
      RETURN
      END
c     
      SUBROUTINE TTTT(N,ALPHA,A,G,NAT,ICO,R,X)
c
C     ICO = 1: Transforms A and G to local origins
C     ICO = 2: Transforms A and G to common origin
C     ICO = 3: Sets A and G to zero 
c     ICO = 4: Delete G
c     ICO = 5: Delete A
C     R supposed to be passed in AU
C
      IMPLICIT none
      INTEGER*4 N,NAT,ICO,I,J,L,IE,II,K,M,IA,IB,IC,ID
      real*8 ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),R(3,N/2),X(3,N/3),
     1S,EPS,SUM
C
      DO 6 I=1,3
      DO 6 J=1,NAT
6     X(I,J)=R(I,J)
C
      IF(ICO.EQ.1)then
       S= 1.0d0
       WRITE(6,*)' A,G transformed into atomic origins'
      elseif(ICO.EQ.2)then
       S=-1.0d0
       WRITE(6,*)' A,G transformed into common origin'
      elseif(ICO.ge.3.and.ICO.LE.5)then
       IF(ICO.eq.3.or.ICO.EQ.4)then
        do 5 II=1,N
        do 5 I=1,3
        do 5 J=1,3
5       G(II,I,J)=0.0d0
        WRITE(6,*)' G set to zero'
       endif
       IF(ICO.eq.3.or.ICO.EQ.5)then
        do 7 II=1,N
        do 7 I=1,3
        do 7 J=1,3
        do 7 K=1,3
7       A(II,I,J,K)=0.0d0
        WRITE(6,*)' A set to zero'
       endif
       return
      else
       write(6,*)ICO
       call report('TTTT unknown instruction')
      endif
C
      DO 1 I=1,3
      DO 1 J=1,3
      DO 1 L=1,NAT
      DO 1 IE=1,3
      II=3*(L-1)+IE
      DO 1 K=1,3
      DO 1 M=1,3
c     G-second index-magnetic
1     G(II,I,J)=G(II,I,J)+S*0.5d0*EPS(J,K,M)*X(K,L)*ALPHA(II,I,M)
C
      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
      II=3*(L-1)+IE
      SUM=0.0d0
      IF(IB.EQ.IC)THEN
       DO 4 ID=1,3
4      SUM=SUM+X(ID,L)*ALPHA(II,IA,ID)
      ENDIF
      A(II,IA,IB,IC)=A(II,IA,IB,IC)+SUM*S
     1-1.5d0*S*(X(IB,L)*ALPHA(II,IA,IC)+X(IC,L)*ALPHA(II,IA,IB))
3     CONTINUE
c
      RETURN
      END
c
      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end

      subroutine ru(NOAT,ic,ni,li)
      implicit none
      integer*4 ic,ni,NOAT,li(3,ni),i,is
      is=0
      open(9,file='FILE.UMA')
      read(9,*)
      read(9,*)
      read(9,*)ni,ni
      if(ic.eq.0)goto 1
      do 2 i=1,NOAT
2     read(9,*)
      do 3 i=1,ni
      read(9,*)li(1,i),li(2,i),li(3,i)
      if(li(1,i).eq.1)is=is+1
      if(li(1,i).eq.4)then
       read(9,*)
       read(9,*)
      endif
3     continue
1     close(9) 
      if(ic.eq.1)write(6,600)ni,is
600   format(i6,' internal coordinates in FILE.UMA,',i6,' stretchings')
      return
      end

      subroutine rb(ni,li,b)
      implicit none
      integer*4 ni,li(3,ni),i,nt,j
      real*8 b(ni,6)
      open(9,file='B.MAT',status='old')
      read(9,*)
      read(9,*)
      read(9,*)
      do 1 i=1,ni
      read(9,*)
      read(9,*)nt
      if(li(1,i).eq.1)then
       if(nt.ne.6)call report('ne <> 6')
       do 2 j=1,nt
2      read(9,*)b(i,j),b(i,j),b(i,j)
      else
       do 3 j=1,nt
3      read(9,*)
      endif
1     continue
      close(9)
      write(6,*)'B.MAT read'
      return
      end

      SUBROUTINE rFF(N,fi)
      IMPLICIT none
      integer*4 N,N1,N3,LN,J
      real*8 fi(N),x
      open(9,file='INTY.FC',status='old')
      N1=1
1     N3=N1+4
      IF(N3.GT.N)N3=N
      read(9,*)
      DO 130 LN=N1,N
      read(9,*)x,(x,J=N1,MIN(LN,N3))
130   if(ln.eq.MIN(LN,N3))fi(LN)=x
c30   read(9,*)x,(xFCAR(LN,J),J=N1,MIN(LN,N3))
      N1=N1+5
      IF(N3.LT.N)GOTO 1
      close(9)
      RETURN
      END

      subroutine wrsp(np,sx,se,an1,an2)
      implicit none
      integer*4 np,i
      real*8 sx(*),se(*),an1,an2
      open(9,file='fit.prn')
      do 1 i=1,np
1     write(9,900)sx(i),se(i)/an2*an1
900   format(f12.2,e12.4)
      close(9)
      return
      end

      subroutine wrfi(ni,N,SCFAC,NM,NOAT,ZIM,CHR)
      implicit none
      integer*4 ni,N,NM(*),NOAT,i
      real*8 SCFAC(*),CHR(*),ZIM(*)
      open(7,file='FTRY.NEW.INP')
      write(7,900)ni
900   format('AUTO   0   1   0   1   0   1   0   1',i4)      
      WRITE(7,333)N
333   FORMAT(20I4)
      WRITE(7,334)(SCFAC(i),i=1,N)
334   FORMAT(6f12.6)
      WRITE(7,333)(NM(i),i=1,N)
      WRITE(7,334)(CHR(I),I=1,NOAT)
      WRITE(7,333)1
      WRITE(7,334)(ZIM(I),I=1,NOAT)
      WRITE(7,*)
      WRITE(7,334)(0.0d0,I=1,ni)
      close(7)
      return
      end