|
DataMuseum.dkPresents historical artifacts from the history of: RC4000/8000/9000 |
This is an automatic "excavation" of a thematic subset of
See our Wiki for more about RC4000/8000/9000 Excavated with: AutoArchaeologist - Free & Open Source Software. |
top - metrics - download
Length: 33792 (0x8400)
Types: TextFile
Names: »cpc11«
└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ. Detaljer om "HC8000" projekt.
└─⟦0364f57e3⟧
└─⟦this⟧ »cpc11«
c ACYA MCHF77. A GENERAL MULTI-CONFIGURATION HARTREE-FOCK PROGRAM.
c C.F. FISCHER.
c REF. IN COMP. PHYS. COMMUN. 14 (1978) 145
C ------------------------------------------------------------------
C M C H F 7 7
C ------------------------------------------------------------------
C
C A MULTI-CONFIGURATION HARTREE FOCK PROGRAM FOR ATOMS FOR THE
C IBM SYSTEM 360/370 (DOUBLE PRECISION) BY
C
C CHARLOTTE FROESE FISCHER
C DEPARTMENT OF COMPUTER SCIENCE
C PENN STATE UNIVERSITY
C UNIVERSITY PARK, PA 16801
C
C
C REVISION OF : MCHF72
C COMPUTER PHYSICS COMMUNICATIONS 4 (1972) 107
C DEVELOPED AT: THE UNIVERSITY OF WATERLOO
C WATERLOO, ONTARIO
C SUPPORTED BY: THE NATIONAL RESEARCH COUNCIL OF CANADA
C
C
C REVISION OF : MULTI-CONFIGURATION HARTREE-FOCK
C COMPUTER PHYSICS COMMUNICATIONS 1 (1969) 151
C DEVELOPED AT: THE UNIVERSITY OF BRITISH COLUMBIA
C VANCOUVER, BC
C SUPPORTED BY: THE NATIONAL RESEARCH COUNCIL OF CANADA
C
C
C
C
C
C
C THE PRESENT PROGRAM WAS DEVELOPED IN PART WHILE AT THE UNIVERSITY
C OF WATERLOO AND SUPPORTED BY A NATIONAL RESEARCH COUNCIL OF
C CANADA GRANT, AND IN PART AT PENN STATE SUPPORTED BY A US ERDA
C GRANT.
C
C ------------------------------------------------------------------
C 1 I N T R O D U C T I O N
C ------------------------------------------------------------------
C
C A MULTI-CONFIGURATION HARTREE-FOCK APPROXIMATION IS A WAVEFUNCTION
C OF THE FORM
C
C W(CLS) = SUM ON I $ WT(I) * ! CONFIG(I) > $
C (I .LE. NCFG)
C
C WHERE CLS SPECIFIES THE STATE AND ! > DESIGNATES A
C A CONFIGURATION STATE FUNCTION DEFINED IN TERMS OF SPIN-ORBITALS
C OF THE FORM
C
C PHI = (1/R)P(NL;R)<SPHERICAL HARMONIC><SPIN FUNCTION>
C
C BOTH THE COEFFICIENTS WT(I) IN THE EXPANSION AND THE RADIAL
C FUNCTIONS P(NL;R) ARE DETERMINED VARIATIONALLY SO AS TO LEAVE
C THE TOTAL ENERGY OF THE SYSTEM STATIONARY WITH RESPECT TO ALL
C PERTURBATIONS SATISFYING CERTAIN ORTHOGONALITY CONSTRAINTS. ALL
C ORBITALS WITHIN A CONFIGURATION STATE FUNCTION ARE ASSUMED TO BE
C ORTHOGONAL BUT SOME FLEXIBILITY IS ALLOWED IN THE ORTHOGONALITY
C CONSTRAINT OF ORBITALS IN DIFFERENT CONFIGURATIONS. THE
C CONSTRAINTS MUST BE SUCH THAT THE MOST GENERAL FORM OF THE ENERGY
C EXPRESSION IS
C
C E(ALS) = SUM ON I $ WT(I)**2 * EAV(I) $
C (I .LE. NCFG)
C
C + SUM ON M $A(M)*WT(J1(M))*WT(J2(M))*FK(M)(I1(M),I2(M))$
C (M .LE. NF)
C
C + SUM ON M $B(M)*WT(J1(M))*WT(J2(M))*GK(M)(I1(M),I2(M))$
C (M .LE. NG)
C
C + SUM ON M $D(M)*WT(J1(M))*WT(J2(M))
C (M .LE. NR)
C *RK(M)(I1(M),I2(M),I3(M),I4(M))
C
C *<IO(M)!JO(M)>**P(M)
C
C + SUM ON M $C(M)*WT(J1(M))*WT(J2(M))
C (M .LE. NL)
C *L(I1(M),I2(M))*<IO(M)!JO(M)>**P(M)
C
C WHERE EAV(I) IS THE AVERAGE ENERGY FOR THE I'TH CONFIGURATION,
C FK(I,J), GK(I,J) AND RK(I1, I2, I3, I4) ARE SLATER INTEGRALS AND
C THE CONTRIBUTION FROM THE ONE-ELECTRON PART OF THE HAMILTONIAN IS
C
C I(NL,N'L) = -(1/2)L(NL,N'L)
C
C
C NOTE THAT CONTRIBUTIONS FROM INTERACTIONS BETWEEN CONFIGURATIONS
C CONTRIBUTE TWICE TO THE ENERGY EXPRESSION AND SO HAVE A FACTOR OF
C TWO INCLUDED IN THE COEFFICIENT. ALSO, BY THE SYMMETRY PROPERTIES
C OF THE RK INTEGRALS, CERTAIN INTERACTION INTEGRALS MAY BE EXPRES-
C SED AS FK OR GK INTEGRALS. THE LATTER ARE TREATED MORE EFFIC-
C IENTLY BY THE PROGRAM AND ARE THE PREFERRED FORM.
C
C THE ABOVE EXPRESSION FOR THE ENERGY ASSUMES THE CONFIGURATION
C STATE FUNCTIONS AND THE RADIAL FUNCTIONS HAVE EACH BEEN ORDERED
C SO THAT THEY CAN BE REFERENCED BY THEIR INDEX IN THE LIST. THE
C INDICES FOR WT ALWAYS REFER TO CONFIGURATIONS WHEREAS THE
C INDICES FOR THE SLATER INTEGRALS, L AND OVERLAP INTEGRALS ARE
C ALWAYS RADIAL FUNCTIONS.
C
C THE STATIONARY CONDITIONS FOR THE RADIAL FUNCTIONS LEAD TO A
C SYSTEM OF COUPLED INTEGRODIFFERENTIAL EQUATIONS THAT DEPEND ON THE
C MIXING COEFFICIENTS WT(I). THESE ARE THE MCHF EQUATIONS FOR THE
C RADIAL FUNCTIONS. AS IN A CONFIGURATION INTERACTION CALCULATION,
C THE STATIONARY CONDITIONS FOR THE MIXING COEFFICIENTS ARE SOLU-
C TIONS OF A SECULAR PROBLEM, (H - E)WT = 0, WHERE H IS THE INTER-
C ACTION MATRIX, AND E THE TOTAL ENERGY. CLEARLY THE ENTRIES IN
C THE MATRIX DEPEND ON THE RADIAL FUNCTIONS. CONSEQUENTLY THE
C MCHF EQUATIONS AND THE SECULAR PROBLEM ARE COUPLED . TOGETHER
C THEY DEFINE THE MCHF PROBLEM. A MORE DETAILED DISCUSSION OF THE
C DERIVATION OF HF EQUATIONS, THEIR PROPERTIES, AND THEIR NUMERICAL
C SOLUTION IS CONTAINED IN THE BOOK, "THE HARTREE-FOCK METHOD FOR
C ATOMS - A NUMERICAL APPROACH", PUBLISHED BY WILEY INTERSCIENCE,
C NEW YORK, 1977. SECTION AND EQUATION REFERENCES IN THE DOCUMENT-
C ATION OF THIS PROGRAM REFER TO THIS BOOK.
C
C THIS PROGRAM DETERMINES THE RADIAL FUNCTIONS AS WELL AS THE
C MIXING COEFFICIENTS FOR A MULTICONFIGURATION APPROXIMATION TO THE
C TOTAL WAVEFUNCTION.
C
C ------------------------------------------------------------------
C 2 I N P U T D A T A
C ------------------------------------------------------------------
C
C THE PROGRAM ALLOWS DATA TO BE READ FROM, AND OUTPUT ROUTED TO,
C A VARIETY OF UNITS. UNLESS INDICATED OTHERWISE, THE INPUT DATA
C HERE WILL BE DESCRIBED ASSUMING ALL DATA IS READ FROM THE SAME
C UNIT, NAMELY INPUT UNIT 1. IF SECTIONS ARE TO BE READ FROM OTHER
C UNITS, THE DATA ON A GIVEN UNIT SHOULD APPEAR IN THE SAME ORDER.
C
CARD 0. IUC, IUD, IUF, IUH, OUC, OUD, OUF, OUH IN FORMAT(4I3)
C
C IUC - INPUT UNIT FOR THE CONFIGURATION CARDS (CARD 2.)
C IUD - INPUT UNIT FOR THE DATA CARDS DEFINING THE ENERGY EXPRESSION
C (CARDS 4., 5., 6., AND 7.) IF ANY.
C IUF - INPUT UNIT FOR THE FUNCTION CARDS (CARD 8.) IF ANY.
C IUH - INPUT UNIT FOR THE HAMILTONIAN MATRIX (CARD 10.) IF ANY.
C THIS UNIT NUMBER SHOULD BE ZERO IF NO MATRIX IS TO BE READ.
C OUC - OUTPUT UNIT FOR A HEADER AND CONFIGURATION CARDS.
C OUD - OUTPUT UNIT FOR VALUES OF THE SLATER INTEGRALS ENTERING INTO
C THE ENERGY EXPRESSION. IF OUD > 0, THE FULL ENERGY EXP-
C RESSION WILL ALSO BE PRINTED; OTHERWISE THIS PRINTING IS
C OMITTED.
C OUF - OUTPUT UNIT FOR THE FUNCTION VALUES
C OUH - OUTPUT UNIT FOR THE HAMILTONIAN MATRIX
C
C IN EACH CASE, THE OUTPUT WILL BE OMITTED IF THE UNIT NUMBER IS SET
C TO ZERO.
C
CARD 1. ATOM, TERM, Z, NO, NWF, NIT, NCFG, NF, NG, NR, NL, ORTHO, OMIT,
C NEW IN THE FORMAT(2A6,F6.0,I6,7I3,2L3,I3)
C
C ATOM - IDENTIFYING LABEL
C TERM - IDENTIFYING LABEL
C Z - ATOMIC NUMBER
C NO - MAXIMUM NUMBER OF POINTS IN THE RANGE OF THE FUNCTION
C VARYING FROM 160 FOR A SMALL ATOM TO 220 FOR A LARGE
C ATOM
C NWF - NUMBER OF FUNCTIONS
C NIT - NUMBER OF FUNCTIONS TO BE MADE SELF-CONSISTENT WITH THE
C REMAINING INNER CORE TO BE KEPT "FROZEN"
C NCFG - NUMBER OF CONFIGURATIONS
C NF - NUMBER OF FK INTEGRALS IN THE EXPRESSION FOR THE ENERGY
C NG - NUMBER OF GK INTEGRALS IN THE EXPRESSION FOR THE ENERGY
C NR - NUMBER OF RK INTEGRALS IN THE EXPRESSION FOR THE ENERGY
C NL - NUMBER OF L INTEGRALS IN THE EXPRESSION FOR THE ENERGY
C ORTHO- LOGICAL VARIABLE: IF .FALSE. ONLY THE ORBITALS WITHIN A
C CONFIGURATION WILL BE MADE ORTHOGONAL, OTHERWISE ALL
C ORBITALS WITH DIFFERENT NL VALUES WILL BE MADE ORTHOGONAL.
C WHEN NON-ORTHOGONAL ORBITALS ARE PRESENT, ORTHO MUST BE
C .TRUE. IF A CORE ORBITAL IS TO BE DETERMINED WHICH SIMULT-
C ANEOUSLY MUST BE ORTHOGONAL TO TWO NON-ORTHOGONAL ORBITALS
C WITH THE SAME NL VALUES. ORTHO MAY BE .FALSE. IF THE
C COMMON CORE IS FIXED. THE PROGRAM CANNOT DEAL WITH CASES
WHERE AN ORBITAL TO BE DETERMINED MUST BE ORTHOGONAL TO
C MORE THAN A SINGLE PAIR OF NON-ORTHOGONAL ORBITALS.
C OMIT - LOGICAL VARIABLE: IF .FALSE. ORTHOGONALIZATION WILL OCCUR
C ONLY AT THE END OF AN SCF CYCLE (WEAK ORTHOGONALITY),
C THOUGH STRICT (STRONG) ORTHOGONALITY WILL BE MAINTAINED
C WITH A FIXED CORE. OTHERWISE, FUNCTIONS WILL BE KEPT
C ORTHOGONAL AT ALL STAGES (STRONG ORTHOGONALITY).
C NEW - THE NUMBER OF NEW CONFIGURATIONS. IF NEW = 0, ALL
C CONFIGURATIONS WILL BE TREATED AS NEW, BUT WHEN NEW > 0,
C AN MXM HAMILTONIAN MATRIX, WHERE M = (NCFG - NEW), WILL BE
C LEFT UNCHANGED OR READ IN AS DATA IF IUH > 0. IT IS
C ASSUMED THE NEW CONFIGURATIONS HAVE BEEN ADDED AT THE END
C OF THE LIST OF CONFIGURATIONS. THIS OPTION ALLOWS A CALC-
C ULATION TO TREAT THE FIRST M CONFIGURATIONS AS PART OF A
C FIXED OR FROZEN CORE. THE RADIAL FUNCTIONS DEFINING THE
C ASSOCIATED ENTRIES IN THE HAMILTONIAN MATRIX SHOULD THEN
C ALSO BE KEPT FIXED.
C
C THE MAXIMUM ALLOWED VALUE OF MANY OF THE ABOVE VARIABLES DEPENDS
C ON THE DIMENSIONS OF ARRAYS. THE MAXIMUM VALUES FOR THE PRESENT
C PROGRAM ARE GIVEN BELOW:
C
C VARIABLE MAXIMUM VALUE
C -------- -------------
C
C NO 220
C NWF 20
C NCFG 40
C NF + NG 200
C NR 300
C NL 20
C
CARD 2. FOR EACH I, I = 1,NCFG A CONFIGURATION CARD WITH
C
C CONFIG1,CONFIG2,CONFIG3,CONFIG4, WT, WTL IN FORMAT( 4A6, F10.8, L1)
C
C CONFIG - IDENTIFYING LABEL FOR THE CONFIGURATION
C WT - THE COEFFICIENT OR WEIGHT FOR THE CONFIGURATION:
C IF OMITTED, ALL CONFIGURATIONS HAVE EQUAL WEIGHT.
C THE WEIGHTS ARE NORMALIZED TO UNITY BY THE PROGRAM.
C IN AN MCHF CALCULATION IT IS IMPORTANT THAT THE INITIAL
C EXPECTED OCCUPATION NUMBER OF ALL FUNCTIONS TO BE
C DETERMINED BE DIFFERENT FROM ZERO, AND THAT FUNCTIONS
C CONSTRAINED BY ORTHOGONALITY HAVE DIFFERING OCCUPATION
C NUMBERS INITIALLY IF THEY WILL DIFFER IN THE FINAL
C ANSWER.
C WTL - LOGICAL VARIABLE: IF .TRUE. THE THE WEIGHT IS TO BE
C LEFT UNCHANGED FROM THE PREVIOUS CASE, OTHERWISE THE
C VALUE OF WT IS TO BE USED.
C
CARD 3. FOR EACH I, I = 1,NWF A CARD WITH
C EL(I),N(I),L(I),S(I),METH(I),ACC(I),IND(I),(QC(I,J),J=1,NCFG)
C IN THE FORMAT(A3, 2I3, F6.2, I3, F3.1, I3, 15F3.0/(24X,15F3.0))
C
C EL - IDENTIFYING LABEL FOR THE ELECTRON
C N - PRINCIPAL QUANTUM NUMBER
C L - ANGULAR QUANTUM NUMBER
C S - SCREENING PARAMETER
C METH- METHOD TO BE USED FOR SOLVING THE DIFFERENTIAL EQUATION
C 1 - METHOD 1 SOLVES A SINGLE BOUNDARY VALUE PROBLEM FOR AN
C ACCEPTABLE SOLUTION WHICH NEED NOT BE NORMALIZED.
C 2 - METHOD 2 SOLVES A SINGLE BOUNDARY VALUE PROBLEM FOR AN
C ACCEPTABLE SOLUTION WHICH IS NORMALIZED TO FIRST ORDER.
C IF THE EXCHANGE FUNCTION IS IDENTICALLY ZERO, THE PROG-
C RAM WILL AUTOMATICALLY SELECT METHOD 2.
C 3 - METHOD 3 IS SIMILAR TO METHOD 1 BUT OMITS ALL CHECKS
C FOR ACCEPTABILITY. IT IS THE PREFERRED METHOD FOR
C VIRTUAL ORBITALS WHOSE OCCUPATION NUMBER IS SMALL.
C ACC - INITIAL ACCELERATING FACTOR, 0 .LE. ACC .LT. 1
C IND - INDICATOR SPECIFYING THE TYPE OF INITIAL ESTIMATES
C -1 - INPUT DATA TO BE READ FROM UNIT IUF
C 0 - SCREENED HYDROGENIC FUNCTIONS
C 1 - SAME AS RESULTS ALREADY IN MEMORY
C QC - NUMBER OF ELECTRONS I IN CONFIGURATION J. A VALUE
C OF -1 INDICATES THAT I'TH ORBITAL IS TO BE TREATED
C AS OCCUPIED IN THE J'TH CONFIGURATION FOR ORTHO-
C GONALITY PURPOSES WHEN ORTHO = .FALSE.
C
CARD 4. FOR EACH M, M = 1,NF (IF NF > 0) A CARD WITH
C A, WW, K, I1, J1, I2, J2 IN THE FORMAT(F12.8,A1,I1,1X,2I2,1X,2I2,1X)
C
C A - COEFFICIENT OF THE FK INTEGRAL
C W - THE CHARACTER 'F'
C K - THE VALUE OF K
C I1, J1 - THE I'TH RADIAL FUNCTION IN THE J'TH CONFIGURATION
C I2, J2
C
CARD 5. FOR EACH M, M = 1,NG (IF NG > 0) A CARD WITH
C B, W, K, I1, J1, I2, J2 IN THE FORMAT(F12.8,A1,I1,1X,2I2,1X,2I2)
C
C B - COEFFICIENT OF THE GK INTEGRAL
C W - THE CHARACTER 'G'
C K - THE VALUE OF K
C I1, J1 - THE I'TH RADIAL FUNCTION IN THE J'TH CONFIGURATION
C I2, J2
C
CARD 6. FOR EACH M, M = 1,NR (IF NR > 0) A CARD WITH
C D, W, K, I1, I2, J1, I3, I4, J2, IO, JO, P
C IN THE FORMAT(F12.8, 1X, I1, 1X, 3I2, 1X, 3I2, 1X,3(1X,I2))
C
C D - COEFFICIENT OF THE RK INTEGRAL
C W - THE CHARACTER 'R'
C K - THE VALUE OF K
C I1,I2,J1 - THE I'TH RADIAL FUNCTION IN THE J'TH CONFIGURATION
C I3,I4,J2
C IO,JO - RADIAL FUNCTIONS IN THE OVERLAP FACTOR <IO!JO>**P
C P - EXPONENT FOR THE OVERLAP FACTOR
C
CARD 7. FOR EACH M, M = 1,NL (IF NL > 0) A CARD WITH
C C, W, K, I1, J1, I2, J2, IO, JO, P
C IN THE FORMAT(F12.8, 2X, 2I2, 1X. 2I2, 1X, 3(1X,I2))
C
C C - COEFFICIENT OF THE L INTEGRAL
C W - THE CHARACTER 'L'
C I1, J1 - THE I'TH RADIAL FUNCTION IN THE J'TH CONFIGURATION
C I2, J2
C IO,JO - RADIAL FUNCTIONS IN THE OVERLAP FACTOR <IO!JO>**P
C P - EXPONENT FOR THE OVERLAP FACTOR.
C
CARD 8. FOR EACH I, I=1,NWF FOR WHICH IND(I) = -1, A RADIAL FUNCTION
C WHICH WAS OUTPUT DURING A PREVIOUS RUN MUST BE READ FROM
C UNIT IUF. IF THIS UNIT IS THE SAME AS THE SYSTEM INPUT UNIT
C FOR THE REST OF THE DATA, THE CARDS MUST BE INSERTED AT THIS
C POINT IN THE ORDER IN WHICH THEY WILL BE READ.
C
CARD 9. PRINT, NSCF, IC, ACFG, ID, CFGTOL, SCFTOL, LD
C IN THE FORMAT(L3, 2I3, F3.1, I3, 2F6.1, L3)
C
C PRINT - A LOGICAL VARIABLE: IF .TRUE. THE RADIAL FUNCTIONS WILL
C BE PRINTED ON UNIT 6, SEVEN FUNCTIONS PER PAGE. THE
C LOGICAL RECORD LENGTH FOR THIS UNIT MUST BE AT LEAST 130.
C NSCF - THE MAXIMUM NUMBER OF CYCLES FOR THE SCF PROCESS. IF
C LEFT BLANK OR ZERO, A DEFAULT OF 12 IS ASSUMED.
C IC - THE NUMBER OF ADDITIONAL IMPROVEMENTS OF RADIAL FUNCTIONS
C SELECTED ON THE BASIS OF GREATEST CHANGE, FOLLOWING A
C SWEEP THROUGH THE SYSTEM OF EQUATIONS BUT PRIOR TO
C REORTHOGONALIZATION, RECOMPUTATION OF OFF-DIAGONAL
C ENERGY PARAMETERS, AND REDIAGONALIZATION OF THE ENERGY
C MATRIX IN A MULTI-CONFIGURATION APPROXIMATION. IF = 0,
C THE DEFAULT VALUE OF 3 + NIT/4 IS ASSUMED IN A SINGLE
C CONFIGURATION APPROXIMATION. THEN A CYCLE CONSISTING ONLY
C OF SWEEPS THROUGH THE SYSTEM OF EQUATIONS REQUIRES THAT IC
C BE SET TO -1. WHEN NCFG > 1, THE ACTUAL VALUE IS USED.
C ACFG - ACCELERATING PARAMETER TO BE APPLIED TO THE WEIGHTS
C WT(I) AFTER AN ENERGY DIAGONALIZATION, 0 .LE. ACFG .LT. 1
C ID - IF NOT = 0, THE ENERGY MATRIX WILL BE COMPUTED BUT NOT
C DIAGONALIZED AND THE WEIGHTS LEFT UNCHANGED.
C CFGTOL- THE ENERGY CONVERGENCE TOLERANCE FOR A MULTI-CONFIGURATION
C CALCULATION. IF = 0, A DEFAULT OF 1.D-10 IS ASSUMED.
C SCFTOL- THE PARAMETER DEFINING THE SELF-CONSISTENCY TOLERANCE
C FOR RADIAL FUNCTIONS. IF = 0, A DEFAULT VALUE OF
C 1.D-7 IS ASSUMED.
C LD - A LOGICAL VARIABLE: IF .TRUE. THE ENERGY MATRIX IS DIAG-
C ONALIZED BEFORE THE SCF ITERATIONS BEGIN. THIS IS RECOM-
C MENDED ONLY WHEN THE RADIAL FUNCTIONS ARE KNOWN TO GOOD
C ACCURACY.
CARD 10. AT THIS POINT, IF IUH > 0, AN MXM HAMILTONIAN MATRIX AS OUTPUT
C FROM SOME PREVIOUS CALCULATION SHOULD BE INSERTED, WHERE
C M = (NCFG - NEW). DURING THE CALCULATION THIS PORTION
C OF THE INTERACTION MATRIX WILL BE LEFT UNCHANGED UNLESS
C DATA CARDS ARE INCLUDED FOR THIS PORTION OF THE MATRIX.
C THE TERMS REPRESENTED BY SUCH DATA CARDS WILL BE ADDED
C TO THE INITIAL MATRIX. WHEN IUH = 0, IT IS ASSUMED THAT
C THIS PORTION OF THE MATRIX IS ALREADY IN MEMORY.
C
CARD 10. END, NEXT, ATOM, ZZ, (ACC(I), I=1,NWF)
C IN THE FORMAT(A1, I2, A6, F6.0, 20F3.1)
C
C END - SPECIAL SYMBOL '*' DENOTING THE END OF A CASE. THE PROGRAM
C ASSUMES THAT THE NEXT CARD, IF ANY, IS A CARD OF TYPE 1.
C (THE ASTERISK IS IMPORTANT ONLY IN THE CASE OF A FAILURE
C TO CONVERGE WHEN THE PROGRAM ATTEMPTS TO FIND THE BEGINNING
C OF THE NEXT CASE.)
C NEXT - A VARIABLE DETERMINING SUBSEQUENT ACTION. IF NEXT = 0, THE
C CASE HAS COMPLETED SUCCESSFULLY AND THE FOLLOWING CARD, IF
C ANY, IS ASSUMED TO BE A CARD OF TYPE 1. IF NEXT = 1, THE
C PROGRAM WILL SCALE PRESENT RESULTS FOR ATOMIC NUMBER Z TO
C THOSE FOR ATOMIC NUMBER ZZ, USING A SCREENED HYDROGENIC
C SCALING LAW, AND REPEAT THE CALCULATION FOR THE NEW CASE.
C THE NEXT CARD WILL BE ASSUMED TO BE OF TYPE 9.
C ZZ - THE NEW ATOMIC NUMBER FOR THE NEXT CASE.
C ACC - THE NEW ACCELERATING PARAMETERS FOR THE NEXT CASE
C
C ALL CALCULATIONS FOR AN ISOELECTRONIC SEQUENCE ARE ASSUMED TO BE
C ONE CASE. SEVERAL EXAMPLES OF INPUT DATA NOW FOLLOW BUT NOTE
C THAT THE DATA CARDS MUST BE REMOVED BEFORE THE PROGRAM CAN BE
C COMPILED.
C
C ------------------------------------------------------------------
C 2-1 E X A M P L E 1
C ------------------------------------------------------------------
C
C THE FIRST EXAMPLE ILLUSTRATES HOW DATA CAN BE SET UP FOR AN
C ISOELECTRONIC SEQUENCE. NOTE THAT FOR EACH ATOM OR ION, THE
C SLATER INTEGRALS AND RADIAL FUNCTIONS ARE PRINTED. THE LATTER
C ALSO ARE OUTPUT, A FUNCTION AT A TIME, FOR FUTURE INPUT.
C
1 1 0 0
6 6 7 0
BE 1S 4. 180 2 2 1
1S2/2S2
1S 1 0 0.0 1 0 0 2
2S 2 0 3.0 1 0 0 2
T 12 4
1 B+ 5.0
T 12 4
1 C+2 6.0
T 12 4
1 O+4 8.0
T 12 4
1 NA+7 11.0
T 12 4
1 P+11 15.0
T 12 4
1 CA+16 20.0
T 12 4
1 MN+21 25.0
T 12 4
1 ZN+26 30.0
T 12 4
*
C
C THOUGH RESULTS ARE MEANINGFUL PHYSICALLY ONLY FOR INTEGRAL
C VALUES OF Z, A USEFUL FEATURE WHEN DOING CALCULATIONS FOR
C NEGATIVE IONS WHERE INITIAL ESTIMATES ARE NOT WELL KNOWN, IS
C THE ABILITY TO TURN OFF THE NUCLEAR CHARGE IN FRACTIONAL STEPS,
C THEREBY APPROACHING THE RELATIVELY UNSTABLE STATE WITH GOOD
C INITIAL ESTIMATES.
C
C ------------------------------------------------------------------
C 2-2 E X A M P L E 2
C ------------------------------------------------------------------
C
C THE SECOND EXAMPLE ILLUSTRATES A FIXED CORE CALCULATION. IN THE
C FIRST CASE, HARTREE-FOCK CALCULATIONS ARE PERFORMED FOR THE NA+
C CORE. THEN, WITH THESE FUNCTIONS FIXED AND STILL IN MEMORY, A
C FIXED CORE CALCULATION IS PERFORMED FIRST FOR 3S AND THEN FOR 3P.
C IN THIS CASE NO CALCULATIONS ARE SUMMARIZED NOR ARE ANY RESULTS
C PUNCHED. NOTE ALSO THE ZERO OCCUPATION NUMBER FOR 3S IN THE
C CALCULATION FOR 3P. THIS IS ALLOWED AS LONG AS THE ASSOCIATED
C RADIAL FUNCTION IS NOT INCLUDED IN THE SET OF FUNCTIONS BEING
C MADE SELF-CONSISTENT.
C
1 1 0 0
0 0 0 0
NA+ 1S 11. 190 3 3 1
2S2/2P6
1S 1 0 0.0 1 0 0 2
2S 2 0 2.0 1 0 0 2
2P 2 1 7.0 1 0 0 6
F 10 10
*
NA 2S 11. 190 4 1 1
2S2/2P6/3S
1S 1 0 0.0 1 0 1 2
2S 2 0 2.0 1 0 1 2
2P 2 1 7.0 1 0 1 6
3S 3 0 10.0 1 0 0 1
F 10 2
*
NA 2P 11. 190 5 1 1
2S2/2P6/3P
1S 1 0 0.0 1 0 1 2
2S 2 0 2.0 1 0 1 2
2P 2 1 7.0 1 0 1 6
3S 3 0 10.0 1 0 1 0
3P 3 1 10.0 1 0 0 1
F 10 2
*
C
C ------------------------------------------------------------------
C 2-3 E X A M P L E 3
C ------------------------------------------------------------------
C
C THIS EXAMPLE ILLUSTRATES THE USE OF MULTIPLE ORBITALS OF THE SAME
C TYPE. THIS FEATURE IS RESTRICTED TO CASES WHERE THE INTERACTION
C BETWEEN CONFIGURATIONS RESULTS IN TERMS WITH AT MOST ONE OVERLAP
C INTEGRAL. THESE ARE ALLOWED ONLY FOR RK AND L INTEGRALS AND SO
C OCCASSIONALLY FK OR GK INTEGRALS MUST BE RE-EXPRESSED AS RK
C INTEGRALS.
C
C
1 1 0 0
0 0 0 0
B 2P 5. 180 4 4 2 1 0 1
2S2/2P 1.0
2P*3 0.3
1S 1 0 0.0 1 0 0 2 2
2S 2 0 2.0 1 0 0 2
2P 2 1 4.0 1 0 0 1
2P* 2 1 3.0 3 0 0 0 3
0.240000000F2( 4 2, 4 2)
-0.94280904R1( 2 2 1, 4 4 2)< 3! 4> 1
F 8 1 0 0 .1E-7 .1E-5
*
C
C
C IN THE ABOVE EXAMPLE THE CFGTOL AND SCFTOL PARAMETERS HAVE BEEN
C RESET TO LARGER VALUES THAN THE DEFAULT VALUES AND NO FUNCTIONS
C ARE PUNCHED. THIS APPROACH IS USEFUL WHEN A TRIAL CALCULATION
C IS BEING PERFORMED.
C
C IN MORE COMPLEX CASES, MORE RESTRICTIONS ARE PRESENT. NORMALLY
C AN OFF-DIAGONAL ENERGY PARAMETER CAN BE COMPUTED DIRECTLY FROM
C A PAIR OF RADIAL EQUATIONS (SEE SEC. 7-3), BUT WHEN NON-ORTHOG-
C ONAL ORBITALS ARE PRESENT, INDICATED BY THE FACT THAT THEIR NL
C VALUES ARE THE SAME, IT MAY BE NECESSARY TO SOLVE A SYSTEM OF
C EQUATIONS FOR SETS OF PARAMETERS. WHEN ORTHO = .FALSE., THE
C PROGRAM ASSUMES OFF-DIAGONAL PARAMETERS CAN BE COMPUTED FROM A
C SINGLE PAIR OF EQUATIONS. WHEN ORTHO = .TRUE., THE PROGRAM
C CHECKS IF OTHER ORBITALS ARE PRESENT, NON-ORTHOGONAL TO EITHER
C ONE OF THE PAIR. AT MOST ONE SUCH ORBITAL FOR EACH MEMBER CAN
C BE ACCOMMODATED. THE IMPLICATION OF THESE RESTRICTIONS CAN BE
C DESCRIBED FOR THE MCHF CALCULATION --
C 2P6.(3S.3P + 3P'.3D )
C IF 2P IS DETERMINED VARIATIONALLY IT MUST BE ORTHOGONAL TO BOTH
C 3P AND 3P' WHICH ARE NON-ORTHOGONAL. HENCE A SYSTEM OF EQUATIONS
C MUST BE SOLVED AND WE NEED ORTHO = .TRUE. ON THE OTHER HAND, IF
C 2P WERE PART OF THE FIXED CORE, THE VALUE OF ORTHO WOULD NOT BE
C IMPORTANT. THE CALCULATION FOR
C 2P6.(3S.3P + 4S.4P + 3P'.3D + 4P'.4D )
C CAN ONLY BE PERFORMED CORRECTLY WITH 2P PART OF A FIXED CORE.
C ORTHO MUST BE .FALSE. OTHERWISE ALL FUNCTIONS WITH DIFFERENT NL
C VALUES WOULD BE MADE ORTHOGONAL. THE -1 OPTION MUST BE USED TO
C INCLUDE 4S,4P IN CONFIGURATION 1, 4P',4D IN CONFIGURATION 3
C FOR ORTHOGONALITY PURPOSES.
C
C
C ------------------------------------------------------------------
C 2-4 E X A M P L E 4
C -------------------------------------------------------------------
C
C THE EXAMPLE BELOW IS A TWO CONFIGURATION PROBLEM WHERE THE INTER-
C ACTION IS SMALL AND ALL FUNCTIONS MUST BE ORTHOGONAL. IN SUCH A
C CALCULATION IT IS IMPORTANT THAT THE FUNCTIONS DEFINING THE DOM-
C INANT CONFIGURATION BE FAIRLY ACCURATE AND SO A ONE CONFIGURATION
C CALCULATION IS USED FOR DETERMINING THE INITIAL ESTIMATES.
C
C THE INTERACTION IN THIS CASE IS SMALL RESULTING IN A SHALLOW
C ENERGY MINIMUM. THIS COMBINED WITH THE MANY ORTHOGONALITY CON-
C STRAINTS MAKES CONVERGENCE DIFFICULT. THE PROGRAM PRINTS OUT
C SEVERAL DIRE MESSAGES BUT THESE DISAPPEAR AS CONVERGENCE SETS IN.
C
1 1 0 0
0 0 0 0
HE 3S 2. 180 2 2 1 0 1 0
1S/2S
1S 1 0 0.0 1 0 0 1
2S 2 0 1.0 1 0 0 1
-.500000000G0( 1 1, 2 1)
F 12 -1
*
HE 3S 2. 180 4 2 2 0 2 2 0 T T
1S/2S 0.9965624
3S/4S -.0016
1S 1 0 0.0 1 0 1 1
2S 2 0 1.0 1 0 1 1
3S 3 0 -5.0 3 0 0 0 1
4S 4 0 -8.7 3 0 0 0 1
-.500000000G0( 1 1, 2 1)
-.5 G0( 3 2, 4 2)
2.0 R0( 1 2 1, 3 4 2)
-2.0 R0( 1 2 1, 4 3 2)
F 15 1 0 0 .1E-7 .1E-6
*
C
C NOTE THAT WHEN NO OVERLAP INTEGRALS ARE PRESENT, THEY CAN SIMPLY
C BE OMITTED. ALSO, NEGATIVE SCREENING PARAMETERS MAY BE USED TO
C TO CONTRACT THE VIRTUAL ORBITALS.
C
C
C ------------------------------------------------------------------
C 2-5 E X A M P L E 5
C -----------------------------------------------------------------
C
C ONE OF THE MORE DIFFICULT CASES, HISTORICALLY WAS 1S2S 1S OF HE.
C THIS IS AN EXAMPLE WHERE LARGE OFF-DIAGONAL ENERGY PARAMETERS
C OCCUR BECAUSE OF THE ORTHOGONALITY CONSTRAINT, THE EFFECT OF
C WHICH IS TO ROTATE THE USUAL 1S, 2S ORBITAL BASIS. ON THE
C OTHER HAND, THE 1S2S 3S CALCULATION IS A PARTICULARLY SIMPLE
C ONE. THE DATA BELOW DESCRIBES A CALCULATION IN WHICH THE
C OUTPUT FOR 3S DEFINES THE INITIAL ESTIMATES FOR THE 1S STATE.
C
C IN THIS CASE, THE 1S STATE IS NOT THE LOWEST STATE OF A GIVEN
C SYMMETRY AND THE HARTREE-FOCK APPROXIMATION IS NOT A PARTICULARLY
C GOOD ONE, THE TOTAL ENERGY BEING CONSIDERABLY BELOW THE OBSERVED.
C AN MCHF APPROXIMATION WHICH INCLUDES 1S2 1S YIELDS A
C MUCH BETTER RESULT. SUCH A CALCULATION FOLLOWS THE SINGLE
C CONFIGURATION CASE BELOW, WITH 3S FUNCTIONS AGAIN SERVING
C AS INITIAL ESTIMATES.
C
C
1 1 0 0
0 0 0 0
HE 3S 2. 180 2 2 1 0 1 0
1S/2S
1S 1 0 0.0 1 0 0 1
2S 2 0 1.0 1 0 0 1
-.500000000G0( 1 1, 2 1)
F 12 -1
*
HE 1S 2. 180 2 2 1 0 1 0
1S/2S
1S 1 0 0.0 1 0 1 1
2S 2 0 1.0 1 0 1 1
1.500000000G0( 1 1, 2 1)
F 12 -1
*
HE 3S 2. 180 2 2 1 0 1 0
1S/2S
1S 1 0 0.0 1 0 1 1
2S 2 0 1.0 1 0 1 1
-.500000000G0( 1 1, 2 1)
F 12 -1
*
HE 1S 2. 180 2 2 2 0 1 1 1 F T
1S/2S 0.9965624
1S2 .114
1S 1 0 0.0 1 0 1 1 2
2S 2 0 1.0 1 0 1 1
1.500000000G0( 1 1, 2 1)
2.82842712 R0( 1 2 1, 1 1 2)
-1.41421356L( 2 1, 1 2)
F 9 -1 .3 0 1.E-7 .1E-6 T
*
C
C
C IN THIS CASE, THE CONFIGURATIONS DIFFER BY EXACTLY ONE ELECTRON
C AND THE INTERACTION NOW ALSO INVOLVES AN L INTEGRAL. THE
C CALCULATIONS HAVE BEEN SET UP IN SUCH A WAY THAT FIRST 1S, THEN
C 2S ARE IMPROVED (IN THE MCHF CASE, 2S IS KEPT STRICTLY
C ORTHOGONAL TO 1S), AFTER WHICH THE FUNCTIONS ARE ROTATED TO
C SATISFY THE STATIONARY CONDITION, AND THE ENERGY MATRIX
C DIAGONALIZED. THE ROTATIONS SEEM TO INTRODUCE OSCILLATIONS
C IN THIS CASE AND SO ACFG WAS SET TO 0.3 IN ORDER TO DAMP OUT
C THESE OSCILLATIONS.
C
C
C ------------------------------------------------------------------
C 2-6 E X A M P L E 6
C -----------------------------------------------------------------
C
C THIS EXAMPLE SHOWS HOW A SIXTH CONFIGURATION CAN BE ADDED TO AN
C INTERACTION MATRIX, WHERE THE MATRIX WAS OUTPUT DURING A PREVIOUS
C CALCULATION. IT IS ASSUMED THAT RADIAL FUNCTIONS CAN BE READ
C FROM UNIT 8.
C
1 1 2 1
6 0 6 6
AR+2 3S 18. 200 6 0 6 2 3 6 0 F F 1
2P5/3P5 0.9920873
2S/3S 0.0364112
3S/3P5/3D(1P) 0.1149757
3S/3P5/3D(3P) -0.0348588
3S/3P5/3D(3P) -0.0015988
2S/3P4/3D -0.018
1S 1 0 0.0 1 0 -1 2 2 2 2 2 2
2S 2 0 3.0 1 0 -1 2 1 2 2 2 1
2P 2 1 7.0 1 0 -1 5 6 5 5 5 6
3S 3 0 10.0 1 0 -1 2 1 1 1 1 2
3P 3 1 10.0 1 0 -1 5 6 5 5 5 4
3D 3 2 12.0 3 0 -1 0 0 1 1 1 1
0.12000000F2( 506, 506)
-0.40000000F2( 506, 606)
0.53333333G1( 506, 606)
-0.08571429G3( 506, 606)
-0.10000000G2( 606, 206)
0.94280904R1( 6 3 6, 5 2 1)
-1.63299316R1( 4 6 6, 5 5 2)
-0.94280904R1( 4 3 6, 5 2 3)
1.41421356R0( 4 3 6, 2 5 3)
-0.81649658R0( 4 3 6, 2 5 4)
-2.30940108R0( 4 3 6, 2 5 5)
f 0 0 0 0 1.e-8 .1e-5
-516.5748431
-0.1543388 -513.0822716
-0.2328091 0.1413721 -514.7085685
0.1012540 -0.2399579 -0.2016125 -514.6445872
0.0062108 -0.0069378 -0.0322798 -0.0086021 -515.0489374
*
C
C NOTICE THAT THE DATA CARDS FOR THE ENERGY EXPRESSION NOW NEED
C REFER ONLY TO THE NEW CONFIGURATION, ALL OTHER ENTRIES OF THE
C INTERACTION MATRIX BEING KEPT FIXED. THIS CAN BE USED TO
C ADVANTAGE IN VERY LARGE PROBLEMS WHERE THE NUMBER OF FK, GK, OR
C RK INTEGRALS MIGHT EXCEED THE MAXIMUM ALLOWED BY THE DIMENSIONS
C OF THE ARRAYS. IN THIS WAY ONE CAN, IN EFFECT, EVALUATE PORTIONS
C OF THE MATRIX AT A TIME.
C
C
C ------------------------------------------------------------------
C 2-7 P E R F O R M A N C E T E S T D A T A
C ------------------------------------------------------------------
C
C
C AN INDICATION OF THE ACCURACY OF THE NUMERICAL PROCEDURE USED IN
C THIS PROGRAM IS GIVEN IN SEC. 6-9 AND ALSO COMPUT. PHYS. COMMUN.
C 4, 107 (1972). THESE RESULTS ALL PERTAIN TO THE HYDROGENIC
C PROBLEM. ADDITIONAL SOURCES OF ERRORS ENTER INTO THE GENERAL HF
C PROBLEM. THE RESULTS FROM THE TEST DATA PROVIDED BELOW, INDICATE
C THE LIMITS OF ACCURACY ON THE IBM 360/370 SERIES, AND ALSO PROVIDE
C INFORMATION ABOUT THE RATE OF CONVERGENCE FOR OUR METHODS. THE
C CASES SELECTED ARE TYPICAL TEST CASES OFTEN USED FOR EVALUATING
C A METHOD.
C
C
1 1 0 0
0 0 0 0
HE 1S 2. 150 1 1 1
1S2
1S 1 0 0.0 1 0 0 2
F 10 10 0 0 .1E-9 .1E-9
*
HE 1P 2. 180 2 2 1 0 1
1S/2P
1S 1 0 0.0 1 0 0 1
2P 2 1 1.0 1 0 0 1
0.50000000G1( 1 1, 2 1)
F 10 10 0 0.1E-12.1E-12
*
HE 3S 2. 180 2 2 1 0 1 0
1S/2S
1S 1 0 0.0 1 0 0 1
2S 2 0 1.0 1 0 0 1
-.500000000G0( 1 1, 2 1)
F 12 -1 0 0 .1E-9 .1E-9
*
HE 1S 2. 180 2 2 1 0 1 0
1S/2S
1S 1 0 0.0 1 0 1 1
2S 2 0 1.0 1 0 1 1
1.500000000G0( 1 1, 2 1)
F 12 -1 0 0 .1E-9 .1E-9
*
LI 4S 3. 200 3 3 1 0 3 0
1S/2S/3S
1S 1 0 0.0 1 0 0 1
2S 2 0 1.0 1 0 0 1
3S 3 0 2.0 1 0 0 1
-.500000000G0( 1 1, 2 1)
-.500000000G0( 1 1, 3 1)
-.500000000G0( 2 1, 3 1)
F 12 -1 0 0 .1E-9 .1E-9
*
W+64 1S 74. 200 3 3 1
1S2/2S2/2P6
1S 1 0 0.0 1 0 0 2
2S 2 0 2.0 1 0 0 2
2P 2 1 7.0 1 0 0 6
F 12 12 0 0.1E-11.1E-11
*
C
C
C THE FOLLOWING TABLE SUMMARIZES THE RESULTS.
C
C T A B L E OF R E S U L T S
C
C -E - TOTAL ENERGY IN ATOMIC UNITS
C PE/KE - RATIO OF POTENTIAL AND KINETIC ENERGY
C DPM - MAXIMUM DIFFERENCE BETWEEN THE FUNCTIONS OF THE
C PREVIOUS AND CURRENT SCF ITERATION
C NXCH - NUMBER OF EXCHANGE FUNCTIONS COMPUTED
C
C
C -E (PE/KE)+2 DPM NXCH
C 1. HE 1S2 1S 2.861679995 4.D-9 1.38D-10 13( 8)
C 2. HE 1S/2P 1P 2.122464215 2.D-9 2.33D-14 19( 7)
C 3. HE 1S/2S 3S 2.174250777 -12.D-9 1.49D-08 18(12)
C 4. HE 1S/2S 1S 2.169854456 -1.D-9 7.37D-09 24(20)
C 5. LI 1S/2S/3S 4S 5.204454130 -16.D-9 4.29D-09 27(21)
C 6. W+64 1S2/2S2/2P6 1S 10318.516188029 1.D-9 2.00D-12 24(12)
C
C THE METHODS USED IN THE PROGRAM RELY ON PROPERTIES OF THE HF
C EQUATIONS, PARTICULARLY WHEN SATISFYING THE ORTHOGONALITY REQUIR-
C EMENT. THE DISCRETIZED APPROXIMATION WILL HAVE THESE PROPERTIES
C ONLY TO LIMITED ACCURACY. FOR EXAMPLE, THE DISCRETIZED PROBLEM
C FOR HE 1S/2S 3S WITH OFF-DIAGONAL ENERGY PARAMTETERS EQUAL TO
C ZERO, HAS SOLUTIONS FOR WHICH <1S!2S> = 1.4D-8. THE ORTHOGON-
C IZATION PROCESS THEN LIMITS THE APPARENT SCF CONVERGENCE. HOW-
C EVER, A SMALLER DPM DOES NOT NECESSARILY QUARANTEE GREATER ACC-
C URACY. A BETTER INDICATION OF ACCURACY IS THE RATIO, (PE/KE)+2,
C WHICH FOR EXACT SOLUTIONS WOULD BE ZERO. IN MOST CASES THE
C MAGNITUDE OF THIS QUANTITY IS SEVERAL UNITS IN THE NINTH DECIMAL
C PLACE.
C
C AN INDICATION OF THE EFFECTIVENESS OF THE PROGRAM CAN BE OBTAINED
C FROM THE NUMBER OF TIMES AN EXCHANGE FUNCTION HAS TO BE COMPUTED.
C (IN LARGE ATOMS MOST OF THE COMPUTATION TIME IS DEVOTED TO THE
C EVALUATION OF EXCHANGE FUNCTIONS, AND SO THIS IS A REASONABLY
C CONVENIENT MACHINE INDEPENDENT MEASURE OF EFFICIENCY, WHERE THE
C MOST EFFICIENT PROGRAM ACHIEVES A GIVEN ACCURACY WITH THE
C LEAST AMOUNT OF COMPUTATION.) THE FIRST NUMBER UNDER THE NXCH
C COLUMN IS THE ACTUAL NUMBER OF CALLS TO XCH MADE BY THE PROGRAM,
C BUT BECAUSE OF THE HIGH ACCURACY REQUESTED, SEVERAL OF THESE DID
C NOT IMPROVE THE DEGREE OF SELF-CONSISTENCY. IN PARENTHESES IS
C THE NUMBER OF CALLS TO XCH REQUIRED FOR AN A POSTERIORI SELF-
C CONSISTENCY OF AT LEAST 8.5D-8, WHERE THE DEGREE OF SELF-CONSIST-
C ENCY IS NOW THE MAXIMUM DIFFERENCE BETWEEN THE CURRENT FUNCTION
C AND ITS NEXT SCF ITERATE. NOTICE THE VERY RAPID RATE OF CONVER-
C GENCE FOR HIGHLY IONIZED SYSTEMS. FOR W+64, AN AVERAGE OF 4.33
C IMPROVEMENTS PER FUNCTION PRODUCED AN A POSTERIORI DEGREE OF
C SELF-CONSISTENCY OF AT LEAST 6.6D-9 . A MORE DETAILED STUDY OF
C THE RATE OF CONVERGENCE FOR THE VARIOUS CASES WILL APPEAR IN THE
C JOURNAL OF COMPUTATIONAL PHYSICS, 1977.
C
▶EOF◀