      SUBROUTINE DRC(STARTV, STARTK)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION STARTV(*), STARTK(*)
************************************************************************
*                                                                      *
*    DRC IS DESIGNED TO FOLLOW A REACTION PATH FROM THE TRANSITION     *
*    STATE.  TWO MODES ARE SUPPORTED, FIRST: GAS PHASE:- AS THE SYSTEM *
*    MOVES FROM THE T/S THE MOMENTUM OF THE ATOMS IS STORED AND THE    *
*    POSITION OF THE ATOMS IS RELATED TO THE OLD POSITION BY (A) THE   *
*    CURRENT VELOCITY OF THE ATOM, AND (B) THE FORCES ACTING ON THAT   *
*    ATOM.  THE SECOND MODE IS CONDENSED PHASE, IN WHICH THE ATOMS MOVE*
*    IN RESPONSE TO THE FORCES ACTING ON THEM. I.E. INFINITELY DAMPED  *
*                                                                      *
************************************************************************
      INCLUDE 'SIZES'
      COMMON /KEYWRD/ KEYWRD
      COMMON /TIMDMP/ TLEFT, TDUMP
      COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
      COMMON /GRADNT/ GRAD(MAXPAR),GNORM
      COMMON /NUMCAL/ NUMCAL
      COMMON /GEOSYM/ NDEP,LOCPAR(MAXPAR),IDEPFN(MAXPAR),LOCDEP(MAXPAR)
      COMMON /GEOM  / GEO(3,NUMATM)
      COMMON /ATMASS/ ATMASS(NUMATM)
      COMMON /GEOVAR/ NVAR, LOC(2,MAXPAR), IDUMY, XPARAM(MAXPAR)
      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,XRACT
      COMMON /DRCCOM/ MCOPRT(2,MAXPAR), NCOPRT, PRTMAX
      CHARACTER KEYWRD*241
      DIMENSION VELREF(MAXPAR), VELO0(MAXPAR), VELO1(MAXPAR),
     1VELO2(MAXPAR), VELO3(MAXPAR), GERROR(MAXPAR),
     2COORD(3,NUMATM), GROLD2(MAXPAR), PAST10(10),
     3GROLD(MAXPAR), PAROLD(MAXPAR), GEOREF(3,NUMATM),
     4 SQRTMS(MAXPAR)
      LOGICAL INT, ADDK, LETOT, LET, VELRED,PRTMAX, IRCDRC
      DATA VELO0/MAXPAR*0.D0/, INT/.TRUE./
      DATA ADDK/.TRUE./
      TNOW=SECOND()
      OLDTIM=SECOND()
      DELOLD=10.D0
      GTOT=0.D0
      OPEN(UNIT=7,STATUS='SCRATCH')
      IF(INDEX(KEYWRD,' PREC').NE.0)THEN
         ACCU=0.25D0
      ELSE
         ACCU=1.D0
      ENDIF
      LPOINT=0
      GNLIM=1.D0
      PAST10(5)=100.D0
      I=INDEX(KEYWRD,'GNORM')
      IF(I.NE.0)GNLIM=READA(KEYWRD,I)
      VELRED=(INDEX(KEYWRD,'VELO').NE.0)
      IF(DOT(STARTV,STARTV,3*NUMAT).GT.0.001D0)THEN
C
C     PRINT OUT INITIAL VELOCITIES
C
         WRITE(6,'(A)')' INITIAL VELOCITY IN DRC'
         WRITE(6,'(3F13.5)')(STARTV(I),I=1,NUMAT*3)
      ENDIF
      LET=(INDEX(KEYWRD,' GEO-OK').NE.0.OR.VELRED)
      IF(INDEX(KEYWRD,' SYM').NE.0)THEN
         WRITE(6,*)'  SYMMETRY SPECIFIED, BUT CANNOT BE USED IN DRC'
         NDEP=0
      ENDIF
C
C      CONVERT TO CARTESIAN COORDINATES, IF NOT ALREADY DONE.
C
      IF(INDEX(KEYWRD,' XYZ').EQ.0)THEN
         NA(1)=0
         CALL GMETRY(GEO,COORD)
         L=0
C
         DO 10 I=1,NUMAT
            LABELS(I)=NAT(I)
            SUM=SQRT(ATMASS(NAT(I)))
            SQRTMS(L+1)=SUM
            SQRTMS(L+2)=SUM
            SQRTMS(L+3)=SUM
            L=L+3
   10    CONTINUE
         DO 30 J=1,3
            DO 20 I=1,NUMAT
               GEO(J,I)=COORD(J,I)
               COORD(J,I)=0.0D0
   20       CONTINUE
   30    CONTINUE
C
         NA(1)=99
      ENDIF
C
C  TRANSFER COORDINATES TO XPARAM AND LOC
C
      IF(INDEX(KEYWRD,' DRC').NE.0)THEN
         PRTMAX=(LOC(1,1).EQ.1)
         IF(PRTMAX)THEN
            J=1
         ELSE
            J=0
         ENDIF
         NVAR=NVAR-J
         DO 40 I=1,NVAR
            MCOPRT(1,I)=LOC(1,I+J)
   40    MCOPRT(2,I)=LOC(2,I+J)
         IF(LOC(1,1).EQ.0)NVAR=0
         NCOPRT=NVAR
      ELSE
         NCOPRT=0
      ENDIF
      L=0
      DO 50 I=1,NUMAT
         LOC(1,L+1)=I
         LOC(2,L+1)=1
         GEOREF(1,I)=GEO(1,I)
         XPARAM(L+1)=GEO(1,I)
C
         LOC(1,L+2)=I
         LOC(2,L+2)=2
         GEOREF(2,I)=GEO(2,I)
         XPARAM(L+2)=GEO(2,I)
C
         LOC(1,L+3)=I
         LOC(2,L+3)=3
         GEOREF(3,I)=GEO(3,I)
         XPARAM(L+3)=GEO(3,I)
C
         L=L+3
   50 CONTINUE
      NVAR=NUMAT*3
C
C DETERMINE DAMPING FACTOR
C
      IRCDRC=(INDEX(KEYWRD,'IRC=').NE.0)
      IF(INDEX(KEYWRD,'DRC=').NE.0) THEN
         HALF=READA(KEYWRD,INDEX(KEYWRD,'DRC='))
         WRITE(6,'(//10X,'' DAMPING FACTOR FOR KINETIC ENERGY ='',F12.6)
     1')HALF
      ELSEIF (INDEX(KEYWRD,'DRC').EQ.0) THEN
         HALF=0.D0
      ELSE
         HALF=1.D6
      ENDIF
C
C  LETOT IS TRUE IF CORRECTIONS ARE NOT TO BE MADE PART WAY INTO
C        THE CALCULATION
C
C  USAGE OF LETOT:
C (1) WHILE LETOT IS FALSE, NO DAMPING WILL BE DONE
C (2) WHEN LETOT IS TURNED TRUE,
C     IF AN IRC, THEN ETOT IS RESET SO THE ERROR IS ZERO.
C     IF A  DRC, EXCESS KINETIC ENERGY USED TO START THE RUN IS REMOVED.
C
      LETOT=(INDEX(KEYWRD,'IRC=').EQ.0 .AND. .NOT. LET)
      HALF=SIGN(MAX(0.000001D0,ABS(HALF)),HALF)
C
C DETERMINE EXCESS KINETIC ENERGY
C
      ISKIN=0
      IF(INDEX(KEYWRD,'KINE').NE.0) THEN
         ISKIN=1
         ADDONK=READA(KEYWRD,INDEX(KEYWRD,'KINE'))
         WRITE(6,'(//10X,'' EXCESS KINETIC ENERGY ENTERED INTO SYSTEM ='
     1',F12.6)')ADDONK
      ELSE
         ADDONK=0.D0
      ENDIF
C
C   LOOP OVER TIME-INTERVALS OF DELTAT SECOND
C
      DELTAT=1.D-16
      QUADR=1.D0
      TOTIME=0.D0
      ONCE=0.D0
      ETOT=0.D0
      ESCF=0.D0
      CONST=1.D0
      IF( INDEX(KEYWRD,'RESTART').NE.0.AND.INDEX(KEYWRD,'IRC=').EQ.0)
     1THEN
C
C  RESTART FROM A PREVIOUS RUN
C
         OPEN(UNIT=9,FILE='FOR009',STATUS='UNKNOWN',FORM='FORMATTED')
         REWIND 9
         OPEN(UNIT=10,FILE='FOR010',STATUS='UNKNOWN',FORM='UNFORMATTED')
         REWIND 10
         READ(9,'(A)')ALPHA
         READ(9,'(3F19.13)')(XPARAM(I),I=1,NVAR)
         READ(9,'(A)')ALPHA
         READ(9,'(3F19.3)')(VELO0(I),I=1,NVAR)
         READ(9,'(A)')ALPHA
         READ(9,*)(GRAD(I),I=1,NVAR)
         READ(9,*)(GROLD(I),I=1,NVAR)
         READ(9,*)(GROLD2(I),I=1,NVAR)
         READ(9,*)ETOT,ESCF,EKIN,DELOLD,DELTAT,DLOLD2,ILOOP,
     1GNORM,LETOT,ELOST1,GTOT
         WRITE(6,'(//10X,''CALCULATION RESTARTED, CURRENT'',
     1'' KINETIC ENERGY='',F10.5,//)')EKIN
         GOTO 100
      ELSE
C                         NOT A RESTART
         ILOOP=1
         IF(INDEX(KEYWRD,'IRC=').NE.0.OR.VELRED)THEN
C
C  GET HOLD OF VELOCITY VECTOR
C
            IF(INDEX(KEYWRD,'IRC=').NE.0)THEN
               K=READA(KEYWRD,INDEX(KEYWRD,'IRC='))
            ELSE
               K=1
            ENDIF
            IF(K.LT.0)THEN
               K=-K
               ONE=-1.D0
            ELSE
               ONE=1.D0
            ENDIF
            KL=(K-1)*NVAR
            SUMM=0.D0
            VELO1(1)=0
            VELO1(2)=0
            VELO1(3)=0
            SUMMAS=0.D0
            I=0
            DO 60 II=1,NUMAT
               AMS=ATMASS(II)
               SUMMAS=SUMMAS+AMS
               VELO0(I+1)=STARTV(KL+I+1)*ONE
               VELREF(I+1)=VELO0(I+1)
               VELO1(1)=VELO1(1)+VELO0(I+1)*AMS
C
               VELO0(I+2)=STARTV(KL+I+2)*ONE
               VELREF(I+2)=VELO0(I+2)
               VELO1(2)=VELO1(2)+VELO0(I+2)*AMS
C
               VELO0(I+3)=STARTV(KL+I+3)*ONE
               VELREF(I+3)=VELO0(I+3)
               VELO1(3)=VELO1(3)+VELO0(I+3)*AMS
C
               I=I+3
   60       CONTINUE
C$DOIT ASIS
            DO 70 I=1,3
   70       VELO1(I)=-VELO1(I)/SUMMAS
            I=0
C$DOIT ASIS
            DO 80 II=1,NUMAT
               AMS=ATMASS(II)
C$DOIT ASIS
               DO 80 I1=1,3
                  I=I+1
                  IF(ADDONK.GT.1.D-5.OR..NOT.VELRED)VELO0(I)=VELO0(I)+VE
     1LO1(I1)
   80       SUMM=SUMM+VELO0(I)**2*AMS
            IF(ADDONK.LT.1.D-5.AND.VELRED)ADDONK=0.5D0*SUMM/4.184D10
            IF(ADDONK.LT.1.D-5.AND..NOT.VELRED)THEN
               IF(ABS(HALF).GT.1.D-3.AND.STARTK(K).GT.105.D0)THEN
                  WRITE(6,'(A,F10.3,A,/,A)')' BY DEFAULT, ONE QUANTUM OF
     1 ENERGY,'//' EQUIVALENT TO',STARTK(K),' CM(-1)',
     2' WILL BE USED TO START THE DRC'
C
C    2.8585086D-3 CONVERTS CM(-1) INTO KCAL/MOLE
C
                  ADDONK=STARTK(K)*2.8585086D-3
                  WRITE(6,'(A,F7.2,A)')' THIS REPRESENTS AN ENERGY OF',A
     1DDONK,' KCALS/MOLE'
               ELSEIF(ABS(HALF).GT.1.D-3)THEN
                  WRITE(6,'(A,F9.2,A)')' THE VIBRATIONAL FREQUENCY (',ST
     1ARTK(K),'CM(-1)) IS TOO SMALL',' FOR ONE QUANTUM TO BE USED'
                  WRITE(6,'(A)')
     1' INSTEAD 0.3KCAL/MOLE WILL BE USED TO START THE IRC'
                  ADDONK=0.3D0
               ELSE
                  ADDONK=0.3D0
               ENDIF
            ENDIF
C
C   AT THIS POINT ADDONK IS IN KCAL/MOLE
C   NORMALIZE SO THAT TOTAL K.E. = ONE QUANTUM (DEFAULT) (DRC ONLY)
C                              OR 0.3KCAL/MOLE (IRC ONLY)
C                              OR ADDONK IF KINETIC=NN SUPPLIED
C
            IF(SUMM.LT.1.D-4) THEN
               WRITE(6,'(A)')' SYSTEM IS APPARENTLY NOT MOVING!'
               RETURN
            ENDIF
C
C  ADDONK IS EXCESS KINETIC ENERGY.  IF THE CALCULATION IS AN IRC,
C  THIS ENERGY MUST BE REMOVED AFTER A SHORT 'TIME'.
C
C  MAKE AN AD-HOC CORRECTION: IF ADDONK IS NON-ZERO AND HALF IS LARGER
C  THAN 0.1, MODIFY ADDONK TO REFLECT ERRORS DUE TO START-UP.
C
            IF(HALF.GT.0.1D0.AND.HALF.LT.10000.D0)
     1ADDONK=ADDONK*(1.D0+0.06972D0/HALF)
C
C  MAKE AN AD-HOC CORRECTION: IF ADDONK IS NON-ZERO AND HALF IS LESS
C  THAN -0.1, MODIFY ADDONK TO REFLECT ERRORS DUE TO START-UP.
C
            IF(HALF.LT.-0.1D0.AND.HALF.GT.-10000.D0)
     1ADDONK=ADDONK*(1.D0+0.06886D0/HALF)
            SUMM=SQRT(ADDONK/(0.5D0*SUMM/4.184D10))
            ADDK=.FALSE.
            IF(SUMM.GT.1.D-10)THEN
               DO  90 I=1,NVAR
   90          VELO0(I)=VELO0(I)*SUMM
C
C  IF IT IS A DRC, DESTROY ADDONK.  THE KINETIC ENERGY USED WILL COME
C  FROM THE VELOCITY ONLY.
C
               IF(HALF.GT.1.D-3)ADDONK=0.D0
            ENDIF
         ENDIF
      ENDIF
  100 CONTINUE
      IUPPER=ILOOP+4999
      ILP=ILOOP
      ONE=0.D0
      IF(INDEX(KEYWRD,'RESTART').NE.0.AND.INDEX(KEYWRD,'IRC=').EQ.0)
     1ONE=1.D0
      DO 190 ILOOP=ILP,IUPPER
C
C  MOVEMENT OF ATOMS WILL BE PROPORTIONAL TO THE AVERAGE VELOCITIES
C  OF THE ATOMS BEFORE AND AFTER TIME INTERVAL
C
C
C  RAPID CHANGE IN GRADIENT IMPLIES SMALL STEP SIZE FOR DELTAT
C
C   KINETIC ENERGY = 1/2 * M * V * V
C                  = 0.5 / (4.184D10) * M * V * V
C   NEW VELOCITY = OLD VELOCITY + GRADIENT * TIME / MASS
C                = KCAL/ANGSTROM*SECOND/(ATOMIC WEIGHT)
C                =4.184*10**10(ERGS)*10**8(PER CM)*DELTAT(SECONDS)
C   NEW POSITION = OLD POSITION - AVERAGE VELOCITY * TIME INTERVAL
C
C
C   ESTABLISH REFERENCE TOTAL ENERGY
C
         ERROR=(ETOT-(EKIN+ESCF))
         IF(ILOOP.GT.2)THEN
            QUADR = 1.D0+ERROR/(EKIN*CONST+0.001D0)*0.5D0
            QUADR = MIN(1.3D0,MAX(0.8D0,QUADR))
         ELSE
            QUADR=1.D0
         ENDIF
         IF((LET.OR.EKIN.GT.0.2).AND.ADDK)THEN
C
C   DUMP IN EXCESS KINETIC ENERGY
C
            ETOT=ETOT+ADDONK
            ADDK=.FALSE.
            ADDONK=0.D0
         ENDIF
         EKOLD=EKIN
C
C  CALCULATE THE DURATION OF THE NEXT STEP.
C  STEP SIZE IS THAT REQUIRED TO PRODUCE A CONSTANT CHANGE IN GEOMETRY
C
C
C  IF DAMPING IS USED, CALCULATE THE NEW TOTAL ENERGY AND
C  THE RATIO FOR REDUCING THE KINETIC ENERGY
C
         CONST=MAX(1.D-36,0.5D0**(DELTAT*1.D15/HALF))
         CONST=SQRT(CONST)
         VELVEC=0.D0
         EKIN=0.D0
         DELTA1=DELOLD+DLOLD2
         ELOST=0.D0
         DO 110 I=1,NVAR
C
C   CALCULATE COMPONENTS OF VELOCITY AS
C   V = V(0) + V'*T + V"*T*T
C   WE NEED ALL THREE TERMS, V(0), V' AND V"
C
            VELO1(I) = 1.D0/ATMASS(LOC(1,I))*GRAD(I)
            IF(ILOOP.GT.3) THEN
               VELO3(I) = 2.D0/ATMASS(LOC(1,I))*
     1(DELTA1*(GROLD(I)-GRAD(I))-DELOLD*(GROLD2(I)-GRAD(I)))/
     2(DELTA1*(DELOLD**2*1.D30)-DELOLD*(DELTA1**2*1.D30))
               VELO2(I)=1.D0/ATMASS(LOC(1,I))*
     1(GRAD(I)-GROLD(I)-0.5D0*VELO3(I)*(1.D30*DELOLD**2))/(DELOLD*1.D15)
            ELSE
               VELO2(I) = 1.D0/ATMASS(LOC(1,I))*
     1                 (GRAD(I)-GROLD(I))/(1.D15*DELOLD)
               VELO3(I)=0.D0
            ENDIF
C
C  MOVE ATOMS THROUGH DISTANCE EQUAL TO VELOCITY * DELTA-TIME, NOTE
C  VELOCITY CHANGES FROM START TO FINISH, THEREFORE AVERAGE.
C
            PAROLD(I)=XPARAM(I)
            XPARAM(I)=XPARAM(I)
     1             -1.D8*(DELTAT*VELO0(I)*ONE
     2             +0.5D0*DELTAT**2*VELO1(I)
     3             +0.16666D0*(DELTAT**2*1.D15)*DELTAT*VELO2(I)
     4             +0.0416666D0*DELTAT**2*(1.D30*DELTAT**2)*VELO3(I))
C
C   CORRECT ERRORS DUE TO CUBIC COMPONENTS IN ENERGY GRADIENT,
C   ALSO TO ADD ON EXCESS ENERGY, IF NECESSARY.
C
            VELVEC=VELVEC+VELO0(I)**2
C
C   MODIFY VELOCITY IN LIGHT OF CURRENT ENERGY GRADIENTS.
C
C   VELOCITY = OLD VELOCITY + (DELTA-T / ATOMIC MASS) * CURRENT GRADIENT
C                           + 1/2 *(DELTA-T * DELTA-T /ATOMIC MASS) *
C                             (SLOPE OF GRADIENT)
C              SLOPE OF GRADIENT = (GRAD(I)-GROLD(I))/DELOLD
C
C
C   THIS EXPRESSION IS ACCURATE TO SECOND ORDER IN TIME.
C
            VELO0(I) = VELO0(I) + DELTAT*VELO1(I) + 0.5D0*DELTAT**2*VELO
     12(I)*1.D15           + 0.166666D0*DELTAT*(1.D30*DELTAT**2)*VELO3(
     2I)
            IF(LET.OR.GNORM.GT.3.D0)THEN
               LET=.TRUE.
               ELOST=ELOST+VELO0(I)**2*ATMASS(LOC(1,I))*(1-CONST**2)
               VELO0(I)=VELO0(I)*CONST*QUADR
            ENDIF
C
C  CALCULATE KINETIC ENERGY (IN 2*ERGS AT THIS POINT)
C
            EKIN=EKIN+VELO0(I)**2*ATMASS(LOC(1,I))
  110    CONTINUE
         ONE=1.D0
         IF(LET.OR.GNORM.GT.3.D0)THEN
            IF(.NOT.LETOT) THEN
               IF(ABS(HALF).LT.1.D-3)THEN
C
C  IT IS AN IRC, SO RESET THE TOTAL ENERGY
C
                  ETOT=ESCF+ELOST1
                  ADDONK=0.D0
                  ELOST1=0.D0
                  ELOST=0.D0
               ELSEIF(ISKIN.EQ.0)THEN
C
C  IT IS A DRC AND KINETIC NOT USED, SO REMOVE EXTRA KINETIC ENERGY
C
                  ETOT=ETOT-ADDONK
               ENDIF
            ENDIF
            LETOT=.TRUE.
         ENDIF
C
C  CONVERT ENERGY INTO KCAL/MOLE
C
         EKIN=0.5*EKIN/4.184D10
C
C  IF IT IS A DAMPED DRC, MODIFY ETOT TO REFLECT LOSS OF KINETIC ENERGY
C
         IF(LETOT.AND.ABS(HALF).GT.0.00001D0)
     1ETOT=ETOT-EKIN/CONST**2+EKIN
         ELOST1=ELOST1+0.5D0*ELOST/4.184D10
C
C STORE OLD GRADIENTS FOR DELTA - VELOCITY CALCULATION
C
         DO 120 I=1,NVAR
            GROLD2(I)=GROLD(I)
            GROLD(I)=GRAD(I)
  120    GRAD(I)=0.D0
C
C   CALCULATE ENERGY AND GRADIENTS
C
         SCFOLD=ESCF
         CALL COMPFG(XPARAM,.TRUE.,ESCF,.TRUE.,GRAD,.TRUE.)
         IF(ILOOP.GT.2)THEN
            GNORM=0.D0
            DO 140 I=1,NVAR,3
               SUM=SQRT(DOT(GRAD(I),GRAD(I),3)/
     1(DOT(VELO0(I),VELO0(I),3)+1.D-20))
               DO 130 J=I,I+2
  130          GERROR(J)=GERROR(J)+GRAD(J)+VELO0(J)*SUM
  140       CONTINUE
            GNORM=SQRT(DOT(GERROR,GERROR,NVAR))
            GTOT=GNORM
         ENDIF
         GNORM=SQRT(DOT(GRAD,GRAD,NVAR))
C
C   CONVERT GRADIENTS INTO ERGS/CM
C
         DO 150 I=1,NVAR
  150    GRAD(I)=GRAD(I)*4.184D18
C
C   SPECIAL TREATMENT FOR FIRST POINT - SET "OLD" GRADIENTS EQUAL TO
C   CURRENT GRADIENTS.
C
         IF(ILOOP.EQ.1) THEN
            DO 160 I=1,NVAR
  160       GROLD(I)=GRAD(I)
         ENDIF
         DLOLD2=DELOLD
         DELOLD=DELTAT
         SUM=0.D0
         DO 170 I=1,NVAR
  170    SUM=SUM + ((GRAD(I)-GROLD(I))/4.184D18)**2
         IF(ABS(HALF).LT.0.001D0)THEN
            DELTAT= DELTAT*
     1MIN(2.D0, (5.D-5*ACCU/(ABS(ESCF+ELOST1-ETOLD)+1.D-20)))**0.25D0
            ETOLD=ESCF+ELOST1
            IF(ILOOP.GT.5.AND.SCFOLD-ESCF.LT.-1.D-3 .OR.
     1      ILOOP.GT.30.AND.SCFOLD-ESCF.LT.0.D0)  THEN
               WRITE(6,'(//,'' IRC CALCULATION COMPLETE '')')
               RETURN
            ENDIF
         ELSE
            DELTAT= DELTAT*MIN(1.05D0, 10.D0*ACCU/(SUM+1.D-4))
            DELTAT=MIN(DELTAT,3.D-15*ACCU)
            PAST10(10)=GNORM
            SUM=0.D0
            DO 180 I=1,9
               SUM=SUM+ABS(PAST10(I)-PAST10(I+1))
  180       PAST10(I)=PAST10(I+1)
            IF(SUM.LT.GNLIM)THEN
               WRITE(6,'(//,A)')' GRADIENT CONSTANT AND SMALL -- ASSUME'
     1//' ALL MOTION STOPPED'
               RETURN
            ENDIF
            DELTAT=MIN(DELTAT,2.D-15)
************************************************************************
*
*         TESTING CODE - REMOVE BEFORE FINAL VERSION ASSEMBLED
C#          (ILOOP/400)*400.EQ.ILOOP)DELTAT=-DELTAT
*
************************************************************************
         ENDIF
         DELTAT=MAX(1.D-16,DELTAT)
         IF(ABS(HALF).LT.0.00001D0)THEN
C
C   FOR THE IRC:
C
C ESCF   = POTENTIAL ENERGY
C ELOST1 = ENERGY LOST (IN DRC, THIS WOULD HAVE BEEN THE KINETIC ENERGY)
C ETOT   = COMPUTED TOTAL ENERGY = STARTING POTENTIAL ENERGY
C
C   IN DRCOUT  'TOTAL' = ESCF + ELOST1
C              'ERROR' = ESCF + ELOST1 - ETOT
C
            CALL PRTDRC(ESCF,DELTAT,XPARAM,GEOREF,
     1ELOST1,GTOT,ETOT,VELO0,NVAR)
         ELSE
C
C   FOR THE DRC:
C
C ESCF   = POTENTIAL ENERGY
C EKIN   = CURRENT KINETIC ENERGY
C ETOT   = COMPUTED TOTAL ENERGY = STARTING POTENTIAL ENERGY -
C          KINETIC ENERGY LOST THROUGH DAMPING, IF PRESENT.
C
C   IN DRCOUT  'TOTAL' = ESCF + EKIN
C              'ERROR' = ESCF + EKIN - ETOT
C
            CALL PRTDRC(ESCF,DELTAT,XPARAM,GEOREF,
     1EKIN,DUMMY,ETOT,VELO0,NVAR)
         ENDIF
         TNOW=SECOND()
         TCYCLE=TNOW-OLDTIM
         OLDTIM=TNOW
         TLEFT=TLEFT-TCYCLE
         IF (ILOOP.EQ.IUPPER.OR.TLEFT.LT.3*TCYCLE) THEN
            IF( INDEX(KEYWRD,'RESTART')+INDEX(KEYWRD,'ISOT').NE.0)
     1CLOSE (9,STATUS='DELETE')
            IF(NUMCAL.NE.1)CLOSE(9)
            OPEN(UNIT=9,FILE='FOR009',STATUS='NEW',FORM='FORMATTED')
            REWIND 9
            OPEN(UNIT=10,FILE='FOR010',STATUS='UNKNOWN',FORM='UNFORMATTE
     1D')
            REWIND 10
            WRITE(9,'(A)')' CARTESIAN GEOMETRY PARAMETERS IN ANGSTROMS'
            WRITE(9,'(3F19.13)')(XPARAM(I),I=1,NVAR)
            WRITE(9,'(A)')' VELOCITY FOR EACH CARTESIAN COORDINATE, IN C
     1M/SEC'
            WRITE(9,'(3F19.3)')(VELO0(I),I=1,NVAR)
            WRITE(9,'(A)')' FIRST, SECOND, AND THIRD-ORDER GRADIENTS, ET
     1C'
            WRITE(9,*)(GRAD(I),I=1,NVAR)
            WRITE(9,*)(GROLD(I),I=1,NVAR)
            WRITE(9,*)(GROLD2(I),I=1,NVAR)
            I=ILOOP+1
            WRITE(9,*)ETOT,ESCF,EKIN,DELOLD,DELTAT,DLOLD2,I,
     1GNORM,LETOT,ELOST1,GTOT
            ESCF=-1.D9
            CALL PRTDRC(ESCF,DELTAT,XPARAM,GEOREF,
     1EKIN,ELOST,ETOT,VELO0,NVAR)
            LINEAR=(NORBS*(NORBS+1))/2
            WRITE(10)(PA(I),I=1,LINEAR)
            IF(NALPHA.NE.0)WRITE(10)(PB(I),I=1,LINEAR)
            WRITE(6,'(//10X,'' RUNNING OUT OF TIME, RESTART FILE WRITTEN
     1'')')
            WRITE(6,'(A)')' GEOMETRY AND VELOCITY ARE IN RESTART FILE'
     1//' IN ASCII'
            RETURN
         ENDIF
  190 CONTINUE
      END
