/* compfg.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 nvar, loc[516]	/* was [2][258] */, idumy;
    doublereal dumy[258];
} geovar_;

#define geovar_1 geovar_

struct {
    integer ndep, locpar[258], idepfn[258], locdep[258];
} geosym_;

#define geosym_1 geosym_

struct {
    doublereal geo[258]	/* was [3][86] */;
} geom_;

#define geom_1 geom_

struct {
    doublereal atheat;
} atheat_;

#define atheat_1 atheat_

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

#define wmatrx_1 wmatrx_

struct {
    doublereal enuclr;
} enuclr_;

#define enuclr_1 enuclr_

struct {
    integer nztype[107], mtype[30], ltype;
} natype_;

#define natype_1 natype_

struct {
    doublereal elect;
} elect_;

#define elect_1 elect_

struct {
    doublereal rxyz[23220], xdumy[43344];
} scrach_;

#define scrach_1 scrach_

struct {
    doublereal h[23220];
} hmatrx_;

#define hmatrx_1 hmatrx_

struct molmec_1_ {
    doublereal htype[4];
    integer nhco[80]	/* was [4][20] */, nnhco, itype;
    logical usemm;
};

#define molmec_1 (*(struct molmec_1_ *) &molmec_)

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 {
    char keywrd[80];
} keywrd_;

#define keywrd_1 keywrd_

/* Initialized data */

struct {
    doublereal e_1[4];
    integer fill_2[83];
    } molmec_ = { 6.1737, 3.3191, 7.1853, 1.7712 };


/* Table of constant values */

static integer c__1 = 1;
static logical c_true = TRUE_;

/* Subroutine */ int compfg_(doublereal *xparam, logical *int__, doublereal *
	escf, logical *fulscf, doublereal *grad, logical *lgrad)
{
    /* Initialized data */

    static logical first = TRUE_;

    /* System generated locals */
    integer i__1, i__2;
    doublereal d__1;
    alist al__1;

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

    /* Local variables */
    extern /* Subroutine */ int iter_(doublereal *, doublereal *, real *, 
	    real *, doublereal *, logical *, logical *);
    static integer i, j, k, l;
    extern /* Subroutine */ int dihed_(doublereal *, integer *, integer *, 
	    integer *, integer *, doublereal *);
    static logical debug;
#define w ((doublereal *)&wmatrx_1)
    static doublereal angle;
    static logical large;
    extern /* Subroutine */ int hcore_(doublereal *, doublereal *, doublereal 
	    *, real *, real *, doublereal *);
    static doublereal coord[258]	/* was [3][86] */;
    extern /* Subroutine */ int deriv_(doublereal *, doublereal *);
    static logical print;
    static doublereal degree[3];
    static logical analyt;
    extern /* Subroutine */ int setupg_(void), gmetry_(doublereal *, 
	    doublereal *), symtry_(void);

    /* Fortran I/O blocks */
    static cilist io___12 = { 0, 6, 0, "(' INTERNAL COORDS',/100(/,3F12.6))", 
	    0 };
    static cilist io___14 = { 0, 6, 0, "(' CARTESIAN COORDS',/100(/,3F12.6))",
	     0 };
    static cilist io___16 = { 0, 6, 0, "(/10X,' HEAT OF FORMATION',G30.17)", 
	    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 */

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

/*   COMPFG CALCULATES (A) THE HEAT OF FORMATION OF THE SYSTEM, AND */
/*                     (B) THE GRADIENTS, IF LGRAD IS .TRUE. */

/*   ON INPUT  XPARAM = ARRAY OF PARAMETERS TO BE USED IN INTERNAL COORDS 
*/
/*             LGRAD  = .TRUE. IF GRADIENTS ARE NEEDED, .FALSE. OTHERWISE 
*/
/*             INT    = NEVER USED, RESERVED FOR FUTURE USE */
/*             FULSCF = .TRUE. IF FULL SCF TO BE DONE, .FALSE. OTHERWISE. 
*/

/*   ON OUTPUT ESCF  = HEAT OF FORMATION. */
/*             GRAD   = ARRAY OF GRADIENTS, IF LGRAD = .TRUE. */

/* ***********************************************************************
 */
    /* Parameter adjustments */
    --grad;
    --xparam;

    /* Function Body */
/*                 MNDO     AM1      PM3      MINDO/3 */
    if (first) {
	first = FALSE_;
/* #         WRITE(6,'(A,F12.4)')' MM CORRECTION:',HTYPE(ITYPE) */
	natype_1.ltype = 0;
	i__1 = molkst_1.numat;
	for (i = 1; i <= i__1; ++i) {
	    if (molkst_1.nat[i - 1] < 99) {
		i__2 = natype_1.ltype;
		for (j = 1; j <= i__2; ++j) {
/* L10: */
		    if (molkst_1.nat[i - 1] == natype_1.mtype[j - 1]) {
			goto L20;
		    }
		}
		++natype_1.ltype;
		natype_1.mtype[natype_1.ltype - 1] = molkst_1.nat[i - 1];
		natype_1.nztype[molkst_1.nat[i - 1] - 1] = natype_1.ltype;

/*       LTYPE = NUMBER OF TYPES OF REAL ATOM PRESENT */
/*       MTYPE = TYPES OF REAL ATOMS PRESENT */
		j = natype_1.ltype;
L20:
		;
	    }
/* L30: */
	}
	analyt = i_indx(keywrd_1.keywrd, "ANALYT", 80L, 6L) != 0;
	if (analyt) {
	    setupg_();
	}
	degree[0] = 1.;
	if (i_indx(keywrd_1.keywrd, " XYZ", 80L, 4L) != 0) {
	    degree[1] = 1.;
	} else {
	    degree[1] = 57.295779531334603;
	}
	degree[2] = degree[1];
	large = i_indx(keywrd_1.keywrd, "LARGE", 80L, 5L) != 0;
	print = i_indx(keywrd_1.keywrd, "COMPFG", 80L, 6L) != 0;
	debug = i_indx(keywrd_1.keywrd, "DEBUG", 80L, 5L) != 0 && print;
    }

/* SET UP COORDINATES FOR CURRENT CALCULATION */

/*       PLACE THE NEW VALUES OF THE VARIABLES IN THE ARRAY GEO. */
/*       MAKE CHANGES IN THE GEOMETRY. */
    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];
/* L40: */
	geom_1.geo[l + k * 3 - 4] = xparam[i];
    }
/* #      WRITE(6,'(3F18.11)')(XPARAM(I),I=1,3) */
/*      IMPOSE THE SYMMETRY CONDITIONS + COMPUTE THE DEPENDENT-PARAMETERS 
*/
    if (geosym_1.ndep != 0) {
	symtry_();
    }
/*      NOW COMPUTE THE ATOMIC COORDINATES. */
    if (debug) {
	if (large) {
	    k = molkst_1.numat;
	} else {
	    k = min(5,molkst_1.numat);
	}
	s_wsfe(&io___12);
	i__1 = k;
	for (i = 1; i <= i__1; ++i) {
	    for (j = 1; j <= 3; ++j) {
		d__1 = geom_1.geo[j + i * 3 - 4] * degree[j - 1];
		do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof(doublereal));
	    }
	}
	e_wsfe();
    }
    gmetry_(geom_1.geo, coord);
    if (debug) {
	if (large) {
	    k = molkst_1.numat;
	} else {
	    k = min(5,molkst_1.numat);
	}
	s_wsfe(&io___14);
	i__1 = k;
	for (i = 1; i <= i__1; ++i) {
	    for (j = 1; j <= 3; ++j) {
		do_fio(&c__1, (char *)&coord[j + i * 3 - 4], (ftnlen)sizeof(
			doublereal));
	    }
	}
	e_wsfe();
    }
    if (analyt) {
	al__1.aerr = 0;
	al__1.aunit = 2;
	f_rew(&al__1);
    }
    hcore_(coord, hmatrx_1.h, w, wmatrx_1.wj, wmatrx_1.wk, &enuclr_1.enuclr);

/* COMPUTE THE HEAT OF FORMATION. */

    if (molkst_1.norbs > 0 && molkst_1.nelecs > 0) {
	iter_(hmatrx_1.h, w, wmatrx_1.wj, wmatrx_1.wk, &elect_1.elect, fulscf,
		 &c_true);
    } else {
	elect_1.elect = 0.;
    }
    *escf = (elect_1.elect + enuclr_1.enuclr) * 23.061 + atheat_1.atheat;
    i__1 = molmec_1.nnhco;
    for (i = 1; i <= i__1; ++i) {
	dihed_(coord, &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);
	*escf += molmec_1.htype[molmec_1.itype - 1] * (d__1 * d__1);
/* L50: */
    }
    if (print) {
	s_wsfe(&io___16);
	do_fio(&c__1, (char *)&(*escf), (ftnlen)sizeof(doublereal));
	e_wsfe();
    }

/* FIND DERIVATIVES IF DESIRED */

    if (*lgrad) {
	deriv_(geom_1.geo, &grad[1]);
    }
    return 0;
} /* compfg_ */

#undef w