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