/* drc.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 {
    doublereal p[23220], pa[23220], pb[23220];
} densty_;

#define densty_1 densty_

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

#define gradnt_1 gradnt_

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 atmass[86];
} atmass_;

#define atmass_1 atmass_

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

#define geovar_1 geovar_

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 xract;
} molkst_;

#define molkst_1 molkst_

struct {
    integer mcoprt[516]	/* was [2][258] */, ncoprt;
    logical prtmax;
} drccom_;

#define drccom_1 drccom_

/* Table of constant values */

static integer c__9 = 9;
static integer c__1 = 1;
static integer c__4 = 4;
static integer c__3 = 3;
static integer c__5 = 5;
static integer c__8 = 8;
static doublereal c_b120 = .5;
static logical c_true = TRUE_;
static doublereal c_b133 = .25;

/* Subroutine */ int drc_(doublereal *startv, doublereal *startk)
{
    /* Initialized data */

    static doublereal velo0[258] = { 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	    0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	    0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	    0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	    0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	    0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	    0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	    0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	    0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	    0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	    0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	    0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	    0.,0.,0. };
    static logical addk = TRUE_;

    /* System generated locals */
    integer i__1, i__2, i__3;
    doublereal d__1, d__2, d__3, d__4;
    olist o__1;
    cllist cl__1;
    alist al__1;

    /* Builtin functions */
    integer f_clos(cllist *), f_open(olist *), i_indx(char *, char *, ftnlen, 
	    ftnlen), s_wsle(cilist *), do_lio(integer *, integer *, char *, 
	    ftnlen), e_wsle(void);
    double sqrt(doublereal);
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
    double d_sign(doublereal *, doublereal *);
    integer s_cmp(char *, char *, ftnlen, ftnlen), f_rew(alist *), s_rsfe(
	    cilist *), e_rsfe(void), s_rsle(cilist *), e_rsle(void);
    double pow_dd(doublereal *, doublereal *);
    /* Subroutine */ int s_stop(char *, ftnlen);
    integer s_wsue(cilist *), do_uio(integer *, char *, ftnlen), e_wsue(void);

    /* Local variables */
    static doublereal half, accu, escf, once, ekin, etot, gtot, summ, tnow, 
	    velo1[258], velo2[258], velo3[258];
    static integer i, j, k, l;
    extern doublereal reada_(char *, integer *, ftnlen);
    static doublereal alpha, ekold, coord[258]	/* was [3][86] */, grold[258],
	     etold, quadr, tleft, const__;
    static logical letot;
    static integer iloop, i1;
    static doublereal error, elost, delta1, dlold2, grold2[258], elost1;
    static char ch[1];
    static integer ii, kl;
    static doublereal addonk, delold, georef[258]	/* was [3][86] */;
    static logical ircdrc;
    extern doublereal second_(void);
    static doublereal deltat;
    extern /* Subroutine */ int getarg_(integer *, char *, ftnlen);
    static doublereal scfold, parold[258];
    static logical velred;
    static doublereal velref[258], velvec, oldtim, gerror[258], totime;
    static integer lpoint;
    extern /* Subroutine */ int gmetry_(doublereal *, doublereal *);
    static doublereal summas;
    static integer iupper;
    extern /* Subroutine */ int compfg_(doublereal *, logical *, doublereal *,
	     logical *, doublereal *, logical *), prtdrc_(doublereal *, 
	    doublereal *, doublereal *, doublereal *, doublereal *, 
	    doublereal *, doublereal *, doublereal *, integer *);
    static doublereal tcycle;
    static integer linear;
    static doublereal sqrtms[258], ams, one;
    static logical let;
    static integer ilp;
    extern doublereal dot_(doublereal *, doublereal *, integer *);
    static doublereal tim, sum;

    /* Fortran I/O blocks */
    static cilist io___16 = { 0, 6, 0, 0, 0 };
    static cilist io___26 = { 0, 6, 0, "(//10X,' DAMPING FACTOR FOR KINETIC "
	    "ENERGY =',F12.6)", 0 };
    static cilist io___29 = { 0, 6, 0, "(//10X,' EXCESS KINETIC ENERGY ENTER"
	    "ED INTO SYSTEM = ',F12.6)", 0 };
    static cilist io___40 = { 0, 9, 0, "(A)", 0 };
    static cilist io___42 = { 0, 9, 0, "(3F19.13)", 0 };
    static cilist io___43 = { 0, 9, 0, "(A)", 0 };
    static cilist io___44 = { 0, 9, 0, "(3F19.3)", 0 };
    static cilist io___45 = { 0, 9, 0, "(A)", 0 };
    static cilist io___46 = { 0, 9, 0, 0, 0 };
    static cilist io___47 = { 0, 9, 0, 0, 0 };
    static cilist io___49 = { 0, 9, 0, 0, 0 };
    static cilist io___51 = { 0, 9, 0, 0, 0 };
    static cilist io___56 = { 0, 6, 0, "(//10X,'CALCULATION RESTARTED, CURRE"
	    "NT',            ' KINETIC ENERGY=',F10.5,//)", 0 };
    static cilist io___67 = { 0, 6, 0, "(A,F10.3,A,/,A)", 0 };
    static cilist io___68 = { 0, 6, 0, "(A,F7.2,A)", 0 };
    static cilist io___69 = { 0, 6, 0, "(A,F7.2,A)", 0 };
    static cilist io___70 = { 0, 6, 0, "(A)", 0 };
    static cilist io___84 = { 0, 6, 0, "(//,' IRC CALCULATION COMPLETE ')", 0 
	    };
    static cilist io___86 = { 0, 9, 0, "(A)", 0 };
    static cilist io___87 = { 0, 9, 0, "(3F19.13)", 0 };
    static cilist io___88 = { 0, 9, 0, "(A)", 0 };
    static cilist io___89 = { 0, 9, 0, "(3F19.3)", 0 };
    static cilist io___90 = { 0, 9, 0, "(A)", 0 };
    static cilist io___91 = { 0, 9, 0, 0, 0 };
    static cilist io___92 = { 0, 9, 0, 0, 0 };
    static cilist io___93 = { 0, 9, 0, 0, 0 };
    static cilist io___94 = { 0, 9, 0, 0, 0 };
    static cilist io___96 = { 0, 10, 0, 0, 0 };
    static cilist io___97 = { 0, 10, 0, 0, 0 };
    static cilist io___98 = { 0, 6, 0, "(//10X,' RUNNING OUT OF TIME, RESTAR"
	    "T FILE WRITTEN')", 0 };


/* ***********************************************************************
 */
/*                                                                      * 
*/
/*    DRC IS DESIGNED TO FOLLOW A REACTION PATH FROM THE TRANSITION     * 
*/
/*    STATE.  TWO MODES ARE SUPPORTED, FIRST: GAS PHASE:- AS THE SYSTEM * 
*/
/*    MOVES FROM THE T/S THE MOMENTUM OF THE ATOMS IS STORED AND THE    * 
*/
/*    POSITION OF THE ATOMS IS RELATED TO THE OLD POSITION BY (A) THE   * 
*/
/*    CURRENT VELOCITY OF THE ATOM, AND (B) THE FORCES ACTING ON THAT   * 
*/
/*    ATOM.  THE SECOND MODE IS CONDENSED PHASE, IN WHICH THE ATOMS MOVE* 
*/
/*    IN RESPONSE TO THE FORCES ACTING ON THEM. I.E. INFINITELY DAMPED  * 
*/
/*                                                                      * 
*/
/* ***********************************************************************
 */
/* 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 */

    /* Parameter adjustments */
    --startk;
    --startv;

    /* Function Body */
    cl__1.cerr = 0;
    cl__1.cunit = 5;
    cl__1.csta = 0;
    f_clos(&cl__1);
    tnow = second_();
    oldtim = second_();
    delold = 10.;
    gtot = 0.;
    o__1.oerr = 0;
    o__1.ounit = 7;
    o__1.ofnm = 0;
    o__1.orl = 0;
    o__1.osta = "SCRATCH";
    o__1.oacc = 0;
    o__1.ofm = 0;
    o__1.oblnk = 0;
    f_open(&o__1);
    if (i_indx(keywrd_1.keywrd, " PREC", 80L, 5L) != 0) {
	accu = .5;
    } else {
	accu = 1.;
    }
    lpoint = 0;
    velred = i_indx(keywrd_1.keywrd, "VELOC", 80L, 5L) != 0;

/*     PRINT OUT INITIAL VELOCITIES */

/* #      WRITE(6,'(A)')' INITIAL VELOCITY IN DRC' */
/* #      WRITE(6,'(3F13.5)')(STARTV(I),I=1,NUMAT*3) */
    let = i_indx(keywrd_1.keywrd, " GEO-OK", 80L, 7L) != 0;
    if (i_indx(keywrd_1.keywrd, " SYMME", 80L, 6L) != 0) {
	s_wsle(&io___16);
	do_lio(&c__9, &c__1, "  SYMMETRY SPECIFIED, BUT CANNOT BE USED IN DRC"
		, 47L);
	e_wsle();
	geosym_1.ndep = 0;
    }

/*      CONVERT TO CARTESIAN COORDINATES, IF NOT ALREADY DONE. */

    if (i_indx(keywrd_1.keywrd, " XYZ", 80L, 4L) == 0) {
	geokst_1.na[0] = 0;
	gmetry_(geom_1.geo, coord);
	l = 0;
	i__1 = molkst_1.numat;
	for (i = 1; i <= i__1; ++i) {
	    geokst_1.labels[i - 1] = molkst_1.nat[i - 1];
	    sum = sqrt(atmass_1.atmass[molkst_1.nat[i - 1] - 1]);
	    for (j = 1; j <= 3; ++j) {
		++l;
		sqrtms[l - 1] = sum;
		geom_1.geo[j + i * 3 - 4] = coord[j + i * 3 - 4];
/* L10: */
		coord[j + i * 3 - 4] = 0.;
	    }
/* L20: */
	}
	geokst_1.na[0] = 99;
    }

/*  TRANSFER COORDINATES TO XPARAM AND LOC */

    if (i_indx(keywrd_1.keywrd, " DRC", 80L, 4L) != 0) {
	drccom_1.prtmax = geovar_1.loc[0] == 1;
	if (drccom_1.prtmax) {
	    j = 1;
	} else {
	    j = 0;
	}
	geovar_1.nvar -= j;
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    drccom_1.mcoprt[(i << 1) - 2] = geovar_1.loc[(i + j << 1) - 2];
/* L30: */
	    drccom_1.mcoprt[(i << 1) - 1] = geovar_1.loc[(i + j << 1) - 1];
	}
	if (geovar_1.loc[0] == 0) {
	    geovar_1.nvar = 0;
	}
	drccom_1.ncoprt = geovar_1.nvar;
    } else {
	drccom_1.ncoprt = 0;
    }
    l = 0;
    i__1 = molkst_1.numat;
    for (i = 1; i <= i__1; ++i) {
	for (j = 1; j <= 3; ++j) {
	    ++l;
	    geovar_1.loc[(l << 1) - 2] = i;
	    geovar_1.loc[(l << 1) - 1] = j;
	    georef[j + i * 3 - 4] = geom_1.geo[j + i * 3 - 4];
/* L40: */
	    geovar_1.xparam[l - 1] = geom_1.geo[j + i * 3 - 4];
	}
    }
    geovar_1.nvar = molkst_1.numat * 3;

/* DETERMINE DAMPING FACTOR */

    ircdrc = i_indx(keywrd_1.keywrd, "IRC=", 80L, 4L) != 0;
    if (i_indx(keywrd_1.keywrd, "DRC=", 80L, 4L) != 0) {
	i__1 = i_indx(keywrd_1.keywrd, "DRC=", 80L, 4L);
	half = reada_(keywrd_1.keywrd, &i__1, 80L);
	s_wsfe(&io___26);
	do_fio(&c__1, (char *)&half, (ftnlen)sizeof(doublereal));
	e_wsfe();
    } else if (i_indx(keywrd_1.keywrd, "DRC", 80L, 3L) == 0) {
	half = 0.;
    } else {
	half = 1e6;
    }
    letot = ! ircdrc && half < 1.;
/* Computing MAX */
    d__2 = 1e-6, d__3 = abs(half);
    d__1 = max(d__2,d__3);
    half = d_sign(&d__1, &half);

/* DETERMINE EXCESS KINETIC ENERGY */

    if (i_indx(keywrd_1.keywrd, "KINE", 80L, 4L) != 0) {
	i__1 = i_indx(keywrd_1.keywrd, "KINE", 80L, 4L);
	addonk = reada_(keywrd_1.keywrd, &i__1, 80L);
	s_wsfe(&io___29);
	do_fio(&c__1, (char *)&addonk, (ftnlen)sizeof(doublereal));
	e_wsfe();
    } else {
	addonk = 0.;
    }

/*   LOOP OVER TIME-INTERVALS OF DELTAT SECOND */

    deltat = 1e-16;
    quadr = 1.;
    totime = 0.;
    once = 0.;
    etot = 0.;
    escf = 0.;
    const__ = 1.;
    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 L60;
	    }
/* L50: */
	}
/*             4 SECONDS TO LOAD IN EXECUTABLE! */
L60:
	tleft = tim - 4;
    } else {
	tleft = 3596.;
    }
    if (i_indx(keywrd_1.keywrd, "REST", 80L, 4L) != 0 && i_indx(
	    keywrd_1.keywrd, "IRC=", 80L, 4L) == 0) {

/*  RESTART FROM A PREVIOUS RUN */

	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 = "FORMATTED";
	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);
	s_rsfe(&io___40);
	do_fio(&c__1, (char *)&alpha, (ftnlen)sizeof(doublereal));
	e_rsfe();
	s_rsfe(&io___42);
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    do_fio(&c__1, (char *)&geovar_1.xparam[i - 1], (ftnlen)sizeof(
		    doublereal));
	}
	e_rsfe();
	s_rsfe(&io___43);
	do_fio(&c__1, (char *)&alpha, (ftnlen)sizeof(doublereal));
	e_rsfe();
	s_rsfe(&io___44);
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    do_fio(&c__1, (char *)&velo0[i - 1], (ftnlen)sizeof(doublereal));
	}
	e_rsfe();
	s_rsfe(&io___45);
	do_fio(&c__1, (char *)&alpha, (ftnlen)sizeof(doublereal));
	e_rsfe();
	s_rsle(&io___46);
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    do_lio(&c__5, &c__1, (char *)&gradnt_1.grad[i - 1], (ftnlen)
		    sizeof(doublereal));
	}
	e_rsle();
	s_rsle(&io___47);
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    do_lio(&c__5, &c__1, (char *)&grold[i - 1], (ftnlen)sizeof(
		    doublereal));
	}
	e_rsle();
	s_rsle(&io___49);
	i__1 = geovar_1.nvar;
	for (i = 1; i <= i__1; ++i) {
	    do_lio(&c__5, &c__1, (char *)&grold2[i - 1], (ftnlen)sizeof(
		    doublereal));
	}
	e_rsle();
	s_rsle(&io___51);
	do_lio(&c__5, &c__1, (char *)&etot, (ftnlen)sizeof(doublereal));
	do_lio(&c__5, &c__1, (char *)&escf, (ftnlen)sizeof(doublereal));
	do_lio(&c__5, &c__1, (char *)&ekin, (ftnlen)sizeof(doublereal));
	do_lio(&c__5, &c__1, (char *)&delold, (ftnlen)sizeof(doublereal));
	do_lio(&c__5, &c__1, (char *)&deltat, (ftnlen)sizeof(doublereal));
	do_lio(&c__5, &c__1, (char *)&dlold2, (ftnlen)sizeof(doublereal));
	do_lio(&c__3, &c__1, (char *)&iloop, (ftnlen)sizeof(integer));
	do_lio(&c__5, &c__1, (char *)&gradnt_1.gnorm, (ftnlen)sizeof(
		doublereal));
	do_lio(&c__8, &c__1, (char *)&letot, (ftnlen)sizeof(logical));
	do_lio(&c__5, &c__1, (char *)&elost1, (ftnlen)sizeof(doublereal));
	do_lio(&c__5, &c__1, (char *)&gtot, (ftnlen)sizeof(doublereal));
	e_rsle();
	s_wsfe(&io___56);
	do_fio(&c__1, (char *)&ekin, (ftnlen)sizeof(doublereal));
	e_wsfe();
	goto L110;
    } else {
	iloop = 1;
	if (i_indx(keywrd_1.keywrd, "IRC=", 80L, 4L) != 0 || velred) {
	    if (i_indx(keywrd_1.keywrd, "IRC=", 80L, 4L) != 0) {
		i__1 = i_indx(keywrd_1.keywrd, "IRC=", 80L, 4L);
		k = (integer) reada_(keywrd_1.keywrd, &i__1, 80L);
	    } else {
		k = 1;
	    }
	    if (k < 0) {
		k = -k;
		one = -1.;
	    } else {
		one = 1.;
	    }
	    kl = (k - 1) * geovar_1.nvar;
	    summ = 0.;
	    velo1[0] = 0.;
	    velo1[1] = 0.;
	    velo1[2] = 0.;
	    summas = 0.;
	    i = 0;
	    i__1 = molkst_1.numat;
	    for (ii = 1; ii <= i__1; ++ii) {
		ams = atmass_1.atmass[ii - 1];
		summas += ams;
		for (i1 = 1; i1 <= 3; ++i1) {
		    ++i;
		    velo0[i - 1] = startv[kl + i] * one;
		    velref[i - 1] = velo0[i - 1];
		    velo1[i1 - 1] += velo0[i - 1] * ams;
/* L70: */
		}
	    }
	    for (i = 1; i <= 3; ++i) {
/* L80: */
		velo1[i - 1] = -velo1[i - 1] / summas;
	    }
	    i = 0;
	    i__1 = molkst_1.numat;
	    for (ii = 1; ii <= i__1; ++ii) {
		ams = atmass_1.atmass[ii - 1];
		for (i1 = 1; i1 <= 3; ++i1) {
		    ++i;
		    if (addonk > 1e-5 || ! velred) {
			velo0[i - 1] += velo1[i1 - 1];
		    }
/* L90: */
/* Computing 2nd power */
		    d__1 = velo0[i - 1];
		    summ += d__1 * d__1 * ams;
		}
	    }
	    if (addonk < 1e-5 && velred) {
		addonk = summ * .5 / 4.184e10;
	    }
	    if (addonk < 1e-5 && ! velred) {
		if (abs(half) > .001 && startk[k] > 105.) {
		    s_wsfe(&io___67);
		    do_fio(&c__1, " BY DEFAULT, ONE QUANTUM OF ENERGY, EQUIV"
			    "ALENT TO", 49L);
		    do_fio(&c__1, (char *)&startk[k], (ftnlen)sizeof(
			    doublereal));
		    do_fio(&c__1, " CM(-1)", 7L);
		    do_fio(&c__1, " WILL BE USED TO START THE DRC", 30L);
		    e_wsfe();

/*    2.8585086D-3 CONVERTS CM(-1) INTO KCAL/MOLE */

		    addonk = startk[k] * .0028585086;
		    s_wsfe(&io___68);
		    do_fio(&c__1, " THIS REPRESENTS AN ENERGY OF", 29L);
		    do_fio(&c__1, (char *)&addonk, (ftnlen)sizeof(doublereal))
			    ;
		    do_fio(&c__1, " KCALS/MOLE", 11L);
		    e_wsfe();
		} else if (abs(half) > .001) {
		    s_wsfe(&io___69);
		    do_fio(&c__1, " THE VIBRATIONAL FREQUENCY (", 28L);
		    do_fio(&c__1, (char *)&startk[k], (ftnlen)sizeof(
			    doublereal));
		    do_fio(&c__1, " IS TOO SMALL FOR ONE QUANTUM TO BE USED", 
			    40L);
		    e_wsfe();
		    s_wsfe(&io___70);
		    do_fio(&c__1, " INSTEAD 0.3KCAL/MOLE WILL BE USED TO STA"
			    "RT THE IRC", 51L);
		    e_wsfe();
		    addonk = .3;
		} else {
		    addonk = .3;
		}
	    }

/*   AT THIS POINT ADDONK IS IN KCAL/MOLE */
/*   NORMALIZE SO THAT TOTAL K.E. = ONE QUANTUM (DEFAULT) (DRC ONL
Y) */
/*                              OR 0.3KCAL/MOLE (IRC ONLY) */
/*                              OR ADDONK IF KINETIC=NN SUPPLIED 
*/

	    summ = sqrt(addonk / (summ * .5 / 4.184e10));
	    addk = FALSE_;
	    if (summ > 1e-10) {
		i__1 = geovar_1.nvar;
		for (i = 1; i <= i__1; ++i) {
/* L100: */
		    velo0[i - 1] *= summ;
		}
	    }
	}
    }
L110:
    iupper = iloop + 4999;
    ilp = iloop;
    one = 0.;
    if (i_indx(keywrd_1.keywrd, "REST", 80L, 4L) != 0 && i_indx(
	    keywrd_1.keywrd, "IRC=", 80L, 4L) == 0) {
	one = 1.;
    }
    i__1 = iupper;
    for (iloop = ilp; iloop <= i__1; ++iloop) {

/*  MOVEMENT OF ATOMS WILL BE PROPORTIONAL TO THE AVERAGE VELOCITIES 
*/
/*  OF THE ATOMS BEFORE AND AFTER TIME INTERVAL */


/*  RAPID CHANGE IN GRADIENT IMPLIES SMALL STEP SIZE FOR DELTAT */

/*   KINETIC ENERGY = 1/2 * M * V * V */
/*                  = 0.5 / (4.184D10) * M * V * V */
/*   NEW VELOCITY = OLD VELOCITY + GRADIENT * TIME / MASS */
/*                = KCAL/ANGSTROM*SECOND/(ATOMIC WEIGHT) */
/*                =4.184*10**10(ERGS)*10**8(PER CM)*DELTAT(SECONDS) */
/*   NEW POSITION = OLD POSITION - AVERAGE VELOCITY * TIME INTERVAL */


/*   ESTABLISH REFERENCE TOTAL ENERGY */

	error = etot - (ekin + escf);
	if (iloop > 2) {
	    quadr = error / (ekin * const__ + .001) * .5 + 1.;
/* Computing MIN */
	    d__1 = 1.3, d__2 = max(.8,quadr);
	    quadr = min(d__1,d__2);
	} else {
	    quadr = 1.;
	}
	if ((let || ekin > .2f) && addk) {

/*   DUMP IN EXCESS KINETIC ENERGY */

	    etot += addonk;
	    addk = FALSE_;
	    addonk = 0.;
	}
	ekold = ekin;

/*  CALCULATE THE DURATION OF THE NEXT STEP. */
/*  STEP SIZE IS THAT REQUIRED TO PRODUCE A CONSTANT CHANGE IN GEOMETR
Y */


/*  IF DAMPING IS USED, CALCULATE THE NEW TOTAL ENERGY AND */
/*  THE RATIO FOR REDUCING THE KINETIC ENERGY */

/* Computing MAX */
	d__3 = deltat * 1e15 / half;
	d__1 = 1e-36, d__2 = pow_dd(&c_b120, &d__3);
	const__ = max(d__1,d__2);
	const__ = sqrt(const__);
	velvec = 0.;
	ekin = 0.;
	delta1 = delold + dlold2;
	elost = 0.;
	i__2 = geovar_1.nvar;
	for (i = 1; i <= i__2; ++i) {

/*   CALCULATE COMPONENTS OF VELOCITY AS */
/*   V = V(0) + V'*T + V"*T*T */
/*   WE NEED ALL THREE TERMS, V(0), V' AND V" */

	    velo1[i - 1] = 1. / atmass_1.atmass[geovar_1.loc[(i << 1) - 2] - 
		    1] * gradnt_1.grad[i - 1];
	    if (iloop > 3) {
/* Computing 2nd power */
		d__1 = delold;
/* Computing 2nd power */
		d__2 = delta1;
		velo3[i - 1] = 2. / atmass_1.atmass[geovar_1.loc[(i << 1) - 2]
			 - 1] * (delta1 * (grold[i - 1] - gradnt_1.grad[i - 1]
			) - delold * (grold2[i - 1] - gradnt_1.grad[i - 1])) /
			 (delta1 * (d__1 * d__1 * 1e30) - delold * (d__2 * 
			d__2 * 1e30));
/* Computing 2nd power */
		d__1 = delold;
		velo2[i - 1] = 1. / atmass_1.atmass[geovar_1.loc[(i << 1) - 2]
			 - 1] * (gradnt_1.grad[i - 1] - grold[i - 1] - velo3[
			i - 1] * .5 * (d__1 * d__1 * 1e30)) / (delold * 1e15);
	    } else {
		velo2[i - 1] = 1. / atmass_1.atmass[geovar_1.loc[(i << 1) - 2]
			 - 1] * (gradnt_1.grad[i - 1] - grold[i - 1]) / (
			delold * 1e15);
		velo3[i - 1] = 0.;
	    }

/*  MOVE ATOMS THROUGH DISTANCE EQUAL TO VELOCITY * DELTA-TIME, NO
TE */
/*  VELOCITY CHANGES FROM START TO FINISH, THEREFORE AVERAGE. */

	    parold[i - 1] = geovar_1.xparam[i - 1];
/* Computing 2nd power */
	    d__1 = deltat;
/* Computing 2nd power */
	    d__2 = deltat;
/* Computing 2nd power */
	    d__3 = deltat;
/* Computing 2nd power */
	    d__4 = deltat;
	    geovar_1.xparam[i - 1] -= (deltat * velo0[i - 1] * one + d__1 * 
		    d__1 * .5 * velo1[i - 1] + d__2 * d__2 * 1e15 * .16666 * 
		    deltat * velo2[i - 1] + d__3 * d__3 * .0416666 * (d__4 * 
		    d__4 * 1e30) * velo3[i - 1]) * 1e8;

/*   CORRECT ERRORS DUE TO CUBIC COMPONENTS IN ENERGY GRADIENT, */
/*   ALSO TO ADD ON EXCESS ENERGY, IF NECESSARY. */

/* Computing 2nd power */
	    d__1 = velo0[i - 1];
	    velvec += d__1 * d__1;

/*   MODIFY VELOCITY IN LIGHT OF CURRENT ENERGY GRADIENTS. */

/*   VELOCITY = OLD VELOCITY + (DELTA-T / ATOMIC MASS) * CURRENT G
RADIENT */
/*                           + 1/2 *(DELTA-T * DELTA-T /ATOMIC MAS
S) * */
/*                             (SLOPE OF GRADIENT) */
/*              SLOPE OF GRADIENT = (GRAD(I)-GROLD(I))/DELOLD */


/*   THIS EXPRESSION IS ACCURATE TO SECOND ORDER IN TIME. */

/* Computing 2nd power */
	    d__1 = deltat;
/* Computing 2nd power */
	    d__2 = deltat;
	    velo0[i - 1] = velo0[i - 1] + deltat * velo1[i - 1] + d__1 * d__1 
		    * .5 * velo2[i - 1] * 1e15 + deltat * .166666 * (d__2 * 
		    d__2 * 1e30) * velo3[i - 1];
	    if (let || gradnt_1.gnorm > 3.) {
		let = TRUE_;
/* Computing 2nd power */
		d__1 = velo0[i - 1];
/* Computing 2nd power */
		d__2 = const__;
		elost += d__1 * d__1 * atmass_1.atmass[geovar_1.loc[(i << 1) 
			- 2] - 1] * (1 - d__2 * d__2);
		velo0[i - 1] = velo0[i - 1] * const__ * quadr;
	    }

/*  CALCULATE KINETIC ENERGY (IN 2*ERGS AT THIS POINT) */

/* Computing 2nd power */
	    d__1 = velo0[i - 1];
	    ekin += d__1 * d__1 * atmass_1.atmass[geovar_1.loc[(i << 1) - 2] 
		    - 1];
/* L120: */
	}
	one = 1.;
	if (let || gradnt_1.gnorm > 3.) {
	    if (! letot) {
		etot = escf + addonk;
		addonk = 0.;
		elost1 = 0.;
		elost = 0.;
	    }
	    letot = TRUE_;
	}

/*  CONVERT ENERGY INTO KCAL/MOLE */

	ekin = ekin * .5f / 4.184e10;
	if (letot && abs(half) > 1e-5) {
/* Computing 2nd power */
	    d__1 = const__;
	    etot = etot - ekin / (d__1 * d__1) + ekin;
	}
	elost1 += elost * .5 / 4.184e10;

/* STORE OLD GRADIENTS FOR DELTA - VELOCITY CALCULATION */

	i__2 = geovar_1.nvar;
	for (i = 1; i <= i__2; ++i) {
	    grold2[i - 1] = grold[i - 1];
	    grold[i - 1] = gradnt_1.grad[i - 1];
/* L130: */
	    gradnt_1.grad[i - 1] = 0.;
	}

/*   CALCULATE ENERGY AND GRADIENTS */

	scfold = escf;
	compfg_(geovar_1.xparam, &c_true, &escf, &c_true, gradnt_1.grad, &
		c_true);
	if (iloop > 2) {
	    gradnt_1.gnorm = 0.;
	    i__2 = geovar_1.nvar;
	    for (i = 1; i <= i__2; i += 3) {
		sum = sqrt(dot_(&gradnt_1.grad[i - 1], &gradnt_1.grad[i - 1], 
			&c__3) / (dot_(&velo0[i - 1], &velo0[i - 1], &c__3) + 
			1e-20));
		i__3 = i + 2;
		for (j = i; j <= i__3; ++j) {
/* L140: */
		    gerror[j - 1] = gerror[j - 1] + gradnt_1.grad[j - 1] + 
			    velo0[j - 1] * sum;
		}
/* L150: */
	    }
	    gradnt_1.gnorm = sqrt(dot_(gerror, gerror, &geovar_1.nvar));
	    gtot = gradnt_1.gnorm;
	}
	gradnt_1.gnorm = sqrt(dot_(gradnt_1.grad, gradnt_1.grad, &
		geovar_1.nvar));

/*   CONVERT GRADIENTS INTO ERGS/CM */

	i__2 = geovar_1.nvar;
	for (i = 1; i <= i__2; ++i) {
/* L160: */
	    gradnt_1.grad[i - 1] *= 4.184e18;
	}

/*   SPECIAL TREATMENT FOR FIRST POINT - SET "OLD" GRADIENTS EQUAL TO 
*/
/*   CURRENT GRADIENTS. */

	if (iloop == 1) {
	    i__2 = geovar_1.nvar;
	    for (i = 1; i <= i__2; ++i) {
/* L170: */
		grold[i - 1] = gradnt_1.grad[i - 1];
	    }
	}
	dlold2 = delold;
	delold = deltat;
	sum = 0.;
	i__2 = geovar_1.nvar;
	for (i = 1; i <= i__2; ++i) {
/* L180: */
/* Computing 2nd power */
	    d__1 = (gradnt_1.grad[i - 1] - grold[i - 1]) / 4.184e18;
	    sum += d__1 * d__1;
	}
	if (abs(half) < .001) {
/* Computing MIN */
	    d__3 = 2., d__4 = accu * 5e-5 / ((d__1 = escf + elost1 - etold, 
		    abs(d__1)) + 1e-20);
	    d__2 = min(d__3,d__4);
	    deltat *= pow_dd(&d__2, &c_b133);
	    etold = escf + elost1;
	    if (iloop > 5 && scfold - escf < -.001 || iloop > 30 && scfold - 
		    escf < 0.) {
		s_wsfe(&io___84);
		e_wsfe();
		s_stop("", 0L);
	    }
	} else {
/* Computing MIN */
	    d__1 = 1.05, d__2 = accu * 10. / (sum + 1e-4);
	    deltat *= min(d__1,d__2);
/* **************************************************************
********* */

/*         TESTING CODE - REMOVE BEFORE FINAL VERSION ASSEMBLED */
/* #          (ILOOP/400)*400.EQ.ILOOP)DELTAT=-DELTAT */

/* **************************************************************
********* */
	}
	deltat = max(1e-16,deltat);
	if (abs(half) < 1e-5) {
	    prtdrc_(&escf, &deltat, geovar_1.xparam, georef, &elost1, &gtot, &
		    etot, velo0, &geovar_1.nvar);
	} else {
	    prtdrc_(&escf, &deltat, geovar_1.xparam, georef, &ekin, &elost, &
		    etot, velo0, &geovar_1.nvar);
	}
	tnow = second_();
	tcycle = tnow - oldtim;
	oldtim = tnow;
	tleft -= tcycle;
	if (iloop == iupper || tleft < tcycle * 3) {
	    if (i_indx(keywrd_1.keywrd, "REST", 80L, 4L) + i_indx(
		    keywrd_1.keywrd, "ISOT", 80L, 4L) != 0) {
		cl__1.cerr = 0;
		cl__1.cunit = 9;
		cl__1.csta = 0;
		f_clos(&cl__1);
	    }
	    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 = "FORMATTED";
	    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 = "UNFORMATTE    D";
	    o__1.oblnk = 0;
	    f_open(&o__1);
	    al__1.aerr = 0;
	    al__1.aunit = 10;
	    f_rew(&al__1);
	    s_wsfe(&io___86);
	    do_fio(&c__1, " CARTESIAN GEOMETRY PARAMETERS IN ANGSTROMS", 43L);
	    e_wsfe();
	    s_wsfe(&io___87);
	    i__2 = geovar_1.nvar;
	    for (i = 1; i <= i__2; ++i) {
		do_fio(&c__1, (char *)&geovar_1.xparam[i - 1], (ftnlen)sizeof(
			doublereal));
	    }
	    e_wsfe();
	    s_wsfe(&io___88);
	    do_fio(&c__1, " VELOCITY FOR EACH CARTESIAN COORDINATE, IN CM/SEC"
		    , 50L);
	    e_wsfe();
	    s_wsfe(&io___89);
	    i__2 = geovar_1.nvar;
	    for (i = 1; i <= i__2; ++i) {
		do_fio(&c__1, (char *)&velo0[i - 1], (ftnlen)sizeof(
			doublereal));
	    }
	    e_wsfe();
	    s_wsfe(&io___90);
	    do_fio(&c__1, " FIRST, SECOND, AND THIRD-ORDER GRADIENTS, ETC", 
		    46L);
	    e_wsfe();
	    s_wsle(&io___91);
	    i__2 = geovar_1.nvar;
	    for (i = 1; i <= i__2; ++i) {
		do_lio(&c__5, &c__1, (char *)&gradnt_1.grad[i - 1], (ftnlen)
			sizeof(doublereal));
	    }
	    e_wsle();
	    s_wsle(&io___92);
	    i__2 = geovar_1.nvar;
	    for (i = 1; i <= i__2; ++i) {
		do_lio(&c__5, &c__1, (char *)&grold[i - 1], (ftnlen)sizeof(
			doublereal));
	    }
	    e_wsle();
	    s_wsle(&io___93);
	    i__2 = geovar_1.nvar;
	    for (i = 1; i <= i__2; ++i) {
		do_lio(&c__5, &c__1, (char *)&grold2[i - 1], (ftnlen)sizeof(
			doublereal));
	    }
	    e_wsle();
	    i = iloop + 1;
	    s_wsle(&io___94);
	    do_lio(&c__5, &c__1, (char *)&etot, (ftnlen)sizeof(doublereal));
	    do_lio(&c__5, &c__1, (char *)&escf, (ftnlen)sizeof(doublereal));
	    do_lio(&c__5, &c__1, (char *)&ekin, (ftnlen)sizeof(doublereal));
	    do_lio(&c__5, &c__1, (char *)&delold, (ftnlen)sizeof(doublereal));
	    do_lio(&c__5, &c__1, (char *)&deltat, (ftnlen)sizeof(doublereal));
	    do_lio(&c__5, &c__1, (char *)&dlold2, (ftnlen)sizeof(doublereal));
	    do_lio(&c__3, &c__1, (char *)&i, (ftnlen)sizeof(integer));
	    do_lio(&c__5, &c__1, (char *)&gradnt_1.gnorm, (ftnlen)sizeof(
		    doublereal));
	    do_lio(&c__8, &c__1, (char *)&letot, (ftnlen)sizeof(logical));
	    do_lio(&c__5, &c__1, (char *)&elost1, (ftnlen)sizeof(doublereal));
	    do_lio(&c__5, &c__1, (char *)&gtot, (ftnlen)sizeof(doublereal));
	    e_wsle();
	    escf = -1e9;
	    prtdrc_(&escf, &deltat, geovar_1.xparam, georef, &ekin, &elost, &
		    etot, velo0, &geovar_1.nvar);
	    linear = molkst_1.norbs * (molkst_1.norbs + 1) / 2;
	    s_wsue(&io___96);
	    i__2 = linear;
	    for (i = 1; i <= i__2; ++i) {
		do_uio(&c__1, (char *)&densty_1.pa[i - 1], (ftnlen)sizeof(
			doublereal));
	    }
	    e_wsue();
	    if (molkst_1.nalpha != 0) {
		s_wsue(&io___97);
		i__2 = linear;
		for (i = 1; i <= i__2; ++i) {
		    do_uio(&c__1, (char *)&densty_1.pb[i - 1], (ftnlen)sizeof(
			    doublereal));
		}
		e_wsue();
	    }
	    s_wsfe(&io___98);
	    e_wsfe();
	    s_stop("", 0L);
	}
/* L190: */
    }
    return 0;
} /* drc_ */