/* powsq.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 iflepo, iscf;
} mesage_;

#define mesage_1 mesage_

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

#define geovar_1 geovar_

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

#define geom_1 geom_

struct {
    integer last;
} last_;

#define last_1 last_

struct {
    char keywrd[80];
} keywrd_;

#define keywrd_1 keywrd_

struct {
    doublereal time0;
} time_;

#define time_1 time_

struct {
    integer nscf;
} numscf_;

#define numscf_1 numscf_

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

#define geosym_1 geosym_

struct {
    doublereal grad[258], gnfina;
} gradnt_;

#define gradnt_1 gradnt_

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

#define geokst_1 geokst_

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 numcal;
} numcal_;

#define numcal_1 numcal_

struct {
    doublereal gnext, amin, anext;
} sigma1_;

#define sigma1_1 sigma1_

struct {
    doublereal gnext1[258], gmin1[258];
} sigma2_;

#define sigma2_1 sigma2_

struct {
    doublereal hess[66564]	/* was [258][258] */, bmat[66564]	/* 
	    was [258][258] */, pmat[66564];
} nllcom_;

#define nllcom_1 nllcom_

struct {
    doublereal pvec[66564];
} scrach_;

#define scrach_1 scrach_

/* Table of constant values */

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

/* Subroutine */ int powsq_(doublereal *xparam, integer *nvar, doublereal *
	funct)
{
    /* Initialized data */

    static char chdot[1] = ".";
    static char zero[1] = "0";
    static char nine[1] = "9";
    static integer icalcn = 0;

    /* Format strings */
    static char fmt_410[] = "(\002  RESTART FILE WRITTEN,  TIME LEFT:\002,f9"
	    ".1,\002 GRAD.:\002,f10.3,\002 HEAT:\002,g14.7)";
    static char fmt_420[] = "(\002 CYCLE:\002,i5,\002 TIME:\002,f6.1,\002 TI"
	    "ME LEFT:\002,f9.1,\002 GRAD.:\002,f10.3,\002 HEAT:\002,g14.7)";

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

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

    /* Local variables */
    static integer indc, icyc;
    static doublereal gmin, xinc;
    static integer ilpr;
    static doublereal pmax, step;
    static integer ipow[9];
    static doublereal work[258], time1, time2;
    static integer i, j, k;
    extern doublereal reada_(char *, integer *, ftnlen);
    static integer l;
    static doublereal p[258], q[258], alpha;
    static logical debug;
    static doublereal glast;
    extern /* Subroutine */ int hqrii_(doublereal *, integer *, integer *, 
	    doublereal *, doublereal *);
    static doublereal tleft, e1[258], e2[258];
    static logical times, rough;
    static doublereal tlast, tdump;
    static integer iloop;
    static doublereal tstep;
    static char ch[1];
    static integer id, if__, ij, ik, il;
    static doublereal sk;
    extern /* Subroutine */ int search_(doublereal *, doublereal *, 
	    doublereal *, integer *, doublereal *, logical *, logical *, 
	    doublereal *);
    extern doublereal second_(void);
    extern /* Subroutine */ int compfg_(doublereal *, logical *, doublereal *,
	     logical *, doublereal *, logical *);
    static doublereal rt;
    static logical resfil;
    static doublereal totime;
    extern /* Subroutine */ int vecprt_(doublereal *, integer *);
    static integer ip1;
    extern /* Subroutine */ int powsav_(doublereal *, doublereal *, 
	    doublereal *, doublereal *, integer *, doublereal *, integer *);
    static logical restrt;
    static doublereal eig[258];
    static logical okc, okf;
    static doublereal sig[258];
    extern doublereal dot_(doublereal *, doublereal *, integer *);
    static doublereal tim, rmu, sum, rmx;
    static logical scf1;
    static integer nat3;
    static doublereal rho2, tol2;

    /* Fortran I/O blocks */
    static cilist io___19 = { 0, 6, 0, "(//10X,' TIME FOR THIS STEP =',F10.2"
	    ",            ' SECONDS')", 0 };
    static cilist io___30 = { 0, 6, 0, "(/,A)", 0 };
    static cilist io___33 = { 0, 6, 0, "(' XPARAM',6F10.6)", 0 };
    static cilist io___34 = { 0, 6, 0, "(//10X,' RESTARTING AT POINT',I3)", 0 
	    };
    static cilist io___35 = { 0, 6, 0, "(//10X,'RESTARTING IN OPTIMISATION',"
	    "                   ' ROUTINES')", 0 };
    static cilist io___36 = { 0, 6, 0, "(' XPARAM')", 0 };
    static cilist io___37 = { 0, 6, 0, "(5(2I3,F10.4))", 0 };
    static cilist io___38 = { 0, 6, 0, "(' STARTING GRADIENTS')", 0 };
    static cilist io___39 = { 0, 6, 0, "(3X,8F9.4)", 0 };
    static cilist io___43 = { 0, 6, 0, "(I3,12(8F9.4,/3X))", 0 };
    static cilist io___46 = { 0, 6, 0, "(' TIME FOR STEP:',F8.2,' LEFT',F8.2)"
	    , 0 };
    static cilist io___47 = { 0, 6, 0, "(//10X,'UN-NORMALIZED HESSIAN MATRIX"
	    "')", 0 };
    static cilist io___48 = { 0, 6, 0, "(8F10.4)", 0 };
    static cilist io___51 = { 0, 6, 0, "(//10X,'HESSIAN MATRIX')", 0 };
    static cilist io___52 = { 0, 6, 0, "(8F10.4)", 0 };
    static cilist io___57 = { 0, 6, 0, "(/10X,'P MATRIX IN POWSQ')", 0 };
    static cilist io___65 = { 0, 6, 0, "(' SEARCH VECTOR')", 0 };
    static cilist io___66 = { 0, 6, 0, "(8F10.5)", 0 };
    static cilist io___78 = { 0, 6, 0, fmt_410, 0 };
    static cilist io___79 = { 0, 6, 0, fmt_420, 0 };
    static cilist io___80 = { 0, 6, 0, "(' TIME FOR STEP:',F8.2,' LEFT',F8.2)"
	    , 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 */

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

/*   POWSQ OPTIMIZES THE GEOMETRY BY MINIMISING THE GRADIENT NORM. */
/*         THUS BOTH GROUND AND TRANSITION STATE GEOMETRIES CAN BE */
/*         CALCULATED. IT IS ROUGHLY EQUIVALENT TO FLEPO, FLEPO MINIMIZES 
*/
/*         THE ENERGY, POWSQ MINIMIZES THE GRADIENT NORM. */

/*  ON ENTRY XPARAM = VALUES OF PARAMETERS TO BE OPTIMIZED. */
/*           NVAR   = NUMBER OF PARAMETERS TO BE OPTIMIZED. */

/*  ON EXIT  XPARAM = OPTIMIZED PARAMETERS. */
/*           FUNCT  = HEAT OF FORMATION IN KCALS. */

/* ********************************************************************* 
*/
/*        *****  ROUTINE PERFORMS  A LEAST SQUARES MINIMIZATION  ***** */
/*        *****  OF A FUNCTION WHICH IS A SUM OF SQUARES.        ***** */
/*        *****  INITIALLY WRITTEN BY J.W. MCIVER JR. AT SUNY/   ***** */
/*        *****  BUFFALO, SUMMER 1971.  REWRITTEN AND MODIFIED   ***** */
/*        *****  BY A.K. AT SUNY BUFFALO AND THE UNIVERSITY OF   ***** */
/*        *****  TEXAS.  DECEMBER 1973                           ***** */

    /* Parameter adjustments */
    --xparam;

    /* Function Body */
    if (icalcn != numcal_1.numcal) {
	icalcn = numcal_1.numcal;
	restrt = i_indx(keywrd_1.keywrd, "RESTART", 80L, 7L) != 0;
	scf1 = i_indx(keywrd_1.keywrd, "1SCF", 80L, 4L) != 0;
	rough = FALSE_;
	time1 = second_();
	time2 = time1;
	icyc = 0;
	times = i_indx(keywrd_1.keywrd, "TIME", 80L, 4L) != 0;
	totime = 3600.;
	i = i_indx(keywrd_1.keywrd, " T=", 80L, 3L);
	if (i != 0) {
	    tim = reada_(keywrd_1.keywrd, &i, 80L);
	    for (j = i + 3; j <= 80; ++j) {
		i__1 = j;
		s_copy(ch, keywrd_1.keywrd + i__1, 1L, j + 1 - i__1);
		if (*(unsigned char *)ch == ' ') {
		    *(unsigned char *)ch = *(unsigned char *)&keywrd_1.keywrd[
			    j - 1];
		    if (*(unsigned char *)ch == 'M') {
			tim *= 60;
		    }
		    if (*(unsigned char *)ch == 'H') {
			tim *= 3600;
		    }
		    if (*(unsigned char *)ch == 'D') {
			tim *= 86400;
		    }
		    goto L20;
		}
/* L10: */
	    }
L20:
	    totime = tim;
	    s_wsfe(&io___19);
	    do_fio(&c__1, (char *)&totime, (ftnlen)sizeof(doublereal));
	    e_wsfe();
	}
	tleft = totime;
	tlast = totime;
	tdump = 3600.;
	i = i_indx(keywrd_1.keywrd, " DUMP", 80L, 5L);
	if (i != 0) {
	    tdump = reada_(keywrd_1.keywrd, &i, 80L);
	    for (j = i + 7; j <= 80; ++j) {
		*(unsigned char *)ch = *(unsigned char *)&keywrd_1.keywrd[j - 
			1];
		if (*(unsigned char *)ch != *(unsigned char *)&chdot[0] && (*(
			unsigned char *)ch < *(unsigned char *)&zero[0] || *(
			unsigned char *)ch > *(unsigned char *)&nine[0])) {
		    if (*(unsigned char *)ch == 'M') {
			tdump *= 60;
		    }
		    goto L40;
		}
/* L30: */
	    }
L40:
	    ;
	}
	resfil = FALSE_;
	step = .02;
	last_1.last = 0;
	iloop = 1;
	nat3 = molkst_1.numat * 3;
	xinc = .00529167;
	rho2 = 1e-4;
	tol2 = .4;
	if (i_indx(keywrd_1.keywrd, "PREC", 80L, 4L) != 0) {
	    tol2 = .01;
	}
	if (i_indx(keywrd_1.keywrd, "GNORM", 80L, 5L) != 0) {
	    i__1 = i_indx(keywrd_1.keywrd, "GNORM", 80L, 5L);
	    tol2 = reada_(keywrd_1.keywrd, &i__1, 80L);
	    if (tol2 < .001f && i_indx(keywrd_1.keywrd, "LET", 80L, 3L) == 0) 
		    {
		s_wsfe(&io___30);
		do_fio(&c__1, "  GNORM HAS BEEN SET TOO LOW, RESET TO 0.001", 
			44L);
		e_wsfe();
		tol2 = .001;
	    }
	}
	debug = i_indx(keywrd_1.keywrd, "POWSQ", 80L, 5L) != 0;
	if (restrt) {

/*   RESTORE STORED DATA */

	    ipow[8] = 0;
	    powsav_(nllcom_1.hess, sigma2_1.gmin1, &xparam[1], nllcom_1.pmat, 
		    &iloop, nllcom_1.bmat, ipow);
	    numscf_1.nscf = ipow[7];
	    i__1 = *nvar;
	    for (i = 1; i <= i__1; ++i) {
		gradnt_1.grad[i - 1] = sigma2_1.gmin1[i - 1];
/* L50: */
		sigma2_1.gnext1[i - 1] = sigma2_1.gmin1[i - 1];
	    }
	    s_wsfe(&io___33);
	    i__1 = *nvar;
	    for (i = 1; i <= i__1; ++i) {
		do_fio(&c__1, (char *)&xparam[i], (ftnlen)sizeof(doublereal));
	    }
	    e_wsfe();
	    if (iloop > 0) {
/* #               ILOOP=ILOOP+1 */
		s_wsfe(&io___34);
		do_fio(&c__1, (char *)&iloop, (ftnlen)sizeof(integer));
		e_wsfe();
	    } else {
		s_wsfe(&io___35);
		e_wsfe();
	    }
	}

/*   DEFINITIONS:   NVAR   = NUMBER OF GEOMETRIC VARIABLES = 3*NUMAT-6
 */

    }
    *nvar = abs(*nvar);
    if (debug) {
	s_wsfe(&io___36);
	e_wsfe();
	s_wsfe(&io___37);
	i__1 = *nvar;
	for (i = 1; i <= i__1; ++i) {
	    do_fio(&c__1, (char *)&geovar_1.loc[(i << 1) - 2], (ftnlen)sizeof(
		    integer));
	    do_fio(&c__1, (char *)&geovar_1.loc[(i << 1) - 1], (ftnlen)sizeof(
		    integer));
	    do_fio(&c__1, (char *)&xparam[i], (ftnlen)sizeof(doublereal));
	}
	e_wsfe();
    }
    if (! restrt) {
	i__1 = *nvar;
	for (i = 1; i <= i__1; ++i) {
/* L60: */
	    gradnt_1.grad[i - 1] = 0.;
	}
	compfg_(&xparam[1], &c_true, funct, &c_true, gradnt_1.grad, &c_true);
    }
    if (debug) {
	s_wsfe(&io___38);
	e_wsfe();
	s_wsfe(&io___39);
	i__1 = *nvar;
	for (i = 1; i <= i__1; ++i) {
	    do_fio(&c__1, (char *)&gradnt_1.grad[i - 1], (ftnlen)sizeof(
		    doublereal));
	}
	e_wsfe();
    }
    gmin = sqrt(dot_(gradnt_1.grad, gradnt_1.grad, nvar));
    glast = gmin;
    i__1 = *nvar;
    for (i = 1; i <= i__1; ++i) {
	sigma2_1.gnext1[i - 1] = gradnt_1.grad[i - 1];
	sigma2_1.gmin1[i - 1] = sigma2_1.gnext1[i - 1];
/* L70: */
    }

/*    NOW TO CALCULATE THE HESSIAN MATRIX. */

    if (iloop < 0) {
	goto L180;
    }

/*   CHECK THAT HESSIAN HAS NOT ALREADY BEEN CALCULATED. */

    ilpr = iloop;
    i__1 = *nvar;
    for (iloop = ilpr; iloop <= i__1; ++iloop) {
	time1 = second_();
	xparam[iloop] += xinc;
	compfg_(&xparam[1], &c_true, funct, &c_true, gradnt_1.grad, &c_true);
	if (scf1) {
	    goto L430;
	}
	if (debug) {
	    s_wsfe(&io___43);
	    do_fio(&c__1, (char *)&iloop, (ftnlen)sizeof(integer));
	    i__2 = *nvar;
	    for (if__ = 1; if__ <= i__2; ++if__) {
		do_fio(&c__1, (char *)&gradnt_1.grad[if__ - 1], (ftnlen)
			sizeof(doublereal));
	    }
	    e_wsfe();
	}
	gradnt_1.grad[iloop - 1] += 1e-5;
	xparam[iloop] -= xinc;
	i__2 = *nvar;
	for (j = 1; j <= i__2; ++j) {
/* L80: */
	    nllcom_1.hess[iloop + j * 258 - 259] = -(gradnt_1.grad[j - 1] - 
		    sigma2_1.gnext1[j - 1]) / xinc;
	}
	time2 = second_();
	tstep = time2 - time1;
	if (times) {
	    s_wsfe(&io___46);
	    do_fio(&c__1, (char *)&tstep, (ftnlen)sizeof(doublereal));
	    do_fio(&c__1, (char *)&tleft, (ftnlen)sizeof(doublereal));
	    e_wsfe();
	}
	if (tlast - tleft > tdump) {
	    tlast = tleft;
	    resfil = TRUE_;
	    ipow[8] = 2;
	    i = iloop;
	    ipow[7] = numscf_1.nscf;
	    powsav_(nllcom_1.hess, sigma2_1.gmin1, &xparam[1], nllcom_1.pmat, 
		    &i, nllcom_1.bmat, ipow);
	}
	if (tleft < tstep * 2.) {

/*  STORE RESULTS TO DATE. */

	    ipow[8] = 1;
	    i = iloop;
	    ipow[7] = numscf_1.nscf;
	    powsav_(nllcom_1.hess, sigma2_1.gmin1, &xparam[1], nllcom_1.pmat, 
		    &i, nllcom_1.bmat, ipow);
	    s_stop("", 0L);
	}
/* L90: */
    }
/*        *****  SCALE -HESSIAN- MATRIX                           ***** */
    if (debug) {
	s_wsfe(&io___47);
	e_wsfe();
	i__1 = *nvar;
	for (i = 1; i <= i__1; ++i) {
/* L100: */
	    s_wsfe(&io___48);
	    i__2 = *nvar;
	    for (j = 1; j <= i__2; ++j) {
		do_fio(&c__1, (char *)&nllcom_1.hess[j + i * 258 - 259], (
			ftnlen)sizeof(doublereal));
	    }
	    e_wsfe();
	}
    }
    i__2 = *nvar;
    for (i = 1; i <= i__2; ++i) {
	sum = 0.;
	i__1 = *nvar;
	for (j = 1; j <= i__1; ++j) {
/* L110: */
/* Computing 2nd power */
	    d__1 = nllcom_1.hess[i + j * 258 - 259];
	    sum += d__1 * d__1;
	}
/* L120: */
	work[i - 1] = 1. / sqrt(sum);
    }
    i__2 = *nvar;
    for (i = 1; i <= i__2; ++i) {
	i__1 = *nvar;
	for (j = 1; j <= i__1; ++j) {
/* L130: */
	    nllcom_1.hess[i + j * 258 - 259] *= work[i - 1];
	}
    }
    if (debug) {
	s_wsfe(&io___51);
	e_wsfe();
	i__1 = *nvar;
	for (i = 1; i <= i__1; ++i) {
/* L140: */
	    s_wsfe(&io___52);
	    i__2 = *nvar;
	    for (j = 1; j <= i__2; ++j) {
		do_fio(&c__1, (char *)&nllcom_1.hess[j + i * 258 - 259], (
			ftnlen)sizeof(doublereal));
	    }
	    e_wsfe();
	}
    }
/*        *****  INITIALIZE B MATIRX                        ***** */
    i__2 = *nvar;
    for (i = 1; i <= i__2; ++i) {
	i__1 = *nvar;
	for (j = 1; j <= i__1; ++j) {
/* L150: */
	    nllcom_1.bmat[i + j * 258 - 259] = 0.;
	}
/* L160: */
	nllcom_1.bmat[i + i * 258 - 259] = work[i - 1] * 2.;
    }
/* ***********************************************************************
 */

/*  THIS IS THE START OF THE BIG LOOP TO OPTIMIZE THE GEOMETRY */

/* ***********************************************************************
 */
    iloop = -99;
    tstep *= 4;
L170:
    if (tlast - tleft > tdump) {
	tlast = tleft;
	resfil = TRUE_;
	ipow[8] = 2;
	i = iloop;
	ipow[7] = numscf_1.nscf;
	powsav_(nllcom_1.hess, sigma2_1.gmin1, &xparam[1], nllcom_1.pmat, &i, 
		nllcom_1.bmat, ipow);
    }
    if (tleft < tstep * 2.) {

/*  STORE RESULTS TO DATE. */

	ipow[8] = 1;
	i = iloop;
	ipow[7] = numscf_1.nscf;
	powsav_(nllcom_1.hess, sigma2_1.gmin1, &xparam[1], nllcom_1.pmat, &i, 
		nllcom_1.bmat, ipow);
	s_stop("", 0L);
    }
L180:
/*        *****  FORM-A- DAGGER-A- IN PA SLONG WITH -P-     ***** */
    ij = 0;
    i__2 = *nvar;
    for (j = 1; j <= i__2; ++j) {
	i__1 = j;
	for (i = 1; i <= i__1; ++i) {
	    ++ij;
	    sum = 0.;
	    i__3 = *nvar;
	    for (k = 1; k <= i__3; ++k) {
/* L190: */
		sum += nllcom_1.hess[i + k * 258 - 259] * nllcom_1.hess[j + k 
			* 258 - 259];
	    }
/* L200: */
	    nllcom_1.pmat[ij - 1] = sum;
	}
    }
    i__1 = *nvar;
    for (i = 1; i <= i__1; ++i) {
	sum = 0.;
	i__2 = *nvar;
	for (k = 1; k <= i__2; ++k) {
/* L210: */
	    sum -= nllcom_1.hess[i + k * 258 - 259] * sigma2_1.gmin1[k - 1];
	}
/* L220: */
	p[i - 1] = -sum;
    }
    l = 0;
    if (debug) {
	s_wsfe(&io___57);
	e_wsfe();
	vecprt_(nllcom_1.pmat, nvar);
    }
    hqrii_(nllcom_1.pmat, nvar, nvar, eig, scrach_1.pvec);
/*        *****  CHECK FOR ZERO EIGENVALUE                  ***** */
/* #      WRITE(6,'(''  EIGS IN POWSQ:'')') */
/* #      WRITE(6,'(6F13.8)')(EIG(I),I=1,NVAR) */
    if (eig[0] < rho2) {
	goto L280;
    }
    indc = 2;
/*        *****  IF MATRIX IS NOT SINGULAR FORM INVERSE     ***** */
/*        *****  BY BACK TRANSFORMING THE EIGENVECTORS      ***** */
    ij = 0;
    i__1 = *nvar;
    for (i = 1; i <= i__1; ++i) {
	i__2 = i;
	for (j = 1; j <= i__2; ++j) {
	    ++ij;
	    sum = 0.;
	    i__3 = *nvar;
	    for (k = 1; k <= i__3; ++k) {
/* L230: */
		sum += scrach_1.pvec[(k - 1) * *nvar + j - 1] * scrach_1.pvec[
			(k - 1) * *nvar + i - 1] / eig[k - 1];
	    }
/* L240: */
	    nllcom_1.pmat[ij - 1] = sum;
	}
    }
/*        *****  FIND -Q- VECTOR                            ***** */
    l = 0;
    il = l + 1;
    l = il + i - 1;
    i__2 = *nvar;
    for (i = 1; i <= i__2; ++i) {
	sum = 0.;
	i__1 = i;
	for (k = 1; k <= i__1; ++k) {
	    ik = i * (i - 1) / 2 + k;
/* L250: */
	    sum += nllcom_1.pmat[ik - 1] * p[k - 1];
	}
	ip1 = i + 1;
	i__1 = *nvar;
	for (k = ip1; k <= i__1; ++k) {
	    ik = k * (k - 1) / 2 + i;
/* L260: */
	    sum += nllcom_1.pmat[ik - 1] * p[k - 1];
	}
/* L270: */
	q[i - 1] = sum;
    }
    goto L300;
L280:
/*        *****  TAKE  -Q- VECTOR AS EIGENVECTOR OF ZERO     ***** */
/*        *****  EIGENVALUE                                 ***** */
    i__2 = *nvar;
    for (i = 1; i <= i__2; ++i) {
/* L290: */
	q[i - 1] = scrach_1.pvec[i - 1];
    }
L300:
/*        *****  FIND SEARCH DIRECTION                      ***** */
    i__2 = *nvar;
    for (i = 1; i <= i__2; ++i) {
	sig[i - 1] = 0.;
	i__1 = *nvar;
	for (j = 1; j <= i__1; ++j) {
/* L310: */
	    sig[i - 1] += q[j - 1] * nllcom_1.bmat[i + j * 258 - 259];
	}
    }
/*        *****  DO A ONE DIMENSIONAL SEARCH                ***** */
    if (debug) {
	s_wsfe(&io___65);
	e_wsfe();
	s_wsfe(&io___66);
	i__1 = *nvar;
	for (i = 1; i <= i__1; ++i) {
	    do_fio(&c__1, (char *)&sig[i - 1], (ftnlen)sizeof(doublereal));
	}
	e_wsfe();
    }
    search_(&xparam[1], &alpha, sig, nvar, &gmin, &okc, &okf, funct);
    if (*nvar == 1) {
	goto L430;
    }

/*  FIRST WE ATTEMPT TO OPTIMIZE GEOMETRY USING SEARCH. */
/*  IF THIS DOES NOT WORK, THEN SWITCH TO LINMIN, WHICH ALWAYS WORKS, */
/*  BUT IS TWICE AS SLOW AS SEARCH. */

    rough = ! okf;
    rmx = 0.;
    i__1 = *nvar;
    for (k = 1; k <= i__1; ++k) {
	rt = (d__1 = sigma2_1.gmin1[k - 1], abs(d__1));
	if (rt > rmx) {
	    rmx = rt;
	}
/* L320: */
    }
    if (rmx < tol2) {
	goto L430;
    }
/*        *****  TWO STEP ESTIMATION OF DERIVATIVES         ***** */
    i__1 = *nvar;
    for (k = 1; k <= i__1; ++k) {
/* L330: */
	e1[k - 1] = (sigma2_1.gmin1[k - 1] - sigma2_1.gnext1[k - 1]) / (
		sigma1_1.amin - sigma1_1.anext);
    }
    rmu = dot_(e1, sigma2_1.gmin1, nvar) / dot_(sigma2_1.gmin1, 
	    sigma2_1.gmin1, nvar);
    i__1 = *nvar;
    for (k = 1; k <= i__1; ++k) {
/* L340: */
	e2[k - 1] = e1[k - 1] - rmu * sigma2_1.gmin1[k - 1];
    }
/*        *****  SCALE -E2- AND -SIG-                       ***** */
    sk = 1. / sqrt(dot_(e2, e2, nvar));
    i__1 = *nvar;
    for (k = 1; k <= i__1; ++k) {
/* L350: */
	sig[k - 1] = sk * sig[k - 1];
    }
    i__1 = *nvar;
    for (k = 1; k <= i__1; ++k) {
/* L360: */
	e2[k - 1] = sk * e2[k - 1];
    }
/*        *****  FIND INDEX OF REPLACEMENT DIRECTION        ***** */
    pmax = -1e20;
    i__1 = *nvar;
    for (i = 1; i <= i__1; ++i) {
	if ((d__1 = p[i - 1] * q[i - 1], abs(d__1)) <= pmax) {
	    goto L370;
	}
	pmax = (d__1 = p[i - 1] * q[i - 1], abs(d__1));
	id = i;
L370:
	;
    }
/*        *****  REPLACE APPROPRIATE DIRECTION AND DERIVATIVE *** */
    i__1 = *nvar;
    for (k = 1; k <= i__1; ++k) {
/* L380: */
	nllcom_1.hess[id + k * 258 - 259] = -e2[k - 1];
    }
/*        *****  REPLACE STARTING POINT                     ***** */
    i__1 = *nvar;
    for (k = 1; k <= i__1; ++k) {
/* L390: */
	nllcom_1.bmat[k + id * 258 - 259] = sig[k - 1] / .529167;
    }
    i__1 = *nvar;
    for (k = 1; k <= i__1; ++k) {
/* L400: */
	sigma2_1.gnext1[k - 1] = sigma2_1.gmin1[k - 1];
    }
    glast = gmin;
    indc = 1;
    time1 = time2;
    time2 = second_();
    tleft = totime - time2 + time_1.time0;
    tstep = time2 - time1;
    ++icyc;
    if (resfil) {
	s_wsfe(&io___78);
	do_fio(&c__1, (char *)&tleft, (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&gmin, (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&(*funct), (ftnlen)sizeof(doublereal));
	e_wsfe();
	resfil = FALSE_;
    } else {
	s_wsfe(&io___79);
	do_fio(&c__1, (char *)&icyc, (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&tstep, (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&tleft, (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&gmin, (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&(*funct), (ftnlen)sizeof(doublereal));
	e_wsfe();
    }
    if (times) {
	s_wsfe(&io___80);
	do_fio(&c__1, (char *)&tstep, (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&tleft, (ftnlen)sizeof(doublereal));
	e_wsfe();
    }
    goto L170;
L430:
    i__1 = *nvar;
    for (i = 1; i <= i__1; ++i) {
/* L440: */
	gradnt_1.grad[i - 1] = 0.;
    }
    last_1.last = 1;
    compfg_(&xparam[1], &c_true, funct, &c_true, gradnt_1.grad, &c_true);
    i__1 = *nvar;
    for (i = 1; i <= i__1; ++i) {
/* L450: */
	gradnt_1.grad[i - 1] = sigma2_1.gmin1[i - 1];
    }
    gradnt_1.gnfina = sqrt(dot_(gradnt_1.grad, gradnt_1.grad, nvar));
    mesage_1.iflepo = 11;
    return 0;
} /* powsq_ */