/* polar.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 koment[80], title[80]; } titles_; #define titles_1 titles_ struct { doublereal polvol[107]; } polvol_; #define polvol_1 polvol_ struct { char keywrd[80]; } keywrd_; #define keywrd_1 keywrd_ struct { integer natoms, labels[86], na[86], nb[86], nc[86]; } geokst_; #define geokst_1 geokst_ struct { integer nvar, loc[516] /* was [2][258] */, idumy; doublereal xparam[258]; } geovar_; #define geovar_1 geovar_ struct { doublereal time0; } time_; #define time_1 time_ struct { char elemnt[214]; } elemts_; #define elemts_1 elemts_ struct { doublereal core[107]; } core_; #define core_1 core_ struct { integer numat, nat[86], nfirst[86], nmidle[86], nlast[86], norbs, nelecs, nalpha, nbeta, nclose, nopen, ndumy; doublereal fract; } molkst_; #define molkst_1 molkst_ struct { integer ndep, locpar[258], idepfn[258], locdep[258]; } geosym_; #define geosym_1 geosym_ struct { doublereal geo[258] /* was [3][86] */; } geom_; #define geom_1 geom_ struct { integer last; } last_; #define last_1 last_ struct { doublereal coord[258] /* was [3][86] */; } coord_; #define coord_1 coord_ struct { doublereal tvec[9] /* was [3][3] */; integer idtvec; } euler_; #define euler_1 euler_ struct { doublereal rxyz[23220], xdumy[43344]; } scrach_; #define scrach_1 scrach_ /* Table of constant values */ static integer c__1 = 1; static logical c_true = TRUE_; static logical c_false = FALSE_; /* Subroutine */ int polar_(void) { /* Format strings */ static char fmt_10[] = "(\0021\002,20(\002*\002),\002 FINITE-FIELD POLAR" "IZABILITIES \002,20(\002*\002),//)"; static char fmt_20[] = "(/\002 ROTATION MATRIX FOR ORIENTATION OF MOLECU" "LE:\002/)"; static char fmt_30[] = "(5x,3f12.6)"; static char fmt_160[] = "(/\002 NEW TRANSLATION VECTOR:\002/,\002 \002,3" "(3f15.5))"; static char fmt_200[] = "(//\002 ENERGY OF \"REORIENTED\" SYSTEM WITHOUT" " FIELD:\002,f20.10)"; /* System generated locals */ integer i__1; doublereal d__1, d__2; /* Builtin functions */ /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen); integer i_indx(char *, char *, ftnlen, ftnlen), s_wsfe(cilist *), e_wsfe( void), do_fio(integer *, char *, ftnlen); /* Local variables */ static doublereal grad[258]; static integer mass; extern /* Subroutine */ int axis_(doublereal *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, integer *, doublereal *) ; static char type[7]; static doublereal sumw, sumx, sumy, sumz, heat0, a, b, c; static integer i, j, k, l; static doublereal atpol, tempv[9] /* was [3][3] */, dipvec[3]; extern /* Subroutine */ int compfg_(doublereal *, logical *, doublereal *, logical *, doublereal *, logical *), ffhpol_(doublereal *, doublereal *, doublereal *); static doublereal rotvec[9] /* was [3][3] */; extern /* Subroutine */ int gmetry_(doublereal *, doublereal *); static doublereal summax; static logical let; static doublereal sum; /* Fortran I/O blocks */ static cilist io___3 = { 0, 6, 0, fmt_10, 0 }; static cilist io___10 = { 0, 6, 0, fmt_20, 0 }; static cilist io___12 = { 0, 6, 0, fmt_30, 0 }; static cilist io___16 = { 0, 6, 0, "(//10X,'CARTESIAN COORDINATES ',/)", 0 }; static cilist io___17 = { 0, 6, 0, "(4X,'NO.',7X,'ATOM',9X,'X', " " 9X,'Y',9X,'Z',/)", 0 }; static cilist io___19 = { 0, 6, 0, "(I6,8X,A2,4X,3F10.4)", 0 }; static cilist io___21 = { 0, 6, 0, fmt_160, 0 }; static cilist io___29 = { 0, 6, 0, fmt_200, 0 }; /* ********************************************************************** */ /* POLAR SETS UP THE CALCULATION OF THE MOLECULAR ELECTRIC RESPONSE */ /* PROPERTIES BY FFHPOL. */ /* ********************************************************************** */ /* 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 */ s_copy(type, " MNDO ", 7L, 7L); let = i_indx(keywrd_1.keywrd, "LET", 80L, 3L) != 0; if (i_indx(keywrd_1.keywrd, "MINDO", 80L, 5L) != 0) { s_copy(type, "MINDO/3", 7L, 7L); } if (i_indx(keywrd_1.keywrd, "AM1", 80L, 3L) != 0) { s_copy(type, " AM1 ", 7L, 7L); } s_wsfe(&io___3); e_wsfe(); gmetry_(geom_1.geo, coord_1.coord); /* Orient the molecule with the moments of inertia. */ /* This is done to ensure a unique, reproduceable set of directions. */ /* If LET is specified, the input orientation will be used. */ if (! let) { mass = 1; axis_(coord_1.coord, &molkst_1.numat, &a, &b, &c, &sumw, &mass, rotvec); s_wsfe(&io___10); e_wsfe(); for (i = 1; i <= 3; ++i) { s_wsfe(&io___12); for (j = 1; j <= 3; ++j) { do_fio(&c__1, (char *)&rotvec[i + j * 3 - 4], (ftnlen)sizeof( doublereal)); } e_wsfe(); /* L40: */ } /* ROTATE ATOMS */ i__1 = molkst_1.numat; for (i = 1; i <= i__1; ++i) { for (j = 1; j <= 3; ++j) { sum = 0.; for (k = 1; k <= 3; ++k) { sum += coord_1.coord[k + i * 3 - 4] * rotvec[k + j * 3 - 4]; /* L50: */ } geom_1.geo[j + i * 3 - 4] = sum; /* L60: */ } /* L70: */ } i__1 = molkst_1.numat; for (i = 1; i <= i__1; ++i) { for (j = 1; j <= 3; ++j) { coord_1.coord[j + i * 3 - 4] = geom_1.geo[j + i * 3 - 4]; /* L80: */ } /* L90: */ } s_wsfe(&io___16); e_wsfe(); s_wsfe(&io___17); e_wsfe(); l = 0; i__1 = molkst_1.numat; for (i = 1; i <= i__1; ++i) { if (molkst_1.nat[i - 1] == 99 || molkst_1.nat[i - 1] == 107) { goto L100; } ++l; s_wsfe(&io___19); do_fio(&c__1, (char *)&l, (ftnlen)sizeof(integer)); do_fio(&c__1, elemts_1.elemnt + (molkst_1.nat[i - 1] - 1 << 1), 2L); for (j = 1; j <= 3; ++j) { do_fio(&c__1, (char *)&coord_1.coord[j + l * 3 - 4], (ftnlen) sizeof(doublereal)); } e_wsfe(); L100: ; } /* IF POLYMER, ROTATE TVEC */ if (euler_1.idtvec > 0) { i__1 = euler_1.idtvec; for (i = 1; i <= i__1; ++i) { for (j = 1; j <= 3; ++j) { sum = 0.; for (k = 1; k <= 3; ++k) { sum += euler_1.tvec[k + i * 3 - 4] * rotvec[k + j * 3 - 4]; /* L110: */ } tempv[j + i * 3 - 4] = sum; /* L120: */ } /* L130: */ } for (i = 1; i <= 3; ++i) { i__1 = euler_1.idtvec; for (j = 1; j <= i__1; ++j) { euler_1.tvec[i + j * 3 - 4] = tempv[i + j * 3 - 4]; /* L140: */ } /* L150: */ } s_wsfe(&io___21); i__1 = euler_1.idtvec; for (i = 1; i <= i__1; ++i) { for (j = 1; j <= 3; ++j) { do_fio(&c__1, (char *)&euler_1.tvec[j + i * 3 - 4], ( ftnlen)sizeof(doublereal)); } } e_wsfe(); } } last_1.last = 1; geokst_1.na[0] = 99; /* SET UP THE VARIABLES IN XPARAM AND LOC, THESE ARE IN CARTESIAN */ /* COORDINATES. */ geosym_1.ndep = 0; molkst_1.numat = 0; sumx = 0.; sumy = 0.; sumz = 0.; i__1 = geokst_1.natoms; for (i = 1; i <= i__1; ++i) { if (geokst_1.labels[i - 1] != 99 && geokst_1.labels[i - 1] != 107) { ++molkst_1.numat; geokst_1.labels[molkst_1.numat - 1] = geokst_1.labels[i - 1]; sumx += coord_1.coord[molkst_1.numat * 3 - 3]; sumy += coord_1.coord[molkst_1.numat * 3 - 2]; sumz += coord_1.coord[molkst_1.numat * 3 - 1]; for (j = 1; j <= 3; ++j) { /* L170: */ geom_1.geo[j + molkst_1.numat * 3 - 4] = coord_1.coord[j + molkst_1.numat * 3 - 4]; } } /* L180: */ } sumx /= molkst_1.numat; sumy /= molkst_1.numat; sumz /= molkst_1.numat; summax = 0.; atpol = 0.; i__1 = molkst_1.numat; for (i = 1; i <= i__1; ++i) { if (geokst_1.labels[i - 1] != 107) { atpol += polvol_1.polvol[geokst_1.labels[i - 1] - 1]; } geom_1.geo[i * 3 - 3] -= sumx; if (summax < (d__1 = geom_1.geo[i * 3 - 3], abs(d__1))) { summax = (d__2 = geom_1.geo[i * 3 - 3], abs(d__2)); } geom_1.geo[i * 3 - 2] -= sumy; if (summax < (d__1 = geom_1.geo[i * 3 - 2], abs(d__1))) { summax = (d__2 = geom_1.geo[i * 3 - 2], abs(d__2)); } geom_1.geo[i * 3 - 1] -= sumz; if (summax < (d__1 = geom_1.geo[i * 3 - 1], abs(d__1))) { summax = (d__2 = geom_1.geo[i * 3 - 1], abs(d__2)); } /* L190: */ } geovar_1.nvar = 0; geokst_1.natoms = molkst_1.numat; compfg_(geom_1.geo, &c_true, &heat0, &c_true, grad, &c_false); s_wsfe(&io___29); do_fio(&c__1, (char *)&heat0, (ftnlen)sizeof(doublereal)); e_wsfe(); ffhpol_(&heat0, &atpol, dipvec); return 0; } /* polar_ */