************************************************************************ * * * 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 ************************************************************************ * * SUBROUTINE ORTHG (IJK) * * * THIS SUBROUTINE CHECKS ORTHOGONALITY. * * * * IJK CONTROLS THE OUTPUT: * * = 0 FINAL PRINT OUT USING (1), (2), (3), (4), (5); * * = 1 CHECK ORTHOGONALITY OF INITIAL WAVE FUNCTIONS USING (3), * * (4), WITH NO PRINT OUT FROM (3). * * * * SUBROUTINES CALLED : RINTC. * * * ************************************************************************ * {DB: IMPLICIT DOUBLEPRECISION (A-H, O-Z)} CHARACTER*2 NH * DIMENSION IT(2025),JT(2025),RT(2025) * COMMON/CONS/ZERO,HALF,TENTH,ONE,TWO,THREE,TEN : /INFORM/IRD,IPD : /ORB1/E({NW}),GAMA({NW}) COMMON/WAVEC/PFC({NPC},{NW}),QFC({NPC},{NW}) : /GRIDC/RGRIDC({NPC}),ECONT,HCONT,MHALF,NCONT,NUMWR,IPHASE : /DEF1/Z,ATW,RNT,H,EPH,ACCY,C,N : /ORB2/IQ({NC},{NW}) : /ORB4/NW,NCF,NP({NW}),NAK({NW}) : /ORB10/NH({NW}) : /WAVE/PF({NP},{NW}),QF({NP},{NW}) : /FIXD/JFIX({NW}) : /KNTR/ITC(20) * * SECTION (3) * 3 IF (NW .EQ. 1) GO TO 8 MT = 0 NTX = 2025 4 IF (MT .EQ. 5) GO TO 17 IF (IJK .NE. 1) WRITE (IPD,310) NT = 0 NWM = NW-1 DO 5 I = 1,NWM IP = I+1 DO 5 J = IP,NW IF (NAK(I) .NE. NAK(J)) GO TO 5 * IF (JFIX(I)+JFIX(J) .EQ. 2) GO TO 5 ORTH = RINTC (I,J,0) IF (IJK .NE. 1) WRITE (IPD,313) NP(I),NH(I),NP(J),NH(J), : ORTH IF (ABS (ORTH) .LT. 1.0D-06) GO TO 5 NT = NT+1 IF (NT .GT. NTX) GO TO 18 IT(NT) = I JT(NT) = J RT(NT) = ORTH 5 CONTINUE 8 CONTINUE RETURN * 17 WRITE (IPD,314) STOP 7 18 WRITE (IPD,315) STOP 5 * 304 FORMAT (1X) 310 FORMAT (/1X,23HORTHOGONALITY INTEGRALS/) 313 FORMAT (1X,I2,A2,2X,I2,A2,3X,1P,E15.7) 314 FORMAT (/1X,22HERROR IN ORTHOGONALITY) 315 FORMAT (/1X,31HDIMENSION ERROR FOR NTX IN PROP) 316 FORMAT (1X,I2,A2,2X,I2,A2,2X,22HSCHMIDT ORTHOGONALISED) END ************************************************************************ * * FUNCTION RINTC (I,J,K) * * * THE VALUE OF THIS FUNCTION IS THE INTEGRAL OF * * THE FUNCTION: * * * * R**K * ( PI(R)*PJ(R) + QI(R)*QJ(R) ) * * * * FOR EXAMPLE, WHEN I = J AND K = 1, THIS GIVES * * THE MEAN RADIUS FOR ORBITAL I. * * * * SUBROUTINES CALLED: SPLINE,SPLINT (SEE FN RINTI), QUAD. * * * ************************************************************************ * {DB: IMPLICIT DOUBLEPRECISION (A-H, O-Z)} * COMMON/WAVEC/PFC({NPC},{NW}),QFC({NPC},{NW}) : /DEF1/Z,ATW,RNT,H,EPH,ACCY,C,N : /TATBC/TAC({NPC}),TBC({NPC}) : /WAVE/PF({NP},{NW}),QF({NP},{NW}) : /TATB/TA({N1}),TB({N1}) : /KNTR/ITC(20) : /ORB2/IQ({NC},{NW}) : /ORB4/NW,NCF,NP({NW}),NAK({NW}) : /GRIDC/RGRIDC({NPC}),ECONT,HCONT,MHALF,NCONT,NUMWR,IPHASE : /INFORM/IRD,IPD * DIMENSION CC({NPC},3) * * TABULATE INTEGRAND AS REQUIRED FOR SUBROUTINE SPLINT. * DO 1 L = 1,NCONT TAC(L) = (RGRIDC(L)**K)*(PFC(L,I)*PFC(L,J)+ : QFC(L,I)*QFC(L,J)) 1 CONTINUE * * CALL CUBIC SPLINES SUBROUTINES FOR INTEGRATION. * CALL SPLINE (RGRIDC,TAC,NCONT,CC) CALL SPLINT (RGRIDC,TAC,1,NCONT,NCONT,CC,RESULT,H,RNT) RINTC = RESULT RETURN END ************************************************************************ * * SUBROUTINE SIMP (W,H,MIN,MAX,SUM,IFLAG) * * * * * THIS SUBROUTINE INTEGRATES ARRAY W WITH INTERVAL H. IT IS USED * * FOR THE CONTINUUM WAVEFUNCTION INTEGRATIONS. * * * * * * NO SUBROUTINES CALLED. * * * ************************************************************************ * {DB: IMPLICIT DOUBLEPRECISION (A-H, O-Z)} * COMMON/GRIDC/RGRIDC({NPC}),ECONT,HCONT,MHALF,NCONT,NUMWR,IPHASE DIMENSION W({NPC}) IF ((IFLAG .EQ. 0) .AND. (MIN .EQ. 1)) THEN DO 100 II=1,MAX W(II)=W(II)*RGRIDC(II) 100 CONTINUE END IF Z1=(W(MIN)-W(MAX))/2.D0 Z2=0.D0 Z3=0.D0 Z4=0.D0 Z5=0.D0 II=MIN+4 IF(MIN.EQ.1)II=II-1 DO 1 M=II,MAX,4 Z1=Z1+W(M) Z2=Z2+W(M-1)+W(M-3) Z3=Z3+W(M-2) 1 CONTINUE SUM=2.D0*H/45.D0*(14.D0*Z1+32.D0*Z2+12.D0*Z3) IF (MIN .EQ. 1) SUM=SUM+(W(3)-3*W(2)+2*W(1))*H*14.D0/45.D0 Z1=(W(MIN)-W(MAX))/2.D0 Z2=0.D0 Z3=0.D0 Z4=0.D0 Z5=0.D0 IMIN=MIN IF (MIN .EQ. 1) IMIN=0 NN=(MAX+IMIN)/2 DO 2 M=II,NN,4 Z1=Z1+W(2*M-IMIN) Z2=Z2+W(2*M-IMIN-2)+W(2*M-IMIN-6) Z3=Z3+W(2*M-IMIN-4) 2 CONTINUE SUM1=4.D0*H*(14.D0*Z1+32.D0*Z2+12.D0*Z3)/45.D0 IF (MIN.EQ.1) SUM1=SUM1+(W(6)-3*W(4)+2*W(2))*2.D0*H*14.D0/45.D0 RETURN END ************************************************************************ * * SUBROUTINE NORMC (IORBIT,C00TEMP,PHASE) * * * THIS SUBROUTINE ENERGY NORMALIZES THE CONTINUUM WAVEFUNCTION. * * FOR CASES WHERE THE CORE IS A NEUTRAL USING EQUATION (15) OF BATES,* * CPC, V.8, PG. 220-235, 1974. A 2-PARAMETER FIT IS PERFORMED, ONE * * FOR THE AMPLITUDE AND THE OTHER FOR THE PHASE. * * * * INPUTS: * * J : SERIAL NUMBER OF ORBITAL * * P,Q : LARGE, SMALL COMPONENTS OF WAVEFUNCTION * * * * OUTPUTS: * * PF,QF,P,Q : AFTER RENORMALIZATION * * * * SUBROUTINES CALLED: SPHJ, SPHN. * * * ************************************************************************ * {DB: IMPLICIT DOUBLEPRECISION (A-H, O-Z)} IMPLICIT INTEGER (I-N) CHARACTER*2 NH * COMMON/CONS/ZERO,HALF,TENTH,ONE,TWO,THREE,TEN COMMON/BESS3/BJ3({NPC},3) : /DEF1/Z,ATW,RNT,H,EPH,ACCY,C,N : /INFORM/IRD,IPD COMMON/WAVEC/PFC({NPC},{NW}),QFC({NPC},{NW}) : /POTEC/YPC({NPC}),YQC({NPC}),XTPC({NPC}),XTQC({NPC}) : /ORB1/E({NW}),GAMADD({NW}) : /ORB4/NW,NCF,NP({NW}),NAK({NW}) : /ORB10/NH({NW}) : /ORB5/NKL({NW}),NKJ({NW}),KAPMAX : /POTE/YP({NP}),YQ({NP}),XP({NP}),XQ({NP}) COMMON/GRIDC/RGRIDC({NPC}),ECONT,HCONT,MHALF,NCONT,NUMWR,IPHASE * DIMENSION CMAT(2),F(200,200),A(200,200),B(200) DIMENSION FTEMP (200),F9(200),P({NPC}),Q({NPC}) * DO 60 I=1,NCONT P(I)=PFC(I,IORBIT) Q(I)=QFC(I,IORBIT) 60 CONTINUE NUMAX = 40 NUMIT = 0 NUMITX = 0 NNORM = 40 ERR = 1.0D-2 IFINAL = 0 M = 2 PI = 4.0D0*ATAN (1.0D0) SIX = TWO*THREE * WRITE (IPD,*) ' ENTERING RENORM SUBROUTINE' * * INTERPOLATE THE DIRECT POTENTIAL * KAPPA = NAK(IORBIT) XMOMEN = SQRT((-E(IORBIT)+C*C)*(-E(IORBIT)+C*C)-C*C*C*C)/C * XMOMEN = SQRT(-TWO*E(IORBIT)) ENERGY = -E(IORBIT)+C*C MSTEP = 10 C00 = ONE C11 = 0.01 * * CALL SPHJ (XMOMEN,NKL(IORBIT)) CALL SPHN (XMOMEN,NKL(IORBIT)) 50 KOUNT = 1 DO 195 J = 1,NNORM JJ = NCONT-(MSTEP+1)*NNORM+1+(J-1)*MSTEP RPT = RGRIDC(JJ) F9(J) = C00*SQRT(TWO*XMOMEN/PI)*RPT*(COS(C11)*BJ3(JJ,1) : -SIN(C11)*BJ3(JJ,2)) F(1,KOUNT) = F9(J)/C00 F(2,KOUNT) = C00*SQRT(TWO*XMOMEN/PI)*RPT*(-SIN(C11)*BJ3(JJ,1) : -COS(C11)*BJ3(JJ,2)) FTEMP(KOUNT) = F9(J) KOUNT = KOUNT + 1 195 CONTINUE DO 5 I = 1,M DO 3 K = 1,I A(K,I) = 0.0E0 DO 2 J = 1,NNORM A(K,I) = A(K,I) + F(I,J)*F(K,J) 2 CONTINUE A(I,K) = A(K,I) 3 CONTINUE 5 CONTINUE DO 6 K = 1,M B(K) = 0.0E0 DO 4 J = 1,NNORM JX = NCONT-(MSTEP+1)*NNORM+1+(J-1)*MSTEP B(K) = B(K)+(P(JX)-FTEMP(J))*F(K,J) 4 CONTINUE 6 CONTINUE DET = A(1,1)*A(2,2)-A(2,1)*A(1,2) CMAT(1) = (B(1)*A(2,2)-B(2)*A(1,2))/DET CMAT(2) = (B(2)*A(1,1)-B(1)*A(2,1))/DET * WRITE (IPD,*) ' C00,C11=',C00,C11 IF (((ABS(CMAT(1)/C00)).LT.1.0E-03).AND. : ((ABS(CMAT(2))).LT.1.0E-03)) THEN C00 = C00 + CMAT(1) C11 = C11 + CMAT(2) ELSE C00 = C00 + (CMAT(1)/2.0) C11 = C11 + (CMAT(2)/2.0) ENDIF * SOME DAMPING IS PROVIDED TO AVOID RAPID OSCILLATIONS IN THE VALUES * OF C00 AND C11 C11 = MOD(C11,TWO*PI) PHASE = C11 * WRITE (IPD,*) ' C00,C11=',C00,C11 IF (ABS(CMAT(1)/C00)+ABS(CMAT(2)) .GT. ERR) THEN NUMIT = NUMIT + 1 IF (NUMIT .EQ. NUMAX) THEN WRITE (IPD,*) ' NUMBER OF ITERATIONS EXCEEDED IN SR NORMC' STOP END IF GO TO 50 ELSE IF (IFINAL .EQ. 0) THEN IFINAL = 1 ERR = 1.0E-4 NUMITX = NUMIT+1 NUMIT = 0 NNORM = 40 GO TO 50 END IF END IF C00TEMP = ONE/ABS(C00) DO 55 I = 1,NCONT P(I) = P(I)*C00TEMP PFC(I,IORBIT) = P(I) Q(I) = Q(I)*C00TEMP QFC(I,IORBIT) = Q(I) 55 CONTINUE * WRITE (IPD,*) ' NUMBER OF ITERATIONS NEEDED IN NORMC=',NUMIT+NUMITX * WRITE (IPD,*) ' ERROR TERMS=',CMAT(1),CMAT(2) rewind 1 do 1301 i=1,ncont,10 write (1,305) RGRIDC(i),pfc(i,iorbit) write (1,306) RGRIDC(i), : SQRT(TWO*XMOMEN/PI)*RGRIDC(i)*(COS(C11)*BJ3(I,1) : -SIN(C11)*BJ3(I,2)) 1301 continue rewind 1 305 format ('1',5P,1e11.4,1e11.4) 306 format ('2',5P,1e11.4,1e11.4) RETURN END *********************************************************************** * * SUBROUTINE SPHJ (XK,NN) * * * ---------------- SECTION 10 SUBPROGRAM ?? ---------------- * * * * THIS ROUTINE EVALUATES THE SPHERICAL BESSEL FUNCTION, jl(xk*r), * * GIVEN THE ELECTRON'S WAVENUMBER, XK, (RELATED TO THE MOMENTUM) * * FOR ORBITAL QUANTUM NUMBER, NN (THE L QUANTUM NUMBER) AND * * STORES THE RESULT IN BJ(I,1), WHERE I=1...N. RECURSION FORMULA * * ARE USED STARTING FROM EXACT EXPRESSIONS FOR j0(xk*r) AND j1(xk*r).* * * * NO SUBROUTINES CALLED. * * * * WRITTEN BY WARREN F PERGER, MICHIGAN TECH UNIVERSITY * * LAST REVISION: 19 DEC 1990 * * * ************************************************************************ * {DB: IMPLICIT DOUBLEPRECISION (A-H, O-Z)} * COMMON/BESS3/BJ3({NPC},3) : /CONS/ZERO,HALF,TENTH,ONE,TWO,THREE,TEN : /GRIDC/RGRIDC({NPC}),ECONT,HCONT,MHALF,NCONT,NUMWR,IPHASE : /INFORM/IRD,IPD COMMON/OPT4/LTC(30) DIMENSION RTEMP({NPC}) DO 10 I=1,NCONT RTEMP(I)=RGRIDC(I)*XK 10 CONTINUE * EPSI = 1.0D-05 * DO 50 I=1,NCONT IF (RTEMP(I) .EQ. ZERO) THEN BJ3(I,1)=ONE ELSE BJ3(I,1)=SIN(RTEMP(I))/RTEMP(I) END IF 50 CONTINUE DO 51 I=1,NCONT IF (RTEMP(I) .EQ. ZERO) THEN BJ3(I,2)=ZERO ELSE BJ3(I,2)=SIN(RTEMP(I))/(RTEMP(I)*RTEMP(I))-COS(RTEMP(I)) : /RTEMP(I) END IF 51 CONTINUE IF (NN .EQ. 0) RETURN IF (NN .EQ. 1) THEN DO 52 I=1,NCONT BJ3(I,1)=BJ3(I,2) 52 CONTINUE RETURN END IF DO 60 I=2,NN DO 65 I2=1,NCONT BJ3(I2,3)=(TWO*DBLE(I)-ONE)*BJ3(I2,2)/RTEMP(I2) - BJ3(I2,1) 65 CONTINUE DO 70 I2=1,NCONT BJ3(I2,1)=BJ3(I2,2) BJ3(I2,2)=BJ3(I2,3) 70 CONTINUE 60 CONTINUE DO 75 I=1,NCONT BJ3(I,1)=BJ3(I,3) 75 CONTINUE * * PRINT OUT BESSEL FUNCTIONS IF (DEBUG) OPTION SET * IF (LTC(16) .EQ. 1) THEN WRITE (IPD,300) NN,(BJ3(II,1),II = 1,NCONT) 23 CONTINUE ENDIF * * ALL DONE * RETURN * 300 FORMAT (/' BESSEL FUNCTION OF ORDER ',I3,/(1P,7D18.10)) * END ************************************************************************ * * SUBROUTINE SPHN (XK,NN) * * * ---------------- SECTION 10 SUBPROGRAM ?? ---------------- * * * * THIS ROUTINE EVALUATES THE SPHERICAL BESSEL FUNCTION, nl(xk*r), * * GIVEN THE ELECTRON'S WAVENUMBER, XK, (RELATED TO THE MOMENTUM) * * FOR ORBITAL QUANTUM NUMBER, NN (THE L QUANTUM NUMBER) AND * * STORES THE RESULT IN BJ(I,2), WHERE I=1...N. RECURSION FORMULA * * ARE USED STARTING FROM EXACT EXPRESSIONS FOR n0(xk*r) AND n1(xk*r).* * * * NO SUBROUTINES CALLED. * * * * WRITTEN BY WARREN F PERGER, MICHIGAN TECH UNIVERSITY * * LAST REVISION: 19 DEC 1990 * * * ************************************************************************ * {DB: IMPLICIT DOUBLEPRECISION (A-H, O-Z)} * COMMON/BESS3/BJ3({NPC},3) : /CONS/ZERO,HALF,TENTH,ONE,TWO,THREE,TEN : /GRIDC/RGRIDC({NPC}),ECONT,HCONT,MHALF,NCONT,NUMWR,IPHASE : /INFORM/IRD,IPD COMMON/OPT4/LTC(30) DIMENSION RTEMP({NPC}),BJ4({NPC},3) DO 10 I=1,NCONT RTEMP(I)=RGRIDC(I)*XK 10 CONTINUE * EPSI = 1.0D-05 * DO 50 I=1,NCONT IF (RTEMP(I) .EQ. ZERO) THEN BJ4(I,1)=ONE ELSE BJ4(I,1)=-COS(RTEMP(I))/RTEMP(I) END IF 50 CONTINUE DO 51 I=1,NCONT IF (RTEMP(I) .EQ. ZERO) THEN BJ4(I,2)=ONE ELSE BJ4(I,2)=-COS(RTEMP(I))/(RTEMP(I)*RTEMP(I))-SIN(RTEMP(I)) : /RTEMP(I) END IF 51 CONTINUE IF (NN .EQ. 0) THEN DO 52 I=1,NCONT BJ3(I,2)=BJ4(I,1) 52 CONTINUE GO TO 11 END IF IF (NN .EQ. 1) THEN DO 53 I=1,NCONT BJ3(I,2)=BJ4(I,2) 53 CONTINUE GO TO 11 END IF DO 60 I=2,NN DO 65 I2=1,NCONT BJ4(I2,3)=(TWO*DBLE(NN)-ONE)*BJ4(I2,2)/RTEMP(I2) - BJ4(I2,1) 65 CONTINUE DO 70 I2=1,NCONT BJ4(I2,1)=BJ4(I2,2) BJ4(I2,2)=BJ4(I2,3) 70 CONTINUE 60 CONTINUE DO 75 I=1,NCONT BJ3(I,2)=BJ4(I,3) 75 CONTINUE * * PRINT OUT BESSEL FUNCTIONS IF (DEBUG) OPTION SET * 11 IF (LTC(16) .EQ. 1) THEN WRITE (IPD,300) NN,(BJ3(II,2),II = 1,NCONT) 23 CONTINUE ENDIF * * ALL DONE * RETURN * 300 FORMAT (/' BESSEL FUNCTION OF ORDER ',I3,/(1P,7D18.10)) * END ************************************************************************ * * SUBROUTINE SLATCO (I3,I1,I4,I2,K,RESULT) * * * * * THIS SUBROUTINE PERFORMS THE SLATER INTEGRAL FOR THE CONTINUUM * * WAVEFUNCTIONS. The inner integral indices are I1 and I2; the * * outer integral indices are I3 and I4. * * * * * * SUBROUTINES CALLED: SIM. * * * ************************************************************************ * IMPLICIT DOUBLE PRECISION (A-H,O-Z) IMPLICIT INTEGER (I-N) * COMMON/WAVEC/PFC({NPC},{NW}),QFC({NPC},{NW}) : /DEF1/Z,ATW,RNT,H,EPH,ACCY,C,N : /GRIDC/RGRIDC({NPC}),ECONT,HCONT,MHALF,NCONT,NUMWR,IPHASE : /INFORM/IRD,IPD : /KNTR/ITC(20) : /ORB4/NW,NCF,NP({NW}),NAK({NW}) * DIMENSION W1({NPC}),W2({NPC}),W3({NPC}),W4({NPC}),GRAL({NPC}), : ORD({NPC}),PI({NPC}),QI({NPC}) * XKA = DBLE(K+1) C1 = SQRT (NAK(I1)**2 - (Z/C)**2) C2 = SQRT (NAK(I2)**2 - (Z/C)**2) MM = MHALF DO 622 M = 1,NCONT PI(M) = PFC(M,I1) QI(M) = QFC(M,I1) GRAL(M) = PFC(M,I2) ORD(M) = QFC(M,I2) 622 CONTINUE DO 40 M=1,NCONT W1(M)=GRAL(M)*PI(M)+ORD(M)*QI(M) 40 CONTINUE DO 50 M=1,MM ORD(M)=RGRIDC(M)**XKA*W1(M) 50 CONTINUE GRAL(1)=ORD(1)/(C1+C2+XKA) GRAL(2)=ORD(2)/(C1+C2+XKA) DO 60 M=3,MM GRAL(M)=GRAL(M-2)+H*(ORD(M-2)+4.D0*ORD(M-1)+ORD(M))/3.D0 60 CONTINUE DO 70 M=1,MM W1(M)=GRAL(M)/RGRIDC(M)**(XKA-1.0) 70 CONTINUE A1=0.D0 A2=0.D0 GRAL(NCONT)=0.D0 M1=MM+1 DO 80 M=M1,NCONT ORD(M)=RGRIDC(M)**(XKA-1.D0)*W1(M) 80 CONTINUE HW=GRAL(MM) OH=ORD(MM)/RGRIDC(MM) GRAL(MM+1)=GRAL(MM)+HCONT*(5.D0*OH+8.D0*ORD(MM+1)-ORD(MM+2))/ : 12.00 M1=MM+2 HW=ORD(MM) ORD(MM)=OH DO 90 M=M1,NCONT GRAL(M)=GRAL(M-2)+HCONT*(ORD(M-2)+4.D0*ORD(M-1)+ORD(M))/3.D0 90 CONTINUE M1=MM+1 DO 100 M=M1,NCONT W1(M)=GRAL(M)/RGRIDC(M)**(XKA+1.D0-1.0) 100 CONTINUE DO 110 M=MM,NCONT ORD(M)=ORD(M)/RGRIDC(M)**(2.D0*XKA-1.0) 110 CONTINUE GRAL(NCONT)=0.D0 GRAL(NCONT-1)=HCONT*(5.D0*ORD(NCONT)+8.D0*ORD(NCONT-1)-ORD(NCONT-2 : ))/12.0D0 M=NCONT-1 120 M=M-1 GRAL(M)=GRAL(M+2)+HCONT*(ORD(M+2)+4.D0*ORD(M+1)+ORD(M))/3.D0 IF (M .GT. MM) GO TO 120 ORD(MM)=HW DO 140 M=1,MM ORD(M)=ORD(M)/RGRIDC(M)**(2.D0*XKA-1.0) 140 CONTINUE OH=ORD(MM) GRAL(MM-1)=GRAL(MM)+H*(5.D0*OH+8.D0*ORD(MM-1)-ORD(MM-2))/12.0 M=MM-1 150 M=M-1 GRAL(M)=GRAL(M+2)+H*(ORD(M+2)+4.D0*ORD(M+1)+ORD(M))/3.D0 IF (M .GT. 1) GO TO 150 DO 160 M=1,MM W4(M)=GRAL(M)*RGRIDC(M)**XKA 160 CONTINUE M1=MM+1 DO 170 M=M1,NCONT W4(M)= RGRIDC(M)**(XKA-1.D0)*GRAL(M) 170 CONTINUE DO 623 M = 1,NCONT PI(M) = PFC(M,I3) QI(M) = QFC(M,I3) 623 CONTINUE DO 624 M = 1,NCONT GRAL(M) = PFC(M,I4) ORD(M) = QFC(M,I4) 624 CONTINUE DO 192 M=1,NCONT W2(M)=GRAL(M)*PI(M)+ORD(M)*QI(M) 192 CONTINUE DO 202 M=1,NCONT W3(M)=(W1(M)+W4(M))*W2(M) 202 CONTINUE A1=0.D0 A2=0.D0 HW=W3(MM) W3(MM)=HW/RGRIDC(MM) CALL SIM(W3,HCONT,MM,NCONT,SUM1,SUM) A1=SUM1 A2=SUM W3(MM)=HW CALL SIM(W3,H,1,MM,SUM1,SUM) A1=A1+SUM1 A2=A2+SUM ST=A2-A1 KL=K * IF (ITC(25) .EQ. 1) WRITE(IPD,270) I1,I2,I3,I4,KL,SUM,A2,ST * 270 FORMAT(1X,5I5,3E15.6) RESULT = A2 RETURN END ************************************************************************ * * SUBROUTINE SIM (W,H,MIN,MAX,SUM1,SUM) * * * * * THIS SUBROUTINE INTEGRATES ARRAY W WITH INTERVAL H. IT IS USED * * FOR THE CONTINUUM WAVEFUNCTION INTEGRATIONS. * * * * * * NO SUBROUTINES CALLED. * * * ************************************************************************ * IMPLICIT DOUBLE PRECISION (A-H,O-Z) IMPLICIT INTEGER (I-N) * DIMENSION W({NPC}) Z1=(W(MIN)-W(MAX))/2.D0 Z2=0.D0 Z3=0.D0 Z4=0.D0 Z5=0.D0 II=MIN+4 IF(MIN.EQ.1)II=II-1 DO 1 M=II,MAX,4 Z1=Z1+W(M) Z2=Z2+W(M-1)+W(M-3) Z3=Z3+W(M-2) 1 CONTINUE SUM=2.D0*H/45.D0*(14.D0*Z1+32.D0*Z2+12.D0*Z3) IF (MIN .EQ. 1) SUM=SUM+(W(3)-3*W(2)+2*W(1))*H*14.D0/45.D0 Z1=(W(MIN)-W(MAX))/2.D0 Z2=0.D0 Z3=0.D0 Z4=0.D0 Z5=0.D0 IMIN=MIN IF (MIN .EQ. 1) IMIN=0 NN=(MAX+IMIN)/2 DO 2 M=II,NN,4 Z1=Z1+W(2*M-IMIN) Z2=Z2+W(2*M-IMIN-2)+W(2*M-IMIN-6) Z3=Z3+W(2*M-IMIN-4) 2 CONTINUE SUM1=4.D0*H*(14.D0*Z1+32.D0*Z2+12.D0*Z3)/45.D0 IF (MIN.EQ.1) SUM1=SUM1+(W(6)-3*W(4)+2*W(2))*2.D0*H*14.D0/45.D0 RETURN END ************************************************************************ ************************************************************************ SUBROUTINE HYDROG (IORBIT) {DB: IMPLICIT DOUBLE PRECISION (A-H,O-Z)} IMPLICIT INTEGER (I-N) DIMENSION F({NPC}),G({NPC}),PP2({NPC}),PQ2({NPC}) COMPLEX*16 GAMMP,FACTR,EIETA,RIP,F1,F2,F3,TX COMPLEX*16 ZTEMP,A,B,Z2,CONHYP,CMPGAMM COMMON/DEF1/Z,ATW,RNT,H,EPH,ACCY,C,N : /CONS/ZERO,HALF,TENTH,ONE,TWO,THREE,TEN : /INFORM/IRD,IPD : /ORB1/E({NW}),GAM({NW}) : /ORB4/NW,NCF,NP({NW}),NAK({NW}) COMMON/WAVEC/PFC({NPC},{NW}),QFC({NPC},{NW}) : /POTEC/YPC({NPC}),YQC({NPC}),XTPC({NPC}),XTQC({NPC}) : /GRIDC/RGRID({NPC}),ECONT,HCONT,MHALF,NCONT,NUMWR,IPHASE * DO 199 I = 1,NCONT PP2(I) = PFC(I,IORBIT) PQ2(I) = QFC(I,IORBIT) 199 CONTINUE 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 WRITE(IPD,*) 'IORBIT= ',IORBIT WRITE(IPD,*) 'E(IORBIT)= ',E(IORBIT) ENERGY = -E(IORBIT)+C*C WRITE(IPD,*) 'ENERGY+M*C**2 =',ENERGY * XMOMEN = SQRT(((ENERGY+C*C)**TWO-C*C*C*C)/(C*C)) WRITE (IPD,*) ' XMOMEN=',XMOMEN ZEFF = YPC(NCONT-20) WRITE (IPD,*) ' NAK(IORBIT),ZEFF,C',NAK(IORBIT),ZEFF,C GAMA = SQRT (NAK(IORBIT)*NAK(IORBIT) - (ZEFF/C)*(ZEFF/C)) WRITE(IPD,*) 'GAMA= ',GAMA Y = ZEFF*ENERGY/(XMOMEN*C*C) WRITE(IPD,*) 'Y= ',Y GAMMP = DCMPLX(GAMA,Y) WRITE(IPD,*) 'GAMMP= ',GAMMP FACTR = DCMPLX(DBLE(NAK(IORBIT)),-ZEFF/XMOMEN) WRITE(IPD,*) 'FACTR= ',FACTR EIETA = SQRT (-FACTR/GAMMP) WRITE(IPD,*) 'EIETA= ',EIETA * XTEMP1 = LGAMMA (TWO*GAMA+ONE) * C1 = EXP(XTEMP1) XTEMP1 = DCMPLX(TWO*GAMA+ONE,ZERO) C1 = CMPGAMM(XTEMP1,0) TX = CMPGAMM(GAMMP,1) C2 = DBLE (TX) PI = 4.0D0*ATAN(1.0D0) C3 = EXP (PI*Y/TWO+C2)/(TWO*SQRT(PI*XMOMEN)*C1) ZTEMP = DCMPLX (TWO*GAMA+ONE,ZERO) WRITE (IPD,*) 'C1= ',C1,' C2=',C2,' C3=',C3 NSTEP = 20 DO 1 JJ = 1,NCONT,NSTEP RIP = DCMPLX(ZERO,XMOMEN*RGRID(JJ)) * IF (RGRID(JJ)*XMOMEN .LE. 9000.0) THEN FACTR = GAMMP*EIETA*CDEXP (-RIP) FCONS = SQRT (ENERGY/(C*C)-ONE)*(TWO*XMOMEN : *RGRID(JJ))**GAMA GCONS = SQRT (ENERGY/(C*C)+ONE)*(TWO*XMOMEN : *RGRID(JJ))**GAMA F3 = DCMPLX(GAMA+ONE,Y) A=F3 B=ZTEMP Z2=TWO*RIP F1 = FACTR*CONHYP(A,B,Z2,0,0) F2 = CONJG (F1) F(JJ) = FCONS*DIMAG ((F1-F2))*C3 G(JJ) = GCONS*DBLE ((F1+F2))*C3 * CALL D2OUT(3,1,RGRID(JJ),G(JJ)) * CALL D2OUT(3,2,RGRID(JJ),PP2(JJ)) * ELSE * DELTA = Y*LOG(2*XMOMEN*RGRID(JJ)) - ATAN(DIMAG(TX)/ * : DBLE(TX)) + DIMAG(LOG(EIETA)) - 0.5*3.14159*GAMA * F(JJ) = -SQRT((ENERGY/(C*C)-1.0D0)/(3.14159*XMOMEN/C)) * : *SIN(XMOMEN*RGRID(JJ)+DELTA)*RGRID(JJ) * G(JJ) = SQRT((ENERGY/(C*C)+1.0D0)/(3.14159*XMOMEN/C)) * : *COS(XMOMEN*RGRID(JJ)+DELTA)*RGRID(JJ) * END IF 1 CONTINUE WRITE(IPD,*) ' R*G(R)= ' DO 7 M = 1,NCONT,NSTEP * WRITE (IPD,1001) RGRID(M),(G(I),I=M,M+9) 7 CONTINUE WRITE(IPD,*) ' R*F(R)= ' DO 9 M = 1,NCONT,NSTEP * WRITE(IPD,1001) RGRID(M),(F(I),I=M,M+9) 9 CONTINUE 1000 FORMAT(1P,10(E12.4)) 1001 FORMAT(1P,11(E16.7)) WRITE(IPD,*) 'RATIO OF PERGER/ROSE= ' WRITE(IPD,*) 'PP2/G = ' DO 17 M = 1,NCONT,NSTEP WRITE(IPD,1001) RGRID(M),PP2(M)/G(M) * WRITE(IPD,1001) RGRID(M),(PP2(I)/G(I),I=M,M+90,10) 17 CONTINUE WRITE(IPD,*) 'PQ2/F = ' DO 18 M = 1,NCONT,NSTEP WRITE(IPD,1001) RGRID(M),PQ2(M)/F(M) * WRITE(IPD,1001) RGRID(M),(PQ2(I)/F(I),I=M,M+90,10) 18 CONTINUE * WRITE (IPD,1030) 30.0/(XMOMEN*2.0) *1030 FORMAT (/, 'SWITCHING OVER TO ASYMPTOTIC EXPANSION AT R= ', * : 1P,1E12.6) RETURN END ************************************************************************ * * FUNCTION HYD (JJ,IORBIT) * * * THIS FUNCTION RETURNS AN APPROXIMATION FOR THE HYDROGENIC * * CONTINUUM WAVEFUNCTION. * * * ************************************************************************ * {DB: IMPLICIT DOUBLE PRECISION (A-H,O-Z)} IMPLICIT INTEGER (I-N) DIMENSION F({NPC}),G({NPC}) COMPLEX*16 HYPER2,GAMMP,FACTR,EIETA,RIP,F1,F2,F3,TX COMPLEX*16 ZTEMP,Z2,CONHYP,CMPGAMM,XTEMP1 COMMON/DEF1/Z,ATW,RNT,H,EPH,ACCY,C,N : /CONS/ZERO,HALF,TENTH,ONE,TWO,THREE,TEN : /INFORM/IRD,IPD : /ORB1/E({NW}),GAM({NW}) : /ORB4/NW,NCF,NP({NW}),NAK({NW}) COMMON/WAVEC/PFC({NPC},{NW}),QFC({NPC},{NW}) : /POTEC/YPC({NPC}),YQC({NPC}),XTPC({NPC}),XTQC({NPC}) : /GRIDC/RGRID({NPC}),ENERGY,HCONT,MHALF,NCONT,NUMWR,IPHASE * 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 * WRITE(IPD,*) 'IORBIT= ',IORBIT * WRITE(IPD,*) 'E(IORBIT)= ',E(IORBIT) ENERGY = -E(IORBIT)+C*C * WRITE(IPD,*) 'ENERGY =',ENERGY * WRITE (IPD,*) ' XMOMEN=',XMOMEN ZEFF = YPC(JJ) * WRITE (IPD,*) ' NAK(IORBIT),ZEFF,C',NAK(IORBIT),ZEFF,C GAMA = SQRT (NAK(IORBIT)*NAK(IORBIT) - (ZEFF/C)*(ZEFF/C)) * WRITE(IPD,*) 'GAMA= ',GAMA Y = ZEFF*ENERGY/(XMOMEN*C*C) * WRITE(IPD,*) 'Y= ',Y GAMMP = DCMPLX(GAMA,Y) * WRITE(IPD,*) 'GAMMP= ',GAMMP FACTR = DCMPLX(DBLE(NAK(IORBIT)),-ZEFF/XMOMEN) * WRITE(IPD,*) 'FACTR= ',FACTR EIETA = SQRT (-FACTR/GAMMP) * WRITE(IPD,*) 'EIETA= ',EIETA * XTEMP1 = LGAMMA (TWO*GAMA+ONE) * C1 = EXP(XTEMP1) XTEMP1 = DCMPLX(TWO*GAMA+ONE,ZERO) C1 = CMPGAMM(XTEMP1,0) TX = CMPGAMM(GAMMP,1) C2 = DBLE (TX) PI = 4.0D0*DATAN(1.0D0) C3 = EXP (PI*Y/TWO+C2)/(TWO*SQRT(PI*XMOMEN)*C1) ZTEMP = DCMPLX (TWO*GAMA+ONE,ZERO) * WRITE (IPD,*) 'C1= ',C1,' C2=',C2,' C3=',C3 RIP = DCMPLX(ZERO,XMOMEN*RGRID(JJ)) FACTR = GAMMP*EIETA*CDEXP (-RIP) FCONS = SQRT (ENERGY/(C*C)-ONE)*(TWO*XMOMEN*RGRID(JJ))**GAMA GCONS = SQRT (ENERGY/(C*C)+ONE)*(TWO*XMOMEN*RGRID(JJ))**GAMA F3 = DCMPLX(GAMA+ONE,Y) Z2 = TWO*RIP IF (ABS(Z2) .GT. 100.0D0) THEN F1 = FACTR*HYPER2 (F3,ZTEMP,Z2) ELSE F1 = FACTR*CONHYP(F3,ZTEMP,Z2,0,0) END IF F2 = CONJG (F1) F(JJ) = FCONS*DIMAG ((F1-F2))*C3 G(JJ) = GCONS*DBLE ((F1+F2))*C3 HYD = G(JJ) RETURN END ************************************************************************ * * FUNCTION HYPER2 (ALPHA,GAMAX,Z) * * * THIS FUNCTION RETURNS AN APPROXIMATION FOR 1F1 FOR LARGE R. * * * ************************************************************************ * {DB: IMPLICIT DOUBLE PRECISION (A-H,O-Z)} COMPLEX*16 TERM,TERM1,HYPER2,SUM1,SUM2,CMPGAMM COMPLEX*16 ALPHA,Z,GAMAX,XTEMP1,XTEMP2,IM DOUBLE PRECISION IX COMMON/INFORM/IRD,IPD * IM = (0.0D0,1.0D0) PI = 4.0D0*ATAN(1.0D0) GAMA = DBLE(GAMAX) SUM1 = DCMPLX (0.0D0,0.0D0) SUM2 = DCMPLX (0.0D0,0.0D0) TERM1 = DCMPLX(1.0D0,0.0D0) DO 1 I = 1,100 IX = FLOAT (I) TERM = (ALPHA+IX-1.0D0)*(ALPHA-GAMA+IX)/(IX*(-Z**IX)) TERM1 = TERM*TERM1 SUM1 = TERM1 + SUM1 IF ((ABS (TERM1)/ABS (SUM1) .LT. 1.0D-7) .AND. : (I .GE. 20)) GO TO 2 1 CONTINUE WRITE(IPD,*) 'ERROR IN FN HYPER2--DID NOT CONVERGE IN 100 TERMS' stop 2 SUM1 = SUM1 + 1.0D0 XTEMP1 = DCMPLX (GAMA,ZERO) XTEMP2 = GAMA-ALPHA SUM1 = SUM1*CDEXP(IM*PI*ALPHA)*CMPGAMM(XTEMP1,0)/ : (Z**ALPHA*CMPGAMM(XTEMP2,0)) TERM1 = DCMPLX(1.0D0,0.0D0) DO 3 I = 1,100 IX = FLOAT (I) TERM = (GAMA-ALPHA+IX-1.0D0)*(-ALPHA+IX)/ : (IX*(Z**IX)) TERM1 = TERM*TERM1 SUM2 = TERM1 + SUM2 IF ((ABS (TERM1)/ABS (SUM2) .LT. 1.0D-7) .AND. : (I .GE. 20)) GO TO 4 3 CONTINUE WRITE(IPD,*) 'ERROR IN FN HYPER2--DID NOT CONVERGE IN 100 TERMS' stop 4 SUM2 = SUM2 + 1.0D0 XTEMP2 = ALPHA SUM2 = SUM2*CMPGAMM(XTEMP1,0)*CDEXP(Z)*Z**(ALPHA-GAMA)/ : CMPGAMM(XTEMP2,0) HYPER2 = SUM1 + SUM2 RETURN END ************************************************************************ * * FUNCTION CMPGAMM (ARG,LNGAM) * * * THIS SUBROUTINE RETURNS IN RES THE COMPLEX GAMMA FUNCTION OF THE * * COMPLEX ARGUMENT ARG. IF LNGAM = 1, THEN THE NATURAL LOG OF THE * * COMPLEX GAMMA FUNCTION IS RETURNED; IF LNGAM = 0, THEN THE COMPLEX * * GAMMA FUNCTION IS RETURNED. THE LANCZOS APPROXIMATION IS USED. * * * ************************************************************************ * {DB: IMPLICIT DOUBLE PRECISION (A-H,O-Z)} COMPLEX*16 CMPGAMM,ARG,GAMMLN,X,TMP,SER * DOUBLE PRECISION COF(6),STP,FPF COMMON/INFORM/IRD,IPD * * DATA ZERO,HALF,TENTH,ONE,TWO,THREE,TEN/ : 0.0D 00,0.5D 00,0.1D 00, : 1.0D 00,2.0D 00,3.0D 00, : 10.0D 00/ DATA COF,STP/76.18009173D0,-86.50532033D0,24.01409822D0, : -1.231739516D0,.120858003D-2,-.536382D-5, : 2.50662827465D0/ DATA FPF/5.5D0/ DATA FIRST/.TRUE./ LOGICAL FIRST IF (FIRST) THEN PI= 4.0D0*ATAN(ONE) * * SET THE MACHINE-DEPENDENT PARAMETERS: * * TENMAX - MAXIMUM SIZE OF EXPONENT OF 10 * ITNMAX = 1 DNUM = TENTH 11 ITNMAX = ITNMAX+1 DNUM = DNUM*TENTH IF (DNUM .GT. ZERO) GOTO 11 ITNMAX = ITNMAX-1 TENMAX = DBLE (ITNMAX) * * EXPMAX - MAXIMUM SIZE OF EXPONENT OF E * DNUM = TENTH**ITNMAX EXPMAX = -LOG (DNUM) * * PRECIS - MACHINE PRECISION * PRECIS = ONE 22 PRECIS = PRECIS/TWO DNUM = PRECIS+ONE IF (DNUM .GT. ONE) GOTO 22 PRECIS = TWO*PRECIS * HLNTPI = HALF*LOG (TWO*PI) * FIRST = .FALSE. * ENDIF X=ARG-ONE TMP=X+FPF TMP=(X+HALF)*CDLOG(TMP)-TMP SER=DCMPLX(ONE,ZERO) DO 2 J=1,6 X=X+ONE SER=SER+COF(J)/X 2 CONTINUE GAMMLN=TMP+CDLOG(STP*SER) CLNGR=DBLE(GAMMLN) CLNGI=DIMAG(GAMMLN) * * RETURN IF LOG OF THE GAMMA FUNCTION IS REQUESTED. * IF (LNGAM .EQ. 1) THEN CMPGAMM = DCMPLX(CLNGR,CLNGI) RETURN END IF * * NOW EXPONENTIATE THE COMPLEX LOG GAMMA FUNCTION TO GET * THE COMPLEX GAMMA FUNCTION * IF ( (CLNGR .LE. EXPMAX) .AND. : (CLNGR .GE. -EXPMAX) ) THEN FAC = EXP (CLNGR) ELSE WRITE (IPD,300) WRITE (IPD,301) CLNGR STOP '10' ENDIF RESR = FAC*COS (CLNGI) RESI = FAC*SIN (CLNGI) CMPGAMM = DCMPLX(RESR,RESI) RETURN 300 FORMAT (///' ***** ERROR IN SUBROUTINE CMPGAMM *****') 301 FORMAT (' ARGUMENT TO EXPONENTIAL FUNCTION (',1P,1D14.7, : ') OUT OF RANGE.') END C **************************************************************** C * * C * SOLUTION TO THE CONFLUENT HYPERGEOMETRIC FUNCTION * C * * C * by * C * * C * MARK NARDIN, * C * * C * W. F. PERGER and ATUL BHALLA * C * * C * * C * Michigan Technological University, Copyright 1989 * C * * C * * C * Description : A numerical evaluator for the confluent * C * hypergeometric function for complex arguments with large * C * magnitudes using a direct summation of the Kummer series. * C * The method used allows an accuracy of up to thirteen * C * decimal places through the use of large real arrays * C * and a single final division. LNCHF is a variable which * C * selects how the result should be represented. A '0' will * C * return the value in standard exponential form. A '1' * C * will return the natural log of the result. IP is an * C * integer variable that specifies how many array positions * C * are desired (usually 10 is sufficient). Setting IP=0 * C * causes the program to estimate the number of array * C * positions. * C * * C * The confluent hypergeometric function is the solution to * C * the differential equation: * C * * C * z M"(a;b;z) + (b-z) M'(a;b;z) - a M(a;b;z) = 0 * C * * C * Subprograms called: BITS, CHGF * C * * C **************************************************************** FUNCTION CONHYP (A,B,Z,LNCHF,IP) INTEGER LNCHF,I,BITS,IP COMPLEX*16 CHGF,A,B,Z,CONHYP DOUBLE PRECISION NTERM,FX,TERM1,MAX,TERM2,ANG,PI COMMON/CONSTS/ZERO,HALF,ONE,TWO,TEN,EPS DOUBLE PRECISION ZERO,HALF,ONE,TWO,TEN,EPS ZERO=0.0D0 PI=TWO*TWO*ATAN(ONE) IF (CDABS(Z) .NE. ZERO) THEN ANG=DATAN2(DIMAG(Z),DBLE(Z)) ELSE ANG=ONE ENDIF IF (DABS(ANG) .LT. (PI*HALF)) THEN ANG=ONE ELSE ANG=DSIN(DABS(ANG)-(PI*HALF))+ONE ENDIF MAX=ZERO NTERM=ZERO FX=ZERO TERM1=ZERO 10 NTERM=NTERM+ONE TERM2=CDABS((A+NTERM-1)*Z/((B+NTERM-1)*NTERM)) IF (TERM2 .EQ. ZERO) GOTO 20 IF (TERM2 .LT. ONE) THEN IF ((DBLE(A)+NTERM-1) .GT. ONE) THEN IF ((DBLE(B)+NTERM-1) .GT. ONE) THEN IF ((TERM2-TERM1) .LT. ZERO) THEN GOTO 20 ENDIF ENDIF ENDIF ENDIF FX=FX+DLOG(TERM2) IF (FX .GT. MAX) MAX=FX TERM1=TERM2 GOTO 10 20 MAX=MAX*2/(BITS()*6.93147181D-1) I=INT(MAX*ANG)+7 IF (I .LT. 5) I=5 IF (IP .GT. I) I=IP CONHYP=CHGF(A,B,Z,I,LNCHF) RETURN END C **************************************************************** C * * C * FUNCTION BITS * C * * C * * C * Description : Determines the number of significant figures * C * of machine precision to arrive at the size of the array * C * the numbers must must be stored in to get the accuracy * C * of the solution. * C * * C * Subprograms called: none * C * * C **************************************************************** INTEGER FUNCTION BITS() DOUBLE PRECISION BIT,BIT2 INTEGER COUNT COMMON/CONSTS/ZERO,HALF,ONE,TWO,TEN,EPS DOUBLE PRECISION ZERO,HALF,ONE,TWO,TEN,EPS BIT=ONE COUNT=0 10 COUNT=COUNT+1 BIT2=BIT*TWO BIT=BIT2+ONE IF ((BIT-BIT2) .NE. 0.0) GOTO 10 BITS=COUNT-1 RETURN END C **************************************************************** C * * C * FUNCTION CHGF * C * * C * * C * Description : Function that sums the Kummer series and * C * returns the solution of the confluent hypergeometric * C * function. * C * * C * Subprograms called: ARMULT, ARYDIV, BITS, CMPADD, CMPMUL * C * * C **************************************************************** FUNCTION CHGF (A,B,Z,L,LNCHF) PARAMETER (LENGTH=777) INTEGER L,I,BITS,BIT,LNCHF,NMACH,ICOUNT,IXCNT COMPLEX*16 A,B,Z,FINAL,CHGF DOUBLE PRECISION AR,AI,CR,CI,XR,XI,SUMR,SUMI,CNT,SIGFIG,MX1,MX2 DOUBLE PRECISION NUMR,NUMI,DENOMR,DENOMI,RMAX DOUBLE PRECISION QR1,QR2,QI1,QI2,AR2,AI2,CR2,CI2,XR2,XI2 DIMENSION SUMR(-1:LENGTH),SUMI(-1:LENGTH),NUMR(-1:LENGTH) DIMENSION NUMI(-1:LENGTH),DENOMR(-1:LENGTH),DENOMI(-1:LENGTH) DIMENSION QR1(-1:LENGTH),QR2(-1:LENGTH),QI1(-1:LENGTH), : QI2(-1:LENGTH) COMMON/CONSTS/ZERO,HALF,ONE,TWO,TEN,EPS DOUBLE PRECISION ZERO,HALF,ONE,TWO,TEN,EPS BIT=BITS() RMAX=TWO**(BIT/2) SIGFIG=TWO**(BIT/4) * * SET TO ZERO ANY ARGUMENTS WHICH ARE BELOW THE PRECISION OF THE * ALGORITHM. * AR2=DBLE(A)*SIGFIG AR=DINT(AR2) AR2=DNINT((AR2-AR)*RMAX) AI2=DIMAG(A)*SIGFIG AI=DINT(AI2) AI2=DNINT((AI2-AI)*RMAX) CR2=DBLE(B)*SIGFIG CR=DINT(CR2) CR2=DNINT((CR2-CR)*RMAX) CI2=DIMAG(B)*SIGFIG CI=DINT(CI2) CI2=DNINT((CI2-CI)*RMAX) XR2=DBLE(Z)*SIGFIG XR=DINT(XR2) XR2=DNINT((XR2-XR)*RMAX) XI2=DIMAG(Z)*SIGFIG XI=DINT(XI2) XI2=DNINT((XI2-XI)*RMAX) * * WARN THE USER THAT THE INPUT VALUE WAS SO CLOSE TO ZERO THAT IT * WAS SET EQUAL TO ZERO. * IF ((DBLE(A).NE.ZERO) .AND. (AR.EQ.ZERO) .AND. (AR2.EQ.ZERO)) : WRITE (6,*) ' WARNING - REAL PART OF A WAS SET TO ZERO' IF ((DIMAG(A).NE.ZERO) .AND. (AI.EQ.ZERO) .AND. (AI2.EQ.ZERO)) : WRITE (6,*) ' WARNING - IMAG PART OF A WAS SET TO ZERO' IF ((DBLE(B).NE.ZERO) .AND. (CR.EQ.ZERO) .AND. (CR2.EQ.ZERO)) : WRITE (6,*) ' WARNING - REAL PART OF B WAS SET TO ZERO' IF ((DIMAG(B).NE.ZERO) .AND. (CI.EQ.ZERO) .AND. (CI2.EQ.ZERO)) : WRITE (6,*) ' WARNING - IMAG PART OF B WAS SET TO ZERO' IF ((DBLE(Z).NE.ZERO) .AND. (XR.EQ.ZERO) .AND. (XR2.EQ.ZERO)) : WRITE (6,*) ' WARNING - REAL PART OF Z WAS SET TO ZERO' IF ((DIMAG(Z).NE.ZERO) .AND. (XI.EQ.ZERO) .AND. (XI2.EQ.ZERO)) : WRITE (6,*) ' WARNING - IMAG PART OF Z WAS SET TO ZERO' * * SCREENING OF THE CASE WHEN B IS ZERO OR A NEGATIVE INTEGER. * IF ((CR .EQ. ZERO) .AND. (CR2 .EQ. ZERO) .AND. : (CI .EQ. ZERO) .AND. (CI2 .EQ. ZERO)) THEN WRITE (6,*) ' ERROR-- ARGUMENT B WAS EQUAL TO ZERO' STOP END IF NMACH=INT(LOG10(TWO**INT(BITS()))) IF ((CI .EQ. ZERO) .AND. (CI2 .EQ. ZERO) .AND. : (DBLE(B) .LT. ZERO)) THEN IF (ABS(DBLE(B)-DBLE(NINT(DBLE(B)))) .LT. TEN**(-NMACH)) THEN WRITE (6,*) ' ERROR-- ARGUMENT B WAS A NEGATIVE INTEGER' STOP END IF END IF SUMR(-1)=ONE SUMI(-1)=ONE NUMR(-1)=ONE NUMI(-1)=ONE DENOMR(-1)=ONE DENOMI(-1)=ONE DO 100 I=0,L+1 SUMR(I)=ZERO SUMI(I)=ZERO NUMR(I)=ZERO NUMI(I)=ZERO DENOMR(I)=ZERO DENOMI(I)=ZERO 100 CONTINUE SUMR(1)=ONE NUMR(1)=ONE DENOMR(1)=ONE CNT=SIGFIG ICOUNT=-1 IF ((AI .EQ. ZERO) .AND. (AI2 .EQ. ZERO) .AND. : (DBLE(A) .LT. ZERO)) THEN IF (ABS(DBLE(A)-DBLE(NINT(DBLE(A)))) .LT. TEN**(-NMACH)) : ICOUNT=-NINT(DBLE(A)) END IF IXCNT=0 110 IF (SUMR(1) .LT. HALF) THEN MX1=SUMI(L+1) ELSE IF (SUMI(1) .LT. HALF) THEN MX1=SUMR(L+1) ELSE MX1=DMAX1(SUMR(L+1),SUMI(L+1)) ENDIF IF (NUMR(1) .LT. HALF) THEN MX2=NUMI(L+1) ELSE IF (NUMI(1) .LT. HALF) THEN MX2=NUMR(L+1) ELSE MX2=DMAX1(NUMR(L+1),NUMI(L+1)) ENDIF IF (MX1-MX2 .GT. 2.0) THEN IF (CR .GT. ZERO) THEN IF (CDABS(CMPLX(AR,AI)*CMPLX(XR,XI)/(CMPLX(CR,CI)*CNT)) : .LE. ONE) GOTO 190 ENDIF ENDIF IF (IXCNT .EQ. ICOUNT) GO TO 190 IXCNT=IXCNT+1 CALL CMPMUL(SUMR,SUMI,CR,CI,QR1,QI1,L,RMAX) CALL CMPMUL(SUMR,SUMI,CR2,CI2,QR2,QI2,L,RMAX) QR2(L+1)=QR2(L+1)-1 QI2(L+1)=QI2(L+1)-1 CALL CMPADD(QR1,QI1,QR2,QI2,SUMR,SUMI,L,RMAX) CALL ARMULT(SUMR,CNT,SUMR,L,RMAX) CALL ARMULT(SUMI,CNT,SUMI,L,RMAX) CALL CMPMUL(DENOMR,DENOMI,CR,CI,QR1,QI1,L,RMAX) CALL CMPMUL(DENOMR,DENOMI,CR2,CI2,QR2,QI2,L,RMAX) QR2(L+1)=QR2(L+1)-1 QI2(L+1)=QI2(L+1)-1 CALL CMPADD(QR1,QI1,QR2,QI2,DENOMR,DENOMI,L,RMAX) CALL ARMULT(DENOMR,CNT,DENOMR,L,RMAX) CALL ARMULT(DENOMI,CNT,DENOMI,L,RMAX) CALL CMPMUL(NUMR,NUMI,AR,AI,QR1,QI1,L,RMAX) CALL CMPMUL(NUMR,NUMI,AR2,AI2,QR2,QI2,L,RMAX) QR2(L+1)=QR2(L+1)-1 QI2(L+1)=QI2(L+1)-1 CALL CMPADD(QR1,QI1,QR2,QI2,NUMR,NUMI,L,RMAX) CALL CMPMUL(NUMR,NUMI,XR,XI,QR1,QI1,L,RMAX) CALL CMPMUL(NUMR,NUMI,XR2,XI2,QR2,QI2,L,RMAX) QR2(L+1)=QR2(L+1)-1 QI2(L+1)=QI2(L+1)-1 CALL CMPADD(QR1,QI1,QR2,QI2,NUMR,NUMI,L,RMAX) CALL CMPADD(SUMR,SUMI,NUMR,NUMI,SUMR,SUMI,L,RMAX) CNT=CNT+SIGFIG AR=AR+SIGFIG CR=CR+SIGFIG GOTO 110 190 CALL ARYDIV(SUMR,SUMI,DENOMR,DENOMI,FINAL,L,LNCHF,RMAX,BIT) CHGF=FINAL RETURN END C **************************************************************** C * * C * SUBROUTINE ARADD * C * * C * * C * Description : Accepts two arrays of numbers and returns * C * the sum of the array. Each array is holding the value * C * of one number in the series. The parameter L is the * C * size of the array representing the number and RMAX is * C * the actual number of digits needed to give the numbers * C * the desired accuracy. * C * * C * Subprograms called: none * C * * C **************************************************************** SUBROUTINE ARADD(A,B,C,L,RMAX) INTEGER L DOUBLE PRECISION A,B,C,Z,RMAX INTEGER EDIFF,I,J DIMENSION A(-1:*),B(-1:*),C(-1:*),Z(-1:777) COMMON/CONSTS/ZERO,HALF,ONE,TWO,TEN,EPS DOUBLE PRECISION ZERO,HALF,ONE,TWO,TEN,EPS DO 110 I=0,L+1 Z(I)=ZERO 110 CONTINUE EDIFF=DNINT(A(L+1)-B(L+1)) IF (DABS(A(1)) .LT. HALF .OR. EDIFF .LE. -L) GOTO 111 IF (DABS(B(1)) .LT. HALF .OR. EDIFF .GE. L) GOTO 113 GOTO 115 111 DO 112 I=-1,L+1 C(I)=B(I) 112 CONTINUE GOTO 311 113 DO 114 I=-1,L+1 C(I)=A(I) 114 CONTINUE GOTO 311 115 Z(-1)=A(-1) IF (DABS(A(-1)-B(-1)) .LT. HALF) GOTO 200 IF (EDIFF .GT. 0) THEN Z(L+1)=A(L+1) GOTO 233 ENDIF IF (EDIFF .LT. 0) THEN Z(L+1)=B(L+1) Z(-1)=B(-1) GOTO 266 ENDIF DO 120 I=1,L IF (A(I) .GT. B(I)) THEN Z(L+1)=A(L+1) GOTO 233 ENDIF IF (A(I) .LT. B(I)) THEN Z(L+1)=B(L+1) Z(-1)=B(-1) GOTO 266 ENDIF 120 CONTINUE GOTO 300 200 IF (EDIFF .GT. 0) GOTO 203 IF (EDIFF .LT. 0) GOTO 207 Z(L+1)=A(L+1) DO 201 I=L,1,-1 Z(I)=A(I)+B(I)+Z(I) IF (Z(I) .GE. RMAX) THEN Z(I)=Z(I)-RMAX Z(I-1)=ONE ENDIF 201 CONTINUE IF (Z(0) .GT. HALF) THEN DO 202 I=L,1,-1 Z(I)=Z(I-1) 202 CONTINUE Z(L+1)=Z(L+1)+ONE Z(0)=ZERO ENDIF GOTO 300 203 Z(L+1)=A(L+1) DO 204 I=L,1+EDIFF,-1 Z(I)=A(I)+B(I-EDIFF)+Z(I) IF (Z(I) .GE. RMAX) THEN Z(I)=Z(I)-RMAX Z(I-1)=ONE ENDIF 204 CONTINUE DO 205 I=EDIFF,1,-1 Z(I)=A(I)+Z(I) IF (Z(I) .GE. RMAX) THEN Z(I)=Z(I)-RMAX Z(I-1)=ONE ENDIF 205 CONTINUE IF (Z(0) .GT. HALF) THEN DO 206 I=L,1,-1 Z(I)=Z(I-1) 206 CONTINUE Z(L+1)=Z(L+1)+1 Z(0)=ZERO ENDIF GOTO 300 207 Z(L+1)=B(L+1) DO 208 I=L,1-EDIFF,-1 Z(I)=A(I+EDIFF)+B(I)+Z(I) IF (Z(I) .GE. RMAX) THEN Z(I)=Z(I)-RMAX Z(I-1)=ONE ENDIF 208 CONTINUE DO 209 I=0-EDIFF,1,-1 Z(I)=B(I)+Z(I) IF (Z(I) .GE. RMAX) THEN Z(I)=Z(I)-RMAX Z(I-1)=ONE ENDIF 209 CONTINUE IF (Z(0) .GT. HALF) THEN DO 210 I=L,1,-1 Z(I)=Z(I-1) 210 CONTINUE Z(L+1)=Z(L+1)+ONE Z(0)=ZERO ENDIF GOTO 300 233 IF (EDIFF .GT. 0) GOTO 243 DO 234 I=L,1,-1 Z(I)=A(I)-B(I)+Z(I) IF (Z(I) .LT. ZERO) THEN Z(I)=Z(I)+RMAX Z(I-1)=-ONE ENDIF 234 CONTINUE GOTO 290 243 DO 244 I=L,1+EDIFF,-1 Z(I)=A(I)-B(I-EDIFF)+Z(I) IF (Z(I) .LT. ZERO) THEN Z(I)=Z(I)+RMAX Z(I-1)=-ONE ENDIF 244 CONTINUE DO 245 I=EDIFF,1,-1 Z(I)=A(I)+Z(I) IF (Z(I) .LT. ZERO) THEN Z(I)=Z(I)+RMAX Z(I-1)=-ONE ENDIF 245 CONTINUE GOTO 290 266 IF (EDIFF .LT. 0) GOTO 276 DO 267 I=L,1,-1 Z(I)=B(I)-A(I)+Z(I) IF (Z(I) .LT. ZERO) THEN Z(I)=Z(I)+RMAX Z(I-1)=-ONE ENDIF 267 CONTINUE GOTO 290 276 DO 277 I=L,1-EDIFF,-1 Z(I)=B(I)-A(I+EDIFF)+Z(I) IF (Z(I) .LT. ZERO) THEN Z(I)=Z(I)+RMAX Z(I-1)=-ONE ENDIF 277 CONTINUE DO 278 I=0-EDIFF,1,-1 Z(I)=B(I)+Z(I) IF (Z(I) .LT. ZERO) THEN Z(I)=Z(I)+RMAX Z(I-1)=-ONE ENDIF 278 CONTINUE 290 IF (Z(1) .GT. HALF) GOTO 300 I=1 291 I=I+1 IF (Z(I) .LT. HALF .AND. I .LT. L+1) GOTO 291 IF (I .EQ. L+1) THEN Z(-1)=ONE Z(L+1)=ZERO GOTO 300 ENDIF 292 DO 293 J=1,L+1-I Z(J)=Z(J+I-1) 293 CONTINUE DO 294 J=L+2-I,L Z(J)=ZERO 294 CONTINUE Z(L+1)=Z(L+1)-I+1 300 DO 310 I=-1,L+1 C(I)=Z(I) 310 CONTINUE 311 IF (C(1) .LT. HALF) THEN C(-1)=ONE C(L+1)=ZERO ENDIF RETURN END C **************************************************************** C * * C * SUBROUTINE ARSUB * C * * C * * C * Description : Accepts two arrays and subtracts each element * C * in the second array from the element in the first array * C * and returns the solution. The parameters L and RMAX are * C * the size of the array and the number of digits needed for * C * the accuracy, respectively. * C * * C * Subprograms called: ARADD * C * * C **************************************************************** SUBROUTINE ARSUB(A,B,C,L,RMAX) INTEGER L,I DOUBLE PRECISION A,B,C,B2,RMAX DIMENSION A(-1:*),B(-1:*),C(-1:*),B2(-1:777) COMMON/CONSTS/ZERO,HALF,ONE,TWO,TEN,EPS DOUBLE PRECISION ZERO,HALF,ONE,TWO,TEN,EPS DO 100 I=-1,L+1 B2(I)=B(I) 100 CONTINUE B2(-1)=(-ONE)*B2(-1) CALL ARADD(A,B2,C,L,RMAX) RETURN END C **************************************************************** C * * C * SUBROUTINE ARMULT * C * * C * * C * Description : Accepts two arrays and returns the product. * C * L and RMAX are the size of the arrays and the number of * C * digits needed to represent the numbers with the required * C * accuracy. * C * * C * Subprograms called: none * C * * C **************************************************************** SUBROUTINE ARMULT(A,B,C,L,RMAX) INTEGER L DOUBLE PRECISION A,B,C,Z,B2,CARRY,RMAX,RMAX2 DIMENSION A(-1:*),C(-1:*),Z(-1:777) INTEGER I COMMON/CONSTS/ZERO,HALF,ONE,TWO,TEN,EPS DOUBLE PRECISION ZERO,HALF,ONE,TWO,TEN,EPS RMAX2=ONE/RMAX Z(-1)=DSIGN(ONE,B)*A(-1) B2=DABS(B) Z(L+1)=A(L+1) DO 100 I=0,L Z(I)=ZERO 100 CONTINUE IF (B2 .LE. EPS .OR. A(1) .LE. EPS) THEN Z(-1)=ONE Z(L+1)=ZERO GOTO 198 ENDIF DO 110 I=L,1,-1 Z(I)=A(I)*B2+Z(I) IF (Z(I) .GE. RMAX) THEN CARRY=DINT(Z(I)/RMAX) Z(I)=Z(I)-CARRY*RMAX Z(I-1)=CARRY ENDIF 110 CONTINUE IF (Z(0) .LT. HALF) GOTO 150 DO 120 I=L,1,-1 Z(I)=Z(I-1) 120 CONTINUE Z(L+1)=Z(L+1)+ONE Z(0)=ZERO 150 CONTINUE 198 DO 199 I=-1,L+1 C(I)=Z(I) 199 CONTINUE IF (C(1) .LT. HALF) THEN C(-1)=ONE C(L+1)=ZERO ENDIF RETURN END C **************************************************************** C * * C * SUBROUTINE CMPADD * C * * C * * C * Description : Takes two arrays representing one real and * C * one imaginary part, and adds two arrays representing * C * another complex number and returns two array holding the * C * complex sum. * C * (CR,CI) = (AR+BR, AI+BI) * C * * C * Subprograms called: ARADD * C * * C **************************************************************** SUBROUTINE CMPADD(AR,AI,BR,BI,CR,CI,L,RMAX) INTEGER L DOUBLE PRECISION AR,AI,BR,BI,CR,CI,RMAX DIMENSION AR(-1:*),AI(-1:*),BR(-1:*),BI(-1:*) DIMENSION CR(-1:*),CI(-1:*) CALL ARADD(AR,BR,CR,L,RMAX) CALL ARADD(AI,BI,CI,L,RMAX) RETURN END C **************************************************************** C * * C * SUBROUTINE CMPSUB * C * * C * * C * Description : Takes two arrays representing one real and * C * one imaginary part, and subtracts two arrays representing * C * another complex number and returns two array holding the * C * complex sum. * C * (CR,CI) = (AR+BR, AI+BI) * C * * C * Subprograms called: ARADD * C * * C **************************************************************** SUBROUTINE CMPSUB(AR,AI,BR,BI,CR,CI,L,RMAX) INTEGER L DOUBLE PRECISION AR,AI,BR,BI,CR,CI,RMAX DIMENSION AR(-1:*),AI(-1:*),BR(-1:*),BI(-1:*) DIMENSION CR(-1:*),CI(-1:*) CALL ARSUB(AR,BR,CR,L,RMAX) CALL ARSUB(AI,BI,CI,L,RMAX) RETURN END C **************************************************************** C * * C * SUBROUTINE CMPMUL * C * * C * * C * Description : Takes two arrays representing one real and * C * one imaginary part, and multiplies it with two arrays * C * representing another complex number and returns the * C * complex product. * C * * C * Subprograms called: ARMULT, ARSUB, ARADD * C * * C **************************************************************** SUBROUTINE CMPMUL(AR,AI,BR,BI,CR,CI,L,RMAX) INTEGER L DOUBLE PRECISION AR,AI,BR,BI,CR,CI,D1,D2,RMAX DIMENSION AR(-1:*),AI(-1:*),CR(-1:*),CI(-1:*) DIMENSION D1(-1:777),D2(-1:777) CALL ARMULT(AR,BR,D1,L,RMAX) CALL ARMULT(AI,BI,D2,L,RMAX) CALL ARSUB(D1,D2,CR,L,RMAX) CALL ARMULT(AR,BI,D1,L,RMAX) CALL ARMULT(AI,BR,D2,L,RMAX) CALL ARADD(D1,D2,CI,L,RMAX) RETURN END C **************************************************************** C * * C * SUBROUTINE ARYDIV * C * * C * * C * Description : Returns the double precision complex number * C * resulting from the division of four arrays, representing * C * two complex numbers. The number returned will be in one * C * two different forms: either standard scientific or as * C * the natural log of the number. * C * * C * Subprograms called: CONV21, CONV12, EADD, ECPDIV, EMULT * C * * C **************************************************************** SUBROUTINE ARYDIV(AR,AI,BR,BI,C,L,LNCHF,RMAX,BIT) INTEGER L,BIT,REXP,IR10,II10,LNCHF COMPLEX*16 C DOUBLE PRECISION AR,AI,BR,BI,PHI,N1,N2,N3,E1,E2,E3,RR10,RI10,X DOUBLE PRECISION AE,BE,X1,X2,DUM1,DUM2,CE,RMAX DIMENSION AR(-1:*),AI(-1:*),BR(-1:*),BI(-1:*) DIMENSION AE(2,2),BE(2,2),CE(2,2) COMMON/CONSTS/ZERO,HALF,ONE,TWO,TEN,EPS DOUBLE PRECISION ZERO,HALF,ONE,TWO,TEN,EPS REXP=BIT/2 X=REXP*(AR(L+1)-2) RR10=X*DLOG10(TWO)/DLOG10(TEN) IR10=INT(RR10) RR10=RR10-IR10 X=REXP*(AI(L+1)-2) RI10=X*DLOG10(TWO)/DLOG10(TEN) II10=INT(RI10) RI10=RI10-II10 DUM1=DSIGN(AR(1)*RMAX*RMAX+AR(2)*RMAX+AR(3),AR(-1)) DUM2=DSIGN(AI(1)*RMAX*RMAX+AI(2)*RMAX+AI(3),AI(-1)) DUM1=DUM1*10**RR10 DUM2=DUM2*10**RI10 CALL CONV12(DCMPLX(DUM1,DUM2),AE) AE(1,2)=AE(1,2)+IR10 AE(2,2)=AE(2,2)+II10 X=REXP*(BR(L+1)-2) RR10=X*DLOG10(TWO)/DLOG10(TEN) IR10=INT(RR10) RR10=RR10-IR10 X=REXP*(BI(L+1)-2) RI10=X*DLOG10(TWO)/DLOG10(TEN) II10=INT(RI10) RI10=RI10-II10 DUM1=DSIGN(BR(1)*RMAX*RMAX+BR(2)*RMAX+BR(3),BR(-1)) DUM2=DSIGN(BI(1)*RMAX*RMAX+BI(2)*RMAX+BI(3),BI(-1)) DUM1=DUM1*10**RR10 DUM2=DUM2*10**RI10 CALL CONV12(DCMPLX(DUM1,DUM2),BE) BE(1,2)=BE(1,2)+IR10 BE(2,2)=BE(2,2)+II10 CALL ECPDIV(AE,BE,CE) IF (LNCHF .EQ. 0) THEN CALL CONV21(CE,C) ELSE CALL EMULT(CE(1,1),CE(1,2),CE(1,1),CE(1,2),N1,E1) CALL EMULT(CE(2,1),CE(2,2),CE(2,1),CE(2,2),N2,E2) CALL EADD(N1,E1,N2,E2,N3,E3) N1=CE(1,1) E1=CE(1,2)-CE(2,2) X2=CE(2,1) IF (E1 .GT. 74.0D0) THEN X1=1.0D75 ELSEIF (E1 .LT. -74.0D0) THEN X1=0 ELSE X1=N1*(10**E1) ENDIF PHI=DATAN2(X2,X1) C=DCMPLX(HALF*(DLOG(N3)+E3*DLOG(TEN)),PHI) ENDIF RETURN END C **************************************************************** C * * C * SUBROUTINE EMULT * C * * C * * C * Description : Takes one base and exponent and multiplies it * C * by another numbers base and exponent to give the product * C * in the form of base and exponent. * C * * C * Subprograms called: none * C * * C **************************************************************** SUBROUTINE EMULT(N1,E1,N2,E2,NF,EF) DOUBLE PRECISION N1,E1,N2,E2,NF,EF COMMON/CONSTS/ZERO,HALF,ONE,TWO,TEN,EPS DOUBLE PRECISION ZERO,HALF,ONE,TWO,TEN,EPS NF=N1*N2 EF=E1+E2 IF (DABS(NF) .GE. TEN) THEN NF=NF/TEN EF=EF+ONE ENDIF RETURN END C **************************************************************** C * * C * SUBROUTINE EDIV * C * * C * * C * Description : returns the solution in the form of base and * C * exponent of the division of two exponential numbers. * C * * C * Subprograms called: none * C * * C **************************************************************** SUBROUTINE EDIV(N1,E1,N2,E2,NF,EF) DOUBLE PRECISION N1,E1,N2,E2,NF,EF COMMON/CONSTS/ZERO,HALF,ONE,TWO,TEN,EPS DOUBLE PRECISION ZERO,HALF,ONE,TWO,TEN,EPS NF=N1/N2 EF=E1-E2 IF ((DABS(NF) .LT. ONE) .AND. (NF .NE. ZERO)) THEN NF=NF*TEN EF=EF-ONE ENDIF RETURN END C **************************************************************** C * * C * SUBROUTINE EADD * C * * C * * C * Description : Returns the sum of two numbers in the form * C * of a base and an exponent. * C * * C * Subprograms called: none * C * * C **************************************************************** SUBROUTINE EADD(N1,E1,N2,E2,NF,EF) DOUBLE PRECISION N1,E1,N2,E2,NF,EF,EDIFF,THIRSIX COMMON/CONSTS/ZERO,HALF,ONE,TWO,TEN,EPS DOUBLE PRECISION ZERO,HALF,ONE,TWO,TEN,EPS EDIFF=E1-E2 THIRSIX=(ONE+TWO)*(TWO+TEN) IF (EDIFF .GT. THIRSIX) THEN NF=N1 EF=E1 ELSE IF (EDIFF .LT. -THIRSIX) THEN NF=N2 EF=E2 ELSE NF=N1*(TEN**EDIFF)+N2 EF=E2 400 IF (DABS(NF) .LT. TEN) GOTO 410 NF=NF/TEN EF=EF+ONE GOTO 400 410 IF ((DABS(NF) .GE. ONE) .OR. (NF .EQ. ZERO)) GOTO 420 NF=NF*TEN EF=EF-ONE GOTO 410 ENDIF 420 RETURN END C **************************************************************** C * * C * SUBROUTINE ESUB * C * * C * * C * Description : Returns the solution to the subtraction of * C * two numbers in the form of base and exponent. * C * * C * Subprograms called: EADD * C * * C **************************************************************** SUBROUTINE ESUB(N1,E1,N2,E2,NF,EF) DOUBLE PRECISION N1,E1,N2,E2,NF,EF COMMON/CONSTS/ZERO,HALF,ONE,TWO,TEN,EPS DOUBLE PRECISION ZERO,HALF,ONE,TWO,TEN,EPS CALL EADD(N1,E1,N2*(-ONE),E2,NF,EF) RETURN END C **************************************************************** C * * C * SUBROUTINE CONV12 * C * * C * * C * Description : Converts a number from complex notation to a * C * form of a 2x2 real array. * C * * C * Subprograms called: none * C * * C **************************************************************** SUBROUTINE CONV12(CN,CAE) COMPLEX*16 CN DOUBLE PRECISION CAE DIMENSION CAE(2,2) COMMON/CONSTS/ZERO,HALF,ONE,TWO,TEN,EPS DOUBLE PRECISION ZERO,HALF,ONE,TWO,TEN,EPS CAE(1,1)=DBLE(CN) CAE(1,2)=ZERO 300 IF (DABS(CAE(1,1)) .LT. TEN) GOTO 310 CAE(1,1)=CAE(1,1)/TEN CAE(1,2)=CAE(1,2)+ONE GOTO 300 310 IF ((DABS(CAE(1,1)) .GE. ONE) .OR. (CAE(1,1) .EQ. ZERO)) : GOTO 320 CAE(1,1)=CAE(1,1)*TEN CAE(1,2)=CAE(1,2)-ONE GOTO 310 320 CAE(2,1)=DIMAG(CN) CAE(2,2)=ZERO 330 IF (DABS(CAE(2,1)) .LT. TEN) GOTO 340 CAE(2,1)=CAE(2,1)/TEN CAE(2,2)=CAE(2,2)+ONE GOTO 330 340 IF ((DABS(CAE(2,1)) .GE. ONE) .OR. (CAE(2,1) .EQ. ZERO)) : GOTO 350 CAE(2,1)=CAE(2,1)*TEN CAE(2,2)=CAE(2,2)-ONE GOTO 340 350 RETURN END C **************************************************************** C * * C * SUBROUTINE CONV21 * C * * C * * C * Description : Converts a number represented in a 2x2 real * C * array to the form of a complex number. * C * * C * Subprograms called: none * C * * C **************************************************************** SUBROUTINE CONV21(CAE,CN) DOUBLE PRECISION CAE COMPLEX*16 CN DIMENSION CAE(2,2) COMMON/CONSTS/ZERO,HALF,ONE,TWO,TEN,EPS DOUBLE PRECISION ZERO,HALF,ONE,TWO,TEN,EPS IF (CAE(1,2) .GT. 75 .OR. CAE(2,2) .GT. 75) THEN CN=CMPLX(1.0D75,1.0D75) ELSE IF (CAE(2,2) .LT. -75) THEN CN=DCMPLX(CAE(1,1)*(10**CAE(1,2)),ZERO) ELSE CN=DCMPLX(CAE(1,1)*(10**CAE(1,2)),CAE(2,1)*(10**CAE(2,2))) ENDIF RETURN END C **************************************************************** C * * C * SUBROUTINE ECPMUL * C * * C * * C * Description : Multiplies two numbers which are each * C * represented in the form of a two by two array and returns * C * the solution in the same form. * C * * C * Subprograms called: EMULT, ESUB, EADD * C * * C **************************************************************** SUBROUTINE ECPMUL(A,B,C) DOUBLE PRECISION A,B,C,N1,E1,N2,E2,C2 DIMENSION A(2,2),B(2,2),C(2,2),C2(2,2) CALL EMULT(A(1,1),A(1,2),B(1,1),B(1,2),N1,E1) CALL EMULT(A(2,1),A(2,2),B(2,1),B(2,2),N2,E2) CALL ESUB(N1,E1,N2,E2,C2(1,1),C2(1,2)) CALL EMULT(A(1,1),A(1,2),B(2,1),B(2,2),N1,E1) CALL EMULT(A(2,1),A(2,2),B(1,1),B(1,2),N2,E2) CALL EADD(N1,E1,N2,E2,C(2,1),C(2,2)) C(1,1)=C2(1,1) C(1,2)=C2(1,2) RETURN END C **************************************************************** C * * C * SUBROUTINE ECPDIV * C * * C * * C * Description : Divides two numbers and returns the solution. * C * All numbers are represented by a 2x2 array. * C * * C * Subprograms called: EADD, ECPMUL, EDIV, EMULT * C * * C **************************************************************** SUBROUTINE ECPDIV(A,B,C) DOUBLE PRECISION A,B,C,N1,E1,N2,E2,B2,N3,E3,C2 DIMENSION A(2,2),B(2,2),C(2,2),B2(2,2),C2(2,2) COMMON/CONSTS/ZERO,HALF,ONE,TWO,TEN,EPS DOUBLE PRECISION ZERO,HALF,ONE,TWO,TEN,EPS B2(1,1)=B(1,1) B2(1,2)=B(1,2) B2(2,1)=-ONE*B(2,1) B2(2,2)=B(2,2) CALL ECPMUL(A,B2,C2) CALL EMULT(B(1,1),B(1,2),B(1,1),B(1,2),N1,E1) CALL EMULT(B(2,1),B(2,2),B(2,1),B(2,2),N2,E2) CALL EADD(N1,E1,N2,E2,N3,E3) CALL EDIV(C2(1,1),C2(1,2),N3,E3,C(1,1),C(1,2)) CALL EDIV(C2(2,1),C2(2,2),N3,E3,C(2,1),C(2,2)) RETURN END C **************************************************************** C * * C * BLOCK DATA BLDAT1 * C * * C * * C * Description : Sets of frequently used numbers in a common * C * block. This makes it easier to convert the code to a * C * single precision version. * C * * C **************************************************************** BLOCK DATA BLDAT1 C COMMON/CONSTS/ZERO,HALF,ONE,TWO,TEN,EPS DOUBLE PRECISION ZERO,HALF,ONE,TWO,TEN,EPS DATA ZERO,HALF,ONE,TWO,TEN,EPS/0.0D0,0.5D0,1.0D0,2.0D0, : 10.0D0,1.0D-10/ END