/* prtdrc.f -- translated by f2c (version 19950314). You must link the resulting object file with the libraries: -lF77 -lI77 -lm (in that order) */ #include "f2c.h" /* Common Block Declarations */ struct { char argz[512]; } argz_; #define argz_1 argz_ struct { char keywrd[80]; } keywrd_; #define keywrd_1 keywrd_ struct { integer numat, nat[86], nfirst[86], nmidle[86], nlast[86], norbs, nelecs, nalpha, nbeta, nclose, nopen, ndumy; doublereal xract; } molkst_; #define molkst_1 molkst_ struct { integer natoms, labels[86], na[86], nb[86], nc[86]; } geokst_; #define geokst_1 geokst_ struct { integer mcoprt[516] /* was [2][258] */, ncoprt; logical parmax; } drccom_; #define drccom_1 drccom_ struct { doublereal core[107]; } core_; #define core_1 core_ struct { doublereal atmass[86]; } atmass_; #define atmass_1 atmass_ struct { doublereal p[23220], pa[23220], pb[23220]; } densty_; #define densty_1 densty_ struct { doublereal allxyz[774] /* was [3][258] */, allvel[774] /* was [3][ 258] */, parref[258], xyz3[774] /* was [3][258] */, vel3[774] /* was [3][258] */, allgeo[774] /* was [3][258] */, geo3[774] /* was [3][258] */, dummy[28509]; } fmatrx_; #define fmatrx_1 fmatrx_ /* Table of constant values */ static integer c__1 = 1; static integer c__5 = 5; static integer c__3 = 3; static integer c__8 = 8; static doublereal c_b174 = 57.29577951; static integer c__0 = 0; /* Subroutine */ int prtdrc_(doublereal *escf, doublereal *deltt, doublereal * xparam, doublereal *ref, doublereal *ekin, doublereal *gtot, doublereal *etot, doublereal *velo0, integer *nvar) { /* Initialized data */ static logical first = TRUE_; static doublereal refscf = 0.; static char cotype[2*3] = "BL" "BA" "DI"; /* System generated locals */ integer i__1; doublereal d__1, d__2; /* Builtin functions */ double sqrt(doublereal); integer i_indx(char *, char *, ftnlen, ftnlen), s_wsfe(cilist *), do_fio( integer *, char *, ftnlen), e_wsfe(void), s_rsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), e_rsle(void), s_wsle(cilist *), e_wsle(void); double d_sign(doublereal *, doublereal *); /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen); double log(doublereal); /* Local variables */ static logical ldrc; static integer ione; static doublereal time, tref, vref[258], oldt, refx; static logical turn; static doublereal escf0, escf1, escf2, escf3[3], ekin0, ekin1, ekin2, ekin3[3], vref0[258], told1, told2, xold0, xold1, xold2, xold3[3], etot0, etot1, gtot0, etot3[3], gtot1, gtot3[3], etot2, gtot2; static char text1[3], text2[2]; static integer i, j; static doublereal sqrt2; extern doublereal reada_(char *, integer *, ftnlen); static integer l; static doublereal xtot0, xtot1, xtot3[3], xtot2; static integer n, k; extern /* Subroutine */ int chrge_(doublereal *, doublereal *); static doublereal fract, fincr; extern /* Subroutine */ int quadr_(doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *); static integer iloop; static doublereal steph, tlast; static integer jloop; static doublereal c1, stept, stepx, t1, t2, aa, bb, cc, dh; static integer ii; static doublereal charge[86], deltat; static integer nfract; static doublereal totime; extern /* Subroutine */ int drcout_(doublereal *, doublereal *, doublereal *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, integer *, doublereal * , doublereal *, char *, char *, integer *, integer *, ftnlen, ftnlen); static logical goturn; static doublereal tsteps[100]; extern /* Subroutine */ int xyzint_(doublereal *, integer *, integer *, integer *, integer *, doublereal *, doublereal *); static doublereal geo[258]; extern doublereal dot_(doublereal *, doublereal *, integer *); static doublereal sum, sum1; /* Fortran I/O blocks */ static cilist io___21 = { 0, 6, 0, "(/,' TIME PRIORITY, INTERVAL =',F4.1" ", ' FEMTOSECONDS',/)", 0 }; static cilist io___22 = { 0, 6, 0, "(/,' KINETIC ENERGY PRIORITY, STEP =" "',F5.2, ' KCAL/MOLE',/)", 0 }; static cilist io___23 = { 0, 6, 0, "(/,' GEOMETRY PRIORITY, STEP =',F5.2" ", ' ANGSTROMS',/)", 0 }; static cilist io___24 = { 0, 9, 0, 0, 0 }; static cilist io___25 = { 0, 9, 0, 0, 0 }; static cilist io___26 = { 0, 9, 0, 0, 0 }; static cilist io___27 = { 0, 9, 0, 0, 0 }; static cilist io___28 = { 0, 9, 0, 0, 0 }; static cilist io___29 = { 0, 9, 0, 0, 0 }; static cilist io___30 = { 0, 9, 0, 0, 0 }; static cilist io___31 = { 0, 9, 0, 0, 0 }; static cilist io___32 = { 0, 9, 0, 0, 0 }; static cilist io___33 = { 0, 9, 0, 0, 0 }; static cilist io___34 = { 0, 9, 0, 0, 0 }; static cilist io___35 = { 0, 9, 0, 0, 0 }; static cilist io___36 = { 0, 9, 0, 0, 0 }; static cilist io___52 = { 0, 9, 0, 0, 0 }; static cilist io___53 = { 0, 9, 0, 0, 0 }; static cilist io___54 = { 0, 9, 0, 0, 0 }; static cilist io___55 = { 0, 9, 0, 0, 0 }; static cilist io___56 = { 0, 9, 0, 0, 0 }; static cilist io___57 = { 0, 9, 0, 0, 0 }; static cilist io___58 = { 0, 9, 0, 0, 0 }; static cilist io___59 = { 0, 9, 0, 0, 0 }; static cilist io___60 = { 0, 9, 0, 0, 0 }; static cilist io___61 = { 0, 9, 0, 0, 0 }; static cilist io___62 = { 0, 9, 0, 0, 0 }; static cilist io___63 = { 0, 9, 0, 0, 0 }; static cilist io___64 = { 0, 9, 0, 0, 0 }; static cilist io___101 = { 0, 6, 0, "(/,20('****'))", 0 }; static cilist io___103 = { 0, 6, 0, "(/,20('****'))", 0 }; static cilist io___104 = { 0, 6, 0, "(/,A,F8.5,A,F8.5,A,F8.1,A)", 0 }; static cilist io___105 = { 0, 6, 0, "(//,A,F11.3,A)", 0 }; /* ******************************************************************** */ /* PRTDRC PREPARES TO PRINT THE GEOMETRY ETC. FOR POINTS IN A DRC */ /* OR IRC */ /* CALCULATION. */ /* ON INPUT ESCF = HEAT OF FORMATION FOR THE CURRENT POINT */ /* DELTT = CHANGE IN TIME, PREVIOUS TO CURRENT POINT */ /* XPARAM = CURRENT CARTESIAN GEOMETRY */ /* EKIN = CURRENT KINETIC ENERGY */ /* GTOT = TOTAL GRADIENT NORM IN IRC CALC'N. */ /* VELO0 = CURRENT VELOCITY */ /* NVAR = NUMBER OF VARIABLES = 3 * NUMBER OF ATOMS. */ /* ******************************************************************* */ /* COMDECK SIZES */ /************************************************************************ ****/ /* THIS FILE CONTAINS ALL THE ARRAY SIZES FOR USE IN MOPAC. **/ /* **/ /* THERE ARE ONLY PARAMETERS THAT THE PROGRAMMER NEED SET: **/ /* MAXHEV = MAXIMUM NUMBER OF HEAVY ATOMS (HEAVY: NON-HYDROGEN ATOMS) **/ /* MAXLIT = MAXIMUM NUMBER OF HYDROGEN ATOMS. **/ /* MAXTIM = DEFAULT TIME FOR A JOB. (SECONDS) **/ /* MAXDMP = DEFAULT TIME FOR AUTOMATIC RESTART FILE GENERATION (SECS) **/ /* **/ /* **/ /************************************************************************ ****/ /* **/ /* THE FOLLOWING CODE DOES NOT NEED TO BE ALTERED BY THE PROGRAMMER **/ /* **/ /************************************************************************ ****/ /* **/ /* ALL OTHER PARAMETERS ARE DERIVED FUNCTIONS OF THESE TWO PARAMETERS **/ /* **/ /* NAME DEFINITION **/ /* NUMATM MAXIMUM NUMBER OF ATOMS ALLOWED. **/ /* MAXORB MAXIMUM NUMBER OF ORBITALS ALLOWED. **/ /* MAXPAR MAXIMUM NUMBER OF PARAMETERS FOR OPTIMISATION. **/ /* N2ELEC MAXIMUM NUMBER OF TWO ELECTRON INTEGRALS ALLOWED. **/ /* MPACK AREA OF LOWER HALF TRIANGLE OF DENSITY MATRIX. **/ /* MORB2 SQUARE OF THE MAXIMUM NUMBER OF ORBITALS ALLOWED. **/ /* MAXHES AREA OF HESSIAN MATRIX **/ /************************************************************************ ****/ /************************************************************************ ****/ /* FOR SHORT VERSION USE LINE WITH NMECI=1, FOR LONG VERSION USE LINE **/ /* WITH NMECI=10 **/ /************************************************************************ ****/ /* PARAMETER (NMECI=1, NPULAY=1) */ /************************************************************************ ****/ /* DECK MOPAC */ /* next line added for Unix implementation for command line arguments */ /* Parameter adjustments */ --velo0; --ref; --xparam; /* Function Body */ if (first) { sqrt2 = sqrt(2.); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { /* L10: */ fmatrx_1.parref[i - 1] = xparam[i]; } *etot = *escf + *ekin; tlast = 0.; goturn = FALSE_; sum = 0.; i__1 = *nvar; for (i = 1; i <= i__1; ++i) { /* Computing 2nd power */ d__1 = velo0[i]; sum += d__1 * d__1; vref0[i - 1] = velo0[i]; /* L20: */ vref[i - 1] = velo0[i]; } ione = 1; ldrc = sum > 1.; first = FALSE_; iloop = 1; oldt = -100.; told1 = 0.; /* DETERMINE TYPE OF PRINT: TIME, ENERGY OR GEOMETRY PRIORITY */ /* OR PRINT ALL POINTS */ stept = 0.; steph = 0.; stepx = 0.; if (i_indx(keywrd_1.keywrd, " T-PRIO", 80L, 7L) != 0) { if (i_indx(keywrd_1.keywrd, " T-PRIORITY=", 80L, 12L) != 0) { i__1 = i_indx(keywrd_1.keywrd, "T-PRIO", 80L, 6L) + 5; stept = reada_(keywrd_1.keywrd, &i__1, 80L); } else { stept = .1; } tref = -1e-6; s_wsfe(&io___21); do_fio(&c__1, (char *)&stept, (ftnlen)sizeof(doublereal)); e_wsfe(); } else if (i_indx(keywrd_1.keywrd, " H-PRIO", 80L, 7L) != 0) { if (i_indx(keywrd_1.keywrd, " H-PRIORITY=", 80L, 12L) != 0) { i__1 = i_indx(keywrd_1.keywrd, "H-PRIO", 80L, 6L) + 5; steph = reada_(keywrd_1.keywrd, &i__1, 80L); } else { steph = .1; } s_wsfe(&io___22); do_fio(&c__1, (char *)&steph, (ftnlen)sizeof(doublereal)); e_wsfe(); } else if (i_indx(keywrd_1.keywrd, " X-PRIO", 80L, 7L) != 0) { if (i_indx(keywrd_1.keywrd, " X-PRIORITY=", 80L, 12L) != 0) { i__1 = i_indx(keywrd_1.keywrd, "X-PRIO", 80L, 6L) + 5; stepx = reada_(keywrd_1.keywrd, &i__1, 80L); } else { stepx = .05; } s_wsfe(&io___23); do_fio(&c__1, (char *)&stepx, (ftnlen)sizeof(doublereal)); e_wsfe(); } if (i_indx(keywrd_1.keywrd, " REST", 80L, 5L) != 0 && i_indx( keywrd_1.keywrd, "IRC=", 80L, 4L) == 0) { s_rsle(&io___24); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.parref[i - 1], (ftnlen) sizeof(doublereal)); } e_rsle(); s_rsle(&io___25); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&vref0[i - 1], (ftnlen)sizeof( doublereal)); } e_rsle(); s_rsle(&io___26); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&vref[i - 1], (ftnlen)sizeof( doublereal)); } e_rsle(); s_rsle(&io___27); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allgeo[i * 3 - 1], ( ftnlen)sizeof(doublereal)); } e_rsle(); s_rsle(&io___28); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allgeo[i * 3 - 2], ( ftnlen)sizeof(doublereal)); } e_rsle(); s_rsle(&io___29); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allgeo[i * 3 - 3], ( ftnlen)sizeof(doublereal)); } e_rsle(); s_rsle(&io___30); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allvel[i * 3 - 1], ( ftnlen)sizeof(doublereal)); } e_rsle(); s_rsle(&io___31); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allvel[i * 3 - 2], ( ftnlen)sizeof(doublereal)); } e_rsle(); s_rsle(&io___32); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allvel[i * 3 - 3], ( ftnlen)sizeof(doublereal)); } e_rsle(); s_rsle(&io___33); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allxyz[i * 3 - 1], ( ftnlen)sizeof(doublereal)); } e_rsle(); s_rsle(&io___34); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allxyz[i * 3 - 2], ( ftnlen)sizeof(doublereal)); } e_rsle(); s_rsle(&io___35); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allxyz[i * 3 - 3], ( ftnlen)sizeof(doublereal)); } e_rsle(); s_rsle(&io___36); do_lio(&c__3, &c__1, (char *)&iloop, (ftnlen)sizeof(integer)); do_lio(&c__8, &c__1, (char *)&ldrc, (ftnlen)sizeof(logical)); do_lio(&c__3, &c__1, (char *)&ione, (ftnlen)sizeof(integer)); do_lio(&c__5, &c__1, (char *)&etot1, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&etot0, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&escf1, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&escf0, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&ekin1, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&ekin0, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&told2, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&told1, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)>ot1, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)>ot0, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&xold2, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&xold1, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&xold0, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&totime, (ftnlen)sizeof(doublereal)); do_lio(&c__3, &c__1, (char *)&jloop, (ftnlen)sizeof(integer)); do_lio(&c__5, &c__1, (char *)&(*etot), (ftnlen)sizeof(doublereal)) ; do_lio(&c__5, &c__1, (char *)&refx, (ftnlen)sizeof(doublereal)); e_rsle(); } } if (*escf < -1e8) { s_wsle(&io___52); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.parref[i - 1], (ftnlen) sizeof(doublereal)); } e_wsle(); s_wsle(&io___53); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&vref0[i - 1], (ftnlen)sizeof( doublereal)); } e_wsle(); s_wsle(&io___54); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&vref[i - 1], (ftnlen)sizeof( doublereal)); } e_wsle(); s_wsle(&io___55); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allgeo[i * 3 - 1], (ftnlen) sizeof(doublereal)); } e_wsle(); s_wsle(&io___56); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allgeo[i * 3 - 2], (ftnlen) sizeof(doublereal)); } e_wsle(); s_wsle(&io___57); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allgeo[i * 3 - 3], (ftnlen) sizeof(doublereal)); } e_wsle(); s_wsle(&io___58); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allvel[i * 3 - 1], (ftnlen) sizeof(doublereal)); } e_wsle(); s_wsle(&io___59); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allvel[i * 3 - 2], (ftnlen) sizeof(doublereal)); } e_wsle(); s_wsle(&io___60); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allvel[i * 3 - 3], (ftnlen) sizeof(doublereal)); } e_wsle(); s_wsle(&io___61); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allxyz[i * 3 - 1], (ftnlen) sizeof(doublereal)); } e_wsle(); s_wsle(&io___62); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allxyz[i * 3 - 2], (ftnlen) sizeof(doublereal)); } e_wsle(); s_wsle(&io___63); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_lio(&c__5, &c__1, (char *)&fmatrx_1.allxyz[i * 3 - 3], (ftnlen) sizeof(doublereal)); } e_wsle(); s_wsle(&io___64); do_lio(&c__3, &c__1, (char *)&iloop, (ftnlen)sizeof(integer)); do_lio(&c__8, &c__1, (char *)&ldrc, (ftnlen)sizeof(logical)); do_lio(&c__3, &c__1, (char *)&ione, (ftnlen)sizeof(integer)); do_lio(&c__5, &c__1, (char *)&etot1, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&etot0, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&escf1, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&escf0, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&ekin1, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&ekin0, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&told2, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&told1, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)>ot1, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)>ot0, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&xold2, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&xold1, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&xold0, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&totime, (ftnlen)sizeof(doublereal)); do_lio(&c__3, &c__1, (char *)&jloop, (ftnlen)sizeof(integer)); do_lio(&c__5, &c__1, (char *)&(*etot), (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&refx, (ftnlen)sizeof(doublereal)); e_wsle(); return 0; } chrge_(densty_1.p, charge); i__1 = molkst_1.numat; for (i = 1; i <= i__1; ++i) { l = molkst_1.nat[i - 1]; /* L30: */ charge[i - 1] = core_1.core[l - 1] - charge[i - 1]; } deltat = *deltt * 1e15; geokst_1.na[1] = -1; xyzint_(&xparam[1], &molkst_1.numat, geokst_1.na, geokst_1.nb, geokst_1.nc, &c_b174, geo); geokst_1.na[0] = 99; if (iloop == 1) { etot1 = etot0; etot0 = *etot; escf1 = *escf; escf0 = *escf; ekin1 = *ekin; ekin0 = *ekin; for (j = 1; j <= 3; ++j) { i__1 = *nvar; for (i = 1; i <= i__1; ++i) { fmatrx_1.allgeo[j + i * 3 - 4] = geo[i - 1]; fmatrx_1.allxyz[j + i * 3 - 4] = xparam[i]; /* L40: */ fmatrx_1.allvel[j + i * 3 - 4] = velo0[i]; } } } else { i__1 = *nvar; for (i = 1; i <= i__1; ++i) { fmatrx_1.allgeo[i * 3 - 1] = fmatrx_1.allgeo[i * 3 - 2]; fmatrx_1.allgeo[i * 3 - 2] = fmatrx_1.allgeo[i * 3 - 3]; fmatrx_1.allgeo[i * 3 - 3] = geo[i - 1]; fmatrx_1.allxyz[i * 3 - 1] = fmatrx_1.allxyz[i * 3 - 2]; fmatrx_1.allxyz[i * 3 - 2] = fmatrx_1.allxyz[i * 3 - 3]; fmatrx_1.allxyz[i * 3 - 3] = xparam[i]; fmatrx_1.allvel[i * 3 - 1] = fmatrx_1.allvel[i * 3 - 2]; fmatrx_1.allvel[i * 3 - 2] = fmatrx_1.allvel[i * 3 - 3]; /* L50: */ fmatrx_1.allvel[i * 3 - 3] = velo0[i]; } } /* FORM QUADRATIC EXPRESSION FOR POSITION AND VELOCITY W.R.T. TIME. */ t1 = max(told2,.02); t2 = max(told1,.02) + t1; i__1 = *nvar; for (i = 1; i <= i__1; ++i) { quadr_(&fmatrx_1.allgeo[i * 3 - 1], &fmatrx_1.allgeo[i * 3 - 2], & fmatrx_1.allgeo[i * 3 - 3], &t1, &t2, &fmatrx_1.geo3[i * 3 - 3], &fmatrx_1.geo3[i * 3 - 2], &fmatrx_1.geo3[i * 3 - 1]); /* *************************************************** */ /* * */ /* QUADR CALCULATES THE A, B AND C IN THE EQUNS. * */ /* * */ /* A = F0 * */ /* A + B.X0 + C.X0**2 = F1 * */ /* A + B.X2 + C.X2**2 = F2 * */ /* GIVEN THE ARGUMENT LIST (F0,F1,F2, X1,X2, A,B,C) * */ /* * */ /* *************************************************** */ quadr_(&fmatrx_1.allxyz[i * 3 - 1], &fmatrx_1.allxyz[i * 3 - 2], & fmatrx_1.allxyz[i * 3 - 3], &t1, &t2, &fmatrx_1.xyz3[i * 3 - 3], &fmatrx_1.xyz3[i * 3 - 2], &fmatrx_1.xyz3[i * 3 - 1]); quadr_(&fmatrx_1.allvel[i * 3 - 1], &fmatrx_1.allvel[i * 3 - 2], & fmatrx_1.allvel[i * 3 - 3], &t1, &t2, &fmatrx_1.vel3[i * 3 - 3], &fmatrx_1.vel3[i * 3 - 2], &fmatrx_1.vel3[i * 3 - 1]); /* L60: */ } etot2 = etot1; etot1 = etot0; etot0 = *etot; quadr_(&etot2, &etot1, &etot0, &t1, &t2, etot3, &etot3[1], &etot3[2]); ekin2 = ekin1; ekin1 = ekin0; ekin0 = *ekin; quadr_(&ekin2, &ekin1, &ekin0, &t1, &t2, ekin3, &ekin3[1], &ekin3[2]); escf2 = escf1; escf1 = escf0; escf0 = *escf; quadr_(&escf2, &escf1, &escf0, &t1, &t2, escf3, &escf3[1], &escf3[2]); gtot2 = gtot1; gtot1 = gtot0; gtot0 = *gtot; quadr_(>ot2, >ot1, >ot0, &t1, &t2, gtot3, >ot3[1], >ot3[2]); xtot2 = xtot1; xtot1 = xtot0; xold2 += xold1; xold1 = xold0; /* CALCULATE CHANGE IN GEOMETRY */ xold0 = 0.; l = 0; xtot0 = 0.; sum1 = 0.; i__1 = molkst_1.numat; for (i = 1; i <= i__1; ++i) { sum = 0.; for (j = 1; j <= 3; ++j) { ++l; /* Computing 2nd power */ d__1 = fmatrx_1.allxyz[l * 3 - 3] - ref[l]; sum1 += d__1 * d__1 * atmass_1.atmass[i - 1]; /* L70: */ /* Computing 2nd power */ d__1 = fmatrx_1.allxyz[l * 3 - 2] - fmatrx_1.allxyz[l * 3 - 3]; sum += d__1 * d__1; } /* L80: */ xold0 += sqrt(sum); } xtot0 += sqrt(sum1) / sqrt2; quadr_(&xtot2, &xtot1, &xtot0, &t1, &t2, xtot3, &xtot3[1], &xtot3[2]); d__1 = xold2 + xold1; d__2 = xold2 + xold1 + xold0; quadr_(&xold2, &d__1, &d__2, &t1, &t2, xold3, &xold3[1], &xold3[2]); /* ********************************************************************** */ /* GO THROUGH THE CRITERIA FOR DECIDING WHETHER OR NOT TO PRINT THIS * */ /* POINT. IF YES, THEN ALSO CALCULATE THE EXACT POINT AS A FRACTION * */ /* BETWEEN THE LAST POINT AND THE CURRENT POINT * */ /* ********************************************************************** */ /* NFRACT IS THE NUMBER OF POINTS TO BE PRINTED IN THE CURRENT DOMAIN */ /* ********************************************************************** */ if (iloop < 3) { goto L170; } fract = -10.; nfract = 1; if (steph != 0.) { /* CRITERION FOR PRINTING RESULTS IS A CHANGE IN HEAT OF FORMATION = */ /* -CHANGE IN KINETIC ENERGY */ if (refscf == 0.) { i = (integer) (escf2 / steph); refscf = i * steph; } dh = (d__1 = escf1 - refscf, abs(d__1)); if (dh > steph) { d__1 = escf1 - refscf; steph = d_sign(&steph, &d__1); nfract = (d__1 = dh / steph, (integer) abs(d__1)); cc = escf3[0]; bb = escf3[1]; aa = escf3[2]; /* *********************************************** */ /* PROGRAMMERS! - BE VERY CAREFUL IF YOU CHANGE * */ /* THIS FOLLOWING SECTION. THERE IS NUMERICAL * */ /* INSTABILITY IF ABS(BB/AA) IS VERY LARGE. NEAR* */ /* INFLECTION POINTS AA CHANGES SIGN. JJPS* */ /* *********************************************** */ if ((d__1 = bb / aa, abs(d__1)) > 30.) { /* USE LINEAR INTERPOLATION */ i__1 = nfract; for (i = 1; i <= i__1; ++i) { /* L90: */ tsteps[i - 1] = -(cc - (refscf + i * steph)) / bb; } } else { /* USE QUADRATIC INTERPOLATION */ i__1 = nfract; for (i = 1; i <= i__1; ++i) { c1 = cc - (refscf + i * steph); /* L100: */ d__1 = sqrt(bb * bb - aa * c1 * 4.); tsteps[i - 1] = (-bb + d_sign(&d__1, &bb)) / (aa * 2.); } } fract = -.1f; refscf += nfract * steph; } } else if (stept != 0.) { /* CRITERION FOR PRINTING RESULTS IS A CHANGE IN TIME. */ if ((d__1 = totime + told2 - tref, abs(d__1)) > stept) { fincr = stept; i = (integer) (totime / stept); fract = i * stept - totime; i = (integer) ((told2 + totime) / stept); j = (integer) (totime / stept); nfract = i - j + ione; ione = 0; i__1 = nfract; for (i = 1; i <= i__1; ++i) { /* L110: */ tsteps[i - 1] = fract + i * stept; } tref += nfract * stept; } } else if (stepx != 0.) { /* CRITERION FOR PRINTING RESULTS IS A CHANGE IN GEOMETRY. */ if (xold2 + xold1 - refx > stepx) { nfract = (integer) ((xold2 + xold1 - refx) / stepx); cc = xold3[0]; bb = xold3[1]; aa = xold3[2]; if ((d__1 = bb / aa, abs(d__1)) > 30.) { /* USE LINEAR INTERPOLATION */ i__1 = nfract; for (i = 1; i <= i__1; ++i) { /* L120: */ tsteps[i - 1] = -(cc - (refx + i * stepx)) / bb; } } else { /* USE QUADRATIC INTERPOLATION */ i__1 = nfract; for (i = 1; i <= i__1; ++i) { c1 = cc - (refx + i * stepx); /* L130: */ d__1 = sqrt(bb * bb - aa * c1 * 4.); tsteps[i - 1] = (-bb + d_sign(&d__1, &bb)) / (aa * 2.); } } refx += nfract * stepx; fract = -.1f; } } else { /* PRINT EVERY POINT. */ fract = 0.f; } if (fract < -9.) { goto L170; } turn = turn || (d__1 = fract - 1., abs(d__1)) > 1e-6; /* LOOP OVER ALL POINTS IN CURRENT DOMAIN */ if (fract == 0. && nfract == 1) { s_copy(text1, " ", 3L, 1L); s_copy(text2, " ", 2L, 1L); ii = 0; drcout_(fmatrx_1.xyz3, fmatrx_1.geo3, fmatrx_1.vel3, nvar, &totime, escf3, ekin3, gtot3, etot3, xtot3, &iloop, charge, &fract, text1, text2, &ii, &jloop, 3L, 2L); n = 0; i__1 = drccom_1.ncoprt; for (i = 1; i <= i__1; ++i) { k = drccom_1.mcoprt[(i << 1) - 2]; j = drccom_1.mcoprt[(i << 1) - 1]; l = k * 3 - 3 + j; if ((d__1 = fmatrx_1.geo3[l * 3 - 1], abs(d__1)) > 1e-20) { fract = -fmatrx_1.geo3[l * 3 - 2] / (fmatrx_1.geo3[l * 3 - 1] * 2.); } if (fract > 0. && fract < told2) { if (fmatrx_1.geo3[l * 3 - 1] > 0.) { s_copy(text1, "MIN", 3L, 3L); } if (fmatrx_1.geo3[l * 3 - 1] < 0.) { s_copy(text1, "MAX", 3L, 3L); } s_copy(text2, cotype + (j - 1 << 1), 2L, 2L); if (n == 0) { ++n; s_wsfe(&io___101); e_wsfe(); } time = totime + fract; drcout_(fmatrx_1.xyz3, fmatrx_1.geo3, fmatrx_1.vel3, nvar, & time, escf3, ekin3, gtot3, etot3, xtot3, &iloop, charge, &fract, text1, text2, &k, &jloop, 3L, 2L); } /* L140: */ } if (n != 0) { s_wsfe(&io___103); e_wsfe(); } if (abs(escf3[2]) > 1e-20) { fract = -escf3[1] / (escf3[2] * 2.); } if (! goturn && fract > 0. && fract < told2 * 1.04 && drccom_1.parmax) { goturn = TRUE_; time = fract + totime; if (escf3[2] > 0.) { s_copy(text1, "MIN", 3L, 3L); if (ldrc) { /* Computing 2nd power */ d__1 = dot_(&velo0[1], vref, nvar); sum = d__1 * d__1 / (dot_(&velo0[1], &velo0[1], nvar) * dot_(vref, vref, nvar) + 1e-10); /* Computing 2nd power */ d__1 = dot_(&velo0[1], vref0, nvar); sum1 = d__1 * d__1 / (dot_(&velo0[1], &velo0[1], nvar) * dot_(vref0, vref0, nvar) + 1e-10); if (sum1 > .1) { s_wsfe(&io___104); do_fio(&c__1, " COEF. OF V(0) =", 29L); do_fio(&c__1, (char *)&sum1, (ftnlen)sizeof( doublereal)); do_fio(&c__1, " LAST V(0)", 12L); do_fio(&c__1, (char *)&sum, (ftnlen)sizeof(doublereal) ); do_fio(&c__1, " HALF-LIFE =", 14L); d__1 = time * -.6931472 / log(sum1); do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof( doublereal)); do_fio(&c__1, " FEMTOSECS", 10L); e_wsfe(); } } s_wsfe(&io___105); do_fio(&c__1, " HALF-CYCLE TIME =", 18L); d__1 = time - tlast; do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof(doublereal)); do_fio(&c__1, " FEMTOSECONDS", 13L); e_wsfe(); tlast = time; i__1 = *nvar; for (i = 1; i <= i__1; ++i) { /* L150: */ vref[i - 1] = velo0[i]; } } if (escf3[2] < 0.) { s_copy(text1, "MAX", 3L, 3L); } s_copy(text2, " ", 2L, 1L); drcout_(fmatrx_1.xyz3, fmatrx_1.geo3, fmatrx_1.vel3, nvar, &time, escf3, ekin3, gtot3, etot3, xtot3, &iloop, charge, &fract, text1, text2, &c__0, &jloop, 3L, 2L); } else { goturn = FALSE_; } } else { i__1 = nfract; for (i = 1; i <= i__1; ++i) { time = totime + tsteps[i - 1]; s_copy(text1, " ", 3L, 1L); s_copy(text2, " ", 2L, 1L); /* # WRITE(6,'(A,4F12.4)')' KINETIC ENERGY, POINT',EKIN3 ,TSTEPS( */ drcout_(fmatrx_1.xyz3, fmatrx_1.geo3, fmatrx_1.vel3, nvar, &time, escf3, ekin3, gtot3, etot3, xtot3, &iloop, charge, & tsteps[i - 1], text1, text2, &c__0, &jloop, 3L, 2L); /* L160: */ } } L170: /* BUFFER TOTAL TIME TO 3 POINTS BACK! */ totime += told2; told2 = told1; told1 = deltat; ++iloop; return 0; } /* prtdrc_ */