/* dipole.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 {
    char keywrd[80];
} keywrd_;

#define keywrd_1 keywrd_

struct {
    doublereal ams[107];
} istope_;

#define istope_1 istope_

struct {
    doublereal dd[107], qq[107], am[107], ad[107], aq[107];
} multip_;

#define multip_1 multip_

/* Table of constant values */

static integer c__1 = 1;

doublereal dipole_(doublereal *p, doublereal *q, doublereal *coord, 
	doublereal *dipvec)
{
    /* Initialized data */

    static struct {
	doublereal e_1;
	doublereal fill_2[106];
	doublereal e_3;
	doublereal fill_4[3];
	doublereal e_5[5];
	doublereal fill_6[4];
	doublereal e_7[4];
	doublereal fill_8[90];
	} equiv_17 = { 0., {0}, 0., {0}, 6.520587, 4.253676, 2.947501, 
		2.139793, 2.2210719, {0}, 6.663059, 5.657623, 6.345552, 
		2.522964 };

#define hyf ((doublereal *)&equiv_17)

    static logical first = TRUE_;

    /* Format strings */
    static char fmt_120[] = "(\002 DIPOLE\002,11x,\002X \002,8x,\002Y \002,8"
	    "x,\002Z \002,6x,\002TOTAL\002,/,\002 POINT-CHG.\002,4f10.3/,\002"
	    " HYBRID\002,4x,4f10.3/,\002 SUM\002,7x,4f10.3)";

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

    /* Builtin functions */
    integer i_indx(char *, char *, ftnlen, ftnlen);
    double sqrt(doublereal);
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);

    /* Local variables */
    static integer i, j, k, l;
    static logical force;
    static integer itype;
    static doublereal wtmol;
    static integer ia, ni;
    static logical chargd;
    static doublereal center[3], dip[12]	/* was [4][3] */, sum;

    /* Fortran I/O blocks */
    static cilist io___16 = { 0, 6, 0, fmt_120, 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 */

/* ***********************************************************************
 */
/*     DIPOLE CALCULATES DIPOLE MOMENTS */

/*  ON INPUT P     = DENSITY MATRIX */
/*           Q     = TOTAL ATOMIC CHARGES, (NUCLEAR + ELECTRONIC) */
/*           NUMAT = NUMBER OF ATOMS IN MOLECULE */
/*           NAT   = ATOMIC NUMBERS OF ATOMS */
/*           NFIRST= START OF ATOM ORBITAL COUNTERS */
/*           COORD = COORDINATES OF ATOMS */

/*  OUTPUT  DIPOLE = DIPOLE MOMENT */
/* ***********************************************************************
 */

/*     IN THE ZDO APPROXIMATION, ONLY TWO TERMS ARE RETAINED IN THE */
/*     CALCULATION OF DIPOLE MOMENTS. */
/*     1. THE POINT CHARGE TERM (INDEPENDENT OF PARAMETERIZATION). */
/*     2. THE ONE-CENTER HYBRIDIZATION TERM, WHICH ARISES FROM MATRIX */
/*     ELEMENTS OF THE FORM <NS/R/NP>. THIS TERM IS A FUNCTION OF */
/*     THE SLATER EXPONENTS (ZS,ZP) AND IS THUS DEPENDENT ON PARAMETER- */
/*     IZATION. THE HYBRIDIZATION FACTORS (HYF(I)) USED IN THIS SUB- */
/*     ROUTINE ARE CALCULATED FROM THE FOLLOWING FORMULAE. */
/*     FOR SECOND ROW ELEMENTS <2S/R/2P> */
/*     HYF(I)= 469.56193322*(SQRT(((ZS(I)**5)*(ZP(I)**5)))/ */
/*           ((ZS(I) + ZP(I))**6)) */
/*     FOR THIRD ROW ELEMENTS <3S/R/3P> */
/*     HYF(I)=2629.107682607*(SQRT(((ZS(I)**7)*(ZP(I)**7)))/ */
/*           ((ZS(I) + ZP(I))**8)) */
/*     FOR FOURTH ROW ELEMENTS AND UP : */
/*     HYF(I)=2*(2.10716)*DD(I) */
/*     WHERE DD(I) IS THE CHARGE SEPARATION IN ATOMIC UNITS */


/*     REFERENCES: */
/*     J.A.POPLE & D.L.BEVERIDGE: APPROXIMATE M.O. THEORY */
/*     S.P.MCGLYNN, ET AL: APPLIED QUANTUM CHEMISTRY */

    /* Parameter adjustments */
    --dipvec;
    coord -= 4;
    --q;
    --p;

    /* Function Body */
    if (first) {
	for (i = 2; i <= 107; ++i) {
/* L10: */
	    hyf[i - 1] = multip_1.dd[i - 1] * 5.0832f;
	}
	wtmol = 0.;
	sum = 0.;
	i__1 = molkst_1.numat;
	for (i = 1; i <= i__1; ++i) {
	    wtmol += istope_1.ams[molkst_1.nat[i - 1] - 1];
/* L20: */
	    sum += q[i];
	}
	chargd = abs(sum) > .5;
	first = FALSE_;
	force = i_indx(keywrd_1.keywrd, "FORCE", 80L, 5L) + i_indx(
		keywrd_1.keywrd, "IRC", 80L, 3L) != 0;
	itype = 1;
	if (i_indx(keywrd_1.keywrd, "MINDO", 80L, 5L) != 0) {
	    itype = 2;
	}
    }
    if (chargd) {

/*   NEED TO RESET ION'S POSITION SO THAT THE CENTER OF MASS IS AT THE
 */
/*   ORIGIN. */

	for (i = 1; i <= 3; ++i) {
/* L30: */
	    center[i - 1] = 0.;
	}
	for (i = 1; i <= 3; ++i) {
	    i__1 = molkst_1.numat;
	    for (j = 1; j <= i__1; ++j) {
/* L40: */
		center[i - 1] += istope_1.ams[molkst_1.nat[j - 1] - 1] * 
			coord[i + j * 3];
	    }
	}
	for (i = 1; i <= 3; ++i) {
/* L50: */
	    center[i - 1] /= wtmol;
	}
	for (i = 1; i <= 3; ++i) {
	    i__1 = molkst_1.numat;
	    for (j = 1; j <= i__1; ++j) {
/* L60: */
		coord[i + j * 3] -= center[i - 1];
	    }
	}
    }
    for (i = 1; i <= 4; ++i) {
	for (j = 1; j <= 3; ++j) {
/* L70: */
	    dip[i + (j << 2) - 5] = 0.;
	}
    }
    i__1 = molkst_1.numat;
    for (i = 1; i <= i__1; ++i) {
	ni = molkst_1.nat[i - 1];
	ia = molkst_1.nfirst[i - 1];
	l = molkst_1.nlast[i - 1] - ia;
	i__2 = l;
	for (j = 1; j <= i__2; ++j) {
	    k = (ia + j) * (ia + j - 1) / 2 + ia;
/* L80: */
	    dip[j + 3] -= hyf[ni + itype * 107 - 108] * p[k];
	}
	for (j = 1; j <= 3; ++j) {
/* L90: */
	    dip[j - 1] += q[i] * 4.803 * coord[j + i * 3];
	}
    }
    for (j = 1; j <= 3; ++j) {
/* L100: */
	dip[j + 7] = dip[j + 3] + dip[j - 1];
    }
    for (j = 1; j <= 3; ++j) {
/* L110: */
/* Computing 2nd power */
	d__1 = dip[(j << 2) - 4];
/* Computing 2nd power */
	d__2 = dip[(j << 2) - 3];
/* Computing 2nd power */
	d__3 = dip[(j << 2) - 2];
	dip[(j << 2) - 1] = sqrt(d__1 * d__1 + d__2 * d__2 + d__3 * d__3);
    }
    if (force) {
	dipvec[1] = dip[8];
	dipvec[2] = dip[9];
	dipvec[3] = dip[10];
    } else {
	s_wsfe(&io___16);
	for (j = 1; j <= 3; ++j) {
	    for (i = 1; i <= 4; ++i) {
		do_fio(&c__1, (char *)&dip[i + (j << 2) - 5], (ftnlen)sizeof(
			doublereal));
	    }
	}
	e_wsfe();
    }
    ret_val = dip[11];
    return ret_val;


} /* dipole_ */

#undef hyf