/* mullik.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], norbx, nelecs, 
	    nalpha, nbeta, nclose, nopen, ndumy;
    doublereal fract;
} molkst_;

#define molkst_1 molkst_

struct {
    char keywrd[80];
} keywrd_;

#define keywrd_1 keywrd_

struct {
    doublereal betas[107], betap[107], betad[107];
} betas_;

#define betas_1 betas_

struct {
    doublereal geo[258]	/* was [3][86] */;
} geom_;

#define geom_1 geom_

struct {
    doublereal zs[107], zp[107], zd[107];
} expont_;

#define expont_1 expont_

/* Table of constant values */

static integer c__5 = 5;
static integer c__1 = 1;
static integer c__2 = 2;

/* Subroutine */ int mullik_(doublereal *c, doublereal *cbeta, logical *uhf, 
	doublereal *h, doublereal *f, integer *norbs, doublereal *vecs, 
	doublereal *store)
{
    /* System generated locals */
    integer i__1, i__2, i__3, i__4;
    doublereal d__1;
    olist o__1;
    cllist cl__1;

    /* Builtin functions */
    integer i_indx(char *, char *, ftnlen, ftnlen);
    double sqrt(doublereal);
    integer f_open(olist *), s_wsue(cilist *), do_uio(integer *, char *, 
	    ftnlen), e_wsue(void), f_clos(cllist *);

    /* Local variables */
    static doublereal eigs[215];
    extern /* Subroutine */ int mult_(doublereal *, doublereal *, doublereal *
	    , integer *);
    static doublereal summ;
    static integer i, j, k, ifact[215];
    static logical graph;
    static doublereal bi, bj;
    static integer if__, jf, ii, ij, jj, il, jl;
    extern /* Subroutine */ int getarg_(integer *, char *, ftnlen);
    static integer linear;
    extern /* Subroutine */ int densit_(doublereal *, integer *, integer *, 
	    integer *, integer *, doublereal *, doublereal *, integer *), 
	    vecprt_(doublereal *, integer *);
    static integer im1;
    extern /* Subroutine */ int gmetry_(doublereal *, doublereal *), rsp_(
	    doublereal *, integer *, integer *, doublereal *, doublereal *);
    static doublereal sum, xyz[258]	/* was [3][86] */;

    /* Fortran I/O blocks */
    static cilist io___19 = { 0, 13, 0, 0, 0 };
    static cilist io___20 = { 0, 13, 0, 0, 0 };
    static cilist io___21 = { 0, 13, 0, 0, 0 };
    static cilist io___23 = { 0, 13, 0, 0, 0 };
    static cilist io___24 = { 0, 13, 0, 0, 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 */

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

/*  MULLIK DOES A MULLIKEN POPULATION ANALYSIS */
/* ON INPUT     C      =  SQUARE ARRAY OF EIGENVECTORS. */
/*              H      =  PACKED ARRAY OF ONE-ELECTRON MATRIX */
/*              F      =  WORKSTORE OF SIZE AT LEAST NORBS*NORBS */
/*              VECS   =  WORKSTORE OF SIZE AT LEAST NORBS*NORBS */
/*              STORE  =  WORKSTORE OF SIZE AT LEAST (NORBS*(NORBS+1))/2 
*/

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

/*  FIRST, RE-CALCULATE THE OVERLAP MATRIX */

/* ********************************************************************* 
*/
    /* Parameter adjustments */
    --store;
    --vecs;
    --f;
    --h;
    --cbeta;
    --c;

    /* Function Body */
    graph = i_indx(keywrd_1.keywrd, "GRAPH", 80L, 5L) != 0;
    i__1 = *norbs;
    for (i = 1; i <= i__1; ++i) {
/* L10: */
	ifact[i - 1] = i * (i - 1) / 2;
    }
    ifact[*norbs] = *norbs * (*norbs + 1) / 2;
    i__1 = molkst_1.numat;
    for (i = 1; i <= i__1; ++i) {
	if__ = molkst_1.nfirst[i - 1];
	il = molkst_1.nlast[i - 1];
	im1 = i - 1;
	bi = betas_1.betas[molkst_1.nat[i - 1] - 1];
	i__2 = il;
	for (k = if__; k <= i__2; ++k) {
	    ii = k * (k - 1) / 2;
	    i__3 = im1;
	    for (j = 1; j <= i__3; ++j) {
		jf = molkst_1.nfirst[j - 1];
		jl = molkst_1.nlast[j - 1];
		bj = betas_1.betas[molkst_1.nat[j - 1] - 1];
		i__4 = jl;
		for (jj = jf; jj <= i__4; ++jj) {
		    ij = ii + jj;
		    h[ij] = h[ij] * 2. / (bi + bj) + 1e-14;
/*  THE  +1.D-14 IS TO PREVENT POSSIBLE ERRORS IN THE DIAG
ONALIZATION. */
		    store[ij] = h[ij];
/* L20: */
		    bj = betas_1.betap[molkst_1.nat[j - 1] - 1];
		}
/* L30: */
	    }
	    i__3 = k;
	    for (jj = if__; jj <= i__3; ++jj) {
		ij = ii + jj;
		store[ij] = 0.;
/* L40: */
		h[ij] = 0.;
	    }
/* L50: */
	    bi = betas_1.betap[molkst_1.nat[i - 1] - 1];
	}
    }
    i__2 = *norbs;
    for (i = 1; i <= i__2; ++i) {
	store[ifact[i]] = 1.;
/* L60: */
	h[ifact[i]] = 1.;
    }
    rsp_(&h[1], norbs, norbs, eigs, &vecs[1]);
    i__2 = *norbs;
    for (i = 1; i <= i__2; ++i) {
/* L70: */
	eigs[i - 1] = 1. / sqrt((d__1 = eigs[i - 1], abs(d__1)));
    }
    ij = 0;
    i__2 = *norbs;
    for (i = 1; i <= i__2; ++i) {
	i__1 = i;
	for (j = 1; j <= i__1; ++j) {
	    ++ij;
	    sum = 0.;
	    i__3 = *norbs;
	    for (k = 1; k <= i__3; ++k) {
/* L80: */
		sum += vecs[i + (k - 1) * *norbs] * eigs[k - 1] * vecs[j + (k 
			- 1) * *norbs];
	    }
	    f[i + (j - 1) * *norbs] = sum;
/* L90: */
	    f[j + (i - 1) * *norbs] = sum;
	}
    }
    if (graph) {
	gmetry_(geom_1.geo, xyz);

/* WRITE TO DISK THE FOLLOWING DATA FOR GRAPHICS CALCULATION, IN ORDER
: */

/*      NUMBER OF ATOMS, ORBITAL, ELECTRONS */
/*      ALL ATOMIC COORDINATES */
/*      ORBITAL COUNTERS */
/*      ORBITAL EXPONENTS, S, P, AND D, AND ATOMIC NUMBERS */
/*      EIGENVECTORS (M.O.S NOT RE-NORMALIZED) */
/*      INVERSE-SQUARE ROOT OF THE OVERLAP MATRIX. */

	getarg_(&c__5, argz_1.argz, 512L);
	o__1.oerr = 0;
	o__1.ounit = 13;
	o__1.ofnmlen = 512;
	o__1.ofnm = argz_1.argz;
	o__1.orl = 0;
	o__1.osta = 0;
	o__1.oacc = 0;
	o__1.ofm = "UNFORMATTED";
	o__1.oblnk = 0;
	f_open(&o__1);
	s_wsue(&io___19);
	do_uio(&c__1, (char *)&molkst_1.numat, (ftnlen)sizeof(integer));
	do_uio(&c__1, (char *)&(*norbs), (ftnlen)sizeof(integer));
	do_uio(&c__1, (char *)&molkst_1.nelecs, (ftnlen)sizeof(integer));
	for (i = 1; i <= 3; ++i) {
	    i__1 = molkst_1.numat;
	    for (j = 1; j <= i__1; ++j) {
		do_uio(&c__1, (char *)&xyz[i + j * 3 - 4], (ftnlen)sizeof(
			doublereal));
	    }
	}
	e_wsue();
	s_wsue(&io___20);
	i__1 = molkst_1.numat;
	for (i = 1; i <= i__1; ++i) {
	    do_uio(&c__1, (char *)&molkst_1.nlast[i - 1], (ftnlen)sizeof(
		    integer));
	    do_uio(&c__1, (char *)&molkst_1.nfirst[i - 1], (ftnlen)sizeof(
		    integer));
	}
	e_wsue();
	s_wsue(&io___21);
	i__1 = molkst_1.numat;
	for (i = 1; i <= i__1; ++i) {
	    do_uio(&c__1, (char *)&expont_1.zs[molkst_1.nat[i - 1] - 1], (
		    ftnlen)sizeof(doublereal));
	}
	i__2 = molkst_1.numat;
	for (i = 1; i <= i__2; ++i) {
	    do_uio(&c__1, (char *)&expont_1.zp[molkst_1.nat[i - 1] - 1], (
		    ftnlen)sizeof(doublereal));
	}
	i__3 = molkst_1.numat;
	for (i = 1; i <= i__3; ++i) {
	    do_uio(&c__1, (char *)&expont_1.zd[molkst_1.nat[i - 1] - 1], (
		    ftnlen)sizeof(doublereal));
	}
	i__4 = molkst_1.numat;
	for (i = 1; i <= i__4; ++i) {
	    do_uio(&c__1, (char *)&molkst_1.nat[i - 1], (ftnlen)sizeof(
		    integer));
	}
	e_wsue();
	linear = *norbs * *norbs;
	s_wsue(&io___23);
	i__1 = linear;
	for (i = 1; i <= i__1; ++i) {
	    do_uio(&c__1, (char *)&c[i], (ftnlen)sizeof(doublereal));
	}
	e_wsue();
	s_wsue(&io___24);
	i__1 = linear;
	for (i = 1; i <= i__1; ++i) {
	    do_uio(&c__1, (char *)&f[i], (ftnlen)sizeof(doublereal));
	}
	e_wsue();
	cl__1.cerr = 0;
	cl__1.cunit = 13;
	cl__1.csta = 0;
	f_clos(&cl__1);
	if (i_indx(keywrd_1.keywrd, "MULLIK", 80L, 6L) == 0) {
	    return 0;
	}
    }

/* OTHERWISE PERFORM MULLIKEN ANALYSIS */

    mult_(&c[1], &f[1], &vecs[1], norbs);
    i = -1;
    densit_(&vecs[1], norbs, norbs, &molkst_1.nclose, &molkst_1.nopen, &
	    molkst_1.fract, &c[1], &c__2);
    linear = *norbs * (*norbs + 1) / 2;
    i__1 = linear;
    for (i = 1; i <= i__1; ++i) {
/* L100: */
	c[i] *= store[i];
    }
    summ = 0.;
    i__1 = *norbs;
    for (i = 1; i <= i__1; ++i) {
	sum = 0.;
	i__2 = i;
	for (j = 1; j <= i__2; ++j) {
/* L110: */
	    sum += c[ifact[i - 1] + j];
	}
	i__2 = *norbs;
	for (j = i + 1; j <= i__2; ++j) {
/* L120: */
	    sum += c[ifact[j - 1] + i];
	}
	summ += sum;
/* L130: */
	c[ifact[i]] = sum;
    }
    vecprt_(&c[1], norbs);
    return 0;
} /* mullik_ */