/* powsav.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 {
    char elemnt[214];
} elemts_;

#define elemts_1 elemts_

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

#define geosym_1 geosym_

struct {
    char koment[80], title[80];
} titles_;

#define titles_1 titles_

struct {
    integer natoms, labels[86], na[86], nb[86], nc[86];
} geokst_;

#define geokst_1 geokst_

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

#define geom_1 geom_

struct {
    integer locvar[516]	/* was [2][258] */;
} locvar_;

#define locvar_1 locvar_

struct {
    char keywrd[80];
} keywrd_;

#define keywrd_1 keywrd_

struct {
    doublereal valvar[258];
    integer numvar;
} valvar_;

#define valvar_1 valvar_

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

#define densty_1 densty_

struct {
    doublereal alparm[774]	/* was [3][258] */, x0, x1, x2;
    integer jloop;
} alparm_;

#define alparm_1 alparm_

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 latom, lparam;
    doublereal react[200];
} path_;

#define path_1 path_

/* Table of constant values */

static integer c__4 = 4;
static integer c__3 = 3;
static integer c__1 = 1;
static doublereal c_b15 = 1.;
static integer c__9 = 9;

/* Subroutine */ int powsav_(doublereal *hess, doublereal *grad, doublereal *
	xparam, doublereal *pmat, integer *iloop, doublereal *bmat, integer *
	ipow)
{
    /* System generated locals */
    integer i__1, i__2;
    olist o__1;
    cllist cl__1;
    alist al__1;

    /* Builtin functions */
    integer f_open(olist *), f_rew(alist *), s_wsfe(cilist *), e_wsfe(void);
    double sqrt(doublereal);
    integer do_fio(integer *, char *, ftnlen), s_wsle(cilist *), e_wsle(void),
	     s_wsue(cilist *), do_uio(integer *, char *, ftnlen), e_wsue(void)
	    , f_clos(cllist *), s_rsue(cilist *), e_rsue(void);

    /* Local variables */
    static integer ivar, i, j, k, l;
    static doublereal x, coord[258]	/* was [3][86] */, funct1;
    static integer ir;
    static doublereal degree, qq[3];
    extern /* Subroutine */ int getarg_(integer *, char *, ftnlen);
    static integer linear;
    extern /* Subroutine */ int xyzint_(doublereal *, integer *, integer *, 
	    integer *, integer *, doublereal *, doublereal *);
    extern doublereal dot_(doublereal *, doublereal *, integer *);
    static integer iel1[3];

    /* Fortran I/O blocks */
    static cilist io___2 = { 0, 6, 0, "(//10X,'- - - - - - - TIME UP - - - -"
	    " - - -',//)", 0 };
    static cilist io___3 = { 0, 6, 0, "(//10X,' - THE CALCULATION IS BEING D"
	    "UMPED TO DISK',/10X,'   RESTART IT USING THE KEY-WORD \"RESTAR"
	    "T\"')", 0 };
    static cilist io___5 = { 0, 6, 0, "(//10X,'CURRENT VALUE OF GRADIENT NOR"
	    "M ='          ,F12.6)", 0 };
    static cilist io___9 = { 0, 6, 0, "(/10X,'CURRENT VALUE OF GEOMETRY',/)", 
	    0 };
    static cilist io___14 = { 0, 6, 0, "(A)", 0 };
    static cilist io___17 = { 0, 6, 0, "(2X,A2,3(F12.6,I3),I4,2I3)", 0 };
    static cilist io___19 = { 0, 6, 0, "(I4,3(F12.6,I3),I4,2I3)", 0 };
    static cilist io___20 = { 0, 6, 0, "(3(I4,','))", 0 };
    static cilist io___21 = { 0, 6, 0, 0, 0 };
    static cilist io___22 = { 0, 0, 0, 0, 0 };
    static cilist io___23 = { 0, 0, 0, 0, 0 };
    static cilist io___24 = { 0, 0, 0, 0, 0 };
    static cilist io___25 = { 0, 0, 0, 0, 0 };
    static cilist io___26 = { 0, 0, 0, 0, 0 };
    static cilist io___28 = { 0, 0, 0, 0, 0 };
    static cilist io___29 = { 0, 10, 0, 0, 0 };
    static cilist io___30 = { 0, 10, 0, 0, 0 };
    static cilist io___31 = { 0, 0, 0, 0, 0 };
    static cilist io___32 = { 0, 0, 0, 0, 0 };
    static cilist io___33 = { 0, 6, 0, "(//10X,' RESTORING DATA FROM DISK'/)",
	     0 };
    static cilist io___34 = { 0, 0, 0, 0, 0 };
    static cilist io___35 = { 0, 0, 0, 0, 0 };
    static cilist io___36 = { 0, 0, 0, 0, 0 };
    static cilist io___37 = { 0, 0, 0, 0, 0 };
    static cilist io___38 = { 0, 0, 0, 0, 0 };
    static cilist io___39 = { 0, 6, 0, "(10X,'FUNCTION =',F13.6//)", 0 };
    static cilist io___40 = { 0, 0, 0, 0, 0 };
    static cilist io___41 = { 0, 10, 0, 0, 0 };
    static cilist io___42 = { 0, 10, 0, 0, 0 };
    static cilist io___43 = { 0, 0, 0, 0, 0 };
    static cilist io___44 = { 0, 0, 0, 0, 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 */

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

/* POWSAV STORES AND RESTORES DATA USED IN THE SIGMA GEOMETRY */
/*        OPTIMISATION. */

/*  ON INPUT HESS   = HESSIAN MATRIX, PARTIAL OR WHOLE. */
/*           GRAD   = GRADIENTS. */
/*           XPARAM = CURRENT STATE OF PARAMETERS. */
/*           ILOOP  = INDEX OF HESSIAN, OR FLAG OF POINT REACHED SO-FAR. 
*/
/*           BMAT   = "B" MATRIX! */
/*           IPOW   = INDICES AND FLAGS. */
/*           IPOW(9)= 0 FOR RESTORE, 1 FOR DUMP */

/* ********************************************************************* 
*/
    /* Parameter adjustments */
    --ipow;
    bmat -= 259;
    --pmat;
    --xparam;
    --grad;
    hess -= 259;

    /* Function Body */
    getarg_(&c__4, argz_1.argz, 512L);
    o__1.oerr = 0;
    o__1.ounit = 9;
    o__1.ofnmlen = 512;
    o__1.ofnm = argz_1.argz;
    o__1.orl = 0;
    o__1.osta = "UNKNOWN";
    o__1.oacc = 0;
    o__1.ofm = "UNFORMATTED";
    o__1.oblnk = 0;
    f_open(&o__1);
    al__1.aerr = 0;
    al__1.aunit = 9;
    f_rew(&al__1);
    getarg_(&c__3, argz_1.argz, 512L);
    o__1.oerr = 0;
    o__1.ounit = 10;
    o__1.ofnmlen = 512;
    o__1.ofnm = argz_1.argz;
    o__1.orl = 0;
    o__1.osta = "UNKNOWN";
    o__1.oacc = 0;
    o__1.ofm = "UNFORMATTED";
    o__1.oblnk = 0;
    f_open(&o__1);
    al__1.aerr = 0;
    al__1.aunit = 10;
    f_rew(&al__1);
    ir = 9;
    if (ipow[9] != 0) {
	if (ipow[9] == 1) {
	    s_wsfe(&io___2);
	    e_wsfe();
	    s_wsfe(&io___3);
	    e_wsfe();
	    funct1 = sqrt(dot_(&grad[1], &grad[1], &geovar_1.nvar));
	    s_wsfe(&io___5);
	    do_fio(&c__1, (char *)&funct1, (ftnlen)sizeof(doublereal));
	    e_wsfe();
	    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];
/* L10: */
		geom_1.geo[l + k * 3 - 4] = xparam[i];
	    }
	    s_wsfe(&io___9);
	    e_wsfe();
	    degree = 57.29577951;
	    if (geokst_1.na[0] == 99) {

/*  CONVERT FROM CARTESIAN COORDINATES TO INTERNAL */

		i__1 = geokst_1.natoms;
		for (i = 1; i <= i__1; ++i) {
		    for (j = 1; j <= 3; ++j) {
/* L20: */
			coord[j + i * 3 - 4] = geom_1.geo[j + i * 3 - 4];
		    }
		}
		xyzint_(coord, &molkst_1.numat, geokst_1.na, geokst_1.nb, 
			geokst_1.nc, &c_b15, geom_1.geo);
	    }
	    geom_1.geo[1] = 0.;
	    geom_1.geo[2] = 0.;
	    geom_1.geo[0] = 0.;
	    geom_1.geo[4] = 0.;
	    geom_1.geo[5] = 0.;
	    geom_1.geo[8] = 0.;
	    ivar = 1;
	    geokst_1.na[0] = 0;
	    s_wsfe(&io___14);
	    do_fio(&c__1, keywrd_1.keywrd, 80L);
	    do_fio(&c__1, titles_1.koment, 80L);
	    do_fio(&c__1, titles_1.title, 80L);
	    e_wsfe();
	    i__1 = geokst_1.natoms;
	    for (i = 1; i <= i__1; ++i) {
		for (j = 1; j <= 3; ++j) {
/* L30: */
		    iel1[j - 1] = 0;
		}
L40:
		if (geovar_1.loc[(ivar << 1) - 2] == i) {
		    iel1[geovar_1.loc[(ivar << 1) - 1] - 1] = 1;
		    ++ivar;
		    goto L40;
		}
		if (i < 4) {
		    iel1[2] = 0;
		    if (i < 3) {
			iel1[1] = 0;
			if (i < 2) {
			    iel1[0] = 0;
			}
		    }
		}
		if (i == path_1.latom) {
		    iel1[path_1.lparam - 1] = -1;
		}
		qq[0] = geom_1.geo[i * 3 - 3];
		qq[1] = geom_1.geo[i * 3 - 2] * degree;
		qq[2] = geom_1.geo[i * 3 - 1] * degree;
/* L50: */
		s_wsfe(&io___17);
		do_fio(&c__1, elemts_1.elemnt + (geokst_1.labels[i - 1] - 1 <<
			 1), 2L);
		for (k = 1; k <= 3; ++k) {
		    do_fio(&c__1, (char *)&qq[k - 1], (ftnlen)sizeof(
			    doublereal));
		    do_fio(&c__1, (char *)&iel1[k - 1], (ftnlen)sizeof(
			    integer));
		}
		do_fio(&c__1, (char *)&geokst_1.na[i - 1], (ftnlen)sizeof(
			integer));
		do_fio(&c__1, (char *)&geokst_1.nb[i - 1], (ftnlen)sizeof(
			integer));
		do_fio(&c__1, (char *)&geokst_1.nc[i - 1], (ftnlen)sizeof(
			integer));
		e_wsfe();
	    }
	    i = 0;
	    x = 0.;
	    s_wsfe(&io___19);
	    do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&x, (ftnlen)sizeof(doublereal));
	    do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&x, (ftnlen)sizeof(doublereal));
	    do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&x, (ftnlen)sizeof(doublereal));
	    do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
	    e_wsfe();
	    if (geosym_1.ndep != 0) {
		i__1 = geosym_1.ndep;
		for (i = 1; i <= i__1; ++i) {
/* L60: */
		    s_wsfe(&io___20);
		    do_fio(&c__1, (char *)&geosym_1.locpar[i - 1], (ftnlen)
			    sizeof(integer));
		    do_fio(&c__1, (char *)&geosym_1.idepfn[i - 1], (ftnlen)
			    sizeof(integer));
		    do_fio(&c__1, (char *)&geosym_1.locdep[i - 1], (ftnlen)
			    sizeof(integer));
		    e_wsfe();
		}
		s_wsle(&io___21);
		e_wsle();
	    }
	}
	io___22.ciunit = ir;
	s_wsue(&io___22);
	do_uio(&c__9, (char *)&ipow[1], (ftnlen)sizeof(integer));
	do_uio(&c__1, (char *)&(*iloop), (ftnlen)sizeof(integer));
	e_wsue();
	io___23.ciunit = ir;
	s_wsue(&io___23);
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    do_uio(&c__1, (char *)&xparam[i], (ftnlen)sizeof(doublereal));
	}
	e_wsue();
	io___24.ciunit = ir;
	s_wsue(&io___24);
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    do_uio(&c__1, (char *)&grad[i], (ftnlen)sizeof(doublereal));
	}
	e_wsue();
	io___25.ciunit = ir;
	s_wsue(&io___25);
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    i__2 = geovar_1.nvar;
	    for (j = 1; j <= i__2; ++j) {
		do_uio(&c__1, (char *)&hess[j + i * 258], (ftnlen)sizeof(
			doublereal));
	    }
	}
	e_wsue();
	io___26.ciunit = ir;
	s_wsue(&io___26);
	i__2 = geovar_1.nvar;
	for (i = 1; i <= i__2; ++i) {
	    i__1 = geovar_1.nvar;
	    for (j = 1; j <= i__1; ++j) {
		do_uio(&c__1, (char *)&bmat[j + i * 258], (ftnlen)sizeof(
			doublereal));
	    }
	}
	e_wsue();
	linear = geovar_1.nvar * (geovar_1.nvar + 1) / 2;
	io___28.ciunit = ir;
	s_wsue(&io___28);
	i__1 = linear;
	for (i = 1; i <= i__1; ++i) {
	    do_uio(&c__1, (char *)&pmat[i], (ftnlen)sizeof(doublereal));
	}
	e_wsue();
	linear = molkst_1.norbs * (molkst_1.norbs + 1) / 2;
	s_wsue(&io___29);
	i__1 = linear;
	for (i = 1; i <= i__1; ++i) {
	    do_uio(&c__1, (char *)&densty_1.pa[i - 1], (ftnlen)sizeof(
		    doublereal));
	}
	e_wsue();
	if (molkst_1.nalpha != 0) {
	    s_wsue(&io___30);
	    i__1 = linear;
	    for (i = 1; i <= i__1; ++i) {
		do_uio(&c__1, (char *)&densty_1.pb[i - 1], (ftnlen)sizeof(
			doublereal));
	    }
	    e_wsue();
	}
	if (path_1.latom != 0) {
	    io___31.ciunit = ir;
	    s_wsue(&io___31);
	    i__1 = geovar_1.nvar;
	    for (i = 1; i <= i__1; ++i) {
		for (j = 1; j <= 3; ++j) {
		    do_uio(&c__1, (char *)&alparm_1.alparm[j + i * 3 - 4], (
			    ftnlen)sizeof(doublereal));
		}
	    }
	    e_wsue();
	    io___32.ciunit = ir;
	    s_wsue(&io___32);
	    do_uio(&c__1, (char *)&alparm_1.jloop, (ftnlen)sizeof(integer));
	    do_uio(&c__1, (char *)&alparm_1.x0, (ftnlen)sizeof(doublereal));
	    do_uio(&c__1, (char *)&alparm_1.x1, (ftnlen)sizeof(doublereal));
	    do_uio(&c__1, (char *)&alparm_1.x2, (ftnlen)sizeof(doublereal));
	    e_wsue();
	}
	cl__1.cerr = 0;
	cl__1.cunit = 9;
	cl__1.csta = 0;
	f_clos(&cl__1);
	cl__1.cerr = 0;
	cl__1.cunit = 10;
	cl__1.csta = 0;
	f_clos(&cl__1);
	return 0;
    } else {
	s_wsfe(&io___33);
	e_wsfe();
	io___34.ciunit = ir;
	s_rsue(&io___34);
	do_uio(&c__9, (char *)&ipow[1], (ftnlen)sizeof(integer));
	do_uio(&c__1, (char *)&(*iloop), (ftnlen)sizeof(integer));
	e_rsue();
	io___35.ciunit = ir;
	s_rsue(&io___35);
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    do_uio(&c__1, (char *)&xparam[i], (ftnlen)sizeof(doublereal));
	}
	e_rsue();
	io___36.ciunit = ir;
	s_rsue(&io___36);
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    do_uio(&c__1, (char *)&grad[i], (ftnlen)sizeof(doublereal));
	}
	e_rsue();
	io___37.ciunit = ir;
	s_rsue(&io___37);
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    i__2 = geovar_1.nvar;
	    for (j = 1; j <= i__2; ++j) {
		do_uio(&c__1, (char *)&hess[j + i * 258], (ftnlen)sizeof(
			doublereal));
	    }
	}
	e_rsue();
	io___38.ciunit = ir;
	s_rsue(&io___38);
	i__2 = geovar_1.nvar;
	for (i = 1; i <= i__2; ++i) {
	    i__1 = geovar_1.nvar;
	    for (j = 1; j <= i__1; ++j) {
		do_uio(&c__1, (char *)&bmat[j + i * 258], (ftnlen)sizeof(
			doublereal));
	    }
	}
	e_rsue();
	funct1 = sqrt(dot_(&grad[1], &grad[1], &geovar_1.nvar));
	s_wsfe(&io___39);
	do_fio(&c__1, (char *)&funct1, (ftnlen)sizeof(doublereal));
	e_wsfe();
	linear = geovar_1.nvar * (geovar_1.nvar + 1) / 2;
	io___40.ciunit = ir;
	s_rsue(&io___40);
	i__1 = linear;
	for (i = 1; i <= i__1; ++i) {
	    do_uio(&c__1, (char *)&pmat[i], (ftnlen)sizeof(doublereal));
	}
	e_rsue();
	linear = molkst_1.norbs * (molkst_1.norbs + 1) / 2;
	s_rsue(&io___41);
	i__1 = linear;
	for (i = 1; i <= i__1; ++i) {
	    do_uio(&c__1, (char *)&densty_1.pa[i - 1], (ftnlen)sizeof(
		    doublereal));
	}
	e_rsue();
	if (molkst_1.nalpha != 0) {
	    s_rsue(&io___42);
	    i__1 = linear;
	    for (i = 1; i <= i__1; ++i) {
		do_uio(&c__1, (char *)&densty_1.pb[i - 1], (ftnlen)sizeof(
			doublereal));
	    }
	    e_rsue();
	}
	if (path_1.latom != 0) {
	    io___43.ciunit = ir;
	    s_rsue(&io___43);
	    i__1 = geovar_1.nvar;
	    for (i = 1; i <= i__1; ++i) {
		for (j = 1; j <= 3; ++j) {
		    do_uio(&c__1, (char *)&alparm_1.alparm[j + i * 3 - 4], (
			    ftnlen)sizeof(doublereal));
		}
	    }
	    e_rsue();
	    io___44.ciunit = ir;
	    s_rsue(&io___44);
	    do_uio(&c__1, (char *)&alparm_1.jloop, (ftnlen)sizeof(integer));
	    do_uio(&c__1, (char *)&alparm_1.x0, (ftnlen)sizeof(doublereal));
	    do_uio(&c__1, (char *)&alparm_1.x1, (ftnlen)sizeof(doublereal));
	    do_uio(&c__1, (char *)&alparm_1.x2, (ftnlen)sizeof(doublereal));
	    e_rsue();
	    ++(*iloop);
	}
	++(*iloop);
	return 0;
    }
    return 0;
} /* powsav_ */