      SUBROUTINE LINMIN(XPARAM,ALPHA,PVECT,NVAR,FUNCT,OKF,IC, DOTT)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension XPARAM(NVAR),PVECT(NVAR)
      common /GRAVEC/ COSINE
      common /NUMCAL/ NUMCAL
C*********************************************************************
C
C  LINMIN DOES A LINE MINIMISATION.
C
C  ON INPUT:  XPARAM = STARTING COORDINATE OF SEARCH.
C             ALPHA  = STEP SIZE FOR INITIATING SEARCH.
C             PVECT  = DIRECTION OF SEARCH.
C             NVAR   = NUMBER OF VARIABLES IN XPARAM.
C             FUNCT  = INITIAL VALUE OF THE FUNCTION TO BE MINIMIZED.
C             ISOK   = NOT IMPORTANT.
C             COSINE = COSINE OF ANGLE OF CURRENT AND PREVIOUS GRADIENT.
C
C  ON OUTPUT: XPARAM = COORDINATE OF MINIMUM OF FUNCTI0N.
C             ALPHA  = NEW STEP SIZE, USED IN NEXT call OF LINMIN.
C             FUNCT  = FINAL, MINIMUM VALUE OF THE FUNCTION.
C             OKF    = TRUE IF LINMIN IMPROVED FUNCT, FALSE OTHERWISE.
C
C**********************************************************************
      common /KEYWRD/ KEYWRD
C
C  THE FOLLOWING common IS USED TO FIND OUT IF A NON-VARIATIONALLY
C  OPTIMIZED WAVE-FUNCTION IS BEING USED.
C
      common /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      character KEYWRD*241
      dimension PHI(3), VT(4), POLD(MAXPAR)
      dimension XSTOR(MAXPAR), XPAREF(MAXPAR)
      INTEGER LEFT,RIGHT,CENTER
      logical PRINT,OKF, HALFE, DIIS
      save ICALCN, PRINT, HALFE, XMAXM, I
      save MAXLIN, ANGLE, YMAXST
      data ICALCN /0/
      if (ICALCN.NE.NUMCAL) then
         HALFE =(INDEX(KEYWRD,'C.I.') .NE. 0 .OR. NCLOSE.NE.NOPEN)
         if (INDEX(KEYWRD,'GNORM') .NE. 0)
     1         DROP=DROP*MIN(READA(KEYWRD,INDEX(KEYWRD,'GNORM')),1.D0)
         XMAXM  = 0.4D0
         DELTA2 = 0.001D0
         if (INDEX(KEYWRD,'NOTH') .EQ. 0) then
            DELTA1 = 0.5D0
         else
            DELTA1 = 0.1D0
         end if
         ALPHA  = 1.D0
         MAXLIN = 15
         if (NVAR.EQ.1) then
            PVECT(1)=0.01D0
            DROP=0.01D0
            ALPHA=1.D0
            DELTA1 = 0.00005D0
            DELTA2 = 0.00001D0
            if (INDEX(KEYWRD,'PREC') .NE. 0) DELTA1=0.0000005
            MAXLIN=30
         end if
         ANGLE=0.8D0
         if (HALFE) ANGLE=-2.D0
         COSINE=99.99D0
C
C  ANGLE IS USED TO DECIDE IF P IS TO BE UPDATED AS CALCULATION
C        PROCEEDS.
C
         YMAXST  = 0.4D0
         PRINT=(INDEX(KEYWRD,'LINMIN') .NE. 0)
         ICALCN=NUMCAL
      end if
C#         write(MIOB,'(A,F12.6)') ' ANGLE BETWEEN PVECT''S:',
C#     % DOT(POLD,PVECT,NVAR)/
C#     % SQRT(DOT(POLD,POLD,NVAR)*DOT(PVECT,PVECT,NVAR)+1.D-20)
      do 10 I=1,NVAR
         XPAREF(I)=XPARAM(I)
   10 POLD(I)=PVECT(I)
      XMAXM=0.D0
      do 20 I=1,NVAR
         PABS=ABS(PVECT(I))
   20 XMAXM=MAX(XMAXM,PABS)
      XMINM=XMAXM
      XMAXM=YMAXST/XMAXM
      if (NVAR.EQ.1)
     1   call COMPFG(XPARAM, .TRUE., FUNCT,.TRUE.,GRAD,.FALSE.)
      FIN=FUNCT
      SSQLST=FUNCT
      DIIS=IC.EQ.1.AND.NVAR.GT.1
      PHI(1)=FUNCT
      ALPHA=1.D0
      VT(1)=0.0D00
      VT(2)=ALPHA
      if (VT(2).GT.XMAXM) VT(2)=XMAXM
      FMAX=FUNCT
      FMIN=FUNCT
      ALPHA=VT(2)
      do 30 I=1,NVAR
   30 XPARAM(I)=XPAREF(I)+ALPHA*PVECT(I)
      call COMPFG(XPARAM, .TRUE., PHI(2),.TRUE.,GRAD,.FALSE.)
      if (PHI(2).GT.FMAX) FMAX=PHI(2)
      if (PHI(2).LT.FMIN) FMIN=PHI(2)
      call EXCHNG (PHI(2),SQSTOR,ENERGY,ESTOR,XPARAM,XSTOR,
     1             ALPHA,ALFS,NVAR)
      if (DIIS) go to 190
      if (NVAR.GT.1) then
C
C   CALCULATE A NEW ALPHA BASED ON THIEL'S FORMULA
C
         ALPHA=-ALPHA**2*DOTT/(2.D0*(PHI(2)-SSQLST-ALPHA*DOTT))
         if (ALPHA.GT.2.D0) ALPHA=2.D0
      else
         if (PHI(2).LT.PHI(1)) then
            ALPHA=2*ALPHA
         else
            ALPHA=-ALPHA
         end if
      end if
C#      if (PRINT) write(MIOB,'(3(A,F12.6))')
C#     %    ' ESTIMATED DROP:',DOTT*0.5D0,
C#     %    '  ACTUAL: ',PHI(2)-SSQLST, '  PREDICTED ALPHA',ALPHA
      OKF=OKF.OR.PHI(2).LT.SSQLST
      if (DELTA1.GT.0.3D0) then
C
C  THIEL'S TESTS # 18 AND 19
C
         if (OKF.AND.ALPHA.LT.2.D0) go to 190
      end if
      VT(3)=ALPHA
      if (VT(3).LE.1.D0) then
         LEFT=3
         CENTER=1
         RIGHT=2
      else
         LEFT=1
         CENTER=2
         RIGHT=3
      end if
      do 40 I=1,NVAR
   40 XPARAM(I)=XPAREF(I)+ALPHA*PVECT(I)
      call COMPFG (XPARAM, .TRUE., FUNCT,.TRUE.,GRAD,.FALSE.)
      if (FUNCT.GT.FMAX) FMAX=FUNCT
      if (FUNCT.LT.FMIN) FMIN=FUNCT
      if (FUNCT.LT.SQSTOR) call EXCHNG (FUNCT,SQSTOR,ENERGY,
     1                          ESTOR,XPARAM,XSTOR,ALPHA,ALFS,NVAR)
      OKF=(OKF.OR.FUNCT.LT.FIN)
      PHI(3)=FUNCT
      if (PRINT) write(MIOB,50) VT(1),PHI(1),PHI(1)-FIN,
     %                          VT(2),PHI(2),PHI(2)-FIN,
     %                          VT(3),PHI(3),PHI(3)-FIN
   50 format( ' ---QLINMN ',/5X, 'LEFT   ...',F17.8,2F17.11/5X,
     1 'CENTER ...',F17.8,2F17.11,/5X, 'RIGHT  ...',F17.8,2F17.11,/)
      do 180 ICTR=3,MAXLIN
         ALPHA=VT(2)-VT(3)
         BETA=VT(3)-VT(1)
         GAMMA=VT(1)-VT(2)
         if (ABS(ALPHA*BETA*GAMMA) .GT. 1.D-4) then
            ALPHA=-(PHI(1)*ALPHA+PHI(2)*BETA+PHI(3)*GAMMA)/
     1                                      (ALPHA*BETA*GAMMA)
         else
C
C   FINISH BECAUSE TWO POINTS CALCULATED ARE VERY CLOSE TOGETHER
C
            go to 190
         end if
         BETA=((PHI(1)-PHI(2))/GAMMA)-ALPHA*(VT(1)+VT(2))
         if (ALPHA) 60,60,90
   60    if (PHI(RIGHT).GT.PHI(LEFT)) go to 70
         ALPHA=3.0D00*VT(RIGHT)-2.0D00*VT(CENTER)
         go to 80
   70    ALPHA=3.0D00*VT(LEFT)-2.0D00*VT(CENTER)
   80    S=ALPHA-ALPOLD
         if (ABS(S).GT.XMAXM) S=SIGN(XMAXM,S)*(1+0.01*(XMAXM/S))
         ALPHA=S+ALPOLD
         go to 100
   90    ALPHA=-BETA/(2.0D00*ALPHA)
         S=ALPHA-ALPOLD
         XXM=2.0D00*XMAXM
         if (ABS(S).GT.XXM) S=SIGN(XXM,S)*(1+0.01*(XXM/S))
         ALPHA=S+ALPOLD
  100    continue
C
C   FINISH if CALCULATED POINT IS NEAR TO POINT ALREADY CALCULATED
C
         do 110 I=1,3
  110    if (ABS(ALPHA-VT(I)).LT.DELTA1*(1.D0+VT(I)).AND.OKF) go to 190
         do 120 I=1,NVAR
  120    XPARAM(I)=XPAREF(I)+ALPHA*PVECT(I)
         FUNOLD=FUNCT
         call COMPFG (XPARAM, .TRUE., FUNCT,.TRUE.,GRAD,.FALSE.)
         if (FUNCT.GT.FMAX) FMAX=FUNCT
         if (FUNCT.LT.FMIN) FMIN=FUNCT
         if (FUNCT.LT.SQSTOR) call EXCHNG (FUNCT,SQSTOR,ENERGY,ESTOR,
     1                                     XPARAM,XSTOR,ALPHA,ALFS,NVAR)
         OKF=OKF .OR. (FUNCT.LT.FIN)
         if (PRINT) write(MIOB,130) VT(LEFT),PHI(LEFT), PHI(LEFT)-FIN,
     1                           VT(CENTER),PHI(CENTER),PHI(CENTER)-FIN,
     2                           VT(RIGHT),PHI(RIGHT),PHI(RIGHT)-FIN,
     3                           ALPHA,FUNCT,FUNCT-FIN
  130    format(5X,'LEFT    ...',F17.8,2F17.11,/5X,'CENTER  ...',
     1F17.8,2F17.11,/5X,'RIGHT   ...',F17.8,2F17.11,/5X,
     2 'NEW     ...',F17.8,2F17.11,/)
C
C TEST TO EXIT FROM LINMIN IF NOT DROPPING IN VALUE OF FUNCTION FAST.
C
         if (ABS(FUNOLD-FUNCT) .LT. DELTA2 .AND. OKF) go to 190
         ALPOLD=ALPHA
         if ((ALPHA.GT.VT(RIGHT)).OR.(ALPHA.GT.VT(CENTER)
     1        .AND.FUNCT.LT.PHI(CENTER)).OR.(ALPHA.GT.VT(LEFT)
     2        .AND.ALPHA.LT.VT(CENTER).AND.FUNCT.GT.PHI(CENTER)))
     3         go to 140
         VT(RIGHT)=ALPHA
         PHI(RIGHT)=FUNCT
         go to 150
  140    VT(LEFT)=ALPHA
         PHI(LEFT)=FUNCT
  150    if (VT(CENTER).LT.VT(RIGHT)) go to 160
         I=CENTER
         CENTER=RIGHT
         RIGHT=I
  160    if (VT(LEFT).LT.VT(CENTER)) go to 170
         I=LEFT
         LEFT=CENTER
         CENTER=I
  170    if (VT(CENTER).LT.VT(RIGHT)) go to 180
         I=CENTER
         CENTER=RIGHT
         RIGHT=I
  180 continue
  190 continue
C
C  IC=1 if THE LAST POINT CALCULATED WAS THE BEST POINT, IC=2 OTHERWISE
C
      IC=2
      if (ABS(ESTOR-ENERGY).LT.1.D-12) IC=1
      call EXCHNG (SQSTOR,FUNCT,ESTOR,ENERGY,XSTOR,XPARAM,
     1             ALFS,ALPHA,NVAR)
      OKF = (FUNCT.LT.SSQLST.OR.DIIS)
      if (FUNCT.GE.SSQLST) return
      if (ALPHA) 200,220,220
  200 ALPHA=-ALPHA
      do 210 I=1,NVAR
  210 PVECT(I)=-PVECT(I)
  220 continue
      return
C
C
      end
