mopac606_nbo/src/prtdrc.f

   1       SUBROUTINE PRTDRC(ESCF,DELTT,XPARAM,REF,EKIN,GTOT,ETOT,VELO0,NVAR)
   2       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
   3       DIMENSION XPARAM(*), VELO0(*), REF(*)
   4 *********************************************************************
   5 *
   6 *    PRTDRC PREPARES TO PRINT THE GEOMETRY ETC. FOR POINTS IN A DRC
   7 *    OR IRC
   8 *    CALCULATION.
   9 *    ON INPUT  ESCF   = HEAT OF FORMATION FOR THE CURRENT POINT
  10 *              DELTT  = CHANGE IN TIME, PREVIOUS TO CURRENT POINT
  11 *              XPARAM = CURRENT CARTESIAN GEOMETRY
  12 *              EKIN   = CURRENT KINETIC ENERGY
  13 *              GTOT   = TOTAL GRADIENT NORM IN IRC CALC'N.
  14 *              VELO0  = CURRENT VELOCITY
  15 *              NVAR   = NUMBER OF VARIABLES = 3 * NUMBER OF ATOMS.
  16 *
  17 ********************************************************************
  18       INCLUDE 'SIZES'
  19       COMMON /KEYWRD/ KEYWRD
  20       COMMON /NUMCAL/ NUMCAL
  21       COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
  22      1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
  23      2                NCLOSE,NOPEN,NDUMY,XRACT
  24       COMMON /GEOKST/ NATOMS,LABELS(NUMATM),
  25      1                NA(NUMATM),NB(NUMATM),NC(NUMATM)
  26       COMMON /DRCCOM/  MCOPRT(2,MAXPAR), NCOPRT, PARMAX
  27       COMMON /CORE  / CORE(107)
  28       COMMON /ATMASS/ ATMASS(NUMATM)
  29       COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
  30       COMMON /FMATRX/ ALLXYZ(3,MAXPAR),ALLVEL(3,MAXPAR),PARREF(MAXPAR),
  31      1XYZ3(3,MAXPAR),VEL3(3,MAXPAR), ALLGEO(3,MAXPAR), GEO3(3,MAXPAR),
  32      2 DUMMY(MAXPAR**2+1-16*MAXPAR), IDUMY(4)
  33       DIMENSION ESCF3(3),EKIN3(3), GTOT3(3), CHARGE(NUMATM), XOLD3(3),
  34      1GEO(3*NUMATM), VREF(MAXPAR), VREF0(MAXPAR), TSTEPS(100), ETOT3(3),
  35      2XTOT3(3)
  36       SAVE  REFSCF, COTYPE
  37       LOGICAL  TURN, PARMAX, LDRC, GOTURN
  38       CHARACTER*241 KEYWRD, TEXT1*3, TEXT2*2,  COTYPE(3)*2
  39       DATA ICALCN/0/
  40       DATA REFSCF/0.D0/
  41       DATA COTYPE/'BL','BA','DI'/
  42       IF (ICALCN.NE.NUMCAL) THEN
  43          ICALCN=NUMCAL
  44          DO 10 I=1,NVAR
  45    10    PARREF(I)=XPARAM(I)
  46          ETOT=ESCF+EKIN
  47          TLAST=0.D0
  48          GOTURN=.FALSE.
  49          SUM=0.D0
  50          DO 20 I=1,NVAR
  51             SUM=SUM+VELO0(I)**2
  52             VREF0(I)=VELO0(I)
  53    20    VREF(I)=VELO0(I)
  54          IONE=1
  55          LDRC=(SUM.GT.1.D0)
  56          ILOOP=1
  57          TOLD1=0.0D0
  58 C
  59 C       DETERMINE TYPE OF PRINT: TIME, ENERGY OR GEOMETRY PRIORITY
  60 C       OR PRINT ALL POINTS
  61 C
  62          STEPT=0.D0
  63          STEPH=0.D0
  64          STEPX=0.D0
  65          IF(INDEX(KEYWRD,' T-PRIO').NE.0)THEN
  66             IF(INDEX(KEYWRD,' T-PRIORITY=').NE.0)THEN
  67                STEPT=READA(KEYWRD,INDEX(KEYWRD,'T-PRIO')+5)
  68             ELSE
  69                STEPT=0.1D0
  70             ENDIF
  71             TREF=-1.D-6
  72             WRITE(6,'(/,'' TIME PRIORITY, INTERVAL ='',F4.1,
  73      1'' FEMTOSECONDS'',/)')STEPT
  74          ELSEIF(INDEX(KEYWRD,' H-PRIO').NE.0)THEN
  75             IF(INDEX(KEYWRD,' H-PRIORITY=').NE.0)THEN
  76                STEPH=READA(KEYWRD,INDEX(KEYWRD,'H-PRIO')+5)
  77             ELSE
  78                STEPH=0.1D0
  79             ENDIF
  80             WRITE(6,'(/,'' KINETIC ENERGY PRIORITY, STEP ='',F5.2,
  81      1'' KCAL/MOLE'',/)')STEPH
  82          ELSEIF(INDEX(KEYWRD,' X-PRIO').NE.0)THEN
  83             IF(INDEX(KEYWRD,' X-PRIORITY=').NE.0)THEN
  84                STEPX=READA(KEYWRD,INDEX(KEYWRD,'X-PRIO')+5)
  85             ELSE
  86                STEPX=0.05D0
  87             ENDIF
  88             WRITE(6,'(/,'' GEOMETRY PRIORITY, STEP ='',F5.2,
  89      1'' ANGSTROMS'',/)')STEPX
  90          ENDIF
  91          IF(INDEX(KEYWRD,' RESTART').NE.0.AND.INDEX(KEYWRD,'IRC=').EQ.0)
  92      1 THEN
  93             READ(9,*)(PARREF(I),I=1,NVAR)
  94             READ(9,*)(REF(I),I=1,NVAR)
  95             READ(9,*)(VREF0(I),I=1,NVAR)
  96             READ(9,*)(VREF(I),I=1,NVAR)
  97             READ(9,*)(ALLGEO(3,I),I=1,NVAR)
  98             READ(9,*)(ALLGEO(2,I),I=1,NVAR)
  99             READ(9,*)(ALLGEO(1,I),I=1,NVAR)
 100             READ(9,*)(ALLVEL(3,I),I=1,NVAR)
 101             READ(9,*)(ALLVEL(2,I),I=1,NVAR)
 102             READ(9,*)(ALLVEL(1,I),I=1,NVAR)
 103             READ(9,*)(ALLXYZ(3,I),I=1,NVAR)
 104             READ(9,*)(ALLXYZ(2,I),I=1,NVAR)
 105             READ(9,*)(ALLXYZ(1,I),I=1,NVAR)
 106             READ(9,*)ILOOP,LDRC,IONE,ETOT1,ETOT0,ESCF1,ESCF0,EKIN1,EKIN0
 107      1,TOLD2,TOLD1,GTOT1,GTOT0,XOLD2,XOLD1,XOLD0,TOTIME,JLOOP,ETOT,REFX
 108      2,XTOT1,XTOT0
 109          ENDIF
 110       ENDIF
 111       IF(ESCF.LT.-1.D8) THEN
 112          WRITE(9,*)(PARREF(I),I=1,NVAR)
 113          WRITE(9,*)(REF(I),I=1,NVAR)
 114          WRITE(9,*)(VREF0(I),I=1,NVAR)
 115          WRITE(9,*)(VREF(I),I=1,NVAR)
 116          WRITE(9,*)(ALLGEO(3,I),I=1,NVAR)
 117          WRITE(9,*)(ALLGEO(2,I),I=1,NVAR)
 118          WRITE(9,*)(ALLGEO(1,I),I=1,NVAR)
 119          WRITE(9,*)(ALLVEL(3,I),I=1,NVAR)
 120          WRITE(9,*)(ALLVEL(2,I),I=1,NVAR)
 121          WRITE(9,*)(ALLVEL(1,I),I=1,NVAR)
 122          WRITE(9,*)(ALLXYZ(3,I),I=1,NVAR)
 123          WRITE(9,*)(ALLXYZ(2,I),I=1,NVAR)
 124          WRITE(9,*)(ALLXYZ(1,I),I=1,NVAR)
 125          WRITE(9,*)ILOOP,LDRC,IONE,ETOT1,ETOT0,ESCF1,ESCF0,EKIN1,EKIN0,
 126      1TOLD2,TOLD1,GTOT1,GTOT0,XOLD2,XOLD1,XOLD0,TOTIME,JLOOP,ETOT,REFX
 127      2,XTOT1,XTOT0
 128          RETURN
 129       ENDIF
 130       CALL CHRGE(P,CHARGE)
 131       DO 30 I=1,NUMAT
 132          L=NAT(I)
 133    30 CHARGE(I)=CORE(L) - CHARGE(I)
 134       DELTAT=DELTT*1.D15
 135       NA(2)=-1
 136       CALL XYZINT(XPARAM,NUMAT,NA,NB,NC,57.29577951D0,GEO)
 137       NA(1)=99
 138       IF(ILOOP.EQ.1)THEN
 139          ETOT1=ETOT0
 140          ETOT0=ETOT
 141          ESCF1=ESCF
 142          ESCF0=ESCF
 143          EKIN1=EKIN
 144          EKIN0=EKIN
 145          DO 40 J=1,3
 146 C$DOIT VBEST
 147             DO 40 I=1,NVAR
 148                ALLGEO(J,I)=GEO(I)
 149                ALLXYZ(J,I)=XPARAM(I)
 150    40    ALLVEL(J,I)=VELO0(I)
 151       ELSE
 152          DO 50 I=1,NVAR
 153             ALLGEO(3,I)=ALLGEO(2,I)
 154             ALLGEO(2,I)=ALLGEO(1,I)
 155             ALLGEO(1,I)=GEO(I)
 156             ALLXYZ(3,I)=ALLXYZ(2,I)
 157             ALLXYZ(2,I)=ALLXYZ(1,I)
 158             ALLXYZ(1,I)=XPARAM(I)
 159             ALLVEL(3,I)=ALLVEL(2,I)
 160             ALLVEL(2,I)=ALLVEL(1,I)
 161    50    ALLVEL(1,I)=VELO0(I)
 162       ENDIF
 163 C
 164 C  FORM QUADRATIC EXPRESSION FOR POSITION AND VELOCITY W.R.T. TIME.
 165 C
 166       T1=MAX(TOLD2,0.02D0)
 167       T2=MAX(TOLD1,0.02D0)+T1
 168       DO 60 I=1,NVAR
 169          CALL QUADR(ALLGEO(3,I),ALLGEO(2,I),ALLGEO(1,I),T1,T2,
 170      1GEO3(1,I),GEO3(2,I),GEO3(3,I))
 171 C
 172 ****************************************************
 173 *                                                  *
 174 *    QUADR CALCULATES THE A, B AND C IN THE EQUNS. *
 175 *                                                  *
 176 *     A                   =   F0                   *
 177 *     A + B.X0 + C.X0**2  =   F1                   *
 178 *     A + B.X2 + C.X2**2  =   F2                   *
 179 * GIVEN THE ARGUMENT LIST (F0,F1,F2, X1,X2, A,B,C) *
 180 *                                                  *
 181 ****************************************************
 182          CALL QUADR(ALLXYZ(3,I),ALLXYZ(2,I),ALLXYZ(1,I),T1,T2,
 183      1XYZ3(1,I),XYZ3(2,I),XYZ3(3,I))
 184          CALL QUADR(ALLVEL(3,I),ALLVEL(2,I),ALLVEL(1,I),T1,T2,
 185      1VEL3(1,I),VEL3(2,I),VEL3(3,I))
 186    60 CONTINUE
 187       ETOT2=ETOT1
 188       ETOT1=ETOT0
 189       ETOT0=ETOT
 190       CALL QUADR(ETOT2,ETOT1,ETOT0,T1,T2,ETOT3(1),ETOT3(2),
 191      1ETOT3(3))
 192       EKIN2=EKIN1
 193       EKIN1=EKIN0
 194       EKIN0=EKIN
 195       CALL QUADR(EKIN2,EKIN1,EKIN0,T1,T2,EKIN3(1),EKIN3(2),
 196      1EKIN3(3))
 197       ESCF2=ESCF1
 198       ESCF1=ESCF0
 199       ESCF0=ESCF
 200       CALL QUADR(ESCF2,ESCF1,ESCF0,T1,T2,ESCF3(1),ESCF3(2),
 201      1ESCF3(3))
 202       GTOT2=GTOT1
 203       GTOT1=GTOT0
 204       GTOT0=GTOT
 205       CALL QUADR(GTOT2,GTOT1,GTOT0,T1,T2,GTOT3(1),GTOT3(2),
 206      1GTOT3(3))
 207       XTOT2=XTOT1
 208       XTOT1=XTOT0
 209       XOLD2=XOLD2+XOLD1
 210       XOLD1=XOLD0
 211 C
 212 C   CALCULATE CHANGE IN GEOMETRY
 213 C
 214       XOLD0=0.D0
 215       L=0
 216       XTOT0=0.D0
 217       SUM1=0.D0
 218       DO 80 I=1,NUMAT
 219          SUM=0.D0
 220          SUM1=0.D0
 221 C$DOIT ASIS
 222          DO 70 J=1,3
 223             L=L+1
 224             SUM1=SUM1+(ALLXYZ(1,L)-REF(L))**2
 225    70    SUM=SUM+(ALLXYZ(2,L)-ALLXYZ(1,L))**2
 226          XOLD0=XOLD0+SQRT(SUM)
 227    80 XTOT0=XTOT0+SQRT(SUM1)
 228       CALL QUADR(XTOT2,XTOT1,XTOT0,T1,T2,
 229      1XTOT3(1),XTOT3(2),XTOT3(3))
 230       CALL QUADR(XOLD2,XOLD2+XOLD1,XOLD2+XOLD1+XOLD0,T1,T2,
 231      1XOLD3(1),XOLD3(2),XOLD3(3))
 232 ***********************************************************************
 233 *   GO THROUGH THE CRITERIA FOR DECIDING WHETHER OR NOT TO PRINT THIS *
 234 *   POINT.  IF YES, THEN ALSO CALCULATE THE EXACT POINT AS A FRACTION *
 235 *   BETWEEN THE LAST POINT AND THE CURRENT POINT                      *
 236 ***********************************************************************
 237 C   NFRACT IS THE NUMBER OF POINTS TO BE PRINTED IN THE CURRENT DOMAIN
 238 ***********************************************************************
 239       IF(ILOOP.LT.3) GOTO 170
 240       FRACT=-10
 241       NFRACT=1
 242       IF(STEPH.NE.0) THEN
 243 C
 244 C   CRITERION FOR PRINTING RESULTS  IS A CHANGE IN HEAT OF FORMATION =
 245 C   -CHANGE IN KINETIC ENERGY
 246 C
 247          IF(REFSCF.EQ.0.D0) THEN
 248             I=ESCF2/STEPH
 249             REFSCF=I*STEPH
 250          ENDIF
 251          DH=ABS(ESCF1-REFSCF)
 252          IF(DH.GT.STEPH)THEN
 253             STEPH=SIGN(STEPH,ESCF1-REFSCF)
 254             NFRACT=ABS(DH/STEPH)
 255             CC=ESCF3(1)
 256             BB=ESCF3(2)
 257             AA=ESCF3(3)
 258 ************************************************
 259 * PROGRAMMERS! - BE VERY CAREFUL IF YOU CHANGE *
 260 * THIS FOLLOWING SECTION.  THERE IS NUMERICAL  *
 261 * INSTABILITY IF ABS(BB/AA) IS VERY LARGE. NEAR*
 262 * INFLECTION POINTS AA CHANGES SIGN.       JJPS*
 263 ************************************************
 264             IF(ABS(BB/AA).GT.30)THEN
 265 C
 266 C   USE LINEAR INTERPOLATION
 267 C
 268                DO 90 I=1,NFRACT
 269    90          TSTEPS(I)=-(CC-(REFSCF+I*STEPH))/BB
 270             ELSE
 271 C
 272 C  USE QUADRATIC INTERPOLATION
 273 C
 274                DO 100 I=1,NFRACT
 275                   C1=CC-(REFSCF+I*STEPH)
 276   100          TSTEPS(I)=(-BB+SIGN(SQRT(BB*BB-4.D0*(AA*C1)),BB))/(2.D0*A
 277      1A)
 278             ENDIF
 279             FRACT=-.1
 280             REFSCF=REFSCF+NFRACT*STEPH
 281          ENDIF
 282       ELSEIF(STEPT.NE.0.D0) THEN
 283 C
 284 C   CRITERION FOR PRINTING RESULTS IS A CHANGE IN TIME.
 285 C
 286          IF(ABS(TOTIME+TOLD2-TREF).GT.STEPT)THEN
 287             I=TOTIME/STEPT
 288             FRACT=I*STEPT-TOTIME
 289             I=(TOLD2+TOTIME)/STEPT
 290             J=TOTIME/STEPT
 291             NFRACT=I-J+ IONE
 292             IONE=0
 293             DO 110 I=1,NFRACT
 294   110       TSTEPS(I)=FRACT+I*STEPT
 295             TREF=TREF+NFRACT*STEPT
 296          ENDIF
 297       ELSEIF(STEPX.NE.0.D0) THEN
 298 C
 299 C   CRITERION FOR PRINTING RESULTS IS A CHANGE IN GEOMETRY.
 300 C
 301          IF(XOLD2+XOLD1-REFX.GT.STEPX)THEN
 302             NFRACT=(XOLD2+XOLD1-REFX)/STEPX
 303             CC=XOLD3(1)
 304             BB=XOLD3(2)
 305             AA=XOLD3(3)
 306             IF(ABS(BB/AA).GT.30)THEN
 307 C
 308 C   USE LINEAR INTERPOLATION
 309 C
 310                DO 120 I=1,NFRACT
 311   120          TSTEPS(I)=-(CC-(REFX+I*STEPX))/BB
 312             ELSE
 313 C
 314 C  USE QUADRATIC INTERPOLATION
 315 C
 316                DO 130 I=1,NFRACT
 317                   C1=CC-(REFX+I*STEPX)
 318   130          TSTEPS(I)=(-BB+SIGN(SQRT(BB*BB-4.D0*(AA*C1)),BB))/(2.D0*A
 319      1A)
 320             ENDIF
 321             REFX=REFX+NFRACT*STEPX
 322             FRACT=-.1
 323          ENDIF
 324       ELSE
 325 C
 326 C   PRINT EVERY POINT.
 327 C
 328          FRACT=0.0
 329       ENDIF
 330       IF(FRACT.LT.-9.D0)GOTO 170
 331       TURN=(TURN.OR.ABS(FRACT-1.D0).GT.1.D-6)
 332 C
 333 C  LOOP OVER ALL POINTS IN CURRENT DOMAIN
 334 C
 335       IF(FRACT.EQ.0.D0.AND.NFRACT.EQ.1)THEN
 336          TEXT1=' '
 337          TEXT2=' '
 338          II=0
 339          CALL DRCOUT(XYZ3,GEO3,VEL3,NVAR,TOTIME,ESCF3,EKIN3,
 340      1GTOT3,ETOT3,XTOT3,ILOOP,CHARGE,FRACT,TEXT1,TEXT2,II,JLOOP)
 341          N=0
 342          DO 140 I=1,NCOPRT
 343             K=MCOPRT(1,I)
 344             J=MCOPRT(2,I)
 345             L=K*3-3+J
 346             IF(ABS(GEO3(3,L)).GT.1.D-20)FRACT=-GEO3(2,L)/(GEO3(3,L)*2.D0
 347      1)
 348             IF(FRACT.GT.0.D0.AND.FRACT.LT.TOLD2) THEN
 349                IF(GEO3(3,L).GT.0.D0)TEXT1='MIN'
 350                IF(GEO3(3,L).LT.0.D0)TEXT1='MAX'
 351                TEXT2=COTYPE(J)
 352                IF(N.EQ.0)THEN
 353                   N=N+1
 354                   WRITE(6,'(/,20(''****''))')
 355                ENDIF
 356                TIME=TOTIME+FRACT
 357                CALL DRCOUT(XYZ3,GEO3,VEL3,NVAR,TIME,ESCF3,EKIN3,
 358      1GTOT3,ETOT3,XTOT3,ILOOP,CHARGE,FRACT,TEXT1,TEXT2,K,JLOOP)
 359             ENDIF
 360   140    CONTINUE
 361          IF(N.NE.0)WRITE(6,'(/,20(''****''))')
 362          IF(ABS(ESCF3(3)).GT.1.D-20)FRACT=-ESCF3(2)/(ESCF3(3)*2.D0)
 363          IF(.NOT.GOTURN.AND.FRACT.GT.0.D0.AND.FRACT.LT.TOLD2*1.04D0
 364      1.AND. PARMAX) THEN
 365             GOTURN=.TRUE.
 366             TIME=FRACT+TOTIME
 367             IF(ESCF3(3).GT.0.D0) THEN
 368                TEXT1='MIN'
 369                IF(LDRC) THEN
 370                   SUM=DOT(VELO0,VREF,NVAR)**2/(DOT(VELO0,VELO0,NVAR)*
 371      1DOT(VREF,VREF,NVAR)+1.D-10)
 372                   SUM1=DOT(VELO0,VREF0,NVAR)**2/(DOT(VELO0,VELO0,NVAR)*
 373      1DOT(VREF0,VREF0,NVAR)+1.D-10)
 374                   IF(SUM1.GT.0.1D0.AND.ABS(SUM1-1.D0).GT.1.D-6)
 375      1WRITE(6,'(/,A,F8.5,A,F8.5,A,G12.3,A)')' COEF. OF V(0)
 376      2=',SUM1,'   LAST V(0)',SUM,'   HALF-LIFE =',
 377      3-0.6931472D0*TIME/LOG(SUM1),' FEMTOSECS'
 378                ENDIF
 379                WRITE(6,'(//,A,F11.3,A)')' HALF-CYCLE TIME ='
 380      1,TIME-TLAST,' FEMTOSECONDS'
 381                TLAST=TIME
 382                DO 150 I=1,NVAR
 383   150          VREF(I)=VELO0(I)
 384             ENDIF
 385             IF(ESCF3(3).LT.0.D0)TEXT1='MAX'
 386             TEXT2=' '
 387             CALL DRCOUT(XYZ3,GEO3,VEL3,NVAR,TIME,ESCF3,EKIN3,
 388      1GTOT3,ETOT3,XTOT3,ILOOP,CHARGE,FRACT,TEXT1,TEXT2,0,JLOOP)
 389          ELSE
 390             GOTURN=.FALSE.
 391          ENDIF
 392       ELSE
 393          DO 160 I=1,NFRACT
 394             TIME=TOTIME+TSTEPS(I)
 395             TEXT1=' '
 396             TEXT2=' '
 397 C#           WRITE(6,'(A,4F12.4)')' KINETIC ENERGY, POINT',EKIN3,TSTEPS(
 398             CALL DRCOUT(XYZ3,GEO3,VEL3,NVAR,TIME,ESCF3,EKIN3,
 399      1GTOT3,ETOT3,XTOT3,ILOOP,CHARGE,TSTEPS(I),TEXT1,TEXT2,0,JLOOP)
 400   160    CONTINUE
 401       ENDIF
 402   170 CONTINUE
 403 C
 404 C BUFFER TOTAL TIME TO 3 POINTS BACK!
 405 C
 406       TOTIME=TOTIME+TOLD2
 407       TOLD2=TOLD1
 408       TOLD1=DELTAT
 409       ILOOP=ILOOP+1
 410       RETURN
 411       END