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

#define keywrd_1 keywrd_

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

#define geokst_1 geokst_

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

#define geovar_1 geovar_

struct {
    doublereal p[23220], pdumy[46440]	/* was [2][23220] */;
} densty_;

#define densty_1 densty_

struct {
    doublereal atmass[86];
} atmass_;

#define atmass_1 atmass_

struct {
    doublereal time0;
} time_;

#define time_1 time_

struct {
    doublereal core[107];
} core_;

#define core_1 core_

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 coord[258]	/* was [3][86] */;
} coord_;

#define coord_1 coord_

struct {
    doublereal evecs[66564];
} scrach_;

#define scrach_1 scrach_

/* Table of constant values */

static integer c__1 = 1;
static logical c_true = TRUE_;
static integer c__0 = 0;
static logical c_false = FALSE_;

/* Subroutine */ int fmat_(doublereal *fmatrx, integer *nreal, doublereal *
	tscf, doublereal *tder, doublereal *deldip, doublereal *heat)
{
    /* Initialized data */

    static doublereal fact = .00695125;

    /* System generated locals */
    integer i__1, i__2, i__3;
    doublereal d__1, d__2, d__3, d__4, d__5;

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

    /* Local variables */
    static doublereal grad[258], escf, cold[258], eigs[258], time;
    static integer ifor;
    static logical prnt;
    static doublereal g2rad[258], g2old[258], time1, time2, time3;
    static integer i, j, k;
    extern doublereal reada_(char *, integer *, ftnlen);
    static integer l;
    static doublereal q[86], heata, heatb;
    static logical debug;
    extern /* Subroutine */ int chrge_(doublereal *, doublereal *);
    static doublereal delta;
    extern /* Subroutine */ int frame_(doublereal *, integer *, integer *, 
	    doublereal *);
    static doublereal grold[258];
    static logical deriv;
    static doublereal tleft, estim;
    static integer iloop;
    static doublereal tlast, tdump;
    static integer k1;
    static doublereal tstep;
    static integer k2;
    static char ch[1];
    static integer ii;
    static doublereal heataa, heatbb;
    static integer kk, ll, lu, ix;
    extern doublereal second_(void);
    extern /* Subroutine */ int compfg_(doublereal *, logical *, doublereal *,
	     logical *, doublereal *, logical *);
    extern doublereal dipole_(doublereal *, doublereal *, doublereal *, 
	    doublereal *);
#define coordl ((doublereal *)&coord_1)
    static logical resfil, precis;
    static doublereal fconst[258];
    static logical analyt;
    extern /* Subroutine */ int forsav_(doublereal *, doublereal *, integer *,
	     integer *, doublereal *, doublereal *, integer *, doublereal *, 
	    doublereal *, integer *, doublereal *);
    static doublereal totime;
    static integer istart, jstart, kountf;
    static doublereal estime;
    extern /* Subroutine */ int matout_(doublereal *, doublereal *, integer *,
	     integer *, integer *);
    static logical restrt;
    static integer lin;
    extern doublereal dot_(doublereal *, doublereal *, integer *);
    static doublereal tim, sum;
    extern /* Subroutine */ int rsp_(doublereal *, integer *, integer *, 
	    doublereal *, doublereal *);
    static doublereal del2[3];

    /* Fortran I/O blocks */
    static cilist io___16 = { 0, 6, 0, "(//4X,'FIRST DERIVATIVES WILL BE USE"
	    "D IN THE'  ,' CALCULATION OF SECOND DERIVATIVES')", 0 };
    static cilist io___20 = { 0, 6, 0, "(/10X,'TIME DEFINED FOR THIS STEP ='"
	    ",F19.2,     ' SECONDS')", 0 };
    static cilist io___21 = { 0, 6, 0, "(/10X,'DEFAULT TIME OF',F8.2,       "
	    "            ' SECONDS ALLOCATED FOR THIS STEP')", 0 };
    static cilist io___34 = { 0, 6, 0, "(/10X,'ESTIMATED TIME TO COMPLETE CA"
	    "LCULATION ='       ,F9.2,' SECONDS')", 0 };
    static cilist io___35 = { 0, 6, 0, "(/10X,'STARTING AGAIN AT LINE',18X,I"
	    "4)", 0 };
    static cilist io___36 = { 0, 6, 0, "(/10X,'TIME USED UP TO RESTART =',F2"
	    "2.2)", 0 };
    static cilist io___56 = { 0, 6, 0, "(' STEP:',I4,' RESTART FILE WRITTEN,"
	    " INTEGRAL = ',F10.2,' TIME LEFT:',F10.2)", 0 };
    static cilist io___57 = { 0, 6, 0, "(' STEP:',I4,' TIME =',F9.2,' SECS, "
	    "INTEGRAL =',F10.2,' TIME LEFT:',F10.2)", 0 };
    static cilist io___59 = { 0, 6, 0, "(//10X,'- - - - - - - TIME UP - - - "
	    "- - - -',//)", 0 };
    static cilist io___60 = { 0, 6, 0, "(/10X,' POINT REACHED =',I4)", 0 };
    static cilist io___61 = { 0, 6, 0, "(/10X,' RESTART USING KEY-WORD \"RES"
	    "TART\"')", 0 };
    static cilist io___62 = { 0, 6, 0, "(10X,'ESTIMATED TIME FOR THE NEXT ST"
	    "EP =',F8.2,  ' SECONDS')", 0 };
    static cilist io___63 = { 0, 6, 0, "(//10X,'FORCE MATRIX WRITTEN TO DISK"
	    "')", 0 };
    static cilist io___64 = { 0, 6, 0, "(//10X,' STARTING TO CALCULATE FORCE"
	    " CONSTANTS',/)", 0 };
    static cilist io___66 = { 0, 6, 0, "('   EIGENVECTORS FROM FIRST CALCULA"
	    "TION')", 0 };
    static cilist io___73 = { 0, 6, 0, "(' STEP:',I4,' TIME =',F9.2,' SECS, "
	    "INTEGRAL =',F10.2,' TIME LEFT:',F10.2)", 0 };
    static cilist io___76 = { 0, 6, 0, "(//10X,'- - - - -  TIME  LIMIT - - -"
	    " - -')", 0 };
    static cilist io___77 = { 0, 6, 0, "(/10X,' POINT REACHED =',I4)", 0 };
    static cilist io___78 = { 0, 6, 0, "(/10X,' RESTART USING KEY-WORD \"RES"
	    "TART\"')", 0 };
    static cilist io___79 = { 0, 6, 0, "(10X,'ESTIMATED TIME FOR THE NEXT ST"
	    "EP =',F8.2,  ' SECONDS')", 0 };
    static cilist io___80 = { 0, 6, 0, "(A)", 0 };
    static cilist io___81 = { 0, 6, 0, "(A)", 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 */

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

/*  VALUE CALCULATES THE SECOND-ORDER OF THE ENERGY WITH */
/*        RESPECT TO THE CARTESIAN COORDINATES I AND J AND PLACES IT */
/*        IN FMATRX */

/*  ON INPUT NATOMS  = NUMBER OF ATOMS IN THE SYSTEM. */
/*           XPARAM  = INTERNAL COORDINATES OF MOLECULE STORED LINEARLY */

/*  VARIABLES USED */
/*           COORDL  = ARRAY OF CARTESIAN COORDINATES, STORED LINEARLY. */
/*           I       = INDEX OF CARTESIAN COORDINATE. */
/*           J       = INDEX OF CARTESIAN COORDINATE. */

/*  ON OUTPUT FMATRX = SECOND DERIVATIVE OF THE ENERGY WITH RESPECT TO */
/*                    CARTESIAN COORDINATES I AND J. */
/* ********************************************************************** 
*/
    /* Parameter adjustments */
    deldip -= 4;
    --fmatrx;

    /* Function Body */

/*    FACT IS THE CONVERSION FACTOR FROM KCAL/MOLE TO ERGS */

/* SET UP CONSTANTS AND FLAGS */
    geokst_1.na[0] = 99;

/*  SET UP THE VARIABLES IN XPARAM ANDLOC,THESE ARE IN CARTESIAN COORDINA 
*/

    molkst_1.numat = 0;
    i__1 = geokst_1.natoms;
    for (i = 1; i <= i__1; ++i) {
	if (geokst_1.labels[i - 1] != 99 && geokst_1.labels[i - 1] != 107) {
	    ++molkst_1.numat;
	    geokst_1.labels[molkst_1.numat - 1] = geokst_1.labels[i - 1];
	}
/* L10: */
    }
    geokst_1.natoms = molkst_1.numat;

/*   THIS IS A QUICK, IF CLUMSY, WAY TO CALCULATE NUMAT, AND TO REMOVE TH 
*/
/*   DUMMY ATOMS FROM THE ARRAY LABELS. */

    geovar_1.nvar = 0;
    i__1 = molkst_1.numat;
    for (i = 1; i <= i__1; ++i) {
	for (j = 1; j <= 3; ++j) {
	    ++geovar_1.nvar;
	    geovar_1.loc[(geovar_1.nvar << 1) - 2] = i;
	    geovar_1.loc[(geovar_1.nvar << 1) - 1] = j;
/* L20: */
	}
    }
    lin = geovar_1.nvar * (geovar_1.nvar + 1) / 2;
    i__1 = lin;
    for (i = 1; i <= i__1; ++i) {
/* L30: */
	fmatrx[i] = 0.;
    }
    prnt = i_indx(keywrd_1.keywrd, "IRC=", 80L, 4L) == 0;
    precis = i_indx(keywrd_1.keywrd, "PRECIS", 80L, 6L) != 0;
    analyt = i_indx(keywrd_1.keywrd, "ANALYT", 80L, 6L) != 0;
    restrt = i_indx(keywrd_1.keywrd, "RESTART", 80L, 7L) != 0;
    if (i_indx(keywrd_1.keywrd, "NLLSQ", 80L, 5L) != 0) {
	restrt = FALSE_;
    }
    debug = i_indx(keywrd_1.keywrd, "FMAT", 80L, 4L) != 0;
    deriv = molkst_1.nclose == molkst_1.nopen && i_indx(keywrd_1.keywrd, 
	    "C.I.", 80L, 4L) == 0;
    if (prnt) {
	s_wsfe(&io___16);
	e_wsfe();
    }
    time = 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;
	    if (s_cmp(keywrd_1.keywrd + i__1, " ", j + 1 - i__1, 1L) == 0) {
		*(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 L50;
	    }
/* L40: */
	}
L50:
	time = tim;
	if (prnt) {
	    s_wsfe(&io___20);
	    do_fio(&c__1, (char *)&time, (ftnlen)sizeof(doublereal));
	    e_wsfe();
	}
    } else {
	if (prnt) {
	    s_wsfe(&io___21);
	    do_fio(&c__1, (char *)&time, (ftnlen)sizeof(doublereal));
	    e_wsfe();
	}
    }
    tleft = time;
    tlast = time;
    tdump = 3600.;
    i = i_indx(keywrd_1.keywrd, " DUMP", 80L, 5L);
    if (i != 0) {
	tdump = reada_(keywrd_1.keywrd, &i, 80L);
	for (j = i + 6; j <= 80; ++j) {
	    i__1 = j;
	    if (s_cmp(keywrd_1.keywrd + i__1, " ", j + 1 - i__1, 1L) == 0) {
		*(unsigned char *)ch = *(unsigned char *)&keywrd_1.keywrd[j - 
			1];
		if (*(unsigned char *)ch == 'M') {
		    tdump *= 60;
		}
		if (*(unsigned char *)ch == 'H') {
		    tdump *= 3600;
		}
		if (*(unsigned char *)ch == 'D') {
		    tdump *= 86400;
		}
		goto L70;
	    }
/* L60: */
	}
L70:
	;
    }
    resfil = FALSE_;
    if (restrt) {
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
/* L80: */
	    cold[i - 1] = coordl[i - 1];
	}
	istart = 0;
	i = 0;
	forsav_(&totime, &deldip[4], &istart, &i, &fmatrx[1], coord_1.coord, &
		geovar_1.nvar, heat, scrach_1.evecs, &jstart, fconst);
	kountf = istart * (istart + 1) / 2;
	++istart;
	++jstart;
	time2 = second_();
	if (istart > geovar_1.nvar) {
	    goto L200;
	}
    } else {
	totime = 0.;
	if (*tscf > 0.) {
	    tleft = tleft - *tscf - *tder;
	}
	istart = 1;
    }
/* CALCULATE FMATRX */
    if (istart > 1) {
	estime = (geovar_1.nvar - istart + 1) * totime / (istart - 1.);
    } else {
	estime = geovar_1.nvar * (*tscf + *tder) * 2.;
	if (precis) {
	    estime *= 2.;
	}
    }
    if (*tscf > 0.) {
	s_wsfe(&io___34);
	do_fio(&c__1, (char *)&estime, (ftnlen)sizeof(doublereal));
	e_wsfe();
    }
    if (restrt) {
	if (istart <= geovar_1.nvar) {
	    s_wsfe(&io___35);
	    do_fio(&c__1, (char *)&istart, (ftnlen)sizeof(integer));
	    e_wsfe();
	}
	s_wsfe(&io___36);
	do_fio(&c__1, (char *)&totime, (ftnlen)sizeof(doublereal));
	e_wsfe();
    }
    lu = kountf;
    time1 = second_();
    molkst_1.numat = geovar_1.nvar / 3;
    i__1 = geovar_1.nvar;
    for (i = istart; i <= i__1; ++i) {
	time2 = second_();
	delta = .016666666666666665;
	if (precis) {

/*   DETERMINE A GOOD STEP SIZE */

	    g2old[0] = 100.;
	    coordl[i - 1] += delta;
	    compfg_(coordl, &c_true, &escf, &c_true, g2old, &c_true);
	    coordl[i - 1] -= delta;
	    delta = delta * 10. / sqrt(dot_(g2old, g2old, &geovar_1.nvar));
/* #         WRITE(6,'(A,F12.5)')' DELTA :',DELTA */
	    g2old[0] = 100.;
	    coordl[i - 1] += delta;
	    compfg_(coordl, &c_true, &escf, &c_true, g2old, &c_true);
/* #         WRITE(6,*)' GNORM:',SQRT(DOT(G2OLD,G2OLD,NVAR)) */
	    coordl[i - 1] -= delta * 2.;
	    g2rad[0] = 100.;
	    compfg_(coordl, &c_true, &heataa, &c_true, g2rad, &c_true);
	    coordl[i - 1] += delta;
	}
	coordl[i - 1] += delta * .5;
	grold[0] = 100.;
	compfg_(coordl, &c_true, &escf, &c_true, grold, &c_true);
/* #         WRITE(6,*)' GNORM:',SQRT(DOT(GROLD,GROLD,NVAR)) */
	chrge_(densty_1.p, q);
	i__2 = molkst_1.numat;
	for (ii = 1; ii <= i__2; ++ii) {
/* L90: */
	    q[ii - 1] = core_1.core[geokst_1.labels[ii - 1] - 1] - q[ii - 1];
	}
	sum = dipole_(densty_1.p, q, coordl, &deldip[i * 3 + 1]);
	coordl[i - 1] -= delta;
	grad[0] = 100.;
	compfg_(coordl, &c_true, &heataa, &c_true, grad, &c_true);
	coordl[i - 1] += delta * .5;
	chrge_(densty_1.p, q);
	i__2 = molkst_1.numat;
	for (ii = 1; ii <= i__2; ++ii) {
/* L100: */
	    q[ii - 1] = core_1.core[geokst_1.labels[ii - 1] - 1] - q[ii - 1];
	}
	sum = dipole_(densty_1.p, q, coordl, del2);
	for (ii = 1; ii <= 3; ++ii) {
/* L110: */
	    deldip[ii + i * 3] = (deldip[ii + i * 3] - del2[ii - 1]) * .5 / 
		    delta;
	}
	ll = lu + 1;
	lu = ll + i - 1;
	l = 0;
	if (precis) {
	    i__2 = lu;
	    for (kountf = ll; kountf <= i__2; ++kountf) {
		++l;
		fmatrx[kountf] += ((grold[l - 1] - grad[l - 1]) * 8. - (g2old[
			l - 1] - g2rad[l - 1])) * .25 / delta * fact / 6.;
/* L120: */
	    }
	    --l;
	    i__2 = geovar_1.nvar;
	    for (k = i; k <= i__2; ++k) {
		++l;
		kk = k * (k - 1) / 2 + i;
		fmatrx[kk] += ((grold[l - 1] - grad[l - 1]) * 8. - (g2old[l - 
			1] - g2rad[l - 1])) * .25 / delta * fact / 6.;
/* L130: */
	    }
	} else {
	    i__2 = lu;
	    for (kountf = ll; kountf <= i__2; ++kountf) {
		++l;
		fmatrx[kountf] += (grold[l - 1] - grad[l - 1]) * .25 / delta *
			 fact;
/* L140: */
	    }
	    --l;
	    i__2 = geovar_1.nvar;
	    for (k = i; k <= i__2; ++k) {
		++l;
		kk = k * (k - 1) / 2 + i;
		fmatrx[kk] += (grold[l - 1] - grad[l - 1]) * .25 / delta * 
			fact;
/* L150: */
	    }
	}
	time3 = second_();
	tstep = time3 - time2;
	totime += tstep;
	tleft -= tstep;
	if (resfil) {
	    s_wsfe(&io___56);
	    do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&totime, (ftnlen)sizeof(doublereal));
	    do_fio(&c__1, (char *)&tleft, (ftnlen)sizeof(doublereal));
	    e_wsfe();
	    resfil = FALSE_;
	} else {
	    s_wsfe(&io___57);
	    do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&tstep, (ftnlen)sizeof(doublereal));
	    do_fio(&c__1, (char *)&totime, (ftnlen)sizeof(doublereal));
	    do_fio(&c__1, (char *)&tleft, (ftnlen)sizeof(doublereal));
	    e_wsfe();
	}
	if (deriv) {
	    estim = totime / i;
	} else {
	    estim = totime * 2. / i;
	}
	if (tlast - tleft > tdump) {
	    tlast = tleft;
	    resfil = TRUE_;
	    jstart = 1;
	    ii = i;
	    forsav_(&totime, &deldip[4], &ii, &geovar_1.nvar, &fmatrx[1], 
		    coord_1.coord, &geovar_1.nvar, heat, scrach_1.evecs, &
		    jstart, fconst);
	}
	if (i != geovar_1.nvar && tleft - 10. < estim) {
	    s_wsfe(&io___59);
	    e_wsfe();
	    s_wsfe(&io___60);
	    do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
	    e_wsfe();
	    s_wsfe(&io___61);
	    e_wsfe();
	    s_wsfe(&io___62);
	    do_fio(&c__1, (char *)&estim, (ftnlen)sizeof(doublereal));
	    e_wsfe();
	    jstart = 1;
	    ii = i;
	    forsav_(&totime, &deldip[4], &ii, &geovar_1.nvar, &fmatrx[1], 
		    coord_1.coord, &geovar_1.nvar, heat, scrach_1.evecs, &
		    jstart, fconst);
	    s_wsfe(&io___63);
	    e_wsfe();
	    s_stop("", 0L);
	}
/* L160: */
    }
/* #      CALL FORSAV(TOTIME,DELDIP,NVAR,NVAR,FMATRX, COORD,NVAR,HEAT, */
/* #     +                EVECS,JSTART,FCONST) */
    if (deriv) {
	goto L290;
    }
    s_wsfe(&io___64);
    e_wsfe();
    frame_(&fmatrx[1], &molkst_1.numat, &c__0, eigs);
    rsp_(&fmatrx[1], &geovar_1.nvar, &geovar_1.nvar, eigs, scrach_1.evecs);
    if (debug) {
	s_wsfe(&io___66);
	e_wsfe();
	matout_(scrach_1.evecs, eigs, &geovar_1.nvar, &geovar_1.nvar, &
		geovar_1.nvar);
    }
    l = 0;
    i__1 = geovar_1.nvar;
    for (i = 1; i <= i__1; ++i) {
	i__2 = i;
	for (j = 1; j <= i__2; ++j) {
	    ++l;
	    sum = 0.;
	    i__3 = *nreal;
	    for (k = 1; k <= i__3; ++k) {
		k1 = (k - 1) * geovar_1.nvar + i;
		k2 = (k - 1) * geovar_1.nvar + j;
/* L170: */
		sum += scrach_1.evecs[k1 - 1] * eigs[k - 1] * scrach_1.evecs[
			k2 - 1];
	    }
/* L180: */
	    fmatrx[l] = sum;
	}
    }
    frame_(&fmatrx[1], &molkst_1.numat, &c__0, eigs);
    rsp_(&fmatrx[1], &geovar_1.nvar, &geovar_1.nvar, eigs, scrach_1.evecs);
/* #      CALL MATOUT(EVECS,EIGS,NVAR,NVAR,NVAR) */
    jstart = 1;
    i__2 = geovar_1.nvar;
    for (i = 1; i <= i__2; ++i) {
/* L190: */
	cold[i - 1] = coordl[i - 1];
    }
L200:
    if (deriv) {
	goto L290;
    }
    l = (jstart - 1) * geovar_1.nvar;
    i__2 = geovar_1.nvar;
    for (iloop = jstart; iloop <= i__2; ++iloop) {

/*   MAKE THE STEP-SIZE ROUGHLY INVERSELY PROPORTIONAL TO THE EIGENVAL
UE */

	if (iloop <= *nreal) {
/* Computing MAX */
/* Computing MIN */
	    d__4 = .2, d__5 = .25 / (d__1 = eigs[iloop - 1], abs(d__1));
	    d__2 = .02, d__3 = min(d__4,d__5);
	    delta = max(d__2,d__3);
	} else if (iloop <= *nreal + 3) {
	    delta = .2;
	} else {
	    delta = .1;
	}
/* #      WRITE(6,*)'DELTA:',DELTA */
	j = l;
	compfg_(cold, &c_true, heat, &c_true, grad, &c_false);
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    ++j;
/* L210: */
	    coordl[i - 1] = cold[i - 1] + scrach_1.evecs[j - 1] * delta;
	}
	compfg_(coordl, &c_true, &heata, &c_true, grad, &c_false);
	heata -= *heat;
	j = l;
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    ++j;
/* L220: */
	    coordl[i - 1] = cold[i - 1] - scrach_1.evecs[j - 1] * delta;
	}
	compfg_(coordl, &c_true, &heatb, &c_true, grad, &c_false);
	heatb -= *heat;
	j = l;
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    ++j;
/* L230: */
	    coordl[i - 1] = cold[i - 1] + scrach_1.evecs[j - 1] * delta * 2;
	}
	compfg_(coordl, &c_true, &heataa, &c_true, grad, &c_false);
	heataa -= *heat;
	j = l;
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    ++j;
/* L240: */
	    coordl[i - 1] = cold[i - 1] - scrach_1.evecs[j - 1] * delta * 2;
	}
	compfg_(coordl, &c_true, &heatbb, &c_true, grad, &c_false);
	heatbb -= *heat;
	sum = ((heata + heatb) * 16 - (heataa + heatbb)) / 12. / delta * fact 
		/ delta * .5;
/* #      WRITE(6,'(5F12.6)')HEATAA+HEAT,HEATA+HEAT,HEAT,HEATB+HEAT, 
*/
/* #     +HEATBB+HEAT */
/* #      WRITE(6,'(5F12.6)')HEATAA,HEATA,0.D0,HEATB,HEATBB */
	fconst[iloop - 1] = sum * .5;
	l += geovar_1.nvar;
	time3 = second_();
	tstep = time3 - time2;
	if (tstep > 1e5) {
	    tstep += -1e6;
	    time3 += -1e6;
	    tleft = -1.;
	}
	time2 = time3;
	totime += tstep;
/* Computing MAX */
	d__1 = -1., d__2 = tleft - tstep;
	tleft = max(d__1,d__2);
	s_wsfe(&io___73);
	do_fio(&c__1, (char *)&iloop, (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&tstep, (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&totime, (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&tleft, (ftnlen)sizeof(doublereal));
	e_wsfe();
	estim = tstep * 5.;

/*    5.0 IS A SAFETY FACTOR */

	if (tlast - tleft > tdump) {
	    tlast = tleft;
	    resfil = TRUE_;
	    ifor = iloop;
	    ix = geovar_1.nvar + 2;

/* VALUE OF IX IS NOT IMPORTANT. SHOULD NOT BE 0 OR NVAR */

	    forsav_(&totime, &deldip[4], &ix, &geovar_1.nvar, &fmatrx[1], 
		    coord_1.coord, &geovar_1.nvar, heat, scrach_1.evecs, &
		    ifor, fconst);
	}
	if (iloop != geovar_1.nvar && tleft - 10. < estim) {
	    s_wsfe(&io___76);
	    e_wsfe();
	    s_wsfe(&io___77);
	    do_fio(&c__1, (char *)&iloop, (ftnlen)sizeof(integer));
	    e_wsfe();
	    s_wsfe(&io___78);
	    e_wsfe();
	    s_wsfe(&io___79);
	    do_fio(&c__1, (char *)&estim, (ftnlen)sizeof(doublereal));
	    e_wsfe();
	    ifor = iloop;
	    ix = geovar_1.nvar + 2;

/* VALUE OF IX IS NOT IMPORTANT. SHOULD NOT BE 0 OR NVAR */

	    forsav_(&totime, &deldip[4], &ix, &geovar_1.nvar, &fmatrx[1], 
		    coord_1.coord, &geovar_1.nvar, heat, scrach_1.evecs, &
		    ifor, fconst);
	}
/* L250: */
    }
    l = 0;
    i__2 = geovar_1.nvar;
    for (i = 1; i <= i__2; ++i) {
	i__1 = i;
	for (j = 1; j <= i__1; ++j) {
	    ++l;
	    sum = 0.;
	    i__3 = geovar_1.nvar;
	    for (k = 1; k <= i__3; ++k) {
		k1 = (k - 1) * geovar_1.nvar + i;
		k2 = (k - 1) * geovar_1.nvar + j;
/* L260: */
		sum += scrach_1.evecs[k1 - 1] * fconst[k - 1] * 
			scrach_1.evecs[k2 - 1];
	    }
/* L270: */
	    fmatrx[l] = sum * 2.;
	}
    }
    i__1 = geovar_1.nvar;
    for (i = 1; i <= i__1; ++i) {
/* L280: */
	coordl[i - 1] = cold[i - 1];
    }
L290:
    i__1 = molkst_1.numat;
    for (i = 1; i <= i__1; ++i) {
	if (atmass_1.atmass[i - 1] < 1e-20) {
	    forsav_(&totime, &deldip[4], &geovar_1.nvar, &geovar_1.nvar, &
		    fmatrx[1], coord_1.coord, &geovar_1.nvar, heat, 
		    scrach_1.evecs, &iloop, fconst);
	    s_wsfe(&io___80);
	    do_fio(&c__1, " AT LEAST ONE ATOM HAS A ZERO MASS. A RESTART", 
		    45L);
	    e_wsfe();
	    s_wsfe(&io___81);
	    do_fio(&c__1, " FILE HAS BEEN WRITTEN AND THE JOB STOPPED", 42L);
	    e_wsfe();
	    s_stop("", 0L);
	}
/* L300: */
    }
    if (istart <= geovar_1.nvar && i_indx(keywrd_1.keywrd, "ISOTOPE", 80L, 7L)
	     != 0) {
	forsav_(&totime, &deldip[4], &geovar_1.nvar, &geovar_1.nvar, &fmatrx[
		1], coord_1.coord, &geovar_1.nvar, heat, scrach_1.evecs, &
		iloop, fconst);
    }
    return 0;
} /* fmat_ */

#undef coordl