/* axis.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 { doublereal atmass[86]; } atmass_; #define atmass_1 atmass_ /* Table of constant values */ static integer c__1 = 1; static integer c__3 = 3; /* Subroutine */ int axis_(doublereal *coord, integer *numat, doublereal *a, doublereal *b, doublereal *c, doublereal *sumw, integer *mass, doublereal *evec) { /* Initialized data */ static doublereal t[6] = { 0.,0.,0.,0.,0.,0. }; static logical first = TRUE_; /* System generated locals */ integer i__1; doublereal d__1, d__2; /* Builtin functions */ integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void), i_indx(char *, char *, ftnlen, ftnlen); /* Local variables */ static integer i, j; static doublereal x[86], y[86], z[86], sumwx, sumwy, sumwz, const1, const2, weight, xyzmom[3], eig[3], rot[3]; extern /* Subroutine */ int rsp_(doublereal *, integer *, integer *, doublereal *, doublereal *); static doublereal sum; /* Fortran I/O blocks */ static cilist io___10 = { 0, 6, 0, "(/10X,'MOLECULAR WEIGHT =',F8.2,/)", 0 }; static cilist io___15 = { 0, 6, 0, "(//10X,' PRINCIPAL MOMENTS OF INERTI" "A IN CM(-1)',/)", 0 }; static cilist io___18 = { 0, 6, 0, "(10X,'A =',F12.6,' B =',F12.6, " " ' C =',F12.6,/)", 0 }; static cilist io___19 = { 0, 6, 0, "(//10X,' PRINCIPAL MOMENTS OF INERTI" "A IN ', 'UNITS OF 10**(-40)*GRAM-CM**2',/)", 0 }; static cilist io___20 = { 0, 6, 0, "(10X,'A =',F12.6,' B =',F12.6, " " ' C =',F12.6,/)", 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 */ /* *********************************************************************** */ /* AXIS CALCULATES THE THREE MOMENTS OF INERTIA AND THE MOLECULAR */ /* WEIGHT. THE MOMENTS OF INERTIA ARE RETURNED IN A, B, AND C. */ /* THE MOLECULAR WEIGHT IN SUMW. */ /* THE UNITS OF INERTIA ARE 10**(-40)GRAM-CM**2, */ /* AND MOL.WEIGHT IN ATOMIC-MASS-UNITS. (AMU'S) */ /* *********************************************************************** */ /* Parameter adjustments */ coord -= 4; evec -= 4; /* Function Body */ /* *********************************************************************** */ /* CONST1 = 10**40/(N*A*A) */ /* N = AVERGADRO'S NUMBER */ /* A = CM IN AN ANGSTROM */ /* 10**40 IS TO ALLOW UNITS TO BE 10**(-40)GRAM-CM**2 */ /* *********************************************************************** */ const1 = 1.66053; /* *********************************************************************** */ /* CONST2 = CONVERSION FACTOR FROM ANGSTROM-AMU TO CM**(-1) */ /* = (PLANCK'S CONSTANT*N*10**16)/(8*PI*PI*C) */ /* = 6.62618*10**(-27)[ERG-SEC]*6.02205*10**23*10**16/ */ /* (8*(3.1415926535)**2*2.997925*10**10[CM/SEC]) */ /* *********************************************************************** */ const2 = 16.8576522; /* FIRST WE CENTRE THE MOLECULE ABOUT THE CENTRE OF GRAVITY, */ /* THIS DEPENDS ON THE ISOTOPIC MASSES, AND THE CARTESIAN GEOMETRY. */ *sumw = 1e-20; sumwx = 0.; sumwy = 0.; sumwz = 0.; weight = 1.; i__1 = *numat; for (i = 1; i <= i__1; ++i) { if (*mass > 0) { weight = atmass_1.atmass[i - 1]; } *sumw += weight; sumwx += weight * coord[i * 3 + 1]; sumwy += weight * coord[i * 3 + 2]; /* L10: */ sumwz += weight * coord[i * 3 + 3]; } if (*mass > 0 && first) { s_wsfe(&io___10); d__1 = min(99999.99,*sumw); do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof(doublereal)); e_wsfe(); } sumwx /= *sumw; sumwy /= *sumw; sumwz /= *sumw; i__1 = *numat; for (i = 1; i <= i__1; ++i) { x[i - 1] = coord[i * 3 + 1] - sumwx; y[i - 1] = coord[i * 3 + 2] - sumwy; /* L20: */ z[i - 1] = coord[i * 3 + 3] - sumwz; } /* *********************************************************************** */ /* MATRIX FOR MOMENTS OF INERTIA IS OF FORM */ /* | Y**2+Z**2 | */ /* | -Y*X Z**2+X**2 | -I =0 */ /* | -Z*X -Z*Y X**2+Y**2 | */ /* *********************************************************************** */ for (i = 1; i <= 6; ++i) { /* L30: */ t[i - 1] = i * 1e-10; } i__1 = *numat; for (i = 1; i <= i__1; ++i) { if (*mass > 0) { weight = atmass_1.atmass[i - 1]; } /* Computing 2nd power */ d__1 = y[i - 1]; /* Computing 2nd power */ d__2 = z[i - 1]; t[0] += weight * (d__1 * d__1 + d__2 * d__2); t[1] -= weight * x[i - 1] * y[i - 1]; /* Computing 2nd power */ d__1 = z[i - 1]; /* Computing 2nd power */ d__2 = x[i - 1]; t[2] += weight * (d__1 * d__1 + d__2 * d__2); t[3] -= weight * z[i - 1] * x[i - 1]; t[4] -= weight * y[i - 1] * z[i - 1]; /* L40: */ /* Computing 2nd power */ d__1 = x[i - 1]; /* Computing 2nd power */ d__2 = y[i - 1]; t[5] += weight * (d__1 * d__1 + d__2 * d__2); } rsp_(t, &c__3, &c__3, eig, &evec[4]); if (*mass > 0 && first && i_indx(keywrd_1.keywrd, "RC=", 80L, 3L) == 0) { s_wsfe(&io___15); e_wsfe(); for (i = 1; i <= 3; ++i) { if (eig[i - 1] < 3e-4) { eig[i - 1] = 0.; rot[i - 1] = 0.; } else { rot[i - 1] = const2 / eig[i - 1]; } /* L50: */ xyzmom[i - 1] = eig[i - 1] * const1; } s_wsfe(&io___18); for (i = 1; i <= 3; ++i) { do_fio(&c__1, (char *)&rot[i - 1], (ftnlen)sizeof(doublereal)); } e_wsfe(); if (i_indx(keywrd_1.keywrd, "RC=", 80L, 3L) == 0) { s_wsfe(&io___19); e_wsfe(); } s_wsfe(&io___20); for (i = 1; i <= 3; ++i) { do_fio(&c__1, (char *)&xyzmom[i - 1], (ftnlen)sizeof(doublereal)); } e_wsfe(); *c = rot[0]; *b = rot[1]; *a = rot[2]; } /* NOW TO ORIENT THE MOLECULE SO THE CHIRALITY IS PRESERVED */ sum = evec[4] * (evec[8] * evec[12] - evec[9] * evec[11]) + evec[7] * ( evec[11] * evec[6] - evec[5] * evec[12]) + evec[10] * (evec[5] * evec[9] - evec[8] * evec[6]); if (sum < 0.) { for (j = 1; j <= 3; ++j) { /* L60: */ evec[j + 3] = -evec[j + 3]; } } i__1 = *numat; for (i = 1; i <= i__1; ++i) { coord[i * 3 + 1] = x[i - 1]; coord[i * 3 + 2] = y[i - 1]; coord[i * 3 + 3] = z[i - 1]; /* L70: */ } if (*mass > 0) { first = FALSE_; } return 0; } /* axis_ */