/* force.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 { integer nvar, loc[516] /* was [2][258] */, idumy; doublereal dumy[258]; } geovar_; #define geovar_1 geovar_ struct { integer ndep, locpar[258], idepfn[258], locdep[258]; } geosym_; #define geosym_1 geosym_ struct { integer natoms, labels[86], na[86], nb[86], nc[86]; } geokst_; #define geokst_1 geokst_ struct { doublereal fmatrx[33411]; } fmatrx_; #define fmatrx_1 fmatrx_ struct { char keywrd[80]; } keywrd_; #define keywrd_1 keywrd_ struct { doublereal grad[258], gnorm; } gradnt_; #define gradnt_1 gradnt_ struct { doublereal cnorml[66564], freq[258], dummy[26058]; } vector_; #define vector_1 vector_ struct { char elemnt[214]; } elemts_; #define elemts_1 elemts_ struct { integer last; } last_; #define last_1 last_ struct { doublereal geo[258] /* was [3][86] */; } geom_; #define geom_1 geom_ struct { doublereal coord[258] /* was [3][86] */; } coord_; #define coord_1 coord_ struct { doublereal tvec[9] /* was [3][3] */; integer id; } euler_; #define euler_1 euler_ struct { doublereal store[66564]; } scrach_; #define scrach_1 scrach_ /* Table of constant values */ static logical c_true = TRUE_; static logical c_false = FALSE_; static integer c__1 = 1; static doublereal c_b42 = 1.; static integer c__2 = 2; static integer c__0 = 0; /* Subroutine */ int force_(void) { /* System generated locals */ integer i__1, i__2, i__3, i__4; doublereal d__1, d__2, d__3, d__4, d__5, d__6; /* Builtin functions */ integer i_indx(char *, char *, ftnlen, ftnlen), s_wsfe(cilist *), do_fio( integer *, char *, ftnlen), e_wsfe(void); double sqrt(doublereal); /* Subroutine */ int s_stop(char *, ftnlen); /* Local variables */ static doublereal escf; extern /* Subroutine */ int fmat_(doublereal *, integer *, doublereal *, doublereal *, doublereal *, doublereal *); static integer ilim; static doublereal tder, tscf, dipt[258]; static integer nvib; extern /* Subroutine */ int axis_(doublereal *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, integer *, doublereal *) ; #define keys ((char *)&keywrd_1) static doublereal summ; static logical prnt; static doublereal escf1; static integer iinc1, iinc2; static doublereal time1, time2, time3; static integer nrem6; static doublereal a, b, c; static integer i, j, k, l; extern doublereal reada_(char *, integer *, ftnlen); static logical debug, large; extern /* Subroutine */ int frame_(doublereal *, integer *, integer *, doublereal *), flepo_(doublereal *, integer *, doublereal *); static doublereal shift[6], trdip[774] /* was [3][258] */; static integer numat; static doublereal const__; extern /* Subroutine */ int nllsq_(doublereal *, integer *), write_( doublereal *, doublereal *); static doublereal wtmol; static integer ii, ij, jj, il, jl; #define gr ((doublereal *)&gradnt_1) static integer iu, ju; extern /* Subroutine */ int anavib_(doublereal *, doublereal *, doublereal *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *); static doublereal deldip[774] /* was [3][258] */; static logical bartel, linear; static doublereal redmas[258]; static integer locold[516] /* was [2][258] */; extern doublereal second_(void); extern /* Subroutine */ int compfg_(doublereal *, logical *, doublereal *, logical *, doublereal *, logical *); static integer nvaold; static doublereal xparam[258], travel[258]; extern /* Subroutine */ int freqcy_(doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, logical *), vecprt_( doublereal *, integer *), thermo_(doublereal *, doublereal *, doublereal *, logical *, doublereal *, doublereal *, doublereal *, integer *, doublereal *), gmetry_(doublereal *, doublereal *); static integer im1; extern /* Subroutine */ int matout_(doublereal *, doublereal *, integer *, integer *, integer *); static logical restrt; extern /* Subroutine */ int xyzint_(doublereal *, integer *, integer *, integer *, integer *, doublereal *, doublereal *), drc_( doublereal *, doublereal *); extern doublereal dot_(doublereal *, doublereal *, integer *); static doublereal rot[9] /* was [3][3] */, sum; extern /* Subroutine */ int rsp_(doublereal *, integer *, integer *, doublereal *, doublereal *); static doublereal sym, sum1; static integer nto6; /* Fortran I/O blocks */ static cilist io___20 = { 0, 6, 0, "(//10X,'HEAT OF FORMATION =',F12.6, " " ' KCALS/MOLE')", 0 }; static cilist io___23 = { 0, 6, 0, "(//10X,'INTERNAL COORDINATE DERIVATI" "VES',//3X, 'ATOM AT. NO.',2X,'BOND',9X,'ANGLE',8X,'DIHEDRAL',/)", 0 }; static cilist io___27 = { 0, 6, 0, "(I6,4X,A2,F13.6,2F10.6)", 0 }; static cilist io___28 = { 0, 6, 0, "(//10X,'GRADIENT NORM =',F10.5)", 0 }; static cilist io___29 = { 0, 6, 0, "(///1X,'** GRADIENT IS VERY LARGE, B" "UT SINCE \"LET\"',' IS USED, CALCULATION WILL CONTINUE')", 0 }; static cilist io___30 = { 0, 6, 0, "(///1X,'** GRADIENT IS TOO LARGE TO " "ALLOW ', 'FORCE MATRIX TO BE CALCULATED, (LIMIT=10" ") **',//)", 0 }; static cilist io___31 = { 0, 6, 0, "(//10X,' GEOMETRY WILL BE OPTIMIZED " "FIRST')", 0 }; static cilist io___32 = { 0, 6, 0, "(15X,'USING NLLSQ')", 0 }; static cilist io___33 = { 0, 6, 0, "(15X,'USING FLEPO')", 0 }; static cilist io___34 = { 0, 6, 0, "(//10X,'GRADIENT NORM =',F10.7)", 0 }; static cilist io___35 = { 0, 6, 0, "(//30X,'**** WARNING ****',// " " 10X,' GRADIENT IS VERY LARGE FOR A THERMO CALCUL" "ATION',/ 10X,' RESULTS ARE LIKELY TO BE INACCURATE IF THE" "RE ARE')", 0 }; static cilist io___36 = { 0, 6, 0, "(10X,' ANY LOW-LYING VIBRATIONS (LES" "S THAN ABOUT ' ,'400CM-1)')", 0 }; static cilist io___37 = { 0, 6, 0, "(10X,' GRADIENT NORM SHOULD BE LESS " "THAN ABOUT ', '0.2 FOR THERMO',/10X,' TO GIVE ACCURATE RESULTS" "')", 0 }; static cilist io___38 = { 0, 6, 0, "(//10X,'TIME FOR SCF CALCULATION ='," "F8.2)", 0 }; static cilist io___39 = { 0, 6, 0, "(//10X,'TIME FOR DERIVATIVES ='," "F8.2)", 0 }; static cilist io___40 = { 0, 6, 0, "(//10X,'SYMMETRY WAS SPECIFIED, BUT " "', 'CANNOT BE USED HERE')", 0 }; static cilist io___47 = { 0, 6, 0, "(/9X,'ORIENTATION OF MOLECULE IN FOR" "CE CALCULATION')", 0 }; static cilist io___48 = { 0, 6, 0, "(/,4X,'NO.',7X,'ATOM',9X,'X', " " 9X,'Y',9X,'Z',/)", 0 }; static cilist io___49 = { 0, 6, 0, "(I6,7X,I3,4X,3F10.4)", 0 }; static cilist io___58 = { 0, 6, 0, "(//10X,' FULL FORCE MATRIX, INVOKED " "BY \"DFORCE\"')", 0 }; static cilist io___59 = { 0, 6, 0, "(//10X,' FORCE MATRIX IN MILLIDYNES/" "ANGSTROM')", 0 }; static cilist io___60 = { 0, 6, 0, "(//10X,'HEAT OF FORMATION =',F12.6, " " ' KCALS/MOLE')", 0 }; static cilist io___62 = { 0, 6, 0, "(//10X,'TRIVIAL VIBRATIONS, SHOULD B" "E ZERO')", 0 }; static cilist io___63 = { 0, 6, 0, "(/, F9.4,'=TX',F9.4,'=TY',F9.4,'=TZ'" ", F9.4,'=RX',F9.4,'=RY',F9.4,'=RZ')", 0 }; static cilist io___64 = { 0, 6, 0, "(//10X,'FORCE CONSTANTS IN MILLIDYNE" "S/ANGSTROM' ,' (= 10**5 DYNES/CM)',/)", 0 }; static cilist io___65 = { 0, 6, 0, "(8F10.5)", 0 }; static cilist io___66 = { 0, 6, 0, "(//10X,' ASSOCIATED EIGENVECTORS')", 0 }; static cilist io___70 = { 0, 6, 0, "(//10X,' ZERO POINT ENERGY' " " , F12.3,' KILOCALORIES PER MOLE')", 0 }; static cilist io___76 = { 0, 6, 0, "(//3X,' THE LAST',I2,' VIBRATIONS AR" "E THE', ' TRANSLATION AND ROTATION MODES')", 0 }; static cilist io___77 = { 0, 6, 0, "(3X,' THE FIRST THREE OF THESE BEING" " TRANSLATIONS', ' IN X, Y, AND Z, RESPECTIVELY')", 0 }; static cilist io___78 = { 0, 6, 0, "(//10X,' FREQUENCIES, REDUCED MASSES" " AND ', 'VIBRATIONAL DIPOLES'/)", 0 }; static cilist io___82 = { 0, 6, 0, "(/)", 0 }; static cilist io___84 = { 0, 6, 0, "(3X,'I',10I10)", 0 }; static cilist io___85 = { 0, 6, 0, "(' FREQ(I)',6F10.4,/)", 0 }; static cilist io___86 = { 0, 6, 0, "(' MASS(I)',6F10.5,/)", 0 }; static cilist io___87 = { 0, 6, 0, "(' DIPX(I)',6F10.5)", 0 }; static cilist io___88 = { 0, 6, 0, "(' DIPY(I)',6F10.5)", 0 }; static cilist io___89 = { 0, 6, 0, "(' DIPZ(I)',6F10.5,/)", 0 }; static cilist io___90 = { 0, 6, 0, "(' DIPT(I)',6F10.5)", 0 }; static cilist io___91 = { 0, 6, 0, "(/)", 0 }; static cilist io___92 = { 0, 6, 0, "(3X,'I',10I10)", 0 }; static cilist io___93 = { 0, 6, 0, "(' FREQ(I)',6F10.4)", 0 }; static cilist io___94 = { 0, 6, 0, "(/,' MASS(I)',6F10.5)", 0 }; static cilist io___95 = { 0, 6, 0, "(/,' DIPX(I)',6F10.5)", 0 }; static cilist io___96 = { 0, 6, 0, "(' DIPY(I)',6F10.5)", 0 }; static cilist io___97 = { 0, 6, 0, "(' DIPZ(I)',6F10.5)", 0 }; static cilist io___98 = { 0, 6, 0, "(/,' DIPT(I)',6F10.5)", 0 }; static cilist io___99 = { 0, 6, 0, "(//10X,' NORMAL COORDINATE ANALYSIS')" , 0 }; static cilist io___100 = { 0, 6, 0, "(//10X,' MASS-WEIGHTED COORDINATE A" "NALYSIS')", 0 }; static cilist io___103 = { 0, 6, 0, "(//1X,'THE LOWEST',I3,' VIBRATIONS " "ARE NOT',/,' TO BE USED IN THE THERMO CALCULATION')", 0 }; static cilist io___104 = { 0, 6, 0, "(//10X,'SYSTEM IS A TRANSITION STAT" "E')", 0 }; static cilist io___105 = { 0, 6, 0, "(//10X,'SYSTEM IS A GROUND STATE')", 0 }; /* 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 */ /* ********************************************************************** */ /* FORCE CALCULATES THE FORCE CONSTANTS FOR THE MOLECULE, AND THE */ /* VIBRATIONAL FREQUENCIES. ISOTOPIC SUBSTITUTION IS ALLOWED. */ /* ********************************************************************** */ /* TEST GEOMETRY TO SEE IF IT IS OPTIMIZED */ time2 = -1e9; gmetry_(geom_1.geo, coord_1.coord); nvaold = geovar_1.nvar; i__1 = geovar_1.nvar; for (i = 1; i <= i__1; ++i) { locold[(i << 1) - 2] = geovar_1.loc[(i << 1) - 2]; /* L10: */ locold[(i << 1) - 1] = geovar_1.loc[(i << 1) - 1]; } geovar_1.nvar = 0; numat = 0; if (geokst_1.labels[0] != 99) { numat = 1; } i__1 = geokst_1.natoms; for (i = 2; i <= i__1; ++i) { if (geokst_1.labels[i - 1] == 99 || geokst_1.labels[i - 1] == 107) { goto L30; } ++numat; if (i == 2) { ilim = 1; } if (i == 3) { ilim = 2; } if (i > 3) { ilim = 3; } i__2 = ilim; for (j = 1; j <= i__2; ++j) { ++geovar_1.nvar; geovar_1.loc[(geovar_1.nvar << 1) - 2] = i; geovar_1.loc[(geovar_1.nvar << 1) - 1] = j; /* L20: */ xparam[geovar_1.nvar - 1] = geom_1.geo[j + i * 3 - 4]; } L30: ; } /* IF A RESTART, THEN TSCF AND TDER WILL BE FAULTY, THEREFORE SET TO -1 */ tscf = -1.; tder = -1.; prnt = i_indx(keywrd_1.keywrd, "RC=", 80L, 3L) == 0; debug = i_indx(keywrd_1.keywrd, "DFORCE", 80L, 6L) != 0; large = i_indx(keywrd_1.keywrd, "LARGE", 80L, 5L) != 0; bartel = i_indx(keywrd_1.keywrd, "NLLSQ", 80L, 5L) != 0; restrt = i_indx(keywrd_1.keywrd, "RESTART", 80L, 7L) != 0; time1 = second_(); if (restrt) { /* CHECK TO SEE IF CALCULATION IS IN NLLSQ OR FORCE. */ if (bartel) { goto L50; } /* CALCULATION IS IN FORCE */ goto L70; } compfg_(xparam, &c_true, &escf, &c_true, gradnt_1.grad, &c_false); if (prnt) { s_wsfe(&io___20); do_fio(&c__1, (char *)&escf, (ftnlen)sizeof(doublereal)); e_wsfe(); } time2 = second_(); tscf = time2 - time1; compfg_(xparam, &c_true, &escf1, &c_false, gradnt_1.grad, &c_true); time3 = second_(); tder = time3 - time2; if (prnt) { s_wsfe(&io___23); e_wsfe(); } l = 0; iu = 0; i__1 = geokst_1.natoms; for (i = 1; i <= i__1; ++i) { if (geokst_1.labels[i - 1] == 99) { goto L40; } ++l; il = iu + 1; if (i == 1) { iu = il - 1; } if (i == 2) { iu = il; } if (i == 3) { iu = il + 1; } if (i > 3) { iu = il + 2; } if (prnt) { s_wsfe(&io___27); do_fio(&c__1, (char *)&l, (ftnlen)sizeof(integer)); do_fio(&c__1, elemts_1.elemnt + (geokst_1.labels[i - 1] - 1 << 1), 2L); i__2 = iu; for (j = il; j <= i__2; ++j) { do_fio(&c__1, (char *)&gradnt_1.grad[j - 1], (ftnlen)sizeof( doublereal)); } e_wsfe(); } L40: ; } /* TEST SUM OF GRADIENTS */ gradnt_1.gnorm = sqrt(dot_(gradnt_1.grad, gradnt_1.grad, &geovar_1.nvar)); if (prnt) { s_wsfe(&io___28); do_fio(&c__1, (char *)&gradnt_1.gnorm, (ftnlen)sizeof(doublereal)); e_wsfe(); } if (gradnt_1.gnorm < 10.) { goto L60; } if (i_indx(keywrd_1.keywrd, " LET ", 80L, 5L) != 0) { s_wsfe(&io___29); e_wsfe(); goto L70; } s_wsfe(&io___30); e_wsfe(); L50: s_wsfe(&io___31); e_wsfe(); if (bartel) { s_wsfe(&io___32); e_wsfe(); nllsq_(xparam, &geovar_1.nvar); } else { s_wsfe(&io___33); e_wsfe(); flepo_(xparam, &geovar_1.nvar, &escf); } compfg_(xparam, &c_true, &escf, &c_true, gradnt_1.grad, &c_true); write_(&time1, &escf); s_wsfe(&io___34); do_fio(&c__1, (char *)&gradnt_1.gnorm, (ftnlen)sizeof(doublereal)); e_wsfe(); gmetry_(geom_1.geo, coord_1.coord); L60: /* NOW TO CALCULATE THE FORCE MATRIX */ /* CHECK OUT SYMMETRY */ L70: /* NEED TO ENSURE THAT XYZINT WILL WORK CORRECTLY BEFORE CALL */ /* TO DRC. */ l = 0; i__1 = geokst_1.natoms; for (i = 1; i <= i__1; ++i) { if (geokst_1.labels[i - 1] != 99) { ++l; geokst_1.labels[l - 1] = geokst_1.labels[i - 1]; } /* L80: */ } geokst_1.natoms = numat; xyzint_(coord_1.coord, &numat, geokst_1.na, geokst_1.nb, geokst_1.nc, & c_b42, geom_1.geo); gmetry_(geom_1.geo, coord_1.coord); if (i_indx(keywrd_1.keywrd, "THERMO", 80L, 6L) != 0 && gradnt_1.gnorm > 1.) { s_wsfe(&io___35); e_wsfe(); s_wsfe(&io___36); e_wsfe(); s_wsfe(&io___37); e_wsfe(); } if (tscf > 0.) { s_wsfe(&io___38); do_fio(&c__1, (char *)&tscf, (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___39); do_fio(&c__1, (char *)&tder, (ftnlen)sizeof(doublereal)); e_wsfe(); } if (geosym_1.ndep > 0) { s_wsfe(&io___40); e_wsfe(); geosym_1.ndep = 0; } if (prnt) { axis_(coord_1.coord, &numat, &a, &b, &c, &wtmol, &c__2, rot); } nvib = numat * 3 - 6; if (abs(c) < 1e-20) { ++nvib; } if (euler_1.id != 0) { nvib = numat * 3 - 3; } if (prnt) { s_wsfe(&io___47); e_wsfe(); s_wsfe(&io___48); e_wsfe(); } l = 0; i__1 = geokst_1.natoms; for (i = 1; i <= i__1; ++i) { if (geokst_1.labels[i - 1] == 99) { goto L90; } ++l; if (prnt) { s_wsfe(&io___49); do_fio(&c__1, (char *)&l, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&geokst_1.labels[i - 1], (ftnlen)sizeof( integer)); for (j = 1; j <= 3; ++j) { do_fio(&c__1, (char *)&coord_1.coord[j + l * 3 - 4], (ftnlen) sizeof(doublereal)); } e_wsfe(); } L90: ; } fmat_(fmatrx_1.fmatrx, &nvib, &tscf, &tder, deldip, &escf); /* THE FORCE MATRIX IS PRINTED AS AN ATOM-ATOM MATRIX RATHER THAN */ /* AS A 3N*3N MATRIX, AS THE 3N MATRIX IS VERY CONFUSING! */ ij = 0; iu = 0; i__1 = numat; for (i = 1; i <= i__1; ++i) { il = iu + 1; iu = il + 2; im1 = i - 1; ju = 0; i__2 = im1; for (j = 1; j <= i__2; ++j) { jl = ju + 1; ju = jl + 2; sum = 0.; i__3 = iu; for (ii = il; ii <= i__3; ++ii) { i__4 = ju; for (jj = jl; jj <= i__4; ++jj) { /* L100: */ /* Computing 2nd power */ d__1 = fmatrx_1.fmatrx[ii * (ii - 1) / 2 + jj - 1]; sum += d__1 * d__1; } } ++ij; /* L110: */ scrach_1.store[ij - 1] = sqrt(sum); } ++ij; /* L120: */ /* Computing 2nd power */ d__1 = fmatrx_1.fmatrx[il * (il + 1) / 2 - 1]; /* Computing 2nd power */ d__2 = fmatrx_1.fmatrx[(il + 1) * (il + 2) / 2 - 1]; /* Computing 2nd power */ d__3 = fmatrx_1.fmatrx[(il + 2) * (il + 3) / 2 - 1]; /* Computing 2nd power */ d__4 = fmatrx_1.fmatrx[(il + 1) * (il + 2) / 2 - 2]; /* Computing 2nd power */ d__5 = fmatrx_1.fmatrx[(il + 2) * (il + 3) / 2 - 3]; /* Computing 2nd power */ d__6 = fmatrx_1.fmatrx[(il + 2) * (il + 3) / 2 - 2]; scrach_1.store[ij - 1] = sqrt(d__1 * d__1 + d__2 * d__2 + d__3 * d__3 + (d__4 * d__4 + d__5 * d__5 + d__6 * d__6) * 2.); } if (debug) { s_wsfe(&io___58); e_wsfe(); i = -geovar_1.nvar; vecprt_(fmatrx_1.fmatrx, &i); } if (prnt) { s_wsfe(&io___59); e_wsfe(); vecprt_(scrach_1.store, &numat); } l = geovar_1.nvar * (geovar_1.nvar + 1) / 2; i__1 = l; for (i = 1; i <= i__1; ++i) { /* L130: */ scrach_1.store[i - 1] = fmatrx_1.fmatrx[i - 1]; } if (prnt) { axis_(coord_1.coord, &numat, &a, &b, &c, &sum, &c__0, rot); } if (prnt) { s_wsfe(&io___60); do_fio(&c__1, (char *)&escf, (ftnlen)sizeof(doublereal)); e_wsfe(); } if (large) { frame_(scrach_1.store, &numat, &c__0, shift); rsp_(scrach_1.store, &geovar_1.nvar, &geovar_1.nvar, vector_1.freq, vector_1.cnorml); i__1 = geovar_1.nvar; for (i = nvib + 1; i <= i__1; ++i) { j = (integer) ((vector_1.freq[i - 1] + 50.) * .01); /* L140: */ vector_1.freq[i - 1] -= j * 100; } if (prnt) { s_wsfe(&io___62); e_wsfe(); s_wsfe(&io___63); i__1 = geovar_1.nvar; for (i = nvib + 1; i <= i__1; ++i) { do_fio(&c__1, (char *)&vector_1.freq[i - 1], (ftnlen)sizeof( doublereal)); } e_wsfe(); s_wsfe(&io___64); e_wsfe(); s_wsfe(&io___65); i__1 = nvib; for (i = 1; i <= i__1; ++i) { do_fio(&c__1, (char *)&vector_1.freq[i - 1], (ftnlen)sizeof( doublereal)); } e_wsfe(); /* CONVERT TO WEIGHTED FMAT */ s_wsfe(&io___66); e_wsfe(); i = -geovar_1.nvar; matout_(vector_1.cnorml, vector_1.freq, &nvib, &i, &geovar_1.nvar) ; } } freqcy_(fmatrx_1.fmatrx, vector_1.freq, vector_1.cnorml, redmas, travel, & c_true); /* CALCULATE ZERO POINT ENERGY */ /* THESE CONSTANTS TAKEN FROM HANDBOOK OF CHEMISTRY AND PHYSICS 62ND ED. */ /* N AVOGADRO'S NUMBER = 6.022045*10**23 */ /* H PLANCK'S CONSTANT = 6.626176*10**(-34)JHZ */ /* C SPEED OF LIGHT = 2.99792458*10**10 CM/SEC */ /* CONST=0.5*N*H*C/(1000*4.184) */ const__ = .0014295718; sum = 0.; i__1 = geovar_1.nvar; for (i = 1; i <= i__1; ++i) { /* L150: */ sum += vector_1.freq[i - 1]; } sum *= const__; if (prnt) { s_wsfe(&io___70); do_fio(&c__1, (char *)&sum, (ftnlen)sizeof(doublereal)); e_wsfe(); } summ = 0.; i__1 = geovar_1.nvar; for (i = 1; i <= i__1; ++i) { sum1 = 1e-20; i__2 = geovar_1.nvar; for (j = 1; j <= i__2; ++j) { /* L160: */ /* Computing 2nd power */ d__1 = vector_1.cnorml[j + (i - 1) * geovar_1.nvar - 1]; sum1 += d__1 * d__1; } sum1 = 1. / sqrt(sum1); for (k = 1; k <= 3; ++k) { /* L170: */ gradnt_1.grad[k - 1] = 0.; } for (k = 1; k <= 3; ++k) { sum = 0.; i__2 = geovar_1.nvar; for (j = 1; j <= i__2; ++j) { /* L180: */ sum += vector_1.cnorml[j + (i - 1) * geovar_1.nvar - 1] * deldip[k + j * 3 - 4]; } summ += abs(sum); /* L190: */ trdip[k + i * 3 - 4] = sum * sum1; } /* Computing 2nd power */ d__1 = trdip[i * 3 - 3]; /* Computing 2nd power */ d__2 = trdip[i * 3 - 2]; /* Computing 2nd power */ d__3 = trdip[i * 3 - 1]; dipt[i - 1] = sqrt(d__1 * d__1 + d__2 * d__2 + d__3 * d__3); /* L200: */ } if (prnt) { s_wsfe(&io___76); i__1 = geovar_1.nvar - nvib; do_fio(&c__1, (char *)&i__1, (ftnlen)sizeof(integer)); e_wsfe(); s_wsfe(&io___77); e_wsfe(); } if (prnt && large) { s_wsfe(&io___78); e_wsfe(); nto6 = geovar_1.nvar / 6; nrem6 = geovar_1.nvar - nto6 * 6; iinc1 = -5; if (nto6 < 1) { goto L220; } i__1 = nto6; for (i = 1; i <= i__1; ++i) { s_wsfe(&io___82); e_wsfe(); iinc1 += 6; iinc2 = iinc1 + 5; s_wsfe(&io___84); i__2 = iinc2; for (j = iinc1; j <= i__2; ++j) { do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer)); } e_wsfe(); s_wsfe(&io___85); i__2 = iinc2; for (j = iinc1; j <= i__2; ++j) { do_fio(&c__1, (char *)&vector_1.freq[j - 1], (ftnlen)sizeof( doublereal)); } e_wsfe(); s_wsfe(&io___86); i__2 = iinc2; for (j = iinc1; j <= i__2; ++j) { do_fio(&c__1, (char *)&redmas[j - 1], (ftnlen)sizeof( doublereal)); } e_wsfe(); s_wsfe(&io___87); i__2 = iinc2; for (j = iinc1; j <= i__2; ++j) { do_fio(&c__1, (char *)&trdip[j * 3 - 3], (ftnlen)sizeof( doublereal)); } e_wsfe(); s_wsfe(&io___88); i__2 = iinc2; for (j = iinc1; j <= i__2; ++j) { do_fio(&c__1, (char *)&trdip[j * 3 - 2], (ftnlen)sizeof( doublereal)); } e_wsfe(); s_wsfe(&io___89); i__2 = iinc2; for (j = iinc1; j <= i__2; ++j) { do_fio(&c__1, (char *)&trdip[j * 3 - 1], (ftnlen)sizeof( doublereal)); } e_wsfe(); s_wsfe(&io___90); i__2 = iinc2; for (j = iinc1; j <= i__2; ++j) { do_fio(&c__1, (char *)&dipt[j - 1], (ftnlen)sizeof(doublereal) ); } e_wsfe(); /* L210: */ } L220: if (nrem6 < 1) { goto L230; } s_wsfe(&io___91); e_wsfe(); iinc1 += 6; iinc2 = iinc1 + (nrem6 - 1); s_wsfe(&io___92); i__1 = iinc2; for (j = iinc1; j <= i__1; ++j) { do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer)); } e_wsfe(); s_wsfe(&io___93); i__1 = iinc2; for (j = iinc1; j <= i__1; ++j) { do_fio(&c__1, (char *)&vector_1.freq[j - 1], (ftnlen)sizeof( doublereal)); } e_wsfe(); s_wsfe(&io___94); i__1 = iinc2; for (j = iinc1; j <= i__1; ++j) { do_fio(&c__1, (char *)&redmas[j - 1], (ftnlen)sizeof(doublereal)); } e_wsfe(); s_wsfe(&io___95); i__1 = iinc2; for (j = iinc1; j <= i__1; ++j) { do_fio(&c__1, (char *)&trdip[j * 3 - 3], (ftnlen)sizeof( doublereal)); } e_wsfe(); s_wsfe(&io___96); i__1 = iinc2; for (j = iinc1; j <= i__1; ++j) { do_fio(&c__1, (char *)&trdip[j * 3 - 2], (ftnlen)sizeof( doublereal)); } e_wsfe(); s_wsfe(&io___97); i__1 = iinc2; for (j = iinc1; j <= i__1; ++j) { do_fio(&c__1, (char *)&trdip[j * 3 - 1], (ftnlen)sizeof( doublereal)); } e_wsfe(); s_wsfe(&io___98); i__1 = iinc2; for (j = iinc1; j <= i__1; ++j) { do_fio(&c__1, (char *)&dipt[j - 1], (ftnlen)sizeof(doublereal)); } e_wsfe(); L230: ; } s_wsfe(&io___99); e_wsfe(); i = -geovar_1.nvar; matout_(vector_1.cnorml, vector_1.freq, &geovar_1.nvar, &i, & geovar_1.nvar); /* CARRY OUT IRC IF REQUESTED. */ if (i_indx(keywrd_1.keywrd, "IRC", 80L, 3L) + i_indx(keywrd_1.keywrd, "DRC", 80L, 3L) != 0) { i__1 = geovar_1.nvar; for (i = 1; i <= i__1; ++i) { geovar_1.loc[(i << 1) - 2] = 0; /* L240: */ geovar_1.loc[(i << 1) - 1] = 0; } geovar_1.nvar = nvaold; i__1 = geovar_1.nvar; for (i = 1; i <= i__1; ++i) { geovar_1.loc[(i << 1) - 2] = locold[(i << 1) - 2]; /* L250: */ geovar_1.loc[(i << 1) - 1] = locold[(i << 1) - 1]; } xyzint_(coord_1.coord, &numat, geokst_1.na, geokst_1.nb, geokst_1.nc, &c_b42, geom_1.geo); last_1.last = 1; drc_(vector_1.cnorml, vector_1.freq); s_stop("", 0L); } freqcy_(fmatrx_1.fmatrx, vector_1.freq, vector_1.cnorml, deldip, deldip, & c_false); s_wsfe(&io___100); e_wsfe(); matout_(vector_1.cnorml, vector_1.freq, &geovar_1.nvar, &i, & geovar_1.nvar); anavib_(coord_1.coord, vector_1.freq, dipt, &geovar_1.nvar, vector_1.cnorml, scrach_1.store, fmatrx_1.fmatrx, travel, redmas); if (i_indx(keywrd_1.keywrd, "THERMO", 80L, 6L) != 0) { gmetry_(geom_1.geo, coord_1.coord); i = i_indx(keywrd_1.keywrd, " ROT", 80L, 4L); if (i != 0) { sym = reada_(keywrd_1.keywrd, &i, 80L); } else { sym = 1.; } linear = (d__1 = a * b * c, abs(d__1)) < 1e-10; i = i_indx(keywrd_1.keywrd, " TRANS", 80L, 6L); /* "I" IS GOING TO MARK THE BEGINNING OF THE GENUINE VIBRATIONS. */ if (i != 0) { i = i_indx(keywrd_1.keywrd, " TRANS=", 80L, 7L); if (i != 0) { i = (integer) (reada_(keywrd_1.keywrd, &i, 80L) + 1); j = nvib - i + 1; s_wsfe(&io___103); i__1 = i - 1; do_fio(&c__1, (char *)&i__1, (ftnlen)sizeof(integer)); e_wsfe(); } else { s_wsfe(&io___104); e_wsfe(); i = 2; j = nvib - 1; } } else { s_wsfe(&io___105); e_wsfe(); i = 1; j = nvib; } thermo_(&a, &b, &c, &linear, &sym, &wtmol, &vector_1.freq[i - 1], &j, &escf); } return 0; } /* force_ */ #undef gr #undef keys