/* dipind.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 {
    doublereal core[107];
} core_;

#define core_1 core_

struct {
    doublereal p[23220], pa[23220], pb[23220];
} densty_;

#define densty_1 densty_

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

#define geom_1 geom_

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_

/* Subroutine */ int dipind_(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_20 = { 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_20)

    static logical first = TRUE_;

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

    /* Builtin functions */
    integer i_indx(char *, char *, ftnlen, ftnlen);
    double sqrt(doublereal);

    /* Local variables */
    static integer i, j, k, l;
    static doublereal q[215];
    extern /* Subroutine */ int chrge_(doublereal *, doublereal *);
    static logical force;
    static doublereal coord[258]	/* was [3][86] */;
    static integer itype;
    static doublereal wtmol, q2[215];
    static integer ia, ni;
    static logical chargd;
    static doublereal dipcom, center[3];
    extern /* Subroutine */ int gmetry_(doublereal *, doublereal *);
    static doublereal dip[12]	/* was [4][3] */, sum;

/* ............................................................... */
/*  MODIFICATION OF DIPOLE SUBROUTINE FOR USE IN THE CALUCLATION */
/*  OF THE INDUCED DIPOLES FOR POLARIZABILITIES. */
/* ............................................................... */
/* 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;

    /* Function Body */

/*  SETUP FOR DIPOLE CALCULATION */

    chrge_(densty_1.p, q2);
    i__1 = molkst_1.numat;
    for (i = 1; i <= i__1; ++i) {
	l = molkst_1.nat[i - 1];
	q[i - 1] = core_1.core[l - 1] - q2[i - 1];
/* L10: */
    }
    gmetry_(geom_1.geo, coord);

    if (first) {
	for (i = 2; i <= 107; ++i) {
/* L20: */
	    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];
/* L30: */
	    sum += q[i - 1];
	}
	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) {
/* L40: */
	    center[i - 1] = 0.;
	}
	for (i = 1; i <= 3; ++i) {
	    i__1 = molkst_1.numat;
	    for (j = 1; j <= i__1; ++j) {
/* L50: */
		center[i - 1] += istope_1.ams[molkst_1.nat[j - 1] - 1] * 
			coord[i + j * 3 - 4];
	    }
	}
	for (i = 1; i <= 3; ++i) {
/* L60: */
	    center[i - 1] /= wtmol;
	}
	for (i = 1; i <= 3; ++i) {
	    i__1 = molkst_1.numat;
	    for (j = 1; j <= i__1; ++j) {
/* L70: */
		coord[i + j * 3 - 4] -= center[i - 1];
	    }
	}
    }
    for (i = 1; i <= 4; ++i) {
	for (j = 1; j <= 3; ++j) {
/* L80: */
	    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;
/* L90: */
	    dip[j + 3] -= hyf[ni + itype * 107 - 108] * densty_1.p[k - 1];
	}
	for (j = 1; j <= 3; ++j) {
/* L100: */
	    dip[j - 1] += q[i - 1] * 4.803 * coord[j + i * 3 - 4];
	}
    }
    for (j = 1; j <= 3; ++j) {
/* L110: */
	dip[j + 7] = dip[j + 3] + dip[j - 1];
    }
    for (j = 1; j <= 3; ++j) {
/* L120: */
/* 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);
    }
    dipvec[1] = -dip[8];
    dipvec[2] = -dip[9];
    dipvec[3] = -dip[10];
    dipcom = dip[11];
/*      WRITE (6,60) ((DIP(I,J),I=1,4),J=1,3) */
/*   60 FORMAT (3(4F10.3)) */
    return 0;

} /* dipind_ */

#undef hyf