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

#define euler_1 euler_

struct {
    integer nvar, loc[516]	/* was [2][258] */, idumy;
    doublereal dummy[258];
} geovar_;

#define geovar_1 geovar_

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 {
    integer natoms, labels[86], na[86], nb[86], nc[86];
} geokst_;

#define geokst_1 geokst_

struct {
    doublereal cosine;
} gravec_;

#define gravec_1 gravec_

struct {
    integer ndep, idumys[774]	/* was [258][3] */;
} geosym_;

#define geosym_1 geosym_

struct {
    integer latom, lparam;
    doublereal react[200];
} path_;

#define path_1 path_

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

#define ucell_1 ucell_

struct {
    doublereal dxyz[6966]	/* was [3][2322] */;
} xyzgra_;

#define xyzgra_1 xyzgra_

struct {
    doublereal enuclr;
} enuclr_;

#define enuclr_1 enuclr_

struct {
    integer numcal;
} numcal_;

#define numcal_1 numcal_

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

#define densty_1 densty_

struct {
    real wj[219386], wk[219386];
} wmatrx_;

#define wmatrx_1 wmatrx_

struct {
    doublereal h[23220];
} hmatrx_;

#define hmatrx_1 hmatrx_

struct {
    doublereal atheat;
} atheat_;

#define atheat_1 atheat_

struct {
    char keywrd[80];
} keywrd_;

#define keywrd_1 keywrd_

struct {
    doublereal errfn[258];
} errfn_;

#define errfn_1 errfn_

/* Table of constant values */

static doublereal c_b11 = 10.;
static integer c__1 = 1;
static logical c_true = TRUE_;
static logical c_false = FALSE_;

/* Subroutine */ int deriv_(doublereal *geo, doublereal *grad)
{
    /* Initialized data */

    static integer icalcn = 0;

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

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

    /* Local variables */
    static doublereal cold[6966]	/* was [3][2322] */;
    static integer idlo;
    static logical fast;
    extern /* Subroutine */ int iter_(doublereal *, doublereal *, real *, 
	    real *, doublereal *, logical *, logical *);
    static doublereal sum11, sum22, xjuc[3], sum12, totl, time1, totl1;
    static integer i, j, k;
    extern doublereal reada_(char *, integer *, ftnlen);
    static integer l;
    static logical halfe, debug;
#define w ((doublereal *)&wmatrx_1)
    extern /* Subroutine */ int dcart_(doublereal *, doublereal *), hcore_(
	    doublereal *, doublereal *, doublereal *, real *, real *, 
	    doublereal *);
    static doublereal coord[258]	/* was [3][86] */, grlim;
    static logical times;
    static doublereal gnorm, press, aa, ee;
    static logical ci;
    static integer ij;
    static doublereal change[3];
    static integer ii, il, jl, kl, ll, kk, idelta, linear;
    extern doublereal second_(void);
    static logical precis;
    static doublereal xparam[258], xderiv[3];
    extern /* Subroutine */ int gmetry_(doublereal *, doublereal *);
    static doublereal xstore;
    extern /* Subroutine */ int symtry_(void);
    static logical scf1;

    /* Fortran I/O blocks */
    static cilist io___18 = { 0, 6, 0, "(' GEO AT START OF DERIV')", 0 };
    static cilist io___19 = { 0, 6, 0, "(F19.5,2F12.5)", 0 };
    static cilist io___24 = { 0, 6, 0, "(' DOING FULL SCF''S IN DERIV')", 0 };
    static cilist io___45 = { 0, 6, 0, "(' GRADIENTS')", 0 };
    static cilist io___46 = { 0, 6, 0, "(10F8.3)", 0 };
    static cilist io___47 = { 0, 6, 0, "(' ERROR FUNCTION')", 0 };
    static cilist io___48 = { 0, 6, 0, "(10F8.3)", 0 };
    static cilist io___49 = { 0, 6, 0, "(' COSINE OF SEARCH DIRECTION =',F30"
	    ".6)", 0 };
    static cilist io___50 = { 0, 6, 0, "(' TIME FOR DERIVATIVES',F12.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 */

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

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

/*    THE MAIN ARRAYS IN DERIV ARE: */
/*        LOC    INTEGER ARRAY, LOC(1,I) CONTAINS THE ADDRESS OF THE ATOM 
*/
/*               INTERNAL COORDINATE LOC(2,I) IS TO BE USED IN THE */
/*               DERIVATIVE CALCULATION. */
/*        GEO    ARRAY \GEO\ HOLDS THE INTERNAL COORDINATES. */
/*        GRAD   ON EXIT, CONTAINS THE DERIVATIVES */

/* ***********************************************************************
 */
    /* Parameter adjustments */
    --grad;
    geo -= 4;

    /* Function Body */
    if (icalcn != numcal_1.numcal) {
	i = i_indx(keywrd_1.keywrd, "PRESS", 80L, 5L);
	press = 0.;
	if (i != 0) {
	    press = reada_(keywrd_1.keywrd, &i, 80L) * 1476.8992;
	}
	idlo = geokst_1.natoms + 1;
	if (geokst_1.labels[geokst_1.natoms - 1] == 107) {
	    idlo = geokst_1.natoms;
	    if (geokst_1.labels[geokst_1.natoms - 2] == 107) {
		idlo = geokst_1.natoms - 1;
		if (geokst_1.labels[geokst_1.natoms - 3] == 107) {
		    idlo = geokst_1.natoms - 2;
		}
	    }
	}
	grlim = .01;
	debug = i_indx(keywrd_1.keywrd, "DERIV", 80L, 5L) != 0;
	precis = i_indx(keywrd_1.keywrd, "PRECIS", 80L, 6L) != 0;
	times = i_indx(keywrd_1.keywrd, "TIME", 80L, 4L) != 0;
	ci = i_indx(keywrd_1.keywrd, "C.I.", 80L, 4L) != 0;
	scf1 = i_indx(keywrd_1.keywrd, "1SCF", 80L, 4L) != 0;
	icalcn = numcal_1.numcal;
	if (i_indx(keywrd_1.keywrd, "RESTART", 80L, 7L) == 0) {
	    i__1 = geovar_1.nvar;
	    for (i = 1; i <= i__1; ++i) {
/* L10: */
		errfn_1.errfn[i - 1] = 0.;
	    }
	}
	grlim = .01;
	if (precis) {
	    grlim = 1e-4;
	}
	if (i_indx(keywrd_1.keywrd, "FULSCF", 80L, 6L) > 0) {
	    grlim = 1e9;
	}
	halfe = molkst_1.nopen > molkst_1.nclose || ci;
	idelta = -7;

/*   IDELTA IS A MACHINE-PRECISION DEPENDANT INTEGER */

	if (halfe && precis) {
	    idelta = -3;
	}
	if (halfe && ! precis) {
	    idelta = -3;
	}
	fast = TRUE_;
	change[0] = pow_di(&c_b11, &idelta);
	change[1] = pow_di(&c_b11, &idelta);
	change[2] = pow_di(&c_b11, &idelta);

/*    CHANGE(I) IS THE STEP SIZE USED IN CALCULATING THE DERIVATIVES. 
*/
/*    FOR "CARTESIAN" DERIVATIVES, CALCULATED USING DCART,AN */
/*    INFINITESIMAL STEP, HERE 0.000001, IS ACCEPTABLE. IN THE */
/*    HALF-ELECTRON METHOD A QUITE LARGE STEP IS NEEDED AS FULL SCF */
/*    CALCULATIONS ARE NEEDED, AND THE DIFFERENCE BETWEEN THE TOTAL */
/*    ENERGIES IS USED. THE STEP CANNOT BE VERY LARGE, AS THE SECOND 
*/
/*    DERIVITIVE IN FLEPO IS CALCULATED FROM THE DIFFERENCES OF TWO */
/*    FIRST DERIVATIVES. CHANGE(1) IS FOR CHANGE IN BOND LENGTH, */
/*    (2) FOR ANGLE, AND (3) FOR DIHEDRAL. */

	xderiv[0] = .5 / change[0];
	xderiv[1] = .5 / change[1];
	xderiv[2] = .5 / change[2];
    }
    gnorm = 0.;
    if (geovar_1.nvar == 0) {
	return 0;
    }
    if (debug) {
	s_wsfe(&io___18);
	e_wsfe();
	s_wsfe(&io___19);
	i__1 = geokst_1.natoms;
	for (i = 1; i <= i__1; ++i) {
	    for (j = 1; j <= 3; ++j) {
		do_fio(&c__1, (char *)&geo[j + i * 3], (ftnlen)sizeof(
			doublereal));
	    }
	}
	e_wsfe();
    }
    i__1 = geovar_1.nvar;
    for (i = 1; i <= i__1; ++i) {
	xparam[i - 1] = geo[geovar_1.loc[(i << 1) - 1] + geovar_1.loc[(i << 1)
		 - 2] * 3];
/* L20: */
/* Computing 2nd power */
	d__1 = grad[i];
	gnorm += d__1 * d__1;
    }
    gnorm = sqrt(gnorm);
    fast = gnorm > grlim && ! scf1 || ! halfe;
    time1 = second_();
    if (geosym_1.ndep != 0) {
	symtry_();
    }
    gmetry_(&geo[4], coord);
    if (! fast) {
	if (debug) {
	    s_wsfe(&io___24);
	    e_wsfe();
	}
	hcore_(coord, hmatrx_1.h, w, wmatrx_1.wj, wmatrx_1.wk, &
		enuclr_1.enuclr);
	if (molkst_1.norbs * molkst_1.nelecs > 0) {
	    iter_(hmatrx_1.h, w, wmatrx_1.wj, wmatrx_1.wk, &aa, &c_true, &
		    c_false);
	} else {
	    aa = 0.;
	}
	linear = molkst_1.norbs * (molkst_1.norbs + 1) / 2;
	i__1 = linear;
	for (i = 1; i <= i__1; ++i) {
/* L30: */
	    densty_1.p[i - 1] = densty_1.pa[i - 1] * 2.;
	}
	aa += enuclr_1.enuclr;
    }
    dcart_(coord, xyzgra_1.dxyz);
    if (geosym_1.ndep != 0) {
	symtry_();
    }
    gmetry_(&geo[4], coord);
    ij = 0;
    i__1 = molkst_1.numat;
    for (ii = 1; ii <= i__1; ++ii) {
	i__2 = ucell_1.l1u;
	for (il = ucell_1.l1l; il <= i__2; ++il) {
	    i__3 = ucell_1.l2u;
	    for (jl = ucell_1.l2l; jl <= i__3; ++jl) {
		i__4 = ucell_1.l3u;
		for (kl = ucell_1.l3l; kl <= i__4; ++kl) {
		    for (ll = 1; ll <= 3; ++ll) {
/* L40: */
			xjuc[ll - 1] = coord[ll + ii * 3 - 4] + euler_1.tvec[
				ll - 1] * il + euler_1.tvec[ll + 2] * jl + 
				euler_1.tvec[ll + 5] * kl;
		    }
		    ++ij;
		    for (kk = 1; kk <= 3; ++kk) {
			cold[kk + ij * 3 - 4] = xjuc[kk - 1];
/* L50: */
		    }
/* L60: */
		}
	    }
	}
/* L70: */
    }
    sum11 = 1e-9;
    sum22 = 1e-9;
    sum12 = 1e-9;
    i__1 = geovar_1.nvar;
    for (i = 1; i <= i__1; ++i) {
	k = geovar_1.loc[(i << 1) - 2];
	l = geovar_1.loc[(i << 1) - 1];
	xstore = xparam[i - 1];
	i__4 = geovar_1.nvar;
	for (j = 1; j <= i__4; ++j) {
/* L80: */
	    geo[geovar_1.loc[(j << 1) - 1] + geovar_1.loc[(j << 1) - 2] * 3] =
		     xparam[j - 1];
	}
	geo[l + k * 3] = xstore - change[l - 1];
	if (geosym_1.ndep != 0) {
	    symtry_();
	}
	gmetry_(&geo[4], coord);
/* #         CALL GEOUT */

/*    USE LOOKUP TABLE OF CARTESIAN DERIVATIVES TO WORK OUT INTERNAL 
*/
/*    COORDINATE DERIVATIVE. */

	totl = 0.;
	ij = 0;
	i__4 = molkst_1.numat;
	for (ii = 1; ii <= i__4; ++ii) {
	    if (euler_1.id == 0) {
		for (ll = 1; ll <= 3; ++ll) {
/* L90: */
		    totl += xyzgra_1.dxyz[ll + ii * 3 - 4] * (coord[ll + ii * 
			    3 - 4] - cold[ll + ii * 3 - 4]);
		}
	    } else {
		i__3 = ucell_1.l1u;
		for (il = ucell_1.l1l; il <= i__3; ++il) {
		    i__2 = ucell_1.l2u;
		    for (jl = ucell_1.l2l; jl <= i__2; ++jl) {
			i__5 = ucell_1.l3u;
			for (kl = ucell_1.l3l; kl <= i__5; ++kl) {
			    for (ll = 1; ll <= 3; ++ll) {
/* L100: */
				xjuc[ll - 1] = coord[ll + ii * 3 - 4] + 
					euler_1.tvec[ll - 1] * il + 
					euler_1.tvec[ll + 2] * jl + 
					euler_1.tvec[ll + 5] * kl;
			    }
			    ++ij;
			    for (kk = 1; kk <= 3; ++kk) {
				totl += xyzgra_1.dxyz[kk + ij * 3 - 4] * (
					xjuc[kk - 1] - cold[kk + ij * 3 - 4]);
/* L110: */
			    }
/* L120: */
			}
		    }
		}
	    }
/* L130: */
	}
	totl *= xderiv[l - 1];

/*   IF NEEDED, CALCULATE "EXACT" DERIVITIVES. */

	if (! fast) {
	    hcore_(coord, hmatrx_1.h, w, wmatrx_1.wj, wmatrx_1.wk, &
		    enuclr_1.enuclr);
	    if (molkst_1.norbs * molkst_1.nelecs > 0) {
		iter_(hmatrx_1.h, w, wmatrx_1.wj, wmatrx_1.wk, &ee, &c_true, &
			c_false);
	    } else {
		ee = 0.;
	    }
	    i__4 = linear;
	    for (ii = 1; ii <= i__4; ++ii) {
/* L140: */
		densty_1.p[ii - 1] = densty_1.pa[ii - 1] * 2.;
	    }
	    ee += enuclr_1.enuclr;
	    totl1 = (aa - ee) * 23.061 * xderiv[l - 1] * 2.;
/* #            WRITE(6,*)AA-EE */
	    errfn_1.errfn[i - 1] = totl1 - totl;
	}
	geo[l + k * 3] = xstore;
/* Computing 2nd power */
	d__1 = grad[i];
	sum11 += d__1 * d__1;
/* Computing 2nd power */
	d__1 = totl;
	sum22 += d__1 * d__1;
	sum12 += totl * grad[i];
	grad[i] = totl + errfn_1.errfn[i - 1];
/* L150: */
    }
    if (debug) {
	s_wsfe(&io___45);
	e_wsfe();
	s_wsfe(&io___46);
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    do_fio(&c__1, (char *)&grad[i], (ftnlen)sizeof(doublereal));
	}
	e_wsfe();
	s_wsfe(&io___47);
	e_wsfe();
	s_wsfe(&io___48);
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    do_fio(&c__1, (char *)&errfn_1.errfn[i - 1], (ftnlen)sizeof(
		    doublereal));
	}
	e_wsfe();
    }
    gravec_1.cosine = sum12 / sqrt(sum11 * sum22);
    if (debug) {
	s_wsfe(&io___49);
	do_fio(&c__1, (char *)&gravec_1.cosine, (ftnlen)sizeof(doublereal));
	e_wsfe();
    }
    if (! fast) {
	gravec_1.cosine = 1.;
    }
    if (times) {
	s_wsfe(&io___50);
	d__1 = second_() - time1;
	do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof(doublereal));
	e_wsfe();
    }
    return 0;
} /* deriv_ */

#undef w