      SUBROUTINE DRC(STARTV, STARTK)
      implicit double precision (A-H,O-Z)
      include 'MOP6FILE'
      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, GETNAM*80
      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
c      open(unit=MIOM,status='SCRATCH') ! what to hell it is used for ?! J/M
      if (INDEX(KEYWRD,' PREC').NE.0) then
         ACCU=0.25D0
      else
         ACCU=1.D0
      end if
      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(MIOB,'(A)')       ' INITIAL VELOCITY IN DRC'
         write(MIOB,'(3F13.5)')  (STARTV(I),I=1,NUMAT*3)
      end if
      LET=(INDEX(KEYWRD,' GEO-OK').NE.0.OR.VELRED)
      if (INDEX(KEYWRD,' SYM').NE.0) then
         write(MIOB,*) '  SYMMETRY SPECifIED, BUT CANNOT BE USED IN DRC'
         NDEP=0
      end if
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
      end if
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
         end if
         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
      end if
      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(MIOB,'(//10X,
     %     '' DAMPING FACTOR FOR KINETIC ENERGY ='',F12.6)') HALF
      else if (INDEX(KEYWRD,'DRC').EQ.0) then
         HALF=0.D0
      else
         HALF=1.D6
      end if
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(MIOB,'(//10X,
     %    '' EXCESS KINETIC ENERGY ENTERED INTO SYSTEM ='',F12.6)')
     %   ADDONK
      else
         ADDONK=0.D0
      end if
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=MIOC,file=GETNAM('FOR009'),
     %               status='UNKNOWN',form='FORMATTED')
         rewind MIOC
         open(unit=MIOD,file=GETNAM('FOR010'),
     %               status='UNKNOWN',form='UNFORMATTED')
         rewind MIOD
         read(MIOC,'(A80)') ALPHA
         read(MIOC,'(3F19.13)') (XPARAM(I),I=1,NVAR)
         read(MIOC,'(A80)') ALPHA
         read(MIOC,'(3F19.3)') (VELO0(I),I=1,NVAR)
         read(MIOC,'(A80)') ALPHA
         read(MIOC,*) (GRAD(I),I=1,NVAR)
         read(MIOC,*) (GROLD(I),I=1,NVAR)
         read(MIOC,*) (GROLD2(I),I=1,NVAR)
         read(MIOC,*) ETOT,ESCF,EKIN,DELOLD,DELTAT,DLOLD2,ILOOP,
     1GNORM,LETOT,ELOST1,GTOT
         write(MIOB,'(//10X,''CALCULATION RESTARTED, CURRENT'',
     1                         '' KINETIC ENERGY='',F10.5,//)')  EKIN
         go to 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
            end if
            if (K.LT.0) then
               K=-K
               ONE=-1.D0
            else
               ONE=1.D0
            end if
            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)
     1                VELO0(I) = VELO0(I)+VELO1(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(MIOB,'(A,F10.3,A,/,A)')
     %             ' BY DEFAULT, ONE QUANTUM OF ENERGY,'
     1           //' 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(MIOB,'(A,F7.2,A)')
     %              ' THIS REPRESENTS AN ENERGY OF',ADDONK,' KCALS/MOLE'
               else if (ABS(HALF).GT.1.D-3) then
                  write(MIOB,'(A,F9.2,A)')
     %              ' THE VIBRATIONAL FREQUENCY (',STARTK(K),
     1              'CM(-1)) IS TOO SMALL',' FOR ONE QUANTUM TO BE USED'
                  write(MIOB,'(A)')
     1             ' INSTEAD 0.3KCAL/MOLE WILL BE USED TO START THE IRC'
                  ADDONK=0.3D0
               else
                  ADDONK=0.3D0
               end if
            end if
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(MIOB,'(A)')  ' SYSTEM IS APPARENTLY NOT MOVING!'
               return
            end if
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)
     1               ADDONK=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)
     1               ADDONK=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
            end if
         end if
      end if
  100 continue
      IUPPER=ILOOP+4999
      ILP=ILOOP
      ONE=0.D0
      if (INDEX(KEYWRD,'RESTART').NE.0.AND.INDEX(KEYWRD,'IRC=').EQ.0)
     1       ONE=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
         end if
         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
         end if
         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
            end if
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
            end if
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
               else if (ISKIN.EQ.0) then
C
C  IT IS A DRC AND KINETIC NOT USED, SO REMOVE EXTRA KINETIC ENERGY
C
                  ETOT=ETOT-ADDONK
               end if
            end if
            LETOT=.TRUE.
         end if
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
         end if
         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)
         end if
         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*
     1               MIN(2.D0, (5.D-5*ACCU/
     2               (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(MIOB,'(//,'' IRC CALCULATION COMPLETE '')')
               return
            end if
         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(MIOB,'(//,A)')
     %            ' GRADIENT CONSTANT AND SMALL -- ASSUME'
     1          //' ALL MOTION STOPPED'
               return
            end if
            DELTAT=MIN(DELTAT,2.D-15)
************************************************************************
*
*         TESTING CODE - REMOVE BEFORE FINAL VERSION ASSEMBLED
C#          (ILOOP/400)*400.EQ.ILOOP)DELTAT=-DELTAT
*
************************************************************************
         end if
         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,
     1                  ELOST1,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,
     1                  EKIN,DUMMY,ETOT,VELO0,NVAR)
         end if
         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)
     1                       close (unit=MIOC,status='DELETE')
            if (NUMCAL.NE.1) close (unit=MIOC)
c
            open(unit=MIOC,file=GETNAM('FOR009'),
     %                  status='NEW',form='FORMATTED')
            rewind MIOC
c
            open(unit=MIOD,file=GETNAM('FOR010'),
     %                  status='UNKNOWN',form='UNFORMATTED')
            rewind MIOD
c
            write(MIOC,'(A)')
     %        ' CARTESIAN GEOMETRY PARAMETERS IN ANGSTROMS'
            write(MIOC,'(3F19.13)') (XPARAM(I),I=1,NVAR)
            write(MIOC,'(A)')
     %        ' VELOCITY FOR EACH CARTESIAN COORDINATE, IN CM/SEC'
            write(MIOC,'(3F19.3)') (VELO0(I),I=1,NVAR)
            write(MIOC,'(A)')
     %        ' FIRST, SECOND, AND THIRD-ORDER GRADIENTS, ETC'
            write(MIOC,*) (GRAD(I)  ,I=1,NVAR)
            write(MIOC,*) (GROLD(I) ,I=1,NVAR)
            write(MIOC,*) (GROLD2(I),I=1,NVAR)
            I=ILOOP+1
            write(MIOC,*) ETOT,ESCF,EKIN,DELOLD,DELTAT,DLOLD2,I,
     %                    GNORM,LETOT,ELOST1,GTOT
            ESCF=-1.D9
            call PRTDRC(ESCF,DELTAT,XPARAM,GEOREF,
     1                  EKIN,ELOST,ETOT,VELO0,NVAR)
            LINEAR=(NORBS*(NORBS+1))/2
            write(MIOD) (PA(I),I=1,LINEAR)
            if (NALPHA.NE.0) write(MIOD) (PB(I),I=1,LINEAR)
            write(MIOB,'(//10X,
     1       '' RUNNING OUT OF TIME, RESTART FILE WRITTEN'')')
            write(MIOB,'(A)')
     1       ' GEOMETRY AND VELOCITY ARE IN RESTART FILE IN ASCII'
            return
         end if
  190 continue
      end
