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,M)*PARM(M+3) DM(9)=DM(9)+EPD(3,M)*PARM(M+3) DO 12 L=8,13 12 DM(L+2)=DM(L+2)+EPD(L,M)*PARM(M+3) DM(4)=DM(3) DM(8)=DM(3) DO 13 IR=1,5 DO 13 IS=1,3 K=3*(IR-1)+IS EMT(IR,IS)=-DM(K) 13 EM(IS,IR)=DM(K) RETURN END ./ 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) C ID=IOVPAR(I,J,R2,PARM) DO 10 L=1,3 EPP(L)=X2(L) 10 EPP(L+3)=X(L)*X(L+1) DO 11 L=1,3 HP=EPP(L) 11 DM(L)=HP*PARM(6)+(1.0-HP)*PARM(7) DO 12 L=4,6 12 DM(L)=EPP(L)*(PARM(6)-PARM(7)) DO 13 IR=1,3 DO 13 IS=1,IR II=IR-IS K=3*II-(II*(II-1))/2+IS EM(IS,IR)=DM(K) 13 EM(IR,IS)=DM(K) RETURN END ./ 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) C R3=SQRT(3.0) D4=X2(3)-0.5*(X2(1)+X2(2)) D5=X2(1)-X2(2) ID=IOVPAR(I,J,R2,PARM) C DO 10 L=1,3 10 ES(L)=R3*X(L)*X(L+1) ES(4)=0.5*R3*D5 ES(5)=D4 DO 11 L=1,5 EM(1,L)=ES(L)*PARM(8) 11 EMT(L,1)=EM(1,L) RETURN END ./ 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) C ID=IOVPAR(I,J,R2,PARM) DO 10 L=1,3 EM(1,L)=X(L)*PARM(9) 10 EMT(L,1)=-EM(1,L) RETURN END ./ ADD NAME=SKSS SUBROUTINE SKSS(X,X2,I,J,R2,IOVPAR,EM,NE) DIMENSION X(6),X2(6),PARM(13),EM(1,1) C ID=IOVPAR(I,J,R2,PARM) EM(1,1)=PARM(10) RETURN END ./ ADD NAME=SLKODE SUBROUTINE SLKODE(DUM,I,J,EM,NGBR) DIMENSION DUM(3),EM(5,5) DIMENSION DD(3),X(6),X2(6),DM(15),E(15,3) COMMENT D ELECTRON HAMILTONIAN CALCULATION R2=0.0 DO 1 L=1,3 X(L)=DUM(L) X2(L)=X(L)*X(L) 1 R2=R2+X2(L) IF(R2.LT.1.0E-4)GOTO 3 R2I=1.0/R2 RI=SQRT(R2I) DO 2 L=1,3 X(L)=X(L)*RI X(L+3)=X(L) X2(L)=X2(L)*R2I 2 X2(L+3)=X2(L) CALL SKDD(X,X2,I,J,R2,NGBR,EM,5) RETURN 3 CALL SELFD(I,J,R2,NGBR,EM,5) RETURN END ./ 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) IF(ISTR.EQ.0)RETURN IO=-ISTR IF(ISTR.GT.0)GOTO 71 READ(IO,1)NNC,NNAT 1 FORMAT(2I5) IF(NC-NNC)11,51,31 11 WRITE(6,21)NC,NNC 21 FORMAT(' NC ARGUMENT(',I4,') IS LESS THAN NC INPUT(',I4,')') STOP 31 WRITE(6,41)NC,NNC 41 FORMAT(' WARNING NC INPUT (',I4,' NOT EQUAL TO ARGUMENT (',I4,')') 51 IF(NNAT.LE.NAT)GOTO 81 WRITE(6,61)NAT,NNAT 61 FORMAT(' NAT ARGUMENT(',I4,') LESS THAN NAT INPUT(',I4,')') STOP 71 WRITE(ISTR,1)NC,NAT NNC=NC NNAT=NAT 81 NAT=NNAT CALL RS4(ISTR,CRD,NNC*NAT) CALL IS2(ISTR,IZP,NAT) IF(ISTR.GT.0)GOTO 42 READ(IO,1)NNM IF(NND.GE.NAT)GOTO 22 WRITE(6,12)NND,NAT 12 FORMAT(' ND ARGUMENT(',I4,') IS LESS THAN NO. OF ATOMS(',I4,')') STOP 22 IF(NNM.LE.NM)GOTO 52 WRITE(6,32)NM,NNM 32 FORMAT(' NM ARGUMENT(',I4,') LESS THAN NM INPUT(',I4,')') STOP 42 WRITE(ISTR,1)NM NNM=NM 52 NM=NNM DO 62 I=1,NNM CALL IS2(ISTR,MM(1,I),NAT) 62 CALL IS2(ISTR,NN(1,I),NAT) IF(NE.EQ.0)RETURN IF(ISTR.GT.0)GOTO 43 READ(IO,1)NNP,NNE IF(NNP.EQ.NP)GOTO 23 WRITE(6,13)NP,NNP 13 FORMAT(' NP ARGUMENT(',I4,') IS NOT EQUAL TO NP INPUT(',I4,')') STOP 23 IF(NNE.LE.NE)GOTO 53 WRITE(6,33)NE,NNE 33 FORMAT(' NE ARGUMENT(',I4,') LESS THAN NE INPUT(',I4,')') STOP 43 WRITE(ISTR,1)NP,NE NNE=NE 53 NE=NNE CALL RS4(ISTR,EE,NP*NP*NNE) CALL IS2(ISTR,IW,2*NNE) CALL RS4(ISTR,VEC,3*NNE) RETURN END SUBROUTINE IS2(IO,I2,NI2) INTEGER*2 I2(NI2) IF(IO.GT.0)GOTO 2 IIO=-IO READ(IIO,1)I2 1 FORMAT(40A2) RETURN 2 WRITE(IO,1)I2 RETURN END SUBROUTINE RS4(IO,R4,NR4) DIMENSION R4(NR4) IF(IO.GT.0)GOTO 2 IIO=-IO READ(IIO,1)R4 1 FORMAT(20A4) RETURN 2 WRITE(IO,1)R4 RETURN END ./ ENDUP ./ ADD NAME=CFTREC DOS FOR A SMALL FCC CLUSTER CONTINUED FRACTION COEFFICIENTS A(I),B2(I),I=1, 11 0.000000E+00 0.100000E+01 -0.736283E-02 0.295643E-02 -0.296524E-01 0.178592E-02 -0.126547E-01 0.221181E-02 -0.225989E-01 0.178256E-02 0.296610E-02 0.141733E-02 -0.129493E-01 0.196434E-02 -0.136124E-01 0.141376E-02 -0.170328E-01 0.181703E-02 -0.161680E-01 0.173103E-02 0.000000E+00 0.148088E-02 CONTINUED FRACTION COEFFICIENTS A(I),B2(I),I=1, 11 0.000000E+00 0.100000E+01 -0.736282E-02 0.295643E-02 -0.302165E-01 0.228776E-02 -0.115333E-01 0.194118E-02 -0.137175E-01 0.151184E-02 -0.177218E-01 0.230826E-02 -0.111183E-01 0.203490E-02 -0.575137E-02 0.141818E-02 -0.307770E-02 0.160240E-02 -0.603465E-02 0.169244E-02 0.000000E+00 0.156978E-02 CONTINUED FRACTION COEFFICIENTS A(I),B2(I),I=1, 11 0.000000E+00 0.100000E+01 -0.736282E-02 0.295643E-02 -0.296524E-01 0.178592E-02 -0.148000E-01 0.252397E-02 -0.189815E-01 0.175460E-02 0.333550E-02 0.158708E-02 -0.315587E-02 0.172282E-02 -0.671181E-02 0.106508E-02 -0.387804E-02 0.140534E-02 -0.193543E-01 0.140736E-02 0.000000E+00 0.152908E-02 CONTINUED FRACTION COEFFICIENTS A(I),B2(I),I=1, 11 0.000000E+00 0.100000E+01 -0.262667E-02 0.211155E-02 -0.871496E-02 0.128841E-02 -0.964703E-02 0.142725E-02 -0.956568E-02 0.948351E-03 -0.577315E-02 0.152401E-02 -0.187157E-01 0.137076E-02 -0.262828E-01 0.283629E-02 -0.914314E-02 0.209163E-02 -0.512340E-02 0.158929E-02 0.000000E+00 0.176099E-02 CONTINUED FRACTION COEFFICIENTS A(I),B2(I),I=1, 11 0.000000E+00 0.100000E+01 -0.262666E-02 0.211155E-02 -0.996707E-02 0.129629E-02 -0.700937E-02 0.136906E-02 -0.647660E-02 0.971641E-03 -0.848719E-02 0.122450E-02 -0.101507E-01 0.130648E-02 -0.329977E-01 0.210095E-02 -0.129033E-01 0.254819E-02 0.459147E-03 0.177905E-02 0.000000E+00 0.152031E-02 ./ ADD NAME=EXRECAL RUNFILE FROM %H=ENDJOB PRINT=SAVE FTVSCLR PROGRAM=%H=ENDPROG DD=FT01F001=.CFTREC/W/FB +++ LIBRARY=SSTLIB.TCMLIB.COM+SSTLIB.RECLIB+PUBLIC.CAMLIB INTEGER*2 MM,NN,IZERO,IZP,IW DIMENSION CRD(3,100),VEC(3,20),DUM(3),PSI(5,100),PMN(5,100),A(34) 1,B2(34),IZP(100),IW(2,20) COMMON /BLKREC/NN(100,15),MM(100,15),NAT,NP,EE(5,5,20),IZERO(100) EQUIVALENCE (IZP(1),PMN(1,1)) EQUIVALENCE (PSI(1,1),CRD(1,1)) COMMON VEC,IW EXTERNAL SLKODE,HOP,FCCBND,EQUIV C C THIS PROGRAM COMPUTES THE RECURSION COEFFICIENTS FOR THE LOCAL D.O.S. C FOR THE 5 D-ORBITALS ON A CENTRAL ATOM OF A SMALL F.C.C. CLUSTER OF C ATOMS. THIS SHOULD PROVIDE A 'BLUEPRINT' FOR OTHER APPLICATIONS OF THE C RECURSION LIBRARY; IN PARTICULAR THE ROUTINES HOP AND SET NEED VERY C LITTLE MODIFICATION TO BE APPROPRIATE FOR MANY ELECTRONIC D.O.S C CALCULATIONS, AND CHANGES SHOULD BE ABLE TO BE RESTRICTED TO THE C DIMENSION AND OUTPUT STATEMENTS.AS THE ARRAYS PSI AND PMN ARE REFERRED C TO WITH A SINGLE SUBSCRIPT IN RECAL, THEIR FIRST DIMENSION C SHOULD EQUAL THE NUMBER OF ORBITALS PER SITE, OTHERWISE THE BASIC C PROGRAM IS SUITABLE FOR ANY NUMBER OF ORBITALS PER SITE. C IN THIS EXAMPLE THE ARRAYS CRD AND IZP ARE NOT NEEDED ONCE THE C HAMILTONIAN HAS BEEN SET UP IN /BLKREC/ AND THEY ARE THEREFORE C EQUIVALENCED TO PSI AND PMN TO SAVE SPACE. THE 'TYPE' ARRAY IZP IS C PROVIDED TO ENABLE INTERACTIONS BETWEEN GEOMETRICALLY EQUIVALENT C SITES TO BE DIFFERENT, AS MAY BE REQUIRED IN LATTICES WITH TWO OR MORE C DIFFERENT ATOM TYPES,OR IN ALLOWING SURFACE INTERACTIONS TO BE C DIFFERENT FROM THOSE IN THE BULK. C LL=11 NED=20 NNDIM=100 NNMX=15 NP=5 NTYPE=1 NX=5 NY=5 NZ=6 C C LL IS THE REQUIRED LENGTH OF THE CONTINUED FRACTION C NED IS THE DIMENSION OF ARRAYS IW,VEC AND EE C NNDIM IS THE SECOND DIMENSION OF ARRAYS NN & MM AND PSI,PMN & CRD C NNMX IS THE MAXIMUM NUMBER OF NEIGHBOUR INTERACTIONS C NP IS THE NUMBER OR ORBITALS PER SITE AND FIRST 2 DIMENSIONS OF EE C NTYPE IS THE NUMBER OF 'TYPES' OF ATOM SITE C FCCLAT GENERATES A F.C.C. LATTICE OF SIZE NX X NY X NZ C NAT=NNDIM CALL FCCLAT(CRD,3,IZP,NAT,NX,NY,NZ,NTYPE) IF(NAT.EQ.0)STOP WRITE(6,4) 4 FORMAT(' OUTPUT FROM EXAMPLE RUN OF RECURSION LIBRARY - RECAL', 1 //20H LATTICE COORDINATES) WRITE(6,5)(I,(CRD(J,I),J=1,3),I=1,NAT) 5 FORMAT(4(I5,3F4.1)) CALL CLCK4X C C SETUP GENERATES THE COMMON BLOCK /BLKREC/ DEFINING THE C HAMILTONIAN MATRIX C C THE VARIABLES GENERATED AND PUT IN THE COMMON /BLKREC/ C C NN IS THE NEIGHBOUR MAP GENERATED BY NNCAL C MM IS THE INTERACTION MAP GENERATED BY MMCAL C NAT IS THE NO. OF ATOMS IN THE CLUSTER C NP IS THE NO. OF 'ORBITALS' ON EACH SITE C EE IS THE LIST OF SUBMATRICES C IZERO IS SUCH THAT IZERO(I)=0 IMPLIES PSI(.,I)=0.0 C C FOR A CLUSTER OF MORE THAN 100 ATOMS ONLY THE OBVIOUS DIMENSIONS C ABOVE NEED TO BE CHANGED , LIKEWISE FOR MORE THAN 15 C CONNECTED ATOMS NM=NNMX CALL SETUP(CRD,3,NAT,EQUIV,NTYPE,IZP,MM,NN,NNDIM,NM,SLKODE,FCCBND 1,FCCBND,EE,NP,NED,NE,VEC,IW) CALL CLCK3F(IT) WRITE(6,6)IT 6 FORMAT(14H TIME IN SETUP,I6) CALL OUT(6,1,1,IZP,IW,VEC,NED,NE,NAT,MM,NN,NNDIM,NM,EE,NP) C C CONSTRUCT INITIAL WAVE FUNCTION NONZERO ON ATOM 4,3,3 AND C COMPUTE CONTINUED FRACTION C DUM(1)=4.0 DUM(2)=3.0 DUM(3)=3.0 DO 41 NSTRT=1,NAT IF(CRD(1,NSTRT).EQ.DUM(1).AND.CRD(2,NSTRT).EQ.DUM(2).AND. 1 CRD(3,NSTRT).EQ.DUM(3))GOTO 43 41 CONTINUE WRITE(6,42) 42 FORMAT(' NO SUCH ATOM TO START') STOP 43 CONTINUE WRITE(6,44) NSTRT 44 FORMAT(I6,' ATOM TAKEN AS STARTING POINT'//) C C COMPUTES LOCAL D.O.S. RECURSION CFTS. FOR THE FIVE ORBITALS ON THE C STARTING ATOM. IZERO IS AN ARRAY SET UP TO AVIOD UNNECCESSARY MATRIX C MULTIPLICATIONS: IZERO(I)=0 IMPLIES PSI(.,I)=0.0 AND C THE ROUTINE HOP UPDATES THIS WITH EACH APPLICATION OF THE MATRIX. C IN THIS CASE (BECAUSE OF SYMMETRY) IT WOULD HAVE BEEN SUFFICIENT C TO COMPUTE ONLY THE 2 DISTINCT DENSITY FUNCTIONS. THE NUMERICAL C DISCREPENCIES IN THE ANSWERS ARE CAUSED BY LACK OF FULL SYMMETRY C IN THE CLUSTER AND STARTING STATES.FOR THE LOCAL D.O.S. SUMMED OVER C ALL ORBITALS ON THE STARTING ATOM , THE STARTING STATE (1,1,1,1,1) C ON THAT ATOM COULD HAVE BEEN USED, AND B2(1) SET EQUAL TO 5.0 C WRITE(1,10) 10 FORMAT(' DOS FOR A SMALL FCC CLUSTER') DO 15 III=1,NP CALL CLCK4X DO 11 I=1,NAT IZERO(I)=0 DO 11 J=1,NP PSI(J,I)=0.0 11 PMN(J,I)=0.0 B2(1)=1.0 A(LL)=0.0 PSI(III,NSTRT)=1.0 IZERO(NSTRT)=1 CALL RECAL(HOP,PSI,PMN,NAT*NP,A,B2,LL) WRITE(1,12)LL 12 FORMAT(5X,48H CONTINUED FRACTION COEFFICIENTS A(I),B2(I),I=1,,I3) WRITE(1,13)(A(I),B2(I),I=1,LL) 13 FORMAT(5X,2E14.6) CALL CLCK3F(IT) WRITE(6,14)IT 14 FORMAT(35H TIME TO COMPUTE TRIDIAGONALISATION,I6) 15 CONTINUE STOP END SUBROUTINE HOP(PSI,PMN,AN) INTEGER*2 MM,NN,IZERO,IDUM DOUBLE PRECISION SUM,DUM DIMENSION PSI(5,100),PMN(5,100),DUM(5) COMMON /BLKREC/NN(100,15),MM(100,15),NAT,NP,EE(5,5,20),IZERO(100) COMMON IDUM(100) C C HOP SATISFIES THE SPEC. IN RECAL. UNLABELLED COMMON IS USED FOR C AN INTEGER*2 WORK ARRAY.ALL INFORMATION FOR THE MATRIX DEFINITION C IS HELD IN THE COMMON BLOCK /BLKREC/ WHICH IS ASSIGNED IN SET C WHERE A DESCRIPTION MAY ALSO BE FOUND. C IZERO(I)=0 INDICATES PSI(.,I)=0.0 AND HENCE THAT CERTAIN LINES C MAY BE SKIPPED TO SAVE TIME. C SUM=0.0 DO 4 I=1,NAT C C THE MATRIX MULTIPLICATION PROCEEDS A LINE AT A TIME BY COMPUTING C AT EACH SITE I THE EFFECT OF THE HAMILTONIAN OPERATING AT EACH C OF ITS NEIGHBOURING SITES, WHICH ARE DEFINED BY NN(I,.) C THE PRODUCT VECTOR IS ACCUMULATED IN DUM AND THE INNER PRODUCT C IN SUM. C IDUM(I)=IZERO(I) DO 11 L=1,NP 11 DUM(L)=0.0 JSZ=NN(I,1) IE=MM(I,1) IF(IZERO(I).EQ.0)GOTO 10 C C MULTIPLY PSI(I) BY THE 'SELF ENERGY' MATRIX (OR DIAGONAL SUBMATRIX) C DO 1 M=1,NP DO 1 L=1,NP 1 DUM(L)=DUM(L)+EE(L,M,IE)*PSI(M,I) 10 IF(JSZ.LT.2)GOTO 3 DO 12 J=2,JSZ JJ=NN(I,J) IF(IZERO(JJ).EQ.0)GOTO 12 IE=MM(I,J) C C MULTIPLY PSI(.,JJ) BY THE APPROPRIATE SUBMATRIX TO CALCULATE THE C EFFECT AT SITE I OF THE HAMILTONIAN OPERATING AT SITE JJ.NOTE C THE PSI(.,I) IS NOW NON-ZERO AND THIS IS RECORDED IN IDUM FOR C LATER COPYING TO IZERO C DO 2 M=1,NP DO 2 L=1,NP 2 DUM(L)=DUM(L)+EE(L,M,IE)*PSI(M,JJ) IDUM(I)=1 12 CONTINUE C C ACCUMULATE THE REQUIRED INNER PRODUCT AND OVERWRITE PMN C 3 DO 4 L=1,NP SUM=SUM+DUM(L)*PSI(L,I) 4 PMN(L,I)=DUM(L)+PMN(L,I) AN=SUM DO 14 I=1,NAT 14 IZERO(I)=IDUM(I) RETURN END ENDPROG ENDJOB ./ ADD NAME=OPRECAL OUTPUT FROM EXAMPLE RUN OF RECURSION LIBRARY - RECAL LATTICE COORDINATES 1 2.0 1.0 1.0 2 4.0 1.0 1.0 3 1.0 2.0 1.0 4 3.0 2.0 1.0 5 5.0 2.0 1.0 6 2.0 3.0 1.0 7 4.0 3.0 1.0 8 1.0 4.0 1.0 9 3.0 4.0 1.0 10 5.0 4.0 1.0 11 2.0 5.0 1.0 12 4.0 5.0 1.0 13 1.0 1.0 2.0 14 3.0 1.0 2.0 15 5.0 1.0 2.0 16 2.0 2.0 2.0 17 4.0 2.0 2.0 18 1.0 3.0 2.0 19 3.0 3.0 2.0 20 5.0 3.0 2.0 21 2.0 4.0 2.0 22 4.0 4.0 2.0 23 1.0 5.0 2.0 24 3.0 5.0 2.0 25 5.0 5.0 2.0 26 2.0 1.0 3.0 27 4.0 1.0 3.0 28 1.0 2.0 3.0 29 3.0 2.0 3.0 30 5.0 2.0 3.0 31 2.0 3.0 3.0 32 4.0 3.0 3.0 33 1.0 4.0 3.0 34 3.0 4.0 3.0 35 5.0 4.0 3.0 36 2.0 5.0 3.0 37 4.0 5.0 3.0 38 1.0 1.0 4.0 39 3.0 1.0 4.0 40 5.0 1.0 4.0 41 2.0 2.0 4.0 42 4.0 2.0 4.0 43 1.0 3.0 4.0 44 3.0 3.0 4.0 45 5.0 3.0 4.0 46 2.0 4.0 4.0 47 4.0 4.0 4.0 48 1.0 5.0 4.0 49 3.0 5.0 4.0 50 5.0 5.0 4.0 51 2.0 1.0 5.0 52 4.0 1.0 5.0 53 1.0 2.0 5.0 54 3.0 2.0 5.0 55 5.0 2.0 5.0 56 2.0 3.0 5.0 57 4.0 3.0 5.0 58 1.0 4.0 5.0 59 3.0 4.0 5.0 60 5.0 4.0 5.0 61 2.0 5.0 5.0 62 4.0 5.0 5.0 63 1.0 1.0 6.0 64 3.0 1.0 6.0 65 5.0 1.0 6.0 66 2.0 2.0 6.0 67 4.0 2.0 6.0 68 1.0 3.0 6.0 69 3.0 3.0 6.0 70 5.0 3.0 6.0 71 2.0 4.0 6.0 72 4.0 4.0 6.0 73 1.0 5.0 6.0 74 3.0 5.0 6.0 75 5.0 5.0 6.0 TIME IN SETUP 10 NEAREST NEIGHBOUR MAP ATOM TYPE CONNECTIVITY NEIGHBOURS 1 1 6 3 4 13 14 16 2 1 6 4 5 14 15 17 3 1 6 1 6 13 16 18 4 1 9 1 2 6 7 14 16 17 19 5 1 6 2 7 15 17 20 6 1 9 3 4 8 9 16 18 19 21 7 1 9 4 5 9 10 17 19 20 22 8 1 6 6 11 18 21 23 9 1 9 6 7 11 12 19 21 22 24 10 1 6 7 12 20 22 25 11 1 6 8 9 21 23 24 12 1 6 9 10 22 24 25 13 1 6 1 3 16 26 28 14 1 9 1 2 4 16 17 26 27 29 15 1 6 2 5 17 27 30 16 1 13 1 3 4 6 13 14 18 19 26 28 29 31 17 1 13 2 4 5 7 14 15 19 20 27 29 30 32 18 1 9 3 6 8 16 21 28 31 33 19 1 13 4 6 7 9 16 17 21 22 29 31 32 34 20 1 9 5 7 10 17 22 30 32 35 21 1 13 6 8 9 11 18 19 23 24 31 33 34 36 22 1 13 7 9 10 12 19 20 24 25 32 34 35 37 23 1 6 8 11 21 33 36 24 1 9 9 11 12 21 22 34 36 37 25 1 6 10 12 22 35 37 26 1 9 13 14 16 28 29 38 39 41 27 1 9 14 15 17 29 30 39 40 42 28 1 9 13 16 18 26 31 38 41 43 29 1 13 14 16 17 19 26 27 31 32 39 41 42 44 30 1 9 15 17 20 27 32 40 42 45 31 1 13 16 18 19 21 28 29 33 34 41 43 44 46 32 1 13 17 19 20 22 29 30 34 35 42 44 45 47 33 1 9 18 21 23 31 36 43 46 48 34 1 13 19 21 22 24 31 32 36 37 44 46 47 49 35 1 9 20 22 25 32 37 45 47 50 36 1 9 21 23 24 33 34 46 48 49 37 1 9 22 24 25 34 35 47 49 50 38 1 6 26 28 41 51 53 39 1 9 26 27 29 41 42 51 52 54 40 1 6 27 30 42 52 55 41 1 13 26 28 29 31 38 39 43 44 51 53 54 56 42 1 13 27 29 30 32 39 40 44 45 52 54 55 57 43 1 9 28 31 33 41 46 53 56 58 44 1 13 29 31 32 34 41 42 46 47 54 56 57 59 45 1 9 30 32 35 42 47 55 57 60 46 1 13 31 33 34 36 43 44 48 49 56 58 59 61 47 1 13 32 34 35 37 44 45 49 50 57 59 60 62 48 1 6 33 36 46 58 61 49 1 9 34 36 37 46 47 59 61 62 50 1 6 35 37 47 60 62 51 1 9 38 39 41 53 54 63 64 66 52 1 9 39 40 42 54 55 64 65 67 53 1 9 38 41 43 51 56 63 66 68 54 1 13 39 41 42 44 51 52 56 57 64 66 67 69 55 1 9 40 42 45 52 57 65 67 70 56 1 13 41 43 44 46 53 54 58 59 66 68 69 71 57 1 13 42 44 45 47 54 55 59 60 67 69 70 72 58 1 9 43 46 48 56 61 68 71 73 59 1 13 44 46 47 49 56 57 61 62 69 71 72 74 60 1 9 45 47 50 57 62 70 72 75 61 1 9 46 48 49 58 59 71 73 74 62 1 9 47 49 50 59 60 72 74 75 63 1 4 51 53 66 64 1 6 51 52 54 66 67 65 1 4 52 55 67 66 1 9 51 53 54 56 63 64 68 69 67 1 9 52 54 55 57 64 65 69 70 68 1 6 53 56 58 66 71 69 1 9 54 56 57 59 66 67 71 72 70 1 6 55 57 60 67 72 71 1 9 56 58 59 61 68 69 73 74 72 1 9 57 59 60 62 69 70 74 75 73 1 4 58 61 71 74 1 6 59 61 62 71 72 75 1 4 60 62 72 13 NON-EQUIVALENT VECTORS FOUND ATOM INDEX OF INTERACTION VECTORS 1 13 1 2 3 4 5 2 13 1 2 3 4 5 3 13 6 2 7 4 5 4 13 8 6 1 2 7 3 4 5 5 13 8 1 7 3 5 6 13 8 6 1 2 7 3 4 5 7 13 8 6 1 2 7 3 4 5 8 13 6 2 7 4 5 9 13 8 6 1 2 7 3 4 5 10 13 8 1 7 3 5 11 13 8 6 7 3 4 12 13 8 6 7 3 4 13 13 9 10 2 4 5 14 13 11 9 10 1 2 3 4 5 15 13 11 10 1 3 5 16 13 12 11 9 10 8 6 1 2 7 3 4 5 17 13 12 11 9 10 8 6 1 2 7 3 4 5 18 13 12 9 10 6 2 7 4 5 19 13 12 11 9 10 8 6 1 2 7 3 4 5 20 13 12 11 10 8 1 7 3 5 21 13 12 11 9 10 8 6 1 2 7 3 4 5 22 13 12 11 9 10 8 6 1 2 7 3 4 5 23 13 12 9 6 7 4 24 13 12 11 9 8 6 7 3 4 25 13 12 11 8 7 3 26 13 11 9 10 1 2 3 4 5 27 13 11 9 10 1 2 3 4 5 28 13 12 9 10 6 2 7 4 5 29 13 12 11 9 10 8 6 1 2 7 3 4 5 30 13 12 11 10 8 1 7 3 5 31 13 12 11 9 10 8 6 1 2 7 3 4 5 32 13 12 11 9 10 8 6 1 2 7 3 4 5 33 13 12 9 10 6 2 7 4 5 34 13 12 11 9 10 8 6 1 2 7 3 4 5 35 13 12 11 10 8 1 7 3 5 36 13 12 11 9 8 6 7 3 4 37 13 12 11 9 8 6 7 3 4 38 13 9 10 2 4 5 39 13 11 9 10 1 2 3 4 5 40 13 11 10 1 3 5 41 13 12 11 9 10 8 6 1 2 7 3 4 5 42 13 12 11 9 10 8 6 1 2 7 3 4 5 43 13 12 9 10 6 2 7 4 5 44 13 12 11 9 10 8 6 1 2 7 3 4 5 45 13 12 11 10 8 1 7 3 5 46 13 12 11 9 10 8 6 1 2 7 3 4 5 47 13 12 11 9 10 8 6 1 2 7 3 4 5 48 13 12 9 6 7 4 49 13 12 11 9 8 6 7 3 4 50 13 12 11 8 7 3 51 13 11 9 10 1 2 3 4 5 52 13 11 9 10 1 2 3 4 5 53 13 12 9 10 6 2 7 4 5 54 13 12 11 9 10 8 6 1 2 7 3 4 5 55 13 12 11 10 8 1 7 3 5 56 13 12 11 9 10 8 6 1 2 7 3 4 5 57 13 12 11 9 10 8 6 1 2 7 3 4 5 58 13 12 9 10 6 2 7 4 5 59 13 12 11 9 10 8 6 1 2 7 3 4 5 60 13 12 11 10 8 1 7 3 5 61 13 12 11 9 8 6 7 3 4 62 13 12 11 9 8 6 7 3 4 63 13 9 10 2 64 13 11 9 10 1 2 65 13 11 10 1 66 13 12 11 9 10 8 6 1 2 67 13 12 11 9 10 8 6 1 2 68 13 12 9 10 6 2 69 13 12 11 9 10 8 6 1 2 70 13 12 11 10 8 1 71 13 12 11 9 10 8 6 1 2 72 13 12 11 9 10 8 6 1 2 73 13 12 9 6 74 13 12 11 9 8 6 75 13 12 11 8 INDEX LATTICE TYPES VECTOR MATRIX 1 1 1 -1.00 1.00 0.00 -0.021 0.000 0.000 0.000 -0.011 0.000 0.005 -0.007 0.000 0.000 0.000 -0.007 0.005 0.000 0.000 0.000 0.000 0.000 0.013 0.000 -0.011 0.000 0.000 0.000 -0.008 2 1 1 1.00 1.00 0.00 -0.021 0.000 0.000 0.000 0.011 0.000 0.005 0.007 0.000 0.000 0.000 0.007 0.005 0.000 0.000 0.000 0.000 0.000 0.013 0.000 0.011 0.000 0.000 0.000 -0.008 3 1 1 -1.00 0.00 1.00 0.005 -0.007 0.000 0.000 0.000 -0.007 0.005 0.000 0.000 0.000 0.000 0.000 -0.021 0.010 0.006 0.000 0.000 0.010 -0.003 -0.009 0.000 0.000 0.006 -0.009 0.007 4 1 1 1.00 0.00 1.00 0.005 0.007 0.000 0.000 0.000 0.007 0.005 0.000 0.000 0.000 0.000 0.000 -0.021 -0.010 -0.006 0.000 0.000 -0.010 -0.003 -0.009 0.000 0.000 -0.006 -0.009 0.007 5 1 1 0.00 1.00 1.00 0.005 0.000 0.007 0.000 0.000 0.000 -0.021 0.000 0.010 -0.006 0.007 0.000 0.005 0.000 0.000 0.000 0.010 0.000 -0.003 0.009 0.000 -0.006 0.000 0.009 0.007 6 1 1 1.00 -1.00 0.00 -0.021 0.000 0.000 0.000 -0.011 0.000 0.005 -0.007 0.000 0.000 0.000 -0.007 0.005 0.000 0.000 0.000 0.000 0.000 0.013 0.000 -0.011 0.000 0.000 0.000 -0.008 7 1 1 0.00 -1.00 1.00 0.005 0.000 -0.007 0.000 0.000 0.000 -0.021 0.000 -0.010 0.006 -0.007 0.000 0.005 0.000 0.000 0.000 -0.010 0.000 -0.003 0.009 0.000 0.006 0.000 0.009 0.007 8 1 1 -1.00 -1.00 0.00 -0.021 0.000 0.000 0.000 0.011 0.000 0.005 0.007 0.000 0.000 0.000 0.007 0.005 0.000 0.000 0.000 0.000 0.000 0.013 0.000 0.011 0.000 0.000 0.000 -0.008 9 1 1 1.00 0.00 -1.00 0.005 -0.007 0.000 0.000 0.000 -0.007 0.005 0.000 0.000 0.000 0.000 0.000 -0.021 0.010 0.006 0.000 0.000 0.010 -0.003 -0.009 0.000 0.000 0.006 -0.009 0.007 10 1 1 0.00 1.00 -1.00 0.005 0.000 -0.007 0.000 0.000 0.000 -0.021 0.000 -0.010 0.006 -0.007 0.000 0.005 0.000 0.000 0.000 -0.010 0.000 -0.003 0.009 0.000 0.006 0.000 0.009 0.007 11 1 1 -1.00 0.00 -1.00 0.005 0.007 0.000 0.000 0.000 0.007 0.005 0.000 0.000 0.000 0.000 0.000 -0.021 -0.010 -0.006 0.000 0.000 -0.010 -0.003 -0.009 0.000 0.000 -0.006 -0.009 0.007 12 1 1 0.00 -1.00 -1.00 0.005 0.000 0.007 0.000 0.000 0.000 -0.021 0.000 0.010 -0.006 0.007 0.000 0.005 0.000 0.000 0.000 0.010 0.000 -0.003 0.009 0.000 -0.006 0.000 0.009 0.007 13 1 1 0.00 0.00 0.00 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 32 ATOM TAKEN AS STARTING POINT TIME TO COMPUTE TRIDIAGONALISATION 22 TIME TO COMPUTE TRIDIAGONALISATION 22 TIME TO COMPUTE TRIDIAGONALISATION 22 TIME TO COMPUTE TRIDIAGONALISATION 23 TIME TO COMPUTE TRIDIAGONALISATION 22 ./ ENDUP ./ ADD NAME=CFTREC DOS FOR A SMALL FCC CLUSTER CONTINUED FRACTION COEFFICIENTS A(I),B2(I),I=1, 11 0.000000E+00 0.100000E+01 -0.736283E-02 0.295643E-02 -0.296524E-01 0.178592E-02 -0.126547E-01 0.221181E-02 -0.225989E-01 0.178256E-02 0.296610E-02 0.141733E-02 -0.129493E-01 0.196434E-02 -0.136124E-01 0.141376E-02 -0.170328E-01 0.181703E-02 -0.161680E-01 0.173103E-02 0.000000E+00 0.148088E-02 CONTINUED FRACTION COEFFICIENTS A(I),B2(I),I=1, 11 0.000000E+00 0.100000E+01 -0.736282E-02 0.295643E-02 -0.302165E-01 0.228776E-02 -0.115333E-01 0.194118E-02 -0.137175E-01 0.151184E-02 -0.177218E-01 0.230826E-02 -0.111183E-01 0.203490E-02 -0.575137E-02 0.141818E-02 -0.307770E-02 0.160240E-02 -0.603465E-02 0.169244E-02 0.000000E+00 0.156978E-02 CONTINUED FRACTION COEFFICIENTS A(I),B2(I),I=1, 11 0.000000E+00 0.100000E+01 -0.736282E-02 0.295643E-02 -0.296524E-01 0.178592E-02 -0.148000E-01 0.252397E-02 -0.189815E-01 0.175460E-02 0.333550E-02 0.158708E-02 -0.315587E-02 0.172282E-02 -0.671181E-02 0.106508E-02 -0.387804E-02 0.140534E-02 -0.193543E-01 0.140736E-02 0.000000E+00 0.152908E-02 CONTINUED FRACTION COEFFICIENTS A(I),B2(I),I=1, 11 0.000000E+00 0.100000E+01 -0.262667E-02 0.211155E-02 -0.871496E-02 0.128841E-02 -0.964703E-02 0.142725E-02 -0.956568E-02 0.948351E-03 -0.577315E-02 0.152401E-02 -0.187157E-01 0.137076E-02 -0.262828E-01 0.283629E-02 -0.914314E-02 0.209163E-02 -0.512340E-02 0.158929E-02 0.000000E+00 0.176099E-02 CONTINUED FRACTION COEFFICIENTS A(I),B2(I),I=1, 11 0.000000E+00 0.100000E+01 -0.262666E-02 0.211155E-02 -0.996707E-02 0.129629E-02 -0.700937E-02 0.136906E-02 -0.647660E-02 0.971641E-03 -0.848719E-02 0.122450E-02 -0.101507E-01 0.130648E-02 -0.329977E-01 0.210095E-02 -0.129033E-01 0.254819E-02 0.459147E-03 0.177905E-02 0.000000E+00 0.152031E-02 ./ ADD NAME=EXPROC RUNFILE FROM %H=ENDJOB PRINT=SAVE FILE SSTLIB.CAMREC84:CFTREC TO &DAT/W/FB FTVSCLR PROGRAM=%H=ENDPROG DD=FT05F001=&DAT LIBRARY=SSTLIB.TCMLIB.DCOM IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(15,5),B2(15,5),AA(15),BB2(15),WORK(15,6) 1,ATR(15),BTR2(15),TABS(201,5),IWK(15),EB(2),IWK2(15) NPTS=51 ALP=0.5 LL=6 ELO=-0.15 EHI=0.15 EPS=5.0E-12 ACC=1.0E-4 NP=5 DE=(EHI-ELO)/FLOAT(NPTS-1) C C LL IS THE LENGTH OF THE C.F. , NP IS THE NUMBER OF DIFFERENT C.F.S C ELO-EHI IS THE RANGE OF ENERGIES OF INTEREST, C NPTS IS THE NUMBER OF POINTS TO BE USED IN PLOTTING C EPS IS THE ACCURACY REQUIRED C READ(5,21)IWK 21 FORMAT(15A4) WRITE(6,21)IWK DO 1 I=1,NP 22 FORMAT(53X,I3) READ(5,22)LLIN 1 READ(5,2)(A(J,I),B2(J,I),J=1,LLIN) 2 FORMAT(5X,2E14.6) 3 FORMAT(' INPUT CONTINUED FRACTION CFTS.',50(/1X,10F7.3)) WRITE(6,3)((A(I,J),B2(I,J),J=1,NP),I=1,LLIN) C C COMPUTE THE SUM OF THE INPUT DENSITIES AND STORE THE RESULTING C COEFFICIENTS IN AA AND BB2 C CALL RECSUM(A,B2,15,LL,NP,AA,BB2,EPS,TABS,375) WRITE(6,2) WRITE(6,12)LL WRITE(6,2)(AA(J),BB2(J),J=1,LL) C C GENERATE TERMINATOR FOR SUMMED DENSITIES C ERR=ACC NBP1=2 LTR=LL CALL TERMGN(AA,BB2,LTR,EPS,ERR,30,EDGE,WIDTH,WEIGHT,NBP1,ATR, 1 BTR2,IWK,WORK,15,TABS,201,IWK2) WRITE(6,99)EDGE,WIDTH,WEIGHT,ERR,NBP1 99 FORMAT(' EDGE',E13.5,' WIDTH',E13.5,' WEIGHT',E13.5, 1/' ERROR IN MATCH',E13.5,' CODE',I4) WRITE(6,98)(ATR(I),BTR2(I),I=1,LTR) 98 FORMAT(' TERMINATOR COEFFICIENTS :',20(/2E16.8)) C C TABULATE D.O.S. ETC. USING C.P.C. ROUTINES C WRITE(6,2) WRITE(6,5) 5 FORMAT(5X,' ENERGY',5X,' TERM.DOS',5X,'QUAD.DOS',4X, 1' INTEGRATED DOS "STRUCTURAL ENERGY"') DO 55 I=1,NPTS TABS(I,1)=ELO+FLOAT(I-1)*DE TABS(I,5)=DENCRS(TABS(I,1),AA,BB2,LL,EDGE,WIDTH,WEIGHT,1, 1 ATR,BTR2) TABS(I,2)=DENQD(TABS(I,1),TABS(I,1),AA,BB2,LL,ALP,EPS,WORK,15, 1 NQ,NE,IWK) IF(NE.LE.0)GOTO 13 TABS(I,3)=WORK(NE,2)*(ALP-1.0) TABS(I,4)=WORK(NE,1)*TABS(I,3) DO 4 J=1,NE TABS(I,3)=TABS(I,3)+WORK(J,2) 4 TABS(I,4)=TABS(I,4)+WORK(J,2)*WORK(J,1) 55 WRITE(6,6)TABS(I,1),TABS(I,5),(TABS(I,J),J=2,4) 6 FORMAT(5E14.5) C C FIND THE FERMI ENERGY FOR INTEGER NUMBERS OF ELECTRONS C THE FIRST CALL FINDS THE F.E. OF THE SUMMED DENSITIES , C AND THE SECOND FINDS IT FROM THE SUM OF THE ORIGINAL DENSITIES. C THE LATTER SHOULD BE MORE ACCURATE, GIVING SOME IDEA OF THE ACCURACY C OF THE APPROXIMATIONS BEING MADE. C WRITE(6,2) WRITE(6,11) EB(1)=TABS(1,1) DO 9 M=1,4 DUM=EB(1) EB(2)=TABS(NPTS,1) AN=FLOAT(M) IFT=60 EF=FENVAL(AN,AA,BB2,15,LL,1,ACC,EPS,EB,WORK,15,IWK,IFT) IFT=60 EB(2)=TABS(NPTS,1) EB(1)=DUM EFAC=FENVAL(AN,A,B2,15,LL,NP,ACC,EPS,EB,WORK,15,IWK,IFT) WRITE(6,10)M,EF,EFAC EB(1)=EF 9 CONTINUE 10 FORMAT(I12,6X,2E18.8) 11 FORMAT(' NO. OF ELECTRONS',5X,'FERMI ENERGY MORE ACCURATE F.E.') 12 FORMAT(' RESULTANT CONTINUED FRACTION COEFFICIENTS',/I3 1 ,' LEVELS USED') STOP 13 WRITE(6,14)NE,TABS(I,1) 14 FORMAT(' DENQD FAILED WITH NE=',I6,' AT E=',E13.5) END ENDPROG ENDJOB ./ ADD NAME=OPPROC DOS FOR A SMALL FCC CLUSTER INPUT CONTINUED FRACTION CFTS. 0.000 1.000 0.000 1.000 0.000 1.000 0.000 1.000 0.000 1.000 -0.007 0.003 -0.007 0.003 -0.007 0.003 -0.003 0.002 -0.003 0.002 -0.030 0.002 -0.030 0.002 -0.030 0.002 -0.009 0.001 -0.010 0.001 -0.013 0.002 -0.012 0.002 -0.015 0.003 -0.010 0.001 -0.007 0.001 -0.023 0.002 -0.014 0.002 -0.019 0.002 -0.010 0.001 -0.006 0.001 0.003 0.001 -0.018 0.002 0.003 0.002 -0.006 0.002 -0.008 0.001 -0.013 0.002 -0.011 0.002 -0.003 0.002 -0.019 0.001 -0.010 0.001 -0.014 0.001 -0.006 0.001 -0.007 0.001 -0.026 0.003 -0.033 0.002 -0.017 0.002 -0.003 0.002 -0.004 0.001 -0.009 0.002 -0.013 0.003 -0.016 0.002 -0.006 0.002 -0.019 0.001 -0.005 0.002 0.000 0.002 0.000 0.001 0.000 0.002 0.000 0.002 0.000 0.002 0.000 0.002 RESULTANT CONTINUED FRACTION COEFFICIENTS 6 LEVELS USED 0.555112E-17 0.500000E+01 -0.583512E-02 0.261848E-02 -0.257177E-01 0.181036E-02 -0.177830E-01 0.230713E-02 -0.188918E-01 0.186159E-02 -0.969331E-02 0.192997E-02 EDGE -0.11241E+00 WIDTH 0.19276E+00 WEIGHT 0.50000E+01 ERROR IN MATCH 0.00000E+00 CODE 2 TERMINATOR COEFFICIENTS : -0.16024797E-01 0.50000000E+01 -0.16024797E-01 0.23223222E-02 -0.16024797E-01 0.23223222E-02 -0.16024797E-01 0.23223222E-02 -0.16024797E-01 0.23223222E-02 -0.16024797E-01 0.23223222E-02 ENERGY TERM.DOS QUAD.DOS INTEGRATED DOS "STRUCTURAL ENERGY" -0.15000E+00 0.00000E+00 0.10673E-01 0.11822E-03 -0.17733E-04 -0.14400E+00 0.00000E+00 0.20041E-01 0.20710E-03 -0.29822E-04 -0.13800E+00 0.00000E+00 0.39528E-01 0.37848E-03 -0.52230E-04 -0.13200E+00 0.00000E+00 0.82659E-01 0.72742E-03 -0.96020E-04 -0.12600E+00 0.00000E+00 0.18544E+00 0.14858E-02 -0.18721E-03 -0.12000E+00 0.00000E+00 0.45241E+00 0.32685E-02 -0.39222E-03 -0.11400E+00 0.00000E+00 0.12127E+01 0.78592E-02 -0.89594E-03 -0.10800E+00 0.16997E+01 0.35003E+01 0.20770E-01 -0.22432E-02 -0.10200E+00 0.90473E+01 0.92274E+01 0.57069E-01 -0.58210E-02 -0.96000E-01 0.11413E+02 0.12423E+02 0.12790E+00 -0.12278E-01 -0.90000E-01 0.93912E+01 0.10439E+02 0.19135E+00 -0.17505E-01 -0.84000E-01 0.94906E+01 0.12650E+02 0.26117E+00 -0.23244E-01 -0.78000E-01 0.11541E+02 0.15174E+02 0.34367E+00 -0.29489E-01 -0.72000E-01 0.15905E+02 0.20462E+02 0.44897E+00 -0.36597E-01 -0.66000E-01 0.23220E+02 0.28024E+02 0.59415E+00 -0.45302E-01 -0.60000E-01 0.32610E+02 0.32069E+02 0.77836E+00 -0.55023E-01 -0.54000E-01 0.39792E+02 0.29845E+02 0.96419E+00 -0.63757E-01 -0.48000E-01 0.41375E+02 0.32371E+02 0.11489E+01 -0.72493E-01 -0.42000E-01 0.39272E+02 0.33811E+02 0.13494E+01 -0.81598E-01 -0.36000E-01 0.36249E+02 0.32510E+02 0.15491E+01 -0.89509E-01 -0.30000E-01 0.33154E+02 0.30411E+02 0.17381E+01 -0.95696E-01 -0.24000E-01 0.29846E+02 0.27743E+02 0.19128E+01 -0.10033E+00 -0.18000E-01 0.26286E+02 0.25376E+02 0.20715E+01 -0.10389E+00 -0.12000E-01 0.22845E+02 0.24598E+02 0.22209E+01 -0.10690E+00 -0.60000E-02 0.20013E+02 0.23611E+02 0.23661E+01 -0.10921E+00 -0.13878E-16 0.18107E+02 0.21812E+02 0.25025E+01 -0.11046E+00 0.60000E-02 0.17272E+02 0.20509E+02 0.26288E+01 -0.11059E+00 0.12000E-01 0.17636E+02 0.20733E+02 0.27517E+01 -0.10961E+00 0.18000E-01 0.19439E+02 0.22581E+02 0.28810E+01 -0.10741E+00 0.24000E-01 0.23030E+02 0.24568E+02 0.30231E+01 -0.10386E+00 0.30000E-01 0.28270E+02 0.24086E+02 0.31704E+01 -0.99360E-01 0.36000E-01 0.32698E+02 0.25278E+02 0.33153E+01 -0.94837E-01 0.42000E-01 0.32088E+02 0.28896E+02 0.34793E+01 -0.89434E-01 0.48000E-01 0.27939E+02 0.29881E+02 0.36557E+01 -0.82323E-01 0.54000E-01 0.26303E+02 0.36034E+02 0.38480E+01 -0.72736E-01 0.60000E-01 0.35398E+02 0.54362E+02 0.41167E+01 -0.56849E-01 0.66000E-01 0.12190E+03 0.57940E+02 0.44102E+01 -0.38926E-01 0.72000E-01 0.18853E+02 0.44006E+02 0.48210E+01 -0.12890E-01 0.78000E-01 0.12624E+01 0.94578E+01 0.49558E+01 -0.34467E-02 0.84000E-01 0.00000E+00 0.24382E+01 0.49863E+01 -0.11528E-02 0.90000E-01 0.00000E+00 0.77956E+00 0.49949E+01 -0.46141E-03 0.96000E-01 0.00000E+00 0.29175E+00 0.49978E+01 -0.20940E-03 0.10200E+00 0.00000E+00 0.12250E+00 0.49990E+01 -0.10409E-03 0.10800E+00 0.00000E+00 0.56132E-01 0.49995E+01 -0.55459E-04 0.11400E+00 0.00000E+00 0.27553E-01 0.49997E+01 -0.31209E-04 0.12000E+00 0.00000E+00 0.14302E-01 0.49998E+01 -0.18362E-04 0.12600E+00 0.00000E+00 0.77763E-02 0.49999E+01 -0.11212E-04 0.13200E+00 0.00000E+00 0.43979E-02 0.49999E+01 -0.70656E-05 0.13800E+00 0.00000E+00 0.25730E-02 0.50000E+01 -0.45754E-05 0.14400E+00 0.00000E+00 0.15505E-02 0.50000E+01 -0.30341E-05 0.15000E+00 0.00000E+00 0.95905E-03 0.50000E+01 -0.20547E-05 NO. OF ELECTRONS FERMI ENERGY MORE ACCURATE F.E. 1 -0.52800692E-01 -0.52332914E-01 2 -0.20770662E-01 -0.21286762E-01 3 0.23056792E-01 0.23207275E-01 4 0.57708116E-01 0.57849931E-01 ./ ENDUP ./ ADD NAME=CFTPHN EXAMPLE FOR PHONON LDOS 7 LEVELS AT 2 SITES ATOM 1 COORDINATE 1 0.609375000000E+00 0.100000000000E+01 0.329112291000E+00 0.155176640000E+00 0.751473308000E+00 0.255199708000E-01 0.801091731000E+00 0.207205534000E+00 0.111061096000E+01 0.309167027000E+00 0.803737521000E+00 0.237515718000E-01 0.000000000000E+00 0.772026777000E-01 ATOM 1 COORDINATE 2 0.911381662000E+00 0.100000000000E+01 0.740966201000E+00 0.425130069000E+00 0.946415722000E+00 0.220495343000E+00 0.443695009000E+00 0.365933366000E-01 0.590427458000E+00 0.689687729000E-01 0.519128859000E+00 0.800270438000E-01 0.000000000000E+00 0.143052824000E-01 ATOM 1 COORDINATE 3 0.335937500000E+00 0.100000000000E+01 0.126945305000E+01 0.286898315000E+00 0.519608140000E+00 0.177648187000E+00 0.569977880000E+00 0.138027556000E-01 0.605309010000E+00 0.186092615000E+00 0.456537068000E+00 0.234522037000E-01 0.000000000000E+00 0.108300984000E+00 ATOM 2 COORDINATE 1 0.339285672000E+00 0.100000000000E+01 0.583071530000E+00 0.147560537000E+00 0.533853590000E+00 0.229671299000E-01 0.112207508000E+01 0.350949228000E+00 0.843207419000E+00 0.154609919000E+00 0.760353982000E+00 0.382584743000E-01 0.000000000000E+00 0.359955318000E-01 ATOM 2 COORDINATE 2 0.339285612000E+00 0.100000000000E+01 0.832189858000E+00 0.168797076000E+00 0.104482174000E+01 0.236061752000E+00 0.814646065000E+00 0.206150889000E+00 0.359184980000E+00 0.457004346000E-01 0.752830744000E+00 0.673142672000E-01 0.000000000000E+00 0.270433836000E-01 ATOM 2 COORDINATE 3 0.339285612000E+00 0.100000000000E+01 0.832197249000E+00 0.168797076000E+00 0.104480648000E+01 0.236062169000E+00 0.814657211000E+00 0.206150711000E+00 0.359177232000E+00 0.457003415000E-01 0.752839923000E+00 0.673142672000E-01 0.000000000000E+00 0.270433389000E-01 ./ ADD NAME=EXPHNDS RUNFILE FROM %H=ENDJOB FTVSCLR PROGRAM=%H=ENDFOR LIBRARY=SSTLIB.TCMLIB.COM+PUBLIC.CAMLIB +++ DATA=%H=ENDAT DD=FT09F001=.CFTPHN/W/FB DIMENSION A(50),B2(50),COORD(3,),ITYP(),TITLE(6) INTEGER*2 NN COMMON /BLKHAM/DM(3,3,),NDM COMMON /BLKMAP/NND,NNMX,NAT,NN(,5) COMMON /BLKPAR/AA,BB(2),RTMASS(2),NON COMMON /BLKWFN/PSI(3,),PMN(3,) EQUIVALENCE (COORD(1,1),PSI(1,1)),(ITYP(1),PMN(1,1)) EXTERNAL DMGEN,PRNT,DMHOP C C THIS PROGRAM GENERATES THE RECURSION COEFFICIENTS FOR THE C PHONON LOCAL DENSITY OF STATES (LDOS) OF A SMALL CLUSTER OF C ATOMS AS AN EXAMPLE OF THE METHOD THAT MAY BE USED ON MORE REALISTIC C SYSTEMS.MOST OF THE PROGRAM AND SUBROUTINES NEED ONLY SMALL C (IF ANY) CHANGES TO PERFORM THE CALCULATION FOR A WIDE RANGE OF C MATERIALS, AND THESE SHOULD BE MAINLY TO DIMENSION STATEMENTS C AND PARAMETER VALUES.THE OUTPUT MAY BE PROCESSED BY A PROGRAM C SUCH AS PHNPLT TO PRODUCE THE ACTUAL GRAPHS OF PHON LDOS AGAINST C FREQUENCY C C C THE IMPORTANT PARAMETERS OF THE PROGRAM ARE AS FOLLOWS: C C NND FIRST DIMENION OF ARRAY NN C NNMX MAXIMUM NUMBER OF NEAREST NEIGHBOURS C NON REMOTENESS OF NEIGHBOURS TO BE TAKEN ONTO ACCOUNT C I.E. =2 FOR 1 & 2 NEIGHBOURS (BOND-STRETCHING & BENDING) C =1 FOR 1ST. NEIGHBOURS ONLY C NAT NUMBER OF ATOMS IN THE CLUSTER C C SEE SUBROUTINE DMGEN FOR A DESCRIPTION OF VARIABLES IN THE C COMMON BLOCKS C NND= NNMX=5 NON=2 C C READ NUMBER OF ATOMS IN CLUSTER, THEN 'NEIGHBOUR MAP' AND ATOM C TYPES AND COORDINATES C READ(5,21)TITLE 21 FORMAT(6A4) WRITE(6,22)TITLE 22 FORMAT(20X,6A4) WRITE(9,21)TITLE READ(5,1)NAT 1 FORMAT(I4) WRITE(6,2)NAT 2 FORMAT(20X,' CLUSTER OF',I6,' ATOMS') DO 4 I=1,NAT READ(5,3)(NN(I,J),J=1,4),ITYP(I),(COORD(J,I),J=1,3) 3 FORMAT(5I4,3F6.2) 4 CONTINUE DO 6 I=1,NAT NJ=NN(I,1) WRITE(6,5)I,ITYP(I),(COORD(J,I),J=1,3),(NN(I,J),J=1,NJ) 5 FORMAT(20X,2I4,3F8.2,4I4) 6 CONTINUE C C GENERATE DYNAMICAL MATRICES USING SCAN AND DMGEN C AA IS BOND STRETCHING FORCE CONSTANT C BB(J) IS BOND BENDING FORCE CONSTANT AT ATOM OF TYPE J C RTMASS(J) IS THE ROOT INVERSE MASS OF AN ATOM OF TYPE J C AA=1.0 BB(1)=0.5 BB(2)=0.25 RTMASS(1)=0.25 RTMASS(2)=SQRT(1.0/28.0) NDM=0 CALL CLCK4X CALL SCAN(NN,NND,NNMX,1,NAT,NON,DMGEN) CALL CLCK3F(ITIM) WRITE(6,7)NDM 7 FORMAT(13X,I8,' DYNAMICAL MATRICES CALCULATED') WRITE(6,8)ITIM 8 FORMAT(1H+,51X,'IN',I8,' 1/100 SECS') C C PRINT OUT DYNAMICAL MATRICES C NDM=0 CALL SCAN(NN,NND,NNMX,1,NAT,NON,PRNT) NSITES=2 NDISP=3 LL=7 WRITE(9,23)LL,NSITES 23 FORMAT(I4,' LEVELS AT',I4,' SITES') WRITE(6,11)LL 11 FORMAT(13X,I4,' LEVELS') C C GENERATE LOCAL DENSITY OF PHONON MODES AT THE FIRST NSITES C ATOMS FOR NDISP CARTESIAN COORDINATES, WITH LL 'LEVELS' EACH. C DO 13 II=1,NSITES DO 13 JJ=1,NDISP CALL CLCK4X DO 12 I=1,NAT DO 12 J=1,3 PSI(J,I)=0.0 12 PMN(J,I)=0.0 B2(1)=1.0 PSI(JJ,II)=1.0 CALL RECAL(DMHOP,PSI,PMN,3*NAT,A,B2,LL) CALL CLCK3F(ITIM) WRITE(9,24)II,JJ WRITE(9,25)(A(I),B2(I),I=1,LL) 24 FORMAT(' ATOM',I5,' COORDINATE',I2) 25 FORMAT(2E20.12) WRITE(6,14)II,JJ,(A(I),B2(I),I=1,LL) WRITE(6,8)ITIM 13 CONTINUE 14 FORMAT(20X,' ATOM',I5,' COORDINATE',I3,50(/20X,2E13.5)) STOP END SUBROUTINE VECGEN(IA,NA,NOP) DIMENSION IA(NOP),NA(NOP),V(3) COMMON /BLKWFN/COORD(3,),ITYP(),JUNK(2,) COMMON /BLKLOC/AB(3,5,5),NS C C SATISFIES THE SPEC. OF SUB IN SCAN, AND CALCULATES AND STORES IN THE C COMMON BLOCK /BLKLOC/ THE NORMALISED VECTOR BETWEEN ATOMS IA(NOP-1) C AND IA(NOP). THE COORDINATES OF THE ATOMS HAVE BEEN STORED IN COORD C AND THE VECTORS GENERATED ARE PUT IN THE ARRAY AB WITH THE LAST TWO C SUBSCRIPTS REFERENCING THE NEIGHBOUR INDEXES NA(NOP-1),NA(NOP) AND C THE FIRST SUBSCRIPT REFERRING TO THE COORDINATE. IF NOP=1 NO ACTION C IS TAKEN. C IF(NOP.EQ.1)RETURN I=IA(NOP-1) J=IA(NOP) SUM=0.0 DO 1 L=1,3 V(L)=COORD(L,J)-COORD(L,I) 1 SUM=SUM+V(L)*V(L) SUM=1.0/SQRT(SUM) I=NA(NOP-1) J=NA(NOP) DO 2 L=1,3 2 AB(L,I,J)=V(L)*SUM RETURN END SUBROUTINE DMGEN(IA,NA,NOP) DIMENSION IA(NOP),NA(NOP),A(3,3) INTEGER*2 NN COMMON /BLKWFN/COORD(3,),ITYP(),JUNK(2,) COMMON /BLKHAM/DM(3,3,),ND COMMON /BLKLOC/AB(3,5,5),NS COMMON /BLKMAP/NND,NNMX,NAT,NN(,5) COMMON /BLKPAR/AA,BB(2),RTMASS(2),NO EXTERNAL VECGEN C C SATISFIES SPEC. OF SUB IN SCAN AND COMPUTES AND ACCUMULATES THE C DYNAMICAL MATRIX INTERACTING BETWEEN ATOMS IA(1) AND IA(NOP) C ACCORDING TO THE FORMULATION IN THE THESIS OF P.E.MEEK .THE C SELF-INTERACTION IS ACCUMULATED AS -1 * THE SUM OF ALL NEIGHBOUR C INTERACTIONS . THE MATRICES ARE ONLY STORED IF IA(NOP).GE.IA(1), C THOSE WHICH ARE STORED ARE PUT, IN THE ORDER GENERATED BY SCAN, C INTO THE ARRAY DM. C SCAN1 IS CALLED IF NOP=1 TO GENERATE THE NORMALISED VECTORS C NEEDED FOR THE CALCULATION OF ALL INTERACTIONS OF THE CURRENT C ATOM. C C VARIABLES IN THE COMMON BLOCKS: C C COORD LIST OF COORDINATES OF THE ATOM CLUSTER C ITYP LIST OF TYPES OF ATOMS C JUNK ARRAY TO MAKE /BLKWFN/ THE CORRECT LENGTH C C DM LIST OF 3 X 3 DYNAMICAL MATRICES C ND NUMBER OF MATRICES GENERATED SO FAR.( INITIALISED IN MAIN) C C AB LIST OF VECTORS GENERATED BY VECGEN CALLED FROM SCAN1 C FOR THE CURRENT 'INITIAL' ATOM C NS INDEX OF SELF-INTERACTION MATRIX FOR CURRENT 'INITIAL' C ATOM.USED IN ACCUMULATION OF THAT MATRIX. C C NND FIRST DIMENSION OF ARRAY NN C NNMX 1+ MAX. NO. OF NEIGHBOURS OF ANY ATOM C NAT NO. OF ATOMS IN CLUSTER C NN NEAREST NEIGHBOUR MAP (MAY BE DECLARED INTEGER*2 TO SAVE SPACE): C NN(I,1) = 1+NO. OF NEIGHBOURS OF ATOM I C NN(I,J),J=2,..,NN(I,1) LIST OF ATOM NUMBERS OF C NEIGHBOURS OF ATOM I C C AA BOND STRETCHING FORCE CONSTANT C BB LIST OF BOND BENDING FORCE CONSTANTS CORRESPONDING TO C EACH ATOM TYPE. C RTMASS LIST OF (MASS OF ATOM)**(-0.5) FOR EACH TYPE OF ATOM C NO =1 FOR FIRST NEIGHBOUR INTERACTIONS ONLY (BOND STRETCHING) C =2 FOR 1ST. AND 2ND.NEIGHBOUR INTERACTIONS (STRETCHIN & BENDING) C GOTO(1,10,20),NOP C C INITIALISE SELF-INTERACTION MATRIX & CONSRUCT NORMALISED VECTORS C 1 NS=ND+1 ND=NS IF(ND.GT.)GOTO 33 DO 2 L=1,3 DO 2 M=1,3 2 DM(L,M,NS)=0.0 CALL SCAN1(NN,NND,NNMX,IA(1),IA(1),NO,VECGEN) RETURN C C COMPUTE FIRST NEIGHBOUR INTERACTION C 10 N1=NA(NOP) DO 11 L=1,3 CON=-8.0*AA*AB(L,1,N1) DO 11 M=1,3 11 A(L,M)=CON*AB(M,1,N1) I=IA(1) NJ=NN(I,1) C C SUM OVER FIRST NEIGHBOURS OF ATOM IA(1) C IT=ITYP(I) DO 13 JJ=2,NJ IF(JJ.EQ.N1)GOTO 13 DO 12 L=1,3 CON=-BB(IT)*(AB(L,1,N1)+AB(L,1,JJ)) DO 12 M=1,3 12 A(L,M)=A(L,M)+CON*AB(M,1,JJ) 13 CONTINUE C C SUM OVER FIRST NEIGHBOURS OF ATOM IA(NOP) C J=IA(NOP) NK=NN(J,1) IT=ITYP(J) DO 15 KK=2,NK IF(NN(J,KK).EQ.I)GOTO 15 DO 14 L=1,3 CON=-BB(IT)*AB(L,N1,KK) DO 14 M=1,3 14 A(L,M)=A(L,M)+CON*(AB(M,N1,KK)-AB(M,1,N1)) 15 CONTINUE GOTO 30 C C COMPUTE SECOND NEIGHBOUR MATRIX C 20 J=IA(NOP-1) IT=ITYP(J) JJ=NA(NOP-1) KK=NA(NOP) DO 21 L=1,3 CON=-BB(IT)*AB(L,JJ,KK) DO 21 M=1,3 21 A(L,M)=CON*AB(M,1,JJ) C C ACCUMULATE SELF INTERACTION MATRIX AND STORE IF IN C UPPER BLOCK TRIANGLE AND SCALE ACCORDING TO APPROPRIATE MASSES C 30 I=IA(1) I=ITYP(I) J=IA(NOP) J=ITYP(J) CON=RTMASS(I)*RTMASS(I) DO 31 L=1,3 DO 31 M=1,3 31 DM(L,M,NS)=DM(L,M,NS)-A(L,M)*CON IF(IA(NOP).LT.IA(1))RETURN ND=ND+1 IF(ND.GT.)GOTO 33 CON=RTMASS(I)*RTMASS(J) DO 32 L=1,3 DO 32 M=1,3 32 DM(L,M,ND)=A(L,M)*CON RETURN 33 WRITE(6,34)ND 34 FORMAT(20X,' NOT ENOUGH ROOM FOR DYNAMICAL MATRICES',I6) STOP END SUBROUTINE PRNT(IA,NA,NOP) DIMENSION IA(NOP),NA(NOP) COMMON /BLKHAM/DM(3,3,),NDM IF(IA(NOP).LT.IA(1))RETURN NDM=NDM+1 WRITE(6,1)IA(NOP),IA(1),((DM(L,M,NDM),M=1,3),L=1,3) 1 FORMAT(20X,' MATRIX FOR RESULT OF H ON ATOM',I4,' AT ATOM',I4, 1 3(/20X,3E13.5)) RETURN END SUBROUTINE DMHOP(X,Y,A) DIMENSION X(3,),Y(3,) INTEGER*2 NN COMMON /BLKA/AA COMMON /BLKHAM/DM(3,3,),NDM COMMON /BLKMAP/NND,NNMX,NAT,NN(,5) COMMON /BLKPAR/AAA,BB(2),RTMASS(2),NON EXTERNAL MLT C C SATISFIES SPEC. OF HOP IN RECAL, CALLING SCAN TO PERFORM C THE OPERATIONS WITH THE ROUTINE MLT, HAVING INITIALISED AA C AND NDM . THE VECTORS REQUIRED ARE TRANSFERED VIA C THE COMMON BLOCK /BLKWFN/ INSTEAD OF THE ARGUMENTS. THIS C IS POSSIBLE BECAUSE OF THE COMMON STATEMENT IN THE MAIN C PROGRAM AND THE STRUCTURE OF RECAL. C AA=0.0 NDM=0 CALL SCAN(NN,NND,NNMX,1,NAT,NON,MLT) A=AA RETURN END SUBROUTINE MLT(IA,NA,NOP) DIMENSION IA(NOP),NA(NOP),DUM(3) COMMON /BLKA/AA COMMON /BLKHAM/DM(3,3,),NDM COMMON /BLKWFN/PSI(3,),PMN(3,) C C SATISFIES THE SPEC. OF SUB IN SCAN, COMPUTING THE PRODUCT C OF THE DYNAMICAL MATRIX ,H, (CONSTRUCTED IN DMGEN) AND A C VECTOR ,PSI, STORED IN COMMON /BLKWFN/ .THE RESULT OF C PMN + H. PSI IS OVERWRITTEN ONTO PMN AND PSI.H.PSI ONTO AA C N=IA(1) NN=IA(NOP) IF(NN.LT.N)RETURN C C COMPUTE THE EFFECT AT ATOM N OF H OPERATING ON ATOM NN C NDM=NDM+1 DO 1 L=1,3 1 DUM(L)=0.0 DO 2 M=1,3 DO 2 L=1,3 2 DUM(L)=DUM(L)+DM(L,M,NDM)*PSI(M,NN) DO 3 L=1,3 AA=AA+DUM(L)*PSI(L,N) 3 PMN(L,N)=PMN(L,N)+DUM(L) IF(N.EQ.NN)RETURN C C COMPUTE THE EFFECT AT ATOM NN OF H OPERATING ON ATOM N C DO 11 L=1,3 11 DUM(L)=0.0 DO 12 M=1,3 DO 12 L=1,3 12 DUM(L)=DUM(L)+DM(M,L,NDM)*PSI(M,N) DO 13 L=1,3 AA=AA+DUM(L)*PSI(L,NN) 13 PMN(L,NN)=PMN(L,NN)+DUM(L) RETURN END ENDFOR EXAMPLE FOR PHONON LDOS 6 4 2 3 6 1 0.00 0.00 0.00 4 1 3 4 2 0.00 1.00 0.0 4 1 2 5 1 0.0 1.0 1.0 4 2 5 6 1 1.0 1.0 0.0 4 3 4 6 2 1.0 1.0 1.0 4 1 4 5 2 1.0 0.0 0.0 ENDAT ENDJOB ./ ADD NAME=OPPHNDS EXAMPLE FOR PHONON LDOS CLUSTER OF 6 ATOMS 1 1 0.00 0.00 0.00 4 2 3 6 2 2 0.00 1.00 0.00 4 1 3 4 3 1 0.00 1.00 1.00 4 1 2 5 4 1 1.00 1.00 0.00 4 2 5 6 5 2 1.00 1.00 1.00 4 3 4 6 6 2 1.00 0.00 0.00 4 1 4 5 33 DYNAMICAL MATRICES CALCULATED IN 0 1/100 SECS MATRIX FOR RESULT OF H ON ATOM 1 AT ATOM 1 0.60938E+00 0.53347E-01 0.22097E-01 0.53347E-01 0.91138E+00 0.31116E+00 0.22097E-01 0.31116E+00 0.33594E+00 MATRIX FOR RESULT OF H ON ATOM 2 AT ATOM 1 -0.35434E-01 0.11811E-01 0.00000E+00 -0.23623E-01 -0.40648E+00 -0.28515E-01 0.00000E+00 -0.28161E-08 -0.23623E-01 MATRIX FOR RESULT OF H ON ATOM 3 AT ATOM 1 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 -0.15625E-01 0.00000E+00 MATRIX FOR RESULT OF H ON ATOM 4 AT ATOM 1 0.00000E+00 -0.15625E-01 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 MATRIX FOR RESULT OF H ON ATOM 3 AT ATOM 1 -0.62500E-01 0.22097E-01 0.22097E-01 -0.22097E-01 -0.30335E+00 -0.25000E+00 -0.22097E-01 -0.29419E+00 -0.30335E+00 MATRIX FOR RESULT OF H ON ATOM 2 AT ATOM 1 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.16704E-01 0.16704E-01 MATRIX FOR RESULT OF H ON ATOM 5 AT ATOM 1 0.00000E+00 -0.16704E-01 -0.16704E-01 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 MATRIX FOR RESULT OF H ON ATOM 6 AT ATOM 1 -0.37796E+00 -0.40327E-01 -0.16704E-01 0.20163E-01 -0.53151E-01 -0.17717E-01 0.83519E-02 -0.17717E-01 -0.17717E-01 MATRIX FOR RESULT OF H ON ATOM 4 AT ATOM 1 0.00000E+00 0.00000E+00 0.00000E+00 -0.15625E-01 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 MATRIX FOR RESULT OF H ON ATOM 5 AT ATOM 1 0.00000E+00 0.00000E+00 0.00000E+00 -0.83519E-02 0.00000E+00 0.00000E+00 -0.83519E-02 0.00000E+00 0.00000E+00 MATRIX FOR RESULT OF H ON ATOM 2 AT ATOM 2 0.33929E+00 -0.89286E-02 0.89286E-02 -0.89286E-02 0.33929E+00 0.89286E-02 0.89286E-02 0.89286E-02 0.33929E+00 MATRIX FOR RESULT OF H ON ATOM 3 AT ATOM 2 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.16704E-01 0.00000E+00 0.00000E+00 0.16704E-01 0.00000E+00 MATRIX FOR RESULT OF H ON ATOM 6 AT ATOM 2 0.00000E+00 0.17857E-01 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 MATRIX FOR RESULT OF H ON ATOM 3 AT ATOM 2 -0.35434E-01 0.00000E+00 0.23623E-01 0.00000E+00 -0.23623E-01 -0.28515E-01 -0.11811E-01 -0.28161E-08 -0.40648E+00 MATRIX FOR RESULT OF H ON ATOM 5 AT ATOM 2 0.00000E+00 0.00000E+00 -0.17857E-01 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 MATRIX FOR RESULT OF H ON ATOM 4 AT ATOM 2 -0.37796E+00 0.11811E-01 -0.11811E-01 -0.23623E-01 -0.35434E-01 0.00000E+00 0.23623E-01 0.00000E+00 -0.35434E-01 MATRIX FOR RESULT OF H ON ATOM 5 AT ATOM 2 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 -0.17857E-01 0.00000E+00 0.00000E+00 MATRIX FOR RESULT OF H ON ATOM 6 AT ATOM 2 0.00000E+00 0.00000E+00 0.00000E+00 0.17857E-01 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 MATRIX FOR RESULT OF H ON ATOM 3 AT ATOM 3 0.60938E+00 -0.22097E-01 -0.53347E-01 -0.22097E-01 0.33594E+00 0.31116E+00 -0.53347E-01 0.31116E+00 0.91138E+00 MATRIX FOR RESULT OF H ON ATOM 6 AT ATOM 3 0.00000E+00 0.16704E-01 0.16704E-01 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 MATRIX FOR RESULT OF H ON ATOM 4 AT ATOM 3 0.00000E+00 0.00000E+00 0.15625E-01 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 MATRIX FOR RESULT OF H ON ATOM 5 AT ATOM 3 -0.37796E+00 0.16704E-01 0.40327E-01 -0.83519E-02 -0.17717E-01 -0.17717E-01 -0.20163E-01 -0.17717E-01 -0.53151E-01 MATRIX FOR RESULT OF H ON ATOM 4 AT ATOM 3 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.15625E-01 0.00000E+00 0.00000E+00 MATRIX FOR RESULT OF H ON ATOM 6 AT ATOM 3 0.00000E+00 0.00000E+00 0.00000E+00 0.83519E-02 0.00000E+00 0.00000E+00 0.83519E-02 0.00000E+00 0.00000E+00 MATRIX FOR RESULT OF H ON ATOM 4 AT ATOM 4 0.59375E+00 0.31250E-01 -0.31250E-01 0.31250E-01 0.59375E+00 -0.15625E-01 -0.31250E-01 -0.15625E-01 0.59375E+00 MATRIX FOR RESULT OF H ON ATOM 5 AT ATOM 4 -0.35434E-01 0.00000E+00 -0.11811E-01 0.00000E+00 -0.29528E-01 -0.14258E-01 0.23623E-01 0.17717E-01 -0.39222E+00 MATRIX FOR RESULT OF H ON ATOM 6 AT ATOM 4 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.83519E-02 0.00000E+00 0.00000E+00 0.83519E-02 MATRIX FOR RESULT OF H ON ATOM 6 AT ATOM 4 -0.35434E-01 0.11811E-01 0.00000E+00 -0.23623E-01 -0.39222E+00 0.17717E-01 0.00000E+00 -0.14258E-01 -0.29528E-01 MATRIX FOR RESULT OF H ON ATOM 5 AT ATOM 4 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.83519E-02 0.00000E+00 0.00000E+00 0.83519E-02 0.00000E+00 MATRIX FOR RESULT OF H ON ATOM 5 AT ATOM 5 0.33036E+00 0.63135E-02 0.15242E-01 0.63135E-02 0.18750E+00 0.16703E+00 0.15242E-01 0.16703E+00 0.49477E+00 MATRIX FOR RESULT OF H ON ATOM 6 AT ATOM 5 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 -0.17857E-01 0.00000E+00 0.00000E+00 0.00000E+00 MATRIX FOR RESULT OF H ON ATOM 6 AT ATOM 5 -0.17857E-01 0.63135E-02 0.63135E-02 -0.63135E-02 -0.15810E+00 -0.15548E+00 -0.63135E-02 -0.14286E+00 -0.15810E+00 MATRIX FOR RESULT OF H ON ATOM 6 AT ATOM 6 0.33036E+00 -0.15242E-01 -0.63135E-02 -0.15242E-01 0.49477E+00 0.16703E+00 -0.63135E-02 0.16703E+00 0.18750E+00 7 LEVELS ATOM 1 COORDINATE 1 0.60938E+00 0.10000E+01 0.32911E+00 0.15518E+00 0.75147E+00 0.25520E-01 0.80109E+00 0.20721E+00 0.11106E+01 0.30917E+00 0.80374E+00 0.23752E-01 0.00000E+00 0.77203E-01 IN 1 1/100 SECS ATOM 1 COORDINATE 2 0.91138E+00 0.10000E+01 0.74097E+00 0.42513E+00 0.94642E+00 0.22050E+00 0.44370E+00 0.36593E-01 0.59043E+00 0.68969E-01 0.51913E+00 0.80027E-01 0.00000E+00 0.14305E-01 IN 1 1/100 SECS ATOM 1 COORDINATE 3 0.33594E+00 0.10000E+01 0.12695E+01 0.28690E+00 0.51961E+00 0.17765E+00 0.56998E+00 0.13803E-01 0.60531E+00 0.18609E+00 0.45654E+00 0.23452E-01 0.00000E+00 0.10830E+00 IN 1 1/100 SECS ATOM 2 COORDINATE 1 0.33929E+00 0.10000E+01 0.58307E+00 0.14756E+00 0.53385E+00 0.22967E-01 0.11221E+01 0.35095E+00 0.84321E+00 0.15461E+00 0.76035E+00 0.38258E-01 0.00000E+00 0.35996E-01 IN 1 1/100 SECS ATOM 2 COORDINATE 2 0.33929E+00 0.10000E+01 0.83219E+00 0.16880E+00 0.10448E+01 0.23606E+00 0.81465E+00 0.20615E+00 0.35918E+00 0.45700E-01 0.75283E+00 0.67314E-01 0.00000E+00 0.27043E-01 IN 1 1/100 SECS ATOM 2 COORDINATE 3 0.33929E+00 0.10000E+01 0.83220E+00 0.16880E+00 0.10448E+01 0.23606E+00 0.81466E+00 0.20615E+00 0.35918E+00 0.45700E-01 0.75284E+00 0.67314E-01 0.00000E+00 0.27043E-01 IN 1 1/100 SECS ./ ENDUP ./ ADD NAME=CFTPHN EXAMPLE FOR PHONON LDOS 7 LEVELS AT 2 SITES ATOM 1 COORDINATE 1 0.609375000000E+00 0.100000000000E+01 0.329112291000E+00 0.155176640000E+00 0.751473308000E+00 0.255199708000E-01 0.801091731000E+00 0.207205534000E+00 0.111061096000E+01 0.309167027000E+00 0.803737521000E+00 0.237515718000E-01 0.000000000000E+00 0.772026777000E-01 ATOM 1 COORDINATE 2 0.911381662000E+00 0.100000000000E+01 0.740966201000E+00 0.425130069000E+00 0.946415722000E+00 0.220495343000E+00 0.443695009000E+00 0.365933366000E-01 0.590427458000E+00 0.689687729000E-01 0.519128859000E+00 0.800270438000E-01 0.000000000000E+00 0.143052824000E-01 ATOM 1 COORDINATE 3 0.335937500000E+00 0.100000000000E+01 0.126945305000E+01 0.286898315000E+00 0.519608140000E+00 0.177648187000E+00 0.569977880000E+00 0.138027556000E-01 0.605309010000E+00 0.186092615000E+00 0.456537068000E+00 0.234522037000E-01 0.000000000000E+00 0.108300984000E+00 ATOM 2 COORDINATE 1 0.339285672000E+00 0.100000000000E+01 0.583071530000E+00 0.147560537000E+00 0.533853590000E+00 0.229671299000E-01 0.112207508000E+01 0.350949228000E+00 0.843207419000E+00 0.154609919000E+00 0.760353982000E+00 0.382584743000E-01 0.000000000000E+00 0.359955318000E-01 ATOM 2 COORDINATE 2 0.339285612000E+00 0.100000000000E+01 0.832189858000E+00 0.168797076000E+00 0.104482174000E+01 0.236061752000E+00 0.814646065000E+00 0.206150889000E+00 0.359184980000E+00 0.457004346000E-01 0.752830744000E+00 0.673142672000E-01 0.000000000000E+00 0.270433836000E-01 ATOM 2 COORDINATE 3 0.339285612000E+00 0.100000000000E+01 0.832197249000E+00 0.168797076000E+00 0.104480648000E+01 0.236062169000E+00 0.814657211000E+00 0.206150711000E+00 0.359177232000E+00 0.457003415000E-01 0.752839923000E+00 0.673142672000E-01 0.000000000000E+00 0.270433389000E-01 ./ ADD NAME=EXPHNPLT RUNFILE FROM %H=ENDJOB FILE SSTLIB.CAMREC84:CFTPHN TO &DAT/W/FB FTVSCLR PROGRAM=%H=ENDPROG DD=FT05F001=&DAT LIBRARY=SSTLIB.TCMLIB.DCOM IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(100,18),B2(100,18),AA(100),BB2(100),WK(100,18,5) 1,TABS(511,3),IWK(100),TITL(6),NAT(6),W1(100),W2(100),W3(100) NDIM=100 IO=6 NPTS=21 WLO=0.0 WHI=1.5 ALP=0.5 EPS=5.0E-8 READ(5,11)TITL 11 FORMAT(6A4) WRITE(IO,10)TITL 10 FORMAT(10X,6A4) READ(5,12)LL,NS 12 FORMAT(I4,10X,I4) NP=NS*3 C C NDIM IS FIRST DIMENSION OF ARRAYS A AND B2 C NPTS IS NUMBER OF POINTS OF TABULATION C WLO-WHI IS THE RANGE OF FREQUENCIES OF INTEREST, C EPS IS THE ACCURACY REQUIRED C TITL HOLDS THE TITLE TO BE COPIED FROM INPUT TO OUTPUT C LL IS THE LENGTH OF THE CONTINUED FRACTIONS C NS IS THE NUMBER OF ATOMS AT WHICH THE PHONON LDOS IS CALCULATED C NP (=3*NS) IS THUS THE NUMBER OF SETS OF CONTINUED FRACTION CFTS C C READ IN COEFFICIENTS C DO 1 I=1,NP II=(I-1)/3+1 READ(5,14)NAT(II),IC 14 FORMAT(5X,I5,11X,I2) READ(5,2)(A(J,I),B2(J,I),J=1,LL) 1 CONTINUE 2 FORMAT(2E20.12) WRITE(IO,13)LL,NS 13 FORMAT(10X,I4,' LEVELS AT',I4,' SITES') WRITE(IO,23)(NAT(I),I=1,NS) 23 FORMAT(10X,' ATOMS',10I5) DO 4 JJ=1,NP,3 JM=MIN0(JJ+2,NP) WRITE(IO,3)((A(I,J),B2(I,J),J=JJ,JM),I=1,LL) 3 FORMAT(10X,' INPUT CONTINUED FRACTION CFTS.',50(/10X,6E11.4)) 4 CONTINUE C C COMPUTE THE SUM OF THE INPUT DENSITIES AND STORE THE RESULTING C COEFFICIENTS IN AA AND BB2 C CALL RECSUM(A,B2,NDIM,LL,NP,AA,BB2,EPS,WK,5*LL*NP) WRITE(IO,2) WRITE(IO,22) 22 FORMAT(10X,' RESULTANT CONTINUED FRACTION COEFFICIENTS') WRITE(IO,24)(AA(J),BB2(J),J=1,LL) 24 FORMAT(10X,2E20.12) DX=(WHI-WLO)/FLOAT(NPTS-1) WRITE(IO,2) C C TABULATE PHONON LDOS GIVEN BY 2*W*N(W*W) C WHERE N(E) CORRESPONDS TO AN ELECTRON LDOS C TABS(I,1) CONTAINS W C TABS(I,2) CONTAINS THE PHONON LOCAL DENSITY OF STATES C TABS(I,3) CONTAINS THE INTEGRAL OF THE ABOVE C WRITE(IO,5) 5 FORMAT(18X,' OMEGA',5X,' DENSITY OF PHONON MODES INTEGRATED ', 1'D.O.S.') DO 15 I=1,NPTS TABS(I,1)=WLO+FLOAT(I-1)*DX E=TABS(I,1)*TABS(I,1) TABS(I,2)=2.0*TABS(I,1)*DENQD(E,E,AA,BB2,LL,ALP,EPS,WK,NDIM, 1 NQ,NE,IWK) IF(NE.GT.0)GOTO 42 WRITE(IO,41)WK(I,1,1),NE,I 41 FORMAT(' TROUBLE AT E=',' NQ=',I3,' I=',I3) STOP 42 TABS(I,3)=WK(NE,2,1)*(ALP-1.0) DO 43 J=1,NE 43 TABS(I,3)=TABS(I,3)+WK(J,2,1) TABS(I,3)=TABS(I,3)*2.0 WRITE(IO,6)(TABS(I,J),J=1,3) 15 CONTINUE 6 FORMAT(10X,3E20.8) STOP END ENDPROG ENDJOB ./ ADD NAME=OPPHNPLT EXAMPLE FOR PHONON LDOS 7 LEVELS AT 2 SITES ATOMS 1 2 INPUT CONTINUED FRACTION CFTS. 0.6094D+00 0.1000D+01 0.9114D+00 0.1000D+01 0.3359D+00 0.1000D+01 0.3291D+00 0.1552D+00 0.7410D+00 0.4251D+00 0.1269D+01 0.2869D+00 0.7515D+00 0.2552D-01 0.9464D+00 0.2205D+00 0.5196D+00 0.1776D+00 0.8011D+00 0.2072D+00 0.4437D+00 0.3659D-01 0.5700D+00 0.1380D-01 0.1111D+01 0.3092D+00 0.5904D+00 0.6897D-01 0.6053D+00 0.1861D+00 0.8037D+00 0.2375D-01 0.5191D+00 0.8003D-01 0.4565D+00 0.2345D-01 0.0000D+00 0.7720D-01 0.0000D+00 0.1431D-01 0.0000D+00 0.1083D+00 INPUT CONTINUED FRACTION CFTS. 0.3393D+00 0.1000D+01 0.3393D+00 0.1000D+01 0.3393D+00 0.1000D+01 0.5831D+00 0.1476D+00 0.8322D+00 0.1688D+00 0.8322D+00 0.1688D+00 0.5339D+00 0.2297D-01 0.1045D+01 0.2361D+00 0.1045D+01 0.2361D+00 0.1122D+01 0.3509D+00 0.8146D+00 0.2062D+00 0.8147D+00 0.2062D+00 0.8432D+00 0.1546D+00 0.3592D+00 0.4570D-01 0.3592D+00 0.4570D-01 0.7604D+00 0.3826D-01 0.7528D+00 0.6731D-01 0.7528D+00 0.6731D-01 0.0000D+00 0.3600D-01 0.0000D+00 0.2704D-01 0.0000D+00 0.2704D-01 RESULTANT CONTINUED FRACTION COEFFICIENTS 0.479091843000D+00 0.600000000000D+01 0.914120235939D+00 0.272556399572D+00 0.103566428624D+01 0.253600669494D+00 0.557354584960D+00 0.109221130287D+00 0.493633149415D+00 0.572026496089D-01 0.643389244701D+00 0.970442487892D-01 0.670391962047D+00 0.186968054287D-01 OMEGA DENSITY OF PHONON MODES INTEGRATED D.O.S. 0.00000000D+00 0.00000000D+00 0.20061257D+01 0.75000000D-01 0.25993608D+01 0.23528144D+01 0.15000000D+00 0.87579080D+01 0.36216510D+01 0.22500000D+00 0.28957244D+01 0.44002676D+01 0.30000000D+00 0.25745059D+01 0.47668401D+01 0.37500000D+00 0.28424823D+01 0.52378133D+01 0.45000000D+00 0.31774934D+01 0.57207138D+01 0.52500000D+00 0.12332777D+01 0.60182049D+01 0.60000000D+00 0.18157666D+01 0.62126757D+01 0.67500000D+00 0.41979454D+01 0.67598706D+01 0.75000000D+00 0.30157468D+01 0.73864784D+01 0.82500000D+00 0.15131930D+01 0.76560981D+01 0.90000000D+00 0.79309662D+01 0.82050211D+01 0.97500000D+00 0.12143623D+02 0.96894785D+01 0.10500000D+01 0.10205478D+01 0.10918245D+02 0.11250000D+01 0.27673040D-01 0.10953992D+02 0.12000000D+01 0.40566780D-02 0.10955601D+02 0.12750000D+01 0.52478330D+00 0.10960721D+02 0.13500000D+01 0.37333981D-03 0.11999986D+02 0.14250000D+01 0.10540232D-04 0.11999999D+02 0.15000000D+01 0.83387870D-06 0.12000000D+02 ./ ENDUP ./ ADD NAME=CFTPEEL CONTINUED FRACTION COEFFICIENTS A(I),B2(I),I=1, 11 ORBITAL PEELING ON A 76 ATOM FCC CLUSTER 2 POSTIONS 3 ORBITALS PER POSITION ADD ATOM POSITION 0.30000E+01 0.20000E+01 0.00000E+00 ORBITAL 1 0.500000E+00 0.100000E+01 0.472490E+00 0.421663E-01 0.195485E-01 0.148046E-01 -0.474055E-02 0.149503E-02 0.273241E+00 0.494780E-02 0.195113E+00 0.687933E-01 -0.102151E-01 0.334883E-02 -0.147545E-01 0.214047E-02 0.342035E+00 0.637294E-02 0.175120E+00 0.626467E-01 0.000000E+00 0.143511E-01 ORBITAL 2 0.500000E+00 0.100000E+01 -0.210645E-02 0.217033E-02 0.189660E-01 0.164294E-02 0.418913E+00 0.116863E-01 0.243933E-01 0.285045E-01 -0.138919E-01 0.225610E-02 -0.208938E-01 0.223641E-02 0.191802E-01 0.223003E-02 0.424895E+00 0.159673E-01 0.386033E-01 0.232389E-01 0.000000E+00 0.202028E-02 ORBITAL 3 0.500000E+00 0.100000E+01 0.150466E-01 0.376133E-02 0.137119E-01 0.141815E-02 0.442441E+00 0.130259E-01 0.238727E-01 0.165956E-01 0.974502E-02 0.165485E-02 0.414811E+00 0.103798E-01 0.488939E-01 0.332176E-01 -0.186314E-01 0.236497E-02 -0.245824E-01 0.223357E-02 0.000000E+00 0.222820E-02 ADD ATOM POSITION 0.35000E+01 0.20000E+01 0.00000E+00 ORBITAL 1 0.500000E+00 0.100000E+01 0.479266E+00 0.418092E-01 0.149541E-01 0.116840E-01 -0.393662E-02 0.805302E-03 0.323287E+00 0.399280E-02 0.164708E+00 0.588453E-01 -0.130959E-01 0.296249E-02 0.249305E+00 0.391418E-02 0.218901E+00 0.716211E-01 0.293668E-01 0.434358E-02 0.000000E+00 0.267288E-01 ORBITAL 2 0.500000E+00 0.100000E+01 0.121590E-02 0.145422E-02 0.149097E+00 0.220581E-02 0.327785E+00 0.569488E-01 -0.534916E-02 0.639670E-02 -0.176335E-01 0.227237E-02 -0.143229E-01 0.192824E-02 0.305057E-01 0.269906E-02 0.429294E+00 0.240395E-01 0.144589E-01 0.153132E-01 0.000000E+00 0.199830E-02 ORBITAL 3 0.500000E+00 0.100000E+01 -0.607080E-02 0.321656E-02 0.165253E+00 0.198901E-02 0.330382E+00 0.577389E-01 -0.199082E-01 0.466610E-02 0.185080E-01 0.211702E-02 0.413457E+00 0.199067E-01 0.211744E-01 0.232999E-01 -0.108829E-01 0.218336E-02 -0.124163E-01 0.188544E-02 0.000000E+00 0.175913E-02 ./ ADD NAME=EXORPEEL RUNFILE FROM %H=ENDJOB PRINT=SAVE FTVSCLR PROGRAM=%H=ENDPROG DD=FT01F001=.CFTPEEL/W/FB +++ LIBRARY=SSTLIB.TCMLIB.COM+SSTLIB.RECLIB+PUBLIC.CAMLIB INTEGER*2 MM,NN,IZERO,IZP,MEM,IW DIMENSION VEC(3,60),CRD(3,100),IZP(100),MEM(3,15),IW(2,60) DIMENSION PSI(5,100),PMN(5,100),A(34),B2(34) COMMON /BLKREC/NN(100,15),MM(100,15),NAT,NP,EE(5,5,60) 1,IZERO(100) COMMON DUM(3) EXTERNAL SLKODE,HOP,FCCBND,EQUIV,PADDOV,PADBND,PADD C C THIS PROGRAM COMPUTES THE RECURSION COEFFICIENTS FOR THE LOCAL D.O.S. C FOR THE 3 P-ORBITALS ON A SURFACE ATOM OF A SMALL F.C.C. CLUSTER OF C ATOMS INTERACTING VIA THEIR D-ELECTRONS. C NNDIM=100 INTDIM=15 NED=60 LL=11 NORP=3 NP=5 NTYPE=1 NX=5 NY=5 NZ=6 NPOSN=2 C C NNDIM IS THE FIRST DIMENSION OF ARRAYS NN & MM AND LAST OF CRD C INTDIM IS THE SECOND DIMENSION OF ARRAYS NN & MM C NED IS THE LAST DIMENSION OF IW,VEC AND EE C LL IS THE REQUIRED LENGTH OF THE CONTINUED FRACTION C NORP THE NUMBER OF ORBITALS TO BE PEELED C NP IS THE NUMBER OF ORBITALS PER ATOM SITE C NTYPE IS THE NUMBER OF 'TYPES' OF ATOM SITE INITIALLY C FCCLAT GENERATES A F.C.C. LATTICE OF SIZE NX X NY X NZ C NPOSN IS THE NUMBER OF ADD ATOM POSTIONS REQUIRED C NAT=NNDIM CALL FCCLAT(CRD,3,IZP,NAT,NX,NY,NZ,NTYPE) IF(NAT.EQ.0)STOP WRITE(6,4) 4 FORMAT(' EXAMPLE OF RECURSION WITH ORBITAL PEELING', 1 /20H LATTICE COORDINATES) WRITE(6,5)(I,(CRD(J,I),J=1,3),I=1,NAT) 5 FORMAT(4(I5,3F4.1)) CALL CLCK4X C C SETUP ASSEMBLES THE HAMILTONIAN MATRIX USING THE SUPPLIED ROUTINES C EQUIV,SLKODE, AND FCCBFE, AND THE LIBRARY ROUTINES NNCAL,MMCAL . C THE INFORMATION IS STORED IN THE COMMON BLOCK /BLKREC/ FOR USE BY HOP C C NN IS THE NEIGHBOUR MAP GENERATED BY NNCAL C MM IS THE INTERACTION MAP GENERATED BY MMCAL C NAT IS THE NO. OF ATOMS IN THE CLUSTER C NP IS THE NO. OF 'ORBITALS' ON EACH SITE C EE IS THE LIST OF SUBMATRICES C IZERO IS SUCH THAT IZERO(I)=0 IMPLIES PSI(.,I)=0.0 C C FOR A CLUSTER OF MORE THAN 100 ATOMS ONLY THE OBVIOUS DIMENSIONS C ABOVE NEED TO BE CHANGED , LIKEWISE FOR MORE THAN 15 C CONNECTED ATOMS ; THE ONLY OTHER CHANGES ARE THE CORRESPONDING C DIMENSION INFORMATION NNDIM AND INTDIM C C SETUP IS CALLED TO GENERATE THE D-ELECTRON INTERACTIONS OF THE C FCC CLUSTER C NM=INTDIM CALL SETUP(CRD,3,NAT,EQUIV,NTYPE,IZP,MM,NN,NNDIM,NM,SLKODE 1,FCCBND,FCCBND,EE,NP,NED,NE,VEC,IW) CALL CLCK3F(IT) WRITE(6,6)IT 6 FORMAT(' TIME IN SETUP',I6) CALL OUT(6,1,1,IZP,IW,VEC,NED,NE,NAT,MM,NN,NNDIM,INTDIM,EE,NP) C C ADD ONE ATOM TO THE CLUSTER AND COMPUTE THE P-ELECTRON INTERACTIONS C USING THE ROUTINE ADDAT C NAT=NAT+1 CRD(3,NAT)=0.0 CRD(2,NAT)=2.0 CRD(1,NAT)=2.0 IZP(NAT)=-2 CALL CLCK4X NE1=NE NM=INTDIM CALL ADDAT(CRD,3,NAT,EQUIV,IZP,MM,NN,NNDIM,NM,PADBND, 1NE,EE,NP,VEC,IW,NED,PADDOV,PADD) CALL CLCK3F(IT) WRITE(6,7)IT 7 FORMAT(' TIME IN ADDAT',I6) CALL OUT(6,NAT,NE1,IZP,IW,VEC,NED,NE,NAT,MM,NN,NNDIM,INTDIM,EE,NP) IZP(NAT)=-IZP(NAT) WRITE(1,12)LL WRITE(1,10)NAT,NPOSN,NORP 10 FORMAT(' ORBITAL PEELING ON A ',I6,' ATOM FCC CLUSTER',/I4 1,' POSTIONS',I3,' ORBITALS PER POSITION') C C ADD A SECOND ATOM TO THE SURFACE IN VARIOUS PLACES C DO 17 NADD=1,NPOSN NAT=NAT+1 CRD(3,NAT)=0.0 CRD(2,NAT)=2.0 CRD(1,NAT)=FLOAT(NADD)*0.5+2.5 IZP(NAT)=-2 NM=INTDIM CALL ADDAT(CRD,3,NAT,EQUIV,IZP,MM,NN,NNDIM,NM,PADBND, 1NE,EE,NP,VEC,IW,NED,PADDOV,PADD) CALL OUT(6,NAT,0,IZP,IW,VEC,NED,NE,NAT,MM,NN,NNDIM,INTDIM,EE,NP) C C CONSTRUCT INITIAL WAVE FUNCTION NONZERO ON LAST ATOM AND C COMPUTE CONTINUED FRACTION C NSTRT=NAT WRITE(6,44) NSTRT 44 FORMAT(I6,' ATOM TAKEN AS STARTING POINT'//) C C COMPUTES LOCAL D.O.S. RECURSION CFTS. FOR THE THREE ORBITALS ON THE C STARTING ATOM. IZERO IS AN ARRAY SET UP TO AVIOD UNNECCESSARY MATRIX C MULTIPLICATIONS: IZERO(I)=0 IMPLIES PSI(.,I)=0.0 AND C THE ROUTINE HOP UPDATES THIS WITH EACH APPLICATION OF THE MATRIX. C WRITE(1,8)(CRD(K,NSTRT),K=1,3) 8 FORMAT(' ADD ATOM POSITION',3E13.5) 9 FORMAT(' ORBITAL',I3) DO 15 III=1,NORP CALL CLCK4X DO 11 I=1,NAT IZERO(I)=0 DO 11 J=1,NP PSI(J,I)=0.0 11 PMN(J,I)=0.0 B2(1)=1.0 A(LL)=0.0 PSI(III,NSTRT)=1.0 IZERO(NSTRT)=1 CALL RECAL(HOP,PSI,PMN,NAT*NP,A,B2,LL) WRITE(1,9)III 12 FORMAT(48H CONTINUED FRACTION COEFFICIENTS A(I),B2(I),I=1,,I3) WRITE(1,13)(A(I),B2(I),I=1,LL) 13 FORMAT(2E14.6) CALL CLCK3F(IT) WRITE(6,14)IT 14 FORMAT(35H TIME TO COMPUTE TRIDIAGONALISATION,I6) ICODE=III IF(III.EQ.NORP)ICODE=NP C C 'PEEL' THE ORBITAL FOR WHICH WE HAVE COMPUTED THE CONTINUED FRACTION C CALL ORPEEL(NSTRT,III,ICODE,MM,NN,NNDIM,INTDIM,EE,NP,NE,NED,MEM) 15 CONTINUE C C DETACH ATOM READY FOR A RUN WITH A NEW POSITION C MEM(1,1)=1 DO 16 I=1,NAT 16 IZERO(I)=IZP(I) IZERO(NAT)=3 CALL PEEL(CRD,3,NAT,NN,NNDIM,INTDIM,MEM,1,IZP,IZERO,2,PSI) 17 CONTINUE STOP END SUBROUTINE HOP(PSI,PMN,AN) INTEGER*2 MM,NN,IZERO,IDUM DOUBLE PRECISION SUM,DUM DIMENSION PSI(5,100),PMN(5,100),DUM(5) COMMON /BLKREC/NN(100,15),MM(100,15),NAT,NP,EE(5,5,60),IZERO(100) COMMON IDUM(100) C C HOP SATISFIES THE SPEC. IN RECAL. UNLABELLED COMMON IS USED FOR C AN INTEGER*2 WORK ARRAY.ALL INFORMATION FOR THE MATRIX DEFINITION C IS HELD IN THE COMMON BLOCK /BLKREC/ WHICH WAS ASSIGNED BY A CALL TO C SETUP C C IZERO(I)=0 INDICATES PSI(.,I)=0.0 AND HENCE THAT CERTAIN LINES C MAY BE SKIPPED TO SAVE TIME. C SUM=0.0 DO 4 I=1,NAT C C THE MATRIX MULTIPLICATION PROCEEDS A LINE AT A TIME BY COMPUTING C AT EACH SITE I THE EFFECT OF THE HAMILTONIAN OPERATING AT EACH C OF ITS NEIGHBOURING SITES, WHICH ARE DEFINED BY NN(I,.) C THE PRODUCT VECTOR IS ACCUMULATED IN DUM AND THE INNER PRODUCT C IN SUM. C IDUM(I)=IZERO(I) DO 11 L=1,NP 11 DUM(L)=0.0 JSZ=NN(I,1) IE=MM(I,1) IF(IZERO(I).EQ.0)GOTO 10 C C MULTIPLY PSI(I) BY THE 'SELF ENERGY' MATRIX (OR DIAGONAL SUBMATRIX) C DO 1 M=1,NP DO 1 L=1,NP 1 DUM(L)=DUM(L)+EE(L,M,IE)*PSI(M,I) 10 IF(JSZ.LT.2)GOTO 3 DO 12 J=2,JSZ JJ=NN(I,J) IF(IZERO(JJ).EQ.0)GOTO 12 IE=MM(I,J) C C MULTIPLY PSI(.,JJ) BY THE APPROPRIATE SUBMATRIX TO CALCULATE THE C EFFECT AT SITE I OF THE HAMILTONIAN OPERATING AT SITE JJ.NOTE C THE PSI(.,I) IS NOW NON-ZERO AND THIS IS RECORDED IN IDUM FOR C LATER COPYING TO IZERO C DO 2 M=1,NP DO 2 L=1,NP 2 DUM(L)=DUM(L)+EE(L,M,IE)*PSI(M,JJ) IDUM(I)=1 12 CONTINUE C C ACCUMULATE THE REQUIRED INNER PRODUCT AND OVERWRITE PMN C 3 DO 4 L=1,NP SUM=SUM+DUM(L)*PSI(L,I) 4 PMN(L,I)=DUM(L)+PMN(L,I) AN=SUM DO 14 I=1,NAT 14 IZERO(I)=IDUM(I) RETURN END INTEGER FUNCTION PADBND(I,J,R2,DD) C C FUNCTION DEFINING 'NEIGHBOUR' FOR THE ADD ATOM : C TAKES THE VALUE 1 IF I & J ARE NEIGHBOURS , 0 OTHERWISE C PADBND=0 IF(R2.GT.4.0)RETURN PADBND=1 RETURN END SUBROUTINE PADD(DUM,I,J,EM,OVP) DIMENSION DUM(3),X(6),X2(6),EM(5,5),E(5,5) EXTERNAL OVP C C COMPUTES THE OVERLAP MATRICES FOR THE ADD ATOM : P-D AND PP C ELECTRONS . SATISFIES THE SPEC OF HCAL FOR SETUP AND ADDAT. C R2=0.0 DO 1 L=1,3 X(L)=DUM(L) X2(L)=X(L)*X(L) 1 R2=R2+X2(L) DO 6 L=1,5 DO 6 M=1,5 6 EM(L,M)=0.0 IF(R2.LT.1.0E-4)GOTO 5 R2I=1.0/R2 RI=SQRT(R2I) DO 2 L=1,3 X(L)=X(L)*RI X(L+3)=X(L) X2(L)=X2(L)*R2I 2 X2(L+3)=X2(L) IF(I.EQ.J)GOTO 4 IF(J.EQ.1)GOTO 3 CALL SKPD(X,X2,I,J,R2,OVP,E,EM,5) RETURN 3 CALL SKPD(X,X2,I,J,R2,OVP,EM,E,5) RETURN 4 CALL SKPP(X,X2,I,J,R2,OVP,EM,5) RETURN 5 CALL SELFP(I,J,R2,OVP,EM,5) RETURN END INTEGER FUNCTION PADDOV(I,J,R2,PAR) DIMENSION PAR(13) C C RETURNS THE VALUES OF THE OVERLAP PARAMETERS FOR THE ADD ATOMS. C C PAR(4)= P D SIGMA C PAR(5)= P D PI C PAR(6)= P P SIGMA C PAR(7)= P P PI C PAR(12) = SELF ENERGY OF P ELECTRON C PADDOV=1 PAR(4)=0.04 PAR(5)=-0.003 PAR(6)=0.2 PAR(7)=-0.002 PAR(12)=0.5 RETURN END ENDPROG ENDJOB ./ ADD NAME=OPORPEEL EXAMPLE OF RECURSION WITH ORBITAL PEELING LATTICE COORDINATES 1 2.0 1.0 1.0 2 4.0 1.0 1.0 3 1.0 2.0 1.0 4 3.0 2.0 1.0 5 5.0 2.0 1.0 6 2.0 3.0 1.0 7 4.0 3.0 1.0 8 1.0 4.0 1.0 9 3.0 4.0 1.0 10 5.0 4.0 1.0 11 2.0 5.0 1.0 12 4.0 5.0 1.0 13 1.0 1.0 2.0 14 3.0 1.0 2.0 15 5.0 1.0 2.0 16 2.0 2.0 2.0 17 4.0 2.0 2.0 18 1.0 3.0 2.0 19 3.0 3.0 2.0 20 5.0 3.0 2.0 21 2.0 4.0 2.0 22 4.0 4.0 2.0 23 1.0 5.0 2.0 24 3.0 5.0 2.0 25 5.0 5.0 2.0 26 2.0 1.0 3.0 27 4.0 1.0 3.0 28 1.0 2.0 3.0 29 3.0 2.0 3.0 30 5.0 2.0 3.0 31 2.0 3.0 3.0 32 4.0 3.0 3.0 33 1.0 4.0 3.0 34 3.0 4.0 3.0 35 5.0 4.0 3.0 36 2.0 5.0 3.0 37 4.0 5.0 3.0 38 1.0 1.0 4.0 39 3.0 1.0 4.0 40 5.0 1.0 4.0 41 2.0 2.0 4.0 42 4.0 2.0 4.0 43 1.0 3.0 4.0 44 3.0 3.0 4.0 45 5.0 3.0 4.0 46 2.0 4.0 4.0 47 4.0 4.0 4.0 48 1.0 5.0 4.0 49 3.0 5.0 4.0 50 5.0 5.0 4.0 51 2.0 1.0 5.0 52 4.0 1.0 5.0 53 1.0 2.0 5.0 54 3.0 2.0 5.0 55 5.0 2.0 5.0 56 2.0 3.0 5.0 57 4.0 3.0 5.0 58 1.0 4.0 5.0 59 3.0 4.0 5.0 60 5.0 4.0 5.0 61 2.0 5.0 5.0 62 4.0 5.0 5.0 63 1.0 1.0 6.0 64 3.0 1.0 6.0 65 5.0 1.0 6.0 66 2.0 2.0 6.0 67 4.0 2.0 6.0 68 1.0 3.0 6.0 69 3.0 3.0 6.0 70 5.0 3.0 6.0 71 2.0 4.0 6.0 72 4.0 4.0 6.0 73 1.0 5.0 6.0 74 3.0 5.0 6.0 75 5.0 5.0 6.0 TIME IN SETUP 10 NEAREST NEIGHBOUR MAP ATOM TYPE CONNECTIVITY NEIGHBOURS 1 1 6 3 4 13 14 16 2 1 6 4 5 14 15 17 3 1 6 1 6 13 16 18 4 1 9 1 2 6 7 14 16 17 19 5 1 6 2 7 15 17 20 6 1 9 3 4 8 9 16 18 19 21 7 1 9 4 5 9 10 17 19 20 22 8 1 6 6 11 18 21 23 9 1 9 6 7 11 12 19 21 22 24 10 1 6 7 12 20 22 25 11 1 6 8 9 21 23 24 12 1 6 9 10 22 24 25 13 1 6 1 3 16 26 28 14 1 9 1 2 4 16 17 26 27 29 15 1 6 2 5 17 27 30 16 1 13 1 3 4 6 13 14 18 19 26 28 29 31 17 1 13 2 4 5 7 14 15 19 20 27 29 30 32 18 1 9 3 6 8 16 21 28 31 33 19 1 13 4 6 7 9 16 17 21 22 29 31 32 34 20 1 9 5 7 10 17 22 30 32 35 21 1 13 6 8 9 11 18 19 23 24 31 33 34 36 22 1 13 7 9 10 12 19 20 24 25 32 34 35 37 23 1 6 8 11 21 33 36 24 1 9 9 11 12 21 22 34 36 37 25 1 6 10 12 22 35 37 26 1 9 13 14 16 28 29 38 39 41 27 1 9 14 15 17 29 30 39 40 42 28 1 9 13 16 18 26 31 38 41 43 29 1 13 14 16 17 19 26 27 31 32 39 41 42 44 30 1 9 15 17 20 27 32 40 42 45 31 1 13 16 18 19 21 28 29 33 34 41 43 44 46 32 1 13 17 19 20 22 29 30 34 35 42 44 45 47 33 1 9 18 21 23 31 36 43 46 48 34 1 13 19 21 22 24 31 32 36 37 44 46 47 49 35 1 9 20 22 25 32 37 45 47 50 36 1 9 21 23 24 33 34 46 48 49 37 1 9 22 24 25 34 35 47 49 50 38 1 6 26 28 41 51 53 39 1 9 26 27 29 41 42 51 52 54 40 1 6 27 30 42 52 55 41 1 13 26 28 29 31 38 39 43 44 51 53 54 56 42 1 13 27 29 30 32 39 40 44 45 52 54 55 57 43 1 9 28 31 33 41 46 53 56 58 44 1 13 29 31 32 34 41 42 46 47 54 56 57 59 45 1 9 30 32 35 42 47 55 57 60 46 1 13 31 33 34 36 43 44 48 49 56 58 59 61 47 1 13 32 34 35 37 44 45 49 50 57 59 60 62 48 1 6 33 36 46 58 61 49 1 9 34 36 37 46 47 59 61 62 50 1 6 35 37 47 60 62 51 1 9 38 39 41 53 54 63 64 66 52 1 9 39 40 42 54 55 64 65 67 53 1 9 38 41 43 51 56 63 66 68 54 1 13 39 41 42 44 51 52 56 57 64 66 67 69 55 1 9 40 42 45 52 57 65 67 70 56 1 13 41 43 44 46 53 54 58 59 66 68 69 71 57 1 13 42 44 45 47 54 55 59 60 67 69 70 72 58 1 9 43 46 48 56 61 68 71 73 59 1 13 44 46 47 49 56 57 61 62 69 71 72 74 60 1 9 45 47 50 57 62 70 72 75 61 1 9 46 48 49 58 59 71 73 74 62 1 9 47 49 50 59 60 72 74 75 63 1 4 51 53 66 64 1 6 51 52 54 66 67 65 1 4 52 55 67 66 1 9 51 53 54 56 63 64 68 69 67 1 9 52 54 55 57 64 65 69 70 68 1 6 53 56 58 66 71 69 1 9 54 56 57 59 66 67 71 72 70 1 6 55 57 60 67 72 71 1 9 56 58 59 61 68 69 73 74 72 1 9 57 59 60 62 69 70 74 75 73 1 4 58 61 71 74 1 6 59 61 62 71 72 75 1 4 60 62 72 13 NON-EQUIVALENT VECTORS FOUND ATOM INDEX OF INTERACTION VECTORS 1 13 1 2 3 4 5 2 13 1 2 3 4 5 3 13 6 2 7 4 5 4 13 8 6 1 2 7 3 4 5 5 13 8 1 7 3 5 6 13 8 6 1 2 7 3 4 5 7 13 8 6 1 2 7 3 4 5 8 13 6 2 7 4 5 9 13 8 6 1 2 7 3 4 5 10 13 8 1 7 3 5 11 13 8 6 7 3 4 12 13 8 6 7 3 4 13 13 9 10 2 4 5 14 13 11 9 10 1 2 3 4 5 15 13 11 10 1 3 5 16 13 12 11 9 10 8 6 1 2 7 3 4 5 17 13 12 11 9 10 8 6 1 2 7 3 4 5 18 13 12 9 10 6 2 7 4 5 19 13 12 11 9 10 8 6 1 2 7 3 4 5 20 13 12 11 10 8 1 7 3 5 21 13 12 11 9 10 8 6 1 2 7 3 4 5 22 13 12 11 9 10 8 6 1 2 7 3 4 5 23 13 12 9 6 7 4 24 13 12 11 9 8 6 7 3 4 25 13 12 11 8 7 3 26 13 11 9 10 1 2 3 4 5 27 13 11 9 10 1 2 3 4 5 28 13 12 9 10 6 2 7 4 5 29 13 12 11 9 10 8 6 1 2 7 3 4 5 30 13 12 11 10 8 1 7 3 5 31 13 12 11 9 10 8 6 1 2 7 3 4 5 32 13 12 11 9 10 8 6 1 2 7 3 4 5 33 13 12 9 10 6 2 7 4 5 34 13 12 11 9 10 8 6 1 2 7 3 4 5 35 13 12 11 10 8 1 7 3 5 36 13 12 11 9 8 6 7 3 4 37 13 12 11 9 8 6 7 3 4 38 13 9 10 2 4 5 39 13 11 9 10 1 2 3 4 5 40 13 11 10 1 3 5 41 13 12 11 9 10 8 6 1 2 7 3 4 5 42 13 12 11 9 10 8 6 1 2 7 3 4 5 43 13 12 9 10 6 2 7 4 5 44 13 12 11 9 10 8 6 1 2 7 3 4 5 45 13 12 11 10 8 1 7 3 5 46 13 12 11 9 10 8 6 1 2 7 3 4 5 47 13 12 11 9 10 8 6 1 2 7 3 4 5 48 13 12 9 6 7 4 49 13 12 11 9 8 6 7 3 4 50 13 12 11 8 7 3 51 13 11 9 10 1 2 3 4 5 52 13 11 9 10 1 2 3 4 5 53 13 12 9 10 6 2 7 4 5 54 13 12 11 9 10 8 6 1 2 7 3 4 5 55 13 12 11 10 8 1 7 3 5 56 13 12 11 9 10 8 6 1 2 7 3 4 5 57 13 12 11 9 10 8 6 1 2 7 3 4 5 58 13 12 9 10 6 2 7 4 5 59 13 12 11 9 10 8 6 1 2 7 3 4 5 60 13 12 11 10 8 1 7 3 5 61 13 12 11 9 8 6 7 3 4 62 13 12 11 9 8 6 7 3 4 63 13 9 10 2 64 13 11 9 10 1 2 65 13 11 10 1 66 13 12 11 9 10 8 6 1 2 67 13 12 11 9 10 8 6 1 2 68 13 12 9 10 6 2 69 13 12 11 9 10 8 6 1 2 70 13 12 11 10 8 1 71 13 12 11 9 10 8 6 1 2 72 13 12 11 9 10 8 6 1 2 73 13 12 9 6 74 13 12 11 9 8 6 75 13 12 11 8 INDEX LATTICE TYPES VECTOR MATRIX 1 1 1 -1.00 1.00 0.00 -0.021 0.000 0.000 0.000 -0.011 0.000 0.005 -0.007 0.000 0.000 0.000 -0.007 0.005 0.000 0.000 0.000 0.000 0.000 0.013 0.000 -0.011 0.000 0.000 0.000 -0.008 2 1 1 1.00 1.00 0.00 -0.021 0.000 0.000 0.000 0.011 0.000 0.005 0.007 0.000 0.000 0.000 0.007 0.005 0.000 0.000 0.000 0.000 0.000 0.013 0.000 0.011 0.000 0.000 0.000 -0.008 3 1 1 -1.00 0.00 1.00 0.005 -0.007 0.000 0.000 0.000 -0.007 0.005 0.000 0.000 0.000 0.000 0.000 -0.021 0.010 0.006 0.000 0.000 0.010 -0.003 -0.009 0.000 0.000 0.006 -0.009 0.007 4 1 1 1.00 0.00 1.00 0.005 0.007 0.000 0.000 0.000 0.007 0.005 0.000 0.000 0.000 0.000 0.000 -0.021 -0.010 -0.006 0.000 0.000 -0.010 -0.003 -0.009 0.000 0.000 -0.006 -0.009 0.007 5 1 1 0.00 1.00 1.00 0.005 0.000 0.007 0.000 0.000 0.000 -0.021 0.000 0.010 -0.006 0.007 0.000 0.005 0.000 0.000 0.000 0.010 0.000 -0.003 0.009 0.000 -0.006 0.000 0.009 0.007 6 1 1 1.00 -1.00 0.00 -0.021 0.000 0.000 0.000 -0.011 0.000 0.005 -0.007 0.000 0.000 0.000 -0.007 0.005 0.000 0.000 0.000 0.000 0.000 0.013 0.000 -0.011 0.000 0.000 0.000 -0.008 7 1 1 0.00 -1.00 1.00 0.005 0.000 -0.007 0.000 0.000 0.000 -0.021 0.000 -0.010 0.006 -0.007 0.000 0.005 0.000 0.000 0.000 -0.010 0.000 -0.003 0.009 0.000 0.006 0.000 0.009 0.007 8 1 1 -1.00 -1.00 0.00 -0.021 0.000 0.000 0.000 0.011 0.000 0.005 0.007 0.000 0.000 0.000 0.007 0.005 0.000 0.000 0.000 0.000 0.000 0.013 0.000 0.011 0.000 0.000 0.000 -0.008 9 1 1 1.00 0.00 -1.00 0.005 -0.007 0.000 0.000 0.000 -0.007 0.005 0.000 0.000 0.000 0.000 0.000 -0.021 0.010 0.006 0.000 0.000 0.010 -0.003 -0.009 0.000 0.000 0.006 -0.009 0.007 10 1 1 0.00 1.00 -1.00 0.005 0.000 -0.007 0.000 0.000 0.000 -0.021 0.000 -0.010 0.006 -0.007 0.000 0.005 0.000 0.000 0.000 -0.010 0.000 -0.003 0.009 0.000 0.006 0.000 0.009 0.007 11 1 1 -1.00 0.00 -1.00 0.005 0.007 0.000 0.000 0.000 0.007 0.005 0.000 0.000 0.000 0.000 0.000 -0.021 -0.010 -0.006 0.000 0.000 -0.010 -0.003 -0.009 0.000 0.000 -0.006 -0.009 0.007 12 1 1 0.00 -1.00 -1.00 0.005 0.000 0.007 0.000 0.000 0.000 -0.021 0.000 0.010 -0.006 0.007 0.000 0.005 0.000 0.000 0.000 0.010 0.000 -0.003 0.009 0.000 -0.006 0.000 0.009 0.007 13 1 1 0.00 0.00 0.00 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 TIME IN ADDAT 0 NEAREST NEIGHBOUR MAP ATOM TYPE CONNECTIVITY NEIGHBOURS 76 -2 6 1 3 4 6 16 24 NON-EQUIVALENT VECTORS FOUND ATOM INDEX OF INTERACTION VECTORS 76 24 19 20 21 22 23 INDEX LATTICE TYPES VECTOR MATRIX 13 1 1 0.00 0.00 0.00 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 14 1 2 0.00 1.00 -1.00 0.002 0.000 0.000 0.000 0.000 0.000 0.024 -0.024 0.000 0.000 -0.002 0.000 0.000 0.000 0.000 0.000 0.011 -0.013 0.000 0.000 0.000 -0.009 0.005 0.000 0.000 15 1 2 1.00 0.00 -1.00 0.000 0.002 0.000 0.000 0.000 0.000 -0.002 0.000 0.000 0.000 0.024 0.000 -0.024 0.000 0.000 -0.011 0.000 0.013 0.000 0.000 -0.009 0.000 0.005 0.000 0.000 16 1 2 -1.00 0.00 -1.00 0.000 -0.002 0.000 0.000 0.000 0.000 -0.002 0.000 0.000 0.000 0.024 0.000 0.024 0.000 0.000 0.011 0.000 0.013 0.000 0.000 0.009 0.000 0.005 0.000 0.000 17 1 2 0.00 -1.00 -1.00 -0.002 0.000 0.000 0.000 0.000 0.000 0.024 0.024 0.000 0.000 -0.002 0.000 0.000 0.000 0.000 0.000 -0.011 -0.013 0.000 0.000 0.000 0.009 0.005 0.000 0.000 18 1 2 0.00 0.00 -2.00 0.000 0.000 0.000 0.000 0.000 0.000 -0.003 0.000 0.000 0.000 -0.003 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.040 0.000 0.000 19 2 1 0.00 -1.00 1.00 0.002 0.000 -0.002 0.000 0.000 0.000 0.024 0.000 0.011 -0.009 0.000 -0.024 0.000 -0.013 0.005 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 ***** #PRC SUBSYSTEM STREAM RESUMPTION ***** 20 2 1 -1.00 0.00 1.00 0.000 0.000 0.024 -0.011 -0.009 0.002 -0.002 0.000 0.000 0.000 0.000 0.000 -0.024 0.013 0.005 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 21 2 1 1.00 0.00 1.00 0.000 0.000 0.024 0.011 0.009 -0.002 -0.002 0.000 0.000 0.000 0.000 0.000 0.024 0.013 0.005 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 22 2 1 0.00 1.00 1.00 -0.002 0.000 -0.002 0.000 0.000 0.000 0.024 0.000 -0.011 0.009 0.000 0.024 0.000 -0.013 0.005 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 23 2 1 0.00 0.00 2.00 0.000 0.000 -0.003 0.000 0.000 0.000 -0.003 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.040 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 24 2 2 0.00 0.00 0.00 0.500 0.000 0.000 0.000 0.000 0.000 0.500 0.000 0.000 0.000 0.000 0.000 0.500 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 NEAREST NEIGHBOUR MAP ATOM TYPE CONNECTIVITY NEIGHBOURS 77 -2 7 1 2 4 6 7 76 36 NON-EQUIVALENT VECTORS FOUND ATOM INDEX OF INTERACTION VECTORS 77 24 31 32 33 34 35 36 77 ATOM TAKEN AS STARTING POINT TIME TO COMPUTE TRIDIAGONALISATION 20 TIME TO COMPUTE TRIDIAGONALISATION 20 TIME TO COMPUTE TRIDIAGONALISATION 19 NEAREST NEIGHBOUR MAP ATOM TYPE CONNECTIVITY NEIGHBOURS 77 -2 6 2 4 5 7 76 46 NON-EQUIVALENT VECTORS FOUND ATOM INDEX OF INTERACTION VECTORS 77 24 42 43 44 45 46 77 ATOM TAKEN AS STARTING POINT TIME TO COMPUTE TRIDIAGONALISATION 19 TIME TO COMPUTE TRIDIAGONALISATION 19 TIME TO COMPUTE TRIDIAGONALISATION 19 ./ ENDUP ./ ADD NAME=CFTPEEL CONTINUED FRACTION COEFFICIENTS A(I),B2(I),I=1, 11 ORBITAL PEELING ON A 76 ATOM FCC CLUSTER 2 POSTIONS 3 ORBITALS PER POSITION ADD ATOM POSITION 0.30000E+01 0.20000E+01 0.00000E+00 ORBITAL 1 0.500000E+00 0.100000E+01 0.472490E+00 0.421663E-01 0.195485E-01 0.148046E-01 -0.474055E-02 0.149503E-02 0.273241E+00 0.494780E-02 0.195113E+00 0.687933E-01 -0.102151E-01 0.334883E-02 -0.147545E-01 0.214047E-02 0.342035E+00 0.637294E-02 0.175120E+00 0.626467E-01 0.000000E+00 0.143511E-01 ORBITAL 2 0.500000E+00 0.100000E+01 -0.210645E-02 0.217033E-02 0.189660E-01 0.164294E-02 0.418913E+00 0.116863E-01 0.243933E-01 0.285045E-01 -0.138919E-01 0.225610E-02 -0.208938E-01 0.223641E-02 0.191802E-01 0.223003E-02 0.424895E+00 0.159673E-01 0.386033E-01 0.232389E-01 0.000000E+00 0.202028E-02 ORBITAL 3 0.500000E+00 0.100000E+01 0.150466E-01 0.376133E-02 0.137119E-01 0.141815E-02 0.442441E+00 0.130259E-01 0.238727E-01 0.165956E-01 0.974502E-02 0.165485E-02 0.414811E+00 0.103798E-01 0.488939E-01 0.332176E-01 -0.186314E-01 0.236497E-02 -0.245824E-01 0.223357E-02 0.000000E+00 0.222820E-02 ADD ATOM POSITION 0.35000E+01 0.20000E+01 0.00000E+00 ORBITAL 1 0.500000E+00 0.100000E+01 0.479266E+00 0.418092E-01 0.149541E-01 0.116840E-01 -0.393662E-02 0.805302E-03 0.323287E+00 0.399280E-02 0.164708E+00 0.588453E-01 -0.130959E-01 0.296249E-02 0.249305E+00 0.391418E-02 0.218901E+00 0.716211E-01 0.293668E-01 0.434358E-02 0.000000E+00 0.267288E-01 ORBITAL 2 0.500000E+00 0.100000E+01 0.121590E-02 0.145422E-02 0.149097E+00 0.220581E-02 0.327785E+00 0.569488E-01 -0.534916E-02 0.639670E-02 -0.176335E-01 0.227237E-02 -0.143229E-01 0.192824E-02 0.305057E-01 0.269906E-02 0.429294E+00 0.240395E-01 0.144589E-01 0.153132E-01 0.000000E+00 0.199830E-02 ORBITAL 3 0.500000E+00 0.100000E+01 -0.607080E-02 0.321656E-02 0.165253E+00 0.198901E-02 0.330382E+00 0.577389E-01 -0.199082E-01 0.466610E-02 0.185080E-01 0.211702E-02 0.413457E+00 0.199067E-01 0.211744E-01 0.232999E-01 -0.108829E-01 0.218336E-02 -0.124163E-01 0.188544E-02 0.000000E+00 0.175913E-02 ./ ADD NAME=EXPEELED RUNFILE FROM %H=ENDJOB PRINT=SAVE FILE SSTLIB.CAMREC84:CFTPEEL TO &DAT/W/FB FTVSCLR PROGRAM=%H=ENDPROG DD=FT05F001=&DAT LIBRARY=SSTLIB.TCMLIB.DCOM IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(30,3),B2(30,3),WK(30,4),IWK(30,2),HEAD(6),E(201) 1,EN(201,10),NL(2) NPTS=17 ND=30 ELO=-0.2 EHI=0.6 EPS=5.0E-12 LL=7 C C NPTS IS THE NUMBER OF POINTS TO BE USED IN PLOTTING C ELO-EHI IS THE RANGE OF ENERGIES OF INTEREST, C EPS IS THE ACCURACY REQUIRED C LL IS THE LENGTH OF THE C.F. C C INPUT THE RESULTS OF THE ORPEEL EXAMPLE PROGRAM C C C LLIN LENGTH OF C.F.S ON INPUT C NPOS NUMBER OF SETS OF C.F.S C NORP NUMBER OF C.F.S IN EACH SET C READ(5,1)LLIN 1 FORMAT(48X,I3) LL=MIN0(LL,LLIN) WRITE(6,2)LL 2 FORMAT(' EXAMPLE OF ENERGY CALCULATION IN ORBITAL PEELING',/I4, 1' LEVELS USED') READ(5,3)HEAD,NPOS,NORP 3 FORMAT(22X,6A4,/I4,9X,I3) WRITE(6,32)HEAD,NPOS,NORP 32 FORMAT(1X,6A4,I4,' SETS OF',I4,' ORBITALS EACH') C C TABULATE VALUES OF FERMI ENERGY REQUIRED C DE=(EHI-ELO)/FLOAT(NPTS-1) DO 34 I=1,NPTS 34 E(I)=ELO+FLOAT(I-1)*DE C C COMPUTE TOTAL ENERGY DIFFERENCE FOR EACH SET OF COEFFICIENTS C DO 10 IPOS=1,NPOS READ(5,35)(HEAD(I),I=1,4),(WK(1,K),K=1,3) 35 FORMAT(2X,4A4,3E13.5) WRITE(6,36)(HEAD(I),I=1,4),(WK(1,K),K=1,3) 36 FORMAT(///2X,4A4,3F6.2) DO 6 I=1,NORP READ(5,1) READ(5,4)(A(J,I),B2(J,I),J=1,LLIN) 4 FORMAT(2E14.6) 5 FORMAT(' INPUT CONTINUED FRACTION CFTS.',50(/1X,6E11.4)) 6 CONTINUE WRITE(6,5)((A(I,J),B2(I,J),J=1,NORP),I=1,LL) C C COMPUTE THE TOTAL ENERGY DIFFERENCE AT EACH VALUE OF C THE FERMI ENERGY AND STORE THE RESULT IN EN C DO 9 I=1,NPTS EN(I,IPOS)=0.0 DO 8 K=1,NORP 8 EN(I,IPOS)=EN(I,IPOS)+EDIFF(E(I),A(1,K),B2(1,K),LL,EPS,WK,IWK) 9 CONTINUE 10 CONTINUE WRITE(6,105)(I,I=1,NPOS) 105 FORMAT(////' TOTAL ENERGY DIFFERENCES TABULATED AGAINST FERMI ENER 1GY',/7X,'EF',9(5X,I3,' POSN')) DO 11 I=1,NPTS WRITE(6,12)E(I),(EN(I,J),J=1,NPOS) 11 CONTINUE 12 FORMAT(10E13.5) STOP END ENDPROG ENDJOB ./ ADD NAME=OPPEELED EXAMPLE OF ENERGY CALCULATION IN ORBITAL PEELING 7 LEVELS USED 76 ATOM FCC CLUSTER 2 SETS OF 3 ORBITALS EACH DD ATOM POSITION 3.00 2.00 0.00 INPUT CONTINUED FRACTION CFTS. 0.5000D+00 0.1000D+01 0.5000D+00 0.1000D+01 0.5000D+00 0.1000D+01 0.4725D+00 0.4217D-01-0.2106D-02 0.2170D-02 0.1505D-01 0.3761D-02 0.1955D-01 0.1480D-01 0.1897D-01 0.1643D-02 0.1371D-01 0.1418D-02 -0.4741D-02 0.1495D-02 0.4189D+00 0.1169D-01 0.4424D+00 0.1303D-01 0.2732D+00 0.4948D-02 0.2439D-01 0.2850D-01 0.2387D-01 0.1660D-01 0.1951D+00 0.6879D-01-0.1389D-01 0.2256D-02 0.9745D-02 0.1655D-02 -0.1022D-01 0.3349D-02-0.2089D-01 0.2236D-02 0.4148D+00 0.1038D-01 DD ATOM POSITION 3.50 2.00 0.00 INPUT CONTINUED FRACTION CFTS. 0.5000D+00 0.1000D+01 0.5000D+00 0.1000D+01 0.5000D+00 0.1000D+01 0.4793D+00 0.4181D-01 0.1216D-02 0.1454D-02-0.6071D-02 0.3217D-02 0.1495D-01 0.1168D-01 0.1491D+00 0.2206D-02 0.1653D+00 0.1989D-02 -0.3937D-02 0.8053D-03 0.3278D+00 0.5695D-01 0.3304D+00 0.5774D-01 0.3233D+00 0.3993D-02-0.5349D-02 0.6397D-02-0.1991D-01 0.4666D-02 0.1647D+00 0.5885D-01-0.1763D-01 0.2272D-02 0.1851D-01 0.2117D-02 -0.1310D-01 0.2962D-02-0.1432D-01 0.1928D-02 0.4135D+00 0.1991D-01 TOTAL ENERGY DIFFERENCES TABULATED AGAINST FERMI ENERGY EF 1 POSN 2 POSN -0.20000D+00 -0.25485D-06 -0.14520D-06 -0.15000D+00 -0.30990D-05 -0.18740D-05 -0.10000D+00 -0.13333D-03 -0.84809D-04 -0.50000D-01 -0.39395D-02 -0.34169D-02 0.14901D-07 -0.11195D-01 -0.10739D-01 0.50000D-01 -0.17915D-01 -0.13708D-01 0.10000D+00 -0.17916D-01 -0.13758D-01 0.15000D+00 -0.17916D-01 -0.13758D-01 0.20000D+00 -0.17917D-01 -0.13758D-01 0.25000D+00 -0.17919D-01 -0.13759D-01 0.30000D+00 -0.17958D-01 -0.13774D-01 0.35000D+00 -0.61194D-01 -0.58537D-01 0.40000D+00 -0.11119D+00 -0.10854D+00 0.45000D+00 -0.16119D+00 -0.15854D+00 0.50000D+00 -0.21124D+00 -0.20863D+00 0.55000D+00 -0.30248D+00 -0.30262D+00 0.60000D+00 -0.40248D+00 -0.40262D+00 ./ ENDUP