/* flepo.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 cosine; } gravec_; #define gravec_1 gravec_ struct { integer last; } last_; #define last_1 last_ struct { integer latom, lparam; doublereal react[200]; } path_; #define path_1 path_ struct { doublereal grad[258], gnorm; } gradnt_; #define gradnt_1 gradnt_ struct { integer iflepo, iscf; } mesage_; #define mesage_1 mesage_ struct { integer nscf; } numscf_; #define numscf_1 numscf_ struct { doublereal time0; } time_; #define time_1 time_ struct { doublereal hesinv[33411]; } fmatrx_; #define fmatrx_1 fmatrx_ struct { integer numcal; } numcal_; #define numcal_1 numcal_ /* Table of constant values */ static integer c__1 = 1; static logical c_true = TRUE_; static logical c_false = FALSE_; /* Subroutine */ int flepo_(doublereal *xparam, integer *nvar, doublereal * funct1) { /* Initialized data */ static integer icalcn = 0; /* Format strings */ static char fmt_130[] = "(\002 FUNCTION VALUE=\002,f13.7,\002 WILL NOT " "BE REPLACED BY VALUE=\002,f13.7,/10x,\002CALCULATED BY RESTART P" "ROCEDURE\002,/)"; static char fmt_140[] = "(\002 FUNCTION VALUE=\002,f13.7,\002 IS BEING R" "EPLACED BY VALUE=\002,f13.7,/10x,\002 FOUND IN RESTART PROCEDUR" "E\002,/,6x,\002THE CORRESPONDING X VALUES AND GRADIENTS ARE ALSO" " BEING REPLACED\002,/)"; static char fmt_250[] = "(//,5x,\002SINCE COS=\002,f9.3,5x,\002THE PROGR" "AM WILL GO TO RE\002,\002START SECTION\002,/)"; static char fmt_270[] = "(\002 THE CURRENT POINT IS ...\002)"; static char fmt_280[] = "(\002 \002,3x,\002I\002,9x,i3,9(8x,i3))"; static char fmt_290[] = "(\002 \002,1x,\002XPARAM(I)\002,1x,f9.4,2x,9(f9" ".4,2x))"; static char fmt_300[] = "(\002 \002,1x,\002GRAD (I)\002,f10.4,1x,9(f10." "4,1x))"; static char fmt_310[] = "(\002 \002,1x,\002PVECT (I)\002,2x,f10.6,1x,9(f" "10.6,1x))"; static char fmt_340[] = "(//,10x,\002HERBERTS TEST SATISFIED - GEOMETRY " "OPTIMIZED\002)"; static char fmt_370[] = "(\002 \002,///,20x,\002NO POINT LOWER IN ENER" "GY \002,\002THAN THE STARTING POINT\002,/,20x,\002COULD BE FOUND " "\002,\002IN THE LINE MINIMIZATION\002)"; static char fmt_390[] = "(\002 \002,//,20x,\002SINCE COS WAS JUST RESET," "THE SEARCH\002,\002 IS BEING ENDED\002)"; static char fmt_400[] = "(\002 \002,20x,\002COS WILL BE RESET AND ANOTHE" "R \002,\002ATTEMPT MADE\002)"; static char fmt_420[] = "(\0020LINMIN FAILED AT CYCLE\002,i5/,\0020\002)"; static char fmt_440[] = "(/,\002 NUMBER OF COUNTS =\002,i6," "\002 COS =\002,f11.4,/,\002 ABSOLUTE CHANGE IN X " " =\002,f13.6,\002 ALPHA =\002,f11.4,/,\002 PREDICTED CHANGE " "IN F = \002,g11.4,\002 ACTUAL = \002,g11.4,/,\002 GRADIE" "NT NORM = \002,g11.4,//)"; static char fmt_460[] = "(\0020TEST ON X SATISFIED\002)"; static char fmt_470[] = "(\002 HEAT OF FORMATION TEST SATISFIED\002)"; static char fmt_480[] = "(\0020TEST ON GRADIENT SATISFIED\002)"; static char fmt_510[] = "(20x,\002HOWEVER, A COMPONENT OF GRADIENT IS" " \002,\002LARGER THAN\002,f6.2,/)"; static char fmt_520[] = "(7x,\002 THERE HAVE BEEN\002,i2,\002 ATTEMPTS T" "O REDUCE THE \002,\002 GRADIENT.\002,/7x,\002 DURING THESE ATTEM" "PTS THE ENERGY DROPPED\002,\002 BY LESS THAN\002,f7.4,\002 KCAL/" "MOLE\002,/10x,\002 FURTHER CALCULATION IS NOT JUSTIFIED AT THIS " "POINT.\002)"; static char fmt_540[] = "(\002 PETERS TEST SATISFIED \002)"; static char fmt_570[] = "(\002 RESTART FILE WRITTEN, TIME LEFT:\002,f9" ".1,\002 GRAD.:\002,f10.3,\002 HEAT:\002,g13.7)"; static char fmt_580[] = "(\002 CYCLE:\002,i4,\002 TIME:\002,f7.2,\002 TI" "ME LEFT:\002,f9.1,\002 GRAD.:\002,f10.3,\002 HEAT:\002,g13.7)"; static char fmt_590[] = "(20x,\002THERE IS NOT ENOUGH TIME FOR ANOTHER C" "YCLE\002,/,30x,\002NOW GOING TO FINAL\002)"; /* System generated locals */ integer i__1, i__2; doublereal d__1, d__2, d__3; static integer equiv_3[9]; static doublereal equiv_11[9]; /* Builtin functions */ integer i_indx(char *, char *, ftnlen, ftnlen), s_cmp(char *, char *, ftnlen, ftnlen); double sqrt(doublereal); integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void); double d_sign(doublereal *, doublereal *); /* Subroutine */ int s_stop(char *, ftnlen); /* Local variables */ static doublereal einc, dell; #define mdfp (equiv_3) #define jcyc (equiv_3) static doublereal tdel; static logical time; #define xdfp (equiv_11) #define drop (equiv_11 + 3) static integer loop; static doublereal dott; static integer iinc1, iinc2, itry1, i, j, itry2; extern doublereal reada_(char *, integer *, ftnlen); static integer k; static doublereal p, s; #define alpha (equiv_11) static doublereal y; static integer ihdim; static doublereal sfact; #define frepf (equiv_11 + 5) static logical geook; static doublereal glast[258], bsmvf, tleft, funct, pvect[258]; static logical reset; #define cycmx (equiv_11 + 6) extern /* Subroutine */ int geout_(void); static doublereal const__, tlast, pmste, tdump, tolrg, xlast[258]; static logical print; #define pnorm (equiv_11 + 2) static doublereal smval, totim; #define jnrst (equiv_3 + 1) static doublereal rootv, funct2; static char ch[1]; static doublereal gd[258], cncadd, gg[258]; static integer ii, ij; static doublereal tf, xd[258]; static logical saddle; static doublereal deltag, delhof, pt, xn, absmin; extern doublereal second_(void); extern /* Subroutine */ int compfg_(doublereal *, logical *, doublereal *, logical *, doublereal *, logical *); static doublereal tx; extern /* Subroutine */ int dfpsav_(doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, integer *, doublereal *) ; static logical resfil; static doublereal gdnorm, tolerg, tolerf; #define totime (equiv_11 + 7) static integer irepet; static doublereal pnlast; #define ncount (equiv_3 + 2) extern /* Subroutine */ int linmin_(doublereal *, doublereal *, doublereal *, integer *, doublereal *, logical *, logical *); static logical minprt; static doublereal tcycle, tolerx; #define lnstop (equiv_3 + 3) static logical restrt; static doublereal tx1, tx2, ggd; #define del (equiv_11 + 4) static logical dfp, okc, okf; #define cos__ (equiv_11 + 1) extern doublereal dot_(doublereal *, doublereal *, integer *); static doublereal tim, rst, yhy, pty; static integer igg1, nto6; /* Fortran I/O blocks */ static cilist io___49 = { 0, 6, 0, "(/,A)", 0 }; static cilist io___59 = { 0, 6, 0, "(' GRADIENT CRITERION IN FLEPO =',F1" "2.5)", 0 }; static cilist io___62 = { 0, 6, 0, "(//10X,'TOTAL TIME USED SO FAR:', " " F13.2,' SECONDS')", 0 }; static cilist io___67 = { 0, 6, 0, "(' GRADIENTS OF OLD GEOMETRY, GNORM=" "',F13.6)", 0 }; static cilist io___68 = { 0, 6, 0, "(6F12.6)", 0 }; static cilist io___70 = { 0, 6, 0, "(' GRADIENTS OF NEW GEOMETRY, GNORM=" "',F13.6)", 0 }; static cilist io___71 = { 0, 6, 0, "(6F12.6)", 0 }; static cilist io___72 = { 0, 6, 0, "(///20X,'CALCULATION ABANDONED AT TH" "IS POINT!')", 0 }; static cilist io___73 = { 0, 6, 0, "(//10X,' SMALL CHANGES IN INTERNAL C" "OORDINATES ARE',/10X,' CAUSING A LARGE CHANGE IN THE DISTANCE BE" "TWEEN')", 0 }; static cilist io___74 = { 0, 6, 0, "(10X,' CHEMICALLY-BOUND ATOMS. THE O" "PTIMISATION',/ 10X,' PROCEDURE WOULD LIKELY PRODUCE INCORRECT R" "ESULTS')", 0 }; static cilist io___78 = { 0, 6, 0, fmt_130, 0 }; static cilist io___79 = { 0, 6, 0, fmt_140, 0 }; static cilist io___93 = { 0, 6, 0, fmt_250, 0 }; static cilist io___94 = { 0, 6, 0, fmt_270, 0 }; static cilist io___97 = { 0, 6, 0, "(/)", 0 }; static cilist io___99 = { 0, 6, 0, fmt_280, 0 }; static cilist io___100 = { 0, 6, 0, fmt_290, 0 }; static cilist io___101 = { 0, 6, 0, fmt_300, 0 }; static cilist io___102 = { 0, 6, 0, fmt_310, 0 }; static cilist io___103 = { 0, 6, 0, fmt_340, 0 }; static cilist io___107 = { 0, 6, 0, fmt_370, 0 }; static cilist io___108 = { 0, 6, 0, fmt_390, 0 }; static cilist io___109 = { 0, 6, 0, fmt_400, 0 }; static cilist io___110 = { 0, 6, 0, fmt_420, 0 }; static cilist io___114 = { 0, 6, 0, "(//,' HEAT OF FORMATION IS ESSENTIA" "LLY STATIONARY ')", 0 }; static cilist io___115 = { 0, 6, 0, fmt_440, 0 }; static cilist io___116 = { 0, 6, 0, fmt_460, 0 }; static cilist io___117 = { 0, 6, 0, fmt_470, 0 }; static cilist io___118 = { 0, 6, 0, fmt_480, 0 }; static cilist io___119 = { 0, 6, 0, fmt_510, 0 }; static cilist io___120 = { 0, 6, 0, fmt_520, 0 }; static cilist io___121 = { 0, 6, 0, fmt_540, 0 }; static cilist io___125 = { 0, 6, 0, fmt_570, 0 }; static cilist io___126 = { 0, 6, 0, fmt_580, 0 }; static cilist io___127 = { 0, 6, 0, fmt_590, 0 }; /* ************************************************************** */ /* * */ /* THIS SUBROUTINE USES THE BFGS UPDATE TO THE INVERSE HESSIAN * */ /* THE NAME FLEPO IS KEPT IN ORDER TO ALLOW COMPATABILITY WITH * */ /* OTHER PROGRAMS. THE DFP FORMULA IS NOTED IN THE COMMENTS. * */ /* * */ /* ************************************************************** */ /* 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 */ /* * */ /* THIS SUBROUTINE ATTEMPTS TO MINIMIZE A REAL-VALUED FUNCTION OF */ /* THE N-COMPONENT REAL VECTOR XPARAM ACCORDING TO THE */ /* BFGS FORMULA. RELEVANT REFERENCES ARE */ /* BROYDEN, C.G., JOURNAL OF THE INSTITUTE FOR MATHEMATICS AND */ /* APPLICATIONS, VOL. 6 PP 222-231, 1970. */ /* FLETCHER, R., COMPUTER JOURNAL, VOL. 13, PP 317-322, 1970. */ /* GOLDFARB, D. MATHEMATICS OF COMPUTATION, VOL. 24, PP 23-26, 1970. */ /* SHANNO, D.F. MATHEMATICS OF COMPUTATION, VOL. 24, PP 647-656 */ /* 1970. */ /* SEE ALSO SUMMARY IN */ /* SHANNO, D.F., J. OF OPTIMIZATION THEORY AND APPLICATIONS */ /* VOL.46, NO 1 PP 87-94 1985. */ /* THE USER MUST SUPPLY THE SUBROUTINE */ /* COMPFG(NVAR,XPARAM,FUNCT,GRAD,1)WHICH */ /* COMPUTES FUNCTION VALUES FUNCT AND GRADIENTS GRAD AT GIVEN */ /* VALUES FOR THE VARIABLES XPARAM. THE MINIMIZATION PROCEEDS BY ONE */ /* OR MORE BINARY CHOPS WHICH FINDS AN IMPROVED VALUE OF THE */ /* FUNCTION IN THE DIRECTION XPARAM+ALPHA*PVECT, */ /* WHERE XPARAM */ /* IS THE VECTOR OF CURRENT VARIABLE VALUES, ALPHA IS A SCALAR */ /* VARIABLE, AND PVECT IS A SEARCH-DIRECTION VECTOR PROVIDED BY THE */ /* BFGS ALGORITHM. A SEQUENCE OF FUNCT VALUES CONVERGING TO SOME */ /* LOCAL MINIMUM VALUE AND A SEQUENCE OF */ /* XPARAM VECTORS CONVERGING TO THE CORRESPONDING MINIMUM POINT */ /* ARE PRODUCED. */ /* CONVERGENCE TESTS. */ /* HERBERTS TEST: THE ESTIMATED DISTANCE FROM THE CURRENT POINT */ /* POINT TO THE MINIMUM IS LESS THAN TOLERA. */ /* "HERBERTS TEST SATISFIED - GEOMETRY OPTIMIZED" */ /* GRADIENT TEST: THE GRADIENT NORM HAS BECOME LESS THAN TOLERG */ /* TIMES THE SQUARE ROOT OF THE NUMBER OF VARIABLES. */ /* "TEST ON GRADIENT SATISFIED". */ /* XPARAM TEST: THE RELATIVE CHANGE IN XPARAM, MEASURED BY ITS NORM, */ /* OVER ANY TWO SUCCESSIVE ITERATION STEPS DROPS BELOW */ /* TOLERX. */ /* "TEST ON XPARAM SATISFIED". */ /* FUNCTION TEST: THE CALCULATED VALUE OF THE HEAT OF FORMATION */ /* BETWEEN ANY TWO CYCLES IS WITHIN TOLERF OF */ /* EACH OTHER. */ /* "HEAT OF FORMATION TEST SATISFIED" */ /* FOR THE GRADIENT, FUNCTION, AND XPARAM TESTS A FURTHER CONDITION, */ /* THAT NO INDIVIDUAL COMPONENT OF THE GRADIENT IS GREATER */ /* THAN TOLERG, MUST BE SATISFIED, IN WHICH CASE THE */ /* CALCULATION EXITS WITH THE MESSAGE */ /* "PETERS TEST SATISFIED" */ /* WILL BE PRINTED, AND FUNCT AND XPARAM WILL CONTAIN THE LAST */ /* FUNCTION VALUE CUM VARIABLE VALUES REACHED. */ /* SIMILAR UNSUCCESSFUL TERMINATIONS WILL TAKE PLACE IF THE COSINE OF */ /* THE SEARCH DIRECTION TO GRADIENT VECTOR IS LESS THAN RST ON TWO */ /* CONSECUTIVE ITERATIONS. */ /* THE BROYDEN-FLETCHER-GOLDFARB-SHANNO ALGORITHM CHOOSES SEARCH */ /* DIRECTIONS ON THE BASIS OF LOCAL PROPERTIES OF THE FUNCTION. */ /* A MATRIX H, WHICH IN FLEPO IS PRESET WITH THE IDENTITY, IS */ /* MAINTAINED AND UPDATED AT EACH ITERATION STEP. */ /* THE MATRIX DESCRIBES A LOCAL METRIC ON THE SURFACE OF FUNCTION */ /* VALUES ABOVE THE POINT XPARAM. THE SEARCH-DIRECTION VECTOR */ /* PVECT IS SIMPLY A TRANSFORMATION OF THE GRADIENT GRAD */ /* BY THE MATRIX H. THE USER THROWS OUT H AND STARTS AGAIN */ /* WHENEVER THE COSINE OF THE ANGLE BETWEEN GRAD AND PVECT BECOMES */ /* LESS THAN RST. */ /* Parameter adjustments */ --xparam; /* Function Body */ if (icalcn != numcal_1.numcal) { /* THE FOLLOWING CONSTANTS SHOULD BE SET BY THE USER. */ rst = .05; tdel = 6.; sfact = 1.5f; pmste = .1; dell = .01; einc = .3; igg1 = 3; *del = dell; /* THESE CONSTANTS SHOULD BE SET BY THE PROGRAM. */ restrt = i_indx(keywrd_1.keywrd, "RESTART", 80L, 7L) != 0; geook = i_indx(keywrd_1.keywrd, "GEO-OK", 80L, 6L) != 0; time = i_indx(keywrd_1.keywrd, "TIME", 80L, 4L) != 0; /* THE DAVIDON-FLETCHER-POWELL METHOD IS NOT RECOMMENDED */ /* BUT CAN BE INVOKED BY USING THE KEY-WORD 'DFP' */ dfp = i_indx(keywrd_1.keywrd, "DFP", 80L, 3L) != 0; tleft = 3600.; tolerg = 1.; const__ = 1.; minprt = i_indx(keywrd_1.keywrd, "SADDLE", 80L, 6L) + i_indx( keywrd_1.keywrd, "SADDLE", 80L, 6L) == 0; saddle = i_indx(keywrd_1.keywrd, "SADDLE", 80L, 6L) != 0; if (! minprt) { minprt = i_indx(keywrd_1.keywrd, "DEBUG", 80L, 5L) != 0; } 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 L20; } /* L10: */ } /* LIMIT JOB TIME TO MAX. OF ONE YEAR, LARGE JOBTIMES STOP */ /* DUMPS WORKING CORRECTLY AS TLEFT-(CYCLE TIME) = TLEFT! */ tim = min(31557600.,tim); L20: tleft = tim; } tlast = tleft; tdump = 3600.; i = i_indx(keywrd_1.keywrd, " DUMP=", 80L, 6L); 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 L40; } /* L30: */ } L40: ; } tx2 = second_(); tleft = tleft - tx2 + time_1.time0; if (i_indx(keywrd_1.keywrd, "GNORM=", 80L, 6L) != 0) { rootv = 1.; const__ = 1e-20; } else { rootv = sqrt(*nvar + 1e-5); } print = i_indx(keywrd_1.keywrd, "FLEPO", 80L, 5L) != 0; tolerx = const__ * 1e-4; delhof = const__ * .001; tolerf = const__ * .002; tolrg = tolerg; if (i_indx(keywrd_1.keywrd, "FORCE", 80L, 5L) != 0) { tolerx = 1e-5; tolerf = 2e-4; tolerg = .1; delhof = 1e-4; } if (i_indx(keywrd_1.keywrd, "PREC", 80L, 4L) != 0) { tolerx *= .01; delhof *= .01; tolerf *= .01; tolerg *= .1; einc *= .01f; } if (i_indx(keywrd_1.keywrd, "GNORM=", 80L, 6L) != 0) { i__1 = i_indx(keywrd_1.keywrd, "GNORM=", 80L, 6L); tolerg = reada_(keywrd_1.keywrd, &i__1, 80L); if (! geook && i_indx(keywrd_1.keywrd, "LET", 80L, 3L) == 0 && tolerg < 1e-4) { s_wsfe(&io___49); do_fio(&c__1, " GNORM HAS BEEN SET TOO LOW, RESET TO 0.0001", 45L); e_wsfe(); tolerg = 1e-4; tolrg = tolerg; } } else { tolerg /= rootv; } } /* THE FOLLOWING CONSTANTS SHOULD BE SET TO SOME ARBITARY LARGE VALUE. */ *drop = 1e15; *frepf = 1e15; /* AND FINALLY, THE FOLLOWING CONSTANTS ARE CALCULATED. */ ihdim = *nvar * (*nvar + 1) / 2; cncadd = 1. / rootv; if (cncadd > .15) { cncadd = .15; } /* FIRST, WE INITIALIZE THE VARIABLES. */ absmin = 1e6; itry1 = 0; itry2 = 0; okc = TRUE_; okf = TRUE_; *jcyc = 0; *lnstop = 1; irepet = 1; *alpha = 1.; *pnorm = 1.; *jnrst = 0; *cycmx = 0.; *cos__ = 0.; *totime = 0.; *ncount = 1; resfil = FALSE_; if (saddle) { /* WE DON'T NEED HIGH PRECISION DURING A SADDLE-POINT CALCULATION. */ if (*nvar > 0) { gradnt_1.gnorm = sqrt(dot_(gradnt_1.grad, gradnt_1.grad, nvar)) - 3.; } if (gradnt_1.gnorm > 10.) { gradnt_1.gnorm = 10.; } if (gradnt_1.gnorm > 1.) { tolerg = tolrg * gradnt_1.gnorm; } s_wsfe(&io___59); do_fio(&c__1, (char *)&tolerg, (ftnlen)sizeof(doublereal)); e_wsfe(); } if (restrt && icalcn != numcal_1.numcal) { mdfp[8] = 0; dfpsav_(totime, &xparam[1], gd, xlast, funct1, mdfp, xdfp); i = (integer) (*totime / 1e6); *totime -= i * 1000000; s_wsfe(&io___62); do_fio(&c__1, (char *)&(*totime), (ftnlen)sizeof(doublereal)); e_wsfe(); numscf_1.nscf = mdfp[4]; if (i_indx(keywrd_1.keywrd, "1SCF", 80L, 4L) != 0) { compfg_(&xparam[1], &c_true, funct1, &c_true, gradnt_1.grad, & c_true); icalcn = numcal_1.numcal; mesage_1.iflepo = 1; time_1.time0 -= *totime; *totime = 0.; return 0; } } else { *totime = 0.; /* CALCULATE THE VALUE OF THE FUNCTION -> FUNCT1, AND GRADIENTS -> GRA D. */ /* NORMAL SET-UP OF FUNCT1 AND GRAD, DONE ONCE ONLY. */ compfg_(&xparam[1], &c_true, funct1, &c_true, gradnt_1.grad, &c_true); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { /* L50: */ gd[i - 1] = gradnt_1.grad[i - 1]; } } icalcn = numcal_1.numcal; if (*nvar != 0) { gradnt_1.gnorm = sqrt(dot_(gradnt_1.grad, gradnt_1.grad, nvar)); } mesage_1.iflepo = 1; if (i_indx(keywrd_1.keywrd, "1SCF", 80L, 4L) != 0) { return 0; } mesage_1.iflepo = 2; if (gradnt_1.gnorm < tolerg || *nvar == 0) { last_1.last = 1; if (restrt) { compfg_(&xparam[1], &c_true, funct1, &c_true, gradnt_1.grad, & c_true); } else { compfg_(&xparam[1], &c_true, funct1, &c_true, gradnt_1.grad, & c_false); } return 0; } tx1 = second_(); tleft = tleft - tx1 + tx2; /* * */ /* START OF EACH ITERATION CYCLE ... */ /* * */ reset = FALSE_; goto L80; L60: if (*cos__ < rst) { i__1 = *nvar; for (i = 1; i <= i__1; ++i) { /* L70: */ gd[i - 1] = .5; } } L80: ++(*jcyc); ++(*jnrst); if (*lnstop != 1 && *cos__ > rst) { goto L160; } /* * */ /* RESTART SECTION */ /* * */ L90: reset = TRUE_; i__1 = *nvar; for (i = 1; i <= i__1; ++i) { xd[i - 1] = xparam[i] - d_sign(del, &gradnt_1.grad[i - 1]); /* L100: */ } /* THIS CALL OF COMPFG IS USED TO CALCULATE THE SECOND-ORDER MATRIX IN H */ /* IF THE NEW POINT HAPPENS TO IMPROVE THE RESULT, THEN IT IS KEPT. */ /* OTHERWISE IT IS SCRAPPED, BUT STILL THE SECOND-ORDER MATRIX IS O.K. */ compfg_(xd, &c_true, &funct2, &c_true, gd, &c_true); if (! geook && sqrt(dot_(gd, gd, nvar)) / gradnt_1.gnorm > 10. && gradnt_1.gnorm / rootv > 20. && *jcyc > 2) { /* THE GEOMETRY IS BADLY SPECIFIED IN THAT MINOR CHANGES IN INTERNAL */ /* COORDINATES LEAD TO LARGE CHANGES IN CARTESIAN COORDINATES, AND TH ESE */ /* LARGE CHANGES ARE BETWEEN PAIRS OF ATOMS THAT ARE CHEMICALLY BONDE D */ /* TOGETHER. */ s_wsfe(&io___67); do_fio(&c__1, (char *)&gradnt_1.gnorm, (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___68); 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(); gdnorm = sqrt(dot_(gd, gd, nvar)); s_wsfe(&io___70); do_fio(&c__1, (char *)&gdnorm, (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___71); i__1 = *nvar; for (i = 1; i <= i__1; ++i) { do_fio(&c__1, (char *)&gd[i - 1], (ftnlen)sizeof(doublereal)); } e_wsfe(); s_wsfe(&io___72); e_wsfe(); s_wsfe(&io___73); e_wsfe(); s_wsfe(&io___74); e_wsfe(); geout_(); s_stop("", 0L); } ++(*ncount); i__1 = ihdim; for (i = 1; i <= i__1; ++i) { /* L110: */ fmatrx_1.hesinv[i - 1] = 0.; } i__1 = *nvar; for (i = 1; i <= i__1; ++i) { ii = i * (i + 1) / 2; deltag = gradnt_1.grad[i - 1] - gd[i - 1]; if (abs(deltag) < 1e-10) { deltag = 1e-10; } ggd = (d__1 = gradnt_1.grad[i - 1], abs(d__1)); if (funct2 < *funct1) { ggd = (d__1 = gd[i - 1], abs(d__1)); } fmatrx_1.hesinv[ii - 1] = d_sign(del, &gradnt_1.grad[i - 1]) / deltag; if (fmatrx_1.hesinv[ii - 1] < 0. && ggd < 1e-12) { fmatrx_1.hesinv[ii - 1] = .01; } if (fmatrx_1.hesinv[ii - 1] < 0.) { fmatrx_1.hesinv[ii - 1] = *del * 6. / ggd; } /* Computing MIN */ d__2 = fmatrx_1.hesinv[ii - 1], d__3 = (d__1 = pmste / max(1e-12,ggd), abs(d__1)); fmatrx_1.hesinv[ii - 1] = min(d__2,d__3); /* L120: */ } *jnrst = 0; if (*jcyc < 2) { gravec_1.cosine = 1.; } if (funct2 >= *funct1) { if (print) { s_wsfe(&io___78); do_fio(&c__1, (char *)&(*funct1), (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&funct2, (ftnlen)sizeof(doublereal)); e_wsfe(); } compfg_(&xparam[1], &c_true, funct1, &c_true, gd, &c_false); gravec_1.cosine = 1.; } else { if (print) { s_wsfe(&io___79); do_fio(&c__1, (char *)&(*funct1), (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&funct2, (ftnlen)sizeof(doublereal)); e_wsfe(); } *funct1 = funct2; gradnt_1.gnorm = 0.; i__1 = *nvar; for (i = 1; i <= i__1; ++i) { xparam[i] = xd[i - 1]; gradnt_1.grad[i - 1] = gd[i - 1]; /* L150: */ /* Computing 2nd power */ d__1 = gradnt_1.grad[i - 1]; gradnt_1.gnorm += d__1 * d__1; } gradnt_1.gnorm = sqrt(gradnt_1.gnorm); } goto L210; /* * */ /* UPDATE VARIABLE-METRIC MATRIX */ /* * */ L160: pty = 0.; ++(*jnrst); yhy = 0.; i__1 = *nvar; for (i = 1; i <= i__1; ++i) { s = 0.; i__2 = *nvar; for (j = 1; j <= i__2; ++j) { if (j > i) { k = j * (j - 1) / 2 + i; } else { k = i * (i - 1) / 2 + j; } /* L170: */ s += fmatrx_1.hesinv[k - 1] * (gradnt_1.grad[j - 1] - glast[j - 1] ); } gg[i - 1] = s; y = gradnt_1.grad[i - 1] - glast[i - 1]; yhy += gg[i - 1] * y; /* L180: */ pty += (xparam[i] - xlast[i - 1]) * y; } if (dfp) { i__1 = *nvar; for (i = 1; i <= i__1; ++i) { pt = xparam[i] - xlast[i - 1]; i__2 = *nvar; for (j = i; j <= i__2; ++j) { k = j * (j - 1) / 2 + i; p = xparam[j] - xlast[j - 1]; /* START OF DAVIDON FLETCHER POWELL FORMULA */ fmatrx_1.hesinv[k - 1] = fmatrx_1.hesinv[k - 1] + pt * p / pty - gg[i - 1] * gg[j - 1] / yhy; /* END OF DAVIDON FLETCHER POWELL FORMULA */ /* L190: */ } } } else { /* BFGS FORMULA ON NEXT LINE */ yhy = yhy / pty + 1.; /* BFGS FORMULA ON LAST LINE */ i__2 = *nvar; for (i = 1; i <= i__2; ++i) { pt = xparam[i] - xlast[i - 1]; i__1 = *nvar; for (j = i; j <= i__1; ++j) { p = xparam[j] - xlast[j - 1]; k = j * (j - 1) / 2 + i; /* START OF BFGS FORMULA ON NEXT LINE */ fmatrx_1.hesinv[k - 1] = fmatrx_1.hesinv[k - 1] - (gg[j - 1] * pt + p * gg[i - 1]) / pty + yhy * p * pt / pty; /* END OF BFGS FORMULA */ /* L200: */ } } } /* * */ /* ESTABLISH NEW SEARCH DIRECTION */ /* * */ L210: pnlast = *pnorm; *pnorm = 0.; dott = 0.; i__1 = *nvar; for (k = 1; k <= i__1; ++k) { s = 0.; i__2 = *nvar; for (i = 1; i <= i__2; ++i) { ij = max(i,k); s -= fmatrx_1.hesinv[ij * (ij - 1) / 2 + i + k - ij - 1] * gradnt_1.grad[i - 1]; /* L220: */ } pvect[k - 1] = s; /* Computing 2nd power */ d__1 = pvect[k - 1]; *pnorm += d__1 * d__1; /* L230: */ dott += pvect[k - 1] * gradnt_1.grad[k - 1]; } *pnorm = sqrt(*pnorm); *cos__ = -dott / (*pnorm * gradnt_1.gnorm); if (*jnrst == 0) { goto L260; } if (*cos__ <= cncadd && *drop > 1.) { goto L240; } if (*cos__ <= rst) { goto L240; } goto L260; L240: *pnorm = pnlast; if (print) { s_wsfe(&io___93); do_fio(&c__1, (char *)&(*cos__), (ftnlen)sizeof(doublereal)); e_wsfe(); } goto L90; L260: if (print) { s_wsfe(&io___94); e_wsfe(); nto6 = (*nvar - 1) / 6 + 1; iinc1 = -5; i__1 = nto6; for (i = 1; i <= i__1; ++i) { s_wsfe(&io___97); e_wsfe(); iinc1 += 6; /* Computing MIN */ i__2 = iinc1 + 5; iinc2 = min(i__2,*nvar); s_wsfe(&io___99); i__2 = iinc2; for (j = iinc1; j <= i__2; ++j) { do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer)); } e_wsfe(); s_wsfe(&io___100); i__2 = iinc2; for (j = iinc1; j <= i__2; ++j) { do_fio(&c__1, (char *)&xparam[j], (ftnlen)sizeof(doublereal)); } e_wsfe(); s_wsfe(&io___101); i__2 = iinc2; for (j = iinc1; j <= i__2; ++j) { do_fio(&c__1, (char *)&gradnt_1.grad[j - 1], (ftnlen)sizeof( doublereal)); } e_wsfe(); s_wsfe(&io___102); i__2 = iinc2; for (j = iinc1; j <= i__2; ++j) { d__1 = *alpha * pvect[j - 1]; do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof(doublereal)); } e_wsfe(); /* L320: */ } } *lnstop = 0; *alpha = *alpha * pnlast / *pnorm; i__1 = *nvar; for (i = 1; i <= i__1; ++i) { glast[i - 1] = gradnt_1.grad[i - 1]; /* L330: */ xlast[i - 1] = xparam[i]; } if (*jnrst == 0) { *alpha = 1.; } *drop = (d__1 = *alpha * dott, abs(d__1)); if (*jnrst != 0 && *drop < delhof) { if (minprt) { s_wsfe(&io___103); e_wsfe(); } /* FLEPO IS ENDING PROPERLY. THIS IS IMMEDIATELY BEFORE THE RETURN. */ last_1.last = 1; compfg_(&xparam[1], &c_true, &funct, &c_true, gradnt_1.grad, &c_false) ; mesage_1.iflepo = 3; time_1.time0 -= *totime; *totime = 0.; return 0; } smval = *funct1; if (gradnt_1.gnorm > 1. && *pnorm > .001) { linmin_(&xparam[1], alpha, pvect, nvar, funct1, &okf, &okc); } else { /* DO A BINARY CHOP TO LOCATE THE MINIMUM */ *alpha = 1.; /* SOMETIMES PNORM IS TOO LARGE, THEREFORE TRIM ALPHA BACK */ if (*pnorm > .1) { *alpha = .1 / *pnorm; } loop = 0; L350: i__1 = *nvar; for (i = 1; i <= i__1; ++i) { /* L360: */ xparam[i] = xlast[i - 1] + pvect[i - 1] * *alpha; } compfg_(&xparam[1], &c_true, funct1, &c_true, gradnt_1.grad, &c_false) ; /* 1.D-5 IS TO PREVENT LOOPING */ if (*funct1 > smval + 1e-5 && *pnorm * *alpha > 1e-6 && loop < 10) { ++loop; *alpha *= .5; goto L350; } } ++(*ncount); if (! okf) { *lnstop = 1; if (minprt) { s_wsfe(&io___107); e_wsfe(); } *funct1 = smval; i__1 = *nvar; for (i = 1; i <= i__1; ++i) { gradnt_1.grad[i - 1] = glast[i - 1]; xparam[i] = xlast[i - 1]; /* L380: */ } if (*jnrst == 0) { s_wsfe(&io___108); e_wsfe(); /* FLEPO IS ENDING BADLY. THIS IS IMMEDIATELY BEFORE THE RETURN. */ last_1.last = 1; compfg_(&xparam[1], &c_true, &funct, &c_true, gradnt_1.grad, & c_true); mesage_1.iflepo = 4; time_1.time0 -= *totime; *totime = 0.; return 0; } if (print) { s_wsfe(&io___109); e_wsfe(); } *cos__ = 0.; goto L560; } /* WE WANT ACCURATE DERIVATIVES AT THIS POINT */ /* LINMIN DOES NOT GENERATE ANY DERIVATIVES, THEREFORE COMPFG MUST BE */ /* CALLED TO END THE SEARCH */ if (reset) { i__1 = *nvar; for (j = 1; j <= i__1; ++j) { /* L410: */ gradnt_1.grad[j - 1] = 0.; } reset = FALSE_; } compfg_(&xparam[1], &c_true, funct1, &c_true, gradnt_1.grad, &c_true); gradnt_1.gnorm = sqrt(dot_(gradnt_1.grad, gradnt_1.grad, nvar)); if (! okc && minprt) { s_wsfe(&io___110); do_fio(&c__1, (char *)&(*jcyc), (ftnlen)sizeof(integer)); e_wsfe(); } xn = 0.; i__1 = *nvar; for (k = 1; k <= i__1; ++k) { /* L430: */ /* Computing 2nd power */ d__1 = xparam[k]; xn += d__1 * d__1; } xn = sqrt(xn); tx = (d__1 = *alpha * *pnorm, abs(d__1)); if (xn != 0.) { tx /= xn; } tf = smval - *funct1; if (*alpha < .01 || absmin - smval < 1e-7) { ++itry2; if (*cos__ < rst) { itry2 = 0; } if (itry2 == 5) { /* RE-CALCULATE ALL DERIVATIVES IF HALF-ELECTRON USED. */ *cos__ = -1.; } } else { if (*alpha > .1) { itry2 = 0; } } if (absmin - smval < 1e-7) { ++itry1; if (itry1 > 10) { s_wsfe(&io___114); e_wsfe(); goto L550; } } else { itry1 = 0; absmin = smval; } if (print) { s_wsfe(&io___115); do_fio(&c__1, (char *)&(*ncount), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*cos__), (ftnlen)sizeof(doublereal)); d__1 = tx * xn; do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&(*alpha), (ftnlen)sizeof(doublereal)); d__2 = -(*drop); do_fio(&c__1, (char *)&d__2, (ftnlen)sizeof(doublereal)); d__3 = -tf; do_fio(&c__1, (char *)&d__3, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&gradnt_1.gnorm, (ftnlen)sizeof(doublereal)); e_wsfe(); } /* L450: */ if (tx <= tolerx) { if (minprt) { s_wsfe(&io___116); e_wsfe(); } goto L490; } if (abs(tf) <= tolerf) { if (minprt) { s_wsfe(&io___117); e_wsfe(); } goto L490; } if (gradnt_1.gnorm <= tolerg * rootv) { if (minprt) { s_wsfe(&io___118); e_wsfe(); } goto L490; } goto L560; L490: i__1 = *nvar; for (i = 1; i <= i__1; ++i) { if ((d__1 = gradnt_1.grad[i - 1], abs(d__1)) > tolerg) { ++irepet; if (irepet > 1) { goto L500; } *frepf = *funct1; *cos__ = 0.; L500: if (minprt) { s_wsfe(&io___119); do_fio(&c__1, (char *)&tolerg, (ftnlen)sizeof(doublereal)); e_wsfe(); } if ((d__1 = *funct1 - *frepf, abs(d__1)) > einc) { irepet = 0; } if (irepet > igg1) { s_wsfe(&io___120); do_fio(&c__1, (char *)&igg1, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&einc, (ftnlen)sizeof(doublereal)); e_wsfe(); last_1.last = 1; compfg_(&xparam[1], &c_true, &funct, &c_true, gradnt_1.grad, & c_false); mesage_1.iflepo = 8; time_1.time0 -= *totime; *totime = 0.; return 0; } else { goto L560; } } /* L530: */ } if (minprt) { s_wsfe(&io___121); e_wsfe(); } L550: last_1.last = 1; compfg_(&xparam[1], &c_true, &funct, &c_true, gradnt_1.grad, &c_false); mesage_1.iflepo = 6; time_1.time0 -= *totime; *totime = 0.; return 0; /* ALL TESTS HAVE FAILED, WE NEED TO DO ANOTHER CYCLE. */ L560: bsmvf = (d__1 = smval - *funct1, abs(d__1)); if (bsmvf > 10.) { *cos__ = 0.; } *del = .002; if (bsmvf > 1.) { *del = dell / 2.; } if (bsmvf > 5.) { *del = dell; } tx2 = second_(); tcycle = tx2 - tx1; tx1 = tx2; /* END OF ITERATION LOOP, EVERYTHING IS STILL O.K. SO GO TO */ /* NEXT ITERATION, IF THERE IS ENOUGH TIME LEFT. */ if (tcycle < 1e5) { /* Computing MAX */ d__1 = *cycmx * .8; *cycmx = max(d__1,tcycle); } tleft -= tcycle; if (tleft < 0.) { tleft = -.1; } if (tcycle > 1e5) { tcycle = 0.; } if (tlast - tleft > tdump) { totim = *totime + second_() - time_1.time0; tlast = tleft; mdfp[8] = 2; resfil = TRUE_; mdfp[4] = numscf_1.nscf; dfpsav_(&totim, &xparam[1], gd, xlast, funct1, mdfp, xdfp); } if (resfil) { if (minprt) { s_wsfe(&io___125); d__1 = min(tleft,9999999.9); do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof(doublereal)); d__2 = min(gradnt_1.gnorm,999999.999); do_fio(&c__1, (char *)&d__2, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&(*funct1), (ftnlen)sizeof(doublereal)); e_wsfe(); } resfil = FALSE_; } else { if (minprt) { s_wsfe(&io___126); do_fio(&c__1, (char *)&(*jcyc), (ftnlen)sizeof(integer)); d__1 = min(tcycle,9999.99); do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof(doublereal)); d__2 = min(tleft,9999999.9); do_fio(&c__1, (char *)&d__2, (ftnlen)sizeof(doublereal)); d__3 = min(gradnt_1.gnorm,999999.999); do_fio(&c__1, (char *)&d__3, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&(*funct1), (ftnlen)sizeof(doublereal)); e_wsfe(); } } if (tleft > sfact * *cycmx) { goto L60; } s_wsfe(&io___127); e_wsfe(); *totime = *totime + second_() - time_1.time0; mdfp[8] = 1; mdfp[4] = numscf_1.nscf; dfpsav_(totime, &xparam[1], gd, xlast, funct1, mdfp, xdfp); s_stop("", 0L); return 0; } /* flepo_ */ #undef cos__ #undef del #undef lnstop #undef ncount #undef totime #undef jnrst #undef pnorm #undef cycmx #undef frepf #undef alpha #undef drop #undef xdfp #undef jcyc #undef mdfp