/* dcart.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 p[23220], pa[23220], pb[23220];
} densty_;

#define densty_1 densty_

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 htype[4];
    integer nhco[80]	/* was [4][20] */, nnhco, itype;
    logical usemm;
} molmec_;

#define molmec_1 molmec_

struct {
    integer l1l, l2l, l3l, l1u, l2u, l3u;
} ucell_;

#define ucell_1 ucell_

struct {
    integer k1l, k2l, k3l, k1u, k2u, k3u;
} dcartc_;

#define dcartc_1 dcartc_

struct {
    doublereal uss[107], upp[107], udd[107];
} onelec_;

#define onelec_1 onelec_

struct {
    integer numcal;
} numcal_;

#define numcal_1 numcal_

/* Table of constant values */

static integer c__1 = 1;
static doublereal c_b36 = 100.;
static integer c__2 = 2;

/* Subroutine */ int dcart_(doublereal *coord, doublereal *dxyz)
{
    /* Initialized data */

    static doublereal chnge = 1e-4;
    static doublereal chnge2 = 5e-5;
    static logical first = TRUE_;

    /* System generated locals */
    integer i__1, i__2, i__3, i__4, i__5, i__6, i__7;
    doublereal d__1;
    alist al__1;

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

    /* Local variables */
    static doublereal padi[171], pbdi[171], heat, refh;
    static integer krep, i, j, k, l;
    extern /* Subroutine */ int dihed_(doublereal *, integer *, integer *, 
	    integer *, integer *, doublereal *);
    static logical debug;
    static doublereal angle;
    static logical force, makep;
    static doublereal deriv, aa, ee;
    static integer if__, jf, ii;
#define lstor1 ((integer *)&ucell_1)
#define lstor2 ((integer *)&dcartc_1)
    static integer im, il, jj, jm;
    static logical anader;
    static integer jl, ik, jk, kl, ij, ncells, iofset;
    extern /* Subroutine */ int analyt_(doublereal *, doublereal *, 
	    doublereal *, doublereal *, integer *, integer *, integer *, 
	    integer *, integer *, integer *, doublereal *, integer *);
    static integer im1, numtot;
    extern /* Subroutine */ int dhc_(doublereal *, doublereal *, doublereal *,
	     doublereal *, integer *, integer *, integer *, integer *, 
	    integer *, integer *, integer *, integer *, doublereal *);
    static doublereal cdi[6]	/* was [3][2] */, del, eng[3];
    static integer ndi[2], iii;
    static doublereal pdi[171];
    static integer jjj;
    static doublereal sum;

    /* Fortran I/O blocks */
    static cilist io___47 = { 0, 6, 0, "(//10X,'CARTESIAN COORDINATE DERIVAT"
	    "IVES',//3X,        'ATOM  AT. NO.',5X,'X',12X,'Y',12X,'Z',/)", 0 }
	    ;
    static cilist io___48 = { 0, 6, 0, "(2I6,F13.6,2F13.6)", 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 */

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

/*    DCART CALCULATES THE DERIVATIVES OF THE ENERGY WITH RESPECT TO THE 
*/
/*          CARTESIAN COORDINATES. THIS IS DONE BY FINITE DIFFERENCES. */

/*    THE MAIN ARRAYS IN DCART ARE: */
/*        DXYZ   ON EXIT CONTAINS THE CARTESIAN DERIVATIVES. */

/* ***********************************************************************
 */
    /* Parameter adjustments */
    dxyz -= 4;
    coord -= 4;

    /* Function Body */

/* CHNGE IS A MACHINE-PRECISION DEPENDENT CONSTANT */
/* CHNGE2=CHNGE/2 */

    if (first) {
	anader = i_indx(keywrd_1.keywrd, "ANALYT", 80L, 6L) != 0;
	debug = i_indx(keywrd_1.keywrd, "DCART", 80L, 5L) != 0;
	force = i_indx(keywrd_1.keywrd, "PRECISE", 80L, 7L) + i_indx(
		keywrd_1.keywrd, "FORCE", 80L, 5L) != 0;
	first = FALSE_;
    }
    ncells = (ucell_1.l1u - ucell_1.l1l + 1) * (ucell_1.l2u - ucell_1.l2l + 1)
	     * (ucell_1.l3u - ucell_1.l3l + 1);
    for (i = 1; i <= 6; ++i) {
	lstor2[i - 1] = lstor1[i - 1];
/* L10: */
	lstor1[i - 1] = 0;
    }
    iofset = (ncells + 1) / 2;
    numtot = molkst_1.numat * ncells;
    i__1 = numtot;
    for (i = 1; i <= i__1; ++i) {
	for (j = 1; j <= 3; ++j) {
/* L20: */
	    dxyz[j + i * 3] = 0.;
	}
    }
    if (anader) {
	al__1.aerr = 0;
	al__1.aunit = 2;
	f_rew(&al__1);
    }
    krep = 0;
    i__1 = molkst_1.numat;
    for (ii = 1; ii <= i__1; ++ii) {
	iii = ncells * (ii - 1) + iofset;
	im1 = ii;
	if__ = molkst_1.nfirst[ii - 1];
	im = molkst_1.nmidle[ii - 1];
	il = molkst_1.nlast[ii - 1];
	ndi[1] = molkst_1.nat[ii - 1];
	for (i = 1; i <= 3; ++i) {
/* L30: */
	    cdi[i + 2] = coord[i + ii * 3];
	}
	i__2 = im1;
	for (jj = 1; jj <= i__2; ++jj) {
	    jjj = ncells * (jj - 1);
/*  FORM DIATOMIC MATRICES */
	    jf = molkst_1.nfirst[jj - 1];
	    jm = molkst_1.nmidle[jj - 1];
	    jl = molkst_1.nlast[jj - 1];
/*   GET FIRST ATOM */
	    ndi[0] = molkst_1.nat[jj - 1];
	    makep = TRUE_;
	    i__3 = dcartc_1.k1u;
	    for (ik = dcartc_1.k1l; ik <= i__3; ++ik) {
		i__4 = dcartc_1.k2u;
		for (jk = dcartc_1.k2l; jk <= i__4; ++jk) {
		    i__5 = dcartc_1.k3u;
		    for (kl = dcartc_1.k3l; kl <= i__5; ++kl) {
			++jjj;
			for (l = 1; l <= 3; ++l) {
/* L40: */
			    cdi[l - 1] = coord[l + jj * 3] + euler_1.tvec[l - 
				    1] * ik + euler_1.tvec[l + 2] * jk + 
				    euler_1.tvec[l + 5] * kl;
			}
			if (! makep) {
			    goto L90;
			}
			makep = FALSE_;
			ij = 0;
			i__6 = jl;
			for (i = jf; i <= i__6; ++i) {
			    k = i * (i - 1) / 2 + jf - 1;
			    i__7 = i;
			    for (j = jf; j <= i__7; ++j) {
				++ij;
				++k;
				padi[ij - 1] = densty_1.pa[k - 1];
				pbdi[ij - 1] = densty_1.pb[k - 1];
/* L50: */
				pdi[ij - 1] = densty_1.p[k - 1];
			    }
			}
/* GET SECOND ATOM FIRST ATOM INTERSECTION */
			i__7 = il;
			for (i = if__; i <= i__7; ++i) {
			    l = i * (i - 1) / 2;
			    k = l + jf - 1;
			    i__6 = jl;
			    for (j = jf; j <= i__6; ++j) {
				++ij;
				++k;
				padi[ij - 1] = densty_1.pa[k - 1];
				pbdi[ij - 1] = densty_1.pb[k - 1];
/* L60: */
				pdi[ij - 1] = densty_1.p[k - 1];
			    }
			    k = l + if__ - 1;
			    i__6 = i;
			    for (l = if__; l <= i__6; ++l) {
				++k;
				++ij;
				padi[ij - 1] = densty_1.pa[k - 1];
				pbdi[ij - 1] = densty_1.pb[k - 1];
/* L70: */
				pdi[ij - 1] = densty_1.p[k - 1];
			    }
/* L80: */
			}
L90:
			if (ii == jj) {
			    goto L120;
			}
			if (anader) {
			    analyt_(pdi, padi, pbdi, cdi, ndi, &jf, &jl, &
				    if__, &il, &molkst_1.norbs, eng, &krep);
			    for (k = 1; k <= 3; ++k) {
				dxyz[k + iii * 3] += eng[k - 1];
/* L100: */
				dxyz[k + jjj * 3] -= eng[k - 1];
			    }
			} else {
			    if (! force) {
				cdi[0] += chnge2;
				cdi[1] += chnge2;
				cdi[2] += chnge2;
				dhc_(pdi, padi, pbdi, cdi, ndi, &jf, &jm, &jl,
					 &if__, &im, &il, &molkst_1.norbs, &
					aa);
			    }
			    for (k = 1; k <= 3; ++k) {
				if (force) {
				    cdi[k + 2] -= chnge2;
				    dhc_(pdi, padi, pbdi, cdi, ndi, &jf, &jm, 
					    &jl, &if__, &im, &il, &
					    molkst_1.norbs, &aa);
				}
				cdi[k + 2] += chnge;
				dhc_(pdi, padi, pbdi, cdi, ndi, &jf, &jm, &jl,
					 &if__, &im, &il, &molkst_1.norbs, &
					ee);
				cdi[k + 2] -= chnge2;
				if (! force) {
				    cdi[k + 2] -= chnge2;
				}
				deriv = (aa - ee) * 46.122 / chnge;
				dxyz[k + iii * 3] += deriv;
				dxyz[k + jjj * 3] -= deriv;
/* L110: */
			    }
			}
L120:
			;
		    }
		}
	    }
/* L130: */
	}
    }
    if (molmec_1.usemm) {

/*   NOW ADD IN MOLECULAR-MECHANICS CORRECTION TO THE H-N-C=O TORSION 
*/

	del = 1e-5;
	i__2 = molmec_1.nnhco;
	for (i = 1; i <= i__2; ++i) {
	    for (j = 1; j <= 4; ++j) {
		for (k = 1; k <= 3; ++k) {
		    coord[k + molmec_1.nhco[j + (i << 2) - 5] * 3] -= del;
		    dihed_(&coord[4], &molmec_1.nhco[(i << 2) - 4], &
			    molmec_1.nhco[(i << 2) - 3], &molmec_1.nhco[(i << 
			    2) - 2], &molmec_1.nhco[(i << 2) - 1], &angle);
/* Computing 2nd power */
		    d__1 = sin(angle);
		    refh = molmec_1.htype[molmec_1.itype - 1] * (d__1 * d__1);
		    coord[k + molmec_1.nhco[j + (i << 2) - 5] * 3] += del * 
			    2.;
		    dihed_(&coord[4], &molmec_1.nhco[(i << 2) - 4], &
			    molmec_1.nhco[(i << 2) - 3], &molmec_1.nhco[(i << 
			    2) - 2], &molmec_1.nhco[(i << 2) - 1], &angle);
		    coord[k + molmec_1.nhco[j + (i << 2) - 5] * 3] -= del;
/* Computing 2nd power */
		    d__1 = sin(angle);
		    heat = molmec_1.htype[molmec_1.itype - 1] * (d__1 * d__1);
		    sum = (refh - heat) / del;
		    dxyz[k + molmec_1.nhco[j + (i << 2) - 5] * 3] += sum;
/* L140: */
		}
/* L150: */
	    }
/* L160: */
	}
    }
    for (i = 1; i <= 6; ++i) {
/* L170: */
	lstor1[i - 1] = lstor2[i - 1];
    }
    if (! debug) {
	return 0;
    }
    s_wsfe(&io___47);
    e_wsfe();
    s_wsfe(&io___48);
    i__2 = numtot;
    for (i = 1; i <= i__2; ++i) {
	do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&molkst_1.nat[(i - 1) / ncells], (ftnlen)sizeof(
		integer));
	for (j = 1; j <= 3; ++j) {
	    do_fio(&c__1, (char *)&dxyz[j + i * 3], (ftnlen)sizeof(doublereal)
		    );
	}
    }
    e_wsfe();
    if (anader) {
	al__1.aerr = 0;
	al__1.aunit = 2;
	f_rew(&al__1);
    }
    return 0;
} /* dcart_ */

#undef lstor2
#undef lstor1


/* Subroutine */ int dhc_(doublereal *p, doublereal *pa, doublereal *pb, 
	doublereal *xi, integer *nat, integer *if__, integer *im, integer *il,
	 integer *jf, integer *jm, integer *jl, integer *norbs, doublereal *
	dener)
{
    /* Initialized data */

    static integer icalcn = 0;

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

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

    /* Local variables */
    static doublereal wlim, f[171], h[171];
    static integer i, j, k;
    static doublereal w[100], shmat[81]	/* was [9][9] */;
    static integer nlast[2], i1, j1, i2;
    extern /* Subroutine */ int h1elec_(integer *, integer *, doublereal *, 
	    doublereal *, doublereal *), fock2d_(doublereal *, doublereal *, 
	    doublereal *, doublereal *, doublereal *, doublereal *, integer *,
	     integer *, integer *, integer *);
    static integer ia, ja, jb, jc, ib, ic, nj, ni, jj, kr, jt, ii;
    static doublereal ee, wj[100], wk[100];
    extern doublereal helect_(integer *, doublereal *, doublereal *, 
	    doublereal *);
    static integer nmidle[2], linear;
    static doublereal e1b[10], e2a[10], enuclr;
    extern /* Subroutine */ int rotate_(integer *, integer *, doublereal *, 
	    doublereal *, doublereal *, integer *, doublereal *, doublereal *,
	     doublereal *, doublereal *);
    static integer nfirst[2];
    extern /* Subroutine */ int solrot_(integer *, integer *, doublereal *, 
	    doublereal *, doublereal *, doublereal *, integer *, doublereal *,
	     doublereal *, doublereal *, doublereal *);
    static logical uhf;

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

/*  DHC CALCULATES THE ENERGY CONTRIBUTIONS FROM THOSE PAIRS OF ATOMS */
/*         THAT HAVE BEEN MOVED BY SUBROUTINE DERIV. */

/* ***********************************************************************
 */
    /* Parameter adjustments */
    --nat;
    xi -= 4;
    --pb;
    --pa;
    --p;

    /* Function Body */
    if (icalcn != numcal_1.numcal) {
	icalcn = numcal_1.numcal;
	wlim = 4.;
	if (euler_1.id == 0) {
	    wlim = 0.;
	}
	uhf = i_indx(keywrd_1.keywrd, "UHF", 80L, 3L) != 0;
    }
    nfirst[0] = 1;
    nmidle[0] = *im - *if__ + 1;
    nlast[0] = *il - *if__ + 1;
    nfirst[1] = nlast[0] + 1;
    nmidle[1] = nfirst[1] + *jm - *jf;
    nlast[1] = nfirst[1] + *jl - *jf;
    linear = nlast[1] * (nlast[1] + 1) / 2;
    i__1 = linear;
    for (i = 1; i <= i__1; ++i) {
	f[i - 1] = 0.;
/* L10: */
	h[i - 1] = 0.;
    }
    i__1 = linear;
    for (i = 1; i <= i__1; ++i) {
/* L20: */
	f[i - 1] = h[i - 1];
    }
    ja = nfirst[1];
    jb = nlast[1];
    jc = nmidle[1];
    ia = nfirst[0];
    ib = nlast[0];
    ic = nmidle[0];
    jt = jb * (jb + 1) / 2;
    j = 2;
    i = 1;
    nj = nat[2];
    ni = nat[1];
    h1elec_(&ni, &nj, &xi[4], &xi[7], shmat);
    if (nat[1] == 102 || nat[2] == 102) {
	k = jb * (jb + 1) / 2;
	i__1 = k;
	for (j = 1; j <= i__1; ++j) {
/* L30: */
	    h[j - 1] = 0.;
	}
    } else {
	j1 = 0;
	i__1 = jb;
	for (j = ja; j <= i__1; ++j) {
	    jj = j * (j - 1) / 2;
	    ++j1;
	    i1 = 0;
	    i__2 = ib;
	    for (i = ia; i <= i__2; ++i) {
		++jj;
		++i1;
		h[jj - 1] = shmat[i1 + j1 * 9 - 10];
		f[jj - 1] = shmat[i1 + j1 * 9 - 10];
/* L40: */
	    }
	}
    }
    kr = 1;
    if (euler_1.id == 0) {
	rotate_(&nj, &ni, &xi[7], &xi[4], &w[kr - 1], &kr, e2a, e1b, &enuclr, 
		&c_b36);
    } else {
	solrot_(&nj, &ni, &xi[7], &xi[4], wj, wk, &kr, e2a, e1b, &enuclr, &
		c_b36);
    }
    if (wj[0] < wlim) {
	i__2 = kr - 1;
	for (i = 1; i <= i__2; ++i) {
/* L50: */
	    wk[i - 1] = 0.;
	}
    }

/*    * ENUCLR IS SUMMED OVER CORE-CORE REPULSION INTEGRALS. */

    i2 = 0;
    i__2 = ic;
    for (i1 = ia; i1 <= i__2; ++i1) {
	ii = i1 * (i1 - 1) / 2 + ia - 1;
	i__1 = i1;
	for (j1 = ia; j1 <= i__1; ++j1) {
	    ++ii;
	    ++i2;
	    h[ii - 1] += e1b[i2 - 1];
/* L60: */
	    f[ii - 1] += e1b[i2 - 1];
	}
    }
    i__1 = ib;
    for (i1 = ic + 1; i1 <= i__1; ++i1) {
	ii = i1 * (i1 + 1) / 2;
	f[ii - 1] += e1b[0];
/* L70: */
	h[ii - 1] += e1b[0];
    }
    i2 = 0;
    i__1 = jc;
    for (i1 = ja; i1 <= i__1; ++i1) {
	ii = i1 * (i1 - 1) / 2 + ja - 1;
	i__2 = i1;
	for (j1 = ja; j1 <= i__2; ++j1) {
	    ++ii;
	    ++i2;
	    h[ii - 1] += e2a[i2 - 1];
/* L80: */
	    f[ii - 1] += e2a[i2 - 1];
	}
    }
    i__2 = jb;
    for (i1 = jc + 1; i1 <= i__2; ++i1) {
	ii = i1 * (i1 + 1) / 2;
	f[ii - 1] += e2a[0];
/* L90: */
	h[ii - 1] += e2a[0];
    }
    fock2d_(f, &p[1], &pa[1], w, wj, wk, &c__2, nfirst, nmidle, nlast);
    ee = helect_(&nlast[1], &pa[1], h, f);
    if (uhf) {
	i__2 = linear;
	for (i = 1; i <= i__2; ++i) {
/* L100: */
	    f[i - 1] = h[i - 1];
	}
	fock2d_(f, &p[1], &pb[1], w, wj, wk, &c__2, nfirst, nmidle, nlast);
	ee += helect_(&nlast[1], &pb[1], h, f);
    } else {
	ee *= 2.;
    }
    *dener = ee + enuclr;
    return 0;

} /* dhc_ */