      SUBROUTINE POWSQ(XPARAM, NVAR, FUNCT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION XPARAM(*)
      COMMON /MESAGE/ IFLEPO,ISCF
**********************************************************************
*
*   POWSQ OPTIMIZES THE GEOMETRY BY MINIMISING THE GRADIENT NORM.
*         THUS BOTH GROUND AND TRANSITION STATE GEOMETRIES CAN BE
*         CALCULATED. IT IS ROUGHLY EQUIVALENT TO FLEPO, FLEPO MINIMIZES
*         THE ENERGY, POWSQ MINIMIZES THE GRADIENT NORM.
*
*  ON ENTRY XPARAM = VALUES OF PARAMETERS TO BE OPTIMIZED.
*           NVAR   = NUMBER OF PARAMETERS TO BE OPTIMIZED.
*
*  ON EXIT  XPARAM = OPTIMIZED PARAMETERS.
*           FUNCT  = HEAT OF FORMATION IN KCALS.
*
**********************************************************************
C        *****  ROUTINE PERFORMS  A LEAST SQUARES MINIMIZATION  *****
C        *****  OF A FUNCTION WHICH IS A SUM OF SQUARES.        *****
C        *****  INITIALLY WRITTEN BY J.W. MCIVER JR. AT SUNY/   *****
C        *****  BUFFALO, SUMMER 1971.  REWRITTEN AND MODIFIED   *****
C        *****  BY A.K. AT SUNY BUFFALO AND THE UNIVERSITY OF   *****
C        *****  TEXAS.  DECEMBER 1973                           *****
C
      COMMON /GEOVAR/ NDUM,LOC(2,MAXPAR), IDUMY, XARAM(MAXPAR)
      COMMON /GEOM  / GEO(3,NUMATM)
      COMMON /LAST  / LAST
      COMMON /KEYWRD/ KEYWRD
      COMMON /TIME  / TIME0
      COMMON /NUMSCF/ NSCF
      COMMON /GEOSYM/ NDEP, LOCPAR(MAXPAR), IDEPFN(MAXPAR),
     1                 LOCDEP(MAXPAR)
      COMMON /GRADNT/ GRAD(MAXPAR),GNFINA
      COMMON /TIMDMP/ TLEFT, TDUMP
      COMMON /GEOKST/ NATOMS,LABELS(NUMATM),
     1                NA(NUMATM),NB(NUMATM),NC(NUMATM)
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /NUMCAL/ NUMCAL
      COMMON /SIGMA1/ GNEXT, AMIN, ANEXT
      COMMON /SIGMA2/ GNEXT1(MAXPAR), GMIN1(MAXPAR)
      COMMON /NLLCOM/ HESS(MAXPAR,MAXPAR),BMAT(MAXPAR,MAXPAR),
     1PMAT(MAXPAR*MAXPAR)
      COMMON /SCRACH/ PVEC
      DIMENSION IPOW(9), SIG(MAXPAR),
     1          E1(MAXPAR), E2(MAXPAR),
     2          P(MAXPAR), WORK(MAXPAR),
     3          PVEC(MAXPAR*MAXPAR), EIG(MAXPAR), Q(MAXPAR)
      LOGICAL DEBUG, RESTRT, TIMES, OKF, ROUGH, SCF1, RESFIL, LOG
      CHARACTER*241 KEYWRD
      SAVE ICALCN
      DATA  ICALCN /0/
      IF(ICALCN.NE.NUMCAL) THEN
         ICALCN=NUMCAL
         RESTRT=(INDEX(KEYWRD,'RESTART') .NE. 0)
         LOG=(INDEX(KEYWRD,'NOLOG') .EQ. 0)
         SCF1=(INDEX(KEYWRD,'1SCF') .NE. 0)
         ROUGH=.FALSE.
         TIME1=SECOND()
         TIME2=TIME1
         ICYC=0
         TIMES=(INDEX(KEYWRD,'TIME') .NE. 0)
         TLAST=TLEFT
         RESFIL=.FALSE.
         STEP=0.02D0
         LAST=0
         ILOOP=1
         NAT3=NUMAT*3
         XINC=0.00529167D0
         RHO2=1.D-4
         TOL2=4.D-1
         IF(INDEX(KEYWRD,'PREC') .NE. 0) TOL2=1.D-2
         IF(INDEX(KEYWRD,'GNORM') .NE. 0) THEN
            TOL2=READA(KEYWRD,INDEX(KEYWRD,'GNORM'))
            IF(TOL2.LT.0.01.AND.INDEX(KEYWRD,' LET').EQ.0)THEN
               WRITE(6,'(/,A)')'  GNORM HAS BEEN SET TOO LOW, RESET TO 0
     1.01'
               TOL2=0.01D0
            ENDIF
         ENDIF
         DEBUG = (INDEX(KEYWRD,'POWSQ') .NE. 0)
         IF(RESTRT) THEN
C
C   RESTORE STORED DATA
C
            IPOW(9)=0
            CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,ILOOP,BMAT,IPOW)
            IF(SCF1) GOTO 390
            NSCF=IPOW(8)
            DO 10 I=1,NVAR
               GRAD(I)=GMIN1(I)
   10       GNEXT1(I)=GMIN1(I)
            WRITE(6,'('' XPARAM'',6F10.6)')(XPARAM(I),I=1,NVAR)
            IF(ILOOP .GT. 0) THEN
C#               ILOOP=ILOOP+1
               WRITE(6,'(//10X,'' RESTARTING AT POINT'',I3)')ILOOP
            ELSE
               WRITE(6,'(//10X,''RESTARTING IN OPTIMISATION'',
     1         '' ROUTINES'')')
            ENDIF
         ENDIF
*
*   DEFINITIONS:   NVAR   = NUMBER OF GEOMETRIC VARIABLES = 3*NUMAT-6
*
      ENDIF
      NVAR=ABS(NVAR)
      IF(DEBUG) THEN
         WRITE(6,'('' XPARAM'')')
         WRITE(6,'(5(2I3,F10.4))')(LOC(1,I),LOC(2,I),XPARAM(I),I=1,NVAR)
      ENDIF
      IF( .NOT. RESTRT) THEN
         DO 20 I=1,NVAR
   20    GRAD(I)=0.D0
         CALL COMPFG(XPARAM, .TRUE., FUNCT, .TRUE., GRAD, .TRUE.)
      ENDIF
      IF(DEBUG) THEN
         WRITE(6,'('' STARTING GRADIENTS'')')
         WRITE(6,'(3X,8F9.4)')(GRAD(I),I=1,NVAR)
      ENDIF
      GMIN=SQRT(DOT(GRAD,GRAD,NVAR))
      GLAST=GMIN
      DO 30 I=1,NVAR
         GNEXT1(I)=GRAD(I)
         GMIN1(I)=GNEXT1(I)
   30 CONTINUE
C
C    NOW TO CALCULATE THE HESSIAN MATRIX.
C
      IF(ILOOP.LT.0) GOTO 140
C
C   CHECK THAT HESSIAN HAS NOT ALREADY BEEN CALCULATED.
C
      ILPR=ILOOP
      DO 50 ILOOP=ILPR,NVAR
         TIME1=SECOND()
         XPARAM(ILOOP)=XPARAM(ILOOP) + XINC
         CALL COMPFG(XPARAM, .TRUE., FUNCT, .TRUE., GRAD, .TRUE.)
         IF(SCF1) GOTO 390
         IF(DEBUG)WRITE(6,'(I3,12(8F9.4,/3X))')
     1    ILOOP,(GRAD(IF),IF=1,NVAR)
         GRAD(ILOOP)=GRAD(ILOOP)+1.D-5
         XPARAM(ILOOP)=XPARAM(ILOOP) - XINC
         DO 40 J=1,NVAR
   40    HESS(ILOOP,J)=-(GRAD(J)-GNEXT1(J))/XINC
         TIME2=SECOND()
         TSTEP=TIME2-TIME1
         IF(TIMES)WRITE(6,'('' TIME FOR STEP:'',F8.2,'' LEFT'',F8.2)')
     1    TSTEP, TLEFT
         IF(TLAST-TLEFT.GT.TDUMP)THEN
            TLAST=TLEFT
            RESFIL=.TRUE.
            IPOW(9)=2
            I=ILOOP
            IPOW(8)=NSCF
            CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,I,BMAT,IPOW)
         ENDIF
         IF( TLEFT .LT. TSTEP*2.D0) THEN
C
C  STORE RESULTS TO DATE.
C
            IPOW(9)=1
            I=ILOOP
            IPOW(8)=NSCF
            CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,I,BMAT,IPOW)
            STOP
         ENDIF
   50 CONTINUE
C        *****  SCALE -HESSIAN- MATRIX                           *****
      IF( DEBUG) THEN
         WRITE(6,'(//10X,''UN-NORMALIZED HESSIAN MATRIX'')')
         DO 60 I=1,NVAR
   60    WRITE(6,'(8F10.4)')(HESS(J,I),J=1,NVAR)
      ENDIF
      DO 80 I=1,NVAR
         SUM = 0.0D0
         DO 70 J=1,NVAR
   70    SUM = SUM+HESS(I,J)**2
   80 WORK(I) = 1.0D0/SQRT(SUM)
      DO 90 I=1,NVAR
         DO 90 J=1,NVAR
   90 HESS(I,J) = HESS(I,J)*WORK(I)
      IF( DEBUG) THEN
         WRITE(6,'(//10X,''HESSIAN MATRIX'')')
         DO 100 I=1,NVAR
  100    WRITE(6,'(8F10.4)')(HESS(J,I),J=1,NVAR)
      ENDIF
C        *****  INITIALIZE B MATIRX                        *****
      DO 120 I=1,NVAR
         DO 110 J=1,NVAR
  110    BMAT(I,J) = 0.0D0
  120 BMAT(I,I) = WORK(I)*2.D0
************************************************************************
*
*  THIS IS THE START OF THE BIG LOOP TO OPTIMIZE THE GEOMETRY
*
************************************************************************
      ILOOP=-99
      TSTEP=TSTEP*4
  130 CONTINUE
      IF(TLAST-TLEFT.GT.TDUMP)THEN
         TLAST=TLEFT
         RESFIL=.TRUE.
         IPOW(9)=2
         I=ILOOP
         IPOW(8)=NSCF
         CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,I,BMAT,IPOW)
      ENDIF
      IF( TLEFT .LT. TSTEP*2.D0) THEN
C
C  STORE RESULTS TO DATE.
C
         IPOW(9)=1
         I=ILOOP
         IPOW(8)=NSCF
         CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,I,BMAT,IPOW)
         IFLEPO=-1
         RETURN
      ENDIF
  140 CONTINUE
C        *****  FORM-A- DAGGER-A- IN PA SLONG WITH -P-     *****
      IJ=0
      DO 160 J=1,NVAR
         DO 160 I=1,J
            IJ=IJ+1
            SUM = 0.0D0
            DO 150 K=1,NVAR
  150       SUM = SUM + HESS(I,K)*HESS(J,K)
  160 PMAT(IJ) = SUM
      DO 180 I=1,NVAR
         SUM = 0.0D0
         DO 170 K=1,NVAR
  170    SUM = SUM-HESS(I,K)*GMIN1(K)
  180 P(I) = -SUM
      L=0
      IF(DEBUG) THEN
         WRITE(6,'(/10X,''P MATRIX IN POWSQ'')')
         CALL VECPRT(PMAT,NVAR)
      ENDIF
      CALL RSP(PMAT,NVAR,NVAR,EIG,PVEC)
C        *****  CHECK FOR ZERO EIGENVALUE                  *****
C#      WRITE(6,'(''  EIGS IN POWSQ:'')')
C#      WRITE(6,'(6F13.8)')(EIG(I),I=1,NVAR)
      IF(EIG(1).LT.RHO2) GO TO 240
      INDC = 2
C        *****  IF MATRIX IS NOT SINGULAR FORM INVERSE     *****
C        *****  BY BACK TRANSFORMING THE EIGENVECTORS      *****
      IJ=0
      DO 200 I=1,NVAR
         DO 200 J=1,I
            IJ=IJ+1
            SUM = 0.0D0
            DO 190 K=1,NVAR
  190       SUM = SUM+PVEC((K-1)*NVAR+J)*PVEC((K-1)*NVAR+I)/EIG(K)
  200 PMAT(IJ) = SUM
C        *****  FIND -Q- VECTOR                            *****
      L=0
      IL=L+1
      L=IL+I-1
      DO 230 I=1,NVAR
         SUM = 0.0D0
         DO 210 K=1,I
            IK=(I*(I-1))/2+K
  210    SUM = SUM+PMAT(IK)*P(K)
         IP1=I+1
         DO 220 K=IP1,NVAR
            IK=(K*(K-1))/2+I
  220    SUM=SUM+PMAT(IK)*P(K)
  230 Q(I) = SUM
      GO TO 260
  240 CONTINUE
C        *****  TAKE  -Q- VECTOR AS EIGENVECTOR OF ZERO     *****
C        *****  EIGENVALUE                                 *****
      DO 250 I=1,NVAR
  250 Q(I) = PVEC(I)
  260 CONTINUE
C        *****  FIND SEARCH DIRECTION                      *****
      DO 270 I=1,NVAR
         SIG(I) = 0.0D0
         DO 270 J=1,NVAR
  270 SIG(I) = SIG(I) + Q(J)*BMAT(I,J)
C        *****  DO A ONE DIMENSIONAL SEARCH                *****
      IF (DEBUG) THEN
         WRITE(6,'('' SEARCH VECTOR'')')
         WRITE(6,'(8F10.5)')(SIG(I),I=1,NVAR)
      ENDIF
      CALL SEARCH(XPARAM, ALPHA, SIG, NVAR, GMIN, OKF, FUNCT)
      IF( NVAR .EQ. 1) GOTO 390
C
C  FIRST WE ATTEMPT TO OPTIMIZE GEOMETRY USING SEARCH.
C  IF THIS DOES NOT WORK, THEN SWITCH TO LINMIN, WHICH ALWAYS WORKS,
C  BUT IS TWICE AS SLOW AS SEARCH.
C
      ROUGH=  (   .NOT.  OKF)
      RMX = 0.0D0
      DO 280 K=1,NVAR
         RT = ABS(GMIN1(K))
         IF(RT.GT.RMX)RMX = RT
  280 CONTINUE
      IF(RMX.LT.TOL2) GO TO 390
C        *****  TWO STEP ESTIMATION OF DERIVATIVES         *****
      DO 290 K=1,NVAR
  290 E1(K) = (GMIN1(K)-GNEXT1(K))/(AMIN-ANEXT)
      RMU = DOT(E1,GMIN1,NVAR)/DOT(GMIN1,GMIN1,NVAR)
      DO 300 K=1,NVAR
  300 E2(K) = E1(K) - RMU*GMIN1(K)
C        *****  SCALE -E2- AND -SIG-                       *****
      SK = 1.0D0/SQRT(DOT(E2,E2,NVAR))
      DO 310 K=1,NVAR
  310 SIG(K) = SK*SIG(K)
      DO 320 K=1,NVAR
  320 E2(K) = SK*E2(K)
C        *****  FIND INDEX OF REPLACEMENT DIRECTION        *****
      PMAX = -1.0D+20
      DO 330 I=1,NVAR
         IF(ABS(P(I)*Q(I)).LE.PMAX) GO TO 330
         PMAX = ABS(P(I)*Q(I))
         ID = I
  330 CONTINUE
C        *****  REPLACE APPROPRIATE DIRECTION AND DERIVATIVE ***
      DO 340 K=1,NVAR
  340 HESS(ID,K) = -E2(K)
C        *****  REPLACE STARTING POINT                     *****
      DO 350 K=1,NVAR
  350 BMAT(K,ID) = SIG(K)/0.529167D0
      DO 360 K=1,NVAR
  360 GNEXT1(K) = GMIN1(K)
      GLAST = GMIN
      INDC = 1
      TIME1=TIME2
      TIME2=SECOND()
      TLEFT=TLEFT-TIME2+TIME0
      TSTEP=TIME2-TIME1
      ICYC=ICYC+1
      IF(RESFIL)THEN
         WRITE(6,370)MIN(TLEFT,9999999.9D0),
     1MIN(GMIN,999999.999D0),FUNCT
         IF(LOG)WRITE(11,370)MIN(TLEFT,9999999.9D0),
     1MIN(GMIN,999999.999D0),FUNCT
  370    FORMAT('  RESTART FILE WRITTEN,  TIME LEFT:',F9.1,
     1' GRAD.:',F10.3,' HEAT:',G14.7)
         RESFIL=.FALSE.
      ELSE
         WRITE(6,380)ICYC,MIN(TSTEP,9999.99D0),
     1MIN(TLEFT,9999999.9D0),MIN(GMIN,999999.999D0),FUNCT
         IF(LOG)WRITE(11,380)ICYC,MIN(TSTEP,9999.99D0),
     1MIN(TLEFT,9999999.9D0),MIN(GMIN,999999.999D0),FUNCT
  380    FORMAT(' CYCLE:',I5,' TIME:',F6.1,' TIME LEFT:',F9.1,
     1' GRAD.:',F10.3,' HEAT:',G14.7)
      ENDIF
      IF(TIMES)WRITE(6,'('' TIME FOR STEP:'',F8.2,'' LEFT'',F8.2)')
     1TSTEP, TLEFT
      GO TO 130
  390 CONTINUE
      DO 400 I=1,NVAR
  400 GRAD(I)=0.D0
      LAST=1
      CALL COMPFG(XPARAM, .TRUE., FUNCT, .TRUE., GRAD, .TRUE.)
      DO 410 I=1,NVAR
  410 GRAD(I)=GMIN1(I)
      GNFINA=SQRT(DOT(GRAD,GRAD,NVAR))
      IFLEPO=11
      IF(SCF1)IFLEPO=13
      RETURN
      END
