C 31 JAN 07 - DGF - TRUNG: DO DK INTEGRALS IF REQUIRED
C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
C  8 MAY 06 - BSD - TRUHAT: CHANGE A PRINT OPTION
C 17 JAN 06 - MWS - CHANGE A FEW FORMATS
C 19 SEP 05 - MWS - ADD TRUE NUCLEAR CHARGE ARRAY TO INFOA COMMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C  1 JUN 05 - BSD - ADDED CONSTRAINTS FOR VB2000 USE, ALLOW 20 PARAMS
C 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL
C  7 APR 04 - BSD - TRUNRG: CORRECT FORMAT STATEMENT
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C  3 SEP 03 - MWS - SYNCH UP RELWFN COMMON BLOCK
C  6 SEP 01 - MWS - ADD DUMMY ARGUMENTS TO NAMEIO CALL
C  1 AUG 01 - DGF - REDO THE UNCONTRACTED BASIS SET FOR OPTMIZ=BASIS
C 12 NOV 98 - GDF - MAKE SYMTRY COMMON CONSISTENT
C 28 SEP 97 - KRG - TRUNRG: FIX CHECK JOBS
C 19 FEB 97 - MWS - TRUNRG: STOP IF ENERGY DOESN'T CONVERGE
C 22 JUN 96 - VAG - CHANGES FOR ADDITION OF CITYP VARIABLE
C  1 APR 96 - MWS - TRUSRC: SAVE ADDED TO TRUSRC TO MAKE AXP HAPPY
C 25 MAR 96 - MWS - TRUNRG: ERASE PREVIOUS CI EIGENVECTOR FILES
C  5 MAR 96 - MWS - CHANGE SHELL SYMMETRY PACKING COMMONS
C 11 NOV 94 - MWS - REMOVE FTNCHEK WARNINGS
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C  1 JUN 94 - MWS - TRUNRG: DON'T ORDER THE AO INTEGRALS FOR CI TRUDGE
C 16 JUL 93 - MWS - INCREASE MAXIMUM CI ROOTS TO 100
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C 10 JAN 92 - TLW - TRUDGX,TRUINP,TRURST: CALL ABRT BEFORE STOP
C 10 JAN 92 - MWS,TLW - CHANGE OPENPK TO SEQOPN
C  8 JAN 92 - TLW - MAKE WRITES PARALLEL; ADD COMMON PAR
C  1 MAR 91 - MWS - ENSURE THAT THIS CODE WORKS IN ANGSTROMS.
C 11 FEB 91 - MWS - TRUDGE CALLS ABRT INSTEAD OF EXIT
C 31 OCT 90 - STE - CLEAN UP TRUSRC SOME
C 31 AUG 90 - MWS - FIX PRINTOUT FOR EXPONENT OPTIMIZATIONS
C  7 AUG 90 - TLW - ADD CF AND CG TO COMMON NSHEL
C  2 JUN 90 - MWS - SOME SMALL DEGREE OF CIVILIZATION
C 12 APR 90 - MK  - CHANGES TO ALLOW FOR NONGRADIENT CI OPTIMIZATION
C  2 APR 90 - MWS - TRUNRG: SAVE STATEMENT, TRUSRC: USE TEN35
C 20 MAR 90 - MK  - NON-GRADIENT OPTIMIZATION OF GEOMETRY --
C                   ANOTHER GEM FROM HONDO 7.0 ...
C*MODULE TRUDGE  *DECK TRUDGX
      SUBROUTINE TRUDGX
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL RSTART,EXACT,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (NNAM=10)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /TRUDIR/ IROOT,JSTART,EXACT
      COMMON /TRUOPT/ OPTMIZ,F,FNOT,TOLF,R(20),RNOT(20),TOLR,
     *                P(20),P0(20),IEX(20),NPAR
      COMMON /TRUPAR/ DR(20),ALPH(20),CURV(20),V(20,20),ALPHA,CURVE,
     *                FNOISE,NDIR,NDIRMX
      COMMON /TRURUN/ TIME,RSTART
C
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
C
      PARAMETER (ZERO=0.0D+00,SMALL=1.0D-12)
      PARAMETER (RTOL=5.0D-02, FTOL=1.0D-03, TOL=5.0D-04)
C
      DATA RNAM/8HTRURST  /
      DATA GEOMET/8HGEOMETRY/
      DATA QNAM/8HKSTART  ,8HJSTART  ,8HTOLF    ,8HTOLR    ,
     *          8HFNOISE  ,8HFNOT    ,8HRNOT    ,8HCURV    ,
     *          8HALPH    ,8HV       /
      DATA KQNAM/1,1,3,3,  3,3,1003,1003,  1003,10003/
C
      KSTART=-1
      JSTART=0
      TOLF=FTOL
      TOLR=RTOL
      FNOISE=TOL
      FNOT=ZERO
      NDIRMX=20
      DO 1 J=1,NDIRMX
      RNOT(J)=ZERO
      CURV(J)=ZERO
      ALPH(J)=ZERO
      DO 1 I=1,NDIRMX
    1 V(I,J)=ZERO
C
      IF (MASWRK) WRITE(IW,9998)
C
C     ----- READ NAMELIST -$TRURST- -----
C
      CALL NAMEIO(IR,KODERR,RNAM,NNAM,QNAM,KQNAM,
     *            KSTART,JSTART,TOLF,TOLR,FNOISE,FNOT,RNOT,CURV,
     *            ALPH,V,
     *            0,0,0,0,
     *   0,0,0,0,0,  0,0,0,0,0,
     *   0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     *   0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0)
      IF(KODERR.EQ.2) THEN
        IF (MASWRK) WRITE(IW,9999)
        CALL ABRT
        STOP
      END IF
        IF (MASWRK) WRITE(IW,9997) KSTART,JSTART,TOLF,TOLR,FNOISE,FNOT
        IX=0
        DO 10 I=1,NDIRMX
        IF(ABS(RNOT(I)).GT.SMALL) IX=IX+1
   10   CONTINUE
      IF((IX.GT.0) .AND. MASWRK) THEN
        WRITE(IW,9996) (RNOT(I),I=1,IX)
        WRITE(IW,9995) (CURV(I),I=1,IX)
        WRITE(IW,9994) (ALPH(I),I=1,IX)
        DO 165 I=1,IX
        IF(I.EQ.1) THEN
          WRITE(IW,9983) (V(J,I),J=1,IX)
        ELSE
          WRITE(IW,9981) (V(J,I),J=1,IX)
        END IF
  165   CONTINUE
      END IF
        IF (MASWRK) WRITE(IW,9992)
      CALL TSECND(TIME)
      IF(TIME.GT.TIMLIM) THEN
         CALL ABRT
         STOP
      END IF
C
C     ----- GET THE PROGRAM STARTED PROPERLY -----
C     ----- THIS SECTION EXECUTES WITH THE INPUT DATA. TO AVOID
C     ----- REPEATING THE SAME CALCULATION, SPECIFY IN THE $TRUDGE
C     ----- GROUP THE P(I) EQUAL TO THE VALUES OF PARAMETERS
C     ----- YOU WANT TRUDGE TO START WITH.
C
C     ----- READ THE $TRUDGE GROUP -----
C
      CALL TRUINP
C
      RSTART=KSTART.GE.0
C
      IF(NDIR.EQ.0) THEN
         CALL ABRT
         STOP
      END IF
C
C           KEEP ON TRUDGING .....
C
      CALL TRUDGE(KSTART)
C
C     ----- PRINT FINAL RESULTS -----
C
      IF(OPTMIZ.EQ.GEOMET) THEN
         CALL TRUHAT
         CALL PRATM(IW,1)
         CALL INTR
      END IF
      IF (MASWRK) WRITE(IW,9990)
      CALL TIMIT(1)
      RETURN
C
 9999 FORMAT(' ERROR IN NAMELIST -TRURST- FOUND. STOP ' )
 9998 FORMAT(/5X,32(1H-)/5X,'NON-GRADIENT ENERGY MINIMIZATION'/
     1           5X,32(1H-))
 9997 FORMAT(' $TRURST KSTART=',I2,'  JSTART=',I2,'  TOLF=',F10.6,
     1       ' TOLR=',F10.6/9X,'FNOISE=',F10.6,'  FNOT=',F14.9)
 9996 FORMAT(9X,'RNOT(1)=',5F10.6/(17X,5F10.6)
     &               /(17X,5F10.6)/(17X,5F10.6))
 9995 FORMAT(9X,'CURV(1)=',5F10.6/(17X,5F10.6)
     &               /(17X,5F10.6)/(17X,5F10.6))
 9994 FORMAT(9X,'ALPH(1)=',5F10.6/(17X,5F10.6)
     &               /(17X,5F10.6)/(17X,5F10.6))
 9983 FORMAT(9X,'   V(1)=',5F10.6/(17X,5F10.6)
     &               /(17X,5F10.6)/(17X,5F10.6))
 9981 FORMAT(17X,5F10.6)
 9992 FORMAT(' $END')
 9990 FORMAT(' ...... END OF NON-GRADIENT ENERGY MINIMIZATION ......')
      END
C*MODULE TRUDGE  *DECK TRUNRG
      SUBROUTINE TRUNRG
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
C
      DIMENSION PPP(10)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      COMMON /FMCOM / XX(1)
      COMMON /FUNCT / E,EG(3,MXATM)
      COMMON /INTNAL/ BOALBE(MXATM,3),SIGN(MXATM),ICONX(MXATM),
     *                IATCON(3,MXATM),INATOM(MXATM),NATIN,IUNT
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /TRUOPT/ OPTMIZ,FUNC,FUNC0,TOLF,R(20),R0(20),TOLR,
     *                P(20),P0(20),IEX(20),NPAR
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
C
      SAVE NSTEP
C
      PARAMETER (ZERO=0.0D+00, PI2=6.28318530717958D+00,
     *           DEGREE=360.0D+00)
C
      DATA NSTEP/0/
      DATA BASIS,GEOMET/8HBASIS   ,8HGEOMETRY/
      DATA CHECK/8HCHECK   /,RNONE/8HNONE    /
C
C     ----- EVALUATE THE ENERGY (NRG) FOR TRUDGE OPTIMIZATIONS -----
C
      IDIR=0
      DO 110 IPAR=1,NPAR
         IF(IEX(IPAR).EQ.0) GO TO 110
         IDIR=IDIR+1
C
         IF(OPTMIZ.EQ.BASIS) THEN
           P(IPAR)=P0(IPAR)*EXP(R(IDIR))
           EX(IEX(IPAR))=P(IPAR)
         ELSE
           IATOM=IEX(IPAR)/10
           IHINT=IEX(IPAR) - 10*IATOM
           P(IPAR)=P0(IPAR)*EXP(R(IDIR))
           BOALBE(IATOM,IHINT)=P(IPAR)
         END IF
  110 CONTINUE
C
      IF(RMETHOD.NE.RNONE  .AND.  MOD(MODQR,2).EQ.1) THEN
         MAXGTOT = MXGTOT
         MAXSH   = MXSH
         CALL VALFM(LOADFM)
         LEXU     = LOADFM   + 1
         LCSU     = LEXU     + MAXGTOT
         LCPU     = LCSU     + MAXGTOT
         LCDU     = LCPU     + MAXGTOT
         LCFU     = LCDU     + MAXGTOT
         LCGU     = LCFU     + MAXGTOT
         LKSTARTU = LCGU     + MAXGTOT
         LKATOMU  = LKSTARTU + MAXSH/NWDVAR + 1
         LKTYPEU  = LKATOMU  + MAXSH/NWDVAR + 1
         LKNGU    = LKTYPEU  + MAXSH/NWDVAR + 1
         LKLOCU   = LKNGU    + MAXSH/NWDVAR + 1
         LKMINU   = LKLOCU   + MAXSH/NWDVAR + 1
         LKMAXU   = LKMINU   + MAXSH/NWDVAR + 1
         LMAPLE   = LKMAXU   + MAXSH/NWDVAR + 1
         LLARCH   = LMAPLE   + MAXGTOT/NWDVAR + 1
         LELM     = LLARCH   + MAXGTOT/NWDVAR + 1
         LAST     = LELM     + MAXGTOT/NWDVAR + 1
         NEED = LAST - LOADFM - 1
         CALL GETFM(NEED)
         CALL BUNCON(MAXGTOT,MAXSH,XX(LEXU),XX(LCSU),XX(LCPU),
     *               XX(LCDU),XX(LCFU),XX(LCGU),XX(LKSTARTU),
     *               XX(LKATOMU),XX(LKTYPEU),XX(LKNGU),XX(LKLOCU),
     *               XX(LKMINU),XX(LKMAXU),XX(LMAPLE),XX(LLARCH),
     *               XX(LELM))
         CALL RETFM(NEED)
      END IF
C
      IF(OPTMIZ.EQ.GEOMET) THEN
        CALL TRUATO
        CALL TRUHAT
      END IF
C
C        PRINT A RESTART VERSION OF $TRUDGE
C
      DO 200 IPAR=1,NPAR
         IF(OPTMIZ.EQ.GEOMET) THEN
            IF(IEX(IPAR).GT.0) THEN
               IATOM=IEX(IPAR)/10
               IHINT=IEX(IPAR) - 10*IATOM
               IF(IHINT.EQ.1) THEN
                  PPP(IPAR)=P(IPAR)
               ELSE
                  PPP(IPAR)=P(IPAR)*DEGREE/PI2
               END IF
            END IF
         ELSE
            PPP(IPAR) = P(IPAR)
         END IF
  200 CONTINUE
      IF (MASWRK) THEN
         WRITE(IW,9996) NSTEP
         WRITE(IW,9995) OPTMIZ,NPAR
         WRITE(IW,9994) (IEX(I),I=1,NPAR)
         WRITE(IW,9993) (PPP(I),I=1,NPAR)
         WRITE(IW,9992)
      END IF
C
C     ----- GET THE ENERGY OF THE WAVEFUNCTION -----
C
      CALL ENERGX
C
      FUNC = E
      IF(FUNC.GE.ZERO .AND. EXETYP.NE.CHECK) THEN
         IF(MASWRK) WRITE(IW,*) 'SCF DID NOT CONVERGE, TRUDGE ABORTED'
         CALL ABRT
      END IF
      IF(MASWRK) WRITE(IW,9080) NSTEP,FUNC
      NSTEP = NSTEP+1
C
C     ----- TURN OFF MOST PRINTOUT -----
C
      NPRINT=-5
      RETURN
C
 9996 FORMAT(1X,'----- TRUDGE RESTART DATA AT NSTEP',I4,' --------')
 9995 FORMAT(1X,'$TRUDGE OPTMIZ=',A8,'  NPAR=',I2)
 9994 FORMAT(1X,'        IEX(1)=',10(I5,1H,)/16X,10(I5,1H,))
 9993 FORMAT(1X,'        P(1)=',4(F14.6,1H,)/1X,13X,4(F14.6,1H,)
     &                  /1X,13X,4(F14.6,1H,)/1X,13X,4(F14.6,1H,)/
     &                  /1X,13X,4(F14.6,1H,))
 9992 FORMAT(1X,'$END')
 9080 FORMAT(/1X,'TRUDGE ENERGY VALUE AT NSTEP=',I4,' IS',F20.10)
      END
C*MODULE TRUDGE  *DECK TRUINP
      SUBROUTINE TRUINP
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
      PARAMETER (NNAM=4)
C
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
C
      COMMON /INTNAL/ BOALBE(MXATM,3),SIGN(MXATM),ICONX(MXATM),
     *                IATCON(3,MXATM),INATOM(MXATM),NATIN,IUNT
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /TRUOPT/ OPTMIZ,FUNC,FUNC0,TOLF,R(20),R0(20),TOLR,
     *                P(20),P0(20),IEX(20),NPAR
      COMMON /TRUPAR/ DUMMY(463),NDIR,NDIRMX
C
      DATA RNAM/8HTRUDGE  /
      DATA QNAM/8HOPTMIZ  ,8HNPAR    ,8HIEX     ,8HP       /
      DATA KQNAM/5,1,201,203/
C
      DATA ZERO/0.0D+00/
      DATA PI2,DEGREE /6.28318530717958D+00,360.0D+00/
      DATA BASIS,GEOMET/8HBASIS   ,8HGEOMETRY/
C
      OPTMIZ=GEOMET
      NPAR=20
      DO 100 I=1,NPAR
         IEX(I) = 0
         P(I) = ZERO
  100 CONTINUE
C
      CALL NAMEIO(IR,KODERR,RNAM,NNAM,QNAM,KQNAM,
     *            OPTMIZ,NPAR,IEX,P,
     *   0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     *   0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     *   0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0)
      IF(KODERR.EQ.1) THEN
         IF (MASWRK) WRITE(IW,9996)
         CALL ABRT
         STOP
      END IF
      IF(KODERR.EQ.2) THEN
         IF (MASWRK) WRITE(IW,9998)
         CALL ABRT
         STOP
      END IF
C
C         COMMON /INTNAL/ EXPECTS TO HAVE DISTANCES IN ANGSTROMS,
C         AND ANGLES IN RADIANS.  CONVERT THE INPUT ANGLES.
C
      IF(OPTMIZ.EQ.GEOMET) THEN
         DO 210 IPAR=1,NPAR
            IF(IEX(IPAR).EQ.0) GO TO 210
            IATOM=IEX(IPAR)/10
            IHINT=IEX(IPAR) - 10*IATOM
            IF(IHINT.NE.1) THEN
               P(IPAR)=P(IPAR)*PI2/DEGREE
            END IF
  210    CONTINUE
      END IF
C
      IDIR=0
      DO 310 IPAR=1,NPAR
         IF(IEX(IPAR).EQ.0) GO TO 310
         IF(OPTMIZ.EQ.BASIS) THEN
            IF(P(IPAR).EQ.ZERO) P (IPAR)=EX(IEX(IPAR))
                                P0(IPAR)=EX(IEX(IPAR))
         ELSE
            IATOM=IEX(IPAR)/10
            IHINT=IEX(IPAR) - 10*IATOM
            IF(P(IPAR).EQ.ZERO)  P(IPAR)=BOALBE(IATOM,IHINT)
                                P0(IPAR)=BOALBE(IATOM,IHINT)
         END IF
         IDIR=IDIR+1
         R(IDIR)= LOG(P(IPAR)/P0(IPAR))
  310 CONTINUE
      NDIR=IDIR
      RETURN
C
 9998 FORMAT(1X,'ERRORS IN NAMELIST $TRUDGE FOUND.')
 9996 FORMAT(1X,'NO NAMELIST $TRUDGE FOUND.')
      END
C*MODULE TRUDGE  *DECK TRUDGE
      SUBROUTINE TRUDGE(KDIR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL RSTART,EXACT,GOPARR,DSKWRK,MASWRK
      DIMENSION U(10),P(10)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /TRUDIR/ IROOT,JDIR,EXACT
      COMMON /TRUEIG/ G(20,20),Q(20,20),SCRAP1(20),SCRAP2(20),SCRAP3(20)
      COMMON /TRUOPT/ OPTMIZ,F,FNOT,TOLF,R(20),RNOT(20),TOLR,
     *                PP(20),P0(20),IEX(20),NPAR
      COMMON /TRUPAR/ DR(20),ALPH(20),CURV(20),V(20,20),ALPHA,CURVE,
     *                FNOISE,NDIR,NDIRMX
      COMMON /TRURUN/ TIME,RSTART
C
      PARAMETER (SMALL=1.0D-12)
C
      DATA CHECK/8HCHECK   /
      DATA ZERO,ONE,TWO,FOUR,EIGHT /0.0D+00,1.0D+00,2.0D+00,4.0D+00,
     1 8.0D+00/
      DATA PT1,PT25,PT5,PT7 /0.1D+00,0.25D+00,0.5D+00,0.7D+00/
C
C                  ****  VERSION  FEB. 1980  ****
C
C   PROGRAM MINIMIZES A FUNCTION OF NDIR VARIABLES BY A MODIFIED POWELL
C   METHOD OF SEARCHES ALONG CONJUGATE DIRECTIONS. COLUMNS OF V MATRIX
C   ARE CURRENT SEARCH DIRECTIONS.  COLUMNS 1 TO KDIR ARE (APPROXIMATELY
C   CONJUGATE.
C
C   WHEN SUBROUTINE IS ENTERED WITH RSTART=.TRUE. THEN VALUES OF KDIR
C   AND OF ALL QUANTITIES IN COMMON  SHOULD HAVE BEEN SET PREVIOUSLY
C   BY THE CALLING PROGRAM.
C   WHEN TRUDGE IS CALLED WITH KDIR NEGATIVE THEN ONLY REQUIRED VALUES
C   ARE TIMLIM, TOLF, TOLR, NOISE, NDIR, AND R. (NOISE NEEDED FOR SEARCH
C   IN THIS CASE INITIAL VALUES OF KDIR,JDIR,RSTART,ALPH,V ARE IGNORED
C   IF TRUDGE IS CALLED WITH NON-NEGATIVE KDIR AND RSTART=.FALSE. THEN
C   VALUES OF TIMLIM, TOLF, TOLR, NOISE, NDIR, R, AND V ARE REQUIRED.
C   IN THIS CASE KDIR IS RESET TO ZERO AFTER FIRST PASS.
C
      IF(NDIR.GT.NDIRMX) THEN
         IF (MASWRK) WRITE(IW,*) 'TRUDGE: NDIR,NDIRMX=',NDIR,NDIRMX
         CALL ABRT
      END IF
      FN=NDIR
      IF(KDIR.LT.0) GO TO 15
      IF(.NOT.RSTART) CALL TRUNRG
      DO 10 I=1,NDIR
   10 IF(ALPH(I).LE.ZERO) ALPH(I)=PT1
      IF(JDIR.LE.0.OR.JDIR.GT.NDIR) JDIR=NDIR
      RSTART=RSTART.AND.KDIR.NE.NDIR
      GO TO 100
C
C     IF K NEGATIVE SET INITIAL V MATRIX TO IDENTITY, STANDARD VALUES.
C
   15 CONTINUE
      KDIR=0
      RSTART=.FALSE.
      CALL TRUNRG
      IF(EXETYP.EQ.CHECK) RETURN
      DO 30 I=1,NDIR
      ALPH(I)=PT1
      DO 20 J=1,NDIR
   20 V(I,J)=ZERO
      CURV(I)=ZERO
   30 V(I,I)=ONE
      JDIR=NDIR
C
C     SEARCH ALONG JTH DIRECTION FOR LOCAL MINIMUM.
C
  100 CONTINUE
      DO 110 I=1,NDIR
  110 DR(I)=V(I,JDIR)
      EXACT=JDIR.LE.KDIR
      ALPHA=ALPH(JDIR)
      ALFSET=ALPHA
      CALL TRUSRC(LEAVE)
      IF (MASWRK) THEN
         WRITE(IW,9999) CURVE,ALPHA
         IF(LEAVE.NE.0) WRITE(IW,9997)
      END IF
 9999 FORMAT(17H  NEW CURVATURE =,E13.5,10H   ALPHA =,F12.7)
 9998 FORMAT(17H  OLD CURVATURE =,E13.5,10H   ALFSET=,F12.7)
 9997 FORMAT(20X,30H ABNORMAL RETURN FROM SEARCH  )
      ALPHA=ABS(ALPHA)
      IF(ALPHA.LT. PT5*ALFSET) ALFSET=PT5*ALFSET
      IF(ALPHA.LT.PT25*ALFSET) ALFSET=PT5*ALFSET
      IF(ALPHA.GT. TWO*ALFSET) ALFSET=TWO*ALFSET
      IF(ALPHA.GT.FOUR*ALFSET) ALFSET=TWO*ALFSET
      IF (MASWRK) WRITE(IW,9998) CURV(JDIR),ALFSET
      CURV(JDIR)=CURVE
      ALPH(JDIR)=ALFSET
      IF(JDIR.EQ.1) GO TO 175
      JDIR=JDIR-1
      CALL TSECND(TIME)
      IF(TIME.LE.TIMLIM) GO TO 100
C
C     CALL FUNC(-1) CAUSES ALL RSTART DATA TO BE OUTPUT,
C     THIS INCLUDES ALL DATA REQUIRED FOR A FUTURE CALL FUNC(1).
C
C     END SECTION: ENTER HERE IF
C        (1) ALL DONE
C        (2) TIME EXCEEDED LIMIT
C
  150 CONTINUE
      IF (MASWRK) WRITE(IW,9995) KDIR,JDIR,F,FNOT
 9995 FORMAT(/,' CONJUGATE DIRECTION = ',I5,/,' PARAMETER = ',I5,
     1 /,' VALUE OF THE FUNCTION = ',F20.10,/,17X,' FNOT = ',F20.10)
      KSTART=KDIR
      JSTART=JDIR
      IF (MASWRK) WRITE(IW,9987) KSTART,JSTART,TOLF,TOLR,FNOISE,FNOT
      IX=0
      DO 160 I=1,NDIRMX
        IF(ABS(RNOT(I)).GT.SMALL) IX=IX+1
  160 CONTINUE
      IF((IX.GT.0) .AND. MASWRK) THEN
        WRITE(IW,9986) (RNOT(I),I=1,IX)
        WRITE(IW,9985) (CURV(I),I=1,IX)
        WRITE(IW,9984) (ALPH(I),I=1,IX)
        DO 165 I=1,IX
        IF(I.EQ.1) THEN
          WRITE(IW,9983) (V(J,I),J=1,IX)
        ELSE
          WRITE(IW,9981) (V(J,I),J=1,IX)
        END IF
  165   CONTINUE
      END IF
      IF (MASWRK) WRITE(IW,9982)
C
 9987 FORMAT(' $TRURST KSTART=',I2,'  JSTART=',I2,'  TOLF=',F10.6,
     1       ' TOLR=',F10.6/9X,'FNOISE=',F10.6,'  FNOT=',F14.9)
 9986 FORMAT(9X,'RNOT(1)=',5F10.6/(17X,5F10.6)
     &               /(17X,5F10.6)/(17X,5F10.6))
 9985 FORMAT(9X,'CURV(1)=',5F10.6/(17X,5F10.6)
     &               /(17X,5F10.6)/(17X,5F10.6))
 9984 FORMAT(9X,'ALPH(1)=',5F10.6/(17X,5F10.6)
     &               /(17X,5F10.6)/(17X,5F10.6))
 9983 FORMAT(9X,'   V(1)=',5F10.6/(17X,5F10.6)
     &               /(17X,5F10.6)/(17X,5F10.6))
 9981 FORMAT(17X,5F10.6)
 9982 FORMAT(' $END')
C
      IF(TIME.GE.TIMLIM) THEN
         IF(MASWRK) WRITE(IW,*) '*** THIS TRUDGE DID NOT CONVERGE ***'
         IF(MASWRK) WRITE(IW,*) 'RESTART WITH $TRUDGE + $TRURST ABOVE.'
      END IF
      RETURN
C
C     END OF SEQUENCE OF SEARCHES ALONG ALL DIRECTIONS IN V MATRIX.
C
  175 CONTINUE
      IF(RSTART) GO TO 185
      DO 180 I=1,NDIR
  180 RNOT(I)=R(I)
      RSTART=.TRUE.
      FNOT=F
      KDIR=0
      JDIR=NDIR
      GO TO 100
C
C     GENERATE NEW CONJUGATE DIRECTION IF NOT CONVERGED.
C
  185 CONTINUE
      DRDR=ZERO
      DO 190 I=1,NDIR
      DUM=R(I)-RNOT(I)
      DR(I)=DUM
  190 DRDR=DUM*DUM+DRDR
      ALFSET=SQRT(DRDR)
      ALPHA=ALFSET
      IF(ALPHA.EQ.ZERO) GO TO 200
      DO 195 I=1,NDIR
      U(I)=DR(I)/ALFSET
  195 DR(I)=U(I)
      EXACT=.TRUE.
      CALL TRUSRC(LEAVE)
      IF (MASWRK) THEN
         IF(LEAVE.NE.0) WRITE(IW,9997)
         WRITE(IW,9999) CURVE,ALPHA
      END IF
      ALPHA=ABS(ALPHA)
  200 DELF=FNOT-F
      FNOT=F
      DRDR=ZERO
      DO 210 I=1,NDIR
      DRDR=(R(I)-RNOT(I))**2+DRDR
  210 RNOT(I)=R(I)
      DELR=SQRT(DRDR/FN)
      IF (MASWRK) WRITE(IW,9991) F,DELF,DELR
 9991 FORMAT(' F=',F20.10,' DELTA F=',E11.5,' DELTA R=',E11.5)
      CALL TRUNRG
C
      IF(DELF.GT.TOLF.OR.DELR.GT.TOLR) GO TO 220
      IF (MASWRK) WRITE(IW,9990)
 9990 FORMAT(/,' NON-GRADIENT ENERGY MINIMIZATION  ...  CONVERGED')
      GO TO 150
  220 CONTINUE
      IF(KDIR.GT.0) GO TO 230
      KDIR=1
      K1=0
      DO 225 I=1,NDIR
  225 R(I)=DR(I)
      DR(NDIR)=CURVE
      P(1)=ALPHA
      GO TO 375
C
C     CALCULATE K+1 Q VECTORS. Q(I) DOT V(J)=KRONECKER DELTA(I,J), J=1,K
C
  230 CONTINUE
      DO 250 J=1,KDIR
      DOT=ZERO
      DO 240 I=1,NDIR
  240 DOT=DR(I)*V(I,J)+DOT
      P(J)=DOT
      DO 250 I=1,NDIR
  250 DR(I)=DR(I)-DOT*V(I,J)
      DOT=ZERO
      DO 260 I=1,NDIR
  260 DOT=DR(I)*U(I)+DOT
      K1=KDIR
      KDIR=K1+1
      RSTART=KDIR.NE.NDIR
      DO 270 I=1,NDIR
      DUM=DR(I)/DOT
      Q(I,KDIR)=DUM
      DO 270 J=1,K1
  270 Q(I,J)=V(I,J)-P(J)*DUM
C
C     G MATRIX=SUM L=1,K+1 OF Q(L)*CURV(L)*TRANSPOSE Q(L).
C
      DO 280 I=1,NDIR
      DUM=CURVE*Q(I,KDIR)
      DO 280 J=1,NDIR
  280 G(I,J)=DUM*Q(J,KDIR)
      DO 300 L=1,K1
      DO 300 I=1,NDIR
      DUM=CURV(L)*Q(I,L)
      DO 300 J=1,NDIR
  300 G(I,J)=DUM*Q(J,L)+G(I,J)
      CALL TRUDIA(NDIR,DR)
C
C     SET P(J) EQUAL TO CORRESPONDING ALPH VALUE.
C
      DO 350 J=1,KDIR
      JJ=NDIR+1-J
      P(J)=ALPH(J)
      DOT=ZERO
      DO 310 I=1,NDIR
      R(I)=Q(I,JJ)
  310 DOT=R(I)*V(I,J)+DOT
      DOT=ABS(DOT)
      IF(DOT.GT.PT7) GO TO 350
      LMAX=J
      DOTMAX=DOT
      DO 340 L=1,K1
      IF(L.EQ.J) GO TO 340
      DOT=ZERO
      DO 320 I=1,NDIR
  320 DOT=R(I)*V(I,L)+DOT
      DOT=ABS(DOT)
      IF(DOT.LT.PT7) GO TO 330
      P(J)=ALPH(L)
      GO TO 350
  330 IF(DOT.LE.DOTMAX) GO TO 340
      DOTMAX=DOT
      LMAX=L
  340 CONTINUE
      P(J)=ALPH(LMAX)
      DOT=ZERO
      DO 345 I=1,NDIR
  345 DOT=R(I)*U(I)+DOT
      IF(ABS(DOT).GT.DOTMAX) P(J)=ALPHA
  350 CONTINUE
C
C     SET NEW CONJUGATE VECTORS IN V EXCEPT KTH VECTOR IS SAVED IN R.
C
      DO 370 J=1,K1
      JJ=NDIR+1-J
      DO 360 I=1,NDIR
  360 V(I,J)=Q(I,JJ)
      ALPH(J)=P(J)
  370 CURV(J)=DR(JJ)
C
C     ORTHOGONALIZE REMAINING SEARCH DIRECTIONS TO CONJUGATE DIRECTIONS.
C
  375 CONTINUE
      IF(KDIR.EQ.NDIR) GO TO 450
      DO 420 J=KDIR,NDIR
      DOT=ZERO
      DO 380 I=1,NDIR
  380 DOT=V(I,J)*R(I)+DOT
      DO 390 I=1,NDIR
  390 V(I,J)=V(I,J)-DOT*R(I)
      IF(KDIR.EQ.1) GO TO 420
      DO 410 L=1,K1
      DOT=ZERO
      DO 400 I=1,NDIR
  400 DOT=V(I,J)*V(I,L)+DOT
      DO 410 I=1,NDIR
  410 V(I,J)=V(I,J)-DOT*V(I,L)
  420 CONTINUE
C
C     CONSTRUCT G MATRIX FOR REMAINING SEARCH DIRECTIONS.
C
      CMAX=ZERO
      DO 425 L=KDIR,NDIR
  425 IF(CURV(L).GT.CMAX) CMAX=CURV(L)
      CMAX=EIGHT*CMAX
      DO 440 I=1,NDIR
      DO 440 J=1,I
      DUM=-CMAX*R(I)*R(J)
      DO 430 L=KDIR,NDIR
  430 DUM=CURV(L)*V(I,L)*V(J,L)+DUM
      IF(K1.LE.0) GO TO 435
      DO 433 L=1,K1
  433 DUM=-CMAX*V(I,L)*V(J,L)+DUM
  435 G(I,J)=DUM
  440 G(J,I)=DUM
C
C     SET KTH CONJUGATE DIRECTION BACK INTO V MATRIX.
C
  450 CONTINUE
      DO 460 I=1,NDIR
      V(I,KDIR)=R(I)
  460 R(I)=RNOT(I)
      ALPH(KDIR)=P(KDIR)
      CURV(KDIR)=DR(NDIR+1-KDIR)
      IF(KDIR.EQ.NDIR) GO TO 480
      CALL TRUDIA(NDIR,DR)
C
C     SET REMAINING SEARCH DIRECTIONS INTO V MATRIX.
C
      K1=KDIR+1
      DO 470 J=K1,NDIR
      CURV(J)=DR(J)
      DO 470 I=1,NDIR
  470 V(I,J)=Q(I,J)
  480 CONTINUE
      JDIR=NDIR
      CALL TSECND(TIME)
      IF(TIME.LE.TIMLIM) GO TO 100
      IF(MASWRK) WRITE(IW,*) 'THIS RUN HAS EXPIRED ITS TIME LIMIT'
      GO TO 150
      END
C*MODULE TRUDGE  *DECK TRUSRC
      SUBROUTINE TRUSRC(LEAVE)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL SET(9),BACK,EXACT,HOLD,GOPARR,DSKWRK,MASWRK
C
      DIMENSION XNOT(20),F(9)
C
      COMMON /IOFILE/ IN,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /TRUOPT/ OPTMIZ,FUNC,FUNC0,TOLF,X(20),X0(20),TOLX,
     *                P(20),P0(20),IEX(20),NPAR
      COMMON /TRUPAR/ DX(20),ALPH(20),CURV(20),V(20,20),AMIN,CURVE,
     *                FNOISE,NDIR,NDIRMX
      COMMON /TRUDIR/ IROOT,JDIR,EXACT
C
      DATA ZERO,ONE,TWO,THREE,FOUR,FIVE,SIX
     1 /0.0D+00,1.0D+00,2.0D+00,3.0D+00,4.0D+00,5.0D+00,6.0D+00/
      DATA TEN3,TEN35 /1.0D+03,1.0D+35/
      DATA PT1875,PT375,PT125,PT025,PT25,PT5,PT75
     1 /0.1875D+00,0.375D+00,0.125D+00,0.025D+00,0.25D+00,0.5D+00,
     2  0.75D+00/
C
C                  ****  VERSION  FEB. 1980  ****
C
C     SUBROUTINE LOCATES A LOCAL MINIMUM ALONG A GIVEN ONE-DIMENSIONAL
C     SEARCH DIRECTION.
C     F(ALF) = FUNC(X+ALF*DX)  FIND OPTIMUM ALF VALUE.
C     LOOKMX = MAXIMUM NUMBER OF FUNCTION EVALUATIONS BEFORE SEARCH IS
C               TERMINATED.
C     ENTER SUBROUTINE WITH FOLLOWING
C     N     = DIMENSIONALITY OF X SPACE
C     X     = INITIAL POSITION (IN N-DIMENSIONAL PARAMETER SPACE)
C     FUNC  = F(0), I.E. FUNCTION VALUE AT INITIAL X.
C     DX    = SEARCH DIRECTION (IN N-DIMENSIONAL SPACE)
C     AMIN  = ESTIMATED INITIAL STEP SIZE
C     NOISE = ACCURACY OF FUNC VALUES.   VARIATIONS SMALLER THAN NOISE
C             ARE NOT CONSIDERED TO BE SIGNIFICANT.
C     EXACT = .FALSE. SUPPRESSES FINAL CHECK OF QUADRATIC FIT.
C
C     EXIT FROM SUBROUTINE WITH FOLLOWING
C     X     = CALCULATED MINIMUM
C     FUNC  = F(AMIN), I.E. FUNCTION VALUE AT FINAL X.
C     AMIN  = CALCULATED ALF VALUE AT MINIMUM
C     CURVE = CALCULATED CURVATURE ALONG SEARCH DIRECTION
C     LEAVE = 0 IF NORMAL EXIT
C           = 2 IF EXCESSIVE NUMBER OF FUNCTION EVALUATIONS
C
C     NDIR,NDIRMX,FUNC,AND X ARE SHARED WITH SUBROUTINE -TRUNRG-.
C     ALL VARIABLES IN COMMON ARE SHARED WITH PROGRAM CALLING -TRUSRC -
C
      DO 2 I=1,9
         F(I) = ZERO
    2 CONTINUE
      IF(EXACT) THEN
         XFACT=ONE
      ELSE
         XFACT=TWO
      END IF
      DO 5 I=1,NDIR
    5 XNOT(I)=X(I)
      LOOKS=0
      LEAVE=0
      LOOKMX=30
      ALF=AMIN
      AMIN=ZERO
      FNOT=FUNC
      FMIN=FUNC
      FNEXT=TEN35
      IR=4
C
   10 LOOKS=LOOKS+1
      IF(LOOKS.LE.LOOKMX) GO TO 25
         LEAVE=2
         IF (MASWRK) WRITE(IW,9999) LOOKS
 9999    FORMAT(/'  EXCESSIVE NUMBER OF FUNCTION EVALUATIONS =',I3/)
C
C     EXIT FROM SUBROUTINE
C
   15    CONTINUE
         DO 20 I=1,NDIR
   20    X(I)=XNOT(I)+AMIN*DX(I)
         FUNC=FMIN
         RETURN
C
C     EVALUATE F(ALF)=FUNC(XNOT+ALF*DX)
C
   25 CONTINUE
      DO 30 I=1,NDIR
   30 X(I)=XNOT(I)+ALF*DX(I)
      CALL TRUNRG
C
      HOLD=.FALSE.
      IF(FUNC.LT.FMIN) THEN
         FNEXT=FMIN
         FMIN=FUNC
         ANEXT=AMIN
         AMIN=ALF
      ELSE IF(FUNC.LT.FNEXT) THEN
         FNEXT=FUNC
         ANEXT=ALF
      END IF
      IF (MASWRK) WRITE(IW,9998) LOOKS,ALF,FUNC,AMIN,FMIN,ANEXT,FNEXT
 9998 FORMAT(' LOOKS =',I3,'   ALF =',F8.4,'  FUNC =',F20.10/
     *                 11X,'  AMIN =',F8.4,'  FMIN =',F20.10/
     *                 11X,' ANEXT =',F8.4,' FNEXT =',F20.10)
      CALL FLSHBF(IW)
      GO TO (150,180,260,320) IR-3
C
C     F1,F3,F5,F7,F9 ARE FUNC VALUES AT EVENLY SPACED POINTS WITH
C     SEPARATION DELTA.  SET(I)=TRUE INDICATES THAT  F(I) HAS BEEN
C     EVALUATED.  DELTA MAY INCREASE OR DECREASE BY FACTORS OF 2.0.
C     F VALUES SHIFT LEFT OR RIGHT SO THAT F5 IS ALWAYS LOWEST
C
  150 CONTINUE
      DELTA=ABS(AMIN-ANEXT)
      BACK=.FALSE.
      DO 155 I=1,9
  155 SET(I)=.FALSE.
      F(5)=FMIN
      SET(5)=.TRUE.
      K=3
      IF(AMIN.LT.ANEXT) K=7
      F(K)=FNEXT
      SET(K)=.TRUE.
C
  160 IF(SET(7)) GO TO 230
      IF(SET(3)) THEN
         IF(SET(1)) THEN
            ALF=AMIN+DELTA
            IR=5
            GO TO 10
C
         END IF
C
C        EXTRAPOLATE USING F1, F3, AND F5.
C        INCREASE STEP SIZE IF F9 IS PREDICTED TO BE LOWER THAN F7.
C        DECREASE STEP SIZE IF F4 IS PREDICTED TO BE LOWER THAN F5.
C
         BACK=FOUR*F(3).LT.(F(1)+THREE*F(5))
         IF((.NOT.(SET(9).OR.HOLD)) .AND.
     *      (FIVE*F(3).GT.(TWO*F(1)+THREE*F(5)))) GO TO 220
         IF(.NOT.(SET(4) .OR. (SIX*F(3).GT.(F(1)+FIVE*F(5))))) GO TO 305
         ALF=AMIN+DELTA
         IR=5
         GO TO 10
C
      END IF
      IF(BACK) THEN
         ALF=AMIN-DELTA
         IR=6
      ELSE
         ALF=AMIN+DELTA
         IR=5
      END IF
      GO TO 10
C
  180 F(7)=FUNC
      SET(7)=.TRUE.
      IF(ALF.NE.AMIN) GO TO 160
C     SHIFT LEFT
      DO 190 I=1,7
         SET(I)=SET(I+2)
         F(I)=F(I+2)
  190 CONTINUE
      SET(8)=.FALSE.
      SET(9)=.FALSE.
      GO TO 160
C
  220 DELTA=DELTA+DELTA
      F(4)=F(3)
      F(3)=F(1)
      F(6)=F(7)
      F(7)=F(9)
      SET(4)=SET(3)
      SET(3)=SET(1)
      SET(6)=SET(7)
      SET(7)=SET(9)
      SET(1)=.FALSE.
      SET(2)=.FALSE.
      SET(8)=.FALSE.
      SET(9)=.FALSE.
      GO TO 160
C
  230 IF(SET(3)) THEN
C
C     INTERPOLATE USING F3, F5, AND F7.
C     EITHER IMMEDIATELY ATTEMPT QUADRATIC FIT OR FIRST REDUCE STEP SIZE
C     IF AMIN=ZERO STEP SIZE MAY BE  REDUCED BY A FACTOR OF 2,4,8, OR 16
C     IN OTHER CASES STEP SIZE IS REDUCED BY A FACTOR OF 2.
C
      SPREAD=MAX(F(3),F(7))-F(5)
      TOL=MAX((FNOT-FMIN)*PT5,TEN3*FNOISE)*XFACT
      DUM=F(3)-F(5)-F(5)+F(7)
      CURVE=DUM/(TWO*DELTA*DELTA)
      DUM=(F(3)-F(7))/DUM
      BACK=DUM.LT.ZERO
      IF(SPREAD.LT.TOL) THEN
C
C        ATTEMPT QUADRATIC FIT
C
         ALF=AMIN+DUM*DELTA*PT5
         FQUAD=F(5)-(F(3)-F(7))*DUM*PT125
         IR=7
         GO TO 10
C
      END IF
      IF(ABS(DUM).GT.PT75.OR.AMIN.NE.ZERO) GO TO 305
      SET(1)=SET(4)
      F(1)=F(4)
      SET(9)=SET(6)
      F(9)=F(6)
      SET(2)=.FALSE.
      SET(3)=.FALSE.
      SET(4)=.FALSE.
      SET(6)=.FALSE.
      SET(7)=.FALSE.
      SET(8)=.FALSE.
      DELTA=PT25*DELTA
      IF(ABS(DUM).GT.PT375) GO TO 160
      SET(1)=.FALSE.
      SET(9)=.FALSE.
      DELTA=PT5*DELTA
      IF(ABS(DUM).LE.PT1875) DELTA=PT5*DELTA
      GO TO 160
C
      END IF
      IF(SET(9)) THEN
C
C     EXTRAPOLATE USING F5, F7, AND F9.
C     INCREASE STEP SIZE IF F1 IS PREDICTED TO BE LOWER THAN F3.
C     DECREASE STEP SIZE IF F6 IS PREDICTED TO BE LOWER THAN F5.
C
      BACK=FOUR*F(7).GT.(F(9)+THREE*F(5))
      IF((.NOT.(SET(1).OR.HOLD)) .AND.
     *   (FIVE*F(7).GT.(TWO*F(9)+THREE*F(5)))) GO TO 220
      IF((.NOT.SET(6)) .AND.
     *   (SIX*F(7).LE.(F(9)+FIVE*F(5)))) GO TO 305
      END IF
      ALF=AMIN-DELTA
      IR=6
      GO TO 10
C
  260 F(3)=FUNC
      SET(3)=.TRUE.
      IF(ALF.EQ.AMIN) THEN
C     SHIFT RIGHT
      DO 270 I=1,7
         SET(10-I)=SET(8-I)
         F(10-I)=F(8-I)
  270 CONTINUE
      SET(1)=.FALSE.
      SET(2)=.FALSE.
      END IF
      GO TO 160
C
  320 IF(.NOT.EXACT) GO TO 15
      TOL=(FNOT-FMIN)*PT025
      IF(TOL.LT.FNOISE) TOL=FNOISE
      IF(ABS(FUNC-FQUAD).LT.TOL) GO TO 15
      IF(ALF.EQ.AMIN.OR.ALF.EQ.ANEXT) GO TO 150
C
C     CUT STEP SIZE IN HALF
C
  305 DELTA=PT5*DELTA
      HOLD=.TRUE.
      SET(1)=SET(3)
      SET(3)=SET(4)
      SET(9)=SET(7)
      SET(7)=SET(6)
      SET(2)=.FALSE.
      SET(4)=.FALSE.
      SET(6)=.FALSE.
      SET(8)=.FALSE.
      F(1)=F(3)
      F(3)=F(4)
      F(9)=F(7)
      F(7)=F(6)
      GO TO 160
      END
C*MODULE TRUDGE  *DECK TRUDIA
      SUBROUTINE TRUDIA(N,EIG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C                  ****  VERSION  3/20/73  ****
C        -----           DESCRIPTION                       -----
C        -----   THIS IS A GIVENS-HOUSHOLDER ROUTINE       -----
C        -----   RECEIVES A MATRIX-A-OF WHICH ONLY LOWER   -----
C        -----   TRANGLE IS USED AND THIS IS DESTROYED     -----
C        -----   EIGENVALUES ARE RETURNED IN ALGEBRAIC     -----
C        -----   ASCENDING ORDER IN ARRAY-EIG-             -----
C        -----   THE EIGENVECTORS ARE RETURNED STORED BY   -----
C        -----   COLUMNS IN VEC.                           -----
C        -----         PARAMETERS PASSED                   -----
C        -----   RHO=TOLERANCE ON FINAL OFF-DIAG ELEMENTS  -----
C        -----   N  IS THE SIZE OF THE MATRIX              -----
C        -----   THIS ROUTINE IS SLIGHTLY MODIFIED VERSION -----
C        -----   OF THAT IN QCPE PROGRAM CINDO BY DOBOSH   -----
C
      DIMENSION W(20),EIG(1),IPOSV(20),IVPOS(20),P(20),Q(20),IORD(20)
      COMMON /TRUEIG/ A(20,20),VEC(20,20),GAMMA(20),BETA(20),BETASQ(20)
      EQUIVALENCE (IPOSV,GAMMA),(IVPOS,BETA),(IORD,BETASQ),
     1 (IVPOS,P),(IVPOS,Q)
      DATA RHO /1.0D-10/
      DATA ZERO,ONE,TWO,PT5 /0.0D+00,1.0D+00,2.0D+00,0.5D+00/
C
      RHOSQ=RHO*RHO
      IF(N.EQ.0) GO TO 560
      SHIFT=ZERO
      N1=N-1
      N2=N-2
      GAMMA(1)=A(1,1)
      IF(N2.EQ.0) GO TO 270
      IF(N2.LT.0) GO TO 280
      DO 260 NR=1,N2
      B=A(NR+1,NR)
      S=ZERO
      DO 130 I=NR,N2
  130 S=S+A(I+2,NR)**2
C        -----  PREPARE FOR POSSIBLE BYPASS OF TRANSFORMATION  -----
      A(NR+1,NR)=ZERO
      IF(S.LE.ZERO) GO TO 250
      S=S+B*B
      SGN=ONE
      IF(B.GE.ZERO) GO TO 160
      SGN=-ONE
  160 SQRTS=SQRT(S)
      D=SGN/(SQRTS+SQRTS)
      TEMP=SQRT(PT5+B*D)
      W(NR)=TEMP
      A(NR+1,NR)=TEMP
      D=D/TEMP
      B=-SGN*SQRTS
C        -----  D IS FACTOR OF PROPORTIONALITY.            -----
C        -----  COMPUTE AND SAVE W VECTOR
      DO 170 I=NR,N2
      TEMP=D*A(I+2,NR)
      W(I+1)=TEMP
  170 A(I+2,NR)=TEMP
C        -----  PREMULTIPLY VECTOR-W-BY MATRIX -A- TO OBTAIN *****
C        -----  -P- VECTOR.  SIMULTANEOUSLY ACCUMULATE DOT   -----
C        -----  PRODUCT WP,(THE SCALAR K).                   -----
      WTAW=ZERO
      DO 220 I=NR,N1
      SUM=ZERO
      DO 180 J=NR,I
  180 SUM=SUM+A(I+1,J+1)*W(J)
      I1=I+1
      IF((N1-I1).LT.0) GO TO 210
      DO 200 J=I1,N1
  200 SUM=SUM+A(J+1,I+1)*W(J)
  210 P(I)=SUM
  220 WTAW=WTAW+SUM*W(I)
C        -----  P VECTOR AND SCALAR K NOW STORED. NEXT     -----
C        ----- COMPUTE -Q- VECTOR                          -----
      DO 230 I=NR,N1
  230 Q(I)=P(I)-WTAW*W(I)
C        -----  NOW FORM PAP MATRIX, REQUIRED PART         -----
      DO 240 J=NR,N1
      QJ=Q(J)
      WJ=W(J)
      DO 240 I=J,N1
  240 A(I+1,J+1)=A(I+1,J+1)-TWO*(W(I)*QJ+WJ*Q(I))
  250 BETA(NR)=B
      BETASQ(NR)=B*B
  260 GAMMA(NR+1)=A(NR+1,NR+1)
  270 B=A(N,N-1)
      BETA(N-1)=B
      BETASQ(N-1)=B*B
      GAMMA(N)=A(N,N)
  280 BETASQ(N)=ZERO
C        -----  ADJOIN AN IDENTITY MATRIX TO BE            -----
C        -----  POSTMULTIPLIED BY ROTATIONS                -----
      DO 300 I=1,N
      DO 290 J=1,N
  290 VEC(I,J)=ZERO
  300 VEC(I,I)=ONE
      M=N
      SUM=ZERO
      NPAS=1
      GO TO 400
  310 SUM=SUM+SHIFT
      COSA=ONE
      G=GAMMA(1)-SHIFT
      PP=G
      PPBS=PP*PP+BETASQ(1)
      PPBR=SQRT(PPBS)
      DO 370 J=1,M
      COSAP=COSA
      IF(PPBS.NE.ZERO) GO TO 320
      SINA=ZERO
      SINA2=ZERO
      COSA=ONE
      GO TO 350
  320 SINA=BETA(J)/PPBR
      SINA2=BETASQ(J)/PPBS
      COSA=PP/PPBR
C        -----  POSTMULTIPLY IDENTITY BY P-TRANSPOSE       -----
      NT=J+NPAS
      IF(NT.LE.N) GO TO 330
      NT = N
  330 DO 340 I=1,NT
      TEMP=COSA*VEC(I,J)+SINA*VEC(I,J+1)
      VEC(I,J+1)=-SINA*VEC(I,J)+COSA*VEC(I,J+1)
  340 VEC(I,J)=TEMP
  350 DIA=GAMMA(J+1)-SHIFT
      U=SINA2*(G+DIA)
      GAMMA(J)=G+U
      G=DIA-U
      PP=DIA*COSA-SINA*COSAP*BETA(J)
      IF(J.NE.M) GO TO 360
      BETA(J)=SINA*PP
      BETASQ(J)=SINA2*PP*PP
      GO TO 380
  360 PPBS=PP*PP+BETASQ(J+1)
      PPBR=SQRT(PPBS)
      BETA(J)=SINA*PPBR
  370 BETASQ(J)=SINA2*PPBS
  380 GAMMA(M+1)=G
C        -----  TEST FOR CONVERGENCE OF LAST DIAGONAL ELEMENT  -----
      NPAS=NPAS+1
      IF(BETASQ(M).GT.RHOSQ) GO TO 410
  390 EIG(M+1)=GAMMA(M+1)+SUM
  400 BETA(M)=ZERO
      BETASQ(M)=ZERO
      M=M-1
      IF(M.EQ.0) GO TO 430
      IF(BETASQ(M).LE.RHOSQ) GO TO 390
C        -----  TAKE ROOT OF CORNER 2 BY 2 NEAREST TO      -----
C        -----  LOWER DIAGONAL IN VALUE AS ESTIMATE OF     -----
C        -----  EIGENVALUE TO USE FOR SHIFT                -----
  410 A2=GAMMA(M+1)
      R2=PT5*A2
      R1=PT5*GAMMA(M)
      R12=R1+R2
      DIF=R1-R2
      TEMP=SQRT(DIF*DIF+BETASQ(M))
      R1=R12+TEMP
      R2=R12-TEMP
      DIF=ABS(A2-R1)-ABS(A2-R2)
      IF(DIF.LT.ZERO) GO TO 420
      SHIFT=R2
      GO TO 310
  420 SHIFT=R1
      GO TO 310
  430 EIG(1)=GAMMA(1)+SUM
C        -----  INITIALIZE AUXILIARY TABLES REQUIRED FOR   -----
C        -----  REARANGING THE VECTORS                     -----
      DO 440 J=1,N
      IPOSV(J)=J
      IVPOS(J)=J
  440 IORD(J) = J
C        ----- USE A TRANSPOSITON SORT TO ORDER THE EIGENVALUES  -----
      M=N
      GO TO 470
  450 DO 460 J=1,M
      IF(EIG(J).LE.EIG(J+1)) GO TO 460
      TEMP=EIG(J)
      EIG(J)=EIG(J+1)
      EIG(J+1)=TEMP
      ITEMP=IORD(J)
      IORD(J)=IORD(J+1)
      IORD(J+1)=ITEMP
  460 CONTINUE
  470 M=M-1
      IF(M.NE.0) GO TO 450
      IF(N1.EQ.0) GO TO 500
      DO 490 L=1,N1
      NV=IORD(L)
      NP=IPOSV(NV)
      IF(NP.EQ.L) GO TO 490
      LV=IVPOS(L)
      IVPOS(NP)=LV
      IPOSV(LV)=NP
      DO 480 I=1,N
      TEMP=VEC(I,L)
      VEC(I,L)=VEC(I,NP)
  480 VEC(I,NP) = TEMP
  490 CONTINUE
  500 CONTINUE
C        -----  BACK TRANSFORM THE VECTORS OF THE TRIPLE   -----
C        -----  DIAGONAL MATRIX                            -----
      DO 550 NRR=1,N
      K=N1
  510 K=K-1
      IF(K.LE.0) GO TO 540
      SUM=ZERO
      DO 520 I=K,N1
  520 SUM=SUM+VEC(I+1,NRR)*A(I+1,K)
      SUM=SUM+SUM
      DO 530 I=K,N1
  530 VEC(I+1,NRR)=VEC(I+1,NRR)-SUM*A(I+1,K)
      GO TO 510
  540 CONTINUE
  550 CONTINUE
  560 CONTINUE
      RETURN
      END
C*MODULE TRUDGE  *DECK TRUATO
      SUBROUTINE TRUATO
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXSH=5000)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTNAL/ BOND(MXATM),ALPHA(MXATM),BETA(MXATM),SIGN(MXATM),
     *                ICONX(MXATM),IATCON(3,MXATM),INATOM(MXATM),
     *                NATIN,IUNT
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /TRANSF/ XSMAL,YSMAL,ZSMAL,XNEW,YNEW,ZNEW,XP,YP,ZP
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      DATA TOL /1.0D-10/
      DATA ZERO/0.0D+00/
      DATA RNONE/8HNONE    /
C
C     READ IN UNIQUE CENTERS AND ATOMIC BASIS SETS GROUPED
C     IN SHELLS
C     GENERATE ALL NEW CENTERS
C     SET TABLE ( CENTERS VERSUS TRANSFORMATIONS )
C     SET TABLE ( SHELLS VERSUS TRANSFORMATIONS )
C
C     APPLY CONSTRAINTS TO COORDINATES IN VB2000
      IF(VBTYP.NE.RNONE) CALL TRUCON
C
      NATMAX=MXATM
      NAT=0
C
C     ----- RESTORE UNIQUE CENTER -----
C
      DO 3000 NATI=1,NATIN
      CALL TRURST(X,Y,Z,C,NAT,NATMAX,NATI)
      NAT=NAT+1
      C(1,NAT)=X
      C(2,NAT)=Y
      C(3,NAT)=Z
C
C     GENERATE EQUIVALENT CENTERS
C
      CALL LOCAL(X,Y,Z,XS,YS,ZS)
      XSMAL=XS
      YSMAL=YS
      ZSMAL=ZS
      DO 2300 IT=1,NT
      IF(IT.EQ.1) GO TO 2300
      NN=9*(IT-1)
      CALL TRANS(NN)
      CALL ROT
      DO 2100 IAT=1,NAT
      TEST=(XP-C(1,IAT))**2+(YP-C(2,IAT))**2+(ZP-C(3,IAT))**2
      IF(TEST.LE.TOL) GO TO 2300
 2100 CONTINUE
      NUAT=NAT
      NAT=NAT+1
      C(1,NAT)=C(1,NUAT)
      C(2,NAT)=C(2,NUAT)
      C(3,NAT)=C(3,NUAT)
      C(1,NUAT)=XP
      C(2,NUAT)=YP
      C(3,NUAT)=ZP
 2300 CONTINUE
 3000 CONTINUE
      DO 3005 IAT=1,NAT
      DO 3005 I=1,3
      DUM=C(I,IAT)
      IF( ABS(DUM).LT.TOL) DUM=ZERO
 3005 C(I,IAT)=DUM
C
      RETURN
      END
C*MODULE TRUDGE  *DECK TRURST
      SUBROUTINE TRURST(X,Y,Z,C,NAT,NATMAX,NATIN)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXATM4=MXATM+4)
C
      DIMENSION C(3,1),CC(3,MXATM4),XYZ(3),XYZ0(3),T(3,3)
C
      COMMON /INTNAL/ BOND(MXATM),ALPHA(MXATM),BETA(MXATM),SIGN(MXATM),
     *                ICONX(MXATM),IATCON(3,MXATM),INATOM(MXATM),
     *                NATIN1,IUNT
C
      DATA ZERO,ONE,TWO /0.0D+00,1.0D+00,2.0D+00/
      DATA UNIT/0.52917724924D+00/
C
C     ----- RESTORE INTERNAL COORDINATES.
C           CONVERT TO CARTESIAN COORDINATES. -----
C
      JMAX=NATMAX+4
      DO 1300 J=1,JMAX
      DO 1300 I=1,3
 1300 CC(I,J)=ZERO
      CC(1,NATMAX+2)=ONE
      CC(2,NATMAX+3)=ONE
      CC(3,NATMAX+4)=ONE
      IF(NAT.EQ.0) GO TO 1500
      DO 1400 J=1,NAT
      DO 1400 I=1,3
 1400 CC(I,J)=C(I,J)
 1500 CONTINUE
C
C     ----- CALCULATE CARTESIAN COORDINATES -----
C
      KCONX=ICONX(NATIN)
      IAT1=IATCON(1,NATIN)
      IAT2=IATCON(2,NATIN)
      RAB=RIJ(CC,IAT1,IAT2)
      DO 1800 I=1,3
 1800 T(I,1)=(CC(I,IAT2)-CC(I,IAT1))/RAB
      IF(KCONX.GT.1) GO TO 2000
C
C     ----- LINEAR CONNECTION -----
C
      DUM=-SIGN(NATIN)*BOND(NATIN)/UNIT
      DO 1900 I=1,3
 1900 XYZ(I)=DUM*T(I,1)
      GO TO 3000
C
C     ----- TRI-ATOMIC CONNECTION -----
C
 2000 CONTINUE
      IAT3=IATCON(3,NATIN)
      RAC=RIJ(CC,IAT1,IAT3)
      RBC=RIJ(CC,IAT2,IAT3)
C
C     ----- DEFINE LOCAL FRAME -----
C
      DO 2100 I=1,3
 2100 T(I,2)=CC(I,IAT3)-CC(I,IAT1)
      DOT=T(1,1)*T(1,2)+T(2,1)*T(2,2)+T(3,1)*T(3,2)
      DO 2200 I=1,3
 2200 T(I,2)=T(I,2)-DOT*T(I,1)
      DOT=T(1,2)*T(1,2)+T(2,2)*T(2,2)+T(3,2)*T(3,2)
      DOT= SQRT(DOT)
      DO 2300 I=1,3
 2300 T(I,2)=T(I,2)/DOT
      T(1,3)=T(2,1)*T(3,2)-T(3,1)*T(2,2)
      T(2,3)=T(3,1)*T(1,2)-T(1,1)*T(3,2)
      T(3,3)=T(1,1)*T(2,2)-T(2,1)*T(1,2)
C
C     ----- DEFINE POLAR COODINATES -----
C
      ALPH=ALPHA(NATIN)
      PHI=ALPH
      GO TO (2400,2500,2600,2700,2500,2700),KCONX
 2400 CALL ABRT
      STOP
 2500 THETA=ZERO
      GO TO 2800
 2600 BET=BETA(NATIN)
      GAM  =ACOS((RAB*RAB+RAC*RAC-RBC*RBC)/(TWO*RAB*RAC))
      THETA=ACOS((COS(BET)-COS(ALPH)*COS(GAM))/(SIN(ALPH)*SIN(GAM)))
      GO TO 2850
 2700 BET=BETA(NATIN)
      THETA=BET
      GO TO 2850
 2800 CONTINUE
      PHI=PHI*SIGN(NATIN)
 2850 CONTINUE
      THETA=THETA*SIGN(NATIN)
      R=BOND(NATIN)/UNIT
C
C     ----- GET CARTESIAN COORDINATES IN LOCAL FRAME -----
C
      XYZ0(1)=R*COS(PHI)
      XYZ0(2)=R*SIN(PHI)*COS(THETA)
      XYZ0(3)=R*SIN(PHI)*SIN(THETA)
C
C     ----- GET CARTESIAN COORDINATES IN MASTER FRAME -----
C
      DO 2900 I=1,3
      XYZ(I)=ZERO
      DO 2900 J=1,3
      XYZ(I)=XYZ(I)+T(I,J)*XYZ0(J)
 2900 CONTINUE
 3000 CONTINUE
      X=XYZ(1)+CC(1,IAT1)
      Y=XYZ(2)+CC(2,IAT1)
      Z=XYZ(3)+CC(3,IAT1)
      RETURN
      END
C*MODULE TRUDGE  *DECK TRUHAT
      SUBROUTINE TRUHAT
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INTNAL/
     +  BOALBE(MXATM,3),SIGN(MXATM),ICONX(MXATM),
     1 IATCON(3,MXATM),INATOM(MXATM),NATIN,IUNT
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DIMENSION ZOALBE(MXATM,3)
C
      DATA PI2,DEGREE /6.28318530717958D+00,360.0D+00/
C
C     ... PRINT CURRENT INTERNAL COORDINATES
C
      DO 10 I=1,NATIN
      ZOALBE(I,1)=BOALBE(I,1)
      DO 10 J=2,3
   10 ZOALBE(I,J)=BOALBE(I,J)*DEGREE/PI2
C
      IF (MASWRK) THEN
         WRITE(IW,9983)
         WRITE(IW,9982)
         DO 20 I=1,NATIN
         WRITE(IW,9981) I,ZOALBE(I,1),ZOALBE(I,2),ZOALBE(I,3),
     1      SIGN(I),ICONX(I),
     2      IATCON(1,I),IATCON(2,I),IATCON(3,I),
     3      INATOM(I)
   20    CONTINUE
         WRITE(IW,9983)
      END IF
C
 9982 FORMAT(
     1 ' INTERNAL COORDINATES (ANGSTROMS/DEGREE)'
     2/'  I   BOND     ALPHA      BETA    SIGN ICONX   IATCON',
     3 '    INATOM')
 9981 FORMAT(I3,F10.7,2F10.5,F5.1,I6,3I4,I7)
 9983 FORMAT(1X,62('-'))
      RETURN
      END
