/* meci.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 { integer numat, nat[86], nfirst[86], nmidle[86], nlast[86], norbs, nelecs, ndummy[2], nclose, nopen, ndumy; doublereal fract; } molkst_; #define molkst_1 molkst_ struct { integer ispqr[1000] /* was [100][10] */, is, i, k; } spqr_; #define spqr_1 spqr_ struct { doublereal p[23220], pa[23220], pb[23220]; } densty_; #define densty_1 densty_ struct { char keywrd[80]; } keywrd_; #define keywrd_1 keywrd_ struct { doublereal occa[10]; } baseoc_; #define baseoc_1 baseoc_ struct { doublereal xy[10000] /* was [10][10][10][10] */; } xyijkl_; #define xyijkl_1 xyijkl_ /* Table of constant values */ static integer c__1 = 1; static integer c__10 = 10; static integer c__9 = 9; doublereal meci_(doublereal *eigs, doublereal *coeff, doublereal *coeffs, doublereal *eiga, integer *n, integer *nmos, integer *ireset, logical *finish) { /* Initialized data */ static char tspin[8*7] = "SINGLET " "DOUBLET " "TRIPLET " "QUARTET " "QUINTET " "SEXTET " "SEPTET "; static logical first = TRUE_; static logical first1 = TRUE_; /* Format strings */ static char fmt_150[] = "(/,5x,\002EIGENVALUES AFTER REMOVAL OF INTER-EL" "ECTRONIC INTERACTIONS\002,/)"; static char fmt_250[] = "(10x,\002COMPONENT OF SPIN = \002,f4.1)"; static char fmt_260[] = "(/10x,\002 IMPOSSIBLE VALUE OF DELTA S\002/)"; static char fmt_280[] = "(//10x,\002 NO OF CONFIGURATIONS CONSIDERED =" " \002,i4)"; static char fmt_440[] = "(///,\002 STATE ENERGIES \002,\002 EXPECTATION " "VALUE OF S**2 S FROM S**2=S(S+1)\002,//)"; static char fmt_510[] = "(i5,f12.6,3x,a8,f15.5,f10.5)"; /* System generated locals */ integer coeff_dim1, coeff_offset, coeffs_dim1, coeffs_offset, i__1, i__2, i__3, i__4, i__5; doublereal ret_val, d__1, d__2; olist o__1; alist al__1; /* 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_open(olist *), f_rew(alist *), s_rsfe(cilist *), e_rsfe(void), s_cmp(char *, char *, ftnlen, ftnlen); /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen); double sqrt(doublereal); integer i_dnnt(doublereal *), s_wsle(cilist *), e_wsle(void), do_lio( integer *, integer *, char *, ftnlen); /* Local variables */ static doublereal diag[2000]; static integer lima, limb; static doublereal conf[10000]; static integer mdim; static char line[80]; static logical doub; extern /* Subroutine */ int ijkl_(integer *, integer *, integer *, integer *, doublereal *, doublereal *, integer *); static logical sing; extern /* Subroutine */ int perm_(integer *, integer *, integer *, integer *, integer *); static logical quar; static doublereal spin[100]; static logical quin, trip; static doublereal summ; static integer nupp; static logical prnt, sext; static integer ntot; static logical last1, prnt2; static integer j; extern doublereal reada_(char *, integer *, ftnlen); static integer l; extern doublereal diagi_(integer *, integer *, doublereal *, doublereal *, integer *); static logical debug; static integer nelec; static doublereal rjkaa[100] /* was [10][10] */, rjkab[100] /* was [10][10] */; static logical large; static doublereal cimat[10000]; static logical force; static doublereal x; static integer kbeta; static doublereal y; static logical geook; static integer nleft; static logical lspin; static integer ndown, izero, i1, j1, iocca1[10], iocca2[10], ioccb1[10], ioccb2[10], lroot; static doublereal smult; static integer k1, l1, kroot, iroot; extern doublereal aababc_(integer *, integer *, integer *, integer *, integer *), aabacd_(integer *, integer *, integer *, integer *, integer *), aabbcd_(integer *, integer *, integer *, integer *, integer *), babbcd_(integer *, integer *, integer *, integer *, integer *), babbbc_(integer *, integer *, integer *, integer *, integer *); static logical lspin1; static integer ne, ik, ii, ji, in, li, il, iu, kalpha, ix, iy, nalpha[100] , microa[20000] /* was [10][2000] */, microb[20000] /* was [10][2000] */, nperma[600] /* was [10][60] */, npermb[ 600] /* was [10][60] */; static logical bigprt; static doublereal xx; static integer maxvec, iofset, natoms; extern /* Subroutine */ int vecprt_(doublereal *, integer *), matout_( doublereal *, doublereal *, integer *, integer *, integer *); static integer lab; static doublereal eig[100]; static integer nfa[20]; static doublereal gse, ams, one; static integer iuj; extern /* Subroutine */ int rsp_(doublereal *, integer *, integer *, doublereal *, doublereal *); static doublereal sum; /* Fortran I/O blocks */ static cilist io___32 = { 0, 6, 0, "(' INITIAL EIGENVALUES')", 0 }; static cilist io___33 = { 0, 6, 0, "(5F12.6)", 0 }; static cilist io___34 = { 0, 6, 0, "(//10X,'NUMBER OF ELECTRONS IN C.I. " "=',F5.1)", 0 }; static cilist io___35 = { 0, 6, 0, "(///10X,A)", 0 }; static cilist io___36 = { 0, 6, 0, "(10X,A)", 0 }; static cilist io___37 = { 0, 6, 0, "(10X,A)", 0 }; static cilist io___38 = { 0, 6, 0, "(10X,A)", 0 }; static cilist io___39 = { 0, 6, 0, "(10X,A)", 0 }; static cilist io___42 = { 0, 6, 0, "(//10X,'EIGENVECTORS',/)", 0 }; static cilist io___43 = { 0, 6, 0, "(6F12.6)", 0 }; static cilist io___51 = { 0, 6, 0, fmt_150, 0 }; static cilist io___52 = { 0, 6, 0, "(6F12.6)", 0 }; static cilist io___53 = { 0, 6, 0, "(///10X,'TWO-ELECTRON J-INTEGRALS',/)" , 0 }; static cilist io___54 = { 0, 6, 0, "(10F10.4)", 0 }; static cilist io___55 = { 0, 6, 0, "(///10X,'TWO-ELECTRON K-INTEGRALS',/)" , 0 }; static cilist io___56 = { 0, 6, 0, "(10F10.4)", 0 }; static cilist io___59 = { 0, 6, 0, "(' MICROSTATES READ IN')", 0 }; static cilist io___61 = { 0, 5, 0, "(A)", 0 }; static cilist io___63 = { 0, 5, 0, "(A)", 0 }; static cilist io___64 = { 0, 5, 0, "(A)", 0 }; static cilist io___68 = { 0, 6, 0, "(20I6)", 0 }; static cilist io___69 = { 0, 6, 0, "(/,'NUMBER OF ELECTRONS IN C.I. REDE" "FINED TO:',I4,/)", 0 }; static cilist io___73 = { 0, 6, 0, fmt_250, 0 }; static cilist io___74 = { 0, 6, 0, fmt_260, 0 }; static cilist io___77 = { 0, 6, 0, fmt_280, 0 }; static cilist io___81 = { 0, 6, 0, "(//10X,'CONFIGURATIONS CONSIDERED IN" " C.I.',//)", 0 }; static cilist io___83 = { 0, 6, 0, "(//10X,'C.I. IS OF SIZE LESS THAN RO" "OT SPECIFIED')", 0 }; static cilist io___84 = { 0, 6, 0, "(10X,'MODIFY SIZE OF C.I. OR ROOT NU" "MBER')", 0 }; static cilist io___89 = { 0, 6, 0, "(/10X,I4,6X,10I4)", 0 }; static cilist io___90 = { 0, 6, 0, "(20X,10I4)", 0 }; static cilist io___97 = { 0, 6, 0, "(//,' C.I. MATRIX')", 0 }; static cilist io___98 = { 0, 6, 0, "(//,' DIAGONAL OF C.I. MATRIX')", 0 }; static cilist io___99 = { 0, 6, 0, "(5F13.6)", 0 }; static cilist io___103 = { 0, 6, 0, "(//20X,'STATE VECTORS',//)", 0 }; static cilist io___104 = { 0, 6, 0, fmt_440, 0 }; static cilist io___110 = { 0, 6, 0, fmt_510, 0 }; static cilist io___111 = { 0, 6, 0, "(//10X,'THE STATE REQUIRED IS NOT P" "RESENT IN THE')", 0 }; static cilist io___112 = { 0, 6, 0, "(10X, ' SET OF CONFIGURATIONS A" "VAILABLE')", 0 }; static cilist io___113 = { 0, 6, 0, "(/ 4X,'NUMBER OF STATES ACCESSIBLE " "USING CURRENT KEY-WORDS',/)", 0 }; static cilist io___114 = { 0, 6, 0, "((24X,A8,I4))", 0 }; static cilist io___117 = { 0, 6, 0, "(' ESR SPECIFIED FOR AN EVEN-ELEC" "TRON SYSTEM')", 0 }; static cilist io___120 = { 0, 6, 0, "(//,' MICROSTATE CONTRIBUTIONS" " TO ', 'STATE EIGENFUNCTION',I3)", 0 }; static cilist io___121 = { 0, 6, 0, "(5F13.6)", 0 }; static cilist io___122 = { 0, 6, 0, "(/,' SPIN DENSITIES FROM EACH M." "O., ENERGY:' ,F7.3)", 0 }; static cilist io___123 = { 0, 6, 0, "(5F12.6)", 0 }; static cilist io___124 = { 0, 6, 0, 0, 0 }; static cilist io___125 = { 0, 6, 0, 0, 0 }; static cilist io___126 = { 0, 6, 0, 0, 0 }; static cilist io___130 = { 0, 6, 0, "(' ATOM',I4,' SPIN DENSITY ',5" "F10.7)", 0 }; static cilist io___131 = { 0, 6, 0, "(' ATOM',I4,' SPIN DENSITY ',F" "10.7,30X,F10.7)", 0 }; /* ********************************************************************** */ /* PROGRAM MECI */ /* A MULTI-ELECTRON CONFIGURATION INTERACTION CALCULATION */ /* WRITTEN BY JAMES J. P. STEWART, AT THE */ /* FRANK J. SEILER RESEARCH LABORATORY */ /* USAFA, COLORADO SPRINGS, CO 80840 */ /* 1985 */ /* ********************************************************************** */ /* 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 */ /* MATRICES FOR SEC.DET., VECTORS, AND EIGENVALUES. */ /* MATRICES TO HOLD ELECTRON CONFIGURATIONS */ /* MATRICES FOR PERMUTATION WORK */ /* MATRICES FOR ONE AND TWO ELECTRON INTEGRALS, AND M.O.S */ /* SPIN MATRICES */ /* Parameter adjustments */ --eigs; --eiga; coeffs_dim1 = *n; coeffs_offset = coeffs_dim1 + 1; coeffs -= coeffs_offset; coeff_dim1 = *n; coeff_offset = coeff_dim1 + 1; coeff -= coeff_offset; /* Function Body */ if (first) { first = FALSE_; geook = i_indx(keywrd_1.keywrd, "GEO-OK", 80L, 6L) != 0; force = i_indx(keywrd_1.keywrd, "FORCE", 80L, 5L) != 0; last1 = i_indx(keywrd_1.keywrd, "1SCF", 80L, 4L) != 0; lspin1 = i_indx(keywrd_1.keywrd, "ESR", 80L, 3L) != 0; mdim = 100; debug = i_indx(keywrd_1.keywrd, "DEBUG", 80L, 5L) != 0; prnt2 = i_indx(keywrd_1.keywrd, "MECI", 80L, 4L) != 0; debug = debug && prnt2; large = i_indx(keywrd_1.keywrd, "LARGE", 80L, 5L) != 0; lroot = 1; if (i_indx(keywrd_1.keywrd, "EXCI", 80L, 4L) != 0) { lroot = 2; } spqr_1.i = i_indx(keywrd_1.keywrd, "ROOT", 80L, 4L); if (spqr_1.i != 0) { lroot = (integer) reada_(keywrd_1.keywrd, &spqr_1.i, 80L); } /* # WRITE(6,'('' NCLOSE, NOPEN'',2I6)')NCLOSE,NOPEN */ /* Computing MAX */ /* Computing MIN */ i__2 = (molkst_1.nclose + molkst_1.nopen + 1) / 2 - (*nmos - 1) / 2, i__3 = molkst_1.norbs - *nmos + 1; i__1 = min(i__2,i__3); j = max(i__1,1); l = 0; i__1 = molkst_1.nclose; for (spqr_1.i = j; spqr_1.i <= i__1; ++spqr_1.i) { ++l; /* L10: */ baseoc_1.occa[l - 1] = 1.; } i__1 = molkst_1.nopen; for (spqr_1.i = molkst_1.nclose + 1; spqr_1.i <= i__1; ++spqr_1.i) { ++l; /* L20: */ baseoc_1.occa[l - 1] = molkst_1.fract * .5; } i__1 = j + *nmos - 1; for (spqr_1.i = molkst_1.nopen + 1; spqr_1.i <= i__1; ++spqr_1.i) { ++l; /* L30: */ baseoc_1.occa[l - 1] = 0.; } /* # WRITE(6,'('' INITIAL ORBITAL OCCUPANCIES'')') */ /* # WRITE(6,'(6F12.6)')(OCCA(L),L=1,NMOS) */ sing = i_indx(keywrd_1.keywrd, "SINGL", 80L, 5L) + i_indx( keywrd_1.keywrd, "EXCI", 80L, 4L) + i_indx(keywrd_1.keywrd, "BIRAD", 80L, 5L) != 0; doub = i_indx(keywrd_1.keywrd, "DOUBL", 80L, 5L) != 0; trip = i_indx(keywrd_1.keywrd, "TRIPL", 80L, 5L) != 0; quar = i_indx(keywrd_1.keywrd, "QUART", 80L, 5L) != 0; quin = i_indx(keywrd_1.keywrd, "QUINT", 80L, 5L) != 0; sext = i_indx(keywrd_1.keywrd, "SEXTE", 80L, 5L) != 0; smult = -.5; if (sing) { smult = 0.; } if (doub) { smult = .75; } if (trip) { smult = 2.; } if (quar) { smult = 3.75; } if (quin) { smult = 6.; } if (sext) { smult = 8.75; } /* # WRITE(6,'('' ORBITAL COUNTERS, PER ATOM, FIRST LAST'')') */ /* # WRITE(6,'(I40,I6)')(NFIRST(I),NLAST(I),I=1,NUMAT) */ x = 0.; i__1 = *nmos; for (j = 1; j <= i__1; ++j) { /* L40: */ x += baseoc_1.occa[j - 1]; } xx = x + x; ne = (integer) (xx + .5f); nelec = (molkst_1.nelecs - ne + 1) / 2; nleft = molkst_1.norbs - *nmos - nelec; } if (*ireset == 1) { one = -1.; goto L520; } else { one = 1.; } prnt = debug || *finish && prnt2; bigprt = prnt && large; last1 = last1 || *finish; /* TEST TO SEE IF THE SET OF ENERGY LEVELS USED IN MECI IS COMPLETE, */ /* I.E., ALL COMPONENTS OF DEGENERATE IRREDUCIBLE REPRESENTATIONS */ /* ARE USED. IF NOT, THEN RESULTS WILL BE NONSENSE. GIVE USERS A */ /* CHANCE TO REALLY FOUL THINGS UP BY ALLOWING JOB TO CONTINUE IF */ /* 'GEO-OK' IS SPECIFIED. */ i__1 = *nmos; for (spqr_1.i = 1; spqr_1.i <= i__1; ++spqr_1.i) { in = spqr_1.i + nelec; i__2 = molkst_1.norbs; for (j = 1; j <= i__2; ++j) { /* L50: */ coeffs[j + spqr_1.i * coeffs_dim1] = coeff[j + in * coeff_dim1]; } /* L60: */ eiga[spqr_1.i] = eigs[in]; } lspin = lspin1 && *finish; if (bigprt) { s_wsfe(&io___32); e_wsfe(); s_wsfe(&io___33); i__1 = *nmos; for (spqr_1.i = 1; spqr_1.i <= i__1; ++spqr_1.i) { do_fio(&c__1, (char *)&eiga[spqr_1.i], (ftnlen)sizeof(doublereal)) ; } e_wsfe(); s_wsfe(&io___34); do_fio(&c__1, (char *)&xx, (ftnlen)sizeof(doublereal)); e_wsfe(); } if (! geook && nelec > 0) { if ((d__1 = eigs[nelec + 1] - eigs[nelec], abs(d__1)) < .1 || (d__2 = eigs[nelec + 1 + *nmos] - eigs[nelec + *nmos], abs(d__2)) < .1) { s_wsfe(&io___35); do_fio(&c__1, "DEGENERATE ENERGY LEVELS DETECTED IN MECI", 41L); e_wsfe(); s_wsfe(&io___36); do_fio(&c__1, "SOME OF THESE LEVELS WOULD BE TREATED BYMECI,", 45L); e_wsfe(); s_wsfe(&io___37); do_fio(&c__1, "WHILE OTHERS WOULD NOT. THIS WOULD RESULT IN", 45L); e_wsfe(); s_wsfe(&io___38); do_fio(&c__1, "NON-REPRODUCIBLE ELECTRONIC ENERGIES.", 37L); e_wsfe(); s_wsfe(&io___39); do_fio(&c__1, " JOB STOPPED. TO CONTINUE, SPECIFY \"GEO-OK\"", 45L); e_wsfe(); s_stop("", 0L); } } kalpha = (ne + 1) / 2; kbeta = ne - kalpha; /* FOR DEBUGGING PURPOSES, UNCOMMENT THE FOLLOWING LINES */ /* THEY LOCK THE PHASE-FACTOR */ /* # DO 55 I=1,NMOS */ /* # SUM=0.D0 */ /* # DO 56 J=1,NORBS */ /* # 56 IF(ABS(SUM).LT.ABS(COEFFS(J,I)))SUM=COEFFS(J,I) */ /* # IF(SUM.LT.0.D0)THEN */ /* # DO 57 J=1,NORBS */ /* # 57 COEFFS(J,I)=-COEFFS(J,I) */ /* # ENDIF */ /* # 55 CONTINUE */ /* FOR DEBUGGING PURPOSES, UNCOMMENT THE PRECEEDING LINES */ if (bigprt) { s_wsfe(&io___42); e_wsfe(); i__1 = molkst_1.norbs; for (spqr_1.i = 1; spqr_1.i <= i__1; ++spqr_1.i) { /* L70: */ s_wsfe(&io___43); i__2 = *nmos; for (j = 1; j <= i__2; ++j) { do_fio(&c__1, (char *)&coeffs[spqr_1.i + j * coeffs_dim1], ( ftnlen)sizeof(doublereal)); } e_wsfe(); } } i__2 = *nmos; for (spqr_1.i = 1; spqr_1.i <= i__2; ++spqr_1.i) { i__1 = *nmos; for (j = 1; j <= i__1; ++j) { i__3 = *nmos; for (spqr_1.k = 1; spqr_1.k <= i__3; ++spqr_1.k) { i__4 = *nmos; for (l = 1; l <= i__4; ++l) { /* L80: */ xyijkl_1.xy[spqr_1.i + (j + (spqr_1.k + l * 10) * 10) * 10 - 1111] = 100.; } } } } nfa[1] = 1; nfa[0] = 1; for (spqr_1.i = 3; spqr_1.i <= 11; ++spqr_1.i) { /* L90: */ nfa[spqr_1.i - 1] = nfa[spqr_1.i - 2] * (spqr_1.i - 1); } i__4 = *nmos; for (spqr_1.i = 1; spqr_1.i <= i__4; ++spqr_1.i) { i__3 = *nmos; for (j = 1; j <= i__3; ++j) { i__1 = *nmos; for (spqr_1.k = 1; spqr_1.k <= i__1; ++spqr_1.k) { i__2 = *nmos; for (l = 1; l <= i__2; ++l) { /* L100: */ xyijkl_1.xy[spqr_1.i + (j + (spqr_1.k + l * 10) * 10) * 10 - 1111] = 100.; } } } } i__2 = *nmos; for (spqr_1.i = 1; spqr_1.i <= i__2; ++spqr_1.i) { i1 = spqr_1.i; i__1 = spqr_1.i; for (j = 1; j <= i__1; ++j) { j1 = j; i__3 = *nmos; for (spqr_1.k = 1; spqr_1.k <= i__3; ++spqr_1.k) { k1 = spqr_1.k; i__4 = spqr_1.k; for (l = 1; l <= i__4; ++l) { l1 = l; /* L110: */ ijkl_(&i1, &k1, &j1, &l1, &x, &coeffs[coeffs_offset], & molkst_1.norbs); } } } } i__4 = *nmos; for (spqr_1.i = 1; spqr_1.i <= i__4; ++spqr_1.i) { i__3 = *nmos; for (j = 1; j <= i__3; ++j) { rjkaa[spqr_1.i + j * 10 - 11] = xyijkl_1.xy[spqr_1.i + (spqr_1.i + (j + j * 10) * 10) * 10 - 1111] - xyijkl_1.xy[spqr_1.i + (j + (spqr_1.i + j * 10) * 10) * 10 - 1111]; /* L120: */ rjkab[spqr_1.i + j * 10 - 11] = xyijkl_1.xy[spqr_1.i + (spqr_1.i + (j + j * 10) * 10) * 10 - 1111]; } } i__3 = *nmos; for (spqr_1.i = 1; spqr_1.i <= i__3; ++spqr_1.i) { x = 0.f; i__4 = *nmos; for (j = 1; j <= i__4; ++j) { x += (rjkaa[spqr_1.i + j * 10 - 11] + rjkab[spqr_1.i + j * 10 - 11]) * baseoc_1.occa[j - 1]; /* L130: */ } eiga[spqr_1.i] -= x; /* # IF(ABS(OCCA(I)-0.5).LT.1.D-4)EIGA(I)=EIGA(I)+XY(I,I,I,I)*0.2 5D0 */ /* L140: */ } if (bigprt) { s_wsfe(&io___51); e_wsfe(); s_wsfe(&io___52); i__3 = *nmos; for (spqr_1.i = 1; spqr_1.i <= i__3; ++spqr_1.i) { do_fio(&c__1, (char *)&eiga[spqr_1.i], (ftnlen)sizeof(doublereal)) ; } e_wsfe(); s_wsfe(&io___53); e_wsfe(); i__3 = *nmos; for (i1 = 1; i1 <= i__3; ++i1) { /* L160: */ s_wsfe(&io___54); i__4 = *nmos; for (j1 = 1; j1 <= i__4; ++j1) { do_fio(&c__1, (char *)&rjkab[i1 + j1 * 10 - 11], (ftnlen) sizeof(doublereal)); } e_wsfe(); } s_wsfe(&io___55); e_wsfe(); i__4 = *nmos; for (i1 = 1; i1 <= i__4; ++i1) { /* L170: */ s_wsfe(&io___56); i__3 = *nmos; for (j1 = 1; j1 <= i__3; ++j1) { d__1 = rjkab[i1 + j1 * 10 - 11] - rjkaa[i1 + j1 * 10 - 11]; do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof(doublereal)); } e_wsfe(); } } natoms = molkst_1.numat; i__3 = *nmos; for (spqr_1.i = 1; spqr_1.i <= i__3; ++spqr_1.i) { i__4 = *nmos; for (j = 1; j <= i__4; ++j) { rjkaa[spqr_1.i + j * 10 - 11] *= .5; /* L180: */ } } if (first1) { spqr_1.i = i_indx(keywrd_1.keywrd, "MICROS", 80L, 6L); if (spqr_1.i != 0) { spqr_1.k = (integer) reada_(keywrd_1.keywrd, &spqr_1.i, 80L); lab = spqr_1.k; if (prnt) { s_wsfe(&io___59); e_wsfe(); } ntot = (integer) (xx + .5f); o__1.oerr = 0; o__1.ounit = 5; o__1.ofnm = 0; o__1.orl = 0; o__1.osta = "OLD"; o__1.oacc = 0; o__1.ofm = 0; o__1.oblnk = "ZERO"; f_open(&o__1); al__1.aerr = 0; al__1.aunit = 5; f_rew(&al__1); for (spqr_1.i = 1; spqr_1.i <= 3; ++spqr_1.i) { /* L190: */ s_rsfe(&io___61); do_fio(&c__1, line, 80L); e_rsfe(); } for (spqr_1.i = 1; spqr_1.i <= 1000; ++spqr_1.i) { s_rsfe(&io___63); do_fio(&c__1, line, 80L); e_rsfe(); /* L200: */ if (i_indx(line, "MICRO", 80L, 5L) != 0) { goto L210; } } L210: i__4 = lab; for (spqr_1.i = 1; spqr_1.i <= i__4; ++spqr_1.i) { s_rsfe(&io___64); do_fio(&c__1, line, 80L); e_rsfe(); /* Computing MAX */ /* Computing MIN */ i__2 = i_indx(line, "0", 80L, 1L), i__5 = i_indx(line, "1", 80L, 1L); i__3 = 0, i__1 = min(i__2,i__5) - 1; izero = max(i__3,i__1); i__3 = *nmos; for (j = 1; j <= i__3; ++j) { i__1 = j + izero - 1; if (s_cmp(line + i__1, "1", j + izero - i__1, 1L) != 0) { i__2 = j + izero - 1; s_copy(line + i__2, "0", j + izero - i__2, 1L); } i__1 = j + *nmos + izero - 1; if (s_cmp(line + i__1, "1", j + *nmos + izero - i__1, 1L) != 0) { i__2 = j + *nmos + izero - 1; s_copy(line + i__2, "0", j + *nmos + izero - i__2, 1L) ; } i__1 = j + izero - 1; microa[j + spqr_1.i * 10 - 11] = *(unsigned char *)&line[ i__1] - '0'; i__1 = j + *nmos + izero - 1; microb[j + spqr_1.i * 10 - 11] = *(unsigned char *)&line[ i__1] - '0'; /* L220: */ } if (prnt) { s_wsfe(&io___68); i__3 = *nmos; for (j = 1; j <= i__3; ++j) { do_fio(&c__1, (char *)µa[j + spqr_1.i * 10 - 11], (ftnlen)sizeof(integer)); } i__1 = *nmos; for (j = 1; j <= i__1; ++j) { do_fio(&c__1, (char *)µb[j + spqr_1.i * 10 - 11], (ftnlen)sizeof(integer)); } e_wsfe(); } spqr_1.k = 0; i__3 = *nmos; for (j = 1; j <= i__3; ++j) { /* L230: */ spqr_1.k = spqr_1.k + microa[j + spqr_1.i * 10 - 11] + microb[j + spqr_1.i * 10 - 11]; } if (spqr_1.k != ntot) { ntot = spqr_1.k; xx = (doublereal) spqr_1.k; s_wsfe(&io___69); do_fio(&c__1, (char *)&spqr_1.k, (ftnlen)sizeof(integer)); e_wsfe(); } /* L240: */ } first1 = FALSE_; goto L310; } nupp = kalpha; ndown = kbeta; ams = (nupp - ndown) * .5; if (prnt) { s_wsfe(&io___73); do_fio(&c__1, (char *)&ams, (ftnlen)sizeof(doublereal)); e_wsfe(); } if (nupp * ndown >= 0) { goto L270; } s_wsfe(&io___74); e_wsfe(); s_stop("", 0L); L270: lima = nfa[*nmos] / (nfa[nupp] * nfa[*nmos - nupp]); limb = nfa[*nmos] / (nfa[ndown] * nfa[*nmos - ndown]); lab = lima * limb; if (prnt) { s_wsfe(&io___77); do_fio(&c__1, (char *)&lab, (ftnlen)sizeof(integer)); e_wsfe(); } /* # IF(LAB.LT.101) GOTO 240 */ /* # WRITE(6,230) */ /* # 230 FORMAT(10X,24H TOO MANY CONFIGURATIONS/) */ /* # GOTO 160 */ /* # 240 CONTINUE */ perm_(nperma, &nupp, nmos, &c__10, &lima); perm_(npermb, &ndown, nmos, &c__10, &limb); spqr_1.k = 0; i__4 = lima; for (spqr_1.i = 1; spqr_1.i <= i__4; ++spqr_1.i) { i__3 = limb; for (j = 1; j <= i__3; ++j) { ++spqr_1.k; i__1 = *nmos; for (l = 1; l <= i__1; ++l) { microa[l + spqr_1.k * 10 - 11] = nperma[l + spqr_1.i * 10 - 11]; /* L290: */ microb[l + spqr_1.k * 10 - 11] = npermb[l + j * 10 - 11]; } } } /* L300: */ L310: lima = lab; limb = lab; } gse = 0.; i__1 = *nmos; for (spqr_1.i = 1; spqr_1.i <= i__1; ++spqr_1.i) { /* # IF(ABS(OCCA(I)-0.5).LT.0.01)GSE=GSE-0.25D0*XY(I,I,I,I) */ gse += eiga[spqr_1.i] * baseoc_1.occa[spqr_1.i - 1] * 2.; gse += xyijkl_1.xy[spqr_1.i + (spqr_1.i + (spqr_1.i + spqr_1.i * 10) * 10) * 10 - 1111] * baseoc_1.occa[spqr_1.i - 1] * baseoc_1.occa[spqr_1.i - 1]; i__3 = *nmos; for (j = spqr_1.i + 1; j <= i__3; ++j) { /* L320: */ gse += (xyijkl_1.xy[spqr_1.i + (spqr_1.i + (j + j * 10) * 10) * 10 - 1111] * 2. - xyijkl_1.xy[spqr_1.i + (j + (spqr_1.i + j * 10) * 10) * 10 - 1111]) * 2. * baseoc_1.occa[spqr_1.i - 1] * baseoc_1.occa[j - 1]; } } /* # IF(PRNT)WRITE(6,'('' GROUND STATE ENERGY:'',F13.6,'' E.V.'')')GSE */ /* .......... */ if (prnt) { s_wsfe(&io___81); e_wsfe(); } j = 0; i__3 = lab; for (spqr_1.i = 1; spqr_1.i <= i__3; ++spqr_1.i) { diag[spqr_1.i - 1] = diagi_(µa[spqr_1.i * 10 - 10], µb[ spqr_1.i * 10 - 10], &eiga[1], xyijkl_1.xy, nmos) - gse; /* L330: */ } L340: if (lab <= mdim) { goto L380; } x = -100.; i__3 = lab; for (spqr_1.i = 1; spqr_1.i <= i__3; ++spqr_1.i) { if (diag[spqr_1.i - 1] > x) { x = diag[spqr_1.i - 1]; j = spqr_1.i; } /* L350: */ } if (j != lab) { i__3 = lab; for (spqr_1.i = j; spqr_1.i <= i__3; ++spqr_1.i) { i1 = spqr_1.i + 1; i__1 = *nmos; for (spqr_1.k = 1; spqr_1.k <= i__1; ++spqr_1.k) { microa[spqr_1.k + spqr_1.i * 10 - 11] = microa[spqr_1.k + i1 * 10 - 11]; /* L360: */ microb[spqr_1.k + spqr_1.i * 10 - 11] = microb[spqr_1.k + i1 * 10 - 11]; } /* L370: */ diag[spqr_1.i - 1] = diag[i1 - 1]; } } --lab; goto L340; L380: /* BEFORE STARTING, CHECK THAT THE ROOT WANTED CAN EXIST */ if (lab < lroot) { s_wsfe(&io___83); e_wsfe(); s_wsfe(&io___84); e_wsfe(); s_stop("", 0L); } /* MAIN LOOP TO FILL SECULAR DETERMINANT */ ik = 0; i__3 = lab; for (spqr_1.i = 1; spqr_1.i <= i__3; ++spqr_1.i) { spqr_1.k = 0; i__1 = *nmos; for (j = 1; j <= i__1; ++j) { ioccb1[j - 1] = microb[j + spqr_1.i * 10 - 11]; iocca1[j - 1] = microa[j + spqr_1.i * 10 - 11]; /* L390: */ spqr_1.k += iocca1[j - 1]; } nalpha[spqr_1.i - 1] = spqr_1.k; if (prnt) { s_wsfe(&io___89); do_fio(&c__1, (char *)&spqr_1.i, (ftnlen)sizeof(integer)); i__1 = *nmos; for (spqr_1.k = 1; spqr_1.k <= i__1; ++spqr_1.k) { do_fio(&c__1, (char *)&iocca1[spqr_1.k - 1], (ftnlen)sizeof( integer)); } e_wsfe(); s_wsfe(&io___90); i__1 = *nmos; for (spqr_1.k = 1; spqr_1.k <= i__1; ++spqr_1.k) { do_fio(&c__1, (char *)&ioccb1[spqr_1.k - 1], (ftnlen)sizeof( integer)); } e_wsfe(); } spqr_1.is = 2; /* INNER LOOP TO FILL SECULAR DETERMINANT */ i__1 = spqr_1.i; for (spqr_1.k = 1; spqr_1.k <= i__1; ++spqr_1.k) { ++ik; cimat[ik - 1] = spqr_1.i * 1e-15 + j * 1e-18; ix = 0; iy = 0; i__4 = *nmos; for (j = 1; j <= i__4; ++j) { ioccb2[j - 1] = microb[j + spqr_1.k * 10 - 11]; iocca2[j - 1] = microa[j + spqr_1.k * 10 - 11]; ix += (i__2 = iocca1[j - 1] - iocca2[j - 1], abs(i__2)); /* L400: */ iy += (i__2 = ioccb1[j - 1] - ioccb2[j - 1], abs(i__2)); } /* CHECK IF MATRIX ELEMENT HAS TO BE ZERO */ if (ix + iy > 4 || nalpha[spqr_1.i - 1] != nalpha[spqr_1.k - 1]) { goto L420; } if (ix + iy == 4) { if (ix == 0) { cimat[ik - 1] = babbcd_(iocca1, ioccb1, iocca2, ioccb2, nmos); } else if (ix == 2) { cimat[ik - 1] = aabbcd_(iocca1, ioccb1, iocca2, ioccb2, nmos); } else { cimat[ik - 1] = aabacd_(iocca1, ioccb1, iocca2, ioccb2, nmos); } } else if (ix == 2) { cimat[ik - 1] = aababc_(iocca1, ioccb1, iocca2, ioccb2, nmos); } else if (iy == 2) { cimat[ik - 1] = babbbc_(iocca1, ioccb1, iocca2, ioccb2, nmos); } else { cimat[ik - 1] = diag[spqr_1.i - 1]; x = 0.; i__2 = *nmos; for (j = 1; j <= i__2; ++j) { /* L410: */ x += iocca1[j - 1] * ioccb1[j - 1]; } spin[spqr_1.i - 1] = -(xx - (nalpha[spqr_1.i - 1] << 1)) * ( xx - (nalpha[spqr_1.i - 1] << 1)) + x * 4.; } L420: ; } spqr_1.ispqr[spqr_1.i - 1] = spqr_1.is - 1; /* L430: */ } if (bigprt) { s_wsfe(&io___97); e_wsfe(); i__3 = -lab; vecprt_(cimat, &i__3); } else { if (prnt) { s_wsfe(&io___98); e_wsfe(); } if (prnt) { s_wsfe(&io___99); i__3 = lab; for (spqr_1.i = 1; spqr_1.i <= i__3; ++spqr_1.i) { do_fio(&c__1, (char *)&cimat[spqr_1.i * (spqr_1.i + 1) / 2 - 1], (ftnlen)sizeof(doublereal)); } e_wsfe(); } } rsp_(cimat, &lab, &lab, eig, conf); /* DECIDE WHICH ROOT TO EXTRACT */ kroot = 0; if (smult < -.1) { ret_val = eig[lroot - 1]; kroot = lroot; } if (bigprt) { s_wsfe(&io___103); e_wsfe(); matout_(conf, eig, &lab, &lab, &lab); } if (prnt) { s_wsfe(&io___104); e_wsfe(); } iroot = 0; for (spqr_1.i = 1; spqr_1.i <= 9; ++spqr_1.i) { /* L450: */ cimat[spqr_1.i - 1] = .1; } i__3 = lab; for (spqr_1.i = 1; spqr_1.i <= i__3; ++spqr_1.i) { x = xx * .5; ii = (spqr_1.i - 1) * lab; i__1 = lab; for (j = 1; j <= i__1; ++j) { ji = j + ii; x -= conf[ji - 1] * conf[ji - 1] * spin[j - 1] * .25; spqr_1.k = spqr_1.ispqr[j - 1]; if (spqr_1.k == 1) { goto L470; } i__2 = spqr_1.k; for (spqr_1.k = 2; spqr_1.k <= i__2; ++spqr_1.k) { li = spqr_1.ispqr[j + spqr_1.k * 100 - 101] + ii; /* L460: */ x += conf[ji - 1] * conf[li - 1] * 2.; } L470: /* L480: */ ; } y = (sqrt(x * 4. + 1.) - 1.) * .5; if ((d__1 = smult - x, abs(d__1)) < .01f) { ++iroot; if (iroot == lroot) { kroot = spqr_1.i; ret_val = eig[spqr_1.i - 1]; } } j = (integer) (y * 2. + 1.5); cimat[j - 1] += 1; /* L490: */ if (prnt) { s_wsfe(&io___110); do_fio(&c__1, (char *)&spqr_1.i, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&eig[spqr_1.i - 1], (ftnlen)sizeof( doublereal)); do_fio(&c__1, tspin + (j - 1 << 3), 8L); do_fio(&c__1, (char *)&x, (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&y, (ftnlen)sizeof(doublereal)); e_wsfe(); } } if (kroot == 0) { s_wsfe(&io___111); e_wsfe(); s_wsfe(&io___112); e_wsfe(); s_wsfe(&io___113); e_wsfe(); for (spqr_1.i = 1; spqr_1.i <= 7; ++spqr_1.i) { /* L500: */ if (cimat[spqr_1.i - 1] > .5) { s_wsfe(&io___114); do_fio(&c__1, tspin + (spqr_1.i - 1 << 3), 8L); i__3 = i_dnnt(&cimat[spqr_1.i - 1]); do_fio(&c__1, (char *)&i__3, (ftnlen)sizeof(integer)); e_wsfe(); } } s_stop("", 0L); } /* # M=0 */ /* # DO 440 I=1,NMOS */ /* # WRITE(6,*) */ /* # DO 440 J=1,NMOS */ /* # WRITE(6,*) */ /* # DO 440 K=1,NMOS */ /* # 440 WRITE(6,'(4I2,8F12.6)')I,J,K,M,(XY(I,J,K,L),L=1,NMOS) */ L520: if (force || last1) { /* REFORM DENSITY MATRIX */ spqr_1.k = (kroot - 1) * lab; i__3 = *nmos; for (spqr_1.i = 1; spqr_1.i <= i__3; ++spqr_1.i) { sum = 0.; i__1 = lab; for (j = 1; j <= i__1; ++j) { /* L530: */ /* Computing 2nd power */ d__1 = conf[j + spqr_1.k - 1]; sum += (microa[spqr_1.i + j * 10 - 11] + microb[spqr_1.i + j * 10 - 11]) * (d__1 * d__1); } /* L540: */ eiga[spqr_1.i] = sum - baseoc_1.occa[spqr_1.i - 1] * 2.; } l = 0; i__3 = molkst_1.norbs; for (spqr_1.i = 1; spqr_1.i <= i__3; ++spqr_1.i) { i__1 = spqr_1.i; for (j = 1; j <= i__1; ++j) { sum = 0.; i__2 = *nmos; for (spqr_1.k = 1; spqr_1.k <= i__2; ++spqr_1.k) { /* L550: */ sum += eiga[spqr_1.k] * coeffs[spqr_1.i + spqr_1.k * coeffs_dim1] * coeffs[j + spqr_1.k * coeffs_dim1]; } ++l; /* L560: */ densty_1.p[l - 1] += sum * one; } } } maxvec = 0; if (lspin) { maxvec = min(4,lab); } if (lspin && ne / 2 << 1 == ne) { s_wsfe(&io___117); e_wsfe(); } i__1 = *nmos; for (spqr_1.i = 1; spqr_1.i <= i__1; ++spqr_1.i) { i__3 = molkst_1.norbs; for (j = 1; j <= i__3; ++j) { /* L570: */ /* Computing 2nd power */ d__1 = coeffs[j + spqr_1.i * coeffs_dim1]; coeffs[j + spqr_1.i * coeffs_dim1] = d__1 * d__1; } } i__3 = maxvec; for (iuj = 1; iuj <= i__3; ++iuj) { iofset = (iuj - 1) * lab; s_wsfe(&io___120); do_fio(&c__1, (char *)&iuj, (ftnlen)sizeof(integer)); e_wsfe(); s_wsfe(&io___121); i__1 = lab; for (spqr_1.i = 1; spqr_1.i <= i__1; ++spqr_1.i) { do_fio(&c__1, (char *)&conf[spqr_1.i + iofset - 1], (ftnlen) sizeof(doublereal)); } e_wsfe(); i__1 = lab; for (spqr_1.i = 1; spqr_1.i <= i__1; ++spqr_1.i) { /* L580: */ /* Computing 2nd power */ d__1 = conf[spqr_1.i + iofset - 1]; conf[spqr_1.i - 1] = d__1 * d__1; } /* SECOND VECTOR! */ i__1 = *nmos; for (spqr_1.i = 1; spqr_1.i <= i__1; ++spqr_1.i) { sum = 0.; i__2 = lab; for (j = 1; j <= i__2; ++j) { /* L590: */ sum += (microa[spqr_1.i + j * 10 - 11] - microb[spqr_1.i + j * 10 - 11]) * conf[j - 1]; } /* L600: */ eiga[spqr_1.i] = sum; } s_wsfe(&io___122); do_fio(&c__1, (char *)&eig[iuj - 1], (ftnlen)sizeof(doublereal)); e_wsfe(); s_wsfe(&io___123); i__1 = *nmos; for (spqr_1.i = 1; spqr_1.i <= i__1; ++spqr_1.i) { do_fio(&c__1, (char *)&eiga[spqr_1.i], (ftnlen)sizeof(doublereal)) ; } e_wsfe(); s_wsle(&io___124); e_wsle(); s_wsle(&io___125); do_lio(&c__9, &c__1, " SPIN DENSITIES FROM EACH ATOMIC ORBITAL", 44L); e_wsle(); s_wsle(&io___126); do_lio(&c__9, &c__1, " S PX " " PY PZ TOTAL", 74L); e_wsle(); i__1 = natoms; for (spqr_1.i = 1; spqr_1.i <= i__1; ++spqr_1.i) { il = molkst_1.nfirst[spqr_1.i - 1]; iu = molkst_1.nlast[spqr_1.i - 1]; l = 0; summ = 0.; i__2 = iu; for (spqr_1.k = il; spqr_1.k <= i__2; ++spqr_1.k) { ++l; sum = 0.; i__4 = *nmos; for (j = 1; j <= i__4; ++j) { /* L610: */ sum += coeffs[spqr_1.k + j * coeffs_dim1] * eiga[j]; } summ += sum; /* L620: */ eigs[l] = sum; } if (l == 4) { s_wsfe(&io___130); do_fio(&c__1, (char *)&spqr_1.i, (ftnlen)sizeof(integer)); i__2 = l; for (spqr_1.k = 1; spqr_1.k <= i__2; ++spqr_1.k) { do_fio(&c__1, (char *)&eigs[spqr_1.k], (ftnlen)sizeof( doublereal)); } do_fio(&c__1, (char *)&summ, (ftnlen)sizeof(doublereal)); e_wsfe(); } else { s_wsfe(&io___131); do_fio(&c__1, (char *)&spqr_1.i, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&eigs[1], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&summ, (ftnlen)sizeof(doublereal)); e_wsfe(); } /* L630: */ } /* L640: */ } return ret_val; } /* meci_ */