/* spcg.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 gss[107], gsp[107], gpp[107], gp2[107], hsp[107], gsd[107], 
	    gpd[107], gdd[107];
} twoele_;

#define twoele_1 twoele_

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

#define euler_1 euler_

struct {
    char keywrd[80];
} keywrd_;

#define keywrd_1 keywrd_

doublereal spcg_(doublereal *c1, doublereal *c2, doublereal *c3, doublereal *
	c4, doublereal *w, real *wj)
{
    /* Initialized data */

    static integer itype = 1;

    /* System generated locals */
    integer i__1, i__2, i__3, i__4, i__5, i__6;
    doublereal ret_val;

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

    /* Local variables */
    static doublereal wint, temp1;
    static integer i, j, k, l;
    static doublereal atemp, elrep;
    static integer i1, j1, k1, ia, ib, ja, jb, ii, jj, kk, is, kr, ix, iy, iz,
	     iminus, im1, is1;
    static logical lid;
    static integer izn;

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

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

/*     SPCG CALCULATES THE REPULSION BETWEEN ELECTRON 1 IN MOLECULAR */
/*     ORBITALS C1 AND C2 AND ELECTRON 2 IN M.O.S C3 AND C4 FOR THE */
/*     VALENCE SP SHELL AT AN MNDO OR MINDO/3 LEVEL. */

/*                            USAGE */
/*      XJ=SPCG(C(1,I),C(1,J),C(1,K),C(1,L)) */
/*  OR, XJ=<I(1),J(1)/K(2),L(2)> */

/*    ON INPUT C1    THE FIRST COLUMN MOLECULAR ORBITAL OF ELECTRON ONE. 
*/
/*             C2        SECOND */
/*             C3        FIRST                                      TWO. 
*/
/*             C4        SECOND */

/*   ON OUTPUT SPCG   =   <C1(1)*C2(1)/C3(2)*C4(2)> */
/* ********************************************************************* 
*/
    /* Parameter adjustments */
    --wj;
    --w;
    --c4;
    --c3;
    --c2;
    --c1;

    /* Function Body */
L10:
    switch (itype) {
	case 1:  goto L20;
	case 2:  goto L70;
	case 3:  goto L30;
    }
L20:
    lid = euler_1.id == 0;
    if (i_indx(keywrd_1.keywrd, "MINDO", 80L, 5L) != 0) {
	itype = 2;
    } else {
	itype = 3;
    }
    goto L10;
/*                           ****************** */
/*                           *      MNDO      * */
/*                           *     OPTION     * */
/*                           ****************** */
L30:
    ret_val = 0.;
    kk = 0;
    i__1 = molkst_1.numat;
    for (ii = 1; ii <= i__1; ++ii) {
	ia = molkst_1.nfirst[ii - 1];
	ib = molkst_1.nlast[ii - 1];
	iminus = ii - 1;
	i__2 = iminus;
	for (jj = 1; jj <= i__2; ++jj) {
	    ja = molkst_1.nfirst[jj - 1];
	    jb = molkst_1.nlast[jj - 1];
	    i__3 = ib;
	    for (i = ia; i <= i__3; ++i) {
		i__4 = i;
		for (j = ia; j <= i__4; ++j) {
		    i__5 = jb;
		    for (k = ja; k <= i__5; ++k) {
			i__6 = k;
			for (l = ja; l <= i__6; ++l) {
			    ++kk;
			    if (lid) {
				wint = w[kk];
			    } else {
				wint = wj[kk];
			    }
			    ret_val += wint * (c1[i] * c2[j] * c3[k] * c4[l] 
				    + c1[k] * c2[l] * c3[i] * c4[j]);
			    if (i != j) {
				ret_val += wint * (c1[j] * c2[i] * c3[k] * c4[
					l] + c1[k] * c2[l] * c3[j] * c4[i]);
			    }
			    if (k != l) {
				ret_val += wint * (c1[i] * c2[j] * c3[l] * c4[
					k] + c1[l] * c2[k] * c3[i] * c4[j]);
			    }
			    if (i != j && k != l) {
				ret_val += wint * (c1[j] * c2[i] * c3[l] * c4[
					k] + c1[l] * c2[k] * c3[j] * c4[i]);
			    }
/* L40: */
			}
		    }
		}
	    }
/* L50: */
	}
/* L60: */
    }
    goto L110;
/*                           ****************** */
/*                           *     MINDO/3    * */
/*                           *     OPTION     * */
/*                           ****************** */
L70:
    ret_val = 0.;
    kr = 0;
    i__1 = molkst_1.numat;
    for (ii = 1; ii <= i__1; ++ii) {
	ia = molkst_1.nfirst[ii - 1];
	ib = molkst_1.nlast[ii - 1];
	im1 = ii - 1;
	i__2 = im1;
	for (jj = 1; jj <= i__2; ++jj) {
	    ++kr;
	    if (lid) {
		elrep = w[kr];
	    } else {
		elrep = wj[kr];
	    }
	    ja = molkst_1.nfirst[jj - 1];
	    jb = molkst_1.nlast[jj - 1];
	    i__6 = ib;
	    for (i = ia; i <= i__6; ++i) {
		i__5 = jb;
		for (k = ja; k <= i__5; ++k) {
/* L80: */
		    ret_val += elrep * (c1[i] * c2[i] * c3[k] * c4[k] + c1[k] 
			    * c2[k] * c3[i] * c4[i]);
		}
	    }
/* L90: */
	}
/* L100: */
    }
L110:
    atemp = ret_val;
    is1 = 0;
    i__1 = molkst_1.numat;
    for (i1 = 1; i1 <= i__1; ++i1) {
	++is1;
	izn = molkst_1.nat[i1 - 1];

/*      (SS/SS) */

	ret_val += c1[is1] * c2[is1] * c3[is1] * c4[is1] * twoele_1.gss[izn - 
		1];
	if (izn < 3) {
	    goto L150;
	}
	is = is1;
	++is1;
	ix = is1;
	++is1;
	iy = is1;
	++is1;
	iz = is1;
	ret_val += twoele_1.gpp[izn - 1] * (c1[ix] * c2[ix] * c3[ix] * c4[ix] 
		+ c1[iy] * c2[iy] * c3[iy] * c4[iy] + c1[iz] * c2[iz] * c3[iz]
		 * c4[iz]);
	ret_val += twoele_1.gsp[izn - 1] * (c1[is] * c2[is] * c3[ix] * c4[ix] 
		+ c1[is] * c2[is] * c3[iy] * c4[iy] + c1[is] * c2[is] * c3[iz]
		 * c4[iz] + c1[ix] * c2[ix] * c3[is] * c4[is] + c1[iy] * c2[
		iy] * c3[is] * c4[is] + c1[iz] * c2[iz] * c3[is] * c4[is]);
	ret_val += twoele_1.gp2[izn - 1] * (c1[ix] * c2[ix] * c3[iy] * c4[iy] 
		+ c1[ix] * c2[ix] * c3[iz] * c4[iz] + c1[iy] * c2[iy] * c3[iz]
		 * c4[iz] + c1[iy] * c2[iy] * c3[ix] * c4[ix] + c1[iz] * c2[
		iz] * c3[ix] * c4[ix] + c1[iz] * c2[iz] * c3[iy] * c4[iy]);
	temp1 = twoele_1.hsp[izn - 1];
	i__2 = iz;
	for (j1 = ix; j1 <= i__2; ++j1) {
	    ret_val += temp1 * (c1[is] * c2[j1] * c3[j1] * c4[is] + c1[is] * 
		    c2[j1] * c3[is] * c4[j1] + c1[j1] * c2[is] * c3[is] * c4[
		    j1] + c1[j1] * c2[is] * c3[j1] * c4[is]);
/* L120: */
	}
	temp1 = (twoele_1.gpp[izn - 1] - twoele_1.gp2[izn - 1]) * .5;
	i__2 = iz;
	for (j1 = ix; j1 <= i__2; ++j1) {
	    i__5 = iz;
	    for (k1 = ix; k1 <= i__5; ++k1) {
		if (j1 == k1) {
		    goto L130;
		}
		ret_val += temp1 * (c1[j1] * c2[k1] * c3[j1] * c4[k1] + c1[j1]
			 * c2[k1] * c3[k1] * c4[j1]);
L130:
		;
	    }
/* L140: */
	}
L150:
	;
    }
    return ret_val;
} /* spcg_ */