************************************************************************ * * * C O N T W V G * * * * RELATIVISTIC CONTINUUM WAVEFUNCTION SOLVER FOR THE GRASP CODE * * * * REFERENCE: K.G. DYALL, ET. AL., COMPUTER PHYSICS COMMUNICATIONS, * * VOL. 55 (1989), PGS. 425-456. * * * * W.F. PERGER, PHYSICS DEPARTMENT, MICHIGAN TECH UNIVERSITY, * * HOUGHTON, MICHIGAN, 49931, USA * * E-MAIL: WFP@MTU.EDU * * * * CHANGES REQUIRED TO GRASP CODE: * * 1) DELETE SUBROUTINES DATAIN AND SCF FROM THE GRASP CODE. * * 2) ADD CONTWVG (THIS CODE). * * 3) IN SUBROUTINE MCDF, AFTER THE LINE: * * : /SEMI/CHK({NC},{OL}),CCR({NC},{OL}) * * ADD THE LINE: * * COMMON/GRIDC/RGRIDC({NPC}),ECONT,HCONT,MHALF,NCONT,NUMWR,IPHASE * * 4) IN SUBROUTINE MCDF, AFTER THE LINE: * * CALL INIT * * ADD THE FOUR LINES: * * IF(ECONT.GT.ZER0) THEN * * CALL MCPIN * * GO TO 4 * * END IF * * 5) IN SUBROUTINE MCDF, CHANGE THE LINE: * * DO 6 NIT = 1,NITIT * * TO: * * 4 DO 6 NIT = 1,NITIT * * 6) IN SUBROUTINE MCDF, AFTER THE FIRST OCCURRENCE OF THE LINE: * * CALL SCF (ACCY,MX,DOALL) * * ADD THE LINE: * * IF (ECONT.GT.ZERO) RETURN * * * * IT SHOULD BE NOTED THAT CONTWVG MUST BE PREPROCESSED, WITH THE SAME * * PREPROCESSOR AS IS USED TO PREPROCESS THE GRASP CODE (BUT WITH THE * * VARIABLE NPC ADDED THE SYMBOL.DATA FILE). * * * ************************************************************************ ************************************************************************ ********************DATAIN.NUM****************************************** ************************************************************************ * * * SECTION 02 D A T A I N * * * * 8 ROUTINES * * * * 01. DATAIN 02. DATR 03. DATANG * * 04. DATSCF 05. CFOUT 06. CSFOUT * * 07. CARDIN 08. FDWTS * * * ************************************************************************ * * PREPROCESSOR SYMBOLS USED : * * SYMBOL USE * * C3 NC + 3 * MC NUMBER OF NONRELATIVISTIC CONFIGURATIONS * MW NUMBER OF NONRELATIVISTIC ORBITALS * NC NUMBER OF CONFIGURATIONS * NRP NUMBER OF RANK/PARITY COMBINATIONS IN MCT/OSCL * NW NUMBER OF ORBITALS * NUC NUMBER OF NUCLEAR PARAMETERS * OL NUMBER OF LEVELS TO OPTIMIZE ON * OP NUMBER OF OPEN NONRELATIVISTIC SUBSHELLS * ORD MAXIMUM ORDER OF FINITE DIFFERENCE FORMULAE * STR LENGTH OF STRING VARIABLES * 2OP 2 * OP * * DB FLAG FOR DOUBLE PRECISION VERSION * CRAY FLAG FOR CRAY 1S VERSION * CYBER FLAG FOR CYBER 205 VERSION (FTN 200) * IBM FLAG FOR IBM VERSION * LOWER FLAG FOR USE OF LOWER CASE LETTERS AS DATA * VAX FLAG FOR VAX UNDER VMS VERSION 4 * ************************************************************************ ************************************************************************ * * SUBROUTINE DATAIN * * * ---------------- SECTION 02 SUBPROGRAM 01 ---------------- * * * * THIS SUBROUTINE CONTROLS THE READING OF THE INPUT DATA. * * * * SUBROUTINES CALLED: CALEN, CARDIN, CFOUT, CSFOUT, DATNR, DATR, * * DATSCF, MATOUT. * * * * LAST REVISION: 14 DEC 1988 * * * ************************************************************************ * {DB: IMPLICIT DOUBLEPRECISION (A-H, O-Z)} CHARACTER*72 IHED,MESSAG CHARACTER*8 ITIME,IDATE,LABL CHARACTER*{STR} CT CHARACTER*2 NH LOGICAL INUSE,LFILE,LOLD,LNEW,LCOEF,LHEAD,LT * PARAMETER (NLAB = 12,SIXTY = 60.0{DB:D|:E}0,ISIXTY = 60) * DIMENSION RT({C3}),IT({C3}),JT({C3}),LT({C3}),KT({C3}),CT({C3}) DIMENSION LABL(NLAB),LABUSD(NLAB) * COMMON/CARD/ICARD : /CONS/ZERO,HALF,TENTH,ONE,TWO,THREE,TEN : /CUP0/TC({NC},{NC}) : /DEF0/TENMAX,EXPMAX,EPSI,PRECIS : /DEBUG/IBUG(6) : /INFORM/IRD,IPD : /MCTA/KA({NRP}),IOPAR({NRP}),NKA : /MCTB/IO(20) : /MCTR/JTC(20) : /NRD4/IPOS({MW}),NPOS({MC}),KPOS({OP},{MC}),NCFN,NWN : /OPT1/ITC(30) : /OPT3/KTC(30) : /OPT4/LTC(30) : /ORB2/IQ({NC},{NW}) : /ORB4/NW,NCF,NP({NW}),NAK({NW}) : /ORB5/NKL({NW}),NKJ({NW}),KAPMAX : /ORB10/NH({NW}) : /OSC1/NLEV,LEVELS({NC}) : /OSC7/CUTOFF : /PRNT/NVEC,IVEC({NC}) COMMON/SYM/ITJPO({NC}),ISPAR({NC}) : /TIME/TIMECK,TIME1,TIME2 : /TITL/IHED,ITIME,IDATE : /WFAC/WFACT COMMON/GRIDC/RGRIDC({NPC}),ECONT,HCONT,MHALF,NCONT,NUMWR,IPHASE * DATA LABL/'MCP ','MCT ','MCBP ','MCDF ','BENA ', : 'OSCL ','END ','LOWFREQ ','LEVELS ','CUTOFF ', : 'PRINT ','CONT '/ * *----------------------------------------------------------------------* * * * TITLE CARD * * * * FORMAT: 1 STRING (1 TO 72 ALPHANUMERIC CHARACTERS) * * * * STRING : [IHED]: USED AS A HEADING IN ALL I/O * * FILES * * * *----------------------------------------------------------------------* * ECONT = -1.0D0 ICARD = 1 * * THIS READ STATEMENT USES THE END PARAMETER TO FIND THE END OF DATA * THIS IS THE ONLY CARD WHICH IS NOT READ BY CARDIN * READ (IRD,300,END = 99) IHED CALL CALEN WRITE (IPD,301) ICARD,IHED * * SET THE OPTIONS TO 0 AND PERFORM OTHER INITIALIZATION * * IF OPTION J IS TO BE SET THEN AT THE APPROPRIATE PLACE * INPUT THE NUMBER J . THIS WILL SET JTC(J) = 1 FOR EXAMPLE . * IF JTC(J) = 0 THE OPTION IS NOT SET . * * ITC CONTAINS OPTIONS FOR THE MCDF MODULE * JTC(J) ,J = 1,9 CONTAINS OPTIONS FOR THE ANGULAR MODULES * JTC(J) ,J = 10,20 CONTROL WHICH MODULES ARE CALLED * KTC CONTAINS OPTIONS FOR BENA MODULE * LTC CONTAINS OPTIONS FOR OSCL MODULE * DO 1 I = 1,30 ITC(I) = 0 KTC(I) = 0 LTC(I) = 0 IF (I .GT. 20) GO TO 1 IO(I) = 0 JTC(I) = 0 IF (I .GT. NLAB) GO TO 1 LABUSD(I) = 0 IF (I .GT. 6) GO TO 1 IBUG(I) = 0 1 CONTINUE NVEC = 0 * *----------------------------------------------------------------------* * * * CONFIGURATION PARAMETERS AND OPTIONS CARD * * * * FORMAT: 2 TO 3 INTEGERS, UP TO 5 LABELS * * * * 1 INTEGER : [NCF]: NUMBER OF CONFIGURATION * * STATE FUNCTIONS (CSF); SPECIFICATION OF A CSF REQUIRES * * ORBITAL, OCCUPATION, AND COUPLING INFORMATION; THE PRE- * * SENT CODE MAY BE GIVEN THIS INFORMATION IN ONE OF TWO * * SCHEMES THAT MUST BE USED CONSISTENTLY IN THE INPUT FOR * * A GIVEN RUN: THE RELATIVISTIC (R) AND NONRELATIVISTIC * * (NR) SCHEMES; DATA PROVIDED USING THE NR SCHEME ARE * * INTERNALLY CONVERTED TO THE R SCHEME WHICH IS APPROPRI- * * ATE FOR A RELATIVISTIC CODE * * 1 INTEGER : [NW]: NUMBER OF ORBITALS; ORBITALS * * FOR THE R SCHEME ARE DEFINED USING AUGMENTED SPECTRO- * * SCOPIC ANGULAR MOMENTUM LABELS: S, P-, P, D-, D, ... * * (SEE THE COMMENTS IN SUBROUTINE DATSCF); IN THE NR CASE * * THE USUAL LABELS ARE APPLICABLE: S, P, D, ... * * 1 INTEGER : [IOP] IF NO LABELS APPEAR ON CARD; * * [IFILE] OTHERWISE: * * [IOP]: IF ABSENT OR "0", THE R SCHEME IS IMPLIED; * * OTHERWISE THE NR CASE IS ASSUMED: DIFFERENT VALUES OF * * [IOP] THEN CONTROL COMPUTATION AND PRINTOUT: IF [IOP] * * IS A NEGATIVE INTEGER, THE SPHERICAL AVERAGE CASE IS * * IMPLIED: THIS IS THE ONLY CASE FOR WHICH COUPLING DATA * * MUST NOT BE GIVEN; IF [IOP] IS "1", "2", OR "3", JJ- * * COUPLED CSF ARE GENERATED FROM THE NR INPUT; IF "2", AN * * LS-COUPLED BASIS IS ALSO GENERATED, AND THE TRANSFORMA- * * TION MATRIX BETWEEN THE JJ- AND LS-COUPLED BASES IS * * IS COMPUTED, AND THIS MATRIX AS WELL AS THE EIGEN- * * VECTORS IN THE LS BASIS ARE PRINTED OUT * * [IFILE]: THE DATASTREAM NUMBER FOR DUMPING TRANSFOR- * * MATION COEFFICIENTS FOR LS-COUPLED AND OTHER BASES * * DEFINED BY SUBSEQUENT INPUT (SEE COMMENTS IN SUBROUTINE * * DATNR) * * LABEL "NOHEAD" : DO NOT PRINT CSF HEADINGS * * LABEL "HEAD" : PRINT CSF HEADINGS * * LABEL "OLDONLY" : PRINT JJ-COUPLED CSF DATA AND * * NOT THOSE IN OTHER BASIS * * LABEL "NEWONLY" : PRINT CSF DATA IN OTHER BASIS * * AND NOT THOSE IN JJ-COUPLING * * LABEL "NOCSF" DO NOT PRINT CSF DATA IN EITHER * * BASIS * * LABEL "CSFONLY" : PRINT ONLY CSF DATA; * * DO NOT PRINT TRANSFORMATION COEFFICIENTS * * LABEL "NOCOEFF" : DO NOT PRINT TRANSFORMATION * * COEFFICIENTS * * LABEL "COEFF" : PRINT TRANSFORMATION COEFFICIENT * * MATRIX * * LABEL "ALL" : PRINT EVERYTHING * * LABEL "NOPRINT" : DO NOT PRINT ANYTHING * * LABEL "NOFILE" : DO NOT COPY DATA TO DATASTREAM * * WITH GIVEN NUMBER * * LABEL "FILE" : COPY DATA TO FILE WITH GIVEN * * DATASTREAM NUMBER * * * * VALUES OF ALL OPTIONAL DATA REMAIN SET THROUGHOUT A RUN * * UNLESS RESET, SO THAT A NUMBER OF SETS OF RESULTS CAN BE * * COPIED TO ONE FILE; WHEN [IFILE] IS CHANGED, THE OLD FILE * * IS REWOUND * * * *----------------------------------------------------------------------* * IFILE = 0 LFILE = .FALSE. LOLD = .TRUE. LNEW = .TRUE. LCOEF = .FALSE. LHEAD = .TRUE. * CALL CARDIN (NN,RT,IT,JT,LT,CT,KT,{C3}) * IF (NN .LT. 2) THEN GO TO 101 ELSEIF (NN .EQ. 2) THEN IF ((KT(1) .NE. 0) .OR. (KT(2) .NE. 0)) GO TO 105 IOP = 0 ELSE IF ((KT(1) .NE. 0) .OR. (KT(2) .NE. 0) .OR. : (KT(3) .NE. 0)) GO TO 105 IF (NN .EQ. 3) THEN IOP = IT(3) LCOEF = IT(3) .EQ. 2 IF (IOP .EQ. 1) LNEW = .FALSE. ELSE IFILE = IT(3) IOP = 2 INQUIRE (UNIT = IFILE,OPENED = INUSE) IF (.NOT. INUSE) OPEN (UNIT = IFILE, : FORM = 'UNFORMATTED') DO 2 I = 4,NN IF (KT(I) .EQ. 0) GO TO 104 IF (CT(I) .EQ. 'NOHEAD') THEN LHEAD = .FALSE. ELSEIF (CT(I) .EQ. 'HEAD') THEN LHEAD = .TRUE. ELSEIF (CT(I) .EQ. 'OLDONLY') THEN LNEW = .FALSE. LOLD = .TRUE. ELSEIF (CT(I) .EQ. 'NEWONLY') THEN LOLD = .FALSE. LNEW = .TRUE. ELSEIF (CT(I) .EQ. 'NOCSF') THEN LNEW = .FALSE. LOLD = .FALSE. ELSEIF (CT(I) .EQ. 'CSFONLY') THEN LOLD = .TRUE. LNEW = .TRUE. LCOEF = .FALSE. ELSEIF (CT(I) .EQ. 'NOCOEFF') THEN LCOEF = .FALSE. ELSEIF (CT(I) .EQ. 'COEFF') THEN LCOEF = .TRUE. ELSEIF (CT(I) .EQ. 'ALL') THEN LOLD = .TRUE. LNEW = .TRUE. LCOEF = .TRUE. ELSEIF (CT(I) .EQ. 'NOPRINT') THEN LOLD = .FALSE. LNEW = .FALSE. LCOEF = .FALSE. ELSEIF (CT(I) .EQ. 'NOFILE') THEN LFILE = .FALSE. ELSEIF (CT(I) .EQ. 'FILE') THEN LFILE = .TRUE. ELSE GO TO 108 ENDIF 2 CONTINUE ENDIF ENDIF * IF (IOP .NE. 0) THEN NCFN = IT(1) NWN = IT(2) CALL DATNR (RT,IT,JT,LT,CT,KT,{C3},IOP,IFILE,LFILE) IF (IOP .GT. 1) ITC(21) = 1 IF (IOP .LT. 0) JTC(10) = 0 ELSE NMAN = IT(1) NWM = IT(2) NCFN = 0 CALL DATR (RT,IT,JT,LT,CT,KT,NMAN,NWM,{NW},{C3}) ENDIF * * READ A LABEL CARD AND CHECK WHETHER THE LABEL HAS BEEN USED * TRANSFER CONTROL TO THE RELEVANT SECTION FOR ANALYSIS OF * THE DATA AND READING OF FURTHER DATA (IF REQUIRED) * K = 0 3 KLAST = K CALL CARDIN (NN,RT,IT,JT,LT,CT,KT,{C3}) IF (KT(1) .EQ. 0) GO TO 104 IF (LT(1)) GO TO 107 NCH = MAX (KT(1),3) DO 4 K = 1,NLAB IF (CT(1)(1:NCH) .EQ. LABL(K)(1:NCH)) GO TO 5 4 CONTINUE GO TO 102 5 IF (LABUSD(K) .NE. 0) GO TO 103 IF ((IOP .LT. 0) .AND. (K .NE. 4) .AND. (K .NE. 7)) GO TO 109 LABUSD(K) = 1 GO TO (6,11,15,9,17,24,50,20,29,31,21,75), K * *----------------------------------------------------------------------* * * * MCP CARD * * * * FORMAT: STRING OF 3 CHARACTERS, 0 TO 3 INTEGERS * * * * LABEL "MCP" : IF FOUND, A CALL IS MADE TO THE * * MCP SECTION OF THE PROGRAM TO OBTAIN MCP COEFFICIENTS * * REQUIRED BY THE MCDF SECTION OF THE PROGRAM * * 1 INTEGER : [IO(1)] DATASTREAM NUMBER FOR OUTPUT * * OF UNSORTED MCP COEFFICIENTS; SET TO ZERO IF NO FILE IS * * REQUIRED. * * 1 INTEGER : [IO(2)] DATASTREAM NUMBER FOR OUTPUT * * OF SORTED MCP COEFFICIENTS; CARE SHOULD BE TAKEN THAT * * THIS NUMBER AGREES WITH THE DATASTREAM NUMBER USED TO * * READ THE MCP COEFFICIENTS IN THE MCDF PROGRAM IF THE * * SECTIONS ARE RUN TOGETHER; SET TO ZERO IF NO DISC/TAPE * * FILE REQUIRED * * 1 INTEGER : [IO(3)] DATASTREAM NUMBER FOR INPUT * * FILE OF MCP COEFFICIENTS; IF EQUAL TO OUTPUT FILE NUM- * * BER THEN A RESTART FROM AN INCOMPLETE CALCULATION IS * * ASSUMED; IF NONZERO OTHERWISE, INPUT IS TAKEN TO CON- * * TAIN A SUBSET OF THE MCP COEFFICIENTS REQUIRED: THESE * * ARE READ IN AND NOT CALCULATED * * * *----------------------------------------------------------------------* * 6 IF ((NN .LT. 1) .OR. (NN .GT. 4)) GO TO 101 JTC(12) = JTC(11) DO 7 I = 2,NN IF (KT(I) .NE. 0) GO TO 105 IF (LT(I)) GO TO 106 IO(I-1) = IT(I) 7 CONTINUE DO 8 I = 1,9 JTC(I) = 0 8 CONTINUE GO TO 3 * *----------------------------------------------------------------------* * * * MCDF CARD * * * * FORMAT: STRING OF 3 OR 4 CHARACTERS, 0 TO 30 INTEGERS * * * * LABEL "MCDF" (MINIMAL ABBREVIATION "MCD") : IF * * FOUND, A CALL IS MADE TO SUBROUTINE DATSCF TO READ DATA * * REQUIRED BY THE MCDF SECTION OF THE PROGRAM * * 0 TO 30 INTEGERS : OPTIONS; UP TO 30 MAY BE SPE- * * CIFIED: * * * * OPTION MEANING * * * * 1 ADJUST LAGRANGE MULTIPLIERS DURING ITERATION IN SOLVS * * 2 SET LAGRANGE MULTIPLIERS TO ZERO * * 3 SET LAGRANGE MULTIPLIERS BETWEEN ORBITALS WITH THE * * SAME GENERALISED OCCUPATION NUMBERS TO ZERO * * 4 WHEN LOADING FROM A DUMP, DO NOT RESCALE WITH Z * * 6 CHANGE MIXING OF OLD AND NEW ESTIMATES OF WAVEFUNC- * * TIONS ONLY WHEN LAST ITERATION WAS ON THE SAME FUNC- * * TION * * 7 DO NOT REDIAGONALIZE THE HAMILTONIAN DURING AN OL * * CALCULATION * * 9 COPY OLD ORBITALS FILE INTO NEW ORBITALS FILE * * 10 USE SLATER EXCHANGE POTENTIAL; SET AUTOMATICALLY BY * * DHFS CARD * * 11 WHEN LOADING A DUMP; READ ONLY THE WAVE FUNCTIONS AND * * PERFORM LOW-ACCURACY SCF CYCLES * * 11 WHEN LOADING FROM A DUMP, READ WAVEFUNCTIONS ONLY * * 12 STOP IF THE MAXIMUM NUMBER OF ITERATIONS IN SCF IS * * EXCEEDED * * 13 IN CASE OF A FAILURE IN SUBROUTINE SOLV, OPTIONS 14 * * 19 ARE ACTIVATED AND A RERUN OF SOLV IS PERFORMED; * * IF NONE OF THESE OPTIONS ARE SET INITIALLY, *ALL* * * WILL BE SET WHEN A FAILURE OCCURS * * 14 DEBUG PRINTOUTS IN SUBROUTINE OUT * * 15 DEBUG PRINTOUT OF INPUT PARAMETERS IN SUBROUTINE SOLV * * 16 PRINT WAVE FUNCTIONS AFTER EACH ITERATION IN SUBROU- * * TINE SOLV; *WARNING* THIS WILL PRODUCE A LARGE OUTPUT * * 17 DEBUG PRINTOUT IN SOLV DURING ITERATION CYCLES * * 18 DEBUG PRINTOUT IN SUBROUTINE OUT * * 19 PRINT POTENTIALS EACH TIME THEY ARE CALCULATED; *WAR- * * NING* THIS WILL PRODUCE A LARGE OUTPUT * * 20 SUPPRESS CALL TO MCDF ROUTINES * * 21 PRINT EIGENVECTORS IN LS BASIS; SET AUTOMATICALLY IF * * [IOP] > 1 OR LABELS PRESENT IN CONFIGURATION * * PARAMETERS AND OPTIONS CARD * * 22 SUPPRESS PRINTING OF EIGENVECTORS IN JJ-COUPLED BASIS * * 23 SUPPRESS PRINTING EIGENVECTORS IN ANY BASIS * * 24 PRINT ONLY THE FIRST 8 VECTORS IN SUBROUTINE MATOUT * * 25 DO NOT PRINT THE LARGEST 5 CONTRIBUTIONS FROM CSF TO * * ASF * * 26 PRINT COEFFICIENTS OF DIRECT AND EXCHANGE POTENTIALS * * IN MCPIN * * 27 WRITE OUT THE EIGENVALUES AND EIGENVECTORS OF THE * * HAMILTONIAN MATRIX WHENEVER IT IS CALCULATED * * 28 WRITE OUT THE GRID AND WAVEFUNCTIONS AT CONVERGENCE * * 29 WRITE OUT THE HAMILTONIAN MATRIX WHENEVER IT IS * * CALCULATED * * 30 PRINT OUT THE STATIC POTENTIAL FROM YPOTX * * * *----------------------------------------------------------------------* * 9 IF (NN .GT. 31) GO TO 101 JTC(13) = 1 IF (NN .GT. 1) THEN DO 10 I = 2,NN IF (KT(I) .NE. 0) GO TO 105 IF (LT(I)) GO TO 106 J = IT(I) IF ((J .LT. 1) .OR. (J .GT. 30)) GO TO 108 ITC(J) = 1 10 CONTINUE IF (ITC(20) .EQ. 1) JTC(13) = 0 ENDIF CALL DATSCF (RT,IT,JT,LT,CT,KT,{C3},IOP) GO TO 3 * *----------------------------------------------------------------------* * * * CONT CARD * * * * FORMAT : LABEL (COL.1-4) , 0 - 20 NUMBERS * * * * LABEL = 'CONT' - IF THIS LABEL IS FOUND THEN CONTINUUM ORBITALS * * ARE CALCULATED, FIXING THE CORE. * * * * THE LABEL IS FOLLOWED BY THE LOCATION OF WHERE THE GRID SWITCHES * * OVER FROM EXPONENTIAL TO LINEAR (DEFAULT=200), FOLLOWED BY THE * * FILE NUMBER TO WRITE THE CONTINUUM ORBITAL TO (DEFAULT=15). * * * *----------------------------------------------------------------------* * 75 IF (RT(2) .LT. 0.0D0) THEN WRITE (IPD,*) ' ERROR IN DATAIN--CONTINUUM E MUST BE > 0' STOP END IF ECONT = RT(2) MHALF = IT(3) IF (IT(4) .EQ. 0) THEN NUMWR = 15 ELSE NUMWR = IT(4) END IF IPHASE = IT(5) GO TO 3 * *----------------------------------------------------------------------* * * * MCT CARD 1 * * * * FORMAT: STRING OF 3 CHARACTERS, 0 TO 3 INTEGERS * * * * LABEL "MCT" : IF FOUND, A CALL IS MADE TO THE * * MCT SECTION OF THE PROGRAM TO OBTAIN MCT COEFFICIENTS. * * 1 INTEGER : [IO(4)]: FILE NUMBER FOR OUTPUT OF * * UNSORTED MCT COEFFICIENTS TO FILE; SET TO ZERO OR OMIT * * IF NO DISC/TAPE FILE REQUIRED. * * TOGETHER, ENSURE THAT THE FILE NUMBERS ARE COMPATIBLE. * * 1 INTEGER : [IO(5)]: FILE NUMBER FOR OUTPUT OF * * SORTED MCT COEFFICIENTS TO FILE; SET TO ZERO OR OMIT IF * * NO DISC/TAPE FILE REQUIRED; IF RUNNING MCT AND OSCL TO- * * GETHER, ENSURE THAT THE FILE NUMBERS ARE COMPATIBLE. * * 1 INTEGER : [IO(6)]: FILE NUMBER FOR INPUT OF * * UNSORTED MCT COEFFICIENTS TO FILE; SET TO ZERO OR OMIT * * IF NO DISC/TAPE FILE REQUIRED. IF RESTARTING UNFINISHED * * MCT RUN, SHOULD BE EQUAL TO IO(4). * * * *----------------------------------------------------------------------* * 11 IF ((NN .LT. 1) .OR. (NN .GT. 4)) GO TO 101 JTC(14) = JTC(11) DO 12 I = 2,NN IF (KT(I) .NE. 0) GO TO 105 IF (LT(I)) GO TO 106 IO(I+2) = IT(I) 12 CONTINUE * *----------------------------------------------------------------------* * * * MCT CARD 2 * * * * FORMAT: 1 TO NRP INTEGERS * * * * 1 TO NRP INTEGERS: [IRP(I), I = 1, * * NRP] LIST OF COMBINATIONS OF RANK AND PARITY; NEGATIVE * * VALUES INDICATE ODD PARITY, POSITIVE VALUES EVEN PARITY * * E.G., -1 = RANK 1, ODD PARITY; 2 = RANK 2, EVEN PARITY; * * IF BOTH PARITIES ARE REQUIRED FOR ALL TENSORS THE FIRST * * VALUE SHOULD BE 0, FOLLOWED BY A LIST OF RANKS ONLY; * * E.G., 0 1 WILL GIVE RANK 1, BOTH PARITIES * * * *----------------------------------------------------------------------* * CALL CARDIN (NN,RT,IT,JT,LT,CT,KT,{C3}) IF (IT(1) .NE. 0) THEN NKA = NN IF (NKA .GT. {NRP}) GO TO 101 DO 13 I = 1,NN IF (KT(I) .NE. 0) GO TO 105 IF (LT(I)) GO TO 106 IF (IT(I) .EQ. 0) GO TO 108 KA(I) = ABS (IT(I)) IOPAR(I) = SIGN (1,IT(I)) 13 CONTINUE ELSE IF (NN .EQ. 1) GO TO 101 NKA = NN-1 IF (NKA .GT. {NRP}) GO TO 101 DO 14 I = 2,NN IF (KT(I) .NE. 0) GO TO 105 IF (LT(I)) GO TO 106 IF (IT(I) .LE. 0) GO TO 108 KA(I-1) = IT(I) IOPAR(I-1) = 0 14 CONTINUE ENDIF GO TO 3 * *----------------------------------------------------------------------* * * * MCBP CARD * * * * FORMAT: STRING OF 3 OR 4 CHARACTERS, 1 TO 2 INTEGERS * * * * LABEL "MCBP" (MINIMAL ABBREVIATION "MCB") : IF * * FOUND, A CALL IS MADE TO THE MCBP SECTION OF THE PRO- * * GRAM TO OBTAIN MCBP COEFFICIENTS REQUIRED BY THE BENA * * SECTION OF THE PROGRAM * * 1 INTEGER : [IO(7)]: DATASTREAM NUMBER FOR * * OUTPUT OF THE MCBP COEFFICIENTS; SET TO ZERO IF NO DISC * * OR TAPE FILE REQUIRED; CARE SHOULD BE TAKEN THAT THIS * * NUMBER AGREES WITH THE DATASTREAM NUMBER USED TO READ * * THE MCBP COEFFICIENTS IN THE BENA PROGRAM IF THE PRO- * * GRAMS ARE RUN TOGETHER * * 1 INTEGER : [IO(8)] DATASTREAM NUMBER FOR INPUT * * OF MCBP COEFFICIENTS; IF EQUAL TO THE DATASTREAM NUMBER * * FOR OUTPUT, A RESTART FROM AN INCOMPLETE RUN IS ASSUMED * * * *----------------------------------------------------------------------* * 15 IF ((NN .LT. 2) .OR. (NN .GT. 3)) GO TO 101 JTC(15) = JTC(11) DO 16 I = 2,NN IF (KT(I) .NE. 0) GO TO 105 IF (LT(I)) GO TO 106 IO(I+5) = IT(I) 16 CONTINUE GO TO 3 * *----------------------------------------------------------------------* * * * BENA CARD 1 * * * * FORMAT: STRING OF 3 TO 4 CHARACTERS, 0 TO 30 INTEGERS * * * * LABEL "BENA" (MINIMAL ABBREVIATION "BEN") : IF * * FOUND, A CALL IS MADE TO THE BENA SECTION OF THE PRO- * * GRAM TO OBTAIN THE BREIT AND QED CORRECTIONS TO THE * * ENERGY LEVELS * * 0 TO 30 INTEGERS : OPTIONS: UP TO 30 CAN BE SPE- * * CIFIED: * * * * OPTION MEANING * * 1 DIAGONALIZE THE TRANSVERSE INTERACTION MATRIX ON ITS * * OWN * * 5 PRINT OUT CONTRIBUTIONS OF TRANSVERSE INTERACTION * * MATRIX ELEMENTS TO EIGENVALUES * * 7 SUPPLY ROUTINE FOR EVALUATING VACUUM POLARIZATION PO- * * TENTIAL IN SUBROUTINE VACUSR * * 8 USE POINT NUCLEUS VACUUM POLARIZATION POTENTIAL ONLY * * 9 USE SECOND-ORDER VACUUM POLARIZATION POTENTIAL ONLY * * 10 SUPPRESS THE FINAL SUMMARY PRINTING * * 11 GIVE INFORMATION ON BESSEL FUNCTION EVALUATION AND * * THE SOURCE OF THE FUNCTION * * 12 PRINT OUT TRANSVERSE INTERACTION MATRIX * * 13 PRINT RESTART INFORMATION * * 14 PRINT OUT EACH CONTRIBUTION TO THE TRANSVERSE INTER- * * ACTION MATRIX ELEMENTS * * 15 PRINT VACUUM POLARIZATION AND SELF-ENERGY INTEGRALS; * * PRINT QED CONTRIBUTIONS TO CSF * * 18 PRINT OUT MATRIX BEFORE DIAGONALIZATION * * 19 PRINT OUT VACUUM POLARIZATION POTENTIAL * * 20 SUPPRESS THE CALL TO THE BENA PACKAGE * * 21 PRINT EIGENVECTORS IN LS BASIS * * 22 DO NOT PRINT EIGENVECTORS IN JJ BASIS * * 23 DO NOT PRINT EIGENVECTORS IN ANY BASIS * * 24 PRINT EIGENVECTORS IN PREVIOUS ASF BASIS * * 25 DO NOT PRINT EIGENENERGIES * * 26 SUPPRESS PRINTING OF EIGENENERGIES AND EIGENVECTORS * * 27 PRINT SUMMARY IN RYDBERG UNITS * * 28 PRINT SUMMARY IN CM-1 UNITS * * 29 PRINT SUMMARY IN EV UNITS * * 30 SUPPRESS FINAL SUMMARY PRINTING * * * *----------------------------------------------------------------------* * 17 IF (NN .GT. 31) GO TO 101 JTC(16) = 1 IF (IOP .GT. 1) KTC(21) = 1 DO 18 I = 2,NN IF (KT(I) .NE. 0) GO TO 105 IF (LT(I)) GO TO 106 J = IT(I) IF ((J .LT. 1) .OR. (J .GT. 30)) GO TO 108 KTC(J) = 1 18 CONTINUE IF (KTC(20) .EQ. 1) JTC(16) = 0 IF (KTC(1)+KTC(2) .GT. 1) GO TO 114 * *----------------------------------------------------------------------* * * * BENA CARD 2 * * * * FORMAT: 1 TO 4 INTEGERS * * * * 1 INTEGER : [IO(9)]: DATASTREAM NUMBER OF FILE * * CONTAINING MCBP COEFFICIENTS * * 1 INTEGER : [IO(10)]: DATASTREAM NUMBER OF FILE * * CONTAINING WAVE FUNCTIONS * * 1 INTEGER : [IO(11)]: DATASTREAM NUMBER OF FILE * * TO WHICH SORTED INTEGRAL LIST AND INTEGRALS ARE TO BE * * WRITTEN * * 1 INTEGER : [IO(12)]: DATASTREAM NUMBER OF FILE * * FROM WHICH SORTED INTEGRAL LIST AND INTEGRALS ARE TO BE * * READ * * * *----------------------------------------------------------------------* * CALL CARDIN (NN,RT,IT,JT,LT,CT,KT,{C3}) IF ((NN .LT. 1) .OR. (NN .GT. 4)) GO TO 101 DO 19 I = 1,NN IF (KT(I) .NE. 0) GO TO 105 IF (LT(I)) GO TO 106 IO(I+8) = IT(I) 19 CONTINUE WFACT = ONE GO TO 3 * *----------------------------------------------------------------------* * * * LOWFREQ CARD * * * * FORMAT: STRING OF 3 TO 7 CHARACTERS, 1 REAL NUMBER * * * * LABEL "LOWFREQ" (MINIMAL ABBREVIATION "LOW") : * * IF FOUND, DATA FOR LOW FREQUENCY CASE FOLLOWS ON CARD. * * 1 REAL NUMBER : FACTOR MULTIPLYING FREQUENCY * * TO ENABLE THE LOW-FREQUENCY (LONG-WAVELENGTH) CASE TO * * BE OBTAINED; A FACTOR OF 1.E-3 IS USUALLY SUFFICIENT TO * * OBTAIN THE LOW-FREQUENCY LIMIT * * * *----------------------------------------------------------------------* * 20 IF (NN .NE. 2) GO TO 101 IF ((KLAST .NE. 5) .AND. (KLAST .NE. 11)) GO TO 110 IF (KT(2) .NE. 0) GO TO 105 WFACT = RT(2) GO TO 3 * *----------------------------------------------------------------------* * * * PRINT CARD * * * * SEE SUBROUTINE DATSCF FOR A DESCRIPTION. * * * *----------------------------------------------------------------------* * 21 IF ((KLAST .NE. 5) .AND. (KLAST .NE. 8)) GO TO 110 IF (NVEC .NE. 0) GO TO 103 IF ((NN .LT. 1) .OR. (NN .GT. NCF)) GO TO 101 IF (NN .GT. 1) THEN NVEC = NN DO 22 I = 2,NN IF (KT(I) .NE. 0) GO TO 105 IF (LT(I)) GO TO 106 IF ((IT(I) .LT. 1) .OR. (IT(I) .GT. NCF)) GO TO 108 IVEC(I) = IT(I) 22 CONTINUE ELSE NVEC = NCF DO 23 I = 1,NCF IVEC(I) = I 23 CONTINUE ENDIF GO TO 3 * *----------------------------------------------------------------------* * * * OSCL CARD 1 * * * * FORMAT: STRING OF 3 TO 4 CHARACTERS, 0 TO 20 INTEGERS * * * * LABEL "OSCL" (MINIMAL ABBREVIATION "OSC") : IF * * FOUND, A CALL IS MADE TO THE OSCILLATOR-STRENGTH PACK- * * AGE TO OBTAIN LEVEL TRANSITION OSCILLATOR STRENGTHS AND * * MATRIX ELEMENTS * * 0 TO 30 INTEGERS : OPTIONS; UP TO 30 CAN BE SPE- * * CIFIED: * * * * OPTION MEANING * * 1 USE COULOMB EIGENVECTORS AND EIGENVALUES ONLY * * 2 USE FULL HAMILTONIAN EIGENVECTORS AND EIGENVALUES * * ONLY * * 3 USE COULOMB EIGENVECTORS AND EIGENVALUES AND CYCLE * * BACK TO USE FULL HAMILTONIAN EIGENVECTORS AND * * EIGENVALUES * * 4 USE FULL HAMILTONIAN (WITHOUT QED) EIGENVECTORS AND * * EIGENVALUES AND CYCLE BACK TO USE EIGENVECTORS AND * * EIGENVALUES WITH QED * * 5 PRINT TRANSITION PROBABILITIES IN AU * * 6 SORT TRANSITIONS IN ASCENDING ORDER OF ENERGY * * 7 PRINT TRANSITION WAVELENGTHS IN ANGSTROM * * 8 PRINT TRANSITION ENERGIES IN CM-1 * * 9 PRINT TRANSITION ENERGIES IN EV * * 10 PRINT TRANSITION FREQUENCIES IN HZ * * 12 PRINT VALUE OF COULOMB AND GAUGE-DEPENDENT INTEGRALS * * 13 PRINT THE INTEGRAND OF THE GAUGE-DEPENDENT INTEGRALS * * 14 PRINT THE I AND J INTEGRALS * * 15 PRINT THE INTEGRANDS OF THE I AND J INTEGRALS * * 16 PRINT THE BESSEL FUNCTIONS * * 18 PRINT THE MCT COEFFICIENTS * * 19 WRITE OUT THE HAMILTONIAN EIGENVECTORS AND * * EIGENVALUES * * 20 SUPPRESS CALL TO THE OSCL PACKAGE * * * *----------------------------------------------------------------------* * 24 IF (NN .GT. 31) GO TO 101 JTC(17) = 1 DO 25 I = 2,NN IF (KT(I) .NE. 0) GO TO 105 IF (LT(I)) GO TO 106 J = IT(I) IF ((J .LT. 1) .OR. (J .GT. 30)) GO TO 101 LTC(J) = 1 25 CONTINUE IF (LTC(20) .EQ. 1) JTC(17) = 0 IF (LTC(7)+LTC(8)+LTC(9)+LTC(10) .GT. 1) GO TO 114 * *----------------------------------------------------------------------* * * * OSCL CARD 2 * * * * FORMAT: 1 TO 3 INTEGERS * * * * 1 INTEGER : [IO(13)]: DATASTREAM NUMBER OF * * FILE CONTAINING MCT COEFFICIENTS * * 1 INTEGER : [IO(14)]: DATASTREAM NUMBER OF FILE * * CONTAINING WAVEFUNCTIONS * * 1 INTEGER : [IO(15)]: DATASTREAM NUMBER OF FILE * * FOR OUTPUT OF TRANSITION PROBABILITIES * * * *----------------------------------------------------------------------* * CALL CARDIN (NN,RT,IT,JT,LT,CT,KT,{C3}) IF ((NN .LT. 1) .OR. (NN .GT. 3)) GO TO 101 DO 26 I = 1,NN IF (KT(I) .NE. 0) GO TO 105 IF (LT(I)) GO TO 106 IO(I+12) = IT(I) 26 CONTINUE IF (IO(14) .EQ. 0) THEN DO 27 I = 1,4 LTC(I) = 0 27 CONTINUE ENDIF NLEV = NCF DO 28 I = 1,NLEV LEVELS(I) = I 28 CONTINUE CUTOFF = ZERO GO TO 3 * *----------------------------------------------------------------------* * * * LEVELS CARD * * * * FORMAT: STRING OF 3 TO 6 CHARACTERS, 1 TO [NCF]-1 INTEGERS * * * * LABEL "LEVELS" (MINIMAL ABBREVIATION "LEV") : * * IF FOUND, THIS CARD SPECIFIES THE LEVELS BETWEEN WHICH * * OSCL CALCULATES DATA * * 2 TO [NCF]-1 INTEGERS : LIST OF * * LEVELS IN ASCENDING ORDER OF ENERGY * * * *----------------------------------------------------------------------* * 29 NLEV = NN-1 IF ((NN .LT. 3) .OR. (NLEV .GT. NCF)) GO TO 101 IF ((KLAST .NE. 6) .AND. (KLAST .NE. 10)) GO TO 111 DO 30 I = 2,NN IF (KT(I) .NE. 0) GO TO 105 IF (LT(I)) GO TO 106 J = IT(I) IF ((J .LT. 1) .OR. (J .GT. NCF)) GO TO 108 LEVELS(I-1) = J 30 CONTINUE GO TO 3 * *----------------------------------------------------------------------* * * * CUTOFF CARD * * * * FORMAT: STRING OF 3 TO 6 CHARACTERS, 1 REAL NUMBER * * * * LABEL "CUTOFF" (MINIMAL ABBREVIATION "CUT") : * * IF FOUND, CUTOFF VALUE FOR PRINTING OF TRANSITION RATES * * APPEARS ON CARD * * 1 REAL NUMBER : TRANSITIONS WITH RATES LESS * * THAN GIVEN VALUE WILL NOT BE PRINTED * * * *----------------------------------------------------------------------* * 31 IF (NN .NE. 2) GO TO 101 IF ((KLAST .NE. 6) .AND. (KLAST .NE. 9)) GO TO 111 IF (KT(2) .NE. 0) GO TO 105 IF (LT(2)) GO TO 106 CUTOFF = RT(2) GO TO 3 * *----------------------------------------------------------------------* * * * END CARD * * * * FORMAT: STRING OF 3 CHARACTERS, BLANK, STRINGS OF LENGTH LESS * * THAN OR EQUAL TO [STR] * * * * LABEL "END" : DENOTES THE END OF DATA FOR THIS * * CASE * * STRING : COMMENTS * * * *----------------------------------------------------------------------* * 50 IF ((JTC(16) .NE. 0) .AND. (IO(10) .EQ. 0) .AND. (JTC(13) .EQ. 0)) : GO TO 112 IF ((JTC(17) .NE. 0) .AND. (IO(14) .EQ. 0)) THEN IF (JTC(16) .NE. 0) THEN IO(14) = -1-KTC(2) ELSEIF (JTC(13) .EQ. 0) THEN GO TO 113 ENDIF ENDIF * IF (JTC(10) .EQ. 0) RETURN * IF (IOP .GT. 0) THEN IF (LOLD) CALL CFOUT IF (LNEW) CALL CSFOUT (LHEAD) IF (LCOEF) THEN DO 51 I = 1,NCF IT(I) = I 51 CONTINUE CALL MATOUT (TC,IT,NCF,NCF,{NC},{NC},2) ENDIF ELSE CALL CFOUT ENDIF * RETURN 99 SECNDS = MOD (TIME1,SIXTY) ISECND = INT (SECNDS*100.0{DB:D|:E}0) IHNDRD = MOD (ISECND,100) ISECND = INT (ISECND/100) IHOURS = INT (TIME1/SIXTY) MINUTS = MOD (IHOURS,ISIXTY) IHOURS = INT (IHOURS/ISIXTY) WRITE (IPD,302) IHOURS,MINUTS,ISECND,IHNDRD STOP * 101 WRITE (IPD,303) WRITE (IPD,304) STOP '020101' 102 WRITE (IPD,303) WRITE (IPD,305) STOP '020102' 103 WRITE (IPD,303) WRITE (IPD,306) STOP '020103' 104 WRITE (IPD,303) WRITE (IPD,307) STOP '020104' 105 WRITE (IPD,303) WRITE (IPD,308) STOP '020105' 106 WRITE (IPD,303) WRITE (IPD,309) STOP '020106' 107 WRITE (IPD,303) WRITE (IPD,310) STOP '020107' 108 WRITE (IPD,303) WRITE (IPD,311) STOP '020108' 109 WRITE (IPD,303) WRITE (IPD,312) STOP '020109' 110 WRITE (IPD,303) WRITE (IPD,313) STOP '020110' 111 WRITE (IPD,303) WRITE (IPD,314) STOP '020111' 112 WRITE (IPD,315) WRITE (IPD,316) STOP '020112' 113 WRITE (IPD,315) WRITE (IPD,317) STOP '020113' 114 WRITE (IPD,303) WRITE (IPD,318) STOP '020114' * 300 FORMAT (1A72) 301 FORMAT ('1++++++++++ DATAIN ++++++++++',//' CARD INPUT', : //1X,I3,' : ',A) 302 FORMAT (/' ***** E N D O F D A T A *****', : //' TOTAL CPU TIME USED = ',I4.2,2(':',I2.2),'.',I2.2) 303 FORMAT (/' ***** ERROR ON ABOVE DATA CARD *****') 304 FORMAT (/' TOO MANY OR TOO FEW VALUES ON CARD') 305 FORMAT (/' LABEL NOT FOUND') 306 FORMAT (/' LABEL ALREADY USED') 307 FORMAT (/' NUMBER GIVEN WHERE STRING EXPECTED') 308 FORMAT (/' STRING GIVEN WHERE NUMBER EXPECTED') 309 FORMAT (/' LOGICAL VALUE GIVEN WHERE NUMBER EXPECTED') 310 FORMAT (/' LOGICAL VALUE GIVEN WHERE STRING EXPECTED') 311 FORMAT (/' INCORRECT VALUE GIVEN') 312 FORMAT (/' ILLEGAL CALL WITH IOP < 0') 313 FORMAT (/' CARD MUST FOLLOW BENA CARD') 314 FORMAT (/' CARD MUST FOLLOW OSCL CARD') 315 FORMAT (/' ***** ERROR IN DATA INPUT *****') 316 FORMAT (/' BENA CALCULATION MUST BE PRECEDED BY MCDF CALCULATION', : ' WHEN DUMP FILE NUMBER IS ZERO') 317 FORMAT (/' OSCL CALCULATION MUST BE PRECEDED BY MCDF CALCULATION', : ' OR BENA CALCULATION WHEN DUMP FILE NUMBER IS ZERO') 318 FORMAT (/' INCONSISTENT OPTIONS') * END ********************SCF.NUM********************************************* ************************************************************************ * * * SECTION 10 S C F * * * * 15 ROUTINES * * * * 01. SCF 02. YPOT 03. XPOT * * 04. SXPOT 05. EIGEN 06. LAGR * * 07. HPOT 08. LINEQ 09. SOLV * * 10. SOLVS 11. OUT 12. IN * * 13. EXTEND 14. MATCH 15. ENDDIF * * * ************************************************************************ * * PREPROCESSOR SYMBOLS USED : * * SYMBOL USE * * LM NUMBER OF LAGRANGE MULTIPLIERS * NC NUMBER OF CONFIGURATIONS * NJ NUMBER OF ORBITAL J VALUES ALLOWED * NN NUMBER OF NODES ALLOWED IN HPOT * NP NUMBER OF GRID POINTS * NUC NUMBER OF NUCLEAR PARAMETERS * NW NUMBER OF ORBITALS * NX NUMBER OF COEFFICIENTS OF YK FUNCTIONS IN EXCHANGE POTL. * NY NUMBER OF COEFFICIENTS OF YK FUNCTIONS IN DIRECT POTENTIAL * N1 NP+10 * ORD MAXIMUM ORDER OF FINITE DIFFERENCE FORMULAE * V1 NW*(NW-1)/2+1 * W1 NW+1 * * DB FLAG FOR DOUBLE PRECISION VERSION * IMSL FLAG FOR USE OF IMSL SUBROUTINE LEQT1F * NAG FLAG FOR USE OF NAG SUBROUTINE F04ATF * ************************************************************************ ************************************************************************ * * SUBROUTINE SCF (ERAL,MX,DOALL) * * * ---------------- SECTION 10 SUBPROGRAM 01 ---------------- * * * * THIS SUBROUTINE CONTROLS THE MAIN SEQUENCE OF THE CALCULATION TO * * PRODUCE A FULL SET OF SELF-CONSISTENT WAVE FUNCTIONS. * * * * ERAL: MAXIMUM DISCREPANCY BETWEEN ESTIMATED AND CALCULATED * * WAVE FUNCTIONS * * MX : ERROR INDICATOR * * * * SUBROUTINES CALLED: DUMP, EIGEN, EXTAIL, HPOT, LAGR, RINT, QUAD, * * SOLV, SOLVS, SXPOT, TIMER, XPOT, YPOT * * * * LAST REVISION: 19 FEB 1988 * * * ************************************************************************ * {DB: IMPLICIT DOUBLEPRECISION (A-H, O-Z)} CHARACTER*2 NH LOGICAL DOALL * PARAMETER (EPS = 1.0{DB:D|:E}-10) * DIMENSION ISTOR(6) * COMMON/CONS/ZERO,HALF,TENTH,ONE,TWO,THREE,TEN : /DEF1/Z,ATW,RNT,H,EPH,ACCY,C,N : /DEF4/CHECK,NITIT,NSTEP,NSOLV,NSCF : /EXCO/PZ({NW}),QZ({NW}) : /FIXD/JFIX({NW}) : /GRID/R({NP}) : /INFORM/IRD,IPD : /INT2/P({NP}),Q({NP}),PC({NP}),QC({NP}) : /OFF0/ECV({LM}),IECC({LM}),NEC : /OFF1/ECO({NW}),IECO({NW}),NECO : /OPT1/ITC(30) : /ORB1/E({NW}),GAMA({NW}) COMMON/ORB2/IQ({NC},{NW}) : /ORB4/NW,NCF,NP({NW}),NAK({NW}) : /ORB5/NKL({NW}),NKJ({NW}),KAPMAX : /ORB7/XCAMAX({NW}),XCAINC({NW}),ISOLVS({NW}) : /ORB8/MCOW({NW}),COWMAN({NW}) : /ORB10/NH({NW}) : /SCF1/UCF({NW}) : /SCF6/MIX({NW}),DIFEIG({NW}),DIFPQ({NW}) : /SCF9/AP({NW}),AQ({NW}),EP({NW}),EQ({NW}),SP({NW}),SQ({NW}), : KP({NW}),KQ({NW}) COMMON/TIME/TIMECK,TIME1,TIME2 : /TATB/TA({N1}),TB({N1}) : /WAVE/PF({NP},{NW}),QF({NP},{NW}) : /POTE/YP({NP}),YQ({NP}),XP({NP}),XQ({NP}) COMMON/WAVEC/PFC({NPC},{NW}),QFC({NPC},{NW}) : /POTEC/YPC({NPC}),YQC({NPC}),XTPC({NPC}),XTQC({NPC}) : /GRIDC/RGRIDC({NPC}),ECONT,HCONT,MHALF,NCONT,NUMWR,IPHASE : /TATBC/TAC({NPC}),TBC({NPC}) PARAMETER (LTIMES=200) DIMENSION PC1({NP}),QC1({NP}),PC2({NPC}),QC2({NPC}),CC({NPC},3) * DIMENSION ECV1({LM}),ORTH1({LM}),ORTEM({LM},LTIMES) * MX = 0 ITR = NSCF J = 0 IF (ITC(13) .NE. 0) THEN NONE = 0 DO 1 I = 1,6 ISTOR(I) = ITC(I+13) NONE = NONE+ISTOR(I) ITC(I+13) = 0 1 CONTINUE IF (NONE .EQ. 0) THEN DO 2 I = 1,6 ISTOR(I) = 1 2 CONTINUE ENDIF ENDIF IF (ECONT .LT. ZERO) THEN IF (NEC .EQ. 0) THEN WRITE (IPD,300) ELSE WRITE (IPD,301) ENDIF ENDIF * * JL IS THE SERIAL NUMBER OF THE FUNCTION TREATED IN THE * PREVIOUS ITERATION. IN THE FIRST ITERATION IT IS ZERO * * DETERMINE WA, THE LARGEST DISCREPANCY, AND J, WHICH * SPECIFIES WHICH FUNCTION SHOWS THIS LARGEST VALUE * J = 0 4 WA = -ONE JL = J DO 5 I = 1,NW IF ((JFIX(I) .EQ. 0) .AND. (ECONT .GE. ZERO)) NP(I) = -1 IF (JFIX(I) .NE. 0) GO TO 5 IF (DOALL) THEN IF (I .LE. JL) GO TO 5 IF (J .GT. JL) GO TO 5 ELSE WB = ABS (DIFPQ(I)) IF (WB .LE. WA) GO TO 5 WA = WB ENDIF J = I 5 CONTINUE * * IF ECONT > 0, THEN CALCULATE CONTINUUM WAVEFUNCTIONS ASSUMING A * FIXED CORE. * IF (ECONT .GE. ZERO) THEN NELEC1 = 0 DO 36 I1 = 1,NW NELEC1 = NELEC1+IQ(1,I1) 36 CONTINUE * * DETERMINE IF LAGRANGE MULTIPLIERS ARE REQUIRED * NEC = 0 IF (ITC(2) .EQ. 1) GO TO 141 IF (NW .EQ. 1) GO TO 142 NWM = NW-1 NWA = NW+1 DO 139 JA = 1,NWM JP = JA+1 DO 138 JB = JP,NW IF (NAK(JA) .NE. NAK(JB)) GO TO 138 IF (NP(JA) .EQ. NP(JB)) GO TO 138 IF (ABS (UCF(JA)-FLOAT (NKJ(JA))-ONE) .GT. EPS) : GO TO 137 IF (ABS (UCF(JB)-FLOAT (NKJ(JB))-ONE) .GT. EPS) : GO TO 137 GO TO 138 137 IF (JFIX(JA)+JFIX(JB) .EQ. 2) GO TO 138 NEC = NEC+1 IF (NEC .GT. 60) GO TO 143 IECC(NEC) = JA+NWA*JB ECV(NEC) = ZERO 138 CONTINUE 139 CONTINUE IF (NEC .EQ. 0) GO TO 142 * WRITE (IPD,315) DO 140 I = 1,NEC MA = IECC(I)/NWA MB = IECC(I)-NWA*MA * WRITE (IPD,316) NP(MB),NH(MB),NP(MA),NH(MA) 140 CONTINUE GO TO 136 141 WRITE (IPD,317) GO TO 136 142 WRITE (IPD,318) GO TO 136 * 143 WRITE (IPD,319) STOP 5 * 136 WINECV = 0.10D0 DO 130 I=1,NEC ECV1(I) = ZERO ECV(I) = WINECV ORTH1(I) = ZERO 130 CONTINUE * * USE AUTOMATIC GRID FEATURE IF MHALF = ZERO * IF (MHALF .EQ. 0) THEN * * NPTS is the number of points per half cycle. * NPTS = 20 * * Estimate the value of HCONT, the spacing in the linear region * of the grid. * HCONT = (3.14159D0/NPTS)*SQRT(ONE/(TWO*ECONT)) IF (HCONT .GT. 0.200D0) HCONT=0.200D0 MHALF = INT((ONE/H)*LOG(HCONT/(RNT*(EXP(-H)-EXP(-TWO*H))))) * * Now, reset HCONT to make sure it is correct. * HCONT = R(MHALF)-R(MHALF-1) NCONT = INT((R(N)-R(MHALF))/HCONT + MHALF) IF (NCONT .GT. {NPC}) THEN NPC = {NPC} WRITE (IPD,320) NPC,NCONT STOP 5 END IF ELSE NCONT = {NPC} HCONT = (R(MHALF)-R(MHALF-1)) END IF * WRITE (IPD,*) ' HCONT=',HCONT * WRITE (IPD,*) ' R(MHALF)=',R(MHALF) RGRIDC(1) = R(1) DO 700 I = 2,NCONT IF (I .LE. MHALF) THEN RGRIDC(I) = R(I) ELSE RGRIDC(I) = RGRIDC(I-1) + HCONT END IF IF ((RGRIDC(I) .GT. R(N)) .AND. (N .LT. {NP}))N=N+1 700 CONTINUE DO 757 I=1,N IF (R(I) .GE. RGRIDC(NCONT)) GO TO 759 757 CONTINUE 759 N8=I-1 DO 710 I = 1,NW IF (JFIX(I) .EQ. 0) THEN DO 758 I1 = 1,NCONT PFC(I1,I) = ZERO QFC(I1,I) = ZERO 758 CONTINUE GO TO 710 END IF DO 720 I1 = 1,N PC1(I1) = PF(I1,I) QC1(I1) = QF(I1,I) 720 CONTINUE CALL SPLINE (R,PC1,N,CC) CALL SPLIN (R,PC1,N,CC,RGRIDC,PC2,NCONT-1) CALL SPLINE (R,QC1,N,CC) CALL SPLIN (R,QC1,N,CC,RGRIDC,QC2,NCONT-1) DO 730 I1 = 1,NCONT-1 PFC(I1,I) = PC2(I1) QFC(I1,I) = QC2(I1) 730 CONTINUE PFC(NCONT,I) = ZERO QFC(NCONT,I) = ZERO 710 CONTINUE END IF IF (DOALL .AND. (J .EQ. JL)) THEN DOALL = .FALSE. GO TO 4 ENDIF * * NORMAL EXIT IF REQUIRED ACCURACY HAS BEEN ACHIEVED * IF ((WA .LE. ERAL) .AND. .NOT.DOALL .AND. : (ECONT .LT. ZERO)) RETURN * * ITR COUNTS ITERATIONS DOWNWARDS FROM THE SPECIFIED MAXIMUM * CALL TIMER (1) IF (TIME2 .LT. TIMECK) GO TO 98 ITR = ITR-1 IF (ITR .LT. 0) GO TO 34 * * STORE ATOMIC NUMBER IN ZST. IPY INDICATES WHICH STAGE * THE AUTOMATIC RECOVERY PROCESS HAS REACHED. * IPY = 0 NORMAL COMPUTATION OF WAVE FUNCTION * IPY = 1 AFTER INITIAL FAILURE, ATOMIC NUMBER IS * BEING INCREASED OR EXCHANGE SET TO MINIMUM * IPY = 2 SOLUTION HAS BEEN FOUND WITH ARTIFICIALLY * INCREASED ATOMIC NUMBER OR DECREASED EXCHANGE, * WHICH IS NOW BEING REDUCED/INCREASED. * THIS REDUCTION OF THE ATOMIC NUMBER IN THE THIRD STAGE IS * ACHIEVED IN NSTEP STAGES, WHICH ARE COUNTED IN IPW, STEPPING * DOWN FROM NSTEP TO ZERO. * ZST = Z IPW = NSTEP IPY = 0 IPWL = NSTEP STEP = NSTEP XCA = XCAMAX(J) ZFAC = ONE * * XCA IS THE FACTOR MULTIPLYING EXCHANGE TERMS * Z IS THE EFFECTIVE ATOMIC NUMBER * * TABULATE THE DIRECT AND EXCHANGE POTENTIALS * 6 CONTINUE IF (ECONT .LT. ZERO) THEN CALL YPOT (J,VN,ZFAC) CALL XPOT (J,XCA) CALL LAGR (J,XCA) END IF IF ((ECONT .GT. ZERO) .AND. (JFIX(J) .EQ. 0)) THEN CS = C*C E(J) = -ECONT XMOMEN = SQRT(((ECONT+CS)**TWO-CS*CS)/CS) * WRITE (IPD,*) ' ENERGY=',-E(J),' K-VALUE=',XMOMEN WRITE (IPD,*) WRITE (IPD,*) WRITE (IPD,*) ' CALCULATION OF A CONTINUUM WITH SR CONTWVG' WRITE (IPD,*) KOUNT = 0 KOUNT2 = 0 VALUE = ZERO ISTOP = 0 ISTP = 1 ISTP3 = 0 ACCYSTP=1.0D-7 * ACCY2=5.0D-5 ACCY2=TEN*ACCYSTP XMIX = 0.7D0 INNER = 0 IF (XCA .LT. EPS) ITC(2) = 1 770 KOUNT = KOUNT + 1 FACTR = TEN**(TWO-DBLE(KOUNT2)) IF (FACTR .LT. ONE) FACTR=ONE IF (KOUNT .GT. 500) THEN WRITE (IPD,*) ' NUMBER OF ITERATIONS IN SCF EXCEEDED' STOP END IF DO 127 I1=1,NCONT TAC(I1)=PFC(I1,J)*PFC(I1,J)+QFC(I1,J)*QFC(I1,J) 127 CONTINUE CALL SIMP (TAC,H,1,MHALF,VALUE1,0) CALL SIMP (TAC,HCONT,MHALF,NCONT,VALUE2,0) VALUE = VALUE1+VALUE2 IF (KOUNT .EQ. 1) GO TO 743 IF (ABS(VALUE) .LT. 1.0E-20) THEN ISTOP = 1 GO TO 743 END IF IF (ABS((OLDVAL-VALUE)/VALUE) .LT. FACTR*ACCYSTP) ISTOP = 1 IF (ISTP .EQ. 0) THEN ISTOP = 0 ISTP = 1 END IF AG = PZ(J) 743 OLDVAL = VALUE IF (KOUNT .EQ. 2) THEN IF ((NINT(Z)-NELEC1) .NE. -1) THEN CALL NORMC2 (J,XX,PHASE,NELEC1) ELSE CALL NORMC (J,XX,PHASE) END IF AG = PZ(J)*XX END IF DO 783 I1 = 1,NCONT PC2(I1) = PFC(I1,J) QC2(I1) = QFC(I1,J) 783 CONTINUE * IF (ISTOP .EQ. 0) THEN CALL YPOT (J,VN,ZFAC) * * INTERPOLATE THE DIRECT POTENTIAL * CALL SPLINE (R,YP,N,CC) CALL SPLIN (R,YP,N,CC,RGRIDC,YPC,NCONT-1) YPC(NCONT) = YPC(NCONT-1) DO 773 I1 = 1,NCONT YQC(I1) = YPC(I1) 773 CONTINUE DO 650 I1=1,NCONT TAC(I1)=PFC(I1,J) TBC(I1)=QFC(I1,J) 650 CONTINUE CALL SPLINE (RGRIDC,TAC,NCONT,CC) CALL SPLIN (RGRIDC,TAC,NCONT,CC,R,PC1,N8) CALL SPLINE (RGRIDC,TBC,NCONT,CC) CALL SPLIN (RGRIDC,TBC,NCONT,CC,R,QC1,N8) DO 651 I1=1,N IF (I1 .LE. N8) THEN PF(I1,J)=PC1(I1) QF(I1,J)=QC1(I1) ELSE PF(I1,J)=ZERO QF(I1,J)=ZERO END IF 651 CONTINUE CALL XPOT (J,XCA) CALL SPLINE (R,XP,N,CC) CALL SPLIN (R,XP,N,CC,RGRIDC,XTPC,NCONT-1) CALL SPLINE (R,XQ,N,CC) CALL SPLIN (R,XQ,N,CC,RGRIDC,XTQC,NCONT-1) * * CLEAR FOR ACCUMULATION OF SUMS. * * ALX = ZERO * ALY = ZERO * * SUM OVER ALL ORBITALS, AVOIDING SERIAL NUMBER J. * DO 210 L = 1,NW IF ((JFIX(J) .EQ. 1) .AND. (JFIX(L) .EQ. 0)) GO TO 210 IF (L .EQ. J) GOTO 210 * * ADD CONTRIBUTIONS FROM OFF-DIAGONAL PARAMETERS. * IF (NEC .EQ. 0) GO TO 210 KA = NWA*J+L IF (J-L) 206,210,207 206 KA = NWA*L+J 207 DO 209 I = 1,NEC IF (IECC(I) .NE. KA) GOTO 209 WA = ECV(I)*XCA/(UCF(J)*C) IF (ABS (WA) .LT. EPS) GO TO 209 * * ADD CONTRIBUTION TO LEADING COEFFICIENTS. * * ALX = ALX+WA*QZ(L) * ALY = ALY-WA*PZ(L) * * ADD CONTRIBUTIONS TO EXCHANGE TERMS. * DO 220 K = 1,N XP(K) = XP(K)+R(K)*WA*QF(K,L) XQ(K) = XQ(K)-R(K)*WA*PF(K,L) 220 CONTINUE DO 208 K = 1,NCONT XTPC(K) = XTPC(K)+RGRIDC(K)*WA*QFC(K,L) XTQC(K) = XTQC(K)-RGRIDC(K)*WA*PFC(K,L) 208 CONTINUE 209 CONTINUE 210 CONTINUE CALL SOLVC (J,VN,XCA,XCP,XCQ,AG,WA,MF,FACTR*ACCYSTP, : KOUNT3) INNER=KOUNT3+INNER * DO 784 I1 = 1,NCONT PFC(I1,J) = XMIX*PFC(I1,J)+(ONE-XMIX)*PC2(I1) QFC(I1,J) = XMIX*QFC(I1,J)+(ONE-XMIX)*QC2(I1) 784 CONTINUE GO TO 770 END IF IF ((NEC .EQ. 0) .OR. (ITC(2) .EQ. 1) .OR. (ISTP3 .EQ. 1)) : GO TO 120 KOUNT2 = KOUNT2+1 IF (KOUNT2 .GT. LTIMES) THEN WRITE (IPD,*) ' INCREASE LTIMES IN SR SCF' STOP 2 END IF DO 119 LB = 1,NEC JA = IECC(LB)/NWA JB = IECC(LB)-NWA*JA IF (J-JA) 102,101,102 101 JJ = JB GO TO 104 102 IF (J-JB) 119,103,119 103 JJ = JA 104 CONTINUE ORTH = RINTC (J,JJ,0) ORTEM(LB,KOUNT2) = ORTH IF (KOUNT2 .EQ. 1) THEN ORTH1(LB) = ORTH ECV1(LB) = ECV(LB) ECV(LB) = -ECV(LB) GO TO 119 END IF SLOPE=(ORTH1(LB)-ORTH)/(ECV1(LB)-ECV(LB)) BINT=ORTH1(LB)-SLOPE*ECV1(LB) ECV1(LB)=ECV(LB) ECV(LB)=-BINT/SLOPE ORTH1(LB)=ORTH ORTEM(LB,KOUNT2)=ORTH 119 CONTINUE ISTOP = 0 ISTP = 0 IF (KOUNT2 .EQ. 1) GO TO 770 DO 125 I1=1,NEC IF (ABS(ORTEM(I1,KOUNT2)) .GT. ACCY2) GO TO 770 125 CONTINUE GO TO 155 120 CONTINUE * 155 WRITE (IPD,*) ' NUMBER OF ITERATIONS IN OUTER LOOP=',KOUNT-1 WRITE (IPD,*) ' NUMBER OF ITERATIONS IN INNER LOOP=',INNER IF ((NINT(Z)-NELEC1) .NE. -1) THEN CALL NORMC2 (J,XX,PHASE,NELEC1) ELSE CALL NORMC (J,XX,PHASE) END IF WRITE (IPD,*) IF (IPHASE .NE. 0) THEN WRITE (IPD,321) -E(J),HCONT,R(MHALF),XMOMEN,PHASE ELSE WRITE (IPD,322) -E(J),HCONT,R(MHALF),XMOMEN END IF CALL MAXMIN (J,1,PEAK) CALL ORTHG(0) 810 FORMAT (1X,11(1P,1E11.3)) IF (NUMWR .NE. 0) THEN WRITE (IPD,314) NH(J),NUMWR WRITE (NUMWR) NH(J),NP(J),NAK(J),NCONT,Z,H,RNT,E(J), : PZY,QZY,MHALF,HCONT WRITE (NUMWR) (PFC(L,J),L=1,NCONT), : (QFC(L,J),L=1,NCONT) * READ (IRD,*) I1,I2,I3,I4,K * CALL SLATCO(I1,I2,I3,I4,K,RK) * WRITE (IPD,*) ' RK integral:' * write (IPD,370) K,I1,I2,I3,I4,RK * 370 format(1x,' R',1i1,'(',1i2,',',1i2,';',1i2,',',1i2,') =',1p, * : 1d12.5) END IF RETURN END IF IF (XCA .LT. ZERO) CALL SXPOT IF ((MCOW(J) .EQ. 1) .AND. (XCA .NE. ZERO)) : CALL HPOT (J,COWMAN(J)) * * IF IN FINAL STAGE OF RECOVERY FROM PREVIOUS FAILURE, * EXTRAPOLATE ESTIMATES OF E AND PZ FROM PREVIOUS * TWO ITERATIONS * IF ((IPY .GE. 2) .AND. (IPWL .NE. IPW)) THEN AG = ACURR+ACURR-ALAST EG = ECURR+ECURR-ELAST ELSE * * USE ESTIMATE OF PZ FROM PREVIOUS ITERATION * AG = PZ(J) * * USE EIGENVALUE FROM PREVIOUS ITERATION, IF ANY * IF ((XCA .GT. EPS) .AND. (E(J) .GT. ZERO)) THEN EG = E(J) * * COMPUTE ESTIMATE OF EIGENVALUE, IF NO OTHER ESTIMATE IS AVAILABLE, * OR IF EXCHANGE TERMS ARE ZERO, OR IF DAMPING HAS BEEN USED. * IF ESTIMATE IS NOT POSITIVE, USE E = 0.1 * ELSE EG = EIGEN (J) IF (EG .LE. ZERO) EG = TENTH ENDIF ENDIF * * SOLVE DIRAC RADIAL EQUATIONS * IF (ISOLVS(J) .NE. 0) GO TO 29 CALL SOLV (J,VN,XCA,AG,QZE,EG,MF) IF (MF .LT. 0) GO TO 26 * * READJUST NUMBER OF GRID POINTS, IF NECESSARY * 7 IF (MF .EQ. N) GO TO 11 N = N+1 8 N = N-1 DO 9 M = 1,NW IF (ABS (PF(N,M)) .GT. EPS) GO TO 11 IF (ABS (QF(N,M)) .GT. EPS) GO TO 11 9 CONTINUE DO 10 M = 1,NW PF(N,M) = ZERO QF(N,M) = ZERO 10 CONTINUE GO TO 8 * * FIND LARGEST DIFFERENCE BETWEEN CALCULATED AND ESTIMATED * FUNCTIONS * 11 WC = ZERO DO 12 I = 1,MF WB = ABS (P(I)-PF(I,J))+ABS (Q(I)-QF(I,J)) IF (WB .GT. WC) WC = WB 12 CONTINUE DIFPQ(J) = WC WE = EG-E(J) * * AVOID ADJUSTMENT OF MIXING COEFFICIENT IF BETWEEN STEPS * IN RECOVERY PROCESS, OR OPTION 5 SET, OR OPTION 6 SET AND * LAST ORBITAL NOT THE SAME AS CURRENT, OR DIFFERENCE IS ZERO * IF ((IPY .EQ. 1) .OR. (IPWL .NE. IPW)) GO TO 13 IF ((ITC(5) .EQ. 1) .OR. ((ITC(6) .EQ. 1) .AND. (J .NE. JL))) : GO TO 13 IF ((DIFEIG(J) .EQ. ZERO) .OR. (WE .EQ. ZERO)) GO TO 13 * * ADJUST MIXING IF NEW CORRECTION IS MORE THAN HALF OLD CORRECTION: * DECREASE IF OF SAME SIGN * INCREASE IF OF OPPOSITE SIGN * INCREASE MIXING IF EIGENVALUE CHANGES BY MORE THAN 10 PERCENT * MIXJ = MIX(J) WER = TWO IF (ABS (WE/E(J)) .LE. TENTH) THEN WER = -WE/DIFEIG(J) IF (ABS (WER) .LE. HALF) GO TO 13 IF (WER .LT. ZERO) THEN MIXJ = MIXJ-1 IF (ABS (WER) .GE. ONE) MIXJ = MIXJ-1 IF (MIXJ .LT. 0) MIXJ = 0 GO TO 13 ENDIF ENDIF MIXJ = MIXJ+1 IF (ABS (WER) .GE. ONE) MIXJ = MIXJ+1 IF (MIXJ .GT. 9) MIXJ = 9 13 DIFEIG(J) = WE WB = MIXJ MIX(J) = MIXJ WB = WB*TENTH * * PRINT RECORD OF ITERATION RESULTS * IF (NECO .EQ. 0) THEN WRITE (IPD,302) J,NP(J),NH(J),EG,Z,WB,XCA,WC,WE ELSE M = IECO(1) DO 14 I = 1,N TA(I) = R(I)*(P(I)*PF(I,M)+Q(I)*QF(I,M)) 14 CONTINUE CALL QUAD (OVLP) ODLM = ECO(1)*UCF(J) WRITE (IPD,302) J,NP(J),NH(J),EG,Z,WB,XCA,WC,WE,NP(M),NH(M), : ODLM,OVLP DO 16 K = 2,NECO M = IECO(K) DO 15 I = 1,N TA(I) = R(I)*(P(I)*PF(I,M)+Q(I)*QF(I,M)) 15 CONTINUE CALL QUAD (OVLP) ODLM = ECO(K)*UCF(J) WRITE (IPD,303) NP(M),NH(M),ODLM,OVLP 16 CONTINUE ENDIF * * TABULATE NEW ESTIMATE AS COMBINATION OF OLD AND COMPUTED * FUNCTIONS * WC = ONE-WB DO 17 I = 1,N PF(I,J) = WC*P(I)+WB*PF(I,J) QF(I,J) = WC*Q(I)+WB*QF(I,J) 17 CONTINUE E(J) = EG PZ(J) = WB*PZ(J)+WC*AG QZ(J) = WB*QZ(J)+WC*QZE * * RE-ORTHOGONALIZE (IF NECESSARY) AND RE-NORMALIZE SOLUTION * DO 19 K = 1,J-1 IF ((ITC(8) .EQ. 0) .AND. (JFIX(K) .EQ. 0)) GO TO 19 IF (NAK(K) .NE. NAK(J)) GO TO 19 WA = RINT (J,K,0) IF (ABS (WA) .LT. EPS) GO TO 19 PZ(J) = PZ(J)-WA*PZ(K) QZ(J) = QZ(J)-WA*QZ(K) DO 18 I = 1,N PF(I,J) = PF(I,J)-WA*PF(I,K) QF(I,J) = QF(I,J)-WA*QF(I,K) 18 CONTINUE 19 CONTINUE CALL EXTAIL (J,MF) DO 20 I = MF,N PF(I,J) = SIGN (EXP (AP(J)+EP(J)*R(I)),SP(J)) QF(I,J) = SIGN (EXP (AQ(J)+EQ(J)*R(I)),SQ(J)) 20 CONTINUE WA = ONE/SQRT (RINT (J,J,0)) IF (PZ(J) .LT. ZERO) WA = -WA DO 21 I = 1,N PF(I,J) = PF(I,J)*WA QF(I,J) = QF(I,J)*WA 21 CONTINUE PZ(J) = PZ(J)*WA QZ(J) = QZ(J)*WA AP(J) = AP(J)*WA AQ(J) = AQ(J)*WA DO 23 K = J+1,NW IF ((ITC(8) .EQ. 0) .OR. (JFIX(K) .NE. 0)) GO TO 23 IF (NAK(K) .NE. NAK(J)) GO TO 23 WA = RINT (J,K,0) IF (ABS (WA) .LT. EPS) GO TO 23 IF (ITC(8) .EQ. 0) GO TO 23 WB = ONE-WA*WA WB = ONE/SQRT (WB+WB) PZ(K) = (PZ(K)-WA*PZ(J))*WB QZ(K) = (QZ(K)-WA*QZ(J))*WB DO 22 I = 1,N PF(I,K) = (PF(I,K)-WA*PF(I,J))*WB QF(I,K) = (QF(I,K)-WA*QF(I,J))*WB 22 CONTINUE 23 CONTINUE * * IF THIS COMPLETES STAGE 1 OF THE RECOVERY PROCESS, * BEGIN STAGE 2, COUNTING UPWARDS FROM ZERO IN IPW * IPWL = IPW IF (IPY-1) 33,24,25 24 IPY = 2 XINC = (XCAMAX(J)-XCA)/STEP ZDF = Z-ZST IPW = 0 IPWL = 0 ACURR = ZERO ECURR = ZERO 25 IF (IPW .GE. NSTEP) GO TO 33 * * STORE VALUES OF E AND PZ FOR SUBSEQUENT EXTRAPOLATION * ITR = ITR-1 * * REPEAT THIS STEP IF CHANGE WAS SIGNIFICANT * IF (ABS (DIFPQ(J)) .GT. 1.0D-2) GO TO 6 * * ... OTHERWISE REDUCE Z AND REPEAT * IPW = IPW+1 ALAST = ACURR ELAST = ECURR ACURR = PZ(J) ECURR = E(J) XCA = XCA+XINC Z = ZST+ZDF*DBLE (NSTEP-IPW) ZFAC = Z/ZST GO TO 6 * * IN CASE OF FAILURE, REPEAT COMPUTATION WITH EXTRA PRINTING * IF OPTION 13 IS SET * 26 IF (ITC(13) .NE. 0) THEN DO 27 I = 1,6 ITC(I+13) = ISTOR(I) 27 CONTINUE CALL SOLV (J,VN,XCA,AG,QZE,EG,MF) DO 28 I = 1,6 ITC(I+13) = 0 28 CONTINUE ENDIF * * ... AND USE ALTERNATIVE METHOD FOR SOLUTION * IF TWO SOLUTIONS, CHOOSE ONE WHOSE EIGENVALUE IS NEAREST PREVIOUS * 29 CALL SOLVS (J,VN,XCA,AG,QZE,EG,MF) IF (MF .GT. 0) GO TO 7 * * IF THIS FAILS, BEGIN RECOVERY BY DECREASING EXCHANGE TERMS * IF (IPY-1) 30,31,32 30 IPY = 1 JL = J 31 IF (XCA .GT. EPS) THEN XCA = XCA-XCAINC(J) IF (XCA .LT. ZERO) XCA = ZERO * * IF FAILURE WITH IPY = 1 INCREASE Z UNTIL CONVERGENCE * ELSE Z = Z+TENTH ZFAC = Z/ZST IF (Z .GT. ZST+5.0D0) THEN MF = 0 GO TO 32 ENDIF ENDIF GO TO 6 * * IF FAILURE WITH IPY = 2 ABANDON CALCULATION * 32 WRITE (IPD,307) NP(J),NH(J) IF (MF+1) 99,96,97 * * ON FINAL SUCCESS, RESET Z TO TRUE VALUE AND TRY ANOTHER * ITERATION * 33 Z = ZST GO TO 4 34 MX = 1 WRITE (IPD,305) IF (ITC(12) .EQ. 0) RETURN STOP '100101' * 96 WRITE (IPD,308) CALL TIMER (0) STOP '100102' 97 WRITE (IPD,309) CALL TIMER (0) STOP '100103' 98 CALL DUMP (0) WRITE (IPD,304) CALL TIMER (0) STOP '100104' 99 WRITE (IPD,306) CALL TIMER (0) STOP '100105' * 300 FORMAT (/'1++++++++++ SCF ++++++++++'//13X,'EIGENVALUE',8X, : 'Z MIX EXCHANGE WAVEFUNCTION EIGENVALUE',/43X,'FACTOR', : ' DIFFERENCE DIFFERENCE',/) 301 FORMAT (/'1++++++++++ SCF ++++++++++'//13X,'EIGENVALUE',8X, : 'Z MIX EXCHANGE WAVEFUNCTION EIGENVALUE ORB ', : 'LAGRANGE OVERLAP',/43X,'FACTOR DIFFERENCE DIFFERENCE', : 10X,'MULTIPLIER',/) 302 FORMAT (2I3,A2,1PD18.10,1X,0PF7.2,F5.1,1PD11.2,2D13.4, : 4X,I2,A2,2D12.4) 303 FORMAT (80X,I2,A2,1P2D12.4) 304 FORMAT (/' INSUFFICIENT TIME') 305 FORMAT (/' ITERATION LIMIT EXCEEDED') 306 FORMAT (/' INSUFFICIENT GRID POINTS') 307 FORMAT (/' FAILURE IN SCF PROCESS FOR ',I2,A2) 308 FORMAT (/' NUMBER OF ITERATIONS EXCEEDED') 309 FORMAT (/' Z INCREASED BEYOND REASONABLE LIMIT') 314 FORMAT (/1X,' WRITING CONTINUUM ORBITAL',1X,1A2,1X, : 'TO TAPE',1X,1I2,' USING THE FOLLOWING ORDER:'/, : 6X,'WRITE (NUMWR) NH(J),NP(J),NAK(J),NCONT,Z,H,RNT,E(J),',/, : 5X,':','PZY,QZY,MHALF,HCONT',/, : 6X,'WRITE (NUMWR) (PFC(L,J),L=1,NCONT),',/, : 5X,':','(QFC(L,J),L=1,NCONT)',//) 315 FORMAT (//37H INCLUDE LAGRANGE MULTIPLIERS BETWEEN/) 316 FORMAT (16X,I2,A2,2X,I2,A2) 317 FORMAT (//44H LAGRANGE MULTIPLIERS ARE NOT TO BE INCLUDED) 318 FORMAT (//38H LAGRANGE MULTIPLIERS ARE NOT REQUIRED) 319 FORMAT (/1X,23HDIMENSION ERROR FOR NEC) 320 FORMAT (/1X,'ERROR - INCREASE NUMBER OF GRID POINTS BY',/, : ' INCREASING THE VALUE OF NPC IN THE PRE-PROCESSOR FILE', : /,' NPC = ',1I4,2X,' NPC REQUIRED =',1I4) 321 FORMAT (/1X,' *** SYNOPSIS OF RESULTS FROM CONTWVG ***',/,/, : 4X,'ENERGY (A.U.)',6X,'HCONT',8X,'R(MHALF)', : 7X,'K-VALUE',7X,'PHASE(RADIANS)',/,1P, : 5X,1E10.4,5X,1E10.4,5X,1E10.4,5X,1E10.4,7X,1E10.4) 322 FORMAT (/1X,' SYNOPSIS OF RESULTS FROM CONTWVG',/,/, : 4X,'ENERGY (A.U.)',6X,'HCONT',8X,'R(MHALF)', : 7X,'K-VALUE',/,1P, : 5X,1E10.4,5X,1E10.4,5X,1E10.4,5X,1E10.4,5X,1E10.4) * END ************************************************************************ * * SUBROUTINE MAXMIN (J,IOUT,PEAK) * * * THIS SUBROUTINE FINDS THE MAXIMA AND MINIMA OF ARRAY PFC. * * * * INPUTS: * * J : SERIAL NUMBER OF ORBITAL * * IOUT : PRINTING OPTION * * * * OUTPUTS: * * PEAK : MAXIMUM AMPLITUDE OF PFC(I,J) * * * * NO SUBROUTINES CALLED. * * * ************************************************************************ * {DB: IMPLICIT DOUBLEPRECISION (A-H, O-Z)} * COMMON/CONS/ZERO,HALF,TENTH,ONE,TWO,THREE,TEN : /INFORM/IRD,IPD COMMON/WAVEC/PFC({NPC},{NW}),QFC({NPC},{NW}) : /GRIDC/RGRIDC({NPC}),ECONT,HCONT,MHALF,NCONT,NUMWR,IPHASE * IF (IOUT .EQ. 1) WRITE (IPD,300) IBIT = 1 PP2MAX = 0.0 DO 72 M = 2,NCONT-1 IF ((ABS(PFC(M,J)) .GT. PP2MAX) .OR. (ABS(PFC(M-1,J)) .GT. : PP2MAX)) THEN IF (IBIT .EQ. 1) PP2MAX = ABS(PFC(M,J)) IF (PFC(M-1,J)*PFC(M,J) .LT. 0.0) IBIT = 1 ELSE IF ((PFC(M-1,J) .GT. 0.0) .AND. (IBIT .EQ. 1)) THEN IF (IOUT .EQ. 1)WRITE (IPD,301) PFC(M-1,J),RGRIDC(M-1) IBIT = 0 PEAK = PFC(M-1,J) PP2MAX = 0.0 ELSEIF ((PFC(M-1,J) .LT. 0.0) .AND. (IBIT .EQ. 1)) : THEN IF (IOUT .EQ. 1)WRITE (IPD,302) PFC(M-1,J),RGRIDC(M-1) PEAK = ABS(PFC(M-1,J)) IBIT = 0 PP2MAX = 0.0 END IF END IF 72 CONTINUE RETURN 300 FORMAT (/,' MAXIMA',9X,'MINIMA',8X,'R') 301 FORMAT (1X,1P,1E12.6,16X,1E12.6) 302 FORMAT (13X,1P,2E14.6) END ************************************************************************ * * SUBROUTINE SPLINE (X,Y,N,C) * * * THIS ROUTINE IS AN ADAPTATION OF THE CUBIC SPLINES INTERPOLATION * * ROUTINE BY FORSYTHE ET AL., "COMPUTER METHODS FOR MATHEMATICAL * * COMPUTATIONS," (1977) BY W.F. PERGER, MICHIGAN TECH UNIV, DEC. 1988* * * * NO SUBROUTINES CALLED. * * * ************************************************************************ * {DB: IMPLICIT DOUBLEPRECISION (A-H, O-Z)} INTEGER N DIMENSION X(N),Y(N),C(N,3) * * THE COEFFICIENTS C(I,1), C(I,2), AND C(I,3), I=1,2,...N ARE COMPUTED * FOR A CUBIC INTERPOLATING SPLINE * * S(X)= Y(I) + C(I,1)*(X-X(I)) + C(I,2)*(X-X(I))**2 + C(I,3)*(X-X(I))**3 * * FOR X(I) .LE. X .LE. X(I+1) * * INPUT: * * N = THE NUMBER OF DATA POINTS OR KNOTS (N.GE.2) * X = THE ABSCISSAS OF THE KNOTS IN STRICTLY INCREASING ORDER * Y = THE ORDINATES OF THE KNOTS * * OUTPUT: * * C(I,J) = ARRAYS OF SPLINE COEFFICIENTS AS DEFINED ABOVE. * * USING P TO DENOTE DIFFERENTIATION * * Y(I) = S(X(I)) * C(I,1) = SP(X(I)) * C(I,2) = SPP(X(I))/2 * C(I,3) = SPPP(X(I))/6 (DERIVATIVE FROM THE RIGHT) * * THE ACCOMPANYING FUNCTION SUBPROGRAM SPLIN IS USED TO EVALUATE * THE SPLINE * INTEGER NM1,IB,I * NM1 = N-1 IF (N .LT. 2) RETURN IF (N .LT. 3) GO TO 50 * * SET UP TRIDIAGONAL SYSTEM * * C(I,1) = DIAGONAL, C(I,3) = OFFDIAGONAL, C(I,2) = RIGHT HAND SIDE * C(1,3) = X(2) - X(1) C(2,2) = (Y(2) - Y(1))/C(1,3) DO 10 I = 2,NM1 C(I,3) = X(I+1) - X(I) C(I,1) = 2.0D0*(C(I-1,3) + C(I,3)) C(I+1,2) = (Y(I+1) - Y(I))/C(I,3) C(I,2) = C(I+1,2) - C(I,2) 10 CONTINUE * * END CONDITIONS. THIRD DERIVATIVES AT X(1) AND X(N) * OBTAINED FROM DIVIDED DIFFERENCES. * C(1,1) = -C(1,3) C(N,1) = -C(N-1,3) C(1,2) = 0.0D0 C(N,2) = 0.0D0 IF (N .EQ. 3) GO TO 15 C(1,2) = C(3,2)/(X(4)-X(2)) - C(2,2)/(X(3)-X(1)) C(N,2) = C(N-1,2)/(X(N)-X(N-2)) - C(N-2,2)/(X(N-1)-X(N-3)) C(1,2) = C(1,2)*C(1,3)**2.0D0/(X(4)-X(1)) C(N,2) = -C(N,2)*C(N-1,3)**2.0D0/(X(N)-X(N-3)) * * FORWARD ELIMINATION * 15 DO 20 I=2,N T = C(I-1,3)/C(I-1,1) C(I,1) = C(I,1) - T*C(I-1,3) C(I,2) = C(I,2) -T*C(I-1,2) 20 CONTINUE * * BACK SUBSTITUTION * C(N,2) = C(N,2)/C(N,1) DO 30 IB =1,NM1 I = N-IB C(I,2) = (C(I,2)-C(I,3)*C(I+1,2))/C(I,1) 30 CONTINUE * * C(I,2) IS NOW THE SIGMA(I) OF THE TEXT * * COMPUTE POLYNOMIAL COEFFICIENTS * C(N,1) = (Y(N) - Y(NM1))/C(NM1,3)+C(NM1,3)*(C(NM1,2)+2.0D0*C(N,2)) DO 40 I = 1,NM1 C(I,1) = (Y(I+1) - Y(I))/C(I,3)-C(I,3)*(C(I+1,2)+2.0D0*C(I,2)) C(I,3) = (C(I+1,2)-C(I,2))/C(I,3) C(I,2) = 3.0D0*C(I,2) 40 CONTINUE C(N,2) = 3.0D0*C(N,2) C(N,3) = C(N-1,3) RETURN * 50 C(1,1) = (Y(2)-Y(1))/(X(2)-X(1)) C(1,2) = 0.0D0 C(1,3) = 0.0D0 C(2,1) = C(1,1) C(2,2) = 0.0D0 C(2,3) = 0.0D0 RETURN END ************************************************************************ * * SUBROUTINE SPLINT (X,Y,JMIN,JMAX,N,YY,RES,HX,RNT) * * * WRITTEN BY WARREN F. PERGER, MICHIGAN TECHNOLOGICAL UNIVERSITY, * * JULY, 1988 * * * * THIS SUBROUTINE IS AN ADAPTATION OF EQUATION 2.3.2.2, "METHODS OF * * NUMERICAL INTEGRATION," BY PHILIP J. DAVIS AND PHILIP RABINOWITZ, * * (1975). IN ADDITION, IT ALSO ESTIMATES THE INTEGRAL BETWEEN R=0 * * AND R=RNT (THE FIRST GRID POINT) BY ASSUMING THAT THE INTEGRAND HAS * * THE FORM R*F(R)=W*(R**PX), WHERE PX IS DETERMINED BY USING THE FIRST * * TWO POINTS IN THE INTEGRAND. THE PROCEDURE IS STARTS BY WRITING THE * * INTEGRAL FROM 0 TO RNT OF W*R**PX = (W*RNT**(PX+1))/(PX+1). THEN, * * F(I=1) = W*RNT**PX, F(I=2) = W*((RNT**PX)*(EXP(H*PX)). TAKING * * F(I=2)/F(I=1) YIELDS PX = (1/H)*LN(F(I=2)/F(I=1)). THE INTEGRAL * * OF F(R)DR FROM 0 TO RNT THEN BECOMES : * * F(I=1)*RNT/(1+(1/H)*LN(F(I=2)/F(I=1))) * * * * INPUTS: X THE ARRAY OF ABSCISSAS. * * Y THE ARRAY OF THE FUNCTION TO BE INTEGRATED AT THE * * LOCATIONS SPECIFIED BY X. * * YY THE ARRAY OF CUBIC SPLINES, WITH YY(I,3) BEING THE * * SECOND DERIVATIVE OF Y(X). * * JMIN THE LOCATION IN X OF THE BEGINNING OF THE INTERVAL. * * JMAX THE LOCATION IN X OF THE END OF THE INTERVAL OF INTEG. * * N THE NUMBER OF POINTS IN X AND Y. * * HX FROM THE EXPONENTIAL GRID, R(I) = RNT*EXP((I-1)*HX) * * RNT (SEE HX) * * * * OUTPUTS: RES THE RESULT OF THE INTEGRATION FROM JMIN TO JMAX * * * * NO SUBROUTINES CALLED. * * ************************************************************************ * {DB: IMPLICIT DOUBLEPRECISION (A-H, O-Z)} * DIMENSION X(N),Y(N),YY(N,3) RES = 0.0D0 DO 10 J = 1,N IF (ABS(Y(J)) .GT. 1.0D-30) GO TO 4 10 CONTINUE RETURN 4 JJ = J+1 DO 1 I = JMIN,JMAX-1 H = X(I+1)-X(I) RES = RES + (Y(I)+Y(I+1))*H/2.0D0 : - (YY(I,3)+YY(I+1,3))*H*H*H/24.D0 1 CONTINUE COP = Y(JJ)/Y(J) IF ((COP .GT. 1.0D0) .AND. (HX .GT. 1.0E-10)) GO TO 6 5 RPX = 0.0D0 * WRITE (IPD,300) GO TO 7 6 CONTINUE RPX = LOG(COP)/HX 7 RES = RES+Y(J)*RNT/(RPX+1.0D0) * * 300 FORMAT (1X,44HSPLINT : CORRECTION NEAR NUCLEUS SET TO ZERO) * RETURN END ************************************************************************ * * SUBROUTINE SPLIN (XA,YA,N,CC,XB,YB,NB) * * * WRITTEN BY WARREN F. PERGER, MICHIGAN TECHNOLOGICAL UNIVERSITY, * * DECEMBER, 1988 * * * * THIS SUBROUTINE IS AN ADAPTATION OF A ROUTINE FROM THE BOOK BY * * FORSYTHE, ET AL., "COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS." * * * ************************************************************************ * {DB: IMPLICIT DOUBLEPRECISION (A-H, O-Z)} * COMMON/INFORM/IRD,IPD DIMENSION XA(N),YA(N),CC(N,3),XB(NB),YB(NB) DATA I/1/ IF (XA(N) .LT. XB(NB)) WRITE (IPD,300) XA(N),XB(NB) DO 40 J=1,NB IF (I .GE. N) I = 1 X = XB(J) IF (X .LT. XA(I)) GO TO 10 IF (X .LE. XA(I+1)) GO TO 30 * * BINARY SEARCH * 10 I = 1 L = N+1 20 K = (I+L)/2 IF (X .LT. XA(K)) L = K IF (X .GE. XA(K)) I = K IF (L .GT. I+1) GO TO 20 * * EVALUATE SPLINE * 30 DX = X-XA(I) YB(J) = YA(I)+DX*(CC(I,1)+DX*(CC(I,2)+DX*CC(I,3))) 40 CONTINUE RETURN 300 FORMAT (/,1X,1P,'WARNING--ERROR IN SR SPLIN: XA(N)=',1X,1E11.3,/, : ' AND XB(NB)=',1X,1E11.3,//) END ************************************************************************ * * SUBROUTINE SPLIND (XA,YA,N,CC,XB,YB,NB) * * * WRITTEN BY WARREN F. PERGER, MICHIGAN TECHNOLOGICAL UNIVERSITY, * * DECEMBER, 1988 * * * * THIS SUBROUTINE IS AN ADAPTATION OF A ROUTINE FROM THE BOOK BY * * FORSYTHE, ET AL., "COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS." * * THE PURPOSE OF THIS ROUTINE IS TO TAKE THE DERIVATIVE, DYA/DXA, AND * * RETURN THE RESULT IN ARRAY YB USING CUBIC SPLINES IN ARRAY CC FROM A * * PREVIOUS CALL TO SUBROUTINE SPLINE. * * * * NO SUBROUTINES CALLED. * * * ************************************************************************ * {DB: IMPLICIT DOUBLEPRECISION (A-H, O-Z)} DIMENSION XA(N),YA(N),CC(N,3),XB(NB),YB(NB) DATA I/1/ DO 40 J=1,NB IF (I .GE. N) I = 1 X = XB(J) IF (X .LT. XA(I)) GO TO 10 IF (X .LE. XA(I+1)) GO TO 30 * * BINARY SEARCH * 10 I = 1 L = N+1 20 K = (I+L)/2 IF (X .LT. XA(K)) L = K IF (X .GE. XA(K)) I = K IF (L .GT. I+1) GO TO 20 * * EVALUATE DERIVATIVE * 30 DX = X-XA(I) YB(J) = CC(I,1)+DX*(2.0D0*CC(I,2)+3.0D0*DX*CC(I,3)) 40 CONTINUE RETURN END ************************************************************************ * * SUBROUTINE SOLVC (J,VN,XCF,ALX,ALY,AST,EST,MX,ACCYSTP,KOUNT) * * * THIS SUBROUTINE SOLVES A SINGLE PAIR OF DIRAC RADIAL EQUATIONS. * * THE DIRECT POTENTIALS ARE TABULATED IN YP AND YQ . * * THE EXCHANGE POTENTIALS ARE TABULATED IN XTP AND XTQ . * * THESE ARRAYS ARE OVERWRITTEN IN SOLV SO THEY MUST BE SET UP * * FOR THE NEXT CALL TO SOLV . * * * * INPUTS: * * * * J : SERIAL NUMBER OF WAVE FUNCTION * * VN : COEFFICIENT OF R IN SERIES EXPANSION OF DIRECT POTENTIAL * * XCF: SCALE FACTOR MULTIPLYING EXCHANGE TERM * * ALX: LEADING COEFFICIENT IN SERIES EXPANSION OF EXCHANGE * * POTENTIAL FOR P EQUATION * * ALY: LEADING COEFFICIENT IN SERIES EXPANSION OF EXCHANGE * * POTENTIAL FOR Q EQUATION * * AST: ESTIMATE OF A * * EST: ESTIMATE OF EIGENVALUE E * * * * OUTPUTS: * * * * MX : SHOULD BE ZERO. A NON-ZERO VALUE INDICATES * * A FAILURE. * * E(J) THE EIGENVALUE * * * * THE WAVE FUNCTIONS WILL BE TABULATED IN ARRAYS * * P AND Q, AND THE COEFFICIENTS IN THE SERIES * * EXPANSIONS WILL BE STORED IN PZ(J), QZ(J). * * * * SUBROUTINES CALLED: * * * * PRWF, OUT, IN, MATCH, QUAD. * * * ************************************************************************ * {DB: IMPLICIT DOUBLEPRECISION (A-H, O-Z)} CHARACTER*2 NH * COMMON/INFORM/IRD,IPD : /WAVE/PF({NP},{NW}),QF({NP},{NW}) : /CONS/ZERO,HALF,TENTH,ONE,TWO,THREE,TEN : /DEF1/Z,ATW,RNT,H,EPH,ACCY,C,N : /ORB2/IQ({NC},{NW}) : /ORB4/NW,NCF,NP({NW}),NAK({NW}) : /ORB10/NH({NW}) : /ORB1/E({NW}),GAMA({NW}) : /POTE/YP({NP}),YQ({NP}),XTP({NP}),XTQ({NP}) : /TATB/TA({N1}),TB({N1}) : /INT1/TC({NP}),TD({NP}),TE({NP}),TF({NP}), : TG({NP}),TH({NP}),TI({NP}),TJ({NP}) : /INT2/P({NP}),Q({NP}),PC({NP}),QC({NP}) : /INT3/XU({NP}),XV({NP}),XR({NP}),XS({NP}) : /EXCO/PZ({NW}),QZ({NW}) : /KNTR/ITC(20) COMMON/WAVEC/PFC({NPC},{NW}),QFC({NPC},{NW}) : /POTEC/YPC({NPC}),YQC({NPC}),XTPC({NPC}),XTQC({NPC}) : /GRIDC/RGRIDC({NPC}),ECONT,HCONT,MHALF,NCONT,NUMWR,IPHASE : /GRID/R({NP}) : /TATBC/TAC({NPC}),TBC({NPC}) DIMENSION XPTEM({NPC}),XQTEM({NPC}),YPTEM({NPC}),YQTEM({NPC}), : YP1({NP}),YQ1({NP}),PTEM({NP}),QTEM({NP}),PT({NP}),QT({NP}) DIMENSION CC(50,3) * IF (ITC(15) .EQ. 1) WRITE (IPD,303) NP(J),NH(J), : VN,XCF,ALX,ALY,AST,EST,Z,N CON = 0.3D0 MAU = 0 * * BEGIN FROM GIVEN ESTIMATES OF A AND E * ACL = AST ECL = E(J) * ELW = -ONE EHG = -ONE FKJ = FLOAT (NAK(J)) HH = -H*HALF HHCONT = -HCONT*HALF WA = -FKJ*HH WAC = ONE+WA WBC = ONE-WA XMOMEN = SQRT((-E(J)+C*C)*(-E(J)+C*C)-C*C*C*C)/C * * SAVE DIRECT & EXCH POTENTIALS FOR REPETITIVE CALLS TO OUT AND OUTC * DO 41 M = 1,N YP1(M) = YP(M) YQ1(M) = YQ(M) PT(M) = PF(M,J) QT(M) = QF(M,J) 41 CONTINUE DO 40 M = 1,NCONT XPTEM(M) = XTPC(M) XQTEM(M) = XTQC(M) 40 CONTINUE VALUE = ZERO KOUNT = 0 IX = 10 IF (XCF .LE. ZERO) GO TO 3 * * TABULATE MODIFIED EXCHANGE TERMS, INCLUDING DIFFERENCE * CORRECTIONS. * 42 KOUNT = KOUNT + 1 IF (KOUNT .GT. 50) THEN WRITE (IPD,*) ' NUMBER OF ITERATIONS EXCEEDED IN SOLVC' STOP END IF DO 45 M = 1,NCONT XTPC(M) = XPTEM(M) XTQC(M) = XQTEM(M) 45 CONTINUE DO 50 I = 1,MHALF+IX * P(I) = PFC(I,J) * Q(I) = QFC(I,J) PC(I) = PT(I) QC(I) = QT(I) YP(I) = YP1(I) YQ(I) = YQ1(I) * XU(I) = XTPC(I) * XV(I) = XTQC(I) XU(I) = XTP(I) XV(I) = XTQ(I) 50 CONTINUE * DO 55 I = -IX/2,IX/2 * PTEM(I+IX/2+1) = P(MHALF+I) * QTEM(I+IX/2+1) = Q(MHALF+I) * TD(I+IX/2+1) = XTPC(MHALF+I) * TE(I+IX/2+1) = XTQC(MHALF+I) * TH(I+IX/2+1) = RGRIDC(MHALF+I) * 55 CONTINUE * DO 60 I = -1,IX/2 * RINT = R(MHALF+I) * IF (I .LT. 0) RINT = R(MHALF)-HCONT * CALL EXTRAP (TH,PTEM,IX,RINT,PC(MHALF+I)) * CALL EXTRAP (TH,QTEM,IX,RINT,QC(MHALF+I)) * 60 CONTINUE * DO 62 I = 0,IX/2 * RINT = R(MHALF+I) * CALL EXTRAP (TH,TD,IX,RINT,XU(MHALF+I)) * CALL EXTRAP (TH,TE,IX,RINT,XV(MHALF+I)) * 62 CONTINUE 1 DO 2 I = 4,MHALF+IX/2 XU(I-2) = HH*(XU(I-1)+XU(I-2))+((PC(I-3)-PC(I)) : /12.0D0 -(PC(I-2)-PC(I-1))*0.25D0)*XCF XV(I-2) = HH*(XV(I-1)+XV(I-2))+((QC(I-3)-QC(I)) : /12.0D0 -(QC(I-2)-QC(I-1))*0.25D0)*XCF 2 CONTINUE DO 6 I = MHALF+2,NCONT H1 = RGRIDC(I-1) H2 = RGRIDC(I-2) XTPC(I-2) = HHCONT*(XTPC(I-1)/H1+XTPC(I-2)/H2)+((PFC(I-3,J) : -PFC(I,J))/12.0D0 -(PFC(I-2,J)-PFC(I-1,J))*0.25D0)*XCF XTQC(I-2) = HHCONT*(XTQC(I-1)/H1+XTQC(I-2)/H2)+((QFC(I-3,J) : -QFC(I,J))/12.0D0 -(QFC(I-2,J)-QFC(I-1,J))*0.25D0)*XCF 6 CONTINUE DO 142 I = 1,IX TH(I)=RGRIDC(NCONT-IX+I-2) TD(I)=XTPC(NCONT-IX+I-2) TE(I)=XTQC(NCONT-IX+I-2) 142 CONTINUE CALL EXTRAP (TH,TD,IX,RGRIDC(NCONT-1),XTPC(NCONT-1)) CALL EXTRAP (TH,TE,IX,RGRIDC(NCONT-1),XTQC(NCONT-1)) * * TABULATE MODIFIED DIRECT POTENTIALS. * 3 WD = HH/C DO 4 I = 1,MHALF+IX/2 YP(I) = WD*YP(I) YQ(I) = WD*YQ(I) 4 CONTINUE DO 15 I = MHALF+1,NCONT YPTEM(I) = HHCONT*YPC(I)/(RGRIDC(I)*C) YQTEM(I) = HHCONT*YQC(I)/(RGRIDC(I)*C) 15 CONTINUE * * TABULATE FUNCTIONS TF(R) AND TG(R). * WC = HH*ECL/C WB = -H*C-WC WB = WB*RNT WC = WC*RNT WCC = HHCONT*ECL/C WBC = -HCONT*C-WCC * * TABULATE FUNCTIONS TF(R) AND TG(R) FOR THE EXPONENTIAL GRID * DO 11 I = 1,MHALF+IX/2 TF(I) = WB+YP(I) TG(I) = WC-YQ(I) WB = WB*EPH WC = WC*EPH 11 CONTINUE DO 16 I = MHALF+1,NCONT TAC(I) = WBC+YPTEM(I) TBC(I) = WCC-YQTEM(I) 16 CONTINUE * * CHOOSE NEW ESTIMATE OF A, IF NEGATIVE. * IF (ACL) 12,12,13 12 CONTINUE IF (KOUNT .EQ. 1) ACL = 100.0*ABS (FKJ) * * PRINT ESTIMATES OF A AND E, IF OPTION 19 IS SET. * 13 IF (ITC(19) .EQ. 1) WRITE (IPD,302) NP(J),NH(J),ACL,ECL * * SPECIAL SECTION AT END OF SUBROUTINE FOR COMPUTATIONS * WITHOUT EXCHANGE. * * IF (XCF)39,39,14 * * OUTWARD INTEGRATION. * 14 CALL OUT (J,VN,ACL,QZE,ECL,MHALF+IX/2,JP,M,MAU,MAE) DO 22 I = 1,MHALF+IX/2 PFC(I,J) = P(I) QFC(I,J) = Q(I) PT(I) = P(I) QT(I) = Q(I) 22 CONTINUE DO 65 I = -IX/2,IX/2 PTEM(I+IX/2+1) = P(MHALF+I) QTEM(I+IX/2+1) = Q(MHALF+I) TH(I+IX/2+1) = R(MHALF+I) 65 CONTINUE RINT = RGRIDC(MHALF+1) * CALL EXTRAP (TH,PTEM,IX+1,RINT,PFC(MHALF+1,J)) * CALL EXTRAP (TH,QTEM,IX+1,RINT,QFC(MHALF+1,J)) CALL SPLINE(TH,PTEM,IX+1,CC) CALL SPLIN (TH,PTEM,IX+1,CC,RINT,PFC(MHALF+1,J),1) CALL SPLINE(TH,QTEM,IX+1,CC) CALL SPLIN (TH,QTEM,IX+1,CC,RINT,QFC(MHALF+1,J),1) CALL OUTC (J,VN,ACL,QZE,ALX,ALY,ECL,NCONT,JP,M,MAU,MAE) DO 127 I1=1,NCONT TAC(I1)=PFC(I1,J)*PFC(I1,J)+QFC(I1,J)*QFC(I1,J) 127 CONTINUE CALL SIMP (TAC,H,1,MHALF,VALUE1,0) CALL SIMP (TAC,HCONT,MHALF,NCONT,VALUE2,0) OLDVAL = VALUE VALUE = VALUE1+VALUE2 IF (ABS((OLDVAL-VALUE)/VALUE) .GT. ACCYSTP) GO TO 42 * WRITE (IPD,*) ' KOUNT=',KOUNT * * NORMAL RETURN, STORING RESULTS. * PZ(J) = ACL QZ(J) = QZE RETURN * 300 FORMAT (1X,6HSOLV17,2X,I2,A2,2X,1P,8D12.4) 301 FORMAT (1X,6HSOLV18,2X,I2,A2,2X,1P,8D12.4) 302 FORMAT (1X,6HSOLV19,2X,I2,A2,2X,1P,8D12.4) 303 FORMAT (1X,6HSOLV15,2X,I2,A2,2X,1P,7D12.4,I5) END ************************************************************************ * * SUBROUTINE NORMC2 (IORBIT,C0,SHIFT,NELEC) * * * THIS SUBROUTINE ENERGY NORMALIZES THE CONTINUUM WAVEFUNCTION, * * USING THE RELATIVISTIC WKB EXPRESSION OF ONG AND RUSSEK, PHYS. REV.* * V. 17, PG. 120, 1978. THIS EXPRESSION IS EQUATION (13) FOR THE * * CONSTANT A. NOTE THAT THEIR VALUE OF K = -KAPPA AND THEIR * * CONTINUUM STATE ENERGY, E, IS GREATER THAN MC**2. THEREFORE, TO * * CONVERT FROM THE GRASP ENERGY, WHICH IS LESS THAN ZERO FOR CONT- * * INUUM'S AND GREATER THAN ZERO FOR BOUND STATES, USE: * * E(ONG&RUSSEK) = -E(GRASP)+MC**2 * * * * INPUTS: * * IORBIT : SERIAL NUMBER OF ORBITAL * * PFC : LARGE COMPONENT OF WAVEFUNCTION * * YPC : ELECTRIC POTENTIAL * * * * OUTPUTS: * * PFC,QFC : AFTER RENORMALIZATION * * * * SUBROUTINES CALLED: SPLINE, SPLIND. * * * * WRITTEN BY WARREN F. PERGER, MICHIGAN TECHNOLOGICAL UNIVERSITY * * LAST UPDATED: 06 AUG 1992 * * * ************************************************************************ * {DB: IMPLICIT DOUBLEPRECISION (A-H, O-Z)} IMPLICIT INTEGER (I-N) * LOGICAL FIRST COMMON/CONS/ZERO,HALF,TENTH,ONE,TWO,THREE,TEN COMMON/DEF1/Z,ATW,RNT,H,EPH,ACCY,C,N : /INFORM/IRD,IPD : /WAVE/PF({NP},{NW}),QF({NP},{NW}) : /EXCO/PZ({NW}),QZ({NW}) : /INT2/P({NP}),Q({NP}),PC({NP}),QC({NP}) : /ORB1/E({NW}),GAMA({NW}) : /ORB2/IQ({NC},{NW}) : /ORB4/NW,NCF,NP({NW}),NAK({NW}) : /ORB5/NKL({NW}),NKJ({NW}),KAPMAX : /TATB/TA({N1}),TB({N1}) : /POTE/YP({NP}),YQ({NP}),XP({NP}),XQ({NP}) COMMON/WAVEC/PFC({NPC},{NW}),QFC({NPC},{NW}) : /POTEC/YPC({NPC}),YQC({NPC}),XTPC({NPC}),XTQC({NPC}) : /GRIDC/RGRIDC({NPC}),ECONT,HCONT,MHALF,NCONT,NUMWR,IPHASE PARAMETER (NFIT=21) * DIMENSION PC1(NFIT),PC3(NFIT),XA(NFIT) DIMENSION CC(NFIT,3) DIMENSION VR(NFIT),DPHIDR(NFIT) DIMENSION DPHI1(NFIT),DPHI2(NFIT),DVBDR(NFIT) DIMENSION TEMP(NFIT),TEMP1(NFIT),TEMP2(NFIT),TEMP3(NFIT) DATA FIRST/.TRUE./ * * WRITE (IPD,*) * WRITE (IPD,*) ' ENTERING NORMC2 SUBROUTINE' * * SET UP CONSTANTS. * PI = TWO*TWO*ATAN(ONE) KAPPA = NAK(IORBIT) C2 = C*C IF (-E(IORBIT) .GT. TEN) THEN XMOMEN = SQRT((-E(IORBIT)+C2)*(-E(IORBIT)+C2)-C2*C2)/C ELSE XMOMEN = SQRT(-TWO*E(IORBIT)) END IF ENERGY = -E(IORBIT)+C2 * CALL MAXMIN (IORBIT,1,PEAK) * * I IS THE GRID POINT TO USE FOR THE WKB EXPRESSION, * I1 IS THE FIRST GRID POINT OF THE CUBIC SPLINES ROUTINES, * NPTS IS THE NUMBER OF POINTS TO USE IN THE CUBIC SPLINES ROUTINES. * NNORM IS THE NUMBER OF POINTS TO USE IN THE FITTING TO COS AND SIN. * NNORM = 20 NPTS = NFIT-1 I1 = NCONT - 2*NPTS M = 2 PI = 4.0D0*DATAN (ONE) MSTEP = 1 C00 = ONE C11 = 0.01 ERR = 1.0D-3 NUMAX = 40 NUMIT = 0 NUMITX = 0 IFINAL = 0 50 K2 = 1 I = I1 + NPTS - 1 IF (NCONT-I .LT. NPTS) THEN WRITE (IPD,*) ' INCREASE NCONT FOR SR NORMC3' STOP 1 END IF 2 RPT = RGRIDC(I) J = IORBIT X1 = C2*(-DBLE(KAPPA))*(-DBLE(KAPPA)-ONE) * * SET UP THE FIRST ESTIMATE OF THE DERIVATIVE W.R.T. R OF PHI (DPHIDR) * ALSO, SET UP THE POTENTIAL, V(R) IN VR, AND THE ARRAY TEMP WHICH * CONTAINS THE QUANTITY (ENERGY-V+M*C*C) * * write (IPD,*) ' dphidr should be about= ',sqrt(xmomen*xmomen * : +two*ypc(i)/rpt) DO 5 K=-NPTS/2,NPTS/2 VR(K+NPTS/2+1) = -YPC(I+K)/RGRIDC(I+K) IF (-E(IORBIT) .GT. TEN) THEN X2 = (ENERGY-VR(K+NPTS/2+1))**TWO - C2*C2 - : X1/(RGRIDC(I+K)*RGRIDC(I+K)) ELSE X2 = (XMOMEN*XMOMEN - TWO*VR(K+NPTS/2+1))*C2 - : X1/(RGRIDC(I+K)*RGRIDC(I+K)) END IF DPHIDR(K+NPTS/2+1) = SQRT(X2/C2) * write (IPD,*) ' dphidr=',dphidr(k+npts/2+1) TEMP(K+NPTS/2+1) = ENERGY-VR(K+NPTS/2+1)+C2 TEMP1(K+NPTS/2+1) = TEMP(K+NPTS/2+1)**(-HALF) XA(K+NPTS/2+1) = RGRIDC(I+K) 5 CONTINUE * * TAKE FIRST DERIVATIVE W.R.T. R OF (ENERGY-V+MC*C)**-1/2 * CALL SPLINE (XA,TEMP1,NPTS,CC) CALL SPLIND (XA,TEMP1,NPTS,CC,XA,TEMP3,NPTS) * * TAKE SECOND DERIVATIVE W.R.T. R OF (ENERGY-V+MC*C)**-1/2 * CALL SPLINE (XA,TEMP3,NPTS,CC) CALL SPLIND (XA,TEMP3,NPTS,CC,XA,TEMP2,NPTS) * * CALCULATE THE DERIVATIVE OF V(R) W.R.T. R * CALL SPLINE (XA,VR,NPTS,CC) CALL SPLIND (XA,VR,NPTS,CC,XA,DVBDR,NPTS) * write (IPD,*) ' dvbdr=',dvbdr(NPTS/2+1) * write (IPD,*) ' dV/dr should be about=',(Z-nelec+1)/(rpt*rpt) * * BEGIN LOOPING FOR ITERATION (TWO TIMES) OF THE DERIVATIVE OF * PHI W.R.T. R * DO 7 IK=1,2 IF (IK .EQ. 1) THEN * * CALCULATE FIRST DERIVATIVE OF DPHIDR**-1/2 * DO 6 K=-NPTS/2,NPTS/2 PC1(K+NPTS/2+1) = DPHIDR(K+NPTS/2+1)**(-HALF) 6 CONTINUE CALL SPLINE (XA,PC1,NPTS,CC) CALL SPLIND (XA,PC1,NPTS,CC,XA,DPHI1,NPTS) * * CALCULATE SECOND DERIVATIVE OF DPHIDR**-1/2 * CALL SPLINE (XA,DPHI1,NPTS,CC) CALL SPLIND (XA,DPHI1,NPTS,CC,XA,DPHI2,NPTS) END IF DO 9 K=-NPTS/2,NPTS/2 X2 = (ENERGY-VR(K+NPTS/2+1))**TWO - C2*C2 - X1/ : (XA(K+NPTS/2+1)*XA(K+NPTS/2+1)) DPHIDR(K+NPTS/2+1) = SQRT((X2/C2) + SQRT(DPHIDR(K+NPTS/2+1))* : DPHI2(K+NPTS/2+1) - SQRT(TEMP(K+NPTS/2+1))* : TEMP2(K+NPTS/2+1) + (KAPPA/XA(K+NPTS/2+1))* : DVBDR(K+NPTS/2+1)/TEMP(K+NPTS/2+1)) * write (IPD,*) 'dphidr(k+npts/2+1)= ',dphidr(k+npts/2+1) 9 CONTINUE 7 CONTINUE * * CALCULATE THE DERIVATIVE OF P(R) USING CUBIC SPLINES * DO 8 K=-NPTS/2,NPTS/2 PC1(K+NPTS/2+1) = PFC(I+K,J) IF (.NOT. FIRST .AND. (IPHASE .EQ. 1)) : PC3(K+NPTS/2+1) = HYD(I+K,J) 8 CONTINUE HYDI = PC3(NPTS/2+1) CALL SPLINE (XA,PC1,NPTS,CC) CALL SPLIND (XA,PC1,NPTS,CC,RGRIDC(I),DPBDR,1) IF (.NOT. FIRST .AND. (IPHASE .EQ. 1)) THEN CALL SPLINE (XA,PC3,NPTS,CC) CALL SPLIND (XA,PC3,NPTS,CC,RGRIDC(I),DHYDBDR,1) END IF DO 15 K=-NPTS/2,NPTS/2 PC1(K+NPTS/2+1) = DPHIDR(K+NPTS/2+1) 15 CONTINUE CALL SPLINE (XA,PC1,NPTS,CC) CALL SPLIND (XA,PC1,NPTS,CC,RGRIDC(I),DPHI22,1) * * CALCULATE U AT R0 * U = (ONE/DPHIDR(NPTS/2+1))*(-DPBDR - HALF*DVBDR(NPTS/2+1)* : PFC(I,J)/TEMP(NPTS/2+1) : - HALF*DPHI22*PFC(I,J)/DPHIDR(NPTS/2+1)) IF (.NOT. FIRST .AND. (IPHASE .EQ. 1)) THEN U2 = (ONE/DPHIDR(NPTS/2+1))*(-DHYDBDR - HALF*DVBDR(NPTS/2+1)* : HYDI/TEMP(NPTS/2+1) : - HALF*DPHI22*HYDI/DPHIDR(NPTS/2+1)) END IF * write (IPD,*) ' u=',u * * CALCULATE THE RENORMALIZATION CONSTANT, C0. * C0 = C*SQRT(PI*DPHIDR(NPTS/2+1)*(PFC(I,J)*PFC(I,J)+U*U)/ : TEMP(NPTS/2+1)) C0 = ONE/C0 DO 100 K=1,NCONT PFC(K,IORBIT) = PFC(K,IORBIT)*C0 QFC(K,IORBIT) = QFC(K,IORBIT)*C0 100 CONTINUE IF (FIRST .OR. (IPHASE .EQ. 0)) THEN FIRST = .FALSE. RETURN END IF * * CALCULATE THE PHASE. * * C01 = C*SQRT(PI*DPHIDR(NPTS/2+1)*(HYDI*HYDI+U2*U2)/ * : TEMP(NPTS/2+1)) * C01 = ONE/C01 PHIR0 = ACOS(ONE/SQRT(ONE + C0*C0*U*U/(PFC(I,J)*PFC(I,J)))) PHSHYD = ACOS(ONE/SQRT(ONE + U2*U2/(HYDI*HYDI))) PHIR0 = ATAN2 (U,PFC(I,J)/C0) PHSHYD = ATAN2 (U2,HYDI) write (IPD,*) SHIFT = PHIR0-PHSHYD IF (ABS(SHIFT) .GT. PI) THEN IF (SHIFT .LT. ZERO) THEN SHIFT = SHIFT+TWO*PI ELSE SHIFT = SHIFT-TWO*PI END IF END IF IF (NAK(IORBIT) .GT. 0) SHIFT=SHIFT-NINT(SHIFT/PI)*PI * WRITE (IPD,*) ' Phase shift, in degrees =',SHIFT*180.0D0/PI * WRITE (IPD,*) ' Phase shift, in radians=', SHIFT * WRITE (IPD,*) RETURN END ************************************************************************ * * SUBROUTINE EXTRAP(X1,Y1,NEXT,X2,Y2) * * * THIS SUBROUTINE PERFORMS A RATIONAL FUNCTION EXTRAPOLATION. * * * * INPUTS: X1,Y1 (ARRAYS CONTAINING X AND Y); NEXT (NUMBER OF POINTS * * TO USE IN EXTRAPOLATION PROCESS); X2 (INDEPENDENT VARIABLE * * TO EXTRAPOLATE). * * * * OUTPUTS: Y2 (EXTRAPOLATED VALUE AT X2) * * * * NO SUBROUTINES CALLED. * * * ************************************************************************ * {DB: IMPLICIT DOUBLEPRECISION (A-H, O-Z)} DIMENSION X1(NEXT),Y1(NEXT),C({NP}),D({NP}) EPS = 1.0E-25 NN=1 HH=ABS(X2-X1(1)) DO 10 I=1,NEXT H=ABS(X2-X1(I)) IF (H.EQ.0.0D0)THEN Y2=Y1(I) DELY=0.0D0 RETURN ELSE IF (H.LT.HH) THEN NN=I HH=H ENDIF C(I)=Y1(I) D(I)=Y1(I)+EPS 10 CONTINUE Y2=Y1(NN) NN=NN-1 DO 20 M=1,NEXT-1 DO 30 I=1,NEXT-M W=C(I+1)-D(I) H=X1(I+M)-X2 T=(X1(I)-X2)*D(I)/H DELD=T-C(I+1) DELD=W/DELD D(I)=C(I+1)*DELD C(I)=T*DELD 30 CONTINUE IF (2*NN .LT. NEXT-M) THEN DELY=C(NN+1) ELSE DELY=D(NN) NN=NN-1 ENDIF Y2=Y2+DELY 20 CONTINUE RETURN END ************************************************************************ * * SUBROUTINE OUTC (J,VN,AST,QZE,ALX,ALY,ECL,JC,JP,M,MAW,MA) * * * THIS SUBROUTINE CARRIES OUT THE STEP-BY-STEP OUTWARD * * INTEGRATION FOR THE INHOMOGENEOUS PAIR OF DIRAC RADIAL * * EQUATIONS. * * * * INPUTS: * * * * J : SERIAL NUMBER OF FUNCTION TO BE COMPUTED * * VN : COEFFICIENT OF R IN SERIES EXPANSION OF DIRECT POTENTIAL * * ALX: LEADING COEFFICIENTS IN THE SERIES EXPANSIONS * * ALY OF THE TWO EXCHANGE TERMS * * AST: CURRENT ESTIMATE OF THE PARAMETER A * * ECL: CURRENT ESTIMATE OF THE EIGENVALUE E * * JC : IF THIS IS NON-ZERO, IT INDICATES THE POINT * * AT WHICH THE STEP-BY-STEP INTEGRATION IS TO * * TERMINATE. IF JC IS ZERO, THIS POINT IS * * DETERMINED IN THE COURSE OF THE CALCULATION. * * MAW: IF THIS IS NON-ZERO, IT INDICATES THE POINT AT * * WHICH THE USE OF THE SERIES EXPANSION IS TO STOP, * * AND THE STEP-BY-STEP INTEGRATION IS TO BEGIN. * * * * OUTPUTS: * * * * M : IS A FAILURE WARNING, AND SHOULD NORMALLY BE ZERO: * * M = -1 INDICATES THAT P(R) HAS PASSED THROUGH * * TOO MANY ZEROES; * * M = +1 INDICATES THAT EITHER THE END OF THE TABLE * * HAS BEEN REACHED, OR THAT P(R) HAS * * EXCEEDED 10**6 IN MAGNITUDE, OR THAT P(R) * * DID NOT PASS THROUGH ENOUGH ZEROES. * * JP : INDICATES THE POINT AT WHICH THE INTEGRATION * * TERMINATED. * * MA : INDICATES THE LAST POINT AT WHICH THE SERIES * * EXPANSION WAS USED. * * QZE: LEADING COEFFICIENT IN EXPANSION OF Q(R) * * P : TWO COMMON ARRAYS IN WHICH THE REQUIRED * * Q SOLUTIONS ARE TABULATED * * * * NO SUBROUTINES CALLED. * * * ************************************************************************ * {DB: IMPLICIT DOUBLEPRECISION (A-H, O-Z)} CHARACTER*2 NH * COMMON/INFORM/IRD,IPD : /KNTR/ITC(20) : /CONS/ZERO,HALF,TENTH,ONE,TWO,THREE,TEN : /DEF1/Z,ATW,RNT,H,EPH,ACCY,C,N : /ORB2/IQ({NC},{NW}) : /ORB4/NW,NCF,NP({NW}),NAK({NW}) : /ORB10/NH({NW}) : /ORB1/E({NW}),GAMA({NW}) : /INT1/TC({NP}),TD({NP}),TE({NP}),TF({NP}), : TG({NP}),TH({NP}),TI({NP}),TJ({NP}) : /INT2/P({NP}),Q({NP}),PC({NP}),QC({NP}) : /INT3/XU({NP}),XV({NP}),XR({NP}),XS({NP}) COMMON/WAVEC/PFC({NPC},{NW}),QFC({NPC},{NW}) : /GRIDC/RGRIDC({NPC}),ECONT,HCONT,MHALF,NCONT,NUMWR,IPHASE : /POTEC/YPC({NPC}),YQC({NPC}),XTPC({NPC}),XTQC({NPC}) : /TATBC/TAC({NPC}),TBC({NPC}) * QUART = 2.5E-01 KJ = NAK(J) FKJ = FLOAT (KJ) * * * DEBUG PRINT OUT IF OPTION 15 1S SET * IF (ITC(15) .EQ. 0) GO TO 16 WRITE (IPD,300) NP(J),NH(J),WC,VN,ALX,ALY,PZE,QZE,PUN,QUN DO 15 I = 1,MA,10 WRITE (IPD,300) NP(J),NH(J),P(I),Q(I) 15 CONTINUE * 16 WA = HCONT*FKJ*HALF * * TABULATE P(R) AND Q(R) BY STEP-BY-STEP INTEGRATION * I = MHALF+1 19 I = I+1 WA1 = WA/RGRIDC(I-1) WA2 = WA/RGRIDC(I) WC = XTPC(I-1)-(WA1-ONE)*PFC(I-1,J)-TAC(I-1)*QFC(I-1,J) WD = XTQC(I-1)+(WA1+ONE)*QFC(I-1,J)-TBC(I-1)*PFC(I-1,J) WE = TAC(I)*TBC(I)+(WA2+ONE)*(WA2-ONE) PFC(I,J) = (TAC(I)*WD+(WA2-ONE)*WC)/WE QFC(I,J) = (TBC(I)*WC-(WA2+ONE)*WD)/WE * * IF JC IS NON-ZERO, STOP WHEN POINT JC IS REACHED. * IF (JC) 23,23,20 20 M = 0 IF (I-JC) 19,21,21 21 JP = I IF (ITC(15) .EQ. 1) WRITE (IPD,301) NP(J),NH(J),JS,M IF (ITC(14) .EQ. 0) RETURN * * DEBUG PRINT OUT IF OPTION 14 IS SET * II = MA+1 DO 22 I = II,JP,10 WRITE (IPD,302) NP(J),NH(J),P(I),Q(I),TAC(I),TBC(I),XTPC(I-1), : XTQC(I-1) 22 CONTINUE 23 RETURN * 300 FORMAT (1X,5HOUT15,2X,I2,A2,2X,1P,8D12.4) 301 FORMAT (1X,5HOUT15,2X,I2,A2,2X,2I4) 302 FORMAT (1X,5HOUT14,2X,I2,A2,2X,1P,8D12.4) END ******************