/* freqcy.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 numat, nat[86], nfirst[86], nmidle[86], nlast[86], norbs, nelecs, 
	    nalpha, nbeta, nclose, nopen, ndumy;
    doublereal fract;
} molkst_;

#define molkst_1 molkst_

struct {
    doublereal atmass[86];
} atmass_;

#define atmass_1 atmass_

struct {
    doublereal oldf[33411], dummy[33153];
} scrach_;

#define scrach_1 scrach_

/* Table of constant values */

static integer c__1 = 1;

/* Subroutine */ int freqcy_(doublereal *fmatrx, doublereal *freq, doublereal 
	*cnorml, doublereal *redmas, doublereal *travel, logical *eorc)
{
    /* Initialized data */

    static doublereal fact = 6.023e23;

    /* System generated locals */
    integer i__1, i__2, i__3;
    doublereal d__1, d__2;

    /* Builtin functions */
    double sqrt(doublereal), d_sign(doublereal *, doublereal *);

    /* Local variables */
    static integer i, j, k, l;
    extern /* Subroutine */ int frame_(doublereal *, integer *, integer *, 
	    doublereal *);
    static doublereal shift[6];
    static integer j1, n3, ii, ij, jj, linear;
    static doublereal weight, wtmass[258];
    static integer jii;
    extern /* Subroutine */ int rsp_(doublereal *, integer *, integer *, 
	    doublereal *, doublereal *);
    static doublereal sum, c2pi, sum1;

/* 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 */

/* ******************************************************************** */

/*  FRCE CALCULATES THE FORCE CONSTANTS AND VIBRATIONAL FREQUENCIES */
/*       FOR A MOLECULE.  IT USES THE ISOTOPIC MASSES TO WEIGHT THE */
/*       FORCE MATRIX */

/* ON INPUT   FMATRX   =  FORCE MATRIX, OF SIZE NUMAT*3*(NUMAT*3+1)/2. */

/* ******************************************************************** */
    /* Parameter adjustments */
    --travel;
    --redmas;
    --cnorml;
    --freq;
    --fmatrx;

    /* Function Body */

/*    CONVERSION FACTOR FOR SPEED OF LIGHT AND 2 PI. */

    c2pi = 5.3087039056530845e-12;
/* NOW TO CALCULATE THE VIBRATIONAL FREQUENCIES */

/*   FIND CONVERSION CONSTANTS FOR MASS WEIGHTED SYSTEM */
    n3 = molkst_1.numat * 3;
    l = 0;
    i__1 = molkst_1.numat;
    for (i = 1; i <= i__1; ++i) {
	weight = 1.4142136 / sqrt(atmass_1.atmass[i - 1]);
	for (j = 1; j <= 3; ++j) {
	    ++l;
/* L10: */
	    wtmass[l - 1] = weight;
	}
    }
/*    CONVERT TO MASS WEIGHTED FMATRX */
    linear = 0;
    i__1 = n3;
    for (i = 1; i <= i__1; ++i) {
	i__2 = i;
	for (j = 1; j <= i__2; ++j) {
	    ++linear;
	    scrach_1.oldf[linear - 1] = fmatrx[linear] * 1e5;
/* L20: */
	    fmatrx[linear] = fmatrx[linear] * wtmass[i - 1] * wtmass[j - 1];
	}
    }

/*    1.D5 IS TO CONVERT FROM MILLIDYNES/ANGSTROM TO DYNES/CM. */

/*    DIAGONALIZE */
    frame_(&fmatrx[1], &molkst_1.numat, &c__1, shift);
    rsp_(&fmatrx[1], &n3, &n3, &freq[1], &cnorml[1]);
    i__2 = n3;
    for (i = 1; i <= i__2; ++i) {
	j = (integer) ((freq[i] + 50.) * .01);
/* L30: */
	freq[i] -= j * 100;
    }
    i__2 = n3;
    for (i = 1; i <= i__2; ++i) {
/* L40: */
	freq[i] *= 1e5;
    }

/*    CALCULATE REDUCED MASSES, STORE IN REDMAS */

    i__2 = n3;
    for (i = 1; i <= i__2; ++i) {
	ii = (i - 1) * n3;
	sum = 0.;
	i__1 = n3;
	for (j = 1; j <= i__1; ++j) {
	    jii = j + ii;
	    jj = j * (j - 1) / 2;
	    i__3 = j;
	    for (k = 1; k <= i__3; ++k) {
/* L50: */
		sum += cnorml[jii] * scrach_1.oldf[jj + k - 1] * cnorml[k + 
			ii];
	    }
	    i__3 = n3;
	    for (k = j + 1; k <= i__3; ++k) {
/* L60: */
		sum += cnorml[jii] * scrach_1.oldf[k * (k - 1) / 2 + j - 1] * 
			cnorml[k + ii];
	    }
/* L70: */
	}
	sum1 = sum * 2.;
	if ((d__1 = freq[i], abs(d__1)) > abs(sum) * 1e-20) {
	    sum = sum * 1. / freq[i];
	} else {
	    sum = 0.;
	}
	d__2 = sqrt(fact * (d__1 = freq[i], abs(d__1))) * c2pi;
	freq[i] = d_sign(&d__2, &freq[i]);
	if ((d__1 = freq[i], abs(d__1)) < abs(sum1) * 1e20) {
	    sum1 = sqrt((d__1 = freq[i] / (sum1 * 1e-5), abs(d__1)));
	} else {
	    sum1 = 0.;
	}
	if (sum < 0. || sum > 100.) {
	    sum = 0.;
	}

/* 0.0063024=SQRT(2*A*B*C/N) WHERE */
/*         A=1.196D8 = CONVERSION OF CM**(-1) TO (ERGS = DYNE.ANGSTROM
S) */
/*         B=1000.0  = MILLIDYNES TO DYNES */
/*         C=1.D8    = CENTIMETERS TO ANGSTROMS */
/*         N=6.02205D23 = AVOGADRO'S NUMBER */
	travel[i] = sum1 * .0063024;
	if (travel[i] > 1.) {
	    travel[i] = 0.;
	}
/* #      WRITE(6,*)TRAVEL(I) */
/* L80: */
	redmas[i] = sum;
    }
    if (*eorc) {

/*    SWITCH EIGENVALUES TO FREQUENCIES */

/*    CONVERT NORMAL VECTORS TO CARTESIAN COORDINATES */
/*    AND NORMALIZE SO THAT THE TOTAL MOVEMENT IS 1.0 ANGSTROM. */

	ij = 0;
	i__2 = n3;
	for (i = 1; i <= i__2; ++i) {
	    sum = 0.;
	    j = 0;
	    i__1 = molkst_1.numat;
	    for (jj = 1; jj <= i__1; ++jj) {
		sum1 = 0.;
		for (j1 = 1; j1 <= 3; ++j1) {
		    ++j;
		    ++ij;
		    cnorml[ij] *= wtmass[j - 1];
/* L90: */
/* Computing 2nd power */
		    d__1 = cnorml[ij];
		    sum1 += d__1 * d__1;
		}
/* L100: */
		sum += sqrt(sum1);
	    }
	    sum = 1. / sum;
	    ij -= n3;
	    i__1 = n3;
	    for (j = 1; j <= i__1; ++j) {
		++ij;
/* L110: */
		cnorml[ij] *= sum;
	    }
/* L120: */
	}

/*          RETURN HESSIAN IN MILLIDYNES/ANGSTROM IN FMATRX */

	i__2 = linear;
	for (i = 1; i <= i__2; ++i) {
/* L130: */
	    fmatrx[i] = scrach_1.oldf[i - 1] * 1e-5;
	}
    } else {

/*  RETURN HESSIAN AS MASS-WEIGHTED FMATRIX */
	linear = 0;

	i__2 = n3;
	for (i = 1; i <= i__2; ++i) {
	    i__1 = i;
	    for (j = 1; j <= i__1; ++j) {
		++linear;
/* L140: */
		fmatrx[linear] = scrach_1.oldf[linear - 1] * 1e-5 * wtmass[i 
			- 1] * wtmass[j - 1];
	    }
	}
    }
    return 0;
} /* freqcy_ */