/* react1.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 {
    doublereal geo[258]	/* was [3][86] */;
} geom_;

#define geom_1 geom_

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

#define geokst_1 geokst_

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

#define densty_1 densty_

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

#define geosym_1 geosym_

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

#define geovar_1 geovar_

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

#define gradnt_1 gradnt_

struct {
    doublereal ams[107];
} istope_;

#define istope_1 istope_

struct {
    doublereal cosine;
} gravec_;

#define gravec_1 gravec_

struct {
    char keywrd[80];
} keywrd_;

#define keywrd_1 keywrd_

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 step, geoa[258]	/* was [3][86] */, geovec[258]	/* was [3][86]
	     */, calcst;
} reactn_;

#define reactn_1 reactn_

/* Table of constant values */

static integer c__5 = 5;
static integer c__1 = 1;
static doublereal c_b63 = 4.;
static logical c_true = TRUE_;
static logical c_false = FALSE_;

/* Subroutine */ int react1_(doublereal *escf)
{
    /* Initialized data */

    static integer irot[6]	/* was [2][3] */ = { 1,2,1,3,2,3 };

    /* System generated locals */
    integer i__1, i__2;
    doublereal d__1, d__2, d__3;
    cllist cl__1;
    static doublereal equiv_0[258];

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

    /* Local variables */
    static doublereal dell, eold, gold, xold[258], swap, summ, sumx, sumy, 
	    sumz, time0;
    static integer idum1[86], idum2[86], idum3[86];
    static doublereal time1, time2, step0;
    static integer i, j, k;
    extern doublereal reada_(char *, integer *, ftnlen);
    static integer l;
    static doublereal x;
    extern /* Subroutine */ int flepo_(doublereal *, integer *, doublereal *);
#define coord (equiv_0)
    static doublereal grold[258];
    extern /* Subroutine */ int geout_(void);
    static doublereal const__;
    static integer iloop;
    static doublereal c1;
    extern /* Subroutine */ int write_(doublereal *, doublereal *);
    static doublereal c2, funct1, ca;
    static integer numat2;
    static doublereal sa, p1stor[23220], p2stor[23220], p3stor[23220];
    static integer ir, jr, linear;
    extern /* Subroutine */ int getgeo_(integer *, integer *, doublereal *, 
	    integer *, integer *, integer *, integer *, doublereal *, integer 
	    *, logical *);
    extern doublereal second_(void);
    static logical gradnt, finish;
    extern /* Subroutine */ int compfg_(doublereal *, logical *, doublereal *,
	     logical *, doublereal *, logical *);
#define idummy ((integer *)equiv_0)
    extern /* Subroutine */ int gmetry_(doublereal *, doublereal *);
    static integer maxstp;
    static doublereal stepmx, xstore[258];
    extern /* Subroutine */ int symtry_(void);
    static logical gok[2];
    static doublereal one;
    extern doublereal dot_(doublereal *, doublereal *, integer *);
    static logical int__;
    static doublereal sum;
    static logical xyz;

    /* Fortran I/O blocks */
    static cilist io___17 = { 0, 6, 0, "(10X,'ERRORS DETECTED IN CONNECTIVIT"
	    "Y')", 0 };
    static cilist io___18 = { 0, 6, 0, "(A,I3,A,I3,A,I3,A)", 0 };
    static cilist io___19 = { 0, 6, 0, "(10X,'ERRORS DETECTED IN CONNECTIVIT"
	    "Y')", 0 };
    static cilist io___20 = { 0, 6, 0, "(A,I3,A,I3,A,I3,A)", 0 };
    static cilist io___21 = { 0, 6, 0, "(10X,'ERRORS DETECTED IN CONNECTIVIT"
	    "Y')", 0 };
    static cilist io___22 = { 0, 6, 0, "(A,I3,A,I3,A,I3,A)", 0 };
    static cilist io___23 = { 0, 6, 0, "(10X,A)", 0 };
    static cilist io___29 = { 0, 6, 0, "(//10X,' NUMBER OF ATOMS IN SECOND S"
	    "YSTEM IS ',     'INCORRECT',/)", 0 };
    static cilist io___30 = { 0, 6, 0, "(' NUMBER OF ATOMS IN FIRST  SYSTEM "
	    "=',I4)", 0 };
    static cilist io___31 = { 0, 6, 0, "(' NUMBER OF ATOMS IN SECOND SYSTEM "
	    "=',I4)", 0 };
    static cilist io___32 = { 0, 6, 0, "(//10X,' GEOMETRY OF SECOND SYSTEM',"
	    "/)", 0 };
    static cilist io___36 = { 0, 6, 0, "(//,'  CARTESIAN GEOMETRY OF FIRST S"
	    "YSTEM',//)", 0 };
    static cilist io___37 = { 0, 6, 0, "(3F14.5)", 0 };
    static cilist io___45 = { 0, 6, 0, "(//,'  CARTESIAN GEOMETRY OF SECOND "
	    "SYSTEM',//)", 0 };
    static cilist io___46 = { 0, 6, 0, "(3F14.5)", 0 };
    static cilist io___47 = { 0, 6, 0, "(//,'   \"DISTANCE\":',F13.6)", 0 };
    static cilist io___48 = { 0, 6, 0, "(//,'  REACTION COORDINATE VECTOR',/"
	    "/)", 0 };
    static cilist io___49 = { 0, 6, 0, "(3F14.5)", 0 };
    static cilist io___50 = { 0, 6, 0, "(///10X,'THERE ARE NO VARIABLES IN T"
	    "HE SADDLE',     ' CALCULATION!')", 0 };
    static cilist io___61 = { 0, 6, 0, "(' TIME=',F9.2)", 0 };
    static cilist io___62 = { 0, 6, 0, "(' CURRENT BAR, STEPMX, GNORM',3F12."
	    "7)", 0 };
    static cilist io___67 = { 0, 6, 0, "(//10X,'FOR POINT',I3)", 0 };
    static cilist io___68 = { 0, 6, 0, "(' DISTANCE A - B  ',F12.6)", 0 };
    static cilist io___70 = { 0, 6, 0, "('  ACTUAL GRADIENTS OF THIS POINT')",
	     0 };
    static cilist io___71 = { 0, 6, 0, "(8F10.4)", 0 };
    static cilist io___72 = { 0, 6, 0, "(' HEAT            ',F12.6)", 0 };
    static cilist io___73 = { 0, 6, 0, "(' GRADIENT NORM   ',F12.6)", 0 };
    static cilist io___74 = { 0, 6, 0, "(' DIRECTION COSINE',F12.6)", 0 };
    static cilist io___76 = { 0, 6, 0, "(//10X,' BOTH SYSTEMS ARE ON THE SAM"
	    "E SIDE OF THE ','TRANSITION STATE -',/10X,' GEOMETRIES OF THE SY"
	    "STEMS', ' ON EACH SIDE OF THE T.S. ARE AS FOLLOWS')", 0 };
    static cilist io___77 = { 0, 6, 0, "(//10X,' GEOMETRY ON ONE SIDE OF THE"
	    " TRANSITION',' STATE')", 0 };
    static cilist io___78 = { 0, 6, 0, "('  REACTANTS AND PRODUCTS SWAPPED A"
	    "ROUND')", 0 };
    static cilist io___79 = { 0, 6, 0, "(' AT END OF REACTION')", 0 };
    static cilist io___82 = { 0, 6, 0, "(' BEST ESTIMATE GEOMETRY OF THE TRA"
	    "NSITION STATE')", 0 };
    static cilist io___83 = { 0, 6, 0, "(//10X,' C1=',F8.3,'C2=',F8.3)", 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 */

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

/*  REACT1 DETERMINES THE TRANSITION STATE OF A CHEMICAL REACTION. */

/*   REACT WORKS BY USING TWO SYSTEMS SIMULTANEOUSLY, THE HEATS OF */
/*   FORMATION OF BOTH ARE CALCULATED, THEN THE MORE STABLE ONE */
/*   IS MOVED IN THE DIRECTION OF THE OTHER. AFTER A STEP THE */
/*   ENERGIES ARE COMPARED, AND THE NOW LOWER-ENERGY FORM IS MOVED */
/*   IN THE DIRECTION OF THE HIGHER-ENERGY FORM. THIS IS REPEATED */
/*   UNTIL THE SADDLE POINT IS REACHED. */

/*   IF ONE FORM IS MOVED 3 TIMES IN SUCCESSION, THEN THE HIGHER ENERGY */
/*   FORM IS RE-OPTIMIZED WITHOUT SHORTENING THE DISTANCE BETWEEN THE TWO 
*/
/*   FORMS. THIS REDUCES THE CHANCE OF BEING CAUGHT ON THE SIDE OF A */
/*   TRANSITION STATE. */

/* ***********************************************************************
 */
    gold = 0.;
    linear = 0;
    gok[0] = FALSE_;
    gok[1] = FALSE_;
    xyz = i_indx(keywrd_1.keywrd, " XYZ", 80L, 4L) != 0;
    gradnt = i_indx(keywrd_1.keywrd, "GRAD", 80L, 4L) != 0;
    i = i_indx(keywrd_1.keywrd, " BAR", 80L, 4L);
    stepmx = .15;
    if (i != 0) {
	stepmx = reada_(keywrd_1.keywrd, &i, 80L);
    }
    maxstp = 1000;

/*    READ IN THE SECOND GEOMETRY. */

    if (xyz) {
	getgeo_(&c__5, geokst_1.labels, reactn_1.geoa, geovar_1.loc, 
		geokst_1.na, geokst_1.nb, geokst_1.nc, istope_1.ams, &
		geokst_1.natoms, &int__);
    } else {
	getgeo_(&c__5, idum1, reactn_1.geoa, idummy, idum1, idum2, idum3, 
		istope_1.ams, &geokst_1.natoms, &int__);

/*  IF INTERNAL COORDINATES ARE TO BE USED, CHECK THE CONNECTIVITY */

	l = 0;
	i__1 = geokst_1.natoms;
	for (i = 1; i <= i__1; ++i) {
	    if (idum1[i - 1] != geokst_1.na[i - 1]) {
		++l;
		if (l == 1) {
		    s_wsfe(&io___17);
		    e_wsfe();
		}
		s_wsfe(&io___18);
		do_fio(&c__1, " FOR ATOM", 9L);
		do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
		do_fio(&c__1, " THE BOND LABELS ARE DIFFERENT:      ", 37L);
		do_fio(&c__1, (char *)&idum1[i - 1], (ftnlen)sizeof(integer));
		do_fio(&c__1, " AND", 4L);
		do_fio(&c__1, (char *)&geokst_1.na[i - 1], (ftnlen)sizeof(
			integer));
		e_wsfe();
	    }
	    if (idum2[i - 1] != geokst_1.nb[i - 1]) {
		++l;
		if (l == 1) {
		    s_wsfe(&io___19);
		    e_wsfe();
		}
		s_wsfe(&io___20);
		do_fio(&c__1, " FOR ATOM", 9L);
		do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
		do_fio(&c__1, " THE BOND ANGLE LABELS ARE DIFFERENT:", 37L);
		do_fio(&c__1, (char *)&idum2[i - 1], (ftnlen)sizeof(integer));
		do_fio(&c__1, " AND", 4L);
		do_fio(&c__1, (char *)&geokst_1.nb[i - 1], (ftnlen)sizeof(
			integer));
		e_wsfe();
	    }
	    if (idum3[i - 1] != geokst_1.nc[i - 1]) {
		++l;
		if (l == 1) {
		    s_wsfe(&io___21);
		    e_wsfe();
		}
		s_wsfe(&io___22);
		do_fio(&c__1, " FOR ATOM", 9L);
		do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
		do_fio(&c__1, " THE DIHEDRAL LABELS ARE DIFFERENT:  ", 37L);
		do_fio(&c__1, (char *)&idum3[i - 1], (ftnlen)sizeof(integer));
		do_fio(&c__1, " AND", 4L);
		do_fio(&c__1, (char *)&geokst_1.nc[i - 1], (ftnlen)sizeof(
			integer));
		e_wsfe();
	    }
/* L8: */
	}
	if (l != 0) {
	    s_wsfe(&io___23);
	    do_fio(&c__1, " CORRECT BEFORE RESUBMISSION", 28L);
	    e_wsfe();
	}
	if (l != 0) {
	    s_stop("", 0L);
	}
    }
    cl__1.cerr = 0;
    cl__1.cunit = 5;
    cl__1.csta = 0;
    f_clos(&cl__1);
    time0 = second_();

/*  SWAP FIRST AND SECOND GEOMETRIES AROUND */
/*  SO THAT GEOUT CAN OUTPUT DATA ON SECOND GEOMETRY. */

    numat2 = 0;
    i__1 = geokst_1.natoms;
    for (i = 1; i <= i__1; ++i) {
	if (geokst_1.labels[i - 1] != 99) {
	    ++numat2;
	}
	const__ = 1.;
	for (j = 1; j <= 3; ++j) {
	    x = reactn_1.geoa[j + i * 3 - 4] * const__;
	    const__ = .0174532925;
	    reactn_1.geoa[j + i * 3 - 4] = geom_1.geo[j + i * 3 - 4];
	    geom_1.geo[j + i * 3 - 4] = x;
/* L10: */
	}
    }
    if (numat2 != molkst_1.numat) {
	s_wsfe(&io___29);
	e_wsfe();
	s_wsfe(&io___30);
	do_fio(&c__1, (char *)&molkst_1.numat, (ftnlen)sizeof(integer));
	e_wsfe();
	s_wsfe(&io___31);
	do_fio(&c__1, (char *)&numat2, (ftnlen)sizeof(integer));
	e_wsfe();
	goto L270;
    }
    s_wsfe(&io___32);
    e_wsfe();
    if (geosym_1.ndep != 0) {
	symtry_();
    }
    geout_();

/*     CONVERT TO CARTESIAN, IF NECESSARY */

    if (xyz) {
	gmetry_(geom_1.geo, coord);
	sumx = 0.;
	sumy = 0.;
	sumz = 0.;
	i__1 = molkst_1.numat;
	for (j = 1; j <= i__1; ++j) {
	    sumx += coord[j * 3 - 3];
	    sumy += coord[j * 3 - 2];
/* L20: */
	    sumz += coord[j * 3 - 1];
	}
	sumx /= molkst_1.numat;
	sumy /= molkst_1.numat;
	sumz /= molkst_1.numat;
	i__1 = molkst_1.numat;
	for (j = 1; j <= i__1; ++j) {
	    geom_1.geo[j * 3 - 3] = coord[j * 3 - 3] - sumx;
	    geom_1.geo[j * 3 - 2] = coord[j * 3 - 2] - sumy;
/* L30: */
	    geom_1.geo[j * 3 - 1] = coord[j * 3 - 1] - sumz;
	}
	s_wsfe(&io___36);
	e_wsfe();
	s_wsfe(&io___37);
	i__1 = molkst_1.numat;
	for (i = 1; i <= i__1; ++i) {
	    for (j = 1; j <= 3; ++j) {
		do_fio(&c__1, (char *)&geom_1.geo[j + i * 3 - 4], (ftnlen)
			sizeof(doublereal));
	    }
	}
	e_wsfe();
	sumx = 0.;
	sumy = 0.;
	sumz = 0.;
	i__1 = molkst_1.numat;
	for (j = 1; j <= i__1; ++j) {
	    sumx += reactn_1.geoa[j * 3 - 3];
	    sumy += reactn_1.geoa[j * 3 - 2];
/* L40: */
	    sumz += reactn_1.geoa[j * 3 - 1];
	}
	sum = 0.;
	sumx /= molkst_1.numat;
	sumy /= molkst_1.numat;
	sumz /= molkst_1.numat;
	i__1 = molkst_1.numat;
	for (j = 1; j <= i__1; ++j) {
	    reactn_1.geoa[j * 3 - 3] -= sumx;
	    reactn_1.geoa[j * 3 - 2] -= sumy;
	    reactn_1.geoa[j * 3 - 1] -= sumz;
/* Computing 2nd power */
	    d__1 = geom_1.geo[j * 3 - 3] - reactn_1.geoa[j * 3 - 3];
/* Computing 2nd power */
	    d__2 = geom_1.geo[j * 3 - 2] - reactn_1.geoa[j * 3 - 2];
/* Computing 2nd power */
	    d__3 = geom_1.geo[j * 3 - 1] - reactn_1.geoa[j * 3 - 1];
	    sum = sum + d__1 * d__1 + d__2 * d__2 + d__3 * d__3;
/* L50: */
	}
	for (l = 3; l >= 1; --l) {

/*     DOCKING IS DONE IN STEPS OF 16, 4, AND 1 DEGREES AT A TIME.
 */

	    i__1 = l - 1;
	    ca = cos(pow_di(&c_b63, &i__1) * .01745329);
/* Computing 2nd power */
	    d__2 = ca;
	    sa = sqrt((d__1 = 1. - d__2 * d__2, abs(d__1)));
	    for (j = 1; j <= 3; ++j) {
		ir = irot[(j << 1) - 2];
		jr = irot[(j << 1) - 1];
		for (i = 1; i <= 10; ++i) {
		    summ = 0.;
		    i__1 = molkst_1.numat;
		    for (k = 1; k <= i__1; ++k) {
			x = ca * reactn_1.geoa[ir + k * 3 - 4] + sa * 
				reactn_1.geoa[jr + k * 3 - 4];
			reactn_1.geoa[jr + k * 3 - 4] = -sa * reactn_1.geoa[
				ir + k * 3 - 4] + ca * reactn_1.geoa[jr + k * 
				3 - 4];
			reactn_1.geoa[ir + k * 3 - 4] = x;
/* Computing 2nd power */
			d__1 = geom_1.geo[k * 3 - 3] - reactn_1.geoa[k * 3 - 
				3];
/* Computing 2nd power */
			d__2 = geom_1.geo[k * 3 - 2] - reactn_1.geoa[k * 3 - 
				2];
/* Computing 2nd power */
			d__3 = geom_1.geo[k * 3 - 1] - reactn_1.geoa[k * 3 - 
				1];
			summ = summ + d__1 * d__1 + d__2 * d__2 + d__3 * d__3;
/* L60: */
		    }
		    if (summ > sum) {
			if (i > 1) {
			    sa = -sa;
			    i__1 = molkst_1.numat;
			    for (k = 1; k <= i__1; ++k) {
				x = ca * reactn_1.geoa[ir + k * 3 - 4] + sa * 
					reactn_1.geoa[jr + k * 3 - 4];
				reactn_1.geoa[jr + k * 3 - 4] = -sa * 
					reactn_1.geoa[ir + k * 3 - 4] + ca * 
					reactn_1.geoa[jr + k * 3 - 4];
				reactn_1.geoa[ir + k * 3 - 4] = x;
/* L70: */
			    }
			    goto L90;
			}
			sa = -sa;
		    }
/* L80: */
		    sum = summ;
		}
L90:
		;
	    }
/* L100: */
	}
	s_wsfe(&io___45);
	e_wsfe();
	s_wsfe(&io___46);
	i__1 = molkst_1.numat;
	for (i = 1; i <= i__1; ++i) {
	    for (j = 1; j <= 3; ++j) {
		do_fio(&c__1, (char *)&reactn_1.geoa[j + i * 3 - 4], (ftnlen)
			sizeof(doublereal));
	    }
	}
	e_wsfe();
	s_wsfe(&io___47);
	do_fio(&c__1, (char *)&sum, (ftnlen)sizeof(doublereal));
	e_wsfe();
	s_wsfe(&io___48);
	e_wsfe();
	s_wsfe(&io___49);
	i__1 = molkst_1.numat;
	for (i = 1; i <= i__1; ++i) {
	    for (j = 1; j <= 3; ++j) {
		d__1 = reactn_1.geoa[j + i * 3 - 4] - geom_1.geo[j + i * 3 - 
			4];
		do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof(doublereal));
	    }
	}
	e_wsfe();
	geokst_1.na[0] = 99;
	j = 0;
	geovar_1.nvar = 0;
	i__1 = geokst_1.natoms;
	for (i = 1; i <= i__1; ++i) {
	    if (geokst_1.labels[i - 1] != 99) {
		++j;
		for (k = 1; k <= 3; ++k) {
		    ++geovar_1.nvar;
		    geovar_1.loc[(geovar_1.nvar << 1) - 1] = k;
/* L110: */
		    geovar_1.loc[(geovar_1.nvar << 1) - 2] = j;
		}
		geokst_1.labels[j - 1] = geokst_1.labels[i - 1];
	    }
/* L120: */
	}
	geokst_1.natoms = molkst_1.numat;
    }

/*   XPARAM HOLDS THE VARIABLE PARAMETERS FOR GEOMETRY IN GEO */
/*   XOLD   HOLDS THE VARIABLE PARAMETERS FOR GEOMETRY IN GEOA */

    if (geovar_1.nvar == 0) {
	s_wsfe(&io___50);
	e_wsfe();
	s_stop("", 0L);
    }
    sum = 0.;
    i__1 = geovar_1.nvar;
    for (i = 1; i <= i__1; ++i) {
	grold[i - 1] = 1.;
	geovar_1.xparam[i - 1] = geom_1.geo[geovar_1.loc[(i << 1) - 1] + 
		geovar_1.loc[(i << 1) - 2] * 3 - 4];
	xold[i - 1] = reactn_1.geoa[geovar_1.loc[(i << 1) - 1] + geovar_1.loc[
		(i << 1) - 2] * 3 - 4];
/* L130: */
/* Computing 2nd power */
	d__1 = geovar_1.xparam[i - 1] - xold[i - 1];
	sum += d__1 * d__1;
    }
    step0 = sqrt(sum);
    one = 1.;
    dell = .1;
    eold = -2e3;
    time1 = second_();
    swap = 0.;
    i__1 = maxstp;
    for (iloop = 1; iloop <= i__1; ++iloop) {
	time2 = second_();
	s_wsfe(&io___61);
	d__1 = time2 - time1;
	do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof(doublereal));
	e_wsfe();
	time1 = time2;

/*   THIS METHOD OF CALCULATING 'STEP' IS QUITE ARBITARY, AND NEEDS */
/*   TO BE IMPROVED BY INTELLIGENT GUESSWORK! */

	if (gradnt_1.gnorm < .001) {
	    gradnt_1.gnorm = .001;
	}
	s_wsfe(&io___62);
	do_fio(&c__1, (char *)&step0, (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&stepmx, (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&gradnt_1.gnorm, (ftnlen)sizeof(doublereal));
	e_wsfe();
/* Computing MIN */
	d__1 = min(swap,.5), d__2 = 6. / gradnt_1.gnorm, d__1 = min(d__1,d__2)
		, d__1 = min(d__1,dell), d__2 = stepmx * step0 + .005;
	reactn_1.step = min(d__1,d__2);
	swap += 1.;
	dell += .1f;
	step0 -= reactn_1.step;
	if (step0 < .01) {
	    goto L240;
	}
	reactn_1.step = step0;
	i__2 = geovar_1.nvar;
	for (i = 1; i <= i__2; ++i) {
/* L150: */
	    xstore[i - 1] = geovar_1.xparam[i - 1];
	}
	flepo_(geovar_1.xparam, &geovar_1.nvar, escf);
	if (linear == 0) {
	    linear = molkst_1.norbs * (molkst_1.norbs + 1) / 2;
	    i__2 = linear;
	    for (i = 1; i <= i__2; ++i) {
		p1stor[i - 1] = densty_1.p[i - 1];
		p2stor[i - 1] = densty_1.pa[i - 1];
/* L140: */
		p3stor[i - 1] = densty_1.pb[i - 1];
	    }
	}
	i__2 = geovar_1.nvar;
	for (i = 1; i <= i__2; ++i) {
/* L160: */
	    geovar_1.xparam[i - 1] = geom_1.geo[geovar_1.loc[(i << 1) - 1] + 
		    geovar_1.loc[(i << 1) - 2] * 3 - 4];
	}
	s_wsfe(&io___67);
	do_fio(&c__1, (char *)&iloop, (ftnlen)sizeof(integer));
	e_wsfe();
	s_wsfe(&io___68);
	do_fio(&c__1, (char *)&reactn_1.step, (ftnlen)sizeof(doublereal));
	e_wsfe();

/*   NOW TO CALCULATE THE "CORRECT" GRADIENTS, SWITCH OFF 'STEP'. */

	reactn_1.step = 0.;
	i__2 = geovar_1.nvar;
	for (i = 1; i <= i__2; ++i) {
/* L170: */
	    gradnt_1.grad[i - 1] = grold[i - 1];
	}
	compfg_(geovar_1.xparam, &c_true, &funct1, &c_false, gradnt_1.grad, &
		c_true);
	i__2 = geovar_1.nvar;
	for (i = 1; i <= i__2; ++i) {
/* L180: */
	    grold[i - 1] = gradnt_1.grad[i - 1];
	}
	if (gradnt) {
	    s_wsfe(&io___70);
	    e_wsfe();
	    s_wsfe(&io___71);
	    i__2 = geovar_1.nvar;
	    for (i = 1; i <= i__2; ++i) {
		do_fio(&c__1, (char *)&gradnt_1.grad[i - 1], (ftnlen)sizeof(
			doublereal));
	    }
	    e_wsfe();
	}
	s_wsfe(&io___72);
	do_fio(&c__1, (char *)&funct1, (ftnlen)sizeof(doublereal));
	e_wsfe();
	gradnt_1.gnorm = sqrt(dot_(gradnt_1.grad, gradnt_1.grad, &
		geovar_1.nvar));
	s_wsfe(&io___73);
	do_fio(&c__1, (char *)&gradnt_1.gnorm, (ftnlen)sizeof(doublereal));
	e_wsfe();
	gravec_1.cosine *= one;
	s_wsfe(&io___74);
	do_fio(&c__1, (char *)&gravec_1.cosine, (ftnlen)sizeof(doublereal));
	e_wsfe();
	geout_();
	if (swap > 2.9 || iloop > 3 && gravec_1.cosine < 0. || *escf > eold) {
	    if (swap > 2.9) {
		swap = 0.;
	    } else {
		swap = .5;
	    }

/*   SWAP REACTANT AND PRODUCT AROUND */

	    finish = gok[0] && gok[1] && gravec_1.cosine < 0.;
	    if (finish) {
		s_wsfe(&io___76);
		e_wsfe();
		i__2 = geovar_1.nvar;
		for (i = 1; i <= i__2; ++i) {
/* L190: */
		    geovar_1.xparam[i - 1] = xstore[i - 1];
		}
		compfg_(geovar_1.xparam, &c_true, &funct1, &c_true, 
			gradnt_1.grad, &c_true);
		s_wsfe(&io___77);
		e_wsfe();
		write_(&time0, &funct1);
	    }
	    s_wsfe(&io___78);
	    e_wsfe();
	    one = -1.;
	    eold = *escf;
	    sum = gold;
	    gold = gradnt_1.gnorm;
	    i = (integer) (one * .5f + 1.7f);
	    if (gradnt_1.gnorm > 10.) {
		gok[i - 1] = TRUE_;
	    }
	    gradnt_1.gnorm = sum;
	    i__2 = molkst_1.numat;
	    for (i = 1; i <= i__2; ++i) {
		for (j = 1; j <= 3; ++j) {
		    x = geom_1.geo[j + i * 3 - 4];
		    geom_1.geo[j + i * 3 - 4] = reactn_1.geoa[j + i * 3 - 4];
/* L200: */
		    reactn_1.geoa[j + i * 3 - 4] = x;
		}
	    }
	    i__2 = geovar_1.nvar;
	    for (i = 1; i <= i__2; ++i) {
		x = xold[i - 1];
		xold[i - 1] = geovar_1.xparam[i - 1];
/* L210: */
		geovar_1.xparam[i - 1] = x;
	    }


/*    SWAP AROUND THE DENSITY MATRICES. */

	    i__2 = linear;
	    for (i = 1; i <= i__2; ++i) {
		x = p1stor[i - 1];
		p1stor[i - 1] = densty_1.p[i - 1];
		densty_1.p[i - 1] = x;
		x = p2stor[i - 1];
		p2stor[i - 1] = densty_1.pa[i - 1];
		densty_1.pa[i - 1] = x;
		x = p3stor[i - 1];
		p3stor[i - 1] = densty_1.pb[i - 1];
		densty_1.pb[i - 1] = x;
/* L220: */
	    }
	    if (finish) {
		goto L240;
	    }
	} else {
	    one = 1.;
	}
/* L230: */
    }
L240:
    s_wsfe(&io___79);
    e_wsfe();
    gold = sqrt(dot_(gradnt_1.grad, gradnt_1.grad, &geovar_1.nvar));
    compfg_(geovar_1.xparam, &c_true, &funct1, &c_true, gradnt_1.grad, &
	    c_true);
    gradnt_1.gnorm = sqrt(dot_(gradnt_1.grad, gradnt_1.grad, &geovar_1.nvar));
    i__1 = geovar_1.nvar;
    for (i = 1; i <= i__1; ++i) {
/* L250: */
	grold[i - 1] = geovar_1.xparam[i - 1];
    }
    write_(&time0, &funct1);

/* THE GEOMETRIES HAVE (A) BEEN OPTIMIZED CORRECTLY, OR */
/*                     (B) BOTH ENDED UP ON THE SAME SIDE OF THE T.S. */

/*  TRANSITION STATE LIES BETWEEN THE TWO GEOMETRIES */

    c1 = gold / (gold + gradnt_1.gnorm);
    c2 = 1. - c1;
    s_wsfe(&io___82);
    e_wsfe();
    s_wsfe(&io___83);
    do_fio(&c__1, (char *)&c1, (ftnlen)sizeof(doublereal));
    do_fio(&c__1, (char *)&c2, (ftnlen)sizeof(doublereal));
    e_wsfe();
    i__1 = geovar_1.nvar;
    for (i = 1; i <= i__1; ++i) {
/* L260: */
	geovar_1.xparam[i - 1] = c1 * grold[i - 1] + c2 * xold[i - 1];
    }
    reactn_1.step = 0.;
    compfg_(geovar_1.xparam, &c_true, &funct1, &c_true, gradnt_1.grad, &
	    c_true);
    write_(&time0, &funct1);
L270:
    s_stop("", 0L);
    return 0;
} /* react1_ */

#undef idummy
#undef coord