/* fock2d.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 {
    doublereal tvec[9]	/* was [3][3] */;
    integer id;
} euler_;

#define euler_1 euler_

struct {
    char keywrd[80];
} keywrd_;

#define keywrd_1 keywrd_

/* Subroutine */ int fock2d_(doublereal *f, doublereal *ptot, doublereal *p, 
	doublereal *w, doublereal *wj, doublereal *wk, integer *numat, 
	integer *nfirst, integer *nmidle, integer *nlast)
{
    /* Initialized data */

    static integer itype = 1;

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

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

    /* Local variables */
    static integer ione;
    static doublereal drep, dpop[2], a;
    static integer i, j, k, l, ifact[20];
    static doublereal elexc, elrep;
    static integer norbs, i2, j2;
    static doublereal sppop[2];
    static integer i1fact[20], linea1;
    static doublereal aa, bb;
    static integer ia, ib, ic, ja, jb, jc, ka, ii, kb, jj, ij, kk, kc, ik, jk,
	     il, jl, kl;
    static doublereal aj, ak;
    static integer kr, ll, iminus, im1;
    static logical lid;
    static doublereal sum;

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

/* FOCK2 FORMS THE TWO-ELECTRON TWO-CENTER REPULSION PART OF THE FOCK */
/* MATRIX */
/* ON INPUT  PTOT = TOTAL DENSITY MATRIX. */
/*           P    = ALPHA OR BETA DENSITY MATRIX. */
/*           W    = TWO-ELECTRON INTEGRAL MATRIX. */

/*  ON OUTPUT F   = PARTIAL FOCK MATRIX */
/* ***********************************************************************
 */
    /* Parameter adjustments */
    --nlast;
    --nmidle;
    --nfirst;
    --wk;
    --wj;
    --w;
    --p;
    --ptot;
    --f;

    /* Function Body */
L10:
    switch (itype) {
	case 1:  goto L20;
	case 2:  goto L170;
	case 3:  goto L40;
    }
L20:

/*   SET UP ARRAY OF (I*(I-1))/2 */

    for (i = 1; i <= 20; ++i) {
	ifact[i - 1] = i * (i - 1) / 2;
/* L30: */
	i1fact[i - 1] = ifact[i - 1] + i;
    }
    lid = euler_1.id == 0;
    ione = 1;
    if (i_indx(keywrd_1.keywrd, "MINDO", 80L, 5L) != 0) {
	itype = 2;
    } else {
	itype = 3;
    }
    goto L10;
L40:
    kk = 0;
    norbs = nlast[*numat];
    linea1 = norbs * (norbs + 1) / 2 + 1;
    p[linea1] = 0.;
    i__1 = *numat;
    for (ii = 1; ii <= i__1; ++ii) {
	ia = nfirst[ii];
	ib = nlast[ii];
	ic = nmidle[ii];
	sum = 0.;
	i__2 = ic;
	for (i = ia; i <= i__2; ++i) {
/* L50: */
	    sum += ptot[i1fact[i - 1]];
	}
	sppop[ii - 1] = sum;
	sum = 0.;
	i__2 = ib;
	for (i = ic + 1; i <= i__2; ++i) {
/* L60: */
	    sum += ptot[i1fact[i - 1]];
	}
	dpop[ii - 1] = sum;
	iminus = ii - ione;
	i__2 = iminus;
	for (jj = 1; jj <= i__2; ++jj) {
	    ja = nfirst[jj];
	    jb = nlast[jj];
	    jc = nmidle[jj];
	    drep = wj[kk + 1];
	    if (lid) {
		drep = w[kk + 1];
		i__3 = ic;
		for (i = ia; i <= i__3; ++i) {
		    ka = ifact[i - 1];
		    i__4 = i;
		    for (j = ia; j <= i__4; ++j) {
			kb = ifact[j - 1];
			ij = ka + j;
			aa = 2.;
			if (i == j) {
			    aa = 1.;
			}
			i__5 = jc;
			for (k = ja; k <= i__5; ++k) {
			    kc = ifact[k - 1];
			    ik = ka + k;
			    jk = kb + k;
			    i__6 = k;
			    for (l = ja; l <= i__6; ++l) {
				il = ka + l;
				jl = kb + l;
				kl = kc + l;
				bb = 2.;
				if (k == l) {
				    bb = 1.;
				}
				++kk;
				a = w[kk];

/*     A  IS THE REPULSION INTEGRAL (I,J/K,L) 
WHERE ORBITALS I AND J ARE */
/*     ON ATOM II, AND ORBITALS K AND L ARE ON
 ATOM JJ. */
/*     AA AND BB ARE CORRECTION FACTORS SINCE 
*/
/*     (I,J/K,L)=(J,I/K,L)=(I,J/L,K)=(J,I/L,K)
 */
/*     IJ IS THE LOCATION OF THE MATRIX ELEMEN
TS BETWEEN ATOMIC ORBITALS */
/*     I AND J.  SIMILARLY FOR IK ETC. */

/* THIS FORMS THE TWO-ELECTRON TWO-CENTER REPU
LSION PART OF THE FOCK */
/* MATRIX.  THE CODE HERE IS HARD TO FOLLOW, A
ND IMPOSSIBLE TO MODIFY!, */
/* BUT IT WORKS, */
				f[ij] += bb * a * ptot[kl];
				f[kl] += aa * a * ptot[ij];
				a = a * aa * bb * .25;
				f[ik] -= a * p[jl];
				f[il] -= a * p[jk];
				f[jk] -= a * p[il];
				f[jl] -= a * p[ik];
/* L70: */
			    }
			}
		    }
		}
	    } else {
		i__6 = ic;
		for (i = ia; i <= i__6; ++i) {
		    ka = ifact[i - 1];
		    i__5 = i;
		    for (j = ia; j <= i__5; ++j) {
			kb = ifact[j - 1];
			ij = ka + j;
			aa = 2.;
			if (i == j) {
			    aa = 1.;
			}
			i__4 = jc;
			for (k = ja; k <= i__4; ++k) {
			    kc = ifact[k - 1];
			    if (i >= k) {
				ik = ka + k;
			    } else {
				ik = linea1;
			    }
			    if (j >= k) {
				jk = kb + k;
			    } else {
				jk = linea1;
			    }
			    i__3 = k;
			    for (l = ja; l <= i__3; ++l) {
				if (i >= l) {
				    il = ka + l;
				} else {
				    il = linea1;
				}
				if (j >= l) {
				    jl = kb + l;
				} else {
				    jl = linea1;
				}
				kl = kc + l;
				bb = 2.;
				if (k == l) {
				    bb = 1.;
				}
				++kk;
				aj = wj[kk];
				ak = wk[kk];

/*     A  IS THE REPULSION INTEGRAL (I,J/K,L) 
WHERE ORBITALS I AND J ARE */
/*     ON ATOM II, AND ORBITALS K AND L ARE ON
 ATOM JJ. */
/*     AA AND BB ARE CORRECTION FACTORS SINCE 
*/
/*     (I,J/K,L)=(J,I/K,L)=(I,J/L,K)=(J,I/L,K)
 */
/*     IJ IS THE LOCATION OF THE MATRIX ELEMEN
TS BETWEEN ATOMIC ORBITALS */
/*     I AND J.  SIMILARLY FOR IK ETC. */

/* THIS FORMS THE TWO-ELECTRON TWO-CENTER REPU
LSION PART OF THE FOCK */
/* MATRIX.  THE CODE HERE IS HARD TO FOLLOW, A
ND IMPOSSIBLE TO MODIFY!, */
/* BUT IT WORKS, */
				if (kl <= ij) {
				    if (i == k && aa + bb < 2.1) {
					bb *= .5;
					aa *= .5;
					f[ij] += bb * aj * ptot[kl];
					f[kl] += aa * aj * ptot[ij];
				    } else {
					f[ij] += bb * aj * ptot[kl];
					f[kl] += aa * aj * ptot[ij];
					a = ak * aa * bb * .25;
					f[ik] -= a * p[jl];
					f[il] -= a * p[jk];
					f[jk] -= a * p[il];
					f[jl] -= a * p[ik];
				    }
				}
/* L80: */
			    }
			}
		    }
		}
	    }

/*   D-ORBITAL CORRECTION */

	    i__3 = ib;
	    for (i = ic + 1; i <= i__3; ++i) {
		ka = ifact[i - 1];
		i__4 = jb;
		for (j = ja; j <= i__4; ++j) {
		    ij = ka + j;

/*   ATOM J (S, P, AND D (IF PRESENT)) EXCHANGE WITH ATOM 
I (D ONLY) */

/* L90: */
		    f[ij] -= drep * .5 * p[ij];
		}
	    }
	    i__4 = ic;
	    for (i = ia; i <= i__4; ++i) {
		ka = ifact[i - 1];
		i__3 = jb;
		for (j = jc + 1; j <= i__3; ++j) {
		    ij = ka + j;

/*    ATOM J (D(IF PRESENT)) EXCHANGE WITH ATOM I (S AND P
 ONLY) */

/* L100: */
		    f[ij] -= drep * .5 * p[ij];
		}
	    }

/*                      THE COULOMB REPULSION TERMS. */

/*     FIRST, ATOM J (S, P AND D SHELLS) BEING REPELLED BY ATOM I(
DSHELL) */

	    i__3 = jb;
	    for (j = ja; j <= i__3; ++j) {
		j2 = i1fact[j - 1];
/* L110: */
		f[j2] += drep * dpop[ii - 1];
	    }

/*     ATOM J (D SHELL) BEING REPELLED BY ATOM I (S AND P SHELLS) 
*/

	    i__3 = jb;
	    for (j = jc + 1; j <= i__3; ++j) {
		j2 = i1fact[j - 1];
/* L120: */
		f[j2] += drep * sppop[ii - 1];
	    }

/*     ATOM I (S, P AND D SHELLS) BEING REPELLED BY ATOM J (D SHEL
L) */

	    i__3 = ib;
	    for (i = ia; i <= i__3; ++i) {
		i2 = i1fact[i - 1];
/* L130: */
		f[i2] += drep * dpop[jj - 1];
	    }

/*    ATOM I (D SHELL) BEING REPELLED BY ATOM J (S AND P SHELLS) 
*/

	    i__3 = ib;
	    for (i = ic + 1; i <= i__3; ++i) {
		i2 = i1fact[i - 1];
/* L140: */
		f[i2] += drep * sppop[jj - 1];
	    }
/* L150: */
	}
/* L160: */
    }

    return 0;
L170:
    kr = 0;
    i__1 = *numat;
    for (ii = 1; ii <= i__1; ++ii) {
	ia = nfirst[ii];
	ib = nlast[ii];
	im1 = ii - ione;
	i__2 = im1;
	for (jj = 1; jj <= i__2; ++jj) {
	    ++kr;
	    if (lid) {
		elrep = w[kr];
		elexc = elrep;
	    } else {
		elrep = wj[kr];
		elexc = wk[kr];
	    }
	    ja = nfirst[jj];
	    jb = nlast[jj];
	    i__3 = ib;
	    for (i = ia; i <= i__3; ++i) {
		ka = ifact[i - 1];
		kk = ka + i;
		i__4 = jb;
		for (k = ja; k <= i__4; ++k) {
		    ll = i1fact[k - 1];
		    ik = ka + k;
		    f[kk] += ptot[ll] * elrep;
		    f[ll] += ptot[kk] * elrep;
/* L180: */
		    f[ik] -= p[ik] * elexc;
		}
	    }
/* L190: */
	}
/* L200: */
    }
    return 0;
} /* fock2d_ */