/* ffhpol.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 { doublereal core[107]; } core_; #define core_1 core_ struct { doublereal geo[258] /* was [3][86] */; } geom_; #define geom_1 geom_ struct { integer numat, nat[86], nfirst[86], nmidle[86], nlast[86], nors, nelecs, nalpha, nbeta, nclose, nopen, ndumy; doublereal fract; } molkst_; #define molkst_1 molkst_ struct { doublereal coord[258] /* was [3][86] */; } coord_; #define coord_1 coord_ struct { char keywrd[80]; } keywrd_; #define keywrd_1 keywrd_ struct { doublereal efield[3]; } field_; #define field_1 field_ struct { doublereal tvec[9] /* was [3][3] */; integer idtvec; } euler_; #define euler_1 euler_ /* Table of constant values */ static integer c__1 = 1; static logical c_true = TRUE_; static logical c_false = FALSE_; /* Subroutine */ int ffhpol_(doublereal *heat0, doublereal *atpol, doublereal *dipvec) { /* Initialized data */ static doublereal heat3m = 0.; static char axis[1*3] = "X" "Y" "Z"; static doublereal heat3p = 0.; static integer iptbd[6] = { 5,7,4,9,6,8 }; /* Format strings */ static char fmt_20[] = "(//\002 APPLIED ELECTRIC FIELD MAGNITUDE: \002,f" "15.5)"; static char fmt_30[] = "(//\002 ATOMIC CONTRIBUTION TO THE POLARIZABILIT" "Y: \002,f15.6)"; static char fmt_40[] = "(//,\002 ****** \002,a1,\002 DIRECTION *****\002" ",/)"; static char fmt_70[] = "(\002 ENERGIES AT: \002,5x,\002F\002,21x,\0022" "F\002,21x,\0023F\002,/)"; static char fmt_80[] = "(\002 + \002,3(f20.10,3x),/,\002 - \002,3(f2" "0.10,3x))"; static char fmt_120[] = "(/,\002 \002,12x,\002+,+\002,15x,\002+,-\002,15" "x,\002-,+\002,15x,\002-,-\002)"; static char fmt_130[] = "(\002 E \002,4f15.6)"; static char fmt_140[] = "(\002 2E \002,4f15.6)"; static char fmt_170[] = "(//,\002 \002,30(\002*\002),\002 DIPOLE \002," "30(\002*\002),//)"; static char fmt_180[] = "(21x,\002E4\002,13x,\002DIP\002,/)"; static char fmt_190[] = "(5x,a1,7x,2f15.6)"; static char fmt_200[] = "(//\002 MAGNITUDE: \002,2f15.6,\002 (A.U.)" "\002,/,\002 \002,12x,2f15.6,\002 (DEBYE)\002)"; static char fmt_210[] = "(//,\002 \002,30(\002*\002),\002 POLARIZABILI" "TY \002,20(\002*\002),//)"; static char fmt_220[] = "(/\002 E4 POLARIZABILITY TENSOR:\002)"; static char fmt_230[] = "(/\002 DIP POLARIZABILITY TENSOR:\002)"; static char fmt_240[] = "(//,\002 AVERAGE POLARIZABILITY:\002,8x,\002E" "4\002,13x,\002DIP\002,/,\002 \002,24x,2f15.6,\002 A.U.\002,/" ",\002 \002,24x,2f15.6,\002 ANG.**3\002,/,\002 \002,24x,2(1pd15." "6),\002 ESU\002)"; static char fmt_250[] = "(//,\002 \002,30(\002*\002),\002 SECOND-ORDER" " \002,25(\002*\002),//)"; static char fmt_260[] = "(\002 COMPONENT\002,12x,\002E4\002,13x,\002DI" "P\002,/)"; static char fmt_270[] = "(\002 \002,5x,a4,5x,2f15.6)"; static char fmt_280[] = "(/)"; static char fmt_290[] = "(\002 \002,6x,a2,6x,2f15.6)"; static char fmt_300[] = "(\002 \002,4x,\002B(AU)\002,5x,2f15.6,/,\002" " \002,4x,\002B(ESU)\002,4x,2f15.6,3x,\002(X10-30)\002)"; static char fmt_310[] = "(//\002 \002,30(\002*\002),\002 THIRD-ORDER " "\002,25(\002*\002),//)"; static char fmt_320[] = "(\002 \002,17x,\002E4\002,13x,\002DIP\002,/)"; static char fmt_330[] = "(5x,a4,2f15.6)"; static char fmt_340[] = "(//\002 GAMMA = \002,1pd15.6,1pd15.6,\002 A.U" ".\002/,\002 \002,8x,1pd15.6,1pd15.6,\002 ESU (X10-36)\002)"; /* System generated locals */ integer i__1, i__2; doublereal d__1; /* Builtin functions */ integer i_indx(char *, char *, ftnlen, ftnlen), s_wsfe(cilist *), do_fio( integer *, char *, ftnlen), e_wsfe(void); double sqrt(doublereal); /* Local variables */ static doublereal grad, bkii, bjii, bijj, diim, gdip; static integer idip; static doublereal diip, eigs[3], bdmu, hnuc, avga3, dipe4[3], dip1m[3], dip2m[3], dip1p[3], dip2p[3], b4esu; static integer i; static logical debug; static doublereal heate[6] /* was [3][2] */; static logical large; static integer nbdip; static doublereal efval, dipdp[3]; static integer ngdip; static doublereal giijj; static integer nbcnt; static doublereal bdesu, hnucj, autoa; static integer ngcnt; static doublereal aterm, eterm, bterm, gterm; static integer i3; static doublereal betae4[9], avga3d, dipe4d, gamme4[6], heat1m, heat2m, heat1p, heat2p, apole4[6], avgpe4, dipe4t, ae, be, ge; static integer id, jd, kd; static doublereal betadp[9], gamdes, gamdip, dipdpd, gammdp[6]; extern /* Subroutine */ int compfg_(doublereal *, logical *, doublereal *, logical *, doublereal *, logical *), dipind_(doublereal *); static doublereal avgesd, apoldp[6], autodb, avgpdp, dipdpt; static logical poldip; static doublereal autokc, avgesu, gamval, gamesu; extern /* Subroutine */ int vecprt_(doublereal *, integer *); static doublereal vectrs[9]; extern /* Subroutine */ int matout_(doublereal *, doublereal *, integer *, integer *, integer *); static doublereal autovm, bx4, by4, bz4; static integer nbd; static doublereal aki, aij, sfe, dmm, bxd, dpm, hmm, dpp, hpm, dmu, hmp, hpp, dmp, byd; static integer ivl; static doublereal bzd; static integer kvl; extern /* Subroutine */ int rsp_(doublereal *, integer *, integer *, doublereal *, doublereal *); static integer idm1; static doublereal h2mm, h2pm, h2mp, b4mu, h2pp; /* Fortran I/O blocks */ static cilist io___17 = { 0, 6, 0, fmt_20, 0 }; static cilist io___20 = { 0, 6, 0, fmt_30, 0 }; static cilist io___22 = { 0, 6, 0, fmt_40, 0 }; static cilist io___37 = { 0, 6, 0, fmt_70, 0 }; static cilist io___38 = { 0, 6, 0, fmt_80, 0 }; static cilist io___69 = { 0, 6, 0, fmt_120, 0 }; static cilist io___70 = { 0, 6, 0, fmt_130, 0 }; static cilist io___75 = { 0, 6, 0, fmt_140, 0 }; static cilist io___84 = { 0, 6, 0, fmt_170, 0 }; static cilist io___89 = { 0, 6, 0, fmt_180, 0 }; static cilist io___90 = { 0, 6, 0, fmt_190, 0 }; static cilist io___91 = { 0, 6, 0, fmt_190, 0 }; static cilist io___92 = { 0, 6, 0, fmt_190, 0 }; static cilist io___93 = { 0, 6, 0, fmt_200, 0 }; static cilist io___94 = { 0, 6, 0, fmt_210, 0 }; static cilist io___95 = { 0, 6, 0, fmt_220, 0 }; static cilist io___102 = { 0, 6, 0, fmt_230, 0 }; static cilist io___106 = { 0, 6, 0, fmt_240, 0 }; static cilist io___107 = { 0, 6, 0, fmt_250, 0 }; static cilist io___118 = { 0, 6, 0, fmt_260, 0 }; static cilist io___119 = { 0, 6, 0, fmt_270, 0 }; static cilist io___120 = { 0, 6, 0, fmt_270, 0 }; static cilist io___121 = { 0, 6, 0, fmt_270, 0 }; static cilist io___122 = { 0, 6, 0, fmt_270, 0 }; static cilist io___123 = { 0, 6, 0, fmt_270, 0 }; static cilist io___124 = { 0, 6, 0, fmt_270, 0 }; static cilist io___125 = { 0, 6, 0, fmt_270, 0 }; static cilist io___126 = { 0, 6, 0, fmt_270, 0 }; static cilist io___127 = { 0, 6, 0, fmt_270, 0 }; static cilist io___128 = { 0, 6, 0, fmt_280, 0 }; static cilist io___129 = { 0, 6, 0, fmt_290, 0 }; static cilist io___130 = { 0, 6, 0, fmt_290, 0 }; static cilist io___131 = { 0, 6, 0, fmt_290, 0 }; static cilist io___132 = { 0, 6, 0, fmt_280, 0 }; static cilist io___133 = { 0, 6, 0, fmt_300, 0 }; static cilist io___134 = { 0, 6, 0, fmt_310, 0 }; static cilist io___139 = { 0, 6, 0, fmt_320, 0 }; static cilist io___140 = { 0, 6, 0, fmt_330, 0 }; static cilist io___141 = { 0, 6, 0, fmt_330, 0 }; static cilist io___142 = { 0, 6, 0, fmt_330, 0 }; static cilist io___143 = { 0, 6, 0, fmt_330, 0 }; static cilist io___144 = { 0, 6, 0, fmt_330, 0 }; static cilist io___145 = { 0, 6, 0, fmt_330, 0 }; static cilist io___146 = { 0, 6, 0, fmt_340, 0 }; /* *********************************************************************** */ /* SUBROUTINE FOR THE FINITE FIELD CALCULATION OF ELECTRIC RESPONSE */ /* PROPERTIES (DIPOLE MOMENT, POLARIZABILITY, AND 1ST AND 2ND */ /* HYPERPOLARIZABILITY. */ /* HENRY A. KURTZ, DEPARTMENT OF CHEMISTRY */ /* MEMPHIS STATE UNIVERSITY */ /* MEMPHIS, TN 38152 */ /* *********************************************************************** */ /* 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 */ /* DIPE4 AND DIPDP HOLD THE CALCULATED DIPOLE MOMENTS */ /* APOLE4 AND APOLDP HOLD THE POLARIZABILITY TENSOR AS */ /* A PACKED ARRAY XX,XY,YY,XZ,YZ,ZZ */ /* BETAE4 AND BETAEP HOLD THE FIRST HYPERPOLARIZABILITY */ /* 1. XXX */ /* 2. YYY 6. YXX */ /* 3. ZZZ 7. YZZ */ /* 4. XYY 8. ZXX */ /* 5. XZZ 9. ZYY */ /* Parameter adjustments */ --dipvec; /* Function Body */ /* Energy: a.u. to kcal/mole */ autokc = 627.50595269999997; /* Length: a.u. to Angstrom */ autoa = .529177; /* Dipole: a.u. to debye */ autodb = 2.541563; /* Electric Field: a.u. to volt/meter */ autovm = 51.4257; nbdip = 1; ngdip = 4; nbcnt = 4; ngcnt = 4; large = i_indx(keywrd_1.keywrd, "LARGE", 80L, 5L) != 0; debug = i_indx(keywrd_1.keywrd, "DEBUG", 80L, 5L) != 0; /* FIELD STRENGTH IN A.U. */ efval = .001; idip = 1; /* # READ (7,10) EFVAL,IDIP */ /* L10: */ s_wsfe(&io___17); do_fio(&c__1, (char *)&efval, (ftnlen)sizeof(doublereal)); e_wsfe(); poldip = FALSE_; if (idip != 0) { poldip = TRUE_; } sfe = 1. / efval; s_wsfe(&io___20); d__1 = *atpol * 6.74834f; do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof(doublereal)); e_wsfe(); /* ....................................................................... */ /* CALCULATE THE POLARIZABILITY AND HYPERPOLARIZABILITIES ALONG */ /* THE THREE PRINCIPLE AXES. (THESE AXES DEPEND ON YOUR ARBITRARY */ /* ORIENTATION AND MAY NOT BE THE TRUE PRINCIPLE AXES.) */ /* ....................................................................... */ for (id = 1; id <= 3; ++id) { if (debug) { s_wsfe(&io___22); do_fio(&c__1, axis + (id - 1), 1L); e_wsfe(); } /* ZERO THE FIELD */ for (i = 1; i <= 3; ++i) { field_1.efield[i - 1] = 0.; /* L50: */ } hnuc = 0.; i__1 = molkst_1.numat; for (i = 1; i <= i__1; ++i) { hnuc += efval * geom_1.geo[id + i * 3 - 4] * core_1.core[ molkst_1.nat[i - 1] - 1] * autovm; /* L60: */ } hnuc *= 23.061; /* +E(ID) */ field_1.efield[id - 1] = efval; compfg_(geom_1.geo, &c_true, &heat1p, &c_true, &grad, &c_false); dipind_(dip1p); diip = dip1p[id - 1]; /* -E(ID) */ field_1.efield[id - 1] = -efval; compfg_(geom_1.geo, &c_true, &heat1m, &c_true, &grad, &c_false); dipind_(dip1m); diim = dip1m[id - 1]; /* +2E(ID) */ field_1.efield[id - 1] = efval * 2.; compfg_(geom_1.geo, &c_true, &heat2p, &c_true, &grad, &c_false); dipind_(dip2p); /* -2E(ID) */ field_1.efield[id - 1] = efval * -2.; compfg_(geom_1.geo, &c_true, &heat2m, &c_true, &grad, &c_false); dipind_(dip2m); /* CORRECT FOR ELECTRIC FIELD - NUCLEAR INTERACTIONS */ heat1p += hnuc; heate[id - 1] = heat1p; heat1m -= hnuc; heate[id + 2] = heat1m; heat2p += hnuc * 2.; heat2m -= hnuc * 2.; if (debug) { s_wsfe(&io___37); e_wsfe(); s_wsfe(&io___38); do_fio(&c__1, (char *)&heat1p, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&heat2p, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&heat3p, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&heat1m, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&heat2m, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&heat3m, (ftnlen)sizeof(doublereal)); e_wsfe(); } /* DIPOLE */ eterm = (heat2p - heat2m) * .083333333333333321 - (heat1p - heat1m) * .66666666666666661; dipe4[id - 1] = eterm * sfe / autokc; /* ALPHA */ ivl = id * (id + 1) / 2; eterm = *heat0 * 2.5 - (heat1p + heat1m) * 1.3333333333333333 + ( heat2p + heat2m) * .083333333333333321; apole4[ivl - 1] = eterm * sfe * sfe / autokc + *atpol * 6.74834f; /* BETA */ eterm = heat1p - heat1m - (heat2p - heat2m) * .5; betae4[id - 1] = eterm * sfe * sfe * sfe / autokc; /* GAMMA */ eterm = (heat1p + heat1m) * 4. - (heat2p + heat2m) - *heat0 * 6.; gamme4[id - 1] = eterm * sfe * sfe * sfe * sfe / autokc; /* DIPOLE CALCULATIONS */ dmu = (dip1p[id - 1] + dip1m[id - 1]) * .66666666666666661 - (dip2p[ id - 1] + dip2m[id - 1]) * .16666666666666665; dipdp[id - 1] = dmu / autodb; ae = (dip1p[id - 1] - dip1m[id - 1]) * .66666666666666661 - (dip2p[id - 1] - dip2m[id - 1]) * .083333333333333321; apoldp[ivl - 1] = ae * sfe / autodb; be = (dip2p[id - 1] + dip2m[id - 1] - dip1p[id - 1] - dip1m[id - 1]) * .3333333333333333; betadp[id - 1] = be * sfe * sfe / autodb; ge = (dip2p[id - 1] - dip2m[id - 1]) * .5 - (dip1p[id - 1] - dip1m[id - 1]); gammdp[id - 1] = ge * sfe * sfe * sfe / autodb; for (kd = 1; kd <= 3; ++kd) { if (kd < id) { kvl = id * (id - 1) / 2 + kd; aki = (dip1p[kd - 1] - dip1m[kd - 1]) * .66666666666666661 - ( dip2p[kd - 1] - dip2m[kd - 1]) * .083333333333333321; apoldp[kvl - 1] = aki * sfe / autodb; } if (kd != id) { bkii = (dip2p[kd - 1] + dip2m[kd - 1] - dip1p[kd - 1] - dip1m[ kd - 1]) * .3333333333333333; nbd = iptbd[nbdip - 1]; betadp[nbd - 1] = bkii * sfe * sfe / autodb; ++nbdip; } /* L90: */ } /* ................................................................... .... */ /* NOW CALCULATE THE OFF AXIS RESULTS. */ /* ................................................................... .... */ idm1 = id - 1; i__1 = idm1; for (jd = 1; jd <= i__1; ++jd) { hnucj = 0.; i__2 = molkst_1.numat; for (i = 1; i <= i__2; ++i) { hnucj += efval * geom_1.geo[jd + i * 3 - 4] * core_1.core[ molkst_1.nat[i - 1] - 1] * 51.4257f; /* L100: */ } hnucj *= 23.061f; for (i = 1; i <= 3; ++i) { field_1.efield[i - 1] = 0.; /* L110: */ } /* DIAGONAL FIELDS WITH COMPONENTS EQUAL TO EFVAL */ field_1.efield[id - 1] = efval; field_1.efield[jd - 1] = efval; compfg_(geom_1.geo, &c_true, &hpp, &c_true, &grad, &c_false); dipind_(dip1p); dpp = dip1p[id - 1]; field_1.efield[jd - 1] = -efval; compfg_(geom_1.geo, &c_true, &hpm, &c_true, &grad, &c_false); dipind_(dip1p); dpm = dip1p[id - 1]; field_1.efield[id - 1] = -efval; compfg_(geom_1.geo, &c_true, &hmm, &c_true, &grad, &c_false); dipind_(dip1p); dmm = dip1p[id - 1]; field_1.efield[jd - 1] = efval; compfg_(geom_1.geo, &c_true, &hmp, &c_true, &grad, &c_false); dipind_(dip1p); dmp = dip1p[id - 1]; hpp = hpp + hnuc + hnucj; hpm = hpm + hnuc - hnucj; hmm = hmm - hnuc - hnucj; hmp = hmp - hnuc + hnucj; if (debug) { s_wsfe(&io___69); e_wsfe(); s_wsfe(&io___70); do_fio(&c__1, (char *)&hpp, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&hpm, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&hmp, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&hmm, (ftnlen)sizeof(doublereal)); e_wsfe(); } /* DIAGONAL FIELDS WITH COMPONENTS EQUAL TO 2*EFVAL */ field_1.efield[id - 1] = efval * 2.; field_1.efield[jd - 1] = efval * 2.; compfg_(geom_1.geo, &c_true, &h2pp, &c_true, &grad, &c_false); field_1.efield[jd - 1] = -efval * 2.; compfg_(geom_1.geo, &c_true, &h2pm, &c_true, &grad, &c_false); field_1.efield[id - 1] = -efval * 2.; compfg_(geom_1.geo, &c_true, &h2mm, &c_true, &grad, &c_false); field_1.efield[jd - 1] = efval * 2.; compfg_(geom_1.geo, &c_true, &h2mp, &c_true, &grad, &c_false); h2pp += (hnuc + hnucj) * 2.; h2pm += (hnuc - hnucj) * 2.; h2mm -= (hnuc + hnucj) * 2.; h2mp -= (hnuc - hnucj) * 2.; if (debug) { s_wsfe(&io___75); do_fio(&c__1, (char *)&h2pp, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&h2pm, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&h2mp, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&h2mm, (ftnlen)sizeof(doublereal)); e_wsfe(); } aterm = (h2pp - h2pm - h2mp + h2mm) * .02083333333333333 - (hpp - hpm - hmp + hmm) * .3333333333333333; aij = aterm * sfe * sfe / autokc; ivl = id * (id - 1) / 2 + jd; apole4[ivl - 1] = aij; bterm = (hmm - hpp + hpm - hmp) * .5 + heate[jd - 1] - heate[jd + 2]; bjii = bterm * sfe * sfe * sfe / autokc; betae4[nbcnt - 1] = bjii; ++nbcnt; bterm = (hmm - hpp + hmp - hpm) * .5 + heate[id - 1] - heate[id + 2]; bijj = bterm * sfe * sfe * sfe / autokc; betae4[nbcnt - 1] = bijj; ++nbcnt; gterm = -(hpp + hmm + hpm + hmp) - *heat0 * 4. + (heate[id - 1] + heate[id + 2]) * 2. + (heate[jd - 1] + heate[jd + 2]) * 2.; giijj = gterm * sfe * sfe * sfe * sfe / autokc; gamme4[ngcnt - 1] = giijj; gdip = (dpp - dmp + dpm - dmm) * .5 - (diip - diim); gammdp[ngcnt - 1] = gdip * sfe * sfe * sfe / autodb; ++ngcnt; /* L150: */ } /* L160: */ } /* ----------------------------------------------------------------------- */ /* SUMMARIZE THE RESULTS */ /* ----------------------------------------------------------------------- */ s_wsfe(&io___84); e_wsfe(); dipe4t = sqrt(dipe4[0] * dipe4[0] + dipe4[1] * dipe4[1] + dipe4[2] * dipe4[2]); dipe4d = dipe4t * autodb; dipdpt = sqrt(dipdp[0] * dipdp[0] + dipdp[1] * dipdp[1] + dipdp[2] * dipdp[2]); dipdpd = dipdpt * autodb; s_wsfe(&io___89); e_wsfe(); s_wsfe(&io___90); do_fio(&c__1, "X", 1L); do_fio(&c__1, (char *)&dipe4[0], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&dipdp[0], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___91); do_fio(&c__1, "Y", 1L); do_fio(&c__1, (char *)&dipe4[1], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&dipdp[1], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___92); do_fio(&c__1, "X", 1L); do_fio(&c__1, (char *)&dipe4[2], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&dipdp[2], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___93); do_fio(&c__1, (char *)&dipe4t, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&dipdpt, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&dipe4d, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&dipdpd, (ftnlen)sizeof(doublereal)); e_wsfe(); /* FIND EIGENVALUES AND EIGENVECTORS OF POLARIZATION MATRIX. */ s_wsfe(&io___94); e_wsfe(); s_wsfe(&io___95); e_wsfe(); i3 = 3; i__1 = -i3; vecprt_(apole4, &i__1); rsp_(apole4, &i3, &i3, eigs, vectrs); matout_(vectrs, eigs, &i3, &i3, &i3); avgpe4 = (eigs[0] + eigs[1] + eigs[2]) / 3.; avga3 = avgpe4 * .14818; avgesu = avgpe4 * 2.96352e-25; s_wsfe(&io___102); e_wsfe(); i__1 = -i3; vecprt_(apoldp, &i__1); rsp_(apoldp, &i3, &i3, eigs, vectrs); matout_(vectrs, eigs, &i3, &i3, &i3); avgpdp = (eigs[0] + eigs[1] + eigs[2]) / 3.; avga3d = avgpdp * .14818; avgesd = avgpdp * 2.96352e-25; s_wsfe(&io___106); do_fio(&c__1, (char *)&avgpe4, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&avgpdp, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&avga3, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&avga3d, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&avgesu, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&avgesd, (ftnlen)sizeof(doublereal)); e_wsfe(); /* CALCULATE "EXPERIMENTAL" HYPERPOLARIZABILITIES */ /* 8.65710D-33 is a.u. to e.s.u. conversion */ s_wsfe(&io___107); e_wsfe(); bx4 = (betae4[0] + betae4[3] + betae4[5]) * .6; by4 = (betae4[1] + betae4[4] + betae4[7]) * .6; bz4 = (betae4[2] + betae4[6] + betae4[8]) * .6; b4mu = (bx4 * dipe4[0] + by4 * dipe4[1] + bz4 * dipe4[2]) / dipe4t; b4esu = b4mu * .0086571; bxd = (betadp[0] + betadp[3] + betadp[5]) * .6; byd = (betadp[1] + betadp[4] + betadp[7]) * .6; bzd = (betadp[2] + betadp[6] + betadp[8]) * .6; bdmu = (bxd * dipdp[0] + byd * dipdp[1] + bzd * dipdp[2]) / dipdpt; bdesu = bdmu * .0086571; s_wsfe(&io___118); e_wsfe(); s_wsfe(&io___119); do_fio(&c__1, "XXX", 3L); do_fio(&c__1, (char *)&betae4[0], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&betadp[0], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___120); do_fio(&c__1, "XYY", 3L); do_fio(&c__1, (char *)&betae4[3], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&betadp[3], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___121); do_fio(&c__1, "XZZ", 3L); do_fio(&c__1, (char *)&betae4[5], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&betadp[5], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___122); do_fio(&c__1, "YYY", 3L); do_fio(&c__1, (char *)&betae4[1], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&betadp[1], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___123); do_fio(&c__1, "YXX", 3L); do_fio(&c__1, (char *)&betae4[4], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&betadp[4], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___124); do_fio(&c__1, "YZZ", 3L); do_fio(&c__1, (char *)&betae4[7], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&betadp[7], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___125); do_fio(&c__1, "ZZZ", 3L); do_fio(&c__1, (char *)&betae4[2], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&betadp[2], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___126); do_fio(&c__1, "ZXX", 3L); do_fio(&c__1, (char *)&betae4[6], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&betadp[6], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___127); do_fio(&c__1, "ZYY", 3L); do_fio(&c__1, (char *)&betae4[8], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&betadp[8], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___128); e_wsfe(); s_wsfe(&io___129); do_fio(&c__1, "BX", 2L); do_fio(&c__1, (char *)&bx4, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&bxd, (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___130); do_fio(&c__1, "BY", 2L); do_fio(&c__1, (char *)&by4, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&byd, (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___131); do_fio(&c__1, "BZ", 2L); do_fio(&c__1, (char *)&bz4, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&bzd, (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___132); e_wsfe(); s_wsfe(&io___133); do_fio(&c__1, (char *)&b4mu, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&bdmu, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&b4esu, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&bdesu, (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___134); e_wsfe(); gamval = gamme4[0] + gamme4[1] + gamme4[2]; gamval += (gamme4[3] + gamme4[4] + gamme4[5]) * 2.; gamval /= 5.; /* 5.05116D-40 is the a.u. to e.s.u. conversion */ gamesu = gamval * 5.05116e-4; gamdip = gammdp[0] + gammdp[1] + gammdp[2]; gamdip += (gammdp[3] + gammdp[4] + gammdp[5]) * 2.; gamdip /= 5.; gamdes = gamdip * 5.05116e-4; s_wsfe(&io___139); e_wsfe(); s_wsfe(&io___140); do_fio(&c__1, "XXXX", 4L); do_fio(&c__1, (char *)&gamme4[0], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&gammdp[0], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___141); do_fio(&c__1, "YYYY", 4L); do_fio(&c__1, (char *)&gamme4[1], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&gammdp[1], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___142); do_fio(&c__1, "ZZZZ", 4L); do_fio(&c__1, (char *)&gamme4[2], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&gammdp[2], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___143); do_fio(&c__1, "XXYY", 4L); do_fio(&c__1, (char *)&gamme4[3], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&gammdp[3], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___144); do_fio(&c__1, "XXZZ", 4L); do_fio(&c__1, (char *)&gamme4[4], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&gammdp[4], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___145); do_fio(&c__1, "YYZZ", 4L); do_fio(&c__1, (char *)&gamme4[5], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&gammdp[5], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___146); do_fio(&c__1, (char *)&gamval, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&gamdip, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&gamesu, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&gamdes, (ftnlen)sizeof(doublereal)); e_wsfe(); return 0; } /* ffhpol_ */