NOTE THIS LIBRARY OF PROGRAMS FOR THE RECURSION METHOD WAS WRITTEN BY CHRIS NEX AND IS MADE AVAILABLE AT NO COST. CHRIS AND I ONLY ASK THAT USERS ACKNOWLEDGE THE RECURSION LIBRARY WHEN IT HAS MADE A SIGNIFICANT CONTRIBUTION TO THEIR WORK, AND THAT THEY NOTIFY ME OF ANY ERRORS THEY DISCOVER. I STRONGLY ENCOURAGE USERS TO CONSULT WITH ME BEFORE EMBARKING ON LARGE PROJECTS TO AVOID THE PITFALLS OTHERS HAVE ENCOUNTERED AND TO AVOID DUPLICATION OF EFFORT. PLEASE ENJOY THE OPPORTUNITIES THESE PROGRAMS OFFER. ROGER HAYDOCK EUGENE, OREGON 2000 APRIL 5 THE CAMBRIDGE RECURSION LIBRARY (1986) ...... CHRIS NEX CONTENTS 1. A BRIEF SURVEY AND DESCRIPTION OF THE CAMBRIDGE RECURSION LIBRARY (1986 VERSION WHICH INCLUDES SOME RECENT ROUTINES WITHOUT EXAMPLES OF USE) 2. A LIST OF RELEVANT SUBROUTINES FROM THE T.C.M. LIBRARY 3. SPECIFICATIONS OF THOSE SUBROUTINES 4. LISTINGS OF THE LIBRARY SUBROUTINES 5. AN EXAMPLE RUN OF THE BASIC RECURSION PROGRAM (EXRECAL). THE CLUSTER OF ATOMS USED IS SMALLER THAN IS GENERALLY USEFUL IN PRODUCTION RUNS , BUT THE PROGRAM INCLUDES A SET OF TYPICAL BACK-UP SUBROUTINES DEFINING THE HAMILTONIAN. 6. AN EXAMPLE RUN OF A PROGRAM PROCESSING THE RESULTS OF 5. TO OBTAIN A SELECTION OF FREQUENTLY REQUIRED QUANTITIES.(EXPROC) 7. AN EXAMPLE RUN OF A PROGRAM GENERATING THE RECURSION COEFFICIENTS FOR THE PHONON LOCAL DENSITY OF STATES FOR AN AMORPHOUS CLUSTER. (EXPHNDS) 8. AN EXAMPLE RUN OF A PROGRAM TO COMPUTE THE PHONON LDOS FROM THE RECURSION COEFFICIENTS TABULATING IT AS A FUNCTION OF FREQUENCY. (EXPHNPLT) 9. AN EXAMPLE OF A PROGRAM ADDING ATOMS TO THE SURFACE OF A CLUSTER AND PERFORMING 'ORBITAL PEELING' AS DEFINED BY N.R.BURKE. (EXORPEEL) 10. AN EXAMPLE RUN OF A PROGRAM TO PROCESS THE RESULTS OF 9. AND TABULATE TOTAL ENERGY DIFFERENCES. (EXPEELED) THE PDS MEMBERS OF SSTLIB.CAMREC86 CONTAINING THE EXAMPLE PROGRAMS ARE PREFIXED EX-, AND THE CORRESPONDING OUTPUT OP- WHILE THE APPROPRIATE RECURSION COEFFICIENTS FOR RE-INPUT ARE PREFIXED CFT- . 1: THE CAMBRIDGE RECURSION LIBRARY 1986 1.0:PREFACE THE RECURSION SUBROUTINE LIBRARY CONTAINS ROUTINES TO ENABLE THE COMPUTATION OF THE INTEGRATED DENSITY OF STATES AND INTEGRALS OF AN ARBITRARY FUNCTION TIMES THE DENSITY OF STATES TO BE EFFECTED. THE APPROPRIATE HAMILTONIANS ARE THOSE FOR A FINITE (60 - 6000) CLUSTER OF ATOMS WITH N ORBITALS ( 1 - 10 OR MORE) AND 'NEIGHBOUR' INTERACTIONS , GIVING AN EFFECTIVE LARGE , SPARSE MATRIX EIGEN- PROBLEM FOR WHICH IS REQUIRED THE DISTRIBUTION OF EIGENVALUES CORRESPONDING TO THE LOCAL DENSITY OF STATES. THE MATHEMATICS OF THE RECURSION METHOD IS DESCRIBED IN HEINE, HAYDOCK, KELLY & BULLETT (1) AND OF THE COMPUTATION OF INTEGRALS IN NEX (2). FOR A DESCRIPTION OF THE BASIC MATHEMATICAL ROUTINES SEE NEX (3), ALTHOUGH THERE ARE A FEW ADDITIONAL SPECIALISED ROUTINES INCLUDED IN THIS COLLECTION. TWO INDEPENDANT PARAMETERS AT THE DISPOSAL OF THE USER, GIVING A CONTROL ON THE ACCURACY , ARE THE SIZE OF THE CLUSTER AND THE 'NUMBER OF LEVELS' IN THE RECURSION. THE APPROPRIATE VALUES OF THESE PARAMETERS WILL VARY FROM PROBLEM TO PROBLEM DEPENDING ON THE HAMILTONIAN, THE NATURE OF THE LATTICE, AND OF THE PHYSICAL QUANTITY OF INTEREST ULTIMATELY BEING COMPUTED, AND SHOULD BE DETERMINED BY SOME TRIAL RUNS OF THE PROGRAM. FOR FIVE D-ELECTRONS PER SITE IN FCC AND BCC CLUSTERS AND THE CALCULATION OF STRUCTURAL ENERGIES IT HAS BEEN FOUND THAT THE NUMBER OF LEVELS APPROPRIATE IS 2-5 TIMES THE 'RADIUS' OF THE CLUSTER (IN NUMBERS OF ATOMS) IN CUBOID CLUSTERS . A MATHEMATICALLY EXACT ERROR BOUND ON THE INTEGRATED DENSITY OF STATES IS AVAILABLE IF ONLY EXACT COEFFICIENTS ARE USED , BUT , AS IS THE NATURE OF SUCH QUANTITIES, THIS BOUND IS PESSIMISTIC , USUALLY BY AT LEAST AN ORDER OF MAGNITUDE. 1.1:INTRODUCTION THE ROUTINES IN THE LIBRARY FALL INTO TWO CATEGORIES : THOSE WHICH ASSIST IN THE SETTING UP OF THE PROBLEM, AND THOSE WHICH COMPUTE FROM THE CALCULATED COEFFICIENTS THE LOCAL DENSITY OF STATES AND RELATED FUNCTIONS. THESE ARE DESCRIBED SEPARATELY IN THE FOLLOWING TWO SECTIONS. THERE EXISTS IN THE LITERATURE SOME SLIGHT CONFUSION OF NOTATION, BUT THAT USED IN THE LIBRARY IS CONSISTENT WITHIN ITSELF AND WITH THE RECURSION FORMULA : B(N+1).PSI(N+1) = ( H - A(N) ).PSI(N) - B(N).PSI(N-1) ,N=0,1,.. WHERE A(N) = AND B(N+1) IS CHOSEN SUCH THAT =1. PSI(-1) (=0) AND PSI(0) ARE GIVEN. (WITHIN THE FORTRAN THE SUBSCRIPTS ARE INCREASED BY 1 TO CONFORM TO FORTRAN IV CONVENTION, ALTHOUGH ALL ROUTINES HAVE ALSO RUN IN FORTRAN 77 UNDER THE IBM VS COMPILER) 1.2: PROBLEM INITIALISATION ROUTINES THE SYSTEM ADOPTED FOR THE SETTING UP OF THE PROBLEM FOR THE APPLICATION OF THE RECURSION METHOD CONSISTS OF A LIBRARY OF ROUTINES DIVIDED INTO THREE SECTIONS , WITH ONLY A SUBSET BEING USED IN ANY ONE PROBLEM , AND FOR SOME PROBLEMS THE USER SUPPLYING HIS OWN ROUTINES TO REPLACE AND SUPPLEMENT SOME OF THE LIBRARY ONES. THE THREE GROUPS OF ROUTINES ARE AS FOLLOWS: THE FIRST SERVE TO DEFINE A CLUSTER OF 'ATOM SITES' , THE NEXT SET UP A 'NEIGHBOUR MAP ' ( OR MATRIX SPARSITY PATTERN ) , AND THE LAST SET UP THE COMPONENT SUBMATRICES OF THE HAMILTONIAN . THE APPROPRIATE RESULTS OF THESE ROUTINES ARE LOADED BY THE USER INTO A COMMON BLOCK WHICH SHOULD THEN CONTAIN ALL THE INFORMATION NECESSARY TO DEFINE THE MATRIX REPRESENTATION OF THE HAMILTONIAN.THIS IS THEN USED BY THE SUBROUTINE AT THE HEART OF THE RECURSION LIBRARY ( RECAL ) TO OPERATE WITH THE HAMILTONIAN AND PRODUCE THE RECURSION COEFFICIENTS. 1.2.1 SUBROUTINES FOR GENERATING A CLUSTER : BCCLAT AND FCCLAT GENERATE CUBOID CLUSTERS IN BCC AND FCC LATTICES RESPECTIVELY . MORE EXOTIC CLUSTER SHAPES MAY BE OBTAINED BY USE OF THE FOLLOWING ROUTINES BUT THEIR USE NEEDS JUSTIFICATION . ONION WILL CLASSIFY EACH ATOM ACCORDING TO THE NUMBER OF 'HOPS' FROM A GIVEN GROUP OF 'CENTRAL' ATOMS , BUT CAN ONLY BE USED AFTER THE ATOM 'NEIGHBOURS' HAVE BEEN DEFINED (E.G. BY NNCAL - SEE BELOW ) ; PEEL EXCISES ATOMS FROM A CLUSTER AND RENUMBERS THE REMAINING ONES . THE LAST TWO ROUTINES MAY BE USEFUL TO USERS WITH PARTICULAR CLUSTER-SHAPE REQUIREMENTS , BUT IN GENERAL THE CUBOID CLUSTERS HAVE PROVED SATISFACTORY FOR ALMOST ALL APPLICATIONS TO DATE. 1.2.2 THE SETTING-UP OF THE HAMILTONIAN MATRIX. NNCAL IS THE CENTRAL SUBROUTINE USED TO SET UP A 'NEIGHBOUR MAP' OF THE ATOM CLUSTER. THIS USES ANOTHER ROUTINE ('IBONDS') DEFINING THE NEIGHBOUR CRITERION ; THIS MAY BE REPLACED BY THE USER WITH MORE PARTICULAR REQUIREMENTS ( I.E. FIRST AND SECOND NEAREST NEIGHBOURS) OR A DIFFERENT LATTICE . ONCE THE NEIGHBOUR MAP IS SET UP, THE NON-EQUIVALENT INTERACTIONS MAY BE SET UP USING MMCAL. THE 'EQUIVALENCE' OF TWO INTERACTIONS IS DEFINED IN A ROUTINE 'EQUIV' WHICH AGAIN MAY BE REPLACED BY A USER WITH SPECIFIC REQUIREMENTS . MMCAL GENERATES A LIST OF NON-EQUIVALENT VECTORS WHICH MAY THEN BE USED TO SET UP THE 'INTERACTION BLOCKS' OF THE HAMILTONIAN MATRIX, FOR INSTANCE AS IN SLKODE TO GIVE THE SLATER- KOSTER INTERACTION MATRICES FOR THE 5 D-ORBITALS . SUBROUTINES GENERATING THE THE S,P OR D SLATER-KOSTER MATRICES ARE SUPPLIED IN THE LIBRARY AND HAVE NAMES STARTING SK-- . THE ESSENTIAL INFORMATION DEFINING THE HAMILTONIAN MATRIX FOR A REGULAR LATTICE MAY BE SET UP BY THE ROUTINE SETUP AND AMENDMENTS TO THE MATRIX MADE USING ADDAT, ORPEEL OR PEEL (ADDING OR EXCISING ATOMS FROM THE CLUSTER ) . FOR A RANDOM SYSTEM AN ALTERNATIVE FORM OF STORAGE OF THE HAMILTONIAN MATRIX IS PROBABLY MORE EFFICIENT AND ROUTINES SCAN AND SCAN1 ARE SUPPLIED TO FACILITATE THE SETTING UP AND OPERATION OF THE MATRIX OF SUCH A PROBLEM. 1.2.3 GENERATION OF RECURSION COEFFICIENTS. THE APPROPRIATE RESULTS OF THE ABOVE SUBROUTINES SHOULD THEN BE LOADED BY THE USER INTO A COMMON BLOCK AND USED BY A ROUTINE , SUCH AS HOP OF THE EXAMPLE RUN , TO SUPPLY THE HAMILTONIAN OPERATOR REQUIRED BY THE ACTUAL RECURSION SUBROUTINE, RECAL. 1.3: PROCESSING THE RECURSION COEFFICIENTS. TWO METHODS ARE PROVIDED FOR THE CALCULATION OF THE LOCAL DENSITY OF STATES : TERMINOR AND QUADRATURE, AS DESCRIBED IN NEX(3). IF INTEGRALS ARE THE FINAL OBJECT OF THE CALCULATION THEN QUADRATURE (DENQD) IS APPROPRIATE , WHILE TO ESTIMATE THE DENSITY ITSELF WITH SECURE KNOWLEDGE OF THE BAND-GAPS THE ANALYTIC TERMINATOR(DENSQ OR DENCRS) CAN BE USED. IF THE SUM OF SEVERAL DENSITIES OF STATES IS REQUIRED (E.G. OVER SEVERAL ORBITALS), THIS AGAIN MAY BE DONE EFFICIENTLY BY 'SUMMING' THEIR TRIDIAGONALISATIONS USING RECSUM TO PRODUCE A RESULTANT SET OF COEFFICIENTS . IF THE FUNCTIONS SUMMED ARE VERY DIFFERENT IN CHARACTER, THIS MAY INTRODUCE SIGNIFICANT ERROR, IN WHICH CASE THE FUNCTIONS SHOULD BE TABULATED SEPARATELY , BUT AGAIN FOR GRAPHICAL PURPOSES THE RESULTS OF RECSUM ARE USUALLY ADEQUATE. FOR QUANTITATIVE WORK , ONLY THE INTEGRATED DENSITY OF STATES , N(E), FUNCTION VALUES SHOULD BE RELIED ON, AND ANY OTHER FUNCTIONS COMPUTED FROM THIS ONE. IT SHOULD BE NOTED THAT THE APPROXIMATIONS OBTAINED FROM THE QUADRATURE FORMULA ARE NOT ANALYTICALLY RELATED IN THE WAY ONE MIGHT EXPECT: INTEGRAL E N(E) DE DOES NOT EQUAL E N(E) - INTEGRAL N(E) DE WHERE THE L.H.S AND N(E) REPRESENT VALUES OBTAINED DIRECTLY FROM DENQD. SUCH IDENTITIES ARE USUALLY SATISFIED APPROXIMATELY, BUT IF PRECISE ANALYTICITY IS DEMANDED THE APPROXIMATION TO N(E) SHOULD BE TAKEN AND ALL OTHER RESULTS COMPUTED DIRECTLY FROM IT . THE EXCEPTION IS THAT D/DE N(E) = N(E) BY DEFINITION . REFERENCES : REFERENCES 1. HEINE,HAYDOCK,KELLY & BULLETT SOLID STATE PHYSICS VOL.35 (1979) 2. NEX C.M.M. J.PHYS. A VOL.11 653 ET SEQ. (1978) 3. NEX C.M.M. COMP.PHYS.COMM (1985) 2: A LIST OF RELEVANT ROUTINES FROM THE TCM LIBRARY THE DOCUMENTATION MAY BE OBTAINED BY COMMANDING THE JOB IN SSTLIB.CAMREC86:PRINT AND A MAGNETIC TAPE LOADED BY COMMANDING THE MEMBER TAPE. THE LIBRARY ROUTINES ARE FILED AS MEMBERS OF THE P.D.S.S SECTION 1 : SSTLIB.TCMLIB.FOR (FORTRAN)& PRECOMPILED SSTLIB.TCMLIB.COM (SINGLE PRECISION), SSTLIB.TCMLIB.DCOM ( DOUBLE PRECISION) SECTION 2: SSTLIB.RECSET (FORTRAN ) & SSTLIB.RECLIB (SINGLE PRECISION ONLY PRECOMPILED) INDEX SECTION 1 : MATHEMATICAL SUBROUTINES FOR PROCESSING COEFFICIENTS BNDCRD ESTIMATES BAND EDGES FROM TABULATED DATA BNDEST CRUDE BAND EDGE ESTIMATION FROM CONTINUED FRACTION BNDREF REFINEMENT OF BAND EDGE ESTIMATION FROM CONTINUED FRACTION BNDWT ESTIMATES THE WEIGHTS IN DIFFERENT BANDS IN A CONTINUED FRACTION REPRESENTATION OF A DENSITY OF STATES CFGEN FINDS TRIDIAGONALISATION CORRESPONDING TO WEIGHT FUNCTION CFGPFN FINDS TRIDIAGONALISATION CORESPONDING TO THE UNION OF WEIGHT FUNCTIONS CRECAL RECURSION CALCULATION FOR COMPLEX HERMITIAN MATRIX DENCRS EVALUATES THE WEIGHT FUNCTION WITH A GENERAL TERMINATOR DENCRQ EVALUATES THE COMPLEX GREENIAN WITH GENERAL TERMINATOR DENINT ESTIMATES INTEGRATED WEIGHT FUNCTION FROM TRIDIAGONALISATION(S) DENQD ESTIMATES WEIGHT FUNCTION ETC. FROM TRIDIAGONALISATION USING QUADRATURE METHOD DENSQ EVALUATES WEIGHT FUNCTION USING A SQUARE-ROOT TERMINATOR FENVAL ESTIMATES 'FERMI ENERGY' GIVEN TRIDIAGONALISATION(S) EDIFF EVALUATES TOTAL ENERGY DIFFERENCE IN 'ORBITAL PEELING' NUMC COUNTS NO. OF EIGENVALUES OF 'JACOBI' MATRIX >X NUMD COMPUTES DET/DET' OF 'JACOBI' MATRIX PLYVAL COMPUTES ORTHOGONAL POLYNOMIALS RECAL FINDS THE TRIDIAGONALISATION OF A SPARSE SYMMETRIC MATRIX RECNO FINDS THE APPROXIMATE TRIDIAGONALISATION OF THE GENERALISED EIGENPROBLEM (M - E S) X = 0 RECPER COMPUTES A PERTURBATION OF A MATRIX TRIDIAGONALISATION RECQD COMPUTES THE COEFFICIENTS IN A GAUSSIAN QUADRATURE RECRTS FINDS THE EIGENVALUES OF A JACOBI MATRIX RECSUM FINDS THE CONTINUED FRACTION COEFFICIENTS OF A SUM OF CONTINUED FRACTIONS RECWT CALCULATES GAUSSIAN QUADRATURE WEIGHT AT A HYPOTHETICAL NODE TABAN IDENTIFIES EXTREMAL VALUES OF A TABULATED FUNCTION TERMGN GENERATES A MATCHING TERMINATOR FOR A CONTINUED FRACTION WTMIN FINDS MINIMUM IN 'LOCAL' WEIGHT OF A CONTINUED FRACTION WEIGHT FUNCTION SECTION 2 : SUBROUTINES FOR SETTING UP THE INITIAL PROBLEM ADDAT EXTENDS THE HAMILTONIAN MATRIX OF A GIVEN CLUSTER BCCBFE 'IBONDS' SUITABLE FOR A BCC LATTICE GENERATED BY BCCLAT BCCLAT GENERATES A BCC LATTICE OF APPROXIMATELY THE DESIRED SIZE EQUIV FUNCTION FOR USE IN MMCAL DEFINING EQUIVALENCE OF VECTORS FCCBND 'IBONDS' SUITABLE FOR A FCC LATTICE GENERATED BY FCCLAT FCCLAT GENERATES A FCC LATTICE OF APPROXIMATELY THE DESIRED SIZE MMCAL GENERATES THE NON-EQUIVALENT INTERACTION VECTORS FOR A LATTICE NNCAL GENERATES THE 'NEIGHBOUR' MAP FOR A GIVEN LATTICE ONION 'SPHERICISES' A GIVEN CLUSTER ORPEEL DELETES A ROW & COLUMN OF THE HAMILTONIAN MATRIX OR RESTORES IT OUT OUTPUTS TO LINE-PRINTER THE HAMILTONIAN MATRIX INFORMATION PEEL EXCISES THE OUTER LAYERS OF A CLUSTER AFTER PROCESSING BY ONION SCAN GENERATES NEIGHBOURS SYSTEMATICALLY FROM A NEIGHBOUR MAP SCAN1 IDENTICAL IN ALL BUT NAME TO SCAN SELFD GENERATES A DIAGONAL 'SELF INTERACTION' MATRIX FOR D-ELECTRONS SELFP GENERATES A DIAGONAL 'SELF INTERACTION' MATRIX FOR P-ELECTRONS SELFS GENERATES A 'SELF INTERACTION' FOR AN S-ELECTRON SETUP GENERATES THE HAMILTONIAN MATRIX FOR A REGULAR CLUSTER SKDD COMPUTES THE SLATER-KOSTER DD OVERLAP MATRIX SKPD COMPUTES THE SLATER-KOSTER PD OVERLAP MATRIX SKPP COMPUTES THE SLATER-KOSTER PP OVERLAP MATRIX SKSD COMPUTES THE SLATER-KOSTER SD OVERLAP MATRIX SKSP COMPUTES THE SLATER-KOSTER SP OVERLAP MATRIX SKSS COMPUTES THE SLATER-KOSTER SS OVERLAP MATRIX SLKODE COMPUTES THE SLATER-KOSTER OVERLAP MATRICES FOR D-ELECTRONS ./ ADD NAME=BNDCRD SUBROUTINE BNDCRD(ET,IC,WTT,NET,NB,AA,RNG) DIMENSION ET(NET),IC(NET),WTT(NET),AA(NB),RNG(NB) ESTIMATES BAND EDGE POSITIONS FROM TABULAR INFORMATION CODED AS IN TABAN. THIS CONSISTS OF A LIST OF POINTS WHERE THE FUNCTION VALUE EITHER CROSSES A LOWER OR UPPER THRESHOLD OR HAS A LOCAL EXTREMUM BETWEEN THEM AND THE ESTIMATION IS A STRAIGHTFORWARD ANALYSIS OF THE TABULAR DATA AS GENERATED BY TABAN. ARGUMENTS : (* INDICATES AN OVERWRITTEN ARGUMENT) ET LIST OF ORDINATES OF EXTREMA (AS O/P FROM TABAN) IC* CODE CHARACTERISING THE EXTREMA (AS O/P FROM TABAN) ON INPUT: A + SIGN INDICATES INCREASING FUNCTION VALUES AND A - SIGN DECREASING ONES. AT LOCAL EXTREMA THIS EXTENDS TO + INDICATING A MINIMUM AND - A MAXIMUM. THE ABSOLUTE VALUE IS CODED BELOW. 1 FUNCTION VALUE CROSSES LOWER THRESHOLD 2 FUNCTION VALUE CROSSES UPPER THRESHOLD 3 FUNCTION VALUE CROSSES BOTH THRESHOLDS 4 FUNCTION VALUE IS A LOCAL EXTREMUM ON OUTPUT: (ONLY IF -4 ON INPUT) 0 CORRESPONDING LOCAL MINIMUM HAS BEEN USED IN THE CONSTRUCTION OF THE BAND POSITIONS WTT IF IABS(IC(I))=4 THIS GIVES THE LOCAL EXTREMUM VALUE NET NUMBER OF VALUES TABULATED IN ET,IC,WTT NB* MAXIMUM NUMBER OF BANDS EXPECTED ON OUTPUT CONTAINS THE NUMBER OF BANDS IDENTIFIED OR IF 0 THEN A MISMATCH OF LEFT AND RIGHT EXTREMA OR NEGATIVE THEN -NB BANDS HAVE BEEN IDENTIFIED, AND THAT NUMBER IS GREATER THAN THE INPUT NB AA* LIST OF LEFT BAND EXTREMA RNG* LIST OF BAND WIDTHS ./ ADD NAME=BNDEST SUBROUTINE BNDEST(A,B2,LL,EPS,AA,RNG,NB,EV,FEV,IC,NET,DE,WK,NW) DIMENSION A(LL),B2(LL),AA(NB),RNG(NB),EV(NET),FEV(NET),IC(NET) 1, WK(NW,4),P(2,3) ESTIMATES THE BAND EDGE POSITIONS OF A CONTINUED FRACTION REPRESENTATION OF A DENSITY OF STATES, GIVEN THE MAXIMUM NUMBER OF BANDS PRESENT. THIS IS DONE BY TABULATION OF THE LOCAL QUADRATURE WEIGHT FOLLOWED BY AN ANALYSIS OF LOCAL MINIMA AND THRESHOLDS IN THAT TABLE. THE LOWER THRESHOLD IS THE VALUE OF THE LOCAL WEIGHT AT THE EDGE OF A SIMPLE ELLIPTICAL BAND. THE ACCURACY OF THE BAND EDGE LOCATION IS DETERMINED BY THE NUMBER, NW, OF POINTS IN THE TABULATION. THIS ROUTINE CALLS TABAN,BNDCRD,RECWT,WTMIN ARGUMENTS (* INDICATES AN OVERWRITTEN ARGUMENT) A DENOMINATOR CONTINUED FRACTION COEFFICIENTS ; I=1,LL-1 B2 NUMERATOR CONTINUED FRACTION COEFFICIENTS ; I=1,LL LL LENGTH OF CONTINUED FRACTION EPS THRESHOLD ACCURACY IN COMPUTATIONS AA* LIST OF ESTIMATED BAND LEFT EXTREMA RNG* LIST OF ESTIMATED BAND WIDTHS NB* ON INPUT: MAXIMUM NUMBER OF BANDS ON OUTPUT:IF >0 NUMBER OF BANDS INDENTIFIED IF OR= 2*NB ON OUTPUT: IF >0 NUMBER OF MATCHING POINTS IF = SUM I=1 TO N { F(X(I)) * G( X(I)) * W(I) } . THEN, WITH P(X,0) = 1 AND P(X,-1) = 0 , FOR N = 1,2...,LM1 A(N) = / < P(X,N-1) , P(X,N-1) > B2(N) = < P(X,N-1) , P(X,N-1) > / < P(X,N-2) , P(X,N-2) > P(X,N) = ( X-A(N) ) * P(X,N-1) - B2(N) * P(X,N-2) THE POLYNOMIALS P(X,N) ARE REPRESENTED BY THEIR VALUES AT THE POINTS X(I) , WITH AN ARBITRARY NORMALISATION TO PRESERVE NUMERICAL STABILITY. ÕHAYDOCK & NEXå ARGUMENTS : (* INDICATES AN OVERWRITTEN ARGUMENT) X NODES OF THE SUMMATION INNER PRODUCT W WEIGHTS OF THE SUMMATION INNER PRODUCT N NUMBER OF TERMS IN THE SUMMATION INNER PRODUCT EPS A VALUE OF B2(I) LESS THAN EPS TERMINATES ROUTINE A* DIAGONAL COEFFICIENTS IN THE RECURRENCE A(I),I=1,LM1 B2* OFF-DIAGONAL COEFFICIENTS IN THE RECURRENCE B2(I),I=1,LM1 LM1* ON INPUT REQUIRED LENGTH OF RECURRENCE ON OUTPUT LENGTH OF RECURRENCE ACTUALLY GENERATED WK WORK ARRAY OF LENGTH AT LEAST 2*N ./ ADD NAME=CFGPGN SUBROUTINE CFGPGN(AA,RNG,WB,NBP1,IC,EPS,A,B2,LM1,WK,NW) DIMENSION AA(NBP1),RNG(NBP1),WB(NBP1),A(LM1),B2(LM1) DIMENSION WK(NW,2,NBP1) GENERATES COEFFICIENTS OF THE J-MATRIX OR 3-TERM RECURRENCE CORRESPONDING TO A LINEAR COMBINATION OF THE SAME TYPE OF WEIGHT FUNCTION. I.E. TO A SET OF BANDS IN AN DENSITY OF STATES. THE ALGORITHM IS THE SAME AS THAT USED IN CFGEN, BUT THE SUMMATION IN THE INNER PRODUCT NOW RUNS IN ADDITION OVER THE SEPARATE BANDS: = SUM J=1 TO NB WB(J)*ÕSUM I=1 TO N {F(X(I,J))*G(X(I,J))*W(I)} å WHERE X(I,J) = X(I) * RNG(J) + AA(J) THE X(I) AND W(I) MAY BE SET EXPLICITLY BY THE USER, OR THERE ARE PROVIDED IN THE ROUTINE MEANS OF GENERATING THEM FOR WEIGHT FUNCTIONS OF THE FORM SQRT(1-X*X) {CHEBYSHEV POLYNOMIALS OF THE SECOND KIND} OR A CONSTANT { LEGENDRE POLYNOMIALS}. ÕHAYDOCK & NEXå THIS ROUTINE USES RECQD,RECRTS,NUMC,NUMD IF IC=2 ARGUMENTS: (* INDICATES AN OVERWRITTEN ARGUMENT) AA LIST OF LEFT EXTREMA OF THE BANDS RNG LIST OF THE WIDTHS OF EACH BAND WB LIST OF WEIGHTS OF EACH BAND NBP1 1+ NUMBER OF BANDS IC* CODE INDICATING TYPE OF BANDS : 1 SQRT(1.0-E*E) -1 B(N+1).PSI(N+1) = (H-A(N)).PSI(N) - B(N).PSI(N-1) = 1 THIS ROUTINE MUST BE SUPPORTED BY A SUBROUTINE SEGMENT HOP, TO CARRY OUT STEPS OF THIS PROCESS : SUBROUTINE HOP(X,Y,A) COMPLEX X,Y DIMENSION X( ),Y( ) WHICH, WITH INPUT X=X, Y=Y, WILL PRODUCE Y = H.X+Y AND A = X.H.X = ARGUMENTS : (* DENOTES OVERWRITTEN ARGUMENT) HOP NAME OF SUBROUTINE SEGMENT PSI* COMPLEX ARRAY INPUT B(0)*PSI(0) OUTPUT B(LL)*PSI(LL) PMN* COMPLEX ARRAY INPUT -B(0)*PSI(-1) OUTPUT -B(LL)*PSI(LL-1) WHERE THE NORMS OF ALL VECTORS, PSI(I), ARE EQUAL TO UNITY M DIMENSION OF MATRIX A* OUTPUT A(I) = A(I-1), I=1,LL-1 B2* INPUT B2(1) = B(0)**2 OUTPUT B2(I) = B(I-1)**2, I=1,LL LL 1+LENGTH OF TRIDIAGONALISATION REQUIRED >OR= 2 ./ ADD NAME=DENCRQ COMPLEX FUNCTION DENCRQ(E,A,B2,LL,AA,RNG,WB,NB,AM,BM2) DIMENSION A(LL),B2(LL),AA(NB),RNG(NB),WB(NB),AM(LL),BM2(LL) COMPUTES THE VALUE OF A GREENIAN REPRESENTED BY A CONTINUED FRACTION USING A TERMINATOR BASED ON THE NUMBER, WEIGHTS AND POSITIONS OF SEPARATE BANDS USING A GENERAL PRESCRIPTION ÕHAYDOCK AND NEX - TO APPEARå. THE MATCHING CONTINUED FRACTION WITH SQUARE ROOT BAND EDGES MAY BE GENERATED USING CFGPGN OR TERGEN, AND SHOULD BE OF THE SAME LENGTH AS THE ORIGINAL. THE FUNCTION F(E) = SUM 8.0*WB(K)/RNG(K)**2 * Õ E - { AA(K) + RNG(K)*0.5 } - SQRT{E-AA(K)}*SQRT{AA(K)+RNG(K)-E} å IS ASSUMED TO CORRESPOND TO THE SUPPLIED COEFFICIENTS AM(I), BM2(I). T(E) = {S(E,N-1)-F(E)*R(E,N)}/{ S(E,N-2)-F(E)*R(E,N-1)} /BM2(N) N(E) = -1/PI * IM { G(E) } G(E) = ÕQ(E,N-1)-B2(N)*T(E)*Q(E,N-2)å/ÕP(E,N)-B2(N)*T(E)*P(E,N-1)å WHERE N=LL AND P,Q AND R,S ARE THE ORTHOGONAL POLYNOMIALS OF THE FIRST AND SECOND KINDS CORRESPONDING TO A,B2 AND AM,BM2 RESPECTIVELY. THIS ROUTINE USES PLYVAL ARGUMENTS : (* INDICATES AN OVERWRITTEN ARGUMENT) DENCRQ TAKES THE REQUIRED VALUE E ARGUMENT OF CONTINUED FRACTION A DENOMINATOR COEFFICIENTS OF CONTINUED FRACTION I=1,LL-1 B2 NUMERATOR COEFFICIENTS OF CONTINUED FRACTION I=1,LL LL LENGTH OF CONTINUED FRACTION AA LIST OF BAND LEFT EXTREMA RNG LIST OF BAND WIDTHS WB LIST OF WEIGHTS OF BANDS NB NUMBER OF BANDS (GREATER THAN 0) AM LL-1 DENOMINATOR COEFFICIENTS OF MATCHING CONTINUED FRACTION BM2 LL NUMERATOR COEFFICIENTS OF MATCHING CONTINUED FRACTION ./ ADD NAME=DENCRS FUNCTION DENCRS(E,A,B2,LL,AA,RNG,WB,NB,AM,BM2) DIMENSION A(LL),B2(LL),AA(NB),RNG(NB),WB(NB),AM(LL),BM2(LL) COMPUTES THE VALUE OF A CONTINUED FRACTION USING A TERMINATOR BASED ON THE NUMBER, WEIGHTS AND POSITIONS OF SEPARATE BANDS USING A GENERAL PRESCRIPTION ÕHAYDOCK AND NEX - TO APPEARå. THE MATCHING CONTINUED FRACTION WITH SQUARE ROOT BAND EDGES MAY BE GENERATED USING CFGPGN OR TERGEN, AND SHOULD BE OF THE SAME LENGTH AS THE ORIGINAL. THE FUNCTION F(E) = SUM 8.0*WB(K)/RNG(K)**2 * Õ E - { AA(K) + RNG(K)*0.5 } - SQRT{E-AA(K)}*SQRT{AA(K)+RNG(K)-E} å IS ASSUMED TO CORRESPOND TO THE SUPPLIED COEFFICIENTS AM(I), BM2(I). T(E) = {S(E,N-1)-F(E)*R(E,N)}/{ S(E,N-2)-F(E)*R(E,N-1)} /BM2(N) N(E) = -1/PI * IM{ ÕQ(E,N-1)-B2(N)*T(E)*Q(E,N-2)å/ ÕP(E,N)-B2(N)*T(E)*P(E,N-1)å } WHERE N=LL AND P,Q AND R,S ARE THE ORTHOGONAL POLYNOMIALS OF THE FIRST AND SECOND KINDS CORRESPONDING TO A,B2 AND AM,BM2 RESPECTIVELY. THIS ROUTINE USES PLYVAL ARGUMENTS : (* INDICATES AN OVERWRITTEN ARGUMENT) DENCRS TAKES THE REQUIRED VALUE E ARGUMENT OF CONTINUED FRACTION A DENOMINATOR COEFFICIENTS OF CONTINUED FRACTION I=1,LL-1 B2 NUMERATOR COEFFICIENTS OF CONTINUED FRACTION I=1,LL LL LENGTH OF CONTINUED FRACTION AA LIST OF BAND LEFT EXTREMA RNG LIST OF BAND WIDTHS WB LIST OF WEIGHTS OF BANDS NB NUMBER OF BANDS (GREATER THAN 0) AM LL-1 DENOMINATOR COEFFICIENTS OF MATCHING CONTINUED FRACTION BM2 LL NUMERATOR COEFFICIENTS OF MATCHING CONTINUED FRACTION ./ ADD NAME=DENINT FUNCTION DENINT(E,A,B2,NA,NP,LL,ALP,EPS,WK,IWK,ICODE) DIMENSION A(NA,NP),B2(NA,NP),WK(LL,4),IWK(LL) EVALUATES THE INTEGRATED DENSITY OF STATES, N(E), CORRESPONDING TO A GIVEN SUM OF CONTINUED FRACTIONS (J-MATRICES) AT A GIVEN POINT E AND RETURNS THAT VALUE, USING THE 'QUADRATURE' APPROACH. THIS ROUTINE USES DENQD,RECWT,RECRTS,NUMC,NUMD ARGUMENTS : (* INDICATES AN OVERWRITTEN ARGUMENT) DENINT TAKES THE COMPUTED VALUE OF THE INTEGRATED DENSITY OF STATES AT E E VALUE AT WHICH INTEGRATED DENSITY OF STATES REQUIRED A* DIAGONAL J-MATRIX ELEMENTS (A(LL,K) OVERWRITTEN) I=1,LL-1 B2 SQUARES OF SUB-DIAGONAL J-MATRIX ELEMENTS I=2,LL B2(1,K) IS THE WEIGHT IN THE K TH BAND NA FIRST DIMENSION OF ARRAYS A AND B2 >= LL NP NUMBER OF DENSITY OF STATES TO BE SUMMED LL LENGTH OF TRIDIAGONALISATIONS ALP PROPORTION OF WEIGHT AT LAST NODE, 0 E A* DIAGONAL J-MATRIX ELEMENTS (A(LL) OVERWRITTEN) I=1,LL-1 B2 SQUARES OF SUB-DIAGONAL J-MATRIX ELEMENTS I=2,LL B2(1) IS THE TOTAL WEIGHT OF THE DENSITY OF STATES LL LENGTH OF TRIDIAGONALISATION ALP PROPORTION OF WEIGHT AT LAST NODE, 0 0 AN REQUIRED VALUE OF INTEGRATES DENSITY OF STATES BETWEEN ZERO AND TOTAL NORMALISATION A(I,J)* DENOMINATOR CONTINUED FRACTION COEFFICIENTS, I=1,LL-1 : J=1,NC. A(LL,J) ONLY ARE OVERWRITTEN B2(I,J) NUMERATOR CONTINUED FRACTION COEFFICIENTS, I=1,LL: J=1,NC ND FIRST DIMENSION OF ARRAYS A AND B2 LL LENGTH OF CONTINUED FRACTION NC NUMBER OF SETS OF COEFFICIENTS IN THE SUM DEFINING THE LDOS ERR ACCURACY REQUIRED IN DETERMINATION OF FERMI ENERGY EPS MACHINE ACCURACY EB* INTERVAL IN WHICH FERMI ENERGY LIES : ON OUTPUT CONTAINS CONFIDENCE INTERVAL IN THE EVENT OF INADEQUATE ACCURACY WK* WORK ARRAY OF DIMENSION AT LEAST LL*4 NW FIRST DIMENSION OF WK : AT LEAST LL IW* INTEGER WORK ARRAY OF LENGTH AT LEAST LL IFT* ON INPUT MAXIMUM NUMBER OF ITERATIONS ALLOWED ON OUTPUT -1 INDICATES GIVEN RANGE OF ENERGIES (EB) INADEQUATE 0 INDICATES DENQD FAILED >0 NUMBER OF ITERATIONS , INPUT VALUE+1 INDICATING FAILURE TO ACHIEVE DESIRED ACCURACY ./ ADD NAME=NUMC FUNCTION NUMC(A,B2,ALAM,LM1) DIMENSION A(LM1),B2(LM1) EVALUATES THE NUMBER OF EIGENVALUES OF A SYMMETRIC TRIDIAGONAL MATRIX STRICTLY GREATER THAN ALAM. THE STURM SEQUENCE PROPERTY IS USED. NO SUB-DIAGONAL ELEMENT MAY BE ZERO.ÕWILKINSONå ARGUMENTS NUMC TAKES THE INTEGER VALUE REQUIRED A DIAGONAL MATRIX ELEMENTS I=1,LM1 B2 SQUARES OF THE SUBDIAGONAL MATRIX ELEMENTS I=2,LM1 ALAM POINT OF EVALUATION LM1 DIMENSION OF MATRIX ./ ADD NAME=NUMD SUBROUTINE NUMD(A,B2,ALAM,LM1,DE) DIMENSION A(LM1),B2(LM1) EVALUATES THE INVERSE OF THE LOGARITHMIC DERIVATIVE OF THE DETERMINANT (SYMMETRIC TRIDIAGONAL MATRIX - ALAM * IDENTITY) WHEN NO SUB-DIAGONAL MATRIX ELEMENT IS ZERO.ÕWILKINSONå ARGUMENTS : ( * INDICATES AN OVERWRITTEN ARGUMENT) A DIAGONAL MATRIX ELEMENTS I=1,LM1 B2 SQUARES OF SUB-DIAGONAL MATRIX ELEMENTS I=2,LM1 ALAM ARGUMENT OF DETERMINANT ABOVE LM1 DIMENSION OF MATRIX DE* DET(MATRIX - ALAM)/DET'(MATRIX - ALAM) ./ ADD NAME=PLYVAL SUBROUTINE PLYVAL(E,A,B2,LM1,P,Q) DIMENSION A(LM1),B2(LM1),P(2),Q(2) COMPUTES RELATIVE VALUES OF LAST TWO POLYNOMIALS OF EACH KIND OF AN ORTHOGONAL SEQUENCE DEFINED BY A THREE-TERM RECURRENCE RELATION. WITH P(X,-1)=Q(X,-1)=0.0 , P(X,0)=1.0 , AND Q(X,1)=B2(1) P(X,I) = (X-A(I)) * P(X,I-1) - B2(I) * P(X,I-2) Q(X,I-1) = (X-A(I)) * Q(X,I-2) - B2(I) * Q(X,I-3) THE SAME NORMALISATION IS USED FOR P AND Q, BUT ITS VALUE VARIES WITH X IN ORDER TO MAINTAIN ACCURACY IN THE RELATIVE VALUES OF THE FUNCTIONS. ARGUMENTS : E ARGUMENT OF POLYNOMIALS A DIAGONAL RECURRENCE COEFFICIENTS I=1,LM1 B2 OFF-DIAGONAL RECURRENCE COEFFICIENTS I=1,LM1 LM1 MAXIMUM DEGREE OF POLYNOMIALS P* P(I) CONTAINS VALUES OF THE POLYNOMIALS OF THE FIRST KIND P(E,LM1+I-2) ,I=1,2 , WITH ARBITRARY NORMALISATION Q* Q(I) CONTAINS VALUES OF THE POLYNOMIALS OF THE SECOND KIND Q(E,LM1+I-3) ,I=1,2 , WITH ARBITRARY NORMALISATION ./ ADD NAME=RECAL SUBROUTINE RECAL(HOP,PSI,PMN,M,A,B2,LL) DIMENSION PSI(M),PMN(M),A(LL),B2(LL) COMPUTES THE TRIDIAGONALISATION OF A LARGE SPARSE SYMMETRIC MATRIX USING THE RECURSION (OR LANCZOS OR PAIGE'S) METHOD : A(N) = B(N+1).PSI(N+1) = (H-A(N)).PSI(N) - B(N).PSI(N-1) = 1 THIS ROUTINE MUST BE SUPPORTED BY A SUBROUTINE SEGMENT HOP, TO CARRY OUT STEPS OF THIS PROCESS : SUBROUTINE HOP(X,Y,A) DIMENSION X( ),Y( ) WHICH, WITH INPUT X=X, Y=Y, WILL PRODUCE Y = H.X+Y AND A = X.H.X = ARGUMENTS : (* DENOTES OVERWRITTEN ARGUMENT) HOP NAME OF SUBROUTINE SEGMENT PSI* INPUT B(0)*PSI(0) OUTPUT B(LL)*PSI(LL) PMN* INPUT -B(0)*PSI(-1) OUTPUT -B(LL)*PSI(LL-1) WHERE THE NORMS OF ALL VECTORS, PSI(I), ARE EQUAL TO UNITY M DIMENSION OF MATRIX A* OUTPUT A(I) = A(I-1), I=1,LL-1 B2* INPUT B2(1) = B(0)**2 OUTPUT B2(I) = B(I-1)**2, I=1,LL LL 1+LENGTH OF TRIDIAGONALISATION REQUIRED >OR= 2 ./ ADD NAME=RECNO SUBROUTINE RECNO(HOP,SOP,U,M,NIT,LS,LL,A,B2,PSI,PMN,EMACH) DIMENSION U(M),PSI(M),PMN(M),A(LL),B2(LL) DOUBLE PRECISION SUM COMMON /BLKNNN/SUM IMPLEMENTS THE 'NON-ORTHOGONAL BASIS' RECURSION METHOD WHERE THE EIGEN-PROBLEM TAKES THE FORM M X = E S X . S IS ASSUMED TO HAVE UNIT DIAGONAL ELEMENTS AND BOTH M AND S TO BE REAL SYMMETRIC. THE INVERSE OF S TIMES A VECTOR IS ESTIMATED BY TAKING NIT APPLICATIONS OF GAUSS-SEIDEL ITERATION. IF NIT IS SET TO BE ZERO THE ROUTINE PERFORMS THE USUAL RECURSION ASSUMING S IS THE IDENTITY. THE NUMBER OF VECTORS NEEDED IS KEPT TO A MINIMUM (3) BY ASKING THE USER TO WRITE THE S MULTIPLICATION ROUTINE IN A PARTICULAR WAY, SO PLEASE NOTE THE SPECIFICATIONS CAREFULLY. IF S IS THE IDENTITY THEN NIT SHOULD BE SET TO ZERO AND THE VECTORS U AND PSI SET TO REFER TO THE SAME REAL ARRAY. THE GREENIAN CALCULATED MAY BE THOUGHT OF AS U(TRANS) S (S E - M)INVERSE S U WHERE U IS THE STARTING VECTOR. THE VECTORS PSI AND PMN ARE IN FACT S TIMES THE USUAL RECURSION BASIS VECTORS. THE VECTOR S U MAY BE SPECIFIED RATHER THAN U AND THIS INDICATED BY REPLACING NIT BY -NIT AND STORING THAT VECTOR IN PSI ON THE FIRST CALL TO RECNO. THE ROUTINE MAY BE RESTARTED TO EXTEND THE NUMBER OF LEVELS AS U,PSI AND PMN AND THE NORM SUM ARE ALL RETAINED (THE LATTER IN THE COMMON BLOCK /BLKNNN/ THE ACTUAL ALGORITHM IS AS FOLLOWS : PMN(0) = 0 PSI(0) = S U(0) B2(1) = U(0)(TRANSPOSED) PSI(0) UP TO NUMBER OF LEVELS (LL) DO : A(N) = U(N)TRANSPOSED M U(N) / U(N)TRANSPOSED PSI(N) PSI(N+1) = M U(N) -A(N) PSI(N) - B2(N) PSI(N-1) U(N+1) = S(INVERSE) PSI(N+1) HERE S(INVERSE) PSI(N+1) IS CALCULATED BY NIT APPLICATIONS OF THE FORMULA (I IS THE ITERATION NUMBER) : U(N+1)(I) = PSI(N+1) - L U(N+1)(I) - R U(N+1)(I-1) WHERE L AND R ARE THE STRICT LEFT AND RIGHT TRIANGULAR PARTS OF S B2(N+1)= U(N+1)TRANSPOSED PSI(N+1) / U(N)TRANSPOSED PSI(N) WITH A RENORMALISATION OF THE VECTORS FOR NUMERICAL STABILITY. CHRIS NEX JAN.1985 ARGUMENTS : (* INDICATES AN OVERWRITTEN ARGUMENT) HOP NAME OF A SUBROUTINE TO PERFORM THE MATRIX MULTIPLICATION BY M . ITS SPECIFICATIONS MUST BE AS FOLLOWS: SUBROUTINE HOP(X,Y,A) DIMENSION X( ),Y( ) WITH INPUT X=X, Y=Y, WILL PRODUCES Y = M.X+Y AND A = X.M.X = SOP NAME OF A SUBROUTINE TO EVALUATE THE PRODUCT OF THE OFF- DIAGONAL ELEMENTS OF S WITH A VECTOR. THE SPECIFICATION IS AS FOLLOWS SUBROUTINE SOP(U,V,W) DIMENSION U( ),V( ),W( ) CALCULATES W = V - OFF-DIAGONAL(S) U THIS IS CALLED WITH U AND W REFERRING TO THE SAME ARRAY TO ACHIEVE A GAUSS-SEIDEL STEP AND WITH V AND W REFERRING TO THE SAME ARRAY TO PERFORM A MORE USUAL MATRIX MULTIPLICATION. N.B. NOTE THE MINUS SIGN THE IMPLICATION FOR THE USER IS THAT THE ELEMENTS OF THE PRODUCT MUST BE EVALUATED AND OVERWRITTEN IN INCREASING ORDER, NOT BY A GLOBAL ACCUMULATION TECHNIQUE. U* INPUT AS STARTING STATE IF NIT >OR= 0 OUTPUT AS LAST RECURSION BASIS VECTOR M DIMENSION OF MATRIX NIT THE ABSOLUTE VALUE DEFINES THE NUMBER OF ITERATIONS IN THE GAUSS-SEIDEL PROCESS. THIS IS RATHER BETTER THAN NIT TERMS IN THE TAYLOR SERIES FOR THE INVERSE OF S. IF NIT IS ZERO S IS ASSUMED TO BE THE IDENTITY MATRIX AND U AND PSI MAY REFER TO THE SAME PHYSICAL ARRAY. IF NIT IS NEGATIVE IT IS ASSUMED THAT PSI DEFINES THE STARTING STATE, AND U WILL BE CALCULATED BY RECNO USING GAUSS-SEIDEL ITERATION TO ESTIMATE S INVERSE PSI. LS STARTING LEVEL IN THE RECURSION PROCESS. IF LS=1 THE NORMALISATION B2(1) IS CALCULATED AS U(TRANS)S U. IF GREATER THAN 1, IT IS ASSUMED THAT ALL ARGUMENTS ARE UNCHANGED FROM A PREVIOUS CALL TO RECNO LL* ON INPUT LAST REQUIRED LEVEL IN TRIDIAGONALISATION ON OUTPUT, IF NEGATIVE, -LL IS THE NUMBER OF LEVELS REACHED BEFORE B2(LL)0 ). THE METHOD USED IS BISECTION BASED ON THE STURM SEQUENCE PROPERTY FOLLOWED BY NEWTON'S METHOD FOR ISOLATED ROOTS.ÕWILKINSONå THIS ROUTINE USES NUMC,NUMD ARGUMENTS : (* INDICATES AN OVERWRITTEN ARGUMENT) A DIAGONAL MATRIX ELEMENTS I=1,LM1 B2 SQUARES OF SUB-DIAGONAL MATRIX ELEMENTS I=2,LM1 LM1 DIMENSION OF MATRIX EPS ABSOLUTE ACCURACY REQUIRED IN EIGENVALUES XLIM UPPER BOUND ON EIGENVALUES TO BE FOUND (IF RELEVANT) N* IF 0 ON INPUT : ONLY EIGENVALUES LESS THAN XLIM ARE FOUND ON OUTPUT : NUMBER OF DISTINCT EIGENVALUES FOUND X* EIGENVALUES IN ASCENDING ORDER MULT* MULTIPLICITY OF EACH EIGENVALUE IF NEGATIVE THEN THE CORRESPONDING EIGENVALUE WAS FOUND WITH LESS ACCURACY THAN EPS BI* REAL WORK ARRAY OF LENGTH AT LEAST LM1 NI* INTEGER WORK ARRAY OF LENGTH AT LEAST LM1 ./ ADD NAME=RECSUM SUBROUTINE RECSUM(AC,BC,NA,LL,NP,A,B2,EPS,WK,NW) DIMENSION AC(NA,NP),BC(NA,NP),WK(NW),A(NA),B(NA) COMPUTES THE TRIDIAGONALISATION (CONTINUED FRACTION, JACOBI MATRIX) CORRESPONDING TO THE SUM OF NP TRIDIAGONALISATIONS, W(M,X) : SQRT(B(N+1,M)).P(N+1,M,X) = (X-A(N,M)).P(N,M,X) - SQRT(B(N,M)).P(N-1,M,X) W(X) = SIGMA(M:1,NC) ¡B(0,M).W(M,X)¢ NOTE THAT THIS ROUTINE USES RECQD,CFGEN,RECRTS,NUMC,NUMD ARGUMENTS : AC AC(N,M) = A(N-1,M), N=1,M, M=1,NP BC BC(N,M) = B(N-1,M)**2, N=1,LL, M=1,NP NA FIRST DIMENSION OF ARRAYS AC AND BC IN CALLING PROGRAM LL* ON INPUT : THE ABSOLUTE VALUE GIVES LENGTH +1 OF EACH TRIDIAGONALISATION. IF >0 M=LL-1 ; IF <0 M=LL ON OUTPUT: LENGTH OF OUTPUT TRIDIAGONALISATION , IF NEGATIVE THEN RECQD FAILED WITH TOO FEW ROOTS NP NUMBER OF CONTINUED FRACTIONS A* A(I) = A(I-1), I=1,M, IN TRIDIAGONALISATION CORRESPONDING TO W(X) B2* B2(I) = B(I-1)**2, I=1,LL, IN TRIDIAGONALISATION CORRESPONDING TO W(X) EPS ACCURACY REQUIRED IN COMPUTATION WK* REAL WORK ARRAY OF LENGTH AT LEAST 5*LL*NP NW LENGTH OF ARRAY WK ./ ADD NAME=RECWT FUNCTION RECWT(E,A,B2,LL,EPS,N,P,NS) DIMENSION A(LL),B2(LL),P(2,3) COMPUTES THE VALUE OF THE WEIGHT AT THE FIXED POINT IN A 1- FIXED POINT GAUSSIAN QUADRATURE, GIVEN THE CORRESPONDING 3-TERM RECURRENCE RELATION: P(E,J)= (E-A(J)) * P(E,J-1) - B2(J) * P(E,J-2) THIS ROUTINE MAY BE CALLED REPEATEDLY WITH INCREASING NUMBER OF LEVELS SUCH THAT IT DOES NOT RECOMPUTE EARLIER POLYNOMIAL VALUES. IF REQUIRED THE VALUE OF THE LAST COEFFICIENT ,A(LL), MAY BE COMPUTED, OR IT MAY BE ASSUMED THAT THIS HAS ALREADY BEEN DONE AND THAT VALUE USED IN THE CALCULATION OF THE WEIGHT. THE EXPRESSION FOR THE WEIGHT USED IS (WITH N=LL) { P(E,N-1)*Q(E,N-1) - P(E,N)*Q(E,N-2) / { P(E,N-1)*P'(E,N) - P'(E,N-1)*P(E,N) + P(E,N)**2/B2(N) } AS THIS FORM IS INDEPENDENT OF THE NORMALISATION OF THE POLYNOMIALS. P AND Q ARE THE MONIC POLYNOMIALS OF THE FIRST AND SECOND KINDS. ARGUMENTS : (* INDICATES AN OVERWRITTEN ARGUMENT) RECWT TAKES THE VALUE OF THE REQUIRED WEIGHT E REQUIRED FIXED POINT IN QUADRATURE. IT MAY BE A NODE OF THE LL-1 OR LL QUADRATURE IF A(LL) IS APPROPRIATELY DEFINED A* DIAGONAL ELEMENTS OF THE RECURRENCE. IF N IS CHANGED FROM -1 INPUT TO 0 ON OUTPUT THEN A(LL) CONTAINS THE ADJUSTED VALUE TO ACHIEVE A GAUSSIAN NODE AT E, OTHERWISE A IS UNCHANGED. B2 SUB-DIAGONAL ELEMENTS OF THE RECURRENCE LL INDEX OF LAST B2 VALUE TO BE USED EPS RELATIVE THRESHOLD VALUE OF THE POLYNOMIAL BELOW WHICH E WILL BE ACCEPTED AS A ZERO N* CODE : -1 A(LL) TO BE OVERWRITTEN. N CHANGED TO 0 IF SUCCESSFUL, UNCHANGED OTHERWISE 0 A(LL) GIVEN (NOT OVERWRITTEN) 1 A(LL) NOT COMPUTED EXPLICITLY (NOT OVERWRITTEN) P* FINAL POLYNOMIAL VALUES USED IN CALCULATION OF WEIGHT TO BE USED UNCHANGED IF ROUTINE IS RE-ENTERED WITH NS=LL IF N=LL-IABS(N) P(2,1)=P(E,N) P(1,1)=P(E,N-1) P(2,2)=P'(E,N) P(1,2)=P'(E,N-1) P(2,3)=Q(E,N-1) P(1,3)=Q(E,N-2) Q(E,M) IS THE POLYNOMIAL OF THE SECOND KIND OF DEGREE M NS POINT AT WHICH RECURRENCE IS INITIATED . THIS SHOULD BE 1 INITIALLY , BUT FOR A SUBSEQUENT CALL, WITH E UNCHANGED AND LARGER LL, SHOULD BE SET TO THE CURRENT VALUE OF LL ./ ADD NAME=TABAN SUBROUTINE TABAN(E,WT,NPTS,THU,THL,ET,IC,WTT,NET) DIMENSION E(NPTS),WT(NPTS),ET(NET),IC(NET),WTT(NET) IDENTIFIES EXTREMAL VALUES OF A TABULATED FUNCTION, WITHIN A USER-DEFINED RANGE OF FUNCTION VALUES. THIS IS EFFECTED BY SIMPLE COMPARISON OF THE TABULAR VALUES. ARGUMENTS : (* INDICATES AN OVERWRITTEN ARGUMENT) E LIST OF ORDINATE VALUES WT LIST OF ABSCISSAE NPTS NUMBER OF COORDINATE PAIRS THU UPPER VALUE OF 'WINDOW' ON FUNCTION THL LOWER VALUE OF 'WINDOW' ON FUNCTION ET* LIST OF ORDINATES OF EXTREMA WITHIN WINDOW IC* CODE OF TYPE OF EXTREMA : A + SIGN INDICATES INCREASING FUNCTION VALUES AND -, DECREASING ONES. AT LOCAL EXTREMA THIS EXTENDS TO + INDICATING A MINIMUM AND - A MAXIMUM. THE ABSOLUTE VALUE IS CODED BELOW. 1 FUNCTION VALUE CROSSES LOWER THRESHOLD 2 FUNCTION VALUE CROSSES UPPER THRESHOLD 3 FUNCTION VALUE CROSSES BOTH THRESHOLDS 4 FUNCTION VALUE IS A LOCAL EXTREMUM WTT* IF IABS(IC(I))=4 THIS GIVES THE LOCAL EXTREMUM VALUE NET* NUMBER OF VALUES TABULATED IN ET,IC,WTT IF NEGATIVE ON OUTPUT, THEN THERE IS NOT ENOUGH SPACE AND THE ABSOLUTE VALUE GIVES THE INDEX OF THE LAST COORDINATE PAIR EXAMINED. THE RESULTS TO THAT POINT ARE STORED AS ABOVE. ./ ADD NAME=TERMGN SUBROUTINE TERMGN(A,B2,LL,EPS,ERR,ITMX,AA,RNG,WB,NBP1,AM,BM2, 1 IC,WK,NW,BWK,NBD,IWK) DIMENSION A(LL),B2(LL),AA(NBP1),RNG(NBP1),WB(NBP1),AM(LL) 1 ,BM2(LL),IC(NW),WK(NW,2,NBP1),BWK(NBD,4),IWK(LL) GENERATES AN ANALYTIC TERMINATOR TO A GIVEN CONTINUED FRACTION. THE FORM OF THE TERMINATOR IS A SUM OF SQUARE ROOTS OF QUADRATICS , F(E) ,AS IN DENCRS, WITH PARAMETERS TO BE ADJUSTED TO MATCH THE APPARENT BAND GAPS IN THE GIVEN CONTINUED FRACTION. THE LOCAL WEIGHT (AS CALCULATED IN RECWT) OF F(E) IS MATCHED TO THAT OF THE GIVEN CONTINUED FRACTION (A(I),B2(I)) AT E VALUES IN THE NEIGHBOURHOOD OF BAND EDGES AND LOCAL MINIMA. THIS ROUTINE MAY SERVE AS AN EXAMPLE FOR THE MATCHING OF OTHER FORMS OF TERMINATING FUNCTION, OR MATCHING ALGORITHMS. THIS ROUTINE USES TABAN,BNDCRD,BNDEST,BNDREF,CFGPGN,RECWT,RECQD,RECRTS,NUMC,NUMD,WTMIN ARGUMENTS : (* INDICATES AN OVERWRITTEN ARGUMENT) A DIAGONAL RECURSION COEFFICIENTS I=1,LL-1 B2 OFF-DIAGONAL RECURSION COEFFICIENTS I=1,LL LL* LENGTH OF GIVEN RECURSION . ON OUTPUT CONTAINS THE LENGTH OF THE COMPUTED CONTINUED FRACTION WHICH IF DIFFERENT FROM INPUT INDICATES FAILURE OF CFGPGN EPS MACHINE ACCURACY ERR* ACCURACY REQUIRED IN LOCATION OF BAND EDGES , ON OUTPUT THE ESTIMATED ACCURACY, SUBJECT TO ITMX MAXIMUM NUMBER OF ITERATIONS IN LOCATION AA* LIST OF BAND LEFT EXTREMA RNG* LIST OF BAND WIDTHS WB* LIST OF BAND WEIGHTS NBP1* 1+NUMBER OF BANDS ,MAXIMUM ON INPUT AND ON OUTPUT CONTAINS THE 1+NUMBER COMPUTED UNLESS THIS EXCEEDS THE INPUT NUMBER WHEN A NEGATIVE VALUE INDICATES THE NUMBER OF BANDS IDENTIFIED BUT NOT COMPUTED. A ZERO VALUE INDICATES A FAILURE IN THE SEARCH PROCEDURE.(INCREASING NW MAY HELP) AM* DIAGONAL C.F. COEFFICIENTS OF MATCHING FUNCTION BM2* OFF-DIAGONAL C.F. COEFFICIENTS OF MATCHING FUNCTION IC* WORK ARRAY OF LENGTH AT LEAST NW WK* WORK ARRAY OF LENGTH AT LEAST LL*2*MAX(3,NBP1) NW FIRST DIMENSION OF WK. NO.OF POINTS USED IN INITIAL SCAN FOR BAND EXTREMA BWK* WORK ARRAY OF MATCHING POINTS OF DIMENSION AT LEAST 8*NBP1 NBD FIRST DIMENSION OF BWK : AT LEAST 2*NBP1 IWK* INTEGER WORK ARRAY OF LENGTH AT LEAST LL ./ ADD NAME=WTMIN SUBROUTINE WTMIN(AA,BB,A,B2,LL,EPS,ACC,ITMX,EM,FEM) DIMENSION A(LL),B2(LL),P(2,3),E(4),FE(4) FINDS A LOCAL MINIMUM IN THE LOCAL WEIGHT, AS DEFINED IN RECWT, OF A DENSITY OF STATES BY SIMPLE SUBDIVISION OF A GIVEN INTERVAL. THE DENSITY OF STATES IS DEFINED BY ITS CONTINUED FRACTION OR J-MATRIX THIS ROUTINE CALLS RECWT ARGUMENTS (* INDICATES AN OVERWRITTEN ARGUMENT) AA LEFT END OF INTERVAL IN WHICH A MINIMUM IS LOCATED BB RIGHT END OF INTERVAL IN WHICH A MINIMUM IS LOCATED A DENOMINATOR COEFFICIENTS OF CONTINUED FRACTION; I=1,LL-1 B2 NUMERATOR COEFFICIENTS OF CONTINUED FRACTION; I=1,LL LL LENGTH OF CONTINUED FRACTION EPS MACHINE ACCURACY ACC ABSOLUTE ACCURACY REQUIRED IN LOCATION OF MINIMUM SUBJECT TO ITMX* MAXIMUM NUMBER OF FUNCTION EVALUATIONS ON OUTPUT: -2 IF MINIMUM OF THREE EQUIDISTANT POINTS IS AT BB -1 IF MINIMUM OF THREE EQUIDISTANT POINTS IS AT AA 0 IF INSUFFICIENT ACCURACY AFTER ITMX FUNCTION EVALUATIONS OTHERWISE CONTAINS THE NUMBER OF FUNCTION EVALUATIONS USED EM* COMPUTED LOCATION OF MINIMUM FEM* COMPUTED MINIMUM WEIGHT ./ ENDUP ./ ADD NAME=ADDAT SUBROUTINE ADDAT(CRD,ND,NAT,EV,IZP,MM,NN,NND,NM,NGBR,NE,EE,NP 1,VEC,IW,NED,OVPAR,HCAL) DIMENSION CRD(ND,NAT),VEC(3,NED),V(3),EE(NP,NP,NED) INTEGER*2 IZP(NAT),NN(NND,NM),IW(2,NED),MM(NND,NM) LOGICAL EV EXTERNAL EV,NGBR,OVPAR EXTENDS THE HAMILTONIAN MATRIX FROM THE USER SUPPLIED ROUTINES EV, HCAL, NGBR AND IOVPAR, AND THE LIBRARY ROUTINES NNCAL & MMCAL . THIS ASSUMES IT HAS ALREADY BEEN SET UP BY SUBROUTINE SETUP IN THE ARRAYS MM,NN,EE,VEC AND IW . ARGUMENTS OF ADDAT : (* INDICATES OVERWRITTEN BY THE ROUTINE) CRD LATTICE COORDINATES ND FIRST DIMENSION OF CRD NAT NO.OF ATOMS IN THE CLUSTER EV LOGICAL FUNCTION OF 2 ARGUMENTS, BOTH REAL ARRAYS OF LENGTH 3 RETURNING THE VALUE .TRUE. IF THE ARRAYS ARE EQUIVALENT AND .FALSE. IF NOT. IZP THE ABSOLUTE VALUE GIVES 'TYPE' OF EACH ATOM IF THE SIGN IS + THEN THE ATOM IS ASSUMED PART OF THE ORIGINAL CLUSTER IF THE SIGN IS - THEN THE ATOM HAS ITS CONNECTIVITY AND INTERACTIONS COMPUTED MM* IS THE INTERACTION MAP GENERATED BY MMCAL NN* IS THE NEIGHBOUR MAP GENERATED BY NNCAL NND FIRST DIMENSION OF ARRAYS MM & NN NM* MAX NO. OF ATOMS CONNECTED BY INTERACTIONS. ON OUTPUT CONTAINS ACTUAL MAX NO. GENERATED NGBR NAME OF A FUNCTION TO SUPPLY INTERACTION INFORMATION TO NNCAL ARGUMENTS ARE : II 'TYPE' OF ATOM I JJ 'TYPE' OF ATOM J R2 SQUARE OF DISTANCE FROM I TO J DD DUMMY ARGUMENT NGBR TAKES THE VALUE 1 IF I & J ARE NEIGHBOURS AND 0 OTHERWISE NE* NO. OF DISTINCT DISPLACEMENT VECTORS (MATRICES) ALREADY FOUND ON OUTPUT CONTAINS THE NEW TOTAL NUMBER FOUND EE* LIST OF INTERACTION MATRICES NP FIRST 2 DIMENSIONS OF ARRAY EE VEC* LIST OF DISTINCT DISPLACEMENT VECTORS FOUND (POSN. I - POSN.J) IW* LIST OF ATOM TYPES AT THE ENDS OF THE VECTORS IN VEC IW(1,.) IS TYPE OF I IW(2,.) IS TYPE OF J NED LAST DIMENSION OF ARRAYS EE,IW,VEC OVPAR NAME OF A FUNCTION TO SUPPLY OVERLAP PARAMETERS TO HCAL ARGUMENTS ARE II 'TYPE' OF ATOM I JJ 'TYPE' OF ATOM J R2 SQUARE OF THE DISTANCE FROM I TO J DD* OVERLAP PARAMETERS AS REQUIRED BY HCAL THE NOTATION USED IS AS FOLLOWS: DD(1) DD SIGMA DD(2) DD PI DD(3) DD DELTA DD(4) PD SIGMA DD(5) PD PI DD(6) PP SIGMA DD(7) PP PI DD(8) SD SIGMA DD(9) SP SIGMA DD(10) SS SIGMA DD(11) D SELF ENERGY DD(12) P SELF ENERGY DD(13) S SELF ENERGY HCAL NAME OF A SUBROUTINE TO CALCULATE THE INTERACTION BETWEEN TWO ATOMS. ARGUMENTS ARE V VECTOR POSITION(I) - POSTITION(J) II TYPE AT I JJ TYPE AT J E* OUTPUT INTERACTION MATRIX H OPERATING ON PSI(J) EFFECT AT I IOVPAR NAME OF FUNCTION SUPPLING INFORMATION TO HCAL ./ ADD NAME=BCCBFE INTEGER FUNCTION BCCBFE(I,J,R2,DD) DIMENSION DD(13) DETERMINES WHETHER A DISTANCE IS A 'NEAREST NEIGHBOUR' OR 'NEXT NEAREST NEIGHBOUR' DISTANCE IN THE BCC LATTICE GENERATED BY BCCLAT , AND IF SO OUTPUTS THE DD PARAMETERS FOR IRON ACCORDING TO D.G.PETTIFOR . (IN RYDBERGS) ARGUMENTS : I 'TYPE' OF ONE LATTICE SITE J 'TYPE' OF THE OTHER LATTICE R2 SQUARE OF THE DISTANCE BETWEEN THE TWO LATTICE SITES DD* OUTPUT AS THE DD PARAMETERS OF D.G.PETTIFOR (SIGMA,PI,DELTA) AND DD(11)=0.0 OF R2<1.0E-4 AS THE SELF ENERGY BCCBFE TAKES THE VALUE 0 IF THE SITES ARE NOT NEIGHBOURS AND 1 IF THEY ARE NEIGHBOURS ./ ADD NAME=BCCLAT SUBROUTINE BCCLAT(CRD,NDIM,IZP,NAT,NX,NY,NZ,NTYPE) DIMENSION CRD(NDIM,NAT) INTEGER*2 IZP(NAT) GENERATES A BCC LATTICE ON POSITIVE INTEGER GRID, ENCLOSED BY A CUBOID OF A GIVEN SIZE. ARGUMENTS:( * INDICATES AN OVERWRITTEN ARGUMENT) CRD* LATTICE COORDINATES ((I,J),I=1,3),J=1,NAT NDIM FIRST FIRST DIMENSION OF ARRAY COORD >OR= 3 IZP* INTEGER*2 ARRAY RETURNING THE VALUE NTYPE IN EACH ELEMENT NAT* ON INPUT THE MAXIMUM NUMBER OF LATTICE POINTS ALLOWED ON OUTPUT THE ACTUAL NUMBER OF POINTS GENERATED NX,NY,NZ INTEGER DIMENSIONS OF THE CUBOID TO CONTAIN THE LATTICE NTYPE 'TYPE' CODE FOR EACH LATTICE SITE ./ ADD NAME=EQUIV LOGICAL FUNCTION EQUIV(V,W) DIMENSION V(3),W(3) DETERMINES WHETHER TWO VECTORE ARE 'EQUIVALENT' (FOR THE PURPOSES OF THE SUBROUTINE MMCAL). INPUT : V REAL ARRAY OF LENGTH 3 W REAL ARRAY OF LENGTH 3 OUTPUT : EQUIV TAKES THE VALUE .TRUE. IF THE VECTORS ARE 'EQUIVALENT' AND .FALSE. OTHERWISE ./ ADD NAME=FCCBND INTEGER FUNCTION FCCBND(I,J,R2,DD) DIMENSION DD(3) DETERMINES WHETHER A DISTANCE IS A 'NEAREST NEIGHBOUR' DISTANCE IN THE FCC LATTICE GENERATED BY FCCLAT , AND IF SO OUTPUTS THE DD PARAMETERS ACCORDING TO D.G.PETTIFOR . (IN RYDBERGS) ARGUMENTS : I 'TYPE' OF ONE LATTICE SITE J 'TYPE' OF THE OTHER LATTICE R2 SQUARE OF THE DISTANCE BETWEEN THE TWO LATTICE SITES DD* OUTPUT AS THE DD PARAMETERS OF D.G.PETTIFOR (SIGMA,PI,DELTA) AND DD(11) IS OUTPUT AS 0.0 (THE SELF-ENERGY) IF R2<1.0E-4 FCCBND TAKES THE VALUE 0 IF THE SITES ARE NOT NEIGHBOURS AND 1 IF THEY ARE NEIGHBOURS ./ ADD NAME=FCCLAT SUBROUTINE FCCLAT(CRD,NDIM,IZP,NAT,NX,NY,NZ,NTYPE) INTEGER*2 IZP(NAT) DIMENSION CRD(NDIM,NAT) GENERATES A FCC LATTICE ON POSITIVE INTEGER GRID, ENCLOSED BY A CUBOID OF A GIVEN SIZE. ARGUMENTS:( * INDICATES AN OVERWRITTEN ARGUMENT) CRD* LATTICE COORDINATES ((I,J),I=1,3),J=1,NAT NDIM FIRST FIRST DIMENSION OF ARRAY COORD >OR= 3 IZP* INTEGER*2 ARRAY RETURNING THE VALUE NTYPE IN EACH ELEMENT NAT* ON INPUT THE MAXIMUM NUMBER OF LATTICE POINTS ALLOWED ON OUTPUT THE ACTUAL NUMBER OF POINTS GENERATED NX,NY,NZ INTEGER DIMENSIONS OF THE CUBOID TO CONTAIN THE LATTICE NTYPE 'TYPE' CODE FOR EACH LATTICE SITE ./ ADD NAME=MMCAL SUBROUTINE MMCAL(CRD,NDIM,NAT,NN,ND,NM,EV,IZP,NMAT,MM,VEC,IW) INTEGER*2 NN(ND,NM),MM(ND,NM),IZP(NAT),IW(2,NMAT) DIMENSION CRD(NDIM,NAT),VEC(NDIM,NMAT) LOGICAL EV COMMON /BLKNNM/NNMAT COMPUTES AN INDEX OF DISTINCT VECTORS LINKING NEIGHBOURING SITES IN A GIVEN LATTICE. THE VECTORS ARE COMPUTED AND INDEXED ACCORDING TO THE 'TYPE' (AS DEFINED BY IZP) OF THE TERMINAL ATOMS AS WELL AS BY THE VECTOR COMPONENTS. THUS IF THERE ARE 3 TYPES OF ATOMS LINKED IN ALL PAIR COMBINATIONS BY EQUIVALENT VECTORS, ALL COMBINATIONS WILL OCCUR IN THE INDEX. (I.E. 12 ENTRIES INCLUDING BOTH SENSES OF THE VECTOR) IF ANY OF THE 'TYPES' IN IZP ARE NEGATIVE, IT IS ASSUMED THAT MMCAL HAS ALREADY BEEN CALLED FOR A SUBCLUSTER OF THE CURRENT CLUSTER AND THAT THOSE ATOMS WITH NEGATIVE IZP ARE NEW ADDITIONS WHOSE INTERACTIONS ARE TO BE COMPUTED. ARGUMENTS : (* INDICATES OVERWRITING BY THE SUBROUTINE) CRD(I,J) COORDINATES OF THE LATTICE (I=1,NDIM) ,J=1,NAT NDIM FIRST DIMENSION OF ARRAYS CRD AND VEC NAT NUMBER OF SITES IN THE LATTICE NN NEAREST NEIGHBOUR MAP AS CALCULATED BY NNCAL : NN(I,1)=1+NO.OF NEIGHBOURS OF SITE I NN(I,J),J=2,NN(I,1) LISTS THE NEIGHBOURS OF SITE I ND FIRST DIMENSION OF ARRAY NN NM SECOND DIMENSION OF ARRAY NN EV LOGICAL FUNCTION (DECLARED EXTERNAL IN THE CALLING ROUTINE) WITH 2 ARGUMENTS ,EACH A REAL ARRAY OF LENGTH NDIM, RETURNING THE VALUE .TRUE. IF ITS ARGUMENTS ARE THE 'SAME' AND .FALSE. IF NOT. THE ARGUMENTS MUST BE UNCHANGED. IZP IZP(I) ABSOLUTE VALUE GIVES 'TYPE' OF I TH LATTICE SITE IF ATOMS ARE BEING ADDED TO AN EXISTING CLUSTER THE A NEGATIVE SIGN INDICATES AN ADDED ATOM. NMAT* ON A FIRST CALL THE MAXIMUM NUMBER OF DISTINCT VECTORS ALLOWED. SUBSEQUENTLY THE NUMBER PREVIOUSLY CALCULATED(AS O/P) ON OUTPUT THE ACTUAL NUMBER OF VECTORS CALCULATED IF 0 THEN NOT ENOUGH STORE HAS BEEN ALLOWED AND NMAT MUST BE INCREASED. MM* INDEX OF VECTORS LINKING NEIGHBOURING SITES: MM(I,J)= K, THE INDEX OF THE VECTOR STORED IN VEC SUCH THAT VEC(K)=SITE VECTOR(NN(I,J)) - SITE VECTOR(I) ,J=2,NN(I,1) VEC(R,K)* LIST OF DISTINCT VECTORS ,(R=1,NDIM) , K=1,NMAT IW(1,K)* 'TYPE' OF ATOM I AT ONE END OF THE K TH VECTOR IW(2,K)* 'TYPE' OF ATOM J AT THE OTHER END OF THE K TH VECTOR ./ ADD NAME=NNCAL SUBROUTINE NNCAL(CRD,NDIM,NAT,IZP,NN,ND,NM,NGBR) INTEGER*2 IZP(NAT),NN(ND,NM) DIMENSION CRD(NDIM,NAT),DUM(3) CALCULATES THE 'NEAREST NEIGHBOUR' MAP OF A LATTICE, GIVEN A SUBROUTINE DEFINING 'NEIGHBOUR'. IT ALSO EXTENDS A MAP GENERATED BY A PREVIOUS CALL, IN WHICH CASE ADDED ATOMS ARE INDICATED BY A NEGATIVE VALUE OF IZP. ARGUMENTS: (* INDICATES OVERWRITING BY THE SUBROUTINE) CRD(I,J) LATTICE COORDINATES (I=1,3),J=1,NAT NDIM FIRST DIMENSION OF ARRAY CRD >OR= 3 NAT NUMBER OF LATTICE POINTS IZP INTEGER*2 ARRAY LISTING THE 'TYPE' OF EACH SITE (FOR NGBR) IF IZP(I) IS NEGATIVE THE ABSOLUTE VALUE IS TAKEN AND ONLY THOSE ATOMS WITH NEGATIVE IZP ARE CONSIDERED FOR MODIFICATIONS TO NN NN* 'NEAREST NEIGHBOUR MAP' : NN(I,1) = 1+NUMBER OF NEIGHBOURS OF SITE I NN(I,J),J=2,NN(I,1) LIST OF SITES CONNECTED TO SITE I ND FIRST DIMENSION OF ARRAY NN NM* SECOND DIMENSION OF ARRAY NN (MAX. NO. OF NEIGHBOURS +1) ON OUTPUT CONTAINS ACTUAL MAX.NO. OF NEIGHBOURS +1 NGBR THE NAME OF A FUNCTION (DECLARED EXTERNAL) OF 4 ARGUMENTS RETURNING THE VALUE 1 IF THE SITES ARE CONNECTED AND 0 IF NOT. THE ARGUMENTS ARE AS FOLLOWS: II 'TYPE' OF SITE A JJ 'TYPE' OF SITE B R2 THE SQUARE OF THE DISTANCE FROM A TO B DUM A DUMMY ARGUMENT NONE OF THESE ARGUMENTS MUST BE MODIFIED BY NGBR ./ ADD NAME=ONION SUBROUTINE ONION(NN,ND,NM,IZERO,NAT,IST,NNS,IW) INTEGER*2 NN(ND,NM),IZERO(NAT),IST(NNS),IW(NAT) ASSIGNS EACH SITE IN A LATTICE (DEFINED BY A 'CONNECTIVITY MAP') TO A SHELL DEFINED BY A 'TOPOLOGICAL' (NUMBER OF 'HOPS') DISTANCE FROM A GIVEN GROUP OF SITES.THE GIVEN GROUP IS LABELLED SHELL 1. ARGUMENTS :(* INDICATES OVERWRITING BY THE SUBROUTINE) NN NEIGHBOUR MAP AS DEFINED BY NNCAL ND FIRST DIMENSION OF ARRAY NN NM SECOND DIMENSION OF ARRAY NN IZERO* INTEGER*2 ARRAY RETURNING THE SHELL NUMBER OF EACH SITE NAT NUMBER OF LATTICE SITES IST INTEGER*2 ARRAY INDEXING THE 'CENTRAL' SITE(S) NNS NUMBER OF CENTRAL SITES IW INTEGER*2 WORK ARRAY OF LENGTH AT LEAST NAT ./ ADD NAME=ORPEEL SUBROUTINE ORPEEL(NSTRT,NORB,NO,MM,NN,ND,ID,EE,NP,NE,NED,MEM) INTEGER*2 MM(ND,ID),NN(ND,ID),MEM(3,ID) DIMENSION EE(NP,NP,NED) IMPLEMENTS ORBITAL PEELING AS SPECIFIED IN THE PH.D THESIS OF N.R. BURKE .AN EQUIVALENT (FUNCTIONAL) DEFINITION IS THAT THE SUBROUTINE DELETES A ROW AND COLUMN OF A SPARCE MATRIX AS SET UP USING NNCAL AND MMCAL. THE MATRIX IS ASSUMED PARTITIONED INTO NP BY NP BLOCKS, OF WHICH THERE ARE ONLY RELATIVELY FEW DISTINCT ONES IN THE OVERALL MATRIX.TO DELETE A ROW AND COLUMN, THEREFORE , A COPY IS MADE OF THE BLOCKS INVOLVED AND THE LIST OF SUBMATRICES MODIFIED ACCORDINGLY. IT IS ASSUMED THAT THE OVERALL PURPOSE IS TO DELETE ROWS AND COLUMNS DEFINED BY A GIVEN DIAGONAL SUBMATRIX. ARGUMENTS : (* INDICATES AN OVERWRITTEN ARGUMENT) NSTRT THE STARTING ATOM .(DIAGONAL SUBMATRIX TO BE DELETED) NORB ORBITAL TO BE PEELED (ROW & COL. OF SUBMATRX TO BE DELETED) NO CODE : IF = 1 THE INTERACTION MATRICES ARE COPIED AND EE EXTENDED ( I.E.FIRST CALL FOR A GIVEN PEELING SEQUENCE) IF BETWEEN 1 & NP THE COPIED INTERACTION MATRICES ARE MODIFIED BY DELETION OF THE APPROPRIATE ROW OR COLUMN (THE NORBTH) IF = NP THE INTERACTION MATRICES ARE RESTORED TO THOSE ORIGINALLY OPERATIVE.(I.E. THE LAST CALL OF A SEQUENCE) MM* THE INDEX OF SUBMATRICES CORRESPONDING TO NN MM(I,J) INDEX OF INTERACTION MATRIX BETWEEN ATOM NN(I,J) AND ATOM I ; J.NE.1 . IF J=1 THEN = INDEX OF THE SELF INTERACTION MATRIX OF ATOM I. NN THE INDEX OF NEIGHBOURS NN(I,1) = 1+ NO. OF NEIGHBOURS OF ATOM I NN(I,J), J=2,NN(I,1) LIST OF NEIGHBOURS OF ATOM I ND FIRST DIMENSION OF ARRAYS NN & MM ID SECOND DIMENSION OF ARRAYS NN & MM EE* LIST OF INTERACTION MATRICES NP FIRST 2 DIMENSIONS OF ARRAY EE NE* NO. OF INTERACTION MATRICES SO FAR COMPUTED NED MAX NO. OF INTERACTION MATRICES ALLOWED( LAST DIMENSION OF EE) MEM* STORAGE SPACE TO ENABLE RESTORATION OF THE ORIGINAL MATRIX ./ ADD NAME=OUT SUBROUTINE OUT(IO,IA,IE,IZP,IW,VEC,NED,NE,NAT,MM,NN,ND,NM,EE,NP) INTEGER*2 MM(ND,NM),NN(ND,NM),IZP(NAT),IW(2,NED) DIMENSION VEC(3,NED),EE(NP,NP,NED) SUBROUTINE TO OUTPUT HAMILTONIAN INFORMATION AS DEFINED IN SETUP ARGUMENTS: IO STREAM NUMBER FOR OUTPUT IA INDEX OF FIRST ATOM INFORMATION OUTPUT FROM NN & MM IE INDEX OF FIRST INTERACTION MATRIX OUTPUT (EE) IF =0 THEN NO INTERACTION MATRICES OUTPUT IZP 'TYPE' OF EACH ATOM IW 'TYPE' OF ATOMS OF EACH INTERACTION MATRIX VEC DISPLACEMENT VECTOR FOR EACH INTERACTION MATRIX NED LAST DIMENSION OF ARRAYS IW & VEC NE NUMBER OF INTERACTION MATRICES NAT NUMBER OF ATOMS MM INDEX OF INTERACTION MATRICES AS GENERATED BY MMCAL NN NEIGHBOUR MAP AS GENERATED BY NNCAL ND FIRST DIMENSION OF ARRAYS MM & NN NM SECOND DIMENSION OF ARRAYS MM & NN EE LIST OF INTERACTION MATRICES NP DIMENSION OF INTERACTION MATRICES ( FIRST 2 DIMENSIONS OF EE) ./ ADD NAME=PEEL SUBROUTINE PEEL(CRD,NDIM,NAT,NN,ND,NM,IST,NS,IZP,IZERO,NSH,IW) INTEGER*2 NN(ND,NM),IST(NS),IZP(NAT),IZERO(NAT),IW(NAT) DIMENSION CRD(NDIM,NAT) GIVEN A CLUSTER OF 'SITES' AND ITS 'SHELL' STRUCTURE, RETAINS ONLY THOSE 'SITES' WITHIN A GIVEN NUMBER OF 'SHELLS'.THE ACCEPTED SITES ARE THEN RENUMBERED AS ARE OTHER RELEVANT REFERENCE ARRAYS. ARGUMENTS : (* INDICATES AN OVERWRITTEN ARGUMENT) CRD* LIST OF LATTICE COORDINATES NDIM FIRST DIMENSION OF ARRAY CRD NAT* ON INPUT THE NUMBER OF LATTICE SITES ON OUTPUT THE NUMBER OF RETAINED LATTICE SITES NN* INTEGER*2 NEAREST NEIGHBOUR MAP (AS O/P FROM NNCAL) ND FIRST DIMENSION OF ARRAY NN NM SECOND DIMENSION AF ARRAY NN IST* INTEGER*2 ARRAY LISTING THE 'CENTRAL' SITES NS NUMBER OF SITES LISTED IN IST IZP* INTEGER*2 'TYPE' OF EACH SITE IZERO* INTEGER*2 'SHELL NUMBER' OF EACH SITE NSH NUMBER OF SHELLS TO BE RETAINED IW* INTEGER*2 IW(I) =NEW INDEX OF OLD SITE I ./ ADD NAME=SCAN SUBROUTINE SCAN(NN,ND,NNMX,N0,NAT,NON,SUB) INTEGER*2 NN DIMENSION NN(ND,NNMX),IA(3),NA(3) GENERATES ALL NEIGHBOURS (0TH, 1ST, AND 2ND IF REQUIRED) OF A SUBCLUSTER OF ATOMS (CONSECUTIVELY NUMBERED) OF A GIVEN CLUSTER. INPUT IS THE 'NEAREST NEIGHBOUR' MAP OF THE WHOLE CLUSTER AND OUTPUT IS VIA A USER SUPPLIED SUBROUTINE WHICH IS CALLED FOR EACH POSSIBLE NEIGHBOUR. ARGUMENTS : NN NEAREST NEIGHBOUR MAP. (N.B. INTEGER*2 ARRAY) NN(I,1) CONTAINS 1+ NO. OF NEIGHOURS OF ATOM I NN(I,J),J=2,..,NN(I,1) IS THE LIST OF ATOM NUMBERS OF THE NEIGHBOURS OF ATOM I ND FIRST DIMENSION OF ARRAY NN NNMX SECOND DIMENSION OF ARRAY NN N0 FIRST ATOM OF THE SUBCLUSTER WHOSE NEIGHBOURS ARE TO BE GENERATED NAT LAST ATOM OF THAT SUBCLUSTER NON 'ORDER' OF NEIGHBOURS REQUIRED I.E. 1 IF FIRST NEIGHBOURS ONLY 2 IF FIRST & SECOND NEIGHBOURS SUB NAME OF A USER SUPPLIED SUBROUTINE (DECLARED EXTERNAL IN THE CALLING ROUTINE)TO PROCESS THE INFORMATION GENERATED. ITS ARGUMENTS , WHICH MUST NOT BE MODIFIED, ARE : ......... (IA,NA,NOP) DIMENSION IA(NOP),NA(NOP) NOP CONTAINS THE CODE AS FOLLOWS: 1 FOR THE SELF INTERACTION 2 FOR A 1ST. NEIGHBOUR INTERACTION 3 FOR A 2ND. (NEIGHBOUR OF NEIGHBOUR) INTERACTION IA(NOP) IS THE INDEX OF THE NEIGHBOUR GENERATED I.E. IA(1)=I IA(2)=INDEX OF FIRST NEIGHBOUR OF I (IF NOP>OR= 2) IA(3)=INDEX OF 2ND. NEIGHBOUR OF I (VIA ATOM IA(2)) IF NOP=3 NA(I) IS THE SUBSCRIPT IN THE NEIGHBOUR MAP NN OF THE GENERATED NEIGHBOUR. I.E. NA(1)=1 NA(2)=J WHERE IA(2)=NN(I,J) (IF NOP>OR= 2) NA(3)=K WHERE IA(3)=NN(J,K) (IF NOP=3) ./ ADD NAME=SELFD SUBROUTINE SELFD(I,J,R2,IOVPAR,EM,NE) DIMENSION PARM(13),EM(NE,5) RETURNS A DIAGONAL 'SELF ENERGY MATRIX' APPROPRIATE FOR D ELECTRONS WRITTEN BY ESTHER HAINES ARGUMENTS : (* INDICATES OVERWRITTEN BY THE ROUTINE) I 'TYPE' OF ATOM J DUMMY ARGUMENT (=I) R2 DUMMY ARGUMENT (=1.0E-4) IOVPAR NAME OF AN INTEGER FUNCTION RETURNING THE SELF ENERGY ARGUMENTS : I TYPE OF ATOM J TYPE OF ATOM R2 DUMMY ARGUMENT (=1.0E-4) PARM PARM(11)= SELF ENERGY EM* DIAGONAL MATRIX (EM(I,I)=PARM(11) ) NE FIRST DIMENSION OF ARRAY EM ./ ADD NAME=SELFP SUBROUTINE SELFP(I,J,R2,IOVPAR,EM,NE) DIMENSION PARM(13),EM(NE,3) RETURNS A DIAGONAL 'SELF ENERGY MATRIX' APPROPRIATE FOR P ELECTRONS WRITTEN BY ESTHER HAINES ARGUMENTS : (* INDICATES OVERWRITTEN BY THE ROUTINE) I 'TYPE' OF ATOM J DUMMY ARGUMENT (=I) R2 DUMMY ARGUMENT (=1.0E-4) IOVPAR NAME OF AN INTEGER FUNCTION RETURNING THE SELF ENERGY ARGUMENTS : I TYPE OF ATOM J TYPE OF ATOM R2 DUMMY ARGUMENT (=1.0E-4) PARM PARM(12)= SELF ENERGY EM* DIAGONAL MATRIX (EM(I,I)=PARM(12) ) NE FIRST DIMENSION OF ARRAY EM ./ ADD NAME=SELFS SUBROUTINE SELFS(I,J,R2,IOVPAR,EM,NE) DIMENSION PARM(13),EM(NE,1) RETURNS A DIAGONAL 'SELF ENERGY MATRIX' APPROPRIATE FOR S ELECTRONS WRITTEN BY ESTHER HAINES ARGUMENTS : (* INDICATES OVERWRITTEN BY THE ROUTINE) I 'TYPE' OF ATOM J DUMMY ARGUMENT (=I) R2 DUMMY ARGUMENT (=1.0E-4) IOVPAR NAME OF AN INTEGER FUNCTION RETURNING THE SELF ENERGY ARGUMENTS : I TYPE OF ATOM J TYPE OF ATOM R2 DUMMY ARGUMENT (=1.0E-4) PARM PARM(13)= SELF ENERGY EM* DIAGONAL MATRIX (EM(1,1)=PARM(13) ) NE FIRST DIMENSION OF ARRAY EM ./ ADD NAME=SETUP SUBROUTINE SETUP(CRD,ND,NAT,EV,NTYPE,IZP,MM,NN,NND,NM,HCAL,NGBR 1,IOVPAR,EE,NP,NED,NE,VEC,IW) INTEGER*2 MM(NND,NM),NN(NND,NM),IZP(NAT),IW(2,NED) LOGICAL EV DIMENSION CRD(ND,NAT),VEC(3,NED),EE(NP,NP,NED) EXTERNAL NGBR,IOVPAR,EV ASSEMBLES THE HAMILTONIAN MATRIX FROM THE USER SUPPLIED ROUTINES EV, HCAL, NGBR AND IOVPAR, AND THE LIBRARY ROUTINES NNCAL & MMCAL . ARGUMENTS OF SETUP : (* INDICATES OVERWRITTEN BY THE ROUTINE) CRD LATTICE COORDINATES ND FIRST DIMENSION OF CRD NAT NO.OF ATOMS IN THE CLUSTER EV LOGICAL FUNCTION OF 2 ARGUMENTS, BOTH REAL ARRAYS OF LENGTH 3 RETURNING THE VALUE .TRUE. IF THE ARRAYS ARE EQUIVALENT AND .FALSE. IF NOT. NTYPE NO. OF DIFFERENT 'TYPES' OF ATOMS IZP 'TYPE' OF EACH ATOM MM* IS THE INTERACTION MAP GENERATED BY MMCAL NN* IS THE NEIGHBOUR MAP GENERATED BY NNCAL NND FIRST DIMENSION OF ARRAYS MM & NN NM* MAX NO. OF ATOMS CONNECTED BY INTERACTIONS. ON OUTPUT CONTAINS ACTUAL MAX NO. GENERATED HCAL NAME OF A SUBROUTINE TO CALCULATE THE INTERACTION BETWEEN TWO ATOMS. ARGUMENTS ARE V VECTOR POSITION(I) - POSTITION(J) II TYPE AT I JJ TYPE AT J E* OUTPUT INTERACTION MATRIX H OPERATING ON PSI(J) EFFECT AT I IOVPAR NAME OF FUNCTION SUPPLING INFORMATION TO HCAL NGBR NAME OF A FUNCTION TO SUPPLY INTERACTION INFORMATION TO NNCAL ARGUMENTS ARE : II 'TYPE' OF ATOM I JJ 'TYPE' OF ATOM J R2 SQUARE OF DISTANCE FROM I TO J DD DUMMY ARGUMENT NGBR TAKES THE VALUE 1 IF I & J ARE NEIGHBOURS AND 0 OTHERWISE IOVPAR NAME OF A FUNCTION TO SUPPLY OVERLAP PARAMETERS TO HCAL ARGUMENTS ARE II 'TYPE' OF ATOM I JJ 'TYPE' OF ATOM J R2 SQUARE OF THE DISTANCE FROM I TO J DD* OVERLAP PARAMETERS AS REQUIRED BY HCAL THE NOTATION USED IS AS FOLLOWS: DD(1) DD SIGMA DD(2) DD PI DD(3) DD DELTA DD(4) PD SIGMA DD(5) PD PI DD(6) PP SIGMA DD(7) PP PI DD(8) SD SIGMA DD(9) SP SIGMA DD(10) SS SIGMA DD(11) D SELF ENERGY DD(12) P SELF ENERGY DD(13) S SELF ENERGY EE LIST OF INTERACTION MATRICES NP FIRST 2 DIMENSIONS OF ARRAY EE NED LAST DIMENSION OF ARRAYS EE,IW,VEC NE* NO. OF DISTINCT DISPLACEMENT VECTORS (MATRICES) FOUND VEC* LIST OF DISTINCT DISPLACEMENT VECTORS FOUND (POSN. J - POSN.I) IW* LIST OF ATOM TYPES AT THE ENDS OF THE VECTORS IN VEC IW(1,.) IS TYPE OF I IW(2,.) IS TYPE OF J ./ ADD NAME=SKDD SUBROUTINE SKDD(X,X2,I,J,R2,IOVPAR,EM,NE) DIMENSION X(6),X2(6),PARM(13),EM(NE,5),E(15,3),DM(15),DD(3) THIS SUBROUTINE CALCULATES THE SLATER-KOSTER OVERLAP MATRICES FOR DD ELECTRONS. WRITTEN BY ESTHER HAINES ARGUMENTS (* INDICATES AN OVERWRITTEN ARGUMENT) X X(1)- X(3) ARE THE DIRECTION COSINES ATOM JJ - ATOM II AND X(4) - X(6) ARE A SECOND COPY OF THEM X2 THE SQUARES OF X I TYPE OF ATOM II J TYPE OF ATOM JJ R2 THE SQUARE OF THE DISTANCE FROM ATOM II TO JJ IOVPAR THE NAME OF AN INTEGER FUNCTION RETURNING THE OVERLAP PARAMETERS . ITS ARGUMENTS ARE AS FOLLOWS: I TYPE OF ATOM II J TYPE OF ATOM JJ R2 THE SQUARE OF THE DISTANCE FROM II TO JJ DD* THE OVERLAP PARAMETERS : DD(1) DD SIGMA DD(2) DD PI DD(3) DD DELTA EM* THE COMPUTED 5X5 OVERLAP MATRIX (OPERATING ON PSI(JJ)) NE FIRST DIMENSION OF ARRAY EM ./ ADD NAME=SKPD SUBROUTINE SKPD(X,X2,I,J,R2,IOVPAR,EM,EMT,NE) DIMENSION X(6),X2(6),PARM(13),EM(NE,5),EMT(NE,5),EPD(13,2),DM(15) THIS SUBROUTINE CALCULATES THE SLATER-KOSTER OVERLAP MATRICES FOR PD ELECTRONS. WRITTEN BY ESTHER HAINES ARGUMENTS (* INDICATES AN OVERWRITTEN ARGUMENT) X X(1)- X(3) ARE THE DIRECTION COSINES ATOM JJ - ATOM II AND X(4) - X(6) ARE A SECOND COPY OF THEM X2 THE SQUARES OF X I TYPE OF ATOM II J TYPE OF ATOM JJ R2 THE SQUARE OF THE DISTANCE FROM ATOM II TO JJ IOVPAR THE NAME OF AN INTEGER FUNCTION RETURNING THE OVERLAP PARAMETERS . ITS ARGUMENTS ARE AS FOLLOWS: I TYPE OF ATOM II J TYPE OF ATOM JJ R2 THE SQUARE OF THE DISTANCE FROM II TO JJ DD* THE OVERLAP PARAMETERS : DD(4) PD SIGMA DD(5) PD PI EM* THE 3X5 COMPUTED MATRIX (ACTING ON THE D STATES OF PSI(JJ)) EMT* THE 5X3 COMPUTED MATRIX (ACTING ON THE P STATES OF PSI(JJ)) NE FIRST DIMENSION OF ARRAY EM AND EMT ./ ADD NAME=SKPP SUBROUTINE SKPP(X,X2,I,J,R2,IOVPAR,EM,NE) DIMENSION X(6),DM(6),X2(6),PARM(13),EM(NE,3),EPP(6) THIS SUBROUTINE CALCULATES THE SLATER-KOSTER OVERLAP MATRICES FOR PP ELECTRONS. WRITTEN BY ESTHER HAINES ARGUMENTS (* INDICATES AN OVERWRITTEN ARGUMENT) X X(1)- X(3) ARE THE DIRECTION COSINES ATOM JJ - ATOM II AND X(4) - X(6) ARE A SECOND COPY OF THEM X2 THE SQUARES OF X I TYPE OF ATOM II J TYPE OF ATOM JJ R2 THE SQUARE OF THE DISTANCE FROM ATOM II TO JJ IOVPAR THE NAME OF AN INTEGER FUNCTION RETURNING THE OVERLAP PARAMETERS . ITS ARGUMENTS ARE AS FOLLOWS: I TYPE OF ATOM II J TYPE OF ATOM JJ R2 THE SQUARE OF THE DISTANCE FROM II TO JJ DD* THE OVERLAP PARAMETERS : DD(6) PP SIGMA DD(7) PP PI EM* THE COMPUTED 3X3 OVERLAP MATRIX (OPERATING ON PSI(JJ)) NE FIRST DIMENSION OF ARRAY EM ./ ADD NAME=SKSD SUBROUTINE SKSD(X,X2,I,J,R2,IOVPAR,EM,EMT,NE) DIMENSION X(6),X2(6),PARM(13),EM(NE,5),EMT(NE,5),ES(5) THIS SUBROUTINE CALCULATES THE SLATER-KOSTER OVERLAP MATRICES FOR SD ELECTRONS. WRITTEN BY ESTHER HAINES ARGUMENTS (* INDICATES AN OVERWRITTEN ARGUMENT) X X(1)- X(3) ARE THE DIRECTION COSINES ATOM JJ - ATOM II AND X(4) - X(6) ARE A SECOND COPY OF THEM X2 THE SQUARES OF X I TYPE OF ATOM II J TYPE OF ATOM JJ R2 THE SQUARE OF THE DISTANCE FROM ATOM II TO JJ IOVPAR THE NAME OF AN INTEGER FUNCTION RETURNING THE OVERLAP PARAMETERS . ITS ARGUMENTS ARE AS FOLLOWS: I TYPE OF ATOM II J TYPE OF ATOM JJ R2 THE SQUARE OF THE DISTANCE FROM II TO JJ DD* THE OVERLAP PARAMETERS : DD(8) SD SIGMA EM* THE 1X5 COMPUTED MATRIX (ACTING ON THE D STATES OF PSI(JJ)) EMT* THE 5X1 COMPUTED MATRIX (ACTING ON THE S STATES OF PSI(JJ)) NE FIRST DIMENSION OF ARRAY EM AND EMT ./ ADD NAME=SKSP SUBROUTINE SKSP(X,X2,I,J,R2,IOVPAR,EM,EMT,NE) DIMENSION X(6),X2(6),PARM(13),EM(NE,3),EMT(NE,3) THIS SUBROUTINE CALCULATES THE SLATER-KOSTER OVERLAP MATRICES FOR SP ELECTRONS. WRITTEN BY ESTHER HAINES ARGUMENTS (* INDICATES AN OVERWRITTEN ARGUMENT) X X(1)- X(3) ARE THE DIRECTION COSINES ATOM JJ - ATOM II AND X(4) - X(6) ARE A SECOND COPY OF THEM X2 THE SQUARES OF X I TYPE OF ATOM II J TYPE OF ATOM JJ R2 THE SQUARE OF THE DISTANCE FROM ATOM II TO JJ IOVPAR THE NAME OF AN INTEGER FUNCTION RETURNING THE OVERLAP PARAMETERS . ITS ARGUMENTS ARE AS FOLLOWS: I TYPE OF ATOM II J TYPE OF ATOM JJ R2 THE SQUARE OF THE DISTANCE FROM II TO JJ DD* THE OVERLAP PARAMETERS : DD(9) SP SIGMA EM* THE 1X3 COMPUTED MATRIX (ACTING ON THE P STATES OF PSI(JJ)) EMT* THE 3X1 COMPUTED MATRIX (ACTING ON THE S STATES OF PSI(JJ)) NE FIRST DIMENSION OF ARRAY EM AND EMT ./ ADD NAME=SKSS SUBROUTINE SKSS(X,X2,I,J,R2,IOVPAR,EM,NE) DIMENSION X(6),X2(6),PARM(13),EM(1,1) THIS SUBROUTINE CALCULATES THE SLATER-KOSTER OVERLAP MATRICES FOR SS ELECTRONS. WRITTEN BY ESTHER HAINES ARGUMENTS (* INDICATES AN OVERWRITTEN ARGUMENT) X X(1)- X(3) ARE THE DIRECTION COSINES ATOM JJ - ATOM II AND X(4) - X(6) ARE A SECOND COPY OF THEM.(NOT USED) X2 THE SQUARES OF X (NOT USED) I TYPE OF ATOM II J TYPE OF ATOM JJ R2 THE SQUARE OF THE DISTANCE FROM ATOM II TO JJ IOVPAR THE NAME OF AN INTEGER FUNCTION RETURNING THE OVERLAP PARAMETERS . ITS ARGUMENTS ARE AS FOLLOWS: I TYPE OF ATOM II J TYPE OF ATOM JJ R2 THE SQUARE OF THE DISTANCE FROM II TO JJ DD* THE OVERLAP PARAMETERS : DD(10) SS SIGMA EM* THE COMPUTED 1X1 OVERLAP MATRIX (OPERATING ON PSI(JJ)) NE FIRST DIMENSION OF ARRAY EM ./ ADD NAME=SLKODE SUBROUTINE SLKODE(DUM,I,J,EM,IBONDS) DIMENSION DUM(3),EM(5,5) CALCULATES THE 5 X 5 2 CENTRE INTERACTION MATRICES BETWEEN THE D-ELECTRONS ON TWO SITES A GIVEN VECTOR APART.THE PRESCRIPTION USED IS THAT OF SLATER & KOSTER IN PHYS.REV 94 VOL 6 P.1498 ET.SEQ.. THE ORBITALS ARE ORDERED XY,YZ,ZX,X*X-Y*Y,3Z*Z-R*R IF DUM HAS SQUARED MODULUS < 1.0E-4 THEN A DIAGONAL SELF INTERACTION MATRIX IS RETURNED THIS ROUTINE CALLS SKDD AND SELFD ARGUMENTS (* INDICATES OVERWRITTEN BY THE ROUTINE) DUM THE VECTOR B-A I THE 'TYPE' OF ATOM A J THE 'TYPE' OF ATOM B EM* THE 5 X 5 MATRIX SUCH THAT (HY(B))(A) = EM Y(B) IBONDS THE NAME OF AN INTEGER FUNCTION SUPPLIED BY THE CALLING ROUTINE AND DECLARED EXTERNAL WITH THE FOLLOWING ARGUMENTS : I 'TYPE' OF ATOM A J 'TYPE' OF ATOM B R2 SQUARE OF DISTANCE FROM A TO B DD* DD(SIGMA,PI,DELTA) FOR THE ABOVE ARGUMENTS DD(11) CONTAINS THE SELF ENERGY IF R2 < 1.0E-4 ./ ADD NAME=STORE SUBROUTINE STORE(ISTR,CRD,NC,IZP,NAT,MM,NN,NND,NM,VEC,IW,NE,EE,NP) INTEGER*2 IZP(NAT),MM(NND,NM),NN(NND,NM),IW(2,NE) DIMENSION CRD(NC,NAT),VEC(3,NE),EE(NP,NP,NE) STORES OR READS THE HAMILTONIAN MATRIX INFORMATION AS GENERATED BY SETUP AND USED IN THE RECURSION METHOD (I.E. RECAL AND CRECAL) NO OVERWRITING OCCURS WHEN OUTPUT IS INTENDED (I.E. ISTR >0) ARGUMENTS (* INDICATES AN OVERWRITTEN ARGUMENT) ISTR IF >0 THE STREAM NUMBER TO WHICH OUTPUT IS SENT IF =0 NO OUTPUT IS MADE IF <0 -THE STREAM NUMBER FROM WHICH INPUT IS READ CRD* THE ATOM COORDINATES NC FIRST DIMENSION OF ARRAY CRD IZP* 'TYPE' OF EACH ATOM .(INTEGER*2) NAT* NUMBER OF ATOMS MM* HAMILTONIAN MAP AS GENERATED BY SETUP NN* NEIGHBOUR MAP AS GENERATED BY SETUP NND FIRST DIMENSION OF ARRAYS NN & MM NM* MAXIMUM CONNECTIVITY (1+MAX. NO. OF NEIGHBOURS) VEC* LIST OF INTERACTION VECTORS IW* LIST OF INTERACTION 'TYPES' NE* NUMBER OF INTERACTION MATRICES. IF 0 THEN NO MATRICES OUTPUT EE* LIST OF INTERACTION MATRICES NP FIRST 2 DIMENSIONS OF EE (DIMENSION OF INTERACTION MATRICES) ./ ENDUP ./ ADD NAME=BNDCRD SUBROUTINE BNDCRD(ET,IC,WTT,NET,NB,AA,RNG) DIMENSION ET(NET),IC(NET),WTT(NET),AA(NB),RNG(NB) DIMENSION NC(9) DO 31 I=1,9 31 NC(I)=0 WMX=0.0 DO 32 I=1,NET IF(IC(I).EQ.4)WMX=AMAX1(WMX,WTT(I)) ICC=IC(I)+5 32 NC(ICC)=NC(ICC)+1 NL=NC(6)+NC(8) NR=NC(2)+NC(4) IF(NL.NE.NR)GOTO 51 IF(NL.GT.NB)GOTO 52 IF(NL.EQ.NB)GOTO 37 NM=NB-NL IF(NM.LT.NC(9))GOTO 34 DO 33 I=1,NET IF(IC(I).EQ.4)IC(I)=0 33 CONTINUE GOTO 37 34 DO 36 I=1,NM WM=WMX DO 35 J=1,NET IF(IC(J).NE.4)GOTO 35 IF(WTT(J).GT.WM)GOTO 35 WM=WTT(J) JJ=J 35 CONTINUE 36 IC(JJ)=0 37 NB=0 DO 41 II=1,NET ICC=IC(II)+5 GOTO(41,39,41,39,40,38,41,38,41),ICC 38 NB=NB+1 AA(NB)=ET(II) GOTO 41 39 RNG(NB)=ET(II)-AA(NB) GOTO 41 40 RNG(NB)=(ET(II)+ET(II-1))*0.5-AA(NB) NB=NB+1 AA(NB)=(ET(II+1)+ET(II))*0.5 41 CONTINUE RETURN 51 NB=0 RETURN 52 NB=-NL RETURN END ./ ADD NAME=BNDEST SUBROUTINE BNDEST(A,B2,LL,EPS,AA,RNG,NB,EV,FEV,IC,NET,DE,WK,NW) DIMENSION A(LL),B2(LL),AA(NB),RNG(NB),EV(NET),FEV(NET),IC(NET) 1 ,WK(NW,4),P(2,3) NPTS=NW ALP=1.0 LM1=LL-1 EMN=AMIN1(A(1)-B2(2),A(LM1)-1.0) EMX=AMAX1(A(1)+B2(2),A(LM1)+1.0) NID=LM1-1 IF(NID.LT.2)GOTO 2 DO 1 I=2,NID EMN=AMIN1(EMN,A(I)-1.0-B2(I+1)) 1 EMX=AMAX1(EMX,A(I)+1.0+B2(I+1)) 2 DE=EMX-EMN EMN=EMN-DE/FLOAT(LL) EMX=EMX+DE/FLOAT(LL) DE=(EMX-EMN)/FLOAT(NPTS-1) E0=EMN DO 3 NE=1,NPTS WK(NE,1)=E0+FLOAT(NE-1)*DE 3 WK(NE,2)=RECWT(WK(NE,1),A,B2,LL,EPS,1,P,1) THU=ALP*B2(1)/FLOAT(LL) THL=B2(1)/(6.0*FLOAT(LL)**3) DO 4 I=1,5 NET=NW CALL TABAN(WK,WK(1,2),NPTS,THU,THL,WK(1,3),IC,WK(1,4),NET) IF(NET.GT.0)GOTO 5 4 THU=(THU+THL)*0.5 RETURN 5 CALL BNDCRD(WK(1,3),IC,WK(1,4),NET,NB,AA,RNG) IF(NB.LE.0)RETURN IB=0 J=0 DO 7 I=1,NET IF(IABS(IC(I)).GT.1)GOTO 7 J=J+1 IF(IC(I).EQ.0)GOTO 6 IC(J)=-IC(I) EV(J)=WK(I,3) FEV(J)=RECWT(EV(J),A,B2,LL,EPS,1,P,1) GOTO 7 6 IC(J)=0 E=WK(I,3) IT=30 CALL WTMIN(E-DE,E+DE,A,B2,LL,EPS,EPS*10.0,IT,EV(J),FEV(J)) 7 CONTINUE NET=J RETURN END ./ ADD NAME=BNDREF SUBROUTINE BNDREF(DEL,AM,BM2,LL,EPS,AA,RNG,NB,BWK,NBD,IC,NET) DIMENSION AM(LL),BM2(LL),AA(NB),RNG(NB),BWK(NBD,3),IC(NBD) 1 ,P(2,3) J=0 DO 33 I=1,NET J=J+1 IF(IC(I).EQ.0)GOTO 32 VAL=RECWT(BWK(I,1),AM,BM2,LL,EPS,1,P,1) BWK(J,3)=(BWK(I,2)-VAL)*FLOAT(IC(I)) GOTO 33 32 IT=30 ERR=ABS(BWK(J,3)) AE=BWK(I,1)-ERR BE=BWK(I,1)+ERR ERR=ERR*0.01 CALL WTMIN(AE,BE,AM,BM2,LL,EPS,ERR,IT,EM,FEM) BWK(J,3)=BWK(I,1)-EM J=J+1 BWK(J,3)=FEM-BWK(I,2) 33 CONTINUE IB=0 J=0 DO 44 I=1,NET J=J+1 IF(IC(I))41,43,42 41 IB=IB+1 DA=SIGN(DEL,BWK(J,3)) AA(IB)=AA(IB)+DA GOTO 44 42 RNG(IB)=RNG(IB)+SIGN(DEL,BWK(J,3))-DA GOTO 44 43 D1=SIGN(DEL,BWK(J,3)) J=J+1 D2=SIGN(DEL,BWK(J,3))*0.5 RNG(IB)=RNG(IB)+D1-D2-DA IB=IB+1 DA=D1+D2 AA(IB)=AA(IB)+DA 44 CONTINUE RETURN END ./ ADD NAME=BNDWT SUBROUTINE BNDWT(AA,RNG,WB,NB,A,B2,LL,EPS,WK,NW,IWK) DIMENSION AA(NB),RNG(NB),WB(NB),A(LL),B2(LL),WK(NW,4),IWK(NW) IF(NB.LT.2)GOTO 23 DO 22 I=2,NB E=AA(I-1)+RNG(I-1) ID=0 WB(I-1)=0.0 DO 22 K=1,2 DUM=DENQD(E,E,A,B2,LL,0.5,EPS,WK,NW,NQ,NE,IWK) IF(NE.LE.0)GOTO 35 JL=NE-ID IF(JL.LT.1)GOTO 21 DO 9 J=1,JL 9 WB(I-1)=WB(I-1)+WK(J,2) 21 ID=1 E=AA(I) 22 CONTINUE 23 WB(NB)=B2(1) IF(NB.LT.2)RETURN DO 24 II=2,NB I=NB-II+2 WB(I-1)=WB(I-1)*0.5 24 WB(I)=WB(I)-WB(I-1) RETURN 35 NB=0 RETURN END ./ ADD NAME=CFGEN SUBROUTINE CFGEN(X,W,N,EPS,A,B2,LM1,WK) DIMENSION X(N),W(N),A(LM1),B2(LM1),WK(N,2) DOUBLE PRECISION SA,SB,DSQRT,DUM SA=0.0D0 SB=0.0D0 DO 1 I=1,N SB=SB+W(I) SA=SA+X(I)*W(I) WK(I,1)=0.0 1 WK(I,2)=1.0 B2(1)=SB A(1)=SA/SB IF(LM1.LT.2)RETURN DO 3 L=2,LM1 ANRM=DSQRT(SB) SA=0.0D0 SB=0.0D0 DO 2 I=1,N DUM=WK(I,2) WK(I,2)=((X(I)-A(L-1))*WK(I,2)-B2(L-1)*WK(I,1))/ANRM WK(I,1)=DUM/ANRM DUM=WK(I,2)*WK(I,2)*W(I) SB=SB+DUM 2 SA=SA+X(I)*DUM B2(L)=SB IF(SB.LT.EPS)GOTO 4 3 A(L)=SA/SB L=LM1+1 4 LM1=L-1 RETURN END ./ ADD NAME=CFGPGN SUBROUTINE CFGPGN(AA,RNG,WB,NBP1,IC,EPS,A,B2,LM1,WK,NW) DIMENSION AA(NBP1),RNG(NBP1),WB(NBP1),A(LM1),B2(LM1) DIMENSION WK(NW,2,NBP1) DOUBLE PRECISION SA,SB,SAJ,SBJ,DSQRT,DUM PI=3.14159265359 NB=NBP1-1 IF(IC-2)10,20,1 10 DO 11 I=1,LM1 TH=FLOAT(LM1-I+1)*PI/FLOAT(LM1+1) WK(I,1,1)=COS(TH) WK(I,2,1)=2.0*(1.0-WK(I,1,1)*WK(I,1,1))/FLOAT(LM1+1) 11 WK(I,1,1)=(WK(I,1,1)+1.0)*0.5 GOTO 1 20 DO 21 L=1,LM1 A(L)=0.5 AL2=FLOAT((L-1)*(L-1)) 21 B2(L)=AL2/(16.0*AL2-4.0) B2(1)=1.0 CALL RECQD(A,B2,LM1,WK,WK(1,2,1),LO,EPS,WK(1,1,2),WK(1,2,2)) IF(LO.EQ.LM1)GOTO 1 IC=0 RETURN 1 SA=0.0D0 SB=0.0D0 DO 3 J=1,NB SAJ=0.0D0 JP=J+1 DO 2 I=1,LM1 SAJ=SAJ+(WK(I,1,1)*RNG(J)+AA(J))*WK(I,2,1) WK(I,1,JP)=0.0 2 WK(I,2,JP)=1.0 SB=SB+WB(J) 3 SA=SA+SAJ*WB(J) B2(1)=SB A(1)=SA/SB IF(LM1.LE.1)RETURN DO 6 L=2,LM1 ANRM=DSQRT(SB) SA=0.0D0 SB=0.0D0 DO 5 J=1,NB JP=J+1 SAJ=0.0D0 SBJ=0.0D0 DO 4 I=1,LM1 XX=WK(I,1,1)*RNG(J)+AA(J) DUM=WK(I,2,JP) WK(I,2,JP)=((XX-A(L-1))*WK(I,2,JP)-B2(L-1)*WK(I,1,JP))/ANRM WK(I,1,JP)=DUM/ANRM DUM=WK(I,2,JP)*WK(I,2,JP)*WK(I,2,1) SBJ=SBJ+DUM 4 SAJ=SAJ+XX*DUM SA=SA+SAJ*WB(J) 5 SB=SB+SBJ*WB(J) B2(L)=SB IF(SB.LT.EPS)GOTO 7 6 A(L)=SA/SB L=LM1+1 7 LM1=L-1 RETURN END ./ ADD NAME=CRECAL SUBROUTINE CRECAL(HOP,PSI,PMN,M,A,B2,N) COMPLEX PSI,PMN,CONJG,DUM DIMENSION PSI(M),PMN(M),A(N),B2(N) DOUBLE PRECISION SUM,DSQRT SUM=B2(1) NM1=N-1 DO 2 J=1,NM1 B2(J)=SUM S=1.0/DSQRT(SUM) CALL HOP(PSI,PMN,A(J)) A(J)=A(J)/B2(J) SUM=0.0D0 DO 1 I=1,M DUM=PSI(I) PSI(I)=(PMN(I)-A(J)*PSI(I))*S PMN(I)=DUM 1 SUM=SUM+REAL(CONJG(PSI(I))*PSI(I)) S=S*SUM DO 2 I=1,M 2 PMN(I)=-PMN(I)*S B2(N)=SUM RETURN END ./ ADD NAME=DENCRQ COMPLEX FUNCTION DENCRQ(E,A,B2,LL,AA,RNG,WB,NB,AM,BM2) DIMENSION A(LL),B2(LL),AA(NB),RNG(NB),WB(NB),AM(LL),BM2(LL) DIMENSION P(2),Q(2) DENCRQ=(0.0,0.0) DO 2 I=1,NB DISC=(E-AA(I))*(E-AA(I)-RNG(I)) WT=8.0*WB(I)/(RNG(I)*RNG(I)) IF (DISC.LT.0.0)GOTO 1 DISC=SQRT(DISC) IF(E.LT.AA(I))DISC=-DISC DENCRQ=DENCRQ+WT*(E-AA(I)-RNG(I)*0.5-DISC) GOTO 2 1 DISC=SQRT(-DISC) DENCRQ=DENCRQ+CMPLX(WT*(E-AA(I)-RNG(I)*0.5),-WT*DISC) 2 CONTINUE CALL PLYVAL(E,AM,BM2,LL-1,P,Q) DENCRQ=(Q(2)-DENCRQ*P(2))/(BM2(LL)*(Q(1)-DENCRQ*P(1))) CALL PLYVAL(E,A,B2,LL-1,P,Q) DENCRQ=(Q(2)-B2(LL)*DENCRQ*Q(1))/(P(2)-B2(LL)*DENCRQ*P(1)) RETURN END ./ ADD NAME=DENCRS FUNCTION DENCRS(E,A,B2,LL,AA,RNG,WB,NB,AM,BM2) DIMENSION A(LL),B2(LL),AA(NB),RNG(NB),WB(NB),AM(LL),BM2(LL) DIMENSION P(2),Q(2) PI=3.14159265359 AIF=0.0 RF=0.0 DENCRS=0.0 DO 2 I=1,NB DISC=(E-AA(I))*(E-AA(I)-RNG(I)) WT=8.0*WB(I)/(RNG(I)*RNG(I)) IF (DISC.LT.0.0)GOTO 1 DISC=SQRT(DISC) IF(E.LT.AA(I))DISC=-DISC RF=RF+WT*(E-AA(I)-RNG(I)*0.5-DISC) GOTO 2 1 DISC=SQRT(-DISC) RF=RF+WT*(E-AA(I)-RNG(I)*0.5) AIF=AIF-WT*DISC 2 CONTINUE IF(AIF.GE.0.0)RETURN CALL PLYVAL(E,AM,BM2,LL-1,P,Q) RD=Q(1)-RF*P(1) AID=-AIF*P(1) DENOM=BM2(LL)*(RD*RD+AID*AID) RT=((Q(2)-RF*P(2))*RD-AIF*P(2)*AID)/DENOM AIT=(-RD*AIF*P(2)-(Q(2)-RF*P(2))*AID)/DENOM CALL PLYVAL(E,A,B2,LL-1,P,Q) T1=P(2)-RT*P(1)*B2(LL) T2=AIT*P(1)*B2(LL) DENCRS=-AIT*B2(LL)*(Q(2)*P(1)-Q(1)*P(2))/((T1*T1+T2*T2)*PI) RETURN END ./ ADD NAME=DENINT FUNCTION DENINT(E,A,B2,NA,NP,LL,ALP,EPS,WK,IWK,ICODE) DIMENSION A(NA,NP),B2(NA,NP),WK(LL,4),IWK(LL) SUM=0.0 DO 2 II=1,NP DUM=DENQD(E,E,A(1,II),B2(1,II),LL,ALP,EPS,WK,LL,NQ,NE,IWK) IF(NE.LT.1) GOTO 3 SUM=SUM+(ALP-1.0)*WK(NE,2) DO 1 I=1,NE 1 SUM=SUM+WK(I,2) 2 CONTINUE ICODE=0 DENINT=SUM RETURN 3 ICODE=-1 RETURN END ./ ADD NAME=DENQD FUNCTION DENQD(E,EMX,A,B2,LL,ALP,EPS,TB,NT,NQ,NE,IWK) DIMENSION A(LL),B2(LL),TB(NT,4),IWK(LL) DOUBLE PRECISION P(2),PP(2),PPP(2),FI(2),FIP(2),XW,SC,DABS DOUBLE PRECISION P1,P0,DUM,AN N=-1 DUM=RECWT(E,A,B2,LL,EPS,N,TB,1) NN=LL+N NR=0 CALL RECRTS(A,B2,NN,EPS,EMX+EPS,NR,TB,IWK,TB(1,2),TB(1,3)) IC=1 NE=0 DO 5 NQ=1,NR IF(IWK(NQ).NE.1)RETURN IF(TB(NQ,1).GT.EMX+EPS)GOTO 6 XV=TB(NQ,1) P(1)=1.0D0 P(2)=XV-A(1) PP(1)=0.0D0 PP(2)=1.0D0 PPP(1)=0.0D0 PPP(2)=0.0D0 FI(1)=0.0D0 FI(2)=B2(1) FIP(1)=0.0D0 FIP(2)=0.0D0 I=2 DO 4 J=2,NN XW=XV-A(J) I1=3-I P(I1)=XW*P(I)-B2(J)*P(I1) PP(I1)=XW*PP(I)-B2(J)*PP(I1)+P(I) PPP(I1)=XW*PPP(I)-B2(J)*PPP(I1)+2.0*PP(I) FI(I1)=XW*FI(I)-B2(J)*FI(I1) FIP(I1)=XW*FIP(I)-B2(J)*FIP(I1)+FI(I) SC=DABS(P(I1))+DABS(PP(I1))+DABS(PPP(I1))+DABS(FI(I1)) 1 +DABS(FIP(I1)) DO 3 K=1,2 P(K)=P(K)/SC PP(K)=PP(K)/SC PPP(K)=PPP(K)/SC FI(K)=FI(K)/SC 3 FIP(K)=FIP(K)/SC 4 I=I1 I2=3-I TB(NQ,2)=FI(I)/PP(I) P2=PP(I)*PP(I) TB(NQ,4)=P(I2)/PP(I) 5 TB(NQ,3)=(FI(I)*PP(I2)-PP(I)*FI(I2)+TB(NQ,4) 1 *(PP(I)*FIP(I)-FI(I)*PPP(I)))/P2 NQ=NR+1 6 NQ=NQ-1 DEV=ABS(TB(1,1)-E) DEN=0.0 I=2 IF(NQ.EQ.1)GOTO 8 DO 7 I=2,NQ WD=ABS(TB(I,1)-E) IF(DEV.LT.WD)GOTO 8 DEV=WD 7 DEN=DEN+TB(I-1,3) I=NQ+1 8 NE=I-1 DENQD=(DEN+TB(NE,3)*ALP)/TB(NE,4) IF(DEV.GT.10.0*EPS)NE=-NE RETURN END ./ ADD NAME=DENSQ FUNCTION DENSQ(E,A,B2,LL,EI) DIMENSION A(LL),B2(LL),EI(2),P(2),Q(2) PI=3.14159265359 DENSQ=0.0 DISC=(E-EI(1))*(E-EI(2)) IF(DISC.GT.0.0)RETURN BI2=(EI(2)-EI(1))*0.25 BI2=BI2*BI2*2.0 AIT=-SQRT(-DISC)/BI2 RT=(E-(EI(1)+EI(2))*0.5)/BI2 CALL PLYVAL(E,A,B2,LL-1,P,Q) T1=P(2)-RT*P(1)*B2(LL) T2=AIT*P(1)*B2(LL) DENSQ=-AIT*B2(LL)*(Q(2)*P(1)-Q(1)*P(2))/((T1*T1+T2*T2)*PI) RETURN END ./ ADD NAME=EDIFF FUNCTION EDIFF(EF,A,B2,LL,EPS,WK,MU) DIMENSION A(LL),B2(LL),WK(LL,4),MU(LL,2) DIMENSION P(2),PP(2),PPP(2),FI(2),FIP(2) DOUBLE PRECISION P,PP,PPP,FI,FIP,XW,SC,DABS DOUBLE PRECISION P1,P0,DUM,AN NN=LL-1 P1=1.0D0 P0=0.0D0 DO 1 I=2,NN DUM=(EF-A(I))*P1-B2(I)*P0 AN=DABS(DUM)+DABS(P1) P0=P1/AN 1 P1=DUM/AN IF(DABS(P1).LT.EPS*DABS(P0))GOTO 2 A(LL)=EF-B2(LL)*P0/P1 NN=LL 2 CONTINUE NR=0 XLIM=EF+EPS CALL RECRTS(A(2),B2(2),NN-1,EPS,XLIM,NR,WK,MU,WK(1,2),WK(1,3)) DO 5 I=1,NR IF(MU(I,1).NE.1)WRITE(6,98)WK(I,1),MU(I,1),EF 98 FORMAT(' *** WARNING ROOT AT ',E16.8,' IS OF MULTIPLICITY', 1 I6,' IN THE CALL OF EDIFF AT EF=',E16.8/, 2 ' INCREASING THE ACCURACY IN THE CALL TO EDIFF MAY CURE THIS') 5 CONTINUE IC=0 CALL RECRTS(A,B2,NN,EPS,EF-EPS,IC,WK(1,2),MU(1,2),WK(1,3),WK(1,4)) DO 6 I=1,IC IF(MU(I,2).NE.1)WRITE(6,98)WK(I,2),MU(I,2),EF 6 CONTINUE IF(IC.EQ.NR)GOTO 7 WRITE(6,44)NR,IC 44 FORMAT(' TROUBLE IN EDIFF ',I4,' ZEROS AND',I4,' POLES - TRY', 1 ' INCREASING THE ACCURACY . EIGENVALUES AS FOLLOWS') WRITE(6,45)(WK(I,1),I=1,NR) WRITE(6,46)(WK(I,2),I=1,IC) 45 FORMAT(4E26.5) 46 FORMAT(E13.5,3E26.5) STOP 7 EDIFF=0.0 DO 8 J=1,NR 8 EDIFF=EDIFF-WK(J,1)+WK(J,2) RETURN END ./ ADD NAME=FENVAL FUNCTION FENVAL(AN,A,B2,ND,LL,NC,ERR,EPS,EB,WK,NW,IW,IFT) DIMENSION A(ND,NC),B2(ND,NC),EB(2),WK(NW,4),IW(NW),FEB(2) ALP=0.5 ITMX=IFT IFT=0 DO 1 I=1,2 FEB(I)=-AN DO 1 J=1,NC ANT=DENQD(EB(I),EB(I),A(1,J),B2(1,J),LL,ALP,EPS,WK,NW,NQ,NE,IW) IF(NE.LE.0)RETURN FEB(I)=FEB(I)-WK(NE,2)*(1.0-ALP) DO 1 K=1,NE 1 FEB(I)=FEB(I)+WK(K,2) IFT=-1 IF(FEB(1)*FEB(2).GT.0)RETURN IFT=0 EF=(EB(1)+EB(2))*0.5 DO 6 IT=1,ITMX IF(ABS(EB(2)-EB(1)).LT.ERR)GOTO 7 ANTL=-AN DANTL=0.0 DO 3 J=1,NC ANT=DENQD(EF,EF,A(1,J),B2(1,J),LL,ALP,EPS,WK,NW,NQ,NE,IW) IF(NE.LT.0)RETURN DANTL=DANTL+ANT ANTL=ANTL-WK(NE,2)*(1.0-ALP) DO 3 K=1,NE 3 ANTL=ANTL+WK(K,2) IF(ANTL*FEB(1).GT.0.0)GOTO 4 EB(2)=EF FEB(2)=ANTL GOTO 5 4 EB(1)=EF FEB(1)=ANTL 5 IF(ABS(DANTL).LT.EPS*ABS(ANTL))DANTL=EPS*ANTL DE=ANTL/DANTL EF=EF-DE IF(ABS(DE).LT.ERR)GOTO 7 IF(EF.LT.EB(1))EF=(EB(1)+EB(2))*0.5 IF(EF.GT.EB(2))EF=(EB(1)+EB(2))*0.5 6 CONTINUE IT=ITMX+1 7 IFT=IT FENVAL=EF RETURN END ./ ADD NAME=NUMC FUNCTION NUMC(A,B2,ALAM,LM1) DIMENSION A(LM1),B2(LM1) 7 NU=0 P0=0.0 P1=1.0 DO 6 I=1,LM1 P2=(A(I)-ALAM)*P1-B2(I)*P0 SC=ABS(P1)+ABS(P2) IF(P2*P1)5,14,4 14 IF(P0*P2)4,5,5 4 NU=NU+1 5 P0=P1/SC 6 P1=P2/SC NUMC=NU RETURN END ./ ADD NAME=NUMD SUBROUTINE NUMD(A,B2,ALAM,LM1,DE) DIMENSION A(LM1),B2(LM1) P0=0.0 PP0=0.0 P1=1.0 PP1=0.0 DO 6 I=1,LM1 P2=(A(I)-ALAM)*P1-B2(I)*P0 PP2=(A(I)-ALAM)*PP1-B2(I)*PP0-P1 SC=ABS(P1)+ABS(P2) P0=P1/SC PP0=PP1/SC PP1=PP2/SC 6 P1=P2/SC DE=P1/PP1 RETURN END ./ ADD NAME=PLYVAL SUBROUTINE PLYVAL(E,A,B2,LM1,P,Q) DIMENSION A(LM1),B2(LM1),P(2),Q(2) P(1)=1.0 P(2)=E-A(1) Q(1)=0.0 Q(2)=B2(1) IF(LM1.LT.2)RETURN DO 1 L=2,LM1 SC=1.0/(ABS(P(2))+ABS(Q(2))) DUM=P(2) P(2)=((E-A(L))*P(2)-B2(L)*P(1))*SC P(1)=DUM*SC DUM=Q(2) Q(2)=((E-A(L))*Q(2)-B2(L)*Q(1))*SC 1 Q(1)=DUM*SC RETURN END ./ ADD NAME=RECAL SUBROUTINE RECAL(HOP,PSI,PMN,M,A,B2,N) DIMENSION PSI(M),PMN(M),A(N),B2(N) DOUBLE PRECISION SUM,DSQRT SUM=B2(1) NM1=N-1 DO 2 J=1,NM1 B2(J)=SUM S=1.0/DSQRT(SUM) CALL HOP(PSI,PMN,A(J)) A(J)=A(J)/B2(J) SUM=0.0D0 DO 1 I=1,M DUM=PSI(I) PSI(I)=(PMN(I)-A(J)*PSI(I))*S PMN(I)=DUM 1 SUM=SUM+PSI(I)*PSI(I) S=S*SUM DO 2 I=1,M 2 PMN(I)=-PMN(I)*S B2(N)=SUM RETURN END ./ ADD NAME=RECNO SUBROUTINE RECNO(HOP,SOP,U,M,NIT,LS,LL,A,B2,PSI,PMN,EMACH) DIMENSION U(M),PSI(M),PMN(M),A(LL),B2(LL) DOUBLE PRECISION SUM COMMON /BLKNNN/SUM NNIT=IABS(NIT) LM1=LL-1 IF(LS.GT.1)GOTO 3 DO 1 I=1,M 1 PMN(I)=0.0 IF (NIT) 11,17,15 11 DO 12 I=1,M 12 U(I)=PSI(I) DO 13 IT=1,NNIT 13 CALL SOP(U,PSI,U) GOTO 17 15 CALL SOP(U,PMN,PSI) DO 16 I=1,M 16 PSI(I)=U(I)-PSI(I) 17 SUM=0.0D0 DO 2 I=1,M 2 SUM=SUM+U(I)*PSI(I) B2(1)=SUM LS=1 3 DO 8 L=LS,LM1 DO 4 I=1,M 4 PMN(I)=-B2(L)*PMN(I) CALL HOP(U,PMN,A(L)) A(L)=A(L)/SUM ANORM=1.0/DSQRT(SUM) DO 5 I=1,M DUM=PMN(I)*ANORM PMN(I)=PSI(I)*ANORM 5 PSI(I)=DUM-A(L)*PMN(I) IF(NNIT.EQ.0)GOTO 10 DO 6 J=1,NNIT 6 CALL SOP(U,PSI,U) 10 SUM=0.0D0 DO 7 I=1,M 7 SUM=SUM+U(I)*PSI(I) IF(SUM.LT.EMACH) GOTO 9 8 B2(L+1)=SUM RETURN 9 LL=-L RETURN END ./ ADD NAME=RECPER SUBROUTINE RECPER(HOP,VOP,W1,W0,A,B,NW,LLIM,NA,NL,AMAT) DIMENSION W1(NW,LLIM),W0(NW,LLIM),A(NA,LLIM),B(NA,LLIM) 1,AMAT(LLIM,LLIM) DOUBLE PRECISION SUM,DSQRT NLM1=NL-1 FAC=1.0/SQRT(B(1,1)) DO 16 N=1,NLM1 DO 2 LL=1,LLIM L=LLIM-LL+1 DO 2 I=1,NW SUM=0.0 DO 1 M=1,L 1 SUM=SUM+B(N,M)*W0(I,L-M+1) 2 W0(I,L)=SUM CALL HOP(W1,W0,AMAT,NW,LLIM,LLIM) A(N,1)=AMAT(1,1) IF(LLIM.LE.1)GOTO 6 DO 3 L=2,LLIM A(N,L)=AMAT(1,L) DO 3 M=2,L 3 A(N,L)=A(N,L)+AMAT(M,L-M+1) CALL VOP(W1,W0(1,2),AMAT,NW,LLIM,LLIM-1) DO 5 L=2,LLIM LL=L-1 DO 4 M=1,LL 4 A(N,L)=A(N,L)+AMAT(M,L-M) 5 A(N,L)=A(N,L)/B(N,1) 6 A(N,1)=A(N,1)/B(N,1) DO 7 L=1,LLIM DO 7 M=1,L DUM=A(N,M) LL=L-M+1 DO 7 I=1,NW 7 W0(I,L)=W0(I,L)-DUM*W1(I,LL) DO 17 L=1,LLIM DO 17 I=1,NW DUM=W0(I,L)*FAC W0(I,L)=-W1(I,L)*FAC 17 W1(I,L)=DUM DO 9 L=1,LLIM SUM=0.0D0 DO 8 M=1,L LL=L-M+1 DO 8 I=1,NW 8 SUM=SUM+W1(I,M)*W1(I,LL) IF(L.EQ.1)FAC=1.0/DSQRT(SUM) 9 B(N+1,L)=SUM IF(LLIM.LE.1)GOTO 16 DO 11 L=2,LLIM LL=L-1 IF(LL.LT.2)GOTO 11 DO 10 M=2,LL 10 B(N+1,L)=B(N+1,L)-B(N+1,M)*B(N+1,L-M+1)/B(N+1,1) 11 B(N+1,L)=B(N+1,L)*0.5 DO 14 L=2,LLIM LL=L-1 DO 14 M=1,LL DUM=B(N+1,L-M+1)/B(N+1,1) DO 14 I=1,NW 14 W1(I,L)=W1(I,L)-DUM*W1(I,M) 16 CONTINUE RETURN END ./ ADD NAME=RECQD SUBROUTINE RECQD(A,B2,LM1,X,WT,M,EPS,WK,MU) DIMENSIONA(LM1),B2(LM1),X(LM1),WT(LM1),WK(LM1),MU(LM1) DIMENSION P(2,3) M=1 CALL RECRTS(A,B2,LM1,EPS,XLIM,M,X,MU,WK,WT) DO 7 II=1,M WT(II)=RECWT(X(II),A,B2,LM1,EPS,0,P,1) 7 CONTINUE RETURN END ./ ADD NAME=RECRTS SUBROUTINE RECRTS(A,B2,LM1,EPS,XLIM,N,X,MULT,BI,NI) DIMENSION A(LM1),B2(LM1),X(LM1),MULT(LM1),BI(LM1),NI(LM1) IF(LM1.GT.1)GOTO 1 IF(N.EQ.0.AND.A(1).GT.XLIM)RETURN X(1)=A(1) MULT(1)=1 N=1 RETURN 1 MULT(1)=0 C C ESTIMATE LIMITS OF ROOTS C CN=AMIN1(A(1)-B2(2),A(LM1)-1.0) NID=LM1-1 IF(NID.LT.2)GOTO 20 DO 2 I=2,NID 2 CN=AMIN1(CN,A(I)-1.0-B2(I+1)) 20 AA=CN IF(N.EQ.0)GOTO 4 CN=AMAX1(A(1)+B2(2),A(LM1)+1.0) IF(NID.LT.2)GOTO 21 DO 3 I=2,NID 3 CN=AMAX1(CN,A(I)+1.0+B2(I+1)) 21 B=CN 4 IF(N.EQ.0) B=XLIM NMAX=AMAX1(ALOG(EPS/(B-AA))/ALOG(0.5)*10,70.0) NA=NUMC(A,B2,AA,LM1) NB=NUMC(A,B2,B,LM1) 5 NL=NA NM=NB N=0 IF (NL.LE.NM) GOTO 19 NE=NL-NM DO 22 I=1,NE X(I)=0.0 22 MULT(I)=0 I=1 NI(1)=NB BI(1)=B NCOUNT=0 C C START BISECTION SEARCH C 9 C=(AA+B)*0.5 NC=NUMC(A,B2,C,LM1) NCOUNT=NCOUNT+1 IF(NCOUNT.GE.NMAX)GOTO 16 IF(NC.LT.NL) GOTO 11 AA=C NA=NC GOTO 15 11 IF(NC.LE.NM)GOTO 14 IF(NI(I).EQ.NC) GOTO 14 I=I+1 BI(I)=C NI(I)=NC 14 B=C NB=NC 15 IF( (NA-NB).GT.1 ) GOTO 36 C DO A FEW NEWTON-RAPHSON STEPS C=(AA+B)*0.5 DO 34 ITR=1,10 CALL NUMD(A,B2,C,LM1,DE) C=C-DE IF(C.GT.B .OR. C.LT.AA) GOTO 9 IF( ABS(DE).LE.EPS) GOTO 35 34 CONTINUE GOTO 9 35 XX=AMAX1(C-EPS,AA) NXX=NUMC(A,B2,XX,LM1) IF( NXX.NE.NA) GOTO 9 XX=AMIN1(C+EPS,B) NXX=NUMC(A,B2,XX,LM1) IF( NXX.NE.NB) GOTO 9 N=N+1 X(N)=C MULT(N)=1 GOTO 37 36 IF(ABS(AA-B).GT.EPS)GOTO 9 16 N=N+1 X(N)=(AA+B)*0.5 MULT(N)=NA-NB IF(NCOUNT.GE.NMAX)MULT(N)=-MULT(N) 37 IF(NB.LE.NM)RETURN IF(I.LE.1)GOTO 19 18 AA=BI(I) NA=NI(I) I=I-1 B=BI(I) NB=NI(I) NL=NL-IABS(MULT(N)) NCOUNT=0 GOTO 9 19 MULT(1)=-MULT(1) RETURN END ./ ADD NAME=RECSUM SUBROUTINE RECSUM(AC,BC,NA,NN,NC,A,B,EPS,WK,NW) DIMENSION AC(NA,NC),BC(NA,NC),WK(NW),A(NA),B(NA) N=IABS(NN) IF(NN.LT.0)GOTO 12 DO 11 I=1,NC 11 AC(N,I)=AC(N-1,I) 12 L1=NC*N L2=2*L1 L3=L2+N NL=0 IC=0 DO 2 I=1,NC CALL RECQD(AC(1,I),BC(1,I),N,WK(IC+1),WK(L1+IC+1),M, 1EPS,WK(L2+1),WK(L3+1)) IF(M.NE.N)GOTO 5 NN=0 DO 1 J=1,N IF(WK(L1+IC+J).GT.0.0)GOTO 1 WK(L1+IC+J)=0.0 NN=NN+1 1 CONTINUE IC=IC+N 2 NL=MAX0(NL,NN) NN=N-NL CALL CFGEN(WK,WK(L1+1),IC,EPS,A,B,NN,WK(L2+1)) RETURN 5 NN=-M RETURN END ./ ADD NAME=RECWT FUNCTION RECWT(E,A,B2,LL,EPS,N,P,NS) DIMENSION A(LL),B2(LL),P(2,3) RECWT=B2(1) IF(LL.EQ.1)RETURN K=NS IF(K.GE.2)GOTO 1 P(1,1)=1.0 P(2,1)=E-A(1) P(1,2)=0.0 P(2,2)=1.0 P(1,3)=0.0 P(2,3)=B2(1) K=2 1 LLIM=LL-IABS(N) IF(LLIM.LT.K)GOTO 4 DO 3 L=K,LLIM SC=ABS(P(1,1))+ABS(P(1,3)) DO 2 J=1,3 DUM=P(2,J) P(2,J)=((E-A(L))*DUM-B2(L)*P(1,J))/SC 2 P(1,J)=DUM/SC 3 P(2,2)=P(2,2)+P(1,1) 4 IF(N)6,5,6 5 RECWT=P(2,3)/P(2,2) GOTO 7 6 RECWT=(P(1,1)*P(2,3)-P(2,1)*P(1,3))/(P(1,1)*P(2,2)-P(1,2) 1 *P(2,1)+P(2,1)*P(2,1)/B2(LL)) 7 IF(RECWT.LT.0.0)RECWT=0.0 IF(N.GE.0)RETURN IF(ABS(P(2,1)).LE.EPS*ABS(P(1,1)))RETURN N=0 A(LL)=E-B2(LL)*P(1,1)/P(2,1) RETURN END ./ ADD NAME=TABAN SUBROUTINE TABAN(E,WT,NPTS,THU,THL,ET,IC,WTT,NET) DIMENSION E(NPTS),WT(NPTS),ET(NET),IC(NET),WTT(NET) C C CONTROL IS MODERATED THROUGHOUT BY THE VARIABLE IREG WHICH C INDICATES IN WHICH REGION THE PREVIOUS POINT LIES. C IREG=1 FUNCTION VALUE LESS THAN LOWER THRESHOLD, THL C IREG=2 FUNCTION VALUE LIES BETWEEN THL AND THU (IN WINDOW) C IREG=3 FUNCTION VALUE LIES ABOVE UPPER THRESHOLD, THU C I=0 IREG=2 IF(WT(1).GT.THU)IREG=3 IF(WT(1).LT.THL)IREG=1 INCR=0 W0=2.0*WT(1)-WT(2) IF(WT(1).GT.W0)INCR=4 IF(WT(1).LT.W0)INCR=-4 IF(INCR.EQ.0)EFL=E(1) DO 10 II=1,NPTS GOTO(1,3,7),IREG 1 IF(WT(II).LT.THL)GOTO 9 IF(WT(II).LT.THU)GOTO 2 IREG=3 I=I+1 ET(I)=(E(II)+E0)*0.5 IC(I)=3 GOTO 9 2 IREG=2 I=I+1 ET(I)=(E(II)+E0)*0.5 IC(I)=1 INCR=4 GOTO 85 3 NCR=0 IF(WT(II).LT.W0)NCR=-4 IF(WT(II).GT.W0)NCR=4 IF(NCR.EQ.0)GOTO 9 IF(NCR.EQ.INCR)GOTO 4 I=I+1 ET(I)=(E0+EFL)*0.5 WTT(I)=W0 IC(I)=NCR IF(I.GE.NET)GOTO 11 4 INCR=NCR IF(WT(II).LT.THU)GOTO 5 IREG=3 I=I+1 ET(I)=(E(II)+E0)*0.5 IC(I)=2 GOTO 9 5 IF(WT(II).GT.THL)GOTO 85 IREG=1 I=I+1 ET(I)=(E(II)+E0)*0.5 IC(I)=-1 GOTO 9 7 IF(WT(II).GT.THU)GOTO 9 IF(WT(II).GT.THL)GOTO 8 IREG=1 I=I+1 ET(I)=(E(II)+E0)*0.5 IC(I)=-3 GOTO 9 8 IREG=2 I=I+1 ET(I)=(E(II)+E0)*0.5 IC(I)=-2 INCR=-4 85 W0=WT(II) EFL=E(II) 9 IF(I.GE.NET)GOTO 11 E0=E(II) 10 CONTINUE II=-I 11 NET=-II RETURN END ./ ADD NAME=TERMGN SUBROUTINE TERMGN(A,B2,LL,EPS,ERR,ITMX,AA,RNG,WB,NBP1,AM,BM2, 1 IC,WK,NW,BWK,NBD,IWK) DIMENSION A(LL),B2(LL),AA(NBP1),RNG(NBP1),WB(NBP1),AM(LL) 1 ,BM2(LL),IC(NW),WK(NW,2,NBP1),BWK(NBD,4),IWK(LL) 2 ,P(2,3) KD=1 NB=NBP1-1 NET=NB*2 CALL BNDEST(A,B2,LL,EPS,AA,RNG,NB,BWK,BWK(1,2),IC,NET,DE,WK,NW) DEL=DE*0.5 DO 4 ITE=1,ITMX,5 DO 1 I=1,NB ID=I+I BWK(ID-1,3)=DEL BWK(ID,3)=DEL BWK(ID-1,4)=AA(I) 1 BWK(ID,4)=RNG(I) CALL BNDWT(AA,RNG,WB,NB,A,B2,LL,EPS,WK,NW,IWK) DO 2 ITER=1,5 CALL CFGPGN(AA,RNG,WB,NB+1,KD,EPS,AM,BM2,LL,WK,NW) CALL BNDREF(DEL,AM,BM2,LL,EPS,AA,RNG,NB,BWK,NBD,IC,NET) DEL=DEL*0.5 2 CONTINUE DEL=0.0 DO 3 I=1,NB ID=I+I DEL=AMAX1(DEL,ABS(AA(I)-BWK(ID-1,4))) 3 DEL=AMAX1(DEL,ABS(RNG(I)-BWK(ID,4))) IF(DEL.LT.ERR)GOTO 5 DEL=DEL*0.5 4 CONTINUE 5 NBP1=NB+1 CALL CFGPGN(AA,RNG,WB,NBP1,KD,EPS,AM,BM2,LL,WK,NW) ERR=DEL RETURN END ./ ADD NAME=WTMIN SUBROUTINE WTMIN(AA,BB,A,B2,LL,EPS,ACC,ITMX,EM,FEM) DIMENSION A(LL),B2(LL),P(2,3),E(4),FE(4) E(1)=AA E(4)=BB FE(1)=RECWT(AA,A,B2,LL,EPS,1,P,1) FE(4)=RECWT(BB,A,B2,LL,EPS,1,P,1) E(2)=(AA+BB)*0.5 FE(2)=RECWT(E(2),A,B2,LL,EPS,1,P,1) IF(FE(2).GT.FE(1))GOTO 6 IF(FE(2).GT.FE(4))GOTO 6 IC=3 DO 4 I=1,ITMX E(IC)=(E(IC+1)+E(IC-1))*0.5 FE(IC)=RECWT(E(IC),A,B2,LL,EPS,1,P,1) IF(FE(2).GT.FE(3))GOTO 2 E(4)=E(3) FE(4)=FE(3) IF(IC.EQ.2)GOTO 3 E(3)=E(2) FE(3)=FE(2) GOTO 3 2 E(1)=E(2) FE(1)=FE(2) IF(IC.EQ.3)GOTO 3 E(2)=E(3) FE(2)=FE(3) 3 ERR=ABS(E(4)-E(1)) IF(ERR.LT.ACC)GOTO 5 4 IC=5-IC IT=0 5 EM=(E(4)+E(1))*0.5 FEM=RECWT(EM,A,B2,LL,EPS,1,P,1) ITMX=IT RETURN 6 IT=1 IF(FE(1).GT.FE(4))IT=4 EM=E(IT) FEM=FE(IT) ITMX=-2+IT/4 RETURN END ./ ENDUP ./ ADD NAME=ADDAT SUBROUTINE ADDAT(CRD,ND,NAT,EV,IZP,MM,NN,NND,NM,NGBR,NE,EE,NP 1,VEC,IW,NED,OVPAR,HCAL) DIMENSION CRD(ND,NAT),VEC(3,NED),V(3),EE(NP,NP,NED) 1 ,IZP(NAT),NN(NND,NM),IW(2,NED),MM(NND,NM) INTEGER*2 IZP,NN,IW,MM LOGICAL EV EXTERNAL EV,NGBR,OVPAR CALL NNCAL(CRD,ND,NAT,IZP,NN,NND,NM,NGBR) NE1=1+NE CALL MMCAL(CRD,ND,NAT,NN,NND,NM,EV,IZP,NE,MM,VEC,IW) IF(NE1.GT.NE)RETURN DO 1 K=NE1,NE II=IW(1,K) IJ=IW(2,K) 1 CALL HCAL(VEC(1,K),II,IJ,EE(1,1,K),OVPAR) DO 2 I=1,3 2 V(I)=0.0 DO 5 I=1,NAT IF(IZP(I).GT.0)GOTO 5 IT=-IZP(I) DO 3 K=1,NE IF(.NOT.EV(V,VEC(1,K)))GOTO 3 IF(IT.EQ.IW(1,K))GOTO 4 3 CONTINUE NE=NE+1 K=NE DO 6 J=1,3 6 VEC(J,K)=0.0 IW(1,K)=IT IW(2,K)=IT CALL HCAL(VEC(1,K),IT,IT,EE(1,1,K),OVPAR) 4 MM(I,1)=K 5 CONTINUE RETURN END ./ ADD NAME=BCCBFE INTEGER FUNCTION BCCBFE(I,J,R2,DD) C C MICHAEL YOU'S IBONDS FOR A BCC LATTICE C FIRST AND SECOND NEAREST NEIGHBOUR DISTANCES C AND THE PETTIFOR PARAMETERS DD. C DIMENSION DD(13) IF(R2-4.5)1,3,3 3 BCCBFE=0 RETURN 1 BCCBFE=1 IF(R2-3.5)2,4,4 4 DD(1)=-0.03195 DD(2)=0.02130 DD(3)=-0.00537 RETURN 2 IF(R2.LT.1.0E-4)GOTO 5 DD(1)=-0.06560 DD(2)=0.04373 DD(3)=-0.01093 RETURN 5 DD(11)=0.0 RETURN END ./ ADD NAME=BCCLAT SUBROUTINE BCCLAT(CRD,NDIM,IZP,NAT,NX,NY,NZ,NTYPE) DIMENSION CRD(NDIM,NAT),IZP(NAT) INTEGER*2 IZP NA=NAT/2 N=0 DO 1 I=2,NX,2 DO 1 J=2,NY,2 DO 1 K=2,NZ,2 N=N+1 IF(N.GT.NA)GOTO 2 IZP(N)=NTYPE CRD(1,N)=FLOAT(I)-1.0 CRD(2,N)=FLOAT(J)-1.0 CRD(3,N)=FLOAT(K)-1.0 1 CONTINUE NAT=N+N GOTO 4 2 WRITE(6,3)I,J,K NAT=0 3 FORMAT(' INCREASE NAT IN THE CALL TO BCCLAT - LATTICE GENERATED' 1 ,' AS FAR AS ',3I4) 4 DO 5 I=1,N K=I+N IZP(K)=NTYPE DO 5 J=1,3 5 CRD(J,K)=CRD(J,I)+1.0 RETURN END ./ ADD NAME=EQUIV LOGICAL FUNCTION EQUIV(V,W) DIMENSION V(3),W(3) EQUIV=.FALSE. DO 1 I=1,3 IF(ABS(V(I)-W(I)).GT.1.0E-4)RETURN 1 CONTINUE EQUIV=.TRUE. RETURN END ./ ADD NAME=FCCBND INTEGER FUNCTION FCCBND(I,J,R2,DD) DIMENSION DD(13) IF(R2-3.0)1,3,3 3 FCCBND=0 RETURN 1 FCCBND=1 IF(R2.LT.1.0E-4)GOTO 2 DD(1)=-0.027784 DD(2)=0.012535 DD(3)=-0.001554 RETURN 2 DD(11)=0.0 RETURN END ./ ADD NAME=FCCLAT SUBROUTINE FCCLAT(CRD,NDIM,IZP,NAT,NX,NY,NZ,NTYPE) DIMENSION CRD(NDIM,NAT),IZP(NAT) INTEGER*2 IZP N=0 KK=2 DO 5 K=1,NZ II=KK RK=FLOAT(K) DO 2 J=1,NY RJ=FLOAT(J) DO 1 I=II,NX,2 RI=FLOAT(I) N=N+1 IF(N.GT.NAT)GOTO 3 IZP(N)=NTYPE CRD(1,N)=RI CRD(2,N)=RJ 1 CRD(3,N)=RK 2 II=3-II 5 KK=3-KK NAT=N RETURN 3 NAT=0 WRITE(6,4)I,J,K 4 FORMAT(' NAT TOO SMALL IN FCCLAT- LATTICE GENERATED UP TO',3I4) RETURN END ./ ADD NAME=MMCAL SUBROUTINE MMCAL(CRD,NDIM,NAT,NN,ND,NM,EV,IZP,NMAT,MM,VEC,IW) INTEGER*2 NN,MM,IZP,IW DIMENSION CRD(NDIM,NAT),VEC(NDIM,NMAT) 1 ,NN(ND,NM),MM(ND,NM),IZP(NAT),IW(2,NMAT) LOGICAL EV COMMON /BLKNNM/NNMAT IFD=NMAT IADD=1 DO 1 II=1,NAT IF(IZP(II).LT.0)GOTO 2 1 CONTINUE IADD=0 NNMAT=NMAT JL=NN(1,1) DO 13 JJ=2,JL J=NN(1,JJ) DO 11 I=1,NDIM VEC(I,NNMAT)=CRD(I,J)-CRD(I,1) 11 VEC(I,JJ-1)=VEC(I,NNMAT) MM(1,JJ)=JJ-1 IW(1,JJ-1)=IZP(1) 13 IW(2,JJ-1)=IZP(J) IFD=JL-1 2 ISTRT=2-IADD DO 30 I=ISTRT,NAT IAZ=IZP(I) IAZ=IABS(IAZ) JL=NN(I,1) IF(JL.LT.2)GOTO 30 DO 20 JJ=2,JL J=NN(I,JJ) IF(IADD.EQ.0)GOTO 3 IF(IZP(I).GT.0.AND.IZP(J).GT.0)GOTO 20 3 JAZ=IZP(J) JAZ=IABS(JAZ) DO 14 K=1,NDIM 14 VEC(K,NNMAT)=CRD(K,J)-CRD(K,I) DO 16 L=1,IFD IF(IW(1,L).NE.IAZ)GOTO 16 IF(IW(2,L).NE.JAZ)GOTO 16 IF(EV(VEC(1,NNMAT),VEC(1,L)))GOTO 19 16 CONTINUE IFD=IFD+1 IF(IFD.GE.NNMAT)GOTO 34 18 FORMAT(' TOO MANY DIFFERENT VECTORS - INCREASE NMAT') DO 29 M=1,NDIM 29 VEC(M,IFD)=VEC(M,NNMAT) IW(1,IFD)=IAZ IW(2,IFD)=JAZ L=IFD 19 MM(I,JJ)=L 20 CONTINUE 30 CONTINUE NMAT=IFD RETURN 34 WRITE(6,18) WRITE(6,41)I 41 FORMAT(26H HAMILTONIAN MAP AS FAR AS,I6,7HTH SITE) NMAT=0 RETURN END ./ ADD NAME=NNCAL SUBROUTINE NNCAL(CRD,NDIM,NAT,IZP,NN,ND,NM,NGBR) INTEGER*2 IZP,NN DIMENSION CRD(NDIM,NAT),DUM(13),IZP(NAT),NN(ND,NM) NNMAX=0 IADD=1 IF(IZP(1).LT.0)NN(1,1)=1 DO 11 II=1,NAT IF(IZP(II).LT.0)GOTO 12 11 CONTINUE NN(1,1)=1 DO 1 I=1,NAT DO 1 J=2,NM 1 NN(I,J)=0 IADD=0 II=2 12 IF(II.LE.1)II=2 50 FORMAT(' FROM NNCAL') DO 5 I=II,NAT IF(IZP(I)*IADD.GT.0)GOTO 5 NN(I,1)=1 IIP=IZP(I) IIP=IABS(IIP) ILJ=I-1 DO 4 J=1,ILJ JJP=IZP(J) JJP=IABS(JJP) R2=0.0 DO 2 L=1,3 DUM(L)=CRD(L,I)-CRD(L,J) 2 R2=R2+DUM(L)*DUM(L) ID=NGBR(IIP,JJP,R2,DUM) IF(ID.EQ.0)GOTO 4 ID=NN(I,1)+1 NN(I,1)=ID NN(I,ID)=J NNMAX=MAX0(NNMAX,ID) ID=NN(J,1)+1 NN(J,1)=ID NN(J,ID)=I NNMAX=MAX0(NNMAX,ID) IF(NNMAX.GT.NM)GOTO 33 3 FORMAT(20H TOO MANY NEIGHBOURS) 4 CONTINUE 5 CONTINUE NM=NNMAX RETURN 33 WRITE(6,3) WRITE(6,38)I 38 FORMAT(24H NEIGHBOUR MAP AS FAR AS,I6,7HTH SITE) DO 39 L=1,NAT MMM=NN(L,1) 39 WRITE(6,40)L,(NN(L,M),M=1,MMM) 40 FORMAT(21I4) STOP END ./ ADD NAME=ONION SUBROUTINE ONION(NN,ND,NM,IZERO,NAT,IST,NNS,IW) INTEGER*2 NN,IZERO,IST,IW DIMENSION NN(ND,NM),IZERO(NAT),IST(NNS),IW(NAT) NS=NNS DO 1 I=1,NAT 1 IZERO(I)=0 DO 2 I=1,NS II=IST(I) IW(I)=II 2 IZERO(II)=1 NF=1 ISH=2 NCD=NS 3 DO 4 I=NF,NS II=IW(I) JL=NN(II,1) DO 4 J=2,JL JJ=NN(II,J) IF(IZERO(JJ).NE.0)GOTO 4 IZERO(JJ)=ISH NCD=NCD+1 IW(NCD)=JJ IF(NCD.EQ.NAT)GOTO 5 4 CONTINUE ISH=ISH+1 NF=NS+1 NS=NCD GOTO 3 5 CONTINUE RETURN END ./ ADD NAME=ORPEEL SUBROUTINE ORPEEL(NSTRT,NORB,NO,MM,NN,ND,ID,EE,NP,NE,NED,MEM) INTEBGER*2 MM,NN,MEM DIMENSION EE(NP,NP,NED),MM(ND,ID),NN(ND,ID),MEM(3,ID) I=NN(NSTRT,1) IF(NO.EQ.NP)GOTO 8 IF(NO.NE.1)GOTO 5 MEM(2,1)=NE DO 1 L=1,I II=MM(NSTRT,L) NE=NE+1 IF(NE.GT.NED)GOTO 11 MEM(1,L)=II MM(NSTRT,L)=NE DO 1 M=1,NP DO 1 N=1,NP 1 EE(M,N,NE)=EE(M,N,II) DO 4 L=2,I J=NN(NSTRT,L) KK=NN(J,1) DO 2 K=2,KK JJ=NN(J,K) IF(JJ.EQ.NSTRT)GOTO 3 2 CONTINUE GOTO 13 3 II=MM(J,K) NE=NE+1 IF(NE.GT.NED)GOTO 11 MEM(2,L)=II MEM(3,L)=K MM(J,K)=NE DO 4 M=1,NP DO 4 N=1,NP 4 EE(M,N,NE)=EE(M,N,II) 5 II=MEM(2,1)+1 DO 55 M=1,NP 55 EE(M,NORB,II)=0.0 DO 6 L=1,I II=MEM(2,1)+L DO 6 M=1,NP 6 EE(NORB,M,II)=0.0 DO 7 L=2,I II=II+1 DO 7 M=1,NP 7 EE(M,NORB,II)=0.0 RETURN 8 DO 9 L=1,I 9 MM(NSTRT,L)=MEM(1,L) DO 10 L=2,I J=NN(NSTRT,L) K=MEM(3,L) 10 MM(J,K)=MEM(2,L) NE=MEM(2,1) RETURN 11 WRITE(6,12) 12 FORMAT(' TOO MANY INTERACTION MATRICES FOR STORE IN ORPEEL : 1 INCREASE DIMENSION AND NED IN CALLING ROUTINE') STOP 13 WRITE(6,14)NSTRT,J,NSTRT 14 FORMAT(' INCONSISTENCY IN NEIGHBOUR MAP : ATOM',I6, 1' IS NOT A NGBR OF ATOM',I6,' WHICH IS A NGBR OF ATOM',I6) STOP END ./ ADD NAME=OUT SUBROUTINE OUT(IO,IA,IE,IZP,IW,VEC,NED,NE,NAT,MM,NN,ND,NM,EE,NP) INTEGER*2 MM,NN,IZP,IW DIMENSION MM(ND,NM),NN(ND,NM),IZP(NAT),IW(2,NED),VEC(3,NED) 1 ,EE(NP,NP,NED) WRITE(IO,11) 11 FORMAT(' NEAREST NEIGHBOUR MAP'/, 1' ATOM TYPE CONNECTIVITY',5X,'NEIGHBOURS') 12 FORMAT(I5,2I3,16I5) DO 13 I=IA,NAT ID=NN(I,1) ID=MIN0(ID,17) WRITE(IO,12)I,IZP(I),ID,(NN(I,J),J=2,ID) IF(ID.EQ.NN(I,1))GOTO 13 ID=NN(I,1) WRITE(IO,125)(NN(I,J),J=18,ID) 125 FORMAT(12X,16I5) 13 CONTINUE WRITE(IO,14)NE 14 FORMAT(I6,' NON-EQUIVALENT VECTORS FOUND',/6X,'ATOM',5X 1,'INDEX OF INTERACTION VECTORS') DO 15 I=IA,NAT ID=NN(I,1) ID=MIN0(ID,16) WRITE(IO,16)I,(MM(I,J),J=1,ID) IF(ID.EQ.NN(I,1))GOTO 15 ID=NN(I,1) WRITE(IO,145)(MM(I,J),J=17,ID) 145 FORMAT(5X,15I5) 15 CONTINUE 16 FORMAT(17I5) IF(IE.EQ.0)RETURN WRITE(IO,19) DO 17 I=IE,NE WRITE(IO,18)I,IW(1,I),IW(2,I),(VEC(K,I),K=1,3) 1,(EE(1,L,I),L=1,NP) IF(NP.EQ.1)GOTO 17 DO 7 K=2,NP 7 WRITE(IO,21)(EE(K,L,I),L=1,NP) 17 CONTINUE 18 FORMAT(/I8,3X,2I3,3X,3F6.2,9F7.3) 19 FORMAT(/4X,' INDEX LATTICE TYPES VECTOR',30X,'MATRIX') 21 FORMAT(38X,9F7.3) RETURN END ./ ADD NAME=PEEL SUBROUTINE PEEL(CRD,NDIM,NAT,NN,ND,NM,IST,NS,IZP,IZERO,NSH,IW) INTEGER*2 NN,IST,IZP,IZERO,IW DIMENSION NN(ND,NM),IST(NS),IZP(NAT),IZERO(NAT),IW(NAT) 1 ,CRD(NDIM,NAT) NA=0 DO 1 I=1,NAT IW(I)=0 IF(IZERO(I).GT.NSH)GOTO 1 NA=NA+1 IW(I)=NA 1 CONTINUE DO 4 I=1,NAT II=IW(I) IF(II.EQ.0)GOTO 4 DO 2 J=1,3 2 CRD(J,II)=CRD(J,I) IZERO(II)=IZERO(I) IZP(II)=IZP(I) JJ=NN(I,1) KV=1 IF(JJ.GE.2) GOTO 38 JJ=2 GOTO 37 38 DO 3 J=2,JJ K=NN(I,J) KK=IW(K) IF(KK.EQ.0)GOTO 3 KV=KV+1 NN(II,KV)=KK 3 CONTINUE JJ=KV+1 IF(JJ.GT.NM)GOTO 36 37 DO 35 J=JJ,NM 35 NN(II,J)=0 36 NN(II,1)=KV 4 CONTINUE NAT=NA DO 5 I=1,NS J=IST(I) 5 IST(I)=IW(J) RETURN END ./ ADD NAME=SCAN SUBROUTINE SCAN(NN,ND,NNMX,N0,NAT,NON,SUB) INTEGER*2 NN DIMENSION NN(ND,NNMX),IA(3),NA(3) DO 2 I=N0,NAT NA(1)=1 IA(1)=I CALL SUB(IA,NA,1) NJ=NN(I,1) DO 2 JJ=2,NJ NA(2)=JJ J=NN(I,JJ) IA(2)=J CALL SUB(IA,NA,2) IF(NON.EQ.1)GOTO 2 NK=NN(J,1) DO 1 KK=2,NK NA(3)=KK K=NN(J,KK) IA(3)=K IF(K.EQ.I)GOTO 1 CALL SUB(IA,NA,3) 1 CONTINUE 2 CONTINUE RETURN END ./ ADD NAME=SCAN1 SUBROUTINE SCAN1(NN,ND,NNMX,N0,NAT,NON,SUB) INTEGER*2 NN DIMENSION NN(ND,NNMX),IA(3),NA(3) DO 2 I=N0,NAT NA(1)=1 IA(1)=I CALL SUB(IA,NA,1) NJ=NN(I,1) DO 2 JJ=2,NJ NA(2)=JJ J=NN(I,JJ) IA(2)=J CALL SUB(IA,NA,2) IF(NON.EQ.1)GOTO 2 NK=NN(J,1) DO 1 KK=2,NK NA(3)=KK K=NN(J,KK) IA(3)=K IF(K.EQ.I)GOTO 1 CALL SUB(IA,NA,3) 1 CONTINUE 2 CONTINUE RETURN END ./ ADD NAME=SELFD SUBROUTINE SELFD(I,J,R2,IOVPAR,EM,NE) DIMENSION PARM(13),EM(NE,5) ID=IOVPAR(I,J,R2,PARM) DO 11 L=1,5 DO 10 M=1,5 10 EM(L,M)=0.0 11 EM(L,L)=PARM(11) RETURN END ./ ADD NAME=SELFP SUBROUTINE SELFP(I,J,R2,IOVPAR,EM,NE) DIMENSION PARM(13),EM(NE,3) ID=IOVPAR(I,J,R2,PARM) DO 11 L=1,3 DO 10 M=1,3 10 EM(L,M)=0.0 11 EM(L,L)=PARM(12) RETURN END ./ ADD NAME=SELFS SUBROUTINE SELFS(I,J,R2,IOVPAR,EM,NE) DIMENSION PARM(13),EM(NE,1) ID=IOVPAR(I,J,R2,PARM) EM(1,1)=PARM(13) RETURN END ./ ADD NAME=SETUP SUBROUTINE SETUP(CRD,ND,NAT,EV,NTYPE,IZP,MM,NN,NND,NM,HCAL,NGBR 1,IOVPAR,EE,NP,NED,NE,VEC,IW) INTEGER*2 MM(NND,NM),NN(NND,NM),IZP(NAT),IW(2,NED) LOGICAL EV DIMENSION CRD(ND,NAT),VEC(3,NED),EE(NP,NP,NED) EXTERNAL NGBR,IOVPAR,EV C C COMPUTE NEIGHBOUR MAP C CALL NNCAL(CRD,ND,NAT,IZP,NN,NND,NM,NGBR) C C CONSTRUCT HAMILTONIAN MAP MM AND EE C NE=NED CALL MMCAL(CRD,ND,NAT,NN,NND,NM,EV,IZP,NE,MM,VEC,IW) IF(NE.EQ.0)GOTO 4 DO 1 I=1,NAT II=IZP(I)+NE 1 MM(I,1)=II C C LOAD HAMILTONIAN MATRICES GENERATED BY EACH DISTINCT VECTOR C DO 2 K=1,NE II=IW(1,K) IJ=IW(2,K) 2 CALL HCAL(VEC(1,K),II,IJ,EE(1,1,K),IOVPAR) IF(NE+NTYPE.GT.NED)GOTO 4 DO 3 KK=1,NTYPE K=NE+KK C C COPY 'SELF ENERGY' MATRICES C DO 8 I=1,3 8 VEC(I,K)=0.0 IW(1,K)=KK IW(2,K)=KK CALL HCAL(VEC(1,K),KK,KK,EE(1,1,K),IOVPAR) 3 CONTINUE NE=NE+NTYPE RETURN 4 WRITE(6,5)NE,KK 5 FORMAT(' INCREASE DIMENSONS OF EE ETC. AS',I4,'+',I4, 1' INTERATION MATRICES ARE REQUIRED') STOP END ./ ADD NAME=SKDD SUBROUTINE SKDD(X,X2,I,J,R2,IOVPAR,EM,NE) DIMENSION X(6),X2(6),PARM(13),EM(NE,5),E(15,3),DM(15),DD(3) R3=SQRT(3.0) D3=X2(1)+X2(2) D4=X2(3)-0.5*D3 D5=X2(1)-X2(2) ID=IOVPAR(I,J,R2,PARM) DO 3 L=1,3 E(L,1)=X2(L)*X2(L+1) E(L,2)=X2(L)+X2(L+1)-4.0*E(L,1) E(L,3)=X2(L+2)+E(L,1) 3 E(L,1)=3.0*E(L,1) E(4,1)=D5*D5 E(4,2)=D3-E(4,1) E(4,3)=X2(3)+0.25*E(4,1) E(4,1)=0.75*E(4,1) E(5,1)=D4*D4 E(5,2)=3.0*X2(3)*D3 E(5,3)=D3*D3*0.75 DD(1)=X(1)*X(3) DD(2)=X(2)*X(1) DD(3)=X(3)*X(2) DO 4 L=1,2 E(L+5,1)=3.0*X2(L+1)*DD(L) E(L+5,2)=DD(L)*(1.0-4.0*X2(L+1)) 4 E(L+5,3)=DD(L)*(X2(L+1)-1.0) E(8,1)=DD(1)*D5*1.5 E(8,2)=DD(1)*(1.0-2.0*D5) E(8,3)=DD(1)*(0.5*D5-1.0) E(9,1)=D5*0.5*D4*R3 E(9,2)=-D5*X2(3)*R3 E(9,3)=D5*0.25*(1.0+X2(3))*R3 E(10,1)=X2(1)*DD(3)*3.0 E(10,2)=(0.25-X2(1))*DD(3)*4.0 E(10,3)=DD(3)*(X2(1)-1.0) E(11,1)=1.5*DD(3)*D5 E(11,2)=-DD(3)*(1.0+2.0*D5) E(11,3)=DD(3)*(1.0+0.5*D5) E(13,3)=0.5*D5*DD(2) E(13,2)=-2.0*DD(2)*D5 E(13,1)=E(13,3)*3.0 E(12,1)=D4*DD(1)*R3 E(14,1)=D4*DD(3)*R3 E(15,1)=D4*DD(2)*R3 E(15,2)=-2.0*R3*DD(2)*X2(3) E(15,3)=0.5*R3*(1.0+X2(3))*DD(2) E(14,2)=R3*DD(3)*(D3-X2(3)) E(14,3)=-R3*0.5*DD(3)*D3 E(12,2)=R3*DD(1)*(D3-X2(3)) E(12,3)=-R3*0.5*DD(1)*D3 DO 11 L=1,15 DM(L)=0.0 DO 11 M=1,3 11 DM(L)=DM(L)+E(L,M)*PARM(M) DO 12 IR=1,5 DO 12 IS=1,IR II=IR-IS K=5*II-(II*(II-1))/2+IS EM(IR,IS)=DM(K) 12 EM(IS,IR)=DM(K) RETURN END ./ ADD NAME=SKPD SUBROUTINE SKPD(X,X2,I,J,R2,IOVPAR,EM,EMT,NE) DIMENSION X(6),X2(6),PARM(13),EM(NE,5),EMT(NE,5),EPD(13,2),DM(15) C R3=SQRT(3.0) D3=X2(1)+X2(2) D4=X2(3)-0.5*D3 D5=X2(1)-X2(2) D6=X(1)*X(2)*X(3) ID=IOVPAR(I,J,R2,PARM) DO 10 L=1,3 EPD(L,1)=R3*X2(L)*X(L+1) EPD(L,2)=X(L+1)*(1.0-2.0*X2(L)) EPD(L+4,1)=R3*X2(L)*X(L+2) EPD(L+4,2)=X(L+2)*(1.0-2.0*X2(L)) EPD(L+7,1)=0.5*R3*X(L)*D5 10 EPD(L+10,1)=X(L)*D4 EPD(4,1)=R3*D6 EPD(4,2)=-2.0*D6 EPD(8,2)=X(1)*(1.0-D5) EPD(9,2)=-X(2)*(1.0+D5) EPD(10,2)=-X(3)*D5 EPD(11,2)=-R3*X(1)*X2(3) EPD(12,2)=-R3*X(2)*X2(3) EPD(13,2)=R3*X(3)*D3 DO 11 L=1,15 11 DM(L)=0.0 DO 12 M=1,2 DM(1)=DM(1)+EPD(1,M)*PARM(M+3) DM(2)=DM(2)+EPD(6,M)*PARM(M+3) DM(3)=DM(3)+EPD(4,M)*PARM(M+3) DM(5)=DM(5)+EPD(2,M)*PARM(M+3) DM(6)=DM(6)+EPD(7,M)*PARM(M+3) DM(7)=DM(7)+EPD(5,