/* hcore.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 uspd[215], dumy[215];
} molorb_;

#define molorb_1 molorb_

struct {
    char keywrd[80];
} keywrd_;

#define keywrd_1 keywrd_

struct {
    doublereal tvec[9]	/* was [3][3] */;
    integer id;
} euler_;

#define euler_1 euler_

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

#define multip_1 multip_

struct {
    doublereal core[107];
} core_;

#define core_1 core_

struct {
    doublereal efield[3];
} field_;

#define field_1 field_

/* Table of constant values */

static integer c__1 = 1;

/* Subroutine */ int hcore_(doublereal *coord, doublereal *h, doublereal *w, 
	real *wj, real *wk, doublereal *enuclr)
{
    /* Initialized data */

    static logical first = TRUE_;

    /* Format strings */
    static char fmt_120[] = "(10f8.4)";

    /* System generated locals */
    integer i__1, i__2, i__3, i__4;

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

    /* Local variables */
    static doublereal half;
    static integer ione;
    static doublereal fnuc, enuc;
    static integer i, j;
    static logical debug, fldon;
    static integer i1, i2, j1, j2;
    extern /* Subroutine */ int h1elec_(integer *, integer *, doublereal *, 
	    doublereal *, doublereal *);
    static integer ia, ib, ic;
    static doublereal di[81]	/* was [9][9] */;
    static integer ja, jb, jc, ii, jj, ni, nj, kr;
    static doublereal fldcon, hterme, cutoff, e1b[10], e2a[10];
    extern /* Subroutine */ int rotate_(integer *, integer *, doublereal *, 
	    doublereal *, doublereal *, integer *, doublereal *, doublereal *,
	     doublereal *, doublereal *), vecprt_(doublereal *, integer *);
    static integer im1, io1, jo1;
    extern /* Subroutine */ int solrot_(integer *, integer *, doublereal *, 
	    doublereal *, doublereal *, doublereal *, integer *, doublereal *,
	     doublereal *, doublereal *, doublereal *);
    static doublereal wjd[100], wkd[100];
    static integer kro;

    /* Fortran I/O blocks */
    static cilist io___37 = { 0, 6, 0, "(//10X,'ONE-ELECTRON MATRIX FROM HCO"
	    "RE')", 0 };
    static cilist io___38 = { 0, 6, 0, "(//10X,'TWO-ELECTRON MATRIX IN HCORE"
	    "'/)", 0 };
    static cilist io___39 = { 0, 6, 0, fmt_120, 0 };
    static cilist io___40 = { 0, 6, 0, "(//10X,'TWO-ELECTRON J MATRIX IN HCO"
	    "RE'/)", 0 };
    static cilist io___41 = { 0, 6, 0, fmt_120, 0 };
    static cilist io___42 = { 0, 6, 0, "(//10X,'TWO-ELECTRON K MATRIX IN HCO"
	    "RE'/)", 0 };
    static cilist io___43 = { 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 */

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

/*   HCORE GENERATES THE ONE-ELECTRON MATRIX AND TWO ELECTRON INTEGRALS */
/*         FOR A GIVEN MOLECULE WHOSE GEOMETRY IS DEFINED IN CARTESIAN */
/*         COORDINATES. */

/*  ON INPUT  COORD   = COORDINATES OF THE MOLECULE. */

/*  ON OUTPUT  H      = ONE-ELECTRON MATRIX. */
/*             W      = TWO-ELECTRON INTEGRALS. */
/*             ENUCLR = NUCLEAR ENERGY */
/* ***********************************************************************
 */
    /* Parameter adjustments */
    --wk;
    --wj;
    --w;
    --h;
    coord -= 4;

    /* Function Body */
    if (first) {
	ione = 1;
	cutoff = 1e10;
	if (euler_1.id != 0) {
	    cutoff = 60.;
	}
	if (euler_1.id != 0) {
	    ione = 0;
	}
	first = FALSE_;
	debug = i_indx(keywrd_1.keywrd, "HCORE", 80L, 5L) != 0;
    }
    fldon = FALSE_;
    if (field_1.efield[0] != 0. || field_1.efield[1] != 0. || field_1.efield[
	    2] != 0.) {
	fldcon = 51.4257;
	fldon = TRUE_;
    }
    i__1 = molkst_1.norbs * (molkst_1.norbs + 1) / 2;
    for (i = 1; i <= i__1; ++i) {
/* L10: */
	h[i] = 0.;
    }
    *enuclr = 0.;
    kr = 1;
    i__1 = molkst_1.numat;
    for (i = 1; i <= i__1; ++i) {
	ia = molkst_1.nfirst[i - 1];
	ib = molkst_1.nlast[i - 1];
	ic = molkst_1.nmidle[i - 1];
	ni = molkst_1.nat[i - 1];

/* FIRST WE FILL THE DIAGONALS, AND OFF-DIAGONALS ON THE SAME ATOM */

	i__2 = ib;
	for (i1 = ia; i1 <= i__2; ++i1) {
	    i2 = i1 * (i1 - 1) / 2 + ia - 1;
	    i__3 = i1;
	    for (j1 = ia; j1 <= i__3; ++j1) {
		++i2;
		h[i2] = 0.;
		if (fldon) {
		    io1 = i1 - ia;
		    jo1 = j1 - ia;
		    if (jo1 == 0 && io1 == 1) {
			hterme = multip_1.dd[ni - 1] * -.529177 * 
				field_1.efield[0] * fldcon;
			h[i2] = hterme;
		    }
		    if (jo1 == 0 && io1 == 2) {
			hterme = multip_1.dd[ni - 1] * -.529177 * 
				field_1.efield[1] * fldcon;
			h[i2] = hterme;
		    }
		    if (jo1 == 0 && io1 == 3) {
			hterme = multip_1.dd[ni - 1] * -.529177 * 
				field_1.efield[2] * fldcon;
			h[i2] = hterme;
		    }
		}
/* L20: */
	    }
	    h[i2] = molorb_1.uspd[i1 - 1];
	    if (fldon) {
		fnuc = -(field_1.efield[0] * coord[i * 3 + 1] + 
			field_1.efield[1] * coord[i * 3 + 2] + field_1.efield[
			2] * coord[i * 3 + 3]) * fldcon;
		h[i2] += fnuc;
	    }
/* L30: */
	}

/*   FILL THE ATOM-OTHER ATOM ONE-ELECTRON MATRIX<PSI(LAMBDA)|PSI(SIGM
A)> */

	im1 = i - ione;
	i__2 = im1;
	for (j = 1; j <= i__2; ++j) {
	    half = 1.;
	    if (i == j) {
		half = .5;
	    }
	    ja = molkst_1.nfirst[j - 1];
	    jb = molkst_1.nlast[j - 1];
	    jc = molkst_1.nmidle[j - 1];
	    nj = molkst_1.nat[j - 1];
	    h1elec_(&ni, &nj, &coord[i * 3 + 1], &coord[j * 3 + 1], di);
	    i2 = 0;
	    i__3 = ib;
	    for (i1 = ia; i1 <= i__3; ++i1) {
		ii = i1 * (i1 - 1) / 2 + ja - 1;
		++i2;
		j2 = 0;
		jj = min(i1,jb);
		i__4 = jj;
		for (j1 = ja; j1 <= i__4; ++j1) {
		    ++ii;
		    ++j2;
/* L40: */
		    h[ii] += di[i2 + j2 * 9 - 10];
		}
	    }

/*   CALCULATE THE TWO-ELECTRON INTEGRALS, W; THE ELECTRON NUCLEAR
 TERMS */
/*   E1B AND E2A; AND THE NUCLEAR-NUCLEAR TERM ENUC. */

	    if (euler_1.id == 0) {
		rotate_(&ni, &nj, &coord[i * 3 + 1], &coord[j * 3 + 1], &w[kr]
			, &kr, e1b, e2a, &enuc, &cutoff);
	    } else {
		kro = kr;
		solrot_(&ni, &nj, &coord[i * 3 + 1], &coord[j * 3 + 1], wjd, 
			wkd, &kr, e1b, e2a, &enuc, &cutoff);
		jj = 0;
		i__4 = kr - 1;
		for (ii = kro; ii <= i__4; ++ii) {
		    ++jj;
		    wj[ii] = wjd[jj - 1];
/* L50: */
		    wk[ii] = wkd[jj - 1];
		}
	    }
	    *enuclr += enuc;

/*   ADD ON THE ELECTRON-NUCLEAR ATTRACTION TERM FOR ATOM I. */

	    i2 = 0;
	    i__4 = ic;
	    for (i1 = ia; i1 <= i__4; ++i1) {
		ii = i1 * (i1 - 1) / 2 + ia - 1;
		i__3 = i1;
		for (j1 = ia; j1 <= i__3; ++j1) {
		    ++ii;
		    ++i2;
/* L60: */
		    h[ii] += e1b[i2 - 1] * half;
		}
	    }
	    i__3 = ib;
	    for (i1 = ic + 1; i1 <= i__3; ++i1) {
		ii = i1 * (i1 + 1) / 2;
/* L70: */
		h[ii] += e1b[0] * half;
	    }

/*   ADD ON THE ELECTRON-NUCLEAR ATTRACTION TERM FOR ATOM J. */

	    i2 = 0;
	    i__3 = jc;
	    for (i1 = ja; i1 <= i__3; ++i1) {
		ii = i1 * (i1 - 1) / 2 + ja - 1;
		i__4 = i1;
		for (j1 = ja; j1 <= i__4; ++j1) {
		    ++ii;
		    ++i2;
/* L80: */
		    h[ii] += e2a[i2 - 1] * half;
		}
	    }
	    i__4 = jb;
	    for (i1 = jc + 1; i1 <= i__4; ++i1) {
		ii = i1 * (i1 + 1) / 2;
/* L90: */
		h[ii] += e2a[0] * half;
	    }
/* L100: */
	}
/* L110: */
    }
    if (! debug) {
	return 0;
    }
    s_wsfe(&io___37);
    e_wsfe();
    vecprt_(&h[1], &molkst_1.norbs);
    j = min(400,kr);
    if (euler_1.id == 0) {
	s_wsfe(&io___38);
	e_wsfe();
	s_wsfe(&io___39);
	i__1 = j;
	for (i = 1; i <= i__1; ++i) {
	    do_fio(&c__1, (char *)&w[i], (ftnlen)sizeof(doublereal));
	}
	e_wsfe();
    } else {
	s_wsfe(&io___40);
	e_wsfe();
	s_wsfe(&io___41);
	i__1 = j;
	for (i = 1; i <= i__1; ++i) {
	    do_fio(&c__1, (char *)&wj[i], (ftnlen)sizeof(real));
	}
	e_wsfe();
	s_wsfe(&io___42);
	e_wsfe();
	s_wsfe(&io___43);
	i__1 = j;
	for (i = 1; i <= i__1; ++i) {
	    do_fio(&c__1, (char *)&wk[i], (ftnlen)sizeof(real));
	}
	e_wsfe();
    }
    return 0;
} /* hcore_ */