      SUBROUTINE PATHS
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      common /PATH  / LATOM,LPARAM,REACT(200)
      common /GEOVAR/ NVAR, LOC(2,MAXPAR), IDUMY, XPARAM(MAXPAR)
      common /KEYWRD/ KEYWRD
      common /TIME  / TIME0
      common /GEOM  / GEO(3,NUMATM)
      common /ALPARM/ ALPARM(3,MAXPAR),X0, X1, X2, ILOOP
************************************************************************
*
*   PATH FOLLOWS A REACTION COORDINATE.   THE REACTION COORDINATE IS ON
*        ATOM LATOM, AND IS A DISTANCE  IF LPARAM=1,
*                           AN ANGLE    IF LPARAM=2,
*                           AN DIHEDRAL IF LPARAM=3.
*
************************************************************************
      dimension GD(MAXPAR),XLAST(MAXPAR),MDFP(20),XDFP(20)
      character*241 KEYWRD
      character*10 TYPE(3)
      save TYPE, GD, XLAST, MDFP, XDFP
      data TYPE / 'ANGSTROMS ','DEGREES   ','DEGREES   '/
      ILOOP=1
      if (INDEX(KEYWRD,'RESTAR') .NE. 0) then
         MDFP(9)=0
         call DFPSAV(TOTIME,XPARAM,GD,XLAST,FUNCT1,MDFP,XDFP)
         write(MIOB,'(//10X,'' RESTARTING AT POINT'',I3)') ILOOP
      end if
      if (ILOOP.GT.1) go to 10
      write(MIOB,'(''  ABOUT TO ENTER FLEPO FROM PATH'')')
      TIME0=SECOND()
      call FLEPO(XPARAM,NVAR,FUNCT)
      write(MIOB,
     %  '(''  OPTIMIZED VALUES OF PARAMETERS, INITIAL POINT'')')
      call WRITMO(TIME0,FUNCT)
      TIME0=SECOND()
   10 continue
      if (ILOOP.GT.2) go to 40
      GEO(LPARAM,LATOM)=REACT(2)
      if (ILOOP.EQ.1) then
         X0=REACT(1)
         X1=X0
         X2=REACT(2)
c
         if (X2.LT. -100.D0) then
           call filclo  !  J/M
           stop ' in PATHS - REACT(2) < -100.0 '
         end if
c
         do 20 I=1,NVAR
            ALPARM(2,I)=XPARAM(I)
   20    ALPARM(1,I)=XPARAM(I)
         ILOOP=2
      end if
      call FLEPO(XPARAM,NVAR,FUNCT)
      RNORD=REACT(2)
      if (LPARAM.GT.1) RNORD=RNORD*57.29577951D0
      write(MIOB,'(1X,16(''*****'')//17X,''REACTION COORDINATE = ''
     %   ,F12.4,2X,A10,19X//1X,16(''*****''))')  RNORD,TYPE(LPARAM)
      write(*,'(1X,16(''*****'')/17X,''REACTION COORDINATE = ''
     %   ,F12.4,2X,A10,19X/1X,16(''*****''))')  RNORD,TYPE(LPARAM)
      call WRITMO(TIME0,FUNCT)
      TIME0=SECOND()
      do 30 I=1,NVAR
   30 ALPARM(3,I)=XPARAM(I)
C
C   NOW FOR THE MAIN INTERPOLATION ROUTE
C
      if (ILOOP.EQ.2) ILOOP=3
   40 continue
      LPR=ILOOP
      do 110 ILOOP = LPR,100
C
         if (REACT(ILOOP).LT. -100.D0) return
C
         RNORD=REACT(ILOOP)
         if (LPARAM.GT.1) RNORD=RNORD*57.29577951D0
         write(*,'(1X,16(''*****'')//19X,''REACTION COORDINATE = ''
     %       ,F12.4,2X,A10,19X//1X,16(''*****''))') RNORD,TYPE(LPARAM)
         write(MIOB,'(1X,16(''*****'')//19X,''REACTION COORDINATE = ''
     %       ,F12.4,2X,A10,19X//1X,16(''*****''))') RNORD,TYPE(LPARAM)
C
         X3=REACT(ILOOP)
         C3=(X0**2-X1**2)*(X1-X2)-(X1**2-X2**2)*(X0-X1)
C      write(MIOB,'(''   C3:'',F13.7)') C3
         if (ABS(C3) .LT. 1.D-8) then
C
C    WE USE A LINEAR INTERPOLATION
C
            CC1=0.D0
            CC2=0.D0
         else
C    WE do A QUADRATIC INTERPOLATION
C
            CC1=(X1-X2)/C3
            CC2=(X0-X1)/C3
         end if
         CB1=1.D0/(X1-X2)
         CB2=(X1**2-X2**2)*CB1
C
C    NOW TO CALCULATE THE INTERPOLATED COORDINATES
C
         do 50 I=1,NVAR
            DELF0=ALPARM(1,I)-ALPARM(2,I)
            DELF1=ALPARM(2,I)-ALPARM(3,I)
            ACONST = CC1*DELF0-CC2*DELF1
            BCONST = CB1*DELF1-ACONST*CB2
            CCONST = ALPARM(3,I) - BCONST*X2 - ACONST*X2**2
            XPARAM(I)=CCONST+BCONST*X3+ACONST*X3**2
            ALPARM(1,I)=ALPARM(2,I)
   50    ALPARM(2,I)=ALPARM(3,I)
C
C   NOW TO CHECK THAT THE GUESSED GEOMETRY IS NOT TOO ABSURD
C
         do 60 I=1,NVAR
   60    if (ABS(XPARAM(I)-ALPARM(3,I)) .GT. 0.2) go to 70
         go to 90
   70    write(MIOB,
     %     '('' GEOMETRY TOO UNSTABLE FOR EXTRAPOLATION TO BE USED''/
     1       '' - THE LAST GEOMETRY IS BEING USED TO START THE NEXT'',
     2       '' CALCULATION'')')
         do 80 I=1,NVAR
   80    XPARAM(I)=ALPARM(3,I)
   90    continue
         X0=X1
         X1=X2
         X2=X3
         GEO(LPARAM,LATOM)=REACT(ILOOP)
         call FLEPO(XPARAM,NVAR,FUNCT)
         call WRITMO(TIME0,FUNCT)
         TIME0=SECOND()
         do 100 I=1,NVAR
  100    ALPARM(3,I)=XPARAM(I)
  110 continue
      end
