SUBROUTINE POWELL
C
C  This is the Powell-method [1] conjugate gradient minimizer.
C
C  For syntax see main parsing loop "HELP"
C
C  Parts of the code were was extracted from the IMSL [2] library.
C
C  Author: Axel T. Brunger
C  =======================
C
C  References:
C  [1] M.J.D. Powell, "Restart Procedures for the Conjugate Gradient
C      Method", Mathematical Programming 12 (1977), pp. 241 - 254
C  [2] IMSL Inc. library manual (1982), Houston, Texas.
C
C Modified by Yajun Liu to contral output
      IMPLICIT NONE
C input/output
      INCLUDE 'cns.inc'
      INCLUDE 'comand.inc'
      INCLUDE 'heap.inc'
      INCLUDE 'mtf.inc'
      INCLUDE 'funct.inc'
C local
      DOUBLE PRECISION  DFPRED, ACC, TOLG
      INTEGER XREF, YREF, ZREF, XCURR, GCURR, W, WXOPT, WGOPT, WGINIT
      INTEGER WRSDG, WRSDX, IER, NDIM, I, MAXFN
      LOGICAL QDEBUG
C YL
      INTEGER NPRINT
C begin
C
C defaults
      MAXFN=500
      DFPRED=0.001D0
      TOLG=0.0001D0
      QDEBUG=.FALSE.
C YL
      NPRINT = 1
C
C parsing
      CALL PUSEND('POWELL>')
      DO WHILE (.NOT.DONE)
      CALL NEXTWD('POWELL>')
C
      IF (WD(1:4).EQ.'HELP') THEN
C
      CALL CNSHELP('cns-minimize-powell')
C
      ELSE IF (WD(1:4).EQ.'NSTE') THEN
      CALL NEXTI('NSTEp=',MAXFN)
      ELSE IF (WD(1:4).EQ.'STEP'.OR.WD(1:4).EQ.'DROP') THEN
      CALL NEXTF('DROP=',DFPRED)
      ELSE IF (WD(1:4).EQ.'TOLG') THEN
      CALL NEXTF('TOLGradient=',TOLG)
      ELSE IF (WD(1:4).EQ.'DEBU') THEN
      CALL NEXTLO('DEBUg=',QDEBUG)
C YL
      ELSE IF (WD(1:4).EQ.'NPRI') THEN
      CALL NEXTI('NPRInt=',NPRINT)
      ELSE IF (WD(1:4).EQ.'?   ') THEN
      WRITE(6,1200) MAXFN, DFPRED, TOLG
1200  FORMAT(/' POWELL-method :  NSTEp=',I7,
     1       '  DROP=',F12.9,'  TOLGradient=',F12.9,/)
      ELSE
      CALL CHKEND('POWELL>',DONE)
      END IF
      END DO
      DONE=.FALSE.
C
C NDIM is the dimension of the problem
      NDIM=0
      DO I=1,NATOM
      IF (IMOVE(I).EQ.0) NDIM=NDIM+3
      END DO
      IF (NDIM.EQ.0) THEN
      WRITE(6,'(A)') ' %POWELL-ERR: all atoms fixed. No action'
      ELSE
C
      WRITE(6,'(A,I6)')
     &  ' POWELL: number of degrees of freedom=',NDIM
C
      ACC=TOLG**2 *NDIM
C
C dynamic
      XREF=ALLHP(IREAL8(NATOM))
      YREF=ALLHP(IREAL8(NATOM))
      ZREF=ALLHP(IREAL8(NATOM))
      XCURR=ALLHP(IREAL8(NDIM))
      GCURR=ALLHP(IREAL8(NDIM))
      W=ALLHP(IREAL8(NDIM))
      WXOPT=ALLHP(IREAL8(NDIM))
      WGOPT=ALLHP(IREAL8(NDIM))
      WGINIT=ALLHP(IREAL8(NDIM))
      WRSDG=ALLHP(IREAL8(NDIM))
      WRSDX=ALLHP(IREAL8(NDIM))
C
      CALL POWEL2(NDIM,ACC,MAXFN,DFPRED,HEAP(XREF),HEAP(YREF),
     &           HEAP(ZREF),HEAP(XCURR),HEAP(GCURR),HEAP(W),
     &           HEAP(WXOPT),HEAP(WGOPT),HEAP(WGINIT),HEAP(WRSDG),
     &           HEAP(WRSDX),IER,QDEBUG,NPRINT)
C YL
C
      IF (IER.EQ.0) THEN
      WRITE(6,'(A)') ' POWELL: Gradient converged. Normal termination'
      ELSE IF (IER.EQ.129) THEN
      WRITE(6,'(A)')' POWELL: Line search terminated'
      ELSE IF (IER.EQ.130) THEN
      WRITE(6,'(A)')' %POWELL-ERR: Search direction uphill'
      ELSE IF (IER.EQ.131) THEN
      WRITE(6,'(A)') ' POWELL: STEP number limit. Normal termination'
      ELSE IF (IER.EQ.132) THEN
      WRITE(6,'(A)')
     & ' %POWELL-ERR: Two consecutive iterations failed to reduce E'
      END IF
      WRITE(6,'(A)') ' POWELL: Current coordinates set to last minimum'
C
      CALL FREHP(XCURR,IREAL8(NDIM))
      CALL FREHP(GCURR,IREAL8(NDIM))
      CALL FREHP(W,IREAL8(NDIM))
      CALL FREHP(WXOPT,IREAL8(NDIM))
      CALL FREHP(WGOPT,IREAL8(NDIM))
      CALL FREHP(WGINIT,IREAL8(NDIM))
      CALL FREHP(WRSDG,IREAL8(NDIM))
      CALL FREHP(WRSDX,IREAL8(NDIM))
      CALL FREHP(ZREF,IREAL8(NATOM))
      CALL FREHP(YREF,IREAL8(NATOM))
      CALL FREHP(XREF,IREAL8(NATOM))
C
      END IF
      RETURN
      END
C
      SUBROUTINE POWEL2(N,ACC,MAXFN,DFPRED,XREF,YREF,ZREF,
     &                 XCURR,GCURR,W,WXOPT,WGOPT,WGINIT,WRSDG,
     &                 WRSDX,IER,QDEBUG,NPRINT)
C YL
C
C See comments in routine POWELL. Original comments were compressed.
C Modifications are added as FLECS sub-procedures. The work array was
C split into 6 individual arrays. Debugging information was added.
C
C DFPRED - A ROUGH ESTIMATE OF THE EXPECTED REDUCTION IN F, WHICH IS
C USED TO DETERMINE THE SIZE OF THE INITIAL CHANGE TO X.
C NOTE THAT DFPRED IS THE EXPECTED REDUCTION ITSELF, AND DOES NOT
C DEPEND ON ANY RATIOS. A BAD VALUE OF DFPRED CAUSES AN ERROR
C MESSAGE, WITH IER=129, AND A RETURN ON THE FIRST ITERATION.
C IER    - ERROR PARAMETER. IER = 0 IMPLIES THAT CONVERGENCE
C WAS ACHIEVED AND NO ERRORS OCCURRED. TERMINAL ERROR IER = 129 IMPLIES
C THAT THE LINE SEARCH OF AN INTEGRATION WAS ABANDONED. THIS ERROR MAY BE
C CAUSED BY AN ERROR IN THE GRADIENT. IER = 130 IMPLIES THAT THE
C CALCULATION CANNOT CONTINUE BECAUSE THE SEARCH DIRECTION IS UPHILL.
C IER = 131 IMPLIES THAT THE ITERATION WAS TERMINATED BECAUSE MAXFN WAS
C EXCEEDED. IER = 132 IMPLIES THAT THE CALCULATION WAS TERMINATED
C BECAUSE TWO CONSECUTIVE ITERATIONS FAILED TO REDUCE F.
C
C REMARKS:
C 1.  THE ROUTINE INCLUDES NO THOROUGH CHECKS ON THE PART OF
C THE USER PROGRAM THAT CALCULATES THE DERIVATIVES OF THE OBJECTIVE
C FUNCTION. THEREFORE, BECAUSE DERIVATIVE CALCULATION IS A FREQUENT
C SOURCE OF ERROR, THE USER SHOULD VERIFY INDEPENDENTLY THE CORRECTNESS
C OF THE DERIVATIVES THAT ARE GIVEN TO THE ROUTINE.
C 2.  BECAUSE OF THE CLOSE RELATION BETWEEN THE CONJUGATE GRADIENT
C METHOD AND THE METHOD OF STEEPEST DESCENTS, IT IS VERY HELPFUL TO
C CHOOSE THE SCALE OF THE VARIABLES IN A WAY THAT BALANCES THE
C MAGNITUDES OF THE COMPONENTS OF A TYPICAL DERIVATE VECTOR. IT CAN BE
C PARTICULARLY INEFFICIENT IF A FEW COMPONENTS OF THE GRADIENT ARE MUCH
C LARGER THAN THE REST.
C 3.  IF THE VALUE OF THE PARAMETER ACC IN THE ARGUMENT LIST OF THE
C ROUTINE IS SET TO ZERO, THEN THE SUBROUTINE WILL CONTINUE ITS
C CALCULATION UNTIL IT STOPS REDUCING THE OBJECTIVE FUNCTION. IN THIS
C CASE THE USUAL BEHAVIOUR IS THAT CHANGES IN THE OBJECTIVE FUNCTION
C BECOME DOMINATED BY COMPUTER ROUNDING ERRORS BEFORE PRECISION IS LOST
C IN THE GRADIENT VECTOR. THEREFORE, BECAUSE THE POINT OF VIEW HAS BEEN
C TAKEN THAT THE USER REQUIRES THE LEAST POSSIBLE VALUE OF THE
C FUNCTION, A VALUE OF THE OBJECTIVE FUNCTION THAT IS SMALL DUE TO
C COMPUTER ROUNDING ERRORS CAN PREVENT FURTHER PROGRESS. HENCE THE
C PRECISION IN THE FINAL VALUES OF THE VARIABLES MAY BE ONLY ABOUT HALF
C THE NUMBER OF SIGNIFICANT DIGITS IN THE COMPUTER ARITHMETIC, BUT THE
C LEAST VALUE OF F IS USUALLY FOUND TO QUITE HIGH ACCURACY.
C
      IMPLICIT NONE
C
      INCLUDE 'cns.inc'
      INCLUDE 'coord.inc'
      INCLUDE 'mtf.inc'
      INCLUDE 'ener.inc'
      INCLUDE 'deriv.inc'
C
      INTEGER  N
      DOUBLE PRECISION   XREF(*), YREF(*), ZREF(*)
      INTEGER  MAXFN, IER
      DOUBLE PRECISION   ACC, DFPRED, XCURR(*), GCURR(*)
      DOUBLE PRECISION   W(*), WXOPT(*), WGOPT(*), WGINIT(*)
      DOUBLE PRECISION   WRSDG(*), WRSDX(*)
      LOGICAL  QDEBUG
C YL
      INTEGER NPRINT
C
C
      INTEGER  MAXLIN, MXFCON, I, IRETRY, IAT3
      INTEGER  ITERC, ITERFM, ITERRS, NCALLS, NFBEG, NFOPT
      DOUBLE PRECISION   BETA, DDSPLN, DFPR, FCH, FINIT, FMIN, GAMDEN
      DOUBLE PRECISION   GAMA, GINIT, GMIN, GNEW, GSPLN, GSQRD, SBOUND
      DOUBLE PRECISION   STEP, STEPCH, STMIN, SUM, WORK, F, DBPREC
      DOUBLE COMPLEX DBCOMP
C
      PARAMETER (MAXLIN=5,MXFCON=2)
C
      IER = 0
C
C W(*) IS USED FOR THE SEARCH DIRECTION OF AN ITERATION. WRSDX(*) AND
C WRSDG(*) CONTAIN THE INFORMATION THAT IS REQUIRED BY THE CONJUGACY
C CONDITIONS OF THE RESTART PROCEDURE. WGINIT(*) CONTAINS THE GRADIENT
C AT THE START OF AN ITERATION. WXOPT(*)==(XREF,YREF,ZREF) CONTAINS THE
C PARAMETERS THAT GIVE THE LEAST CALCULATED VALUE OF F. WGOPT(*)
C CONTAINS THE GRADIENT VECTOR WHERE F IS LEAST.
C
C SET SOME PARAMETERS TO BEGIN THE CALCULATION. ITERC AND NCALLS COUNT
C THE NUMBER OF ITERATIONS AND CALLS OF ENERGY. ITERFM IS THE NUMBER OF
C THE MOST RECENT ITERATION THAT DECREASES F.
C
C ===========================
C insertion: initialize XCURR
      IAT3=1
      DO I=1,NATOM
      IF (IMOVE(I).EQ.0) THEN
      XCURR(IAT3)=X(I)
      XCURR(IAT3+1)=Y(I)
      XCURR(IAT3+2)=Z(I)
      IAT3=IAT3+3
      END IF
      END DO
C end insertion
C ============================
C ======================================
C insertion: initialize XREF, YREF, ZREF
      DO I=1,NATOM
      XREF(I)=X(I)
      YREF(I)=Y(I)
      ZREF(I)=Z(I)
      END DO
C end insertion
C ======================================
C
      ITERC = 0
      NCALLS = 0
      ITERFM = ITERC
      STEPCH=0.0D0
C
C CALL ENERGY. LET THE INITIAL SEARCH DIRECTION BE MINUS THE GRADIENT
C VECTOR. USUALLY THE PARAMETER ITERRS GIVES THE ITERATION NUMBER OF
C THE MOST RECENT RESTART, BUT IT IS SET TO ZERO WHEN THE STEEPEST
C DESCENT DIRECTION IS USED.
C
    5 NCALLS = NCALLS+1
      DBPREC = NCALLS
      CALL DECLAR( 'MINI_CYCLES', 'DP', ' ', DBCOMP, DBPREC )
      IAT3=1
C =======================================================
C insertion: get energy and gradient gcurr at point xcurr
      DO I=1,NATOM
      IF (IMOVE(I).EQ.0) THEN
      X(I)=XCURR(IAT3)
      Y(I)=XCURR(IAT3+1)
      Z(I)=XCURR(IAT3+2)
      IAT3=IAT3+3
      END IF
      END DO
      CALL ENERGY
C zero forces of fixed atoms
      DO I=1,NATOM
      IF (IMOVE(I).NE.0) THEN
      DX(I)=0.0D0
      DY(I)=0.0D0
      DZ(I)=0.0D0
      END IF
      END DO
      F=RENR(SSENER)
      IAT3=1
      DO I=1,NATOM
      IF (IMOVE(I).EQ.0) THEN
      XCURR(IAT3)=X(I)
      XCURR(IAT3+1)=Y(I)
      XCURR(IAT3+2)=Z(I)
      GCURR(IAT3)=DX(I)
      GCURR(IAT3+1)=DY(I)
      GCURR(IAT3+2)=DZ(I)
      IAT3=IAT3+3
      END IF
      END DO
      IF (QDEBUG) WRITE(6,1800) STEPCH
1800  FORMAT(' POWELL: STEPCH=',G17.10)
C YL
      IF(NPRINT.GT.0) THEN
        IF(MOD(NCALLS,NPRINT).EQ.0) THEN
          WRITE(6,'(A,I6,A,F10.4,A)')
     & ' --------------- cycle=',NCALLS,
     & ' ------ stepsize=',STEPCH,' -----------------------'
      CALL PRINTE(RENR)
      WRITE(6,'(2A)') ' --------------------------------------------',
     & '-----------------------------------'
        ENDIF
      ENDIF
C end insertion
C =========================================================
C
      IF (NCALLS.GE.2) GO TO 20
   10 CONTINUE
      DO I=1,N
      W(I) = -GCURR(I)
      END DO
      ITERRS = 0
      IF (ITERC.GT.0) GO TO 80
C
C SET SUM TO G SQUARED. GMIN AND GNEW ARE THE OLD AND THE NEW
C DIRECTIONAL DERIVATIVES ALONG THE CURRENT SEARCH DIRECTION. LET FCH
C BE THE DIFFERENCE BETWEEN F AND THE PREVIOUS BEST VALUE OF THE
C OBJECTIVE FUNCTION.
C
   20 GNEW = 0.0D0
      SUM = 0.0D0
      DO I=1,N
      GNEW = GNEW+W(I)*GCURR(I)
      SUM = SUM+GCURR(I)**2
      END DO
      IF (QDEBUG) WRITE(6,1050) GNEW,F
1050  FORMAT(' POWELL: New search point. GNEW=',G17.10,' Fnew=',
     1       G17.10)
      IF (NCALLS.EQ.1) GO TO 35
      FCH = F-FMIN
C
C STORE THE VALUES OF XCURR, F AND GCURR, IF  THEY ARE THE BEST THAT
C HAVE BEEN CALCULATED SO FAR, AND NOTE GCURR  SQUARED AND THE VALUE OF
C NCALLS. TEST FOR CONVERGENCE.
C
      IF (FCH) 35,30,50
   30 IF (GNEW/GMIN.LT.-1.0D0) GO TO 45
   35 FMIN = F
      GSQRD = SUM
      NFOPT = NCALLS
      IF (QDEBUG) WRITE(6,1030)
1030  FORMAT(' POWELL: least energy point set to current point.')
      DO I=1,N
      WXOPT(I) = XCURR(I)
      WGOPT(I) = GCURR(I)
      END DO
C insertion: copy wxopt to xref, yref, zref
C ===========================================
      IAT3=1
      DO I=1,NATOM
      IF (IMOVE(I).EQ.0) THEN
      XREF(I)=WXOPT(IAT3)
      YREF(I)=WXOPT(IAT3+1)
      ZREF(I)=WXOPT(IAT3+2)
      IAT3=IAT3+3
      ELSE
      XREF(I)=X(I)
      YREF(I)=Y(I)
      ZREF(I)=Z(I)
      END IF
      END DO
C ===========================================
C end insertion
   45 IF (SUM.LE.ACC) GO TO 9005
C
C CHECK IF THE VALUE OF MAXFN ALLOWS ANOTHER CALL OF ENERGY.
C
   50 IF (NCALLS.NE.MAXFN) GO TO 55
      IER = 131
      GO TO 9000
   55 IF (NCALLS.GT.1) GO TO 100
C
C SET DFPR TO THE ESTIMATE OF THE  REDUCTION IN F GIVEN IN THE
C ARGUMENT LIST, IN ORDER THAT THE  INITIAL CHANGE TO THE PARAMETERS
C IS OF A SUITABLE SIZE. THE VALUE  OF STMIN IS USUALLY THE
C STEP-LENGTH OF THE MOST RECENT  LINE SEARCH THAT GIVES THE LEAST
C CALCULATED VALUE OF F.
C
CCC*** modification.  Original two lines commented out.
CCC*** modification by Jiang Jiansheng, York, UK.  12/1/90, ATB.
CCC**      DFPR = DFPRED
CCC**      STMIN = DFPRED/GSQRD
      DFPR = MAX(DFPRED, SQRT(GSQRD/N))
      STMIN = DFPR/GSQRD
CCC***
C
C BEGIN THE ITERATION
C
   80 ITERC = ITERC+1
C
C STORE THE INITIAL ENERGY VALUE AND GRADIENT, CALCULATE THE INITIAL
C DIRECTIONAL DERIVATIVE, AND BRANCH IF ITS VALUE IS NOT NEGATIVE. SET
C SBOUND TO MINUS ONE TO INDICATE THAT A BOUND ON THE STEP IS NOT KNOWN
C YET, AND SET NFBEG TO THE CURRENT VALUE OF NCALLS. THE PARAMETER
C IRETRY SHOWS THE NUMBER OF ATTEMPTS AT SATISFYING THE BETA CONDITION.
C
      FINIT = F
      GINIT = 0.0D0
      DO I=1,N
      WGINIT(I) = GCURR(I)
      GINIT = GINIT+W(I)*GCURR(I)
      END DO
      IF (GINIT.GE.0.0D0) GO TO 165
      GMIN = GINIT
      SBOUND = -1.0D0
      NFBEG = NCALLS
      IRETRY = -1
C
C SET STEPCH SO THAT THE INITIAL STEP-LENGTH IS CONSISTENT WITH THE
C PREDICTED REDUCTION IN F, SUBJECT TO THE CONDITION THAT IT DOES NOT
C EXCEED THE STEP-LENGTH OF THE PREVIOUS ITERATION. LET STMIN BE THE
C STEP TO THE LEAST CALCULATED VALUE OF F.
C
      STEPCH = MIN(STMIN,ABS(DFPR/GINIT))
      STMIN = 0.0D0
C
C CALL ENERGY AT THE VALUE OF X THAT IS DEFINED BY THE NEW CHANGE TO
C THE STEP-LENGTH, AND LET THE NEW STEP-LENGTH BE STEP. THE VARIABLE
C WORK IS USED AS WORK SPACE.
C
   90 STEP = STMIN+STEPCH
      WORK = 0.0D0
      DO I=1,N
      XCURR(I) = WXOPT(I)+STEPCH*W(I)
      WORK = MAX(WORK,ABS(XCURR(I)-WXOPT(I)))
      END DO
      IF (WORK.GT.0.0D0) GO TO 5
C
C TERMINATE THE LINE SEARCH IF STEPCH IS EFFECTIVELY ZERO.
C
      IF (NCALLS.GT.NFBEG+1) GO TO 115
      IF (ABS(GMIN/GINIT)-0.2D0) 170,170,115
C
C LET SPLN BE THE QUADRATIC SPLINE THAT INTERPOLATES THE CALCULATED
C ENERGY VALUES AND DIRECTIONAL DERIVATIVES AT THE POINTS STMIN AND
C STEP OF THE LINE SEARCH, WHERE THE KNOT OF THE SPLINE IS AT
C 0.5*(STMIN+STEP). REVISE STMIN, GMIN AND SBOUND, AND SET DDSPLN TO
C THE SECOND DERIVATIVE OF SPLN AT THE NEW STMIN. HOWEVER, IF FCH IS
C ZERO, IT IS ASSUMED THAT THE MAXIMUM ACCURACY IS ALMOST ACHIEVED, SO
C DDSPLN IS CALCULATED USING ONLY THE CHANGE IN THE GRADIENT.
C
  100 WORK = (FCH+FCH)/STEPCH-GNEW-GMIN
      DDSPLN = (GNEW-GMIN)/STEPCH
      IF (NCALLS.GT.NFOPT) SBOUND = STEP
      IF (NCALLS.GT.NFOPT) GO TO 105
      IF (GMIN*GNEW.LE.0.0D0) SBOUND = STMIN
      STMIN = STEP
      GMIN = GNEW
      STEPCH = -STEPCH
  105 IF (FCH.NE.0.0D0) DDSPLN = DDSPLN+(WORK+WORK)/STEPCH
C
C TEST FOR CONVERGENCE OF THE LINE SEARCH, BUT FORCE AT LEAST TWO STEPS
C TO BE TAKEN IN ORDER NOT TO LOSE QUADRATIC TERMINATION.
C
      IF (GMIN.EQ.0.0D0) GO TO 170
      IF (NCALLS.LE.NFBEG+1) GO TO 120
      IF (ABS(GMIN/GINIT).LE.0.2D0) GO TO 170
C
C APPLY THE TEST THAT DEPENDS ON THE PARAMETER MAXLIN.
C
  110 IF (NCALLS.LT.NFOPT+MAXLIN) GO TO 120
  115 IER = 129
      GO TO 170
C
C SET STEPCH TO THE GREATEST CHANGE TO THE CURRENT VALUE OF STMIN THAT
C IS ALLOWED BY THE BOUND ON THE LINE SEARCH. SET GSPLN TO THE GRADIENT
C OF THE QUADRATIC SPLINE AT (STMIN+STEPCH). HENCE CALCULATE THE VALUE
C OF STEPCH THAT MINIMIZES THE SPLINE FUNCTION, AND THEN OBTAIN THE NEW
C FUNCTION AND GRADIENT VECTOR, FOR THIS VALUE OF THE CHANGE TO THE
C STEP-LENGTH.
C
  120 STEPCH = 0.5D0*(SBOUND-STMIN)
      IF (SBOUND.LT.-0.5D0) STEPCH = 9.0D0*STMIN
      GSPLN = GMIN+STEPCH*DDSPLN
      IF (GMIN*GSPLN.LT.0.0D0) STEPCH = STEPCH*GMIN/(GMIN-GSPLN)
      GO TO 90
C
C CALCULATE THE VALUE OF BETA THAT OCCURS IN THE NEW SEARCH
C DIRECTION.
C
  125 SUM = 0.0D0
      DO I=1,N
      SUM = SUM+GCURR(I)*WGINIT(I)
      END DO
      BETA = (GSQRD-SUM)/(GMIN-GINIT)
C
C TEST THAT THE NEW SEARCH DIRECTION CAN BE MADE DOWNHILL. IF IT
C CANNOT, THEN MAKE ONE ATTEMPT TO IMPROVE THE ACCURACY OF THE LINE
C SEARCH.
C
      IF (ABS(BETA*GMIN).LE.0.2D0*GSQRD) GO TO 135
      IRETRY = IRETRY+1
      IF (IRETRY.LE.0) GO TO 110
C
C APPLY THE TEST THAT DEPENDS ON THE PARAMETER MXFCON. SET DFPR TO THE
C PREDICTED REDUCTION IN F ON THE NEXT ITERATION.
C
  135 IF (F.LT.FINIT) ITERFM = ITERC
      IF (ITERC.LT.ITERFM+MXFCON) GO TO 140
      IER = 132
      GO TO 9000
  140 DFPR = STMIN*GINIT
C
C BRANCH IF A RESTART PROCEDURE IS REQUIRED DUE TO THE ITERATION NUMBER
C OR DUE TO THE SCALAR PRODUCT OF CONSECUTIVE GRADIENTS.
C
      IF (IRETRY.GT.0) GO TO 10
      IF (ITERRS.EQ.0) GO TO 155
      IF (ITERC-ITERRS.GE.N) GO TO 155
      IF (ABS(SUM).GE.0.2D0*GSQRD) GO TO 155
C
C CALCULATE THE VALUE OF GAMA THAT OCCURS IN THE NEW SEARCH DIRECTION,
C AND SET SUM TO A SCALAR PRODUCT FOR THE TEST BELOW. THE VALUE OF
C GAMDEN IS SET BY THE RESTART PROCEDURE.
C
      GAMA = 0.0D0
      SUM = 0.0D0
      DO I=1,N
      GAMA = GAMA+GCURR(I)*WRSDG(I)
      SUM = SUM+GCURR(I)*WRSDX(I)
      END DO
      GAMA = GAMA/GAMDEN
C
C RESTART IF THE NEW SEARCH DIRECTION IS NOT SUFFICIENTLY DOWNHILL.
C
      IF (ABS(BETA*GMIN+GAMA*SUM).GE.0.2D0*GSQRD) GO TO 155
C
C CALCULATE THE NEW SEARCH DIRECTION.
C
      IF (QDEBUG) WRITE(6,1190) BETA,GAMA
1190  FORMAT(' POWELL: New search direction. BETA=',G17.10,' GAMA=',
     1       G17.10)
      DO I=1,N
      W(I) = -GCURR(I)+BETA*W(I)+GAMA*WRSDX(I)
      END DO
C
C ==================================================
C insertion: project w at point xcurr
      IAT3=1
      DO I=1,NATOM
      IF (IMOVE(I).EQ.0) THEN
      DX(I)=W(IAT3)
      DY(I)=W(IAT3+1)
      DZ(I)=W(IAT3+2)
      IAT3=IAT3+3
      ELSE
      DX(I)=0.0D0
      DY(I)=0.0D0
      DZ(I)=0.0D0
      END IF
      END DO
      IAT3=1
      DO I=1,NATOM
      IF (IMOVE(I).EQ.0) THEN
      W(IAT3)=DX(I)
      W(IAT3+1)=DY(I)
      W(IAT3+2)=DZ(I)
      IAT3=IAT3+3
      END IF
      END DO
C end insertion
C ==================================================
C
      GO TO 80
C
C APPLY THE RESTART PROCEDURE.
C
  155 GAMDEN = GMIN-GINIT
      IF (QDEBUG) WRITE(6,1200) BETA
1200  FORMAT(' POWELL: Applying restart procedure. BETA=',G17.10)
      DO I=1,N
      WRSDX(I) = W(I)
      WRSDG(I) = GCURR(I)-WGINIT(I)
      W(I) = -GCURR(I)+BETA*W(I)
      END DO
C
C =========================================================
C insertion: project w at point xcurr (note: same as above)
      IAT3=1
      DO I=1,NATOM
      IF (IMOVE(I).EQ.0) THEN
      DX(I)=W(IAT3)
      DY(I)=W(IAT3+1)
      DZ(I)=W(IAT3+2)
      IAT3=IAT3+3
      ELSE
      DX(I)=0.0D0
      DY(I)=0.0D0
      DZ(I)=0.0D0
      END IF
      END DO
      IAT3=1
      DO I=1,NATOM
      IF (IMOVE(I).EQ.0) THEN
      W(IAT3)=DX(I)
      W(IAT3+1)=DY(I)
      W(IAT3+2)=DZ(I)
      IAT3=IAT3+3
      END IF
      END DO
C end insertion
C ==================================================
      ITERRS = ITERC
      GO TO 80
C
C SET IER TO INDICATE THAT THE SEARCH DIRECTION IS UPHILL.
C
  165 IER = 130
C
C ENSURE THAT F, X AND G ARE OPTIMAL.
C
  170 IF (NCALLS.EQ.NFOPT) GO TO 180
      IF (QDEBUG) WRITE(6,1300)
1300  FORMAT(' POWELL: Current point set to optimal point.')
      F = FMIN
      DO I=1,N
      XCURR(I) = WXOPT(I)
      GCURR(I) = WGOPT(I)
      END DO
  180 IF (IER.EQ.0) GO TO 125
 9000 CONTINUE
 9005 CONTINUE
C
C================================================
C-insertion: copy least energy point to x, y, z
      IAT3=1
      DO I=1,NATOM
      IF (IMOVE(I).EQ.0) THEN
      X(I)=WXOPT(IAT3)
      Y(I)=WXOPT(IAT3+1)
      Z(I)=WXOPT(IAT3+2)
      DX(I)=WGOPT(IAT3)
      DY(I)=WGOPT(IAT3+1)
      DZ(I)=WGOPT(IAT3+2)
      IAT3=IAT3+3
      ELSE
      DX(I)=0.0D0
      DY(I)=0.0D0
      DZ(I)=0.0D0
      END IF
      END DO
C end insertion
C ================================================
C
      RETURN
      END