/* iter.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 f[23220], fb[23220]; } fokmat_; #define fokmat_1 fokmat_ struct { doublereal p[23220], pa[23220], pb[23220]; } densty_; #define densty_1 densty_ struct { doublereal c[46225], eigs[215], cbeta[46225], eigb[215]; } vector_; #define vector_1 vector_ struct { doublereal dumy[258], gnorm; } gradnt_; #define gradnt_1 gradnt_ struct { integer last; } last_; #define last_1 last_ struct { integer iflepo, iiter; } mesage_; #define mesage_1 mesage_ struct { doublereal atheat; } atheat_; #define atheat_1 atheat_ struct { doublereal enuclr; } enuclr_; #define enuclr_1 enuclr_ struct { doublereal xi, xj, xk; } citerm_; #define citerm_1 citerm_ struct { integer latom, lparam; doublereal react[200]; } path_; #define path_1 path_ struct { integer numcal; } numcal_; #define numcal_1 numcal_ struct { doublereal time0; } time_; #define time_1 time_ 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 dummy[215], pdiag[215]; } molorb_; #define molorb_1 molorb_ struct { char keywrd[80]; } keywrd_; #define keywrd_1 keywrd_ struct { integer nscf; } numscf_; #define numscf_1 numscf_ /* Table of constant values */ static integer c__3 = 3; static integer c__1 = 1; static logical c_false = FALSE_; static doublereal c_b87 = 9999.; static integer c__23220 = 23220; static integer c__0 = 0; /* Subroutine */ int iter_(doublereal *h, doublereal *w, real *wj, real *wk, doublereal *ee, logical *fulscf, logical *rand) { /* Initialized data */ static integer icalcn = 0; static logical debug = FALSE_; static logical prtfok = FALSE_; static logical prteig = FALSE_; static logical prtden = FALSE_; static logical prteng = FALSE_; static logical prt1el = FALSE_; static char abprt[5*3] = " " "ALPHA" " BETA"; /* Format strings */ static char fmt_200[] = "(\002 FOCK MATRIX ON ITERATION\002,i3)"; static char fmt_220[] = "(10x,a,\002 EIGENVALUES ON ITERATION\002,i3,/1" "0(6g13.6,/))"; static char fmt_270[] = "(//10x,\002\"\"\"\"\"\"\"\"\"\"\"\"\"UNABLE TO " "ACHIEVE SELF-CONSISTENCE\002,/)"; static char fmt_280[] = "(//,10x,\002DELTAE= \002,e12.4,5x,\002DELTAP=" " \002,e12.4,///)"; /* System generated locals */ integer i__1, i__2; doublereal d__1, d__2, d__3; olist o__1; alist al__1; /* Builtin functions */ integer i_indx(char *, char *, ftnlen, ftnlen), f_open(olist *), f_rew( alist *), s_rsue(cilist *), do_uio(integer *, char *, ftnlen), e_rsue(void), s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void); double d_sign(doublereal *, doublereal *); /* Subroutine */ int s_stop(char *, ftnlen); /* Local variables */ extern /* Subroutine */ int diag_(doublereal *, doublereal *, integer *, doublereal *, integer *, integer *); static doublereal diff; extern doublereal meci_(doublereal *, doublereal *, doublereal *, doublereal *, integer *, integer *, integer *, logical *); static doublereal escf, eold; static integer ialp, jalp, jbet, ibet, ncis; extern /* Subroutine */ int cnvg_(doublereal *, doublereal *, doublereal * , integer *, integer *, doublereal *); static doublereal pold[23220]; extern /* Subroutine */ int swap_(doublereal *, integer *, integer *, integer *, integer *); static integer nmos; static logical frst; static integer na1el, na2el, nb2el, nb1el; extern /* Subroutine */ int fock2_(doublereal *, doublereal *, doublereal *, doublereal *, real *, real *, integer *, integer *, integer *, integer *), fock1_(doublereal *, doublereal *, doublereal *, doublereal *); static doublereal time1, pold2[23220], pold3[615]; static integer i, j; extern doublereal reada_(char *, integer *, ftnlen); static integer l; static logical halfe, makea, makeb; static integer modea, modeb; static doublereal pbold[23220]; static logical force; static integer ifill; static logical newdg, ready, capps; static integer ihomo; static doublereal shift; static logical bfrst, times; static integer iredy; static doublereal trans; static integer niter; static doublereal scorr; extern /* Subroutine */ int pulay_(doublereal *, doublereal *, integer *, doublereal *, doublereal *, doublereal *, integer *, integer *, integer *, logical *, doublereal *), write_(doublereal *, doublereal *); static logical prtpl; static doublereal w1, w2, pbold2[23220], pbold3[615], titer1, titer2; static logical ci; static doublereal pl; static logical camkin, allcon, excitd; static integer ihomob; extern /* Subroutine */ int epseta_(doublereal *, doublereal *); static logical oknewd, incitr, prtvec; static doublereal ar1[46440], ar2[46440], ar3[46440], ar4[46440], br1[ 46440], br2[46440], br3[46440], br4[46440]; static logical minprt; static doublereal shfmax; static logical okpuly; static integer linear; static doublereal bshift, scfcrt, pltest; static integer itrmax; extern /* Subroutine */ int getarg_(integer *, char *, ftnlen); static doublereal random, selcon; extern doublereal second_(void); extern /* Subroutine */ int vecprt_(doublereal *, integer *); static doublereal shiftb, tenold; extern doublereal helect_(integer *, doublereal *, doublereal *, doublereal *), capcor_(integer *, integer *, integer *, integer *, doublereal *, doublereal *); static doublereal sellim; extern /* Subroutine */ int interp_(integer *, integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *), matout_( doublereal *, doublereal *, integer *, integer *, integer *), densit_(doublereal *, integer *, integer *, integer *, integer *, doublereal *, doublereal *, integer *); static doublereal eta, plb; static logical uhf; static doublereal ten, eps, sum; extern /* Subroutine */ int rsp_(doublereal *, integer *, integer *, doublereal *, doublereal *); static logical scf1; /* Fortran I/O blocks */ static cilist io___48 = { 0, 10, 0, 0, 0 }; static cilist io___49 = { 0, 10, 0, 0, 0 }; static cilist io___57 = { 0, 6, 0, "(' SCF CRITERION =',G14.4)", 0 }; static cilist io___58 = { 0, 6, 0, "(//2X,' THERE IS A RISK OF INFINITE " "LOOPING WITH', ' THE SCFCRT LESS THAN 1.D-12')", 0 }; static cilist io___59 = { 0, 6, 0, "(' SCF CRITERION =',G14.4)", 0 }; static cilist io___64 = { 0, 6, 0, "(' SELCON, GNORM',2G16.7)", 0 }; static cilist io___66 = { 0, 6, 0, "(//10X,'ONE-ELECTRON MATRIX AT ENTRA" "NCE TO ITER')", 0 }; static cilist io___75 = { 0, 6, 0, "(//,' ALL CONVERGERS ARE NOW FORCED " "ON',/ ' SHIFT=10, PULAY ON, CAMP-KING ON',/ " " ' AND ITERATION COUNTER RESET',//)", 0 } ; static cilist io___80 = { 0, 6, 0, fmt_200, 0 }; static cilist io___84 = { 0, 6, 0, "(' ITERATION',I3,' PLS=',2E10.3,' EN" "ERGY ', F14.7,' DELTAE',F13.7)", 0 }; static cilist io___94 = { 0, 6, 0, "(//10X,A, " " ' EIGENVECTORS AND EIGENVALUES ON ITERATION',I" "3)", 0 }; static cilist io___95 = { 0, 6, 0, fmt_220, 0 }; static cilist io___104 = { 0, 6, 0, "(//10X,A,' EIGENVECTORS AND EIGENVA" "LUES ON ', 'ITERATION',I3)", 0 }; static cilist io___105 = { 0, 6, 0, fmt_220, 0 }; static cilist io___106 = { 0, 6, 0, "(' DENSITY MATRIX ON ITERATION',I4)", 0 }; static cilist io___107 = { 0, 6, 0, "(' \"\"\"\"\"\"\"\"\"\"\"\"\"\"\"UN" "ABLE TO ACHIEVE SELF-CONSISTENCE, JOB CONTINUING')", 0 }; static cilist io___108 = { 0, 6, 0, fmt_270, 0 }; static cilist io___109 = { 0, 6, 0, fmt_280, 0 }; static cilist io___110 = { 0, 6, 0, "(27X,'AFTER MECI, ENERGY ',F14.7)", 0 }; static cilist io___112 = { 0, 6, 0, "(' TIME FOR SCF CALCULATION',F8.2, " " ' INTEGRAL',F8.2)", 0 }; static cilist io___113 = { 0, 6, 0, "(' NO. OF ITERATIONS =',I6)", 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 */ /* *********************************************************************** */ /* ITER GENERATES A SCF FIELD AND RETURNS THE ENERGY IN "ENERGY" */ /* THE MAIN ARRAYS USED IN ITER ARE: */ /* P ONLY EVER CONTAINS THE TOTAL DENSITY MATRIX */ /* PA ONLY EVER CONTAINS THE ALPHA DENSITY MATRIX */ /* PB ONLY EVER CONTAINS THE BETA DENSITY MATRIX */ /* C ONLY EVER CONTAINS THE EIGENVECTORS */ /* H ONLY EVER CONTAINS THE ONE-ELECTRON MATRIX */ /* F STARTS OFF CONTAINING THE ONE-ELECTRON MATRIX, */ /* AND IS USED TO HOLD THE FOCK MATRIX */ /* W ONLY EVER CONTAINS THE TWO-ELECTRON MATRIX */ /* THE MAIN INTEGERS CONSTANTS IN ITER ARE: */ /* LINEAR SIZE OF PACKED TRIANGLE = NORBS*(NORBS+1)/2 */ /* THE MAIN INTEGER VARIABLES ARE */ /* NITER NUMBER OF ITERATIONS EXECUTED */ /* PRINCIPAL REFERENCES: */ /* ON MNDO: "GROUND STATES OF MOLECULES. 38. THE MNDO METHOD. */ /* APPROXIMATIONS AND PARAMETERS." */ /* DEWAR, M.J.S., THIEL,W., J. AM. CHEM. SOC.,99,4899,(1977). */ /* ON SHIFT: "THE DYNAMIC 'LEVEL SHIFT' METHOD FOR IMPROVING THE */ /* CONVERGENCE OF THE SCF PROCEDURE", A. V. MITIN, J. COMP. */ /* CHEM. 9, 107-110 (1988) */ /* ON HALF-ELECTRON: "MINDO/3 COMPARISON OF THE GENERALIZED S.C.F. */ /* COUPLING OPERATOR AND "HALF-ELECTRON" METHODS FOR */ /* CALCULATING THE ENERGIES AND GEOMETRIES OF OPEN SHELL */ /* SYSTEMS" */ /* DEWAR, M.J.S., OLIVELLA, S., J. CHEM. SOC. FARA. II, */ /* 75,829,(1979). */ /* ON PULAY'S CONVERGER: "CONVERGANCE ACCELERATION OF ITERATIVE */ /* SEQUENCES. THE CASE OF SCF ITERATION", PULAY, P., */ /* CHEM. PHYS. LETT, 73, 393, (1980). */ /* ON CNVG: IT ENCORPORATES THE IMPROVED ITERATION SCHEME (IIS) BY */ /* PIOTR BADZIAG & FRITZ SOLMS. ACCEPTED FOR PUBLISHING */ /* IN COMPUTERS & CHEMISTRY */ /* ON PSEUDODIAGONALISATION: "FAST SEMIEMPIRICAL CALCULATIONS", */ /* STEWART. J.J.P., CSASZAR, P., PULAY, P., J. COMP. CHEM., */ /* 3, 227, (1982) */ /* *********************************************************************** */ /* *********************************************************************** */ /* * */ /* IF, FOR ANY REASON, ITER FAILS TO GENERATE AN SCF, THEN UNCOMMENT * */ /* ALL LINES THAT BEGIN WITH A ''. THE EFFECT OF DOING THIS IS TO * */ /* INCREASE THE STORAGE REQUIREMENT BY ABOUT 20% AND TO ALLOW PULAY'S * */ /* AND CAMP AND KING'S CONVERGERS TO BE USED. * */ /* * */ /* *********************************************************************** */ /* Parameter adjustments */ --wk; --wj; --w; --h; /* Function Body */ /* INITIALIZE */ ifill = 0; /* Computing MAX */ i__1 = 1, i__2 = molkst_1.nclose + molkst_1.nalpha; ihomo = max(i__1,i__2); /* Computing MAX */ i__1 = 1, i__2 = molkst_1.nclose + molkst_1.nbeta; ihomob = max(i__1,i__2); eold = 100.; ready = FALSE_; if (icalcn != numcal_1.numcal) { epseta_(&eps, &eta); /* ULTIMATE SCF CRITERION: HEAT OF FORMATION CONVERGED WITHIN A FACTO R */ /* OF 10 OF THE LIMITING PRECISION OF THE COMPUTER */ eps = eps * 23.061 * 2.; shift = 0.; icalcn = numcal_1.numcal; shfmax = 20.; linear = molkst_1.norbs * (molkst_1.norbs + 1) / 2; /* DEBUG KEY-WORDS WORKED OUT */ debug = i_indx(keywrd_1.keywrd, "DEBUG", 80L, 5L) != 0; minprt = i_indx(keywrd_1.keywrd, "SADDLE", 80L, 6L) + path_1.latom == 0 || debug; prteig = i_indx(keywrd_1.keywrd, "EIGS", 80L, 4L) != 0; prteng = i_indx(keywrd_1.keywrd, "ENERGY", 80L, 6L) != 0; prtpl = i_indx(keywrd_1.keywrd, " PL ", 80L, 4L) != 0; prt1el = i_indx(keywrd_1.keywrd, "1ELEC", 80L, 5L) != 0 && debug; prtden = i_indx(keywrd_1.keywrd, "DENS", 80L, 4L) != 0 && debug; prtfok = i_indx(keywrd_1.keywrd, "FOCK", 80L, 4L) != 0 && debug; prtvec = i_indx(keywrd_1.keywrd, "VECT", 80L, 4L) != 0 && debug; debug = i_indx(keywrd_1.keywrd, "ITER", 80L, 4L) != 0; /* INITIALIZE SOME LOGICALS AND CONSTANTS */ newdg = FALSE_; pl = 1.; bshift = -15.; shift = 1.; /* SCFCRT AND PLTEST ARE MACHINE-PRECISION DEPENDENT */ scfcrt = 1e-6; pltest = 1e-4; itrmax = 200; nmos = 0; ncis = 0; na2el = molkst_1.nclose; na1el = molkst_1.nalpha + molkst_1.nopen; nb2el = 0; nb1el = molkst_1.nbeta + molkst_1.nopen; /* USE KEY-WORDS TO ASSIGN VARIOUS CONSTANTS */ if (i_indx(keywrd_1.keywrd, "C.I.", 80L, 4L) != 0) { i__1 = i_indx(keywrd_1.keywrd, "C.I.", 80L, 4L) + 5; nmos = (integer) reada_(keywrd_1.keywrd, &i__1, 80L); } if (i_indx(keywrd_1.keywrd, "MICROS", 80L, 6L) != 0) { i__1 = i_indx(keywrd_1.keywrd, "MICROS", 80L, 6L); ncis = (integer) reada_(keywrd_1.keywrd, &i__1, 80L); } if (i_indx(keywrd_1.keywrd, "FILL", 80L, 4L) != 0) { i__1 = i_indx(keywrd_1.keywrd, "FILL", 80L, 4L); ifill = (integer) (-reada_(keywrd_1.keywrd, &i__1, 80L)); } if (i_indx(keywrd_1.keywrd, "SHIFT", 80L, 5L) != 0) { i__1 = i_indx(keywrd_1.keywrd, "SHIFT", 80L, 5L); bshift = -reada_(keywrd_1.keywrd, &i__1, 80L); } if (bshift != 0.) { ten = bshift; } if (i_indx(keywrd_1.keywrd, "ITRY", 80L, 4L) != 0) { i__1 = i_indx(keywrd_1.keywrd, "ITRY", 80L, 4L); itrmax = (integer) reada_(keywrd_1.keywrd, &i__1, 80L); } camkin = i_indx(keywrd_1.keywrd, "KING", 80L, 4L) + i_indx( keywrd_1.keywrd, "CAMP", 80L, 4L) != 0; ci = i_indx(keywrd_1.keywrd, "MICROS", 80L, 6L) + i_indx( keywrd_1.keywrd, "C.I.", 80L, 4L) != 0; okpuly = FALSE_; okpuly = i_indx(keywrd_1.keywrd, "PULAY", 80L, 5L) != 0; uhf = i_indx(keywrd_1.keywrd, "UHF", 80L, 3L) != 0; scf1 = i_indx(keywrd_1.keywrd, "1SCF", 80L, 4L) != 0; oknewd = abs(bshift) < .001; if (camkin && abs(bshift) > 1e-5) { bshift = 4.44; } excitd = i_indx(keywrd_1.keywrd, "EXCITED", 80L, 7L) != 0; times = i_indx(keywrd_1.keywrd, "TIMES", 80L, 5L) != 0; force = i_indx(keywrd_1.keywrd, "FORCE", 80L, 5L) != 0; allcon = okpuly || camkin; /* SET UP C.I. PARAMETERS */ /* NMOS IS NO. OF M.O.S USED IN C.I. */ /* NCIS IS CHANGE IN SPIN, OR NUMBER OF STATES */ if (nmos == 0) { nmos = molkst_1.nopen - molkst_1.nclose; } if (ncis == 0) { if (i_indx(keywrd_1.keywrd, "TRIPLET", 80L, 7L) + i_indx( keywrd_1.keywrd, "QUART", 80L, 5L) != 0) { ncis = 1; } if (i_indx(keywrd_1.keywrd, "QUINTET", 80L, 7L) + i_indx( keywrd_1.keywrd, "SEXTE", 80L, 5L) != 0) { ncis = 2; } } /* DO WE NEED A CAPPED ATOM CORRECTION? */ capps = FALSE_; i__1 = molkst_1.numat; for (i = 1; i <= i__1; ++i) { /* L10: */ if (molkst_1.nat[i - 1] == 102) { capps = TRUE_; } } mesage_1.iiter = 1; trans = .1; if (i_indx(keywrd_1.keywrd, "RESTART", 80L, 7L) + i_indx( keywrd_1.keywrd, "OLDENS", 80L, 6L) != 0) { if (i_indx(keywrd_1.keywrd, "OLDENS", 80L, 6L) != 0) { 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_rsue(&io___48); i__1 = linear; for (i = 1; i <= i__1; ++i) { do_uio(&c__1, (char *)&densty_1.pa[i - 1], (ftnlen)sizeof( doublereal)); } e_rsue(); if (uhf) { s_rsue(&io___49); i__1 = linear; for (i = 1; i <= i__1; ++i) { do_uio(&c__1, (char *)&densty_1.pb[i - 1], (ftnlen)sizeof( doublereal)); } e_rsue(); i__1 = linear; for (i = 1; i <= i__1; ++i) { pold[i - 1] = densty_1.pa[i - 1]; pbold[i - 1] = densty_1.pb[i - 1]; /* L20: */ densty_1.p[i - 1] = densty_1.pa[i - 1] + densty_1.pb[i - 1]; } } else { i__1 = linear; for (i = 1; i <= i__1; ++i) { densty_1.pb[i - 1] = densty_1.pa[i - 1]; pbold[i - 1] = densty_1.pa[i - 1]; pold[i - 1] = densty_1.pa[i - 1]; /* L30: */ densty_1.p[i - 1] = densty_1.pa[i - 1] * 2.; } } } else { numscf_1.nscf = 0; i__1 = linear; for (i = 1; i <= i__1; ++i) { densty_1.p[i - 1] = 0.; densty_1.pa[i - 1] = 0.; /* L40: */ densty_1.pb[i - 1] = 0.; } w1 = na1el / (na1el + 1e-6 + nb1el); w2 = 1. - w1; if (w1 < 1e-6) { w1 = .5; } if (w2 < 1e-6) { w2 = .5; } random = 1.; if (uhf && na1el == nb1el) { random = 1.1; } i__1 = molkst_1.norbs; for (i = 1; i <= i__1; ++i) { j = i * (i + 1) / 2; densty_1.p[j - 1] = molorb_1.pdiag[i - 1]; densty_1.pa[j - 1] = densty_1.p[j - 1] * w1 * random; random = 1. / random; /* L50: */ densty_1.pb[j - 1] = densty_1.p[j - 1] * w2 * random; } i__1 = linear; for (i = 1; i <= i__1; ++i) { pbold[i - 1] = densty_1.pb[i - 1]; /* L60: */ pold[i - 1] = densty_1.pa[i - 1]; } } halfe = molkst_1.nopen != molkst_1.nclose; /* DETERMINE THE SELF-CONSISTENCY CRITERION */ if (halfe || i_indx(keywrd_1.keywrd, "PREC", 80L, 4L) != 0) { scfcrt *= .01; } if (i_indx(keywrd_1.keywrd, "POLAR", 80L, 5L) + i_indx( keywrd_1.keywrd, "NLLSQ", 80L, 5L) + i_indx(keywrd_1.keywrd, "SIGMA", 80L, 5L) != 0) { scfcrt *= .001; } if (force) { scfcrt *= 1e-4; } scfcrt = max(scfcrt,1e-12); /* THE USER CAN STATE THE SCF CRITERION, IF DESIRED. */ i = i_indx(keywrd_1.keywrd, "SCFCRT", 80L, 6L); if (i != 0) { scfcrt = reada_(keywrd_1.keywrd, &i, 80L); s_wsfe(&io___57); do_fio(&c__1, (char *)&scfcrt, (ftnlen)sizeof(doublereal)); e_wsfe(); if (scfcrt < 1e-12) { s_wsfe(&io___58); e_wsfe(); } } else { if (debug) { s_wsfe(&io___59); do_fio(&c__1, (char *)&scfcrt, (ftnlen)sizeof(doublereal)); e_wsfe(); } } last_1.last = 0; /* END OF INITIALIZATION SECTION. */ } else if (force) { /* RESET THE DENSITY MATRIX IF MECI HAS FORMED AN EXCITED STATE */ /* SUM IS NOT USED */ sum = meci_(vector_1.eigs, vector_1.c, vector_1.cbeta, vector_1.eigb, &molkst_1.norbs, &nmos, &c__1, &c_false); } /* INITIALIZATION OPERATIONS DONE EVERY TIME ITER IS CALLED */ makea = TRUE_; makeb = TRUE_; if (newdg) { newdg = abs(bshift) < .001; } if (last_1.last == 1) { newdg = FALSE_; } selcon = scfcrt * 23.061f; if (pltest < scfcrt) { pltest = scfcrt; } if (molkst_1.nalpha != molkst_1.nbeta || ! uhf) { pltest = .001; } if (! force && ! halfe) { if (gradnt_1.gnorm > 5.) { selcon = scfcrt * gradnt_1.gnorm * .2; } if (gradnt_1.gnorm > 200.) { selcon = scfcrt * 50.; } } if (debug) { s_wsfe(&io___64); do_fio(&c__1, (char *)&selcon, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&gradnt_1.gnorm, (ftnlen)sizeof(doublereal)); e_wsfe(); } titer1 = second_(); if (prt1el) { s_wsfe(&io___66); e_wsfe(); vecprt_(&h[1], &molkst_1.norbs); } iredy = 1; L70: niter = 0; time1 = second_(); frst = TRUE_; if (camkin) { modea = 1; modeb = 1; } else { modea = 0; modeb = 0; } bfrst = TRUE_; /* ********************************************************************* */ /* * */ /* * */ /* START THE SCF LOOP HERE * */ /* * */ /* * */ /* ********************************************************************* */ incitr = TRUE_; L80: incitr = modea != 3 && modeb != 3; if (incitr) { ++niter; } if (niter > itrmax - 10 && ! allcon) { /* ****************************************************************** ***** */ /* * */ /* SWITCH ON ALL CONVERGERS * */ /* * */ /* ****************************************************************** ***** */ s_wsfe(&io___75); e_wsfe(); allcon = TRUE_; bshift = 4.44; iredy = -4; eold = 100.; okpuly = TRUE_; newdg = FALSE_; camkin = ! halfe; goto L70; } /* *********************************************************************** */ /* * */ /* MAKE THE ALPHA FOCK MATRIX * */ /* * */ /* *********************************************************************** */ if (abs(shift) > 1e-10 && bshift != 0.) { l = 0; if (niter > 1) { if (newdg && ! (halfe || camkin)) { /* SHIFT WILL APPLY TO THE VIRTUAL ENERGY LEVELS USED IN THE */ /* PSEUDODIAGONALIIZATION. IF DIFF IS -VE, GOOD, THEN LOWER T HE */ /* HOMO-LUMO GAP BY 0.1EV, OTHERWISE INCREASE IT. */ if (diff > 0.) { shift = 1.; /* IF THE PSEUDODIAGONALIZATION APPROXIMATION -- THAT THE WAVEFUNCTION */ /* IS ALMOST STABLE -- IS INVALID, TURN OFF NEWDG */ if (diff > 1.) { newdg = FALSE_; } } else { shift = -.1; } } else { shift = ten + vector_1.eigs[ihomo] - vector_1.eigs[ihomo - 1] + shift; } if (diff > 0.) { if (shift > 4.) { shfmax = 4.5; } if (shift > shfmax) { /* Computing MAX */ d__1 = shfmax - .5; shfmax = max(d__1,0.); } } /* IF SYSTEM GOES UNSTABLE, LIMIT SHIFT TO THE RANGE -INFINITY - SHFMAX */ /* BUT IF SYSTEM IS STABLE, LIMIT SHIFT TO THE RANGE -INFINITY - +20 */ /* Computing MAX */ d__1 = -20., d__2 = min(shfmax,shift); shift = max(d__1,d__2); if ((d__1 = shift - shfmax, abs(d__1)) < 1e-5) { shfmax += .01; } /* THE CAMP-KING AND PULAY CONVERGES NEED A CONSTANT SHIFT. */ /* IF THE SHIFT IS ALLOWED TO VARY, THESE CONVERGERS WILL NOT */ /* WORK PROPERLY. */ if (okpuly || (d__1 = bshift - 4.44, abs(d__1)) < 1e-5) { shift = -8.; if (newdg) { shift = 0.; } } if (uhf) { if (newdg && ! (halfe || camkin)) { shiftb = ten - tenold; } else { shiftb = ten + vector_1.eigs[ihomob] - vector_1.eigs[ ihomob - 1] + shiftb; } if (diff > 0.) { shiftb = min(4.,shiftb); } /* Computing MAX */ d__1 = -20., d__2 = min(shfmax,shiftb); shiftb = max(d__1,d__2); if (okpuly || (d__1 = bshift - 4.44, abs(d__1)) < 1e-5) { shiftb = -8.; if (newdg) { shift = 0.; } } /* # WRITE(6,*)'SHIFT:',SHIFT,SHIFTB */ i__1 = molkst_1.norbs; for (i = ihomob + 1; i <= i__1; ++i) { /* L90: */ vector_1.eigb[i - 1] += shiftb; } } else { /* # WRITE(6,*)'SHIFT:',SHIFT */ } } tenold = ten; i__1 = molkst_1.norbs; for (i = ihomo + 1; i <= i__1; ++i) { /* L100: */ vector_1.eigs[i - 1] += shift; } i__1 = molkst_1.norbs; for (i = 1; i <= i__1; ++i) { i__2 = i; for (j = 1; j <= i__2; ++j) { ++l; /* L110: */ fokmat_1.f[l - 1] = h[l] + shift * densty_1.pa[l - 1]; } /* L120: */ fokmat_1.f[l - 1] -= shift; } } else if (*rand && last_1.last == 0 && niter < 2 && *fulscf) { random = .001; i__1 = linear; for (i = 1; i <= i__1; ++i) { random = -random; /* L130: */ fokmat_1.f[i - 1] = h[i] + random; } } else { i__1 = linear; for (i = 1; i <= i__1; ++i) { /* L140: */ fokmat_1.f[i - 1] = h[i]; } } L150: fock2_(fokmat_1.f, densty_1.p, densty_1.pa, &w[1], &wj[1], &wk[1], & molkst_1.numat, molkst_1.nfirst, molkst_1.nmidle, molkst_1.nlast); fock1_(fokmat_1.f, densty_1.p, densty_1.pa, densty_1.pb); /* *********************************************************************** */ /* * */ /* MAKE THE BETA FOCK MATRIX * */ /* * */ /* *********************************************************************** */ if (uhf) { if (shiftb != 0.) { l = 0; i__1 = molkst_1.norbs; for (i = 1; i <= i__1; ++i) { i__2 = i; for (j = 1; j <= i__2; ++j) { ++l; /* L160: */ fokmat_1.fb[l - 1] = h[l] + shiftb * densty_1.pb[l - 1]; } /* L170: */ fokmat_1.fb[l - 1] -= shiftb; } } else if (*rand && last_1.last == 0 && niter < 2 && *fulscf) { random = .001f; i__1 = linear; for (i = 1; i <= i__1; ++i) { random = -random; /* L180: */ fokmat_1.fb[i - 1] = h[i] + random; } } else { i__1 = linear; for (i = 1; i <= i__1; ++i) { /* L190: */ fokmat_1.fb[i - 1] = h[i]; } } fock2_(fokmat_1.fb, densty_1.p, densty_1.pb, &w[1], &wj[1], &wk[1], & molkst_1.numat, molkst_1.nfirst, molkst_1.nmidle, molkst_1.nlast); fock1_(fokmat_1.fb, densty_1.p, densty_1.pb, densty_1.pa); } if (! (*fulscf)) { goto L290; } if (prtfok) { s_wsfe(&io___80); do_fio(&c__1, (char *)&niter, (ftnlen)sizeof(integer)); e_wsfe(); vecprt_(fokmat_1.f, &molkst_1.norbs); } /* *********************************************************************** */ /* * */ /* CALCULATE THE ENERGY IN KCAL/MOLE * */ /* * */ /* *********************************************************************** */ *ee = helect_(&molkst_1.norbs, densty_1.pa, &h[1], fokmat_1.f); if (uhf) { *ee += helect_(&molkst_1.norbs, densty_1.pb, &h[1], fokmat_1.fb); } else { *ee *= 2.; } if (capps) { *ee += capcor_(molkst_1.nat, molkst_1.nfirst, molkst_1.nlast, & molkst_1.numat, densty_1.p, &h[1]); } scorr = shift * (molkst_1.nopen - molkst_1.nclose) * 23.061 * .25 * ( molkst_1.fract * (2. - molkst_1.fract)); escf = (*ee + enuclr_1.enuclr) * 23.061 + atheat_1.atheat + scorr; /* Computing MAX */ /* Computing MAX */ d__3 = abs(*ee); d__1 = selcon, d__2 = eps * max(d__3,1.); sellim = max(d__1,d__2); if (incitr) { diff = escf - eold; if (diff > 0.) { ten += -1.; } else { ten = ten * .975 + .05; } /* # WRITE(6,'(2F12.6)')TEN,DIFF */ /* Computing MAX */ d__2 = 1., d__3 = abs(*ee); if (niter > 4 && (pl == 0. || pl < pltest && (d__1 = diff / max(d__2, d__3), abs(d__1)) < sellim) && ready) { /* ************************************************************** ********* */ /* * */ /* SELF-CONSISTENCY TEST, EXIT MODE FROM ITERATIONS * */ /* * */ /* ************************************************************** ********* */ if (abs(shift) < 1e-10) { goto L290; } shift = 0.; shiftb = 0.; i__1 = linear; for (i = 1; i <= i__1; ++i) { /* L210: */ fokmat_1.f[i - 1] = h[i]; } makea = TRUE_; makeb = TRUE_; goto L150; } ready = iredy > 0 && (abs(diff) < sellim * 10. || pl == 0.); ++iredy; } if (prtpl || debug && niter > itrmax - 20) { if (abs(escf) > 99999.) { escf = d_sign(&c_b87, &escf); } if (abs(diff) > 9999.) { diff = 0.; } if (incitr) { s_wsfe(&io___84); do_fio(&c__1, (char *)&niter, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&pl, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&plb, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&escf, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&diff, (ftnlen)sizeof(doublereal)); e_wsfe(); } } if (incitr) { eold = escf; } /* *********************************************************************** */ /* * */ /* INVOKE THE CAMP-KING CONVERGER * */ /* * */ /* *********************************************************************** */ if (niter > 2 && camkin && makea) { i__1 = molkst_1.norbs - na1el; d__1 = escf / 23.061; interp_(&molkst_1.norbs, &na1el, &i__1, &modea, &d__1, fokmat_1.f, vector_1.c, ar1, ar2, ar3, ar4, ar1); } makeb = FALSE_; if (modea == 3) { goto L230; } makeb = TRUE_; if (newdg) { /* ****************************************************************** ***** */ /* * */ /* INVOKE PULAY'S CONVERGER * */ /* * */ /* ****************************************************************** ***** */ if (okpuly && makea && iredy > 1) { pulay_(fokmat_1.f, densty_1.pa, &molkst_1.norbs, pold, pold2, pold3, &jalp, &ialp, &c__23220, &frst, &pl); } /* ****************************************************************** ***** */ /* * */ /* DIAGONALIZE THE ALPHA OR RHF SECULAR DETERMINANT * */ /* WHERE POSSIBLE, USE THE PULAY-STEWART METHOD, OTHERWISE USE BEPPU'S * */ /* * */ /* ****************************************************************** ***** */ if (halfe || camkin) { rsp_(fokmat_1.f, &molkst_1.norbs, &molkst_1.norbs, vector_1.eigs, vector_1.c); } else { diag_(fokmat_1.f, vector_1.c, &na1el, vector_1.eigs, & molkst_1.norbs, &molkst_1.norbs); } } else { rsp_(fokmat_1.f, &molkst_1.norbs, &molkst_1.norbs, vector_1.eigs, vector_1.c); } j = 1; if (prtvec) { j = 1; if (uhf) { j = 2; } s_wsfe(&io___94); do_fio(&c__1, abprt + (j - 1) * 5, 5L); do_fio(&c__1, (char *)&niter, (ftnlen)sizeof(integer)); e_wsfe(); matout_(vector_1.c, vector_1.eigs, &molkst_1.norbs, &molkst_1.norbs, & molkst_1.norbs); } else { if (prteig) { s_wsfe(&io___95); do_fio(&c__1, abprt + (j - 1) * 5, 5L); do_fio(&c__1, (char *)&niter, (ftnlen)sizeof(integer)); i__1 = molkst_1.norbs; for (i = 1; i <= i__1; ++i) { do_fio(&c__1, (char *)&vector_1.eigs[i - 1], (ftnlen)sizeof( doublereal)); } e_wsfe(); } } L230: if (ifill != 0) { swap_(vector_1.c, &molkst_1.norbs, &molkst_1.norbs, &na2el, &ifill); } /* *********************************************************************** */ /* * */ /* CALCULATE THE ALPHA OR RHF DENSITY MATRIX * */ /* * */ /* *********************************************************************** */ if (uhf) { densit_(vector_1.c, &molkst_1.norbs, &molkst_1.norbs, &na2el, &na1el, &molkst_1.fract, densty_1.pa, &c__1); if (modea != 3 && ! (newdg && okpuly)) { cnvg_(densty_1.pa, pold, pold2, &molkst_1.norbs, &niter, &pl); } } else { densit_(vector_1.c, &molkst_1.norbs, &molkst_1.norbs, &na2el, &na1el, &molkst_1.fract, densty_1.p, &c__1); if (modea != 3 && ! (newdg && okpuly)) { cnvg_(densty_1.p, pold, pold2, &molkst_1.norbs, &niter, &pl); } } /* *********************************************************************** */ /* * */ /* UHF-SPECIFIC CODE * */ /* * */ /* *********************************************************************** */ if (uhf) { /* ****************************************************************** ***** */ /* * */ /* INVOKE THE CAMP-KING CONVERGER * */ /* * */ /* ****************************************************************** ***** */ if (niter > 2 && camkin && makeb) { i__1 = molkst_1.norbs - nb1el; d__1 = escf / 23.061; interp_(&molkst_1.norbs, &nb1el, &i__1, &modeb, &d__1, fokmat_1.fb, vector_1.cbeta, br1, br2, br3, br4, br1); } makea = FALSE_; if (modeb == 3) { goto L240; } makea = TRUE_; if (newdg) { /* ************************************************************** ********* */ /* * */ /* INVOKE PULAY'S CONVERGER * */ /* * */ /* ************************************************************** ********* */ if (okpuly && makeb && iredy > 1) { pulay_(fokmat_1.fb, densty_1.pb, &molkst_1.norbs, pbold, pbold2, pbold3, &jbet, &ibet, &c__23220, &bfrst, &plb) ; } /* ************************************************************** ********* */ /* * */ /* DIAGONALIZE THE ALPHA OR RHF SECULAR DETERMINANT * */ /* WHERE POSSIBLE, USE THE PULAY-STEWART METHOD, OTHERWISE USE BEP PU'S * */ /* * */ /* ************************************************************** ********* */ if (halfe || camkin) { rsp_(fokmat_1.fb, &molkst_1.norbs, &molkst_1.norbs, vector_1.eigb, vector_1.cbeta); } else { diag_(fokmat_1.fb, vector_1.cbeta, &nb1el, vector_1.eigb, & molkst_1.norbs, &molkst_1.norbs); } } else { rsp_(fokmat_1.fb, &molkst_1.norbs, &molkst_1.norbs, vector_1.eigb, vector_1.cbeta); } if (prtvec) { s_wsfe(&io___104); do_fio(&c__1, abprt + 10, 5L); do_fio(&c__1, (char *)&niter, (ftnlen)sizeof(integer)); e_wsfe(); matout_(vector_1.cbeta, vector_1.eigb, &molkst_1.norbs, & molkst_1.norbs, &molkst_1.norbs); } else { if (prteig) { s_wsfe(&io___105); do_fio(&c__1, abprt + 10, 5L); do_fio(&c__1, (char *)&niter, (ftnlen)sizeof(integer)); i__1 = molkst_1.norbs; for (i = 1; i <= i__1; ++i) { do_fio(&c__1, (char *)&vector_1.eigb[i - 1], (ftnlen) sizeof(doublereal)); } e_wsfe(); } } /* ****************************************************************** ***** */ /* * */ /* CALCULATE THE BETA DENSITY MATRIX * */ /* * */ /* ****************************************************************** ***** */ L240: densit_(vector_1.cbeta, &molkst_1.norbs, &molkst_1.norbs, &nb2el, & nb1el, &molkst_1.fract, densty_1.pb, &c__1); if (! (newdg && okpuly)) { cnvg_(densty_1.pb, pbold, pbold2, &molkst_1.norbs, &niter, &plb); } } /* *********************************************************************** */ /* * */ /* CALCULATE THE TOTAL DENSITY MATRIX * */ /* * */ /* *********************************************************************** */ if (uhf) { i__1 = linear; for (i = 1; i <= i__1; ++i) { /* L250: */ densty_1.p[i - 1] = densty_1.pa[i - 1] + densty_1.pb[i - 1]; } } else { i__1 = linear; for (i = 1; i <= i__1; ++i) { densty_1.pa[i - 1] = densty_1.p[i - 1] * .5; /* L260: */ densty_1.pb[i - 1] = densty_1.pa[i - 1]; } } if (prtden) { s_wsfe(&io___106); do_fio(&c__1, (char *)&niter, (ftnlen)sizeof(integer)); e_wsfe(); vecprt_(densty_1.p, &molkst_1.norbs); } oknewd = pl < sellim || oknewd; newdg = pl < trans && oknewd || newdg; if (pl < trans * .3333) { oknewd = TRUE_; } if (niter >= itrmax) { if (diff < .001 && pl < 1e-4 && ! force) { s_wsfe(&io___107); e_wsfe(); goto L290; } if (minprt) { s_wsfe(&io___108); e_wsfe(); } s_wsfe(&io___109); do_fio(&c__1, (char *)&diff, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&pl, (ftnlen)sizeof(doublereal)); e_wsfe(); mesage_1.iflepo = 9; mesage_1.iiter = 2; write_(&time_1.time0, &escf); s_stop("", 0L); } goto L80; /* ********************************************************************* */ /* * */ /* * */ /* END THE SCF LOOP HERE * */ /* NOW CALCULATE THE ELECTRONIC ENERGY * */ /* * */ /* * */ /* ********************************************************************* */ /* SELF-CONSISTENCE ACHEIVED. */ L290: *ee = helect_(&molkst_1.norbs, densty_1.pa, &h[1], fokmat_1.f); if (uhf) { *ee += helect_(&molkst_1.norbs, densty_1.pb, &h[1], fokmat_1.fb); } else { *ee = *ee * 2. + shift * (molkst_1.nopen - molkst_1.nclose) * 23.061 * .25 * (molkst_1.fract * (2. - molkst_1.fract)); } if (capps) { *ee += capcor_(molkst_1.nat, molkst_1.nfirst, molkst_1.nlast, & molkst_1.numat, densty_1.p, &h[1]); } /* NORMALLY THE EIGENVALUES ARE INCORRECT BECAUSE THE */ /* PSEUDODIAGONALIZATION HAS BEEN USED. IF THIS */ /* IS THE LAST SCF, THEN DO AN EXACT DIAGONALIZATION */ if (numscf_1.nscf == 0 || last_1.last == 1 || ci || halfe) { /* PUT F AND FB INTO POLD IN ORDER TO NOT DESTROY F AND FB */ /* AND DO EXACT DIAGONALISATIONS */ i__1 = linear; for (i = 1; i <= i__1; ++i) { /* L300: */ pold[i - 1] = fokmat_1.f[i - 1]; } rsp_(pold, &molkst_1.norbs, &molkst_1.norbs, vector_1.eigs, vector_1.c); if (uhf) { i__1 = linear; for (i = 1; i <= i__1; ++i) { /* L310: */ pold[i - 1] = fokmat_1.fb[i - 1]; } rsp_(pold, &molkst_1.norbs, &molkst_1.norbs, vector_1.eigb, vector_1.cbeta); } i__1 = linear; for (i = 1; i <= i__1; ++i) { /* L320: */ pold[i - 1] = densty_1.p[i - 1]; } if (ci || halfe) { sum = meci_(vector_1.eigs, vector_1.c, vector_1.cbeta, vector_1.eigb, &molkst_1.norbs, &nmos, &c__0, &c_false); *ee += sum; if (prtpl) { escf = (*ee + enuclr_1.enuclr) * 23.061f + atheat_1.atheat; s_wsfe(&io___110); do_fio(&c__1, (char *)&escf, (ftnlen)sizeof(doublereal)); e_wsfe(); } } } ++numscf_1.nscf; titer2 = second_(); if (times) { s_wsfe(&io___112); d__1 = titer2 - titer1; do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof(doublereal)); d__2 = titer2 - time_1.time0; do_fio(&c__1, (char *)&d__2, (ftnlen)sizeof(doublereal)); e_wsfe(); } if (debug) { s_wsfe(&io___113); do_fio(&c__1, (char *)&niter, (ftnlen)sizeof(integer)); e_wsfe(); } /* IF(FORCE) SCFCRT=1.D-5 */ if (allcon && (d__1 = bshift - 4.44, abs(d__1)) < 1e-7) { camkin = FALSE_; allcon = FALSE_; newdg = FALSE_; bshift = -10.; okpuly = FALSE_; } shift = 1.; return 0; } /* iter_ */