- Email: [email protected]
Rapid evaluation of four-body cluster contributions
COMPUTER PHYSICS COMMUNICATIONS 2 (1971) 214-222. NORTH-HOLLAND PUBLISHING COMPANY
RAPID
EVALUATION
OF
FOUR-BODY
CLUSTER
CONTRIBUTIONS
G. P. MUELLER~ Code 7660, Naval Research Laboratory, Washington, D.C. 20390, USA Received 9 February 1971
PROGRAM
SUMMARY
Title of program (32 characters maximum): FOURBODY Catalogue number: ABIA Computer for which the program is designed and others upon which it is operable Computer: CDC 3800. Installation: Naval Research Laboratory, Washington, D.C. 20390 Operating system or monitor under which the program is executed: SCOPE Programming languages used: FORTRAN High speed store required: 3789 words. No. of bits in a word: 48 Is the program overlaid? No No. of magnetic tapes required: None What other peripherals are used?Card Reader; Line Printer No, of cards in combined program and test deck: 467 Card punching code: BCD Keywords descriptive of probleni and method of solution: Nuclear, Cluster Expansion, Numerical Integration, FourBody Contributions, Haselgrove, Alpha Matter, Monte—Carlo, Equidistributed Sequence. Nature of physical problem Within the framework of the cluster expansion for— malism we present a method for rapidly evaluating the four-body contributions to the ground state energy of an extended system of bosons, specifically, alpha matter. Since the numerical results for alpha matter have been published elsewhere [1], we present here a fictitious example for which the exact results are obtamable, Method of solution These calculations involve the numerical integration of six-dimensional integrals, which we evaluate using a method due to Haselgrove [2]. For the same accuracy, this method seems to require less running time and is easier to program than more conventional methods. ~ NRC—NRL Postdoctoral Resident Research Associate.
Restrictions on the complexity of the problem The program is restricted to the evaluation of six— dimensional integrals, but can easily be modified to evaluate an integral of almost any dimension. Typical running time Approximately 15 seconds are needed on the CDC 3800 for results that are accurate to 5% for the example given here. For a more realistic problem, perhaps twice that much running time would be needed. References [1] G. P. Mueller and J. W. Clark, Nuci. Phys. .4155 (1970) 561. [2] C. B. Haselgrove, Math. Computation 15 (1961) 323.
FOUR-BODY CLUSTER CONTRIBUTIONS
215
LONG WRITE-UP 1. JINTRODUCTION The primary purpose of this paper is to bring to the attention of physicists a method for rapidly evaluating many -dimensional numerical integrals. An example to illustrate its use is drawn from the cluster expansion approach for calculating the ground-state energy per particle of an infinitely extended system of bosons. Specifically, we present a program for evaluating the four-body contributions (involving six-dimensional integrals) that arise in this approach. This integration technique should be useful whenever a manydimensional integral must be evaluated numerically.
where p is the number density of particles and each of the four integrals runs over the volume, ~2, of the system. The effective potential V*(r) is given by
v*(r)
=
V(r)
2f(r)
-
~—[f(ryc.7
-
(17f(r))2]
,
(3)
where V(r) is a phenomenological two-body potential obtained from an analysis of the relevant two-body scattering data. Terms which are negligible in the limit N, ~2—, 00~ p constant, are not included in eq. (2). The factors i~(r)are given by ~(r) =f2(r) - 1 . (4) In order to estimate the contribution E 4 we
2. CLUSTER EXPANSION FOR A BOSON SYSTEM Consider an infinitely extended system of alpha particles interacting through a two-body, central potential. The alpha particles are assumed to be structureless. The ground-state energy per alpha particle can be written [1-3] E=E2+E3+E4+...
(1)
,
where E~is the contribution to the energy per particle arising from n-body configurations. The cluster expansion method is variational in nature; E depends functionally on a trial correlation factorf(r), which, in turn, depends on several variational parameters, some of which may be fixed by auxiliary conditions. The optimal values of the free variational parameters are those which minimize E, or, more particularly, which minimize a truncated form of the series (1) representing E. Several recent applications of the cluster expansion approach [1,4] use the sum of E2 and E3 to represent E. This truncation seems justified by the apparent rapid convergence of the series; at least, E31 was genera.lly much smaller than 1E21. These calculations are obviously improved by a reasonably accurate estimate of E4, which serves both as a correction to E2 + E3 and as a further check on the rate of convergence of the cluster series. For a boson system interacting through central, two-body forces the contribution E4 can be written 3/cl) fdr 2(r E4 = ~(p 1dr2dr3dr4V*(r12)f 12) +
~
2~(r34)77(r41)fl(r24) (2)
need to evaluate integrals of the form 1 fdr 9= ci1dr2dr3dr4g1(r12g2(r23) (5) x g3(r24)g4(r34)g5(r41)g6(r13) Provided that the integrand vanishes sufficiently rapidly as any of the four radial coordinates becomes large, 9 can be written as 00 2 00 dR 00 ~= 2(2~ j dRR 2R~ dR3R~
j
x
1
1
-1
-1
f
2~
f d~12I ~23 I dcog1(R1)g2(R2)
(6)
0
X g3(R3)g4(R23)g5(R13)g6(R12) where .R12 = R~3=
+
-
2R1R2x12
+
-
~1R3x23 (7)
R~1=R~+R~ - 2R3R1[x12x23 - {(i -x~2)(1_x~3)}1/2cosq.~]
3. THE HASELGROVE METHOD In order to evaluate these many-dimensional integrals we use a modified Monte-Carlo method based on the theory of equidistributed sequences (EDS). By way of introduction, consider a MonteCarlo approach to evaluating the integral ~ 1 1
I-I
dx1
fd~...I
dxdF(xl,x2...,xd).
(8)
216
G. P. MUELLER
The simplest Monte-Carlo method involves
of bounded variation, then the sequences S
estimating Iby calculating the sum
and S2(N) will satisfy the following relations for
M
1
‘MC
M1
=
(n) (n) F(X~ , X2
(n) ~.
. . ,
X~
1(N)
N sufficiently large:
)
(9)
,
Each component of the sampling point ,x~f))is a random number lying on the sufficiently interval well [0,1]. behaved One can function show {5] F and that for for large a enough M, I- Ii~vt-’
constant/M
~
1
19
/
(10)
.
The Monte-Carlo method can be improved by the adoption of some results from the theory of equidistributed defined as follows: sequences let z [5,61.An EDS can be 1, z2, ... be a sequence of points in the interval [a,d] and let {b,c] be a subinterval of [a,dj. Of the first N members of the sequence z1, z2, Z3, .., let MNbe the number that lies in the interval [b, C]. Then the sequence z1,z2,z3, . . . is said to be ~quidistributed if MN b -c lim—=---—=-. a -d
(11)
N—~°N
A particularly useful example of an EDS is the following: let (z) represent the non-integer portion of z, i.e., —
z
-
‘12’/
1 1 ‘
‘
equal to{z]z.isThen where the greatest if 0 is aninteger irrational less number, than or the sequence y(fl)
=
(nO)
n
,
=
is equi isLri u e ~
0, 1, 2, 3,
. . .,
(13)
-
~
2N+l S ‘N’ 2” ‘2 -
N
(15)
~
(N+1)
where the set of numbers {o~}, i = 1,2,.. .,d, is any set of irrational numbers that are linearly independent over the rationals. We note that 5 1(N)/(2N+1) is just the average value of the integrand over 2the is aNmodified samplingCesàro points mean and that of the 52(N)/(N÷1) sequence F((n91),.. ~, (nO,j)), n = 0, 1, 2
4. CHANGES OF VARIABLES There are many changes of variables which will bring an integral of the form (6) into the form (8). The primary consideration is to make the integrand as smooth as possible as a function of the new variables. Given the complexity of the integrands, we often can only guess at which new variables would be most useful. For many cluster calculations, V*(r) and q(r) are dominated by exponentials or gaussians, at least for large ticularly, the change r. Forofthe variables exponential case, par1°’ ‘ = 1 2 3 ‘ 3 ‘16 ~‘i‘ — ex ~-R~ E ~‘ ‘ ‘ ‘ ‘ ~‘ in (6) seems a useful one. For the angular integrations we use 2
2
X
~.
Let the sequences S1(.N) and S2(N) be defined
=
12 X12 ~----j—~—~—-~
X23 y5
=
X23
k--T--~—~—’-~‘6 ~
by
(17)
S1(N) =F(0 ‘ 0 I
...
‘
0)
. The purpose of the transformations yielding Y4
N +
F((nB
2 n=1
)
(no ‘
)
and Y5 is to emphasize the region of integration where the magnitudes R1 -Roj and R~-RQ~
(n9..~))
...
,
N
and =
~ S1(n) n=0
.
(14)
Haselgrove [7] demonstrates that if F(x1,x2,. . . ,xd) is continuous and possesses a derivative d’F(x1 , X2,. . . , Xd) ax1ax2.. .axd
,
are small, since the integrand contains factors which depend exponentially on these quantities. These transformations effectively smooth out the
‘
integrand. 5. PROGRAM DESCRIPTION
The detailed workings of the program are described by comment cards incorporated into it. In general terms, the program is split into two parts: a main program called FOURBODY and a
FOUR-BODY CLUSTER CONTRIBUTIONS
two-part subroutine called CECIL. FOURBODY consists of the integration routine; it generates the six-dimensional sampling points, calls on CECIL(2) for the values of the integrands at the sampling points and inserts these values into the sums S 1 and S2. CECIL(1) returns the overall 2/m for factor thatand instance) multiplies also the each valueintegral of the integrands (~p~h for the sampling point (0,0,. . ., 0). FOURBODY can be used with any V(r) and f(r) or, for that matter, for any six-dimensional integral, with only the following changes to be made by the user: the variables NIN and NOUT, which define the input and output units, must be set for the user’s machine and the FORMAT statements # 660 and # 680 must be adapted to the form of the user’s output. Much of the subroutine CECIL must be modified by the user to generate the particular integrands in which. he is interested, though the logical structure of the subroutine should remain intact,
6. NUMERICAL RESULTS
217
arising from the attractive and repulsive portions of the alpha-alpha potential.
7. GENERALIZATIONS AND IMPROVEMENTS We can suggest a number of modifications which could be incorporated into the program to enhance its usefulness. In addition, the program can be modified to evaluate almost any manydimensional integral. 1) If no more than a few tens of thousands of sampling points are evaluated and if the user’s computer carries at least eight significant figures, it is unnecessary to carry the sampling points in double precision and probably unnecessary to retain the running averages (described by comment cards in the program) to reduce round-off errors. 2) At present there is no provision for restarting the program at the N + 1 sampling point after it has run and printed out the results for N points. If the printout included the values of VONE(N) and VTWO(N), then the program could easily be modified to restart at the N + 1 sampling point. 3) With little effort, FOURBODY can be generalized to perform any n-dimensional integral, where n = 4, 5, 6 For n ‘ 3 it would be more efficient (greater accuracy for given running time) to use one of the more common methods, such as gaussian quadrature. We would like to emphasize again that the present program, though it may be useful to some researchers as it stands, is meant primanly as a guide to the use of the Haselgrove method.
The program included with this paper is used to evaluate a fictitious problem. We could have supplied the CECIL subroutine that was used for alpha matter [1], for which V(r) is the difference of two gaussians and f(r) is a third order polynomial in exp (- r), but it would be quite difficult to evaluate the integrals exactly for this case. Instead, we estimate the integrals that arise for the artificial case where 2 I (r)V(r) = ~(r) = exl (- r). (18) The four integrals listed on the computer output correspond to the four terms in brackets in eq. (2); they are labeled RING, DIAGONAL, OPENEND and CONNECTED on the printout. These names arise from a graphical representation of the boson cluster expansion. The running time on a CDC 3800 is approximately 2.7 seconds per one thousand sampling points. An examination of the output indicates that 2000 points is sufficient for 10% accuracy, 4000 for 5%, 10000 for 2% and 20000 for better than 1%. For more realistic cases the accuracy is somewhat reduced, perhaps halved. In the case of alpha matter, for instance, there is con-
[1] F. Iwamoto and M.Yamada, Progr. Theoret. Phys. 17 (1957) 543. [21 3. B. Aviles, Ann. Phys. 5 (1958) 251. [3] G. P. Mueller and 3. W. Clark, Nuci. Phys. A155 (1970) 561. [4] S.-O. Backman et al., Nucl. Phys. A130 (1969) 635. [5] P. J. Davis and P. Rabinowits, Numerical integration (Blaisdell, Waltham, 1967) P. 150. [6]P.3. Davis, Interpolation and approximation (Blaisdeli, New York, 1963) p. 354.
siderable cancelation between the contributions
[7]C. B. Haselgrove, Math. Computation 15 (1961) 323.
REFERENCES
218
G.P.MUELLER TEST TSST MO.
RUN
OUTPUT
c OF tOUaROOY W
POtr.TSs 04 25o120060
CO.
N
RI40
DIAGONAL
LI
RlhG
OI~GONAL
OF
IWEGRALS=
‘
OPENFNO
OPEWNll
CGkNECTEO CGhNEClEO
2TS.@D2
33.062 33.062
33.lw2 33.062
356.RCQ 4?9.566
44.019 49.18?
43.259 46.1hn
7.253 8.208
500
295.1_14 342.347
37.256 41.479
37.962 42.?7?
5.V06 6,986
FL0
3nl.lPV
5.405
3ir.s9n
35.943 38.403
37.396
75n
39.311
6,181
loon
291.5e7
36.976
37.159n
5.748
loon
3114.040
37.371
38.415
5.035
12% 125tl
240.941
36.514
3no.29n
37.004
36.1~7 37.Rnl
5.518 5.739
15nn
273.S13
34.948
35.112
5.188
15on
294.473
36.574
37.111
5.614
1750
2C7.7~0
34.5P4
1750
2PR.271
36.118
34.37n 36.974
5.031 5.473
2non 2060
271.131 2P4.4?2
IS.572 35.957
33.*9?. 36.n~r
5.300 5.432
22% 22%
26n.541 290.476
34.396 35.742
33.40= 35.560
5.093 5.363
2soo 75nn
2C0.910 Z?'I.PPO
34.6A3 35.422
33.765 35.179
5.3Sl 5.321
2750
2r4.407 273.9!2
34.689
33.47n
5.355
35.321
34.911
5.325
3001) znotl
2Lfl.944
34.306
33.417
5.2b3
27P.ZPO
35.164
34.67.I
5.324
3250 3250
255.331 269.941
33.1593 34.991
33.7a7 34.*m
5.134 5.302
3560
2&?.9?b
5.331
2**.5z2
34.150 34.R72
33.515
3500
34.3SL
5.295
3750 3750
2cU.915 267.612
33.826 34.766
33.nn3 34.mll
5.222 5.293
4000 4onn
24~.n24 zc*..*.
33.9134 34.634
33.337 34.002
5.156 5.279
4250 4256
2=9.w3 2h5,q72
34.038 34.557
33.249 33.9c(b
5.143 5.263
4500 l5nn
2ch.7lb
33.642
2rr.790
34.495
33.+m 33.9nR
5.125 5.251
4150 l15n
2C4.614 **1.7:7
33.663 34.41s
33.122 33.130
5.061 S.234
5000
2c5.241 ZrP.wL
33.563 34.330
33.3Rb
5.039
33.774
5.216
.%S.dfi?
2-n 750 500
77so
5000
(EXACT (EXACT
VALUES) VALUESI
FOUR-BODY CLUSTER CONTRIBUTIONS 9000 9000
257,711 255,792
33,444 33.565
33,405 33,549
5.043 5.081
9250 9250
257,043 235,882
23,496 33.560
33,356 33.350
5.045 5.078
9500 9500
256,970 255,979
33.53* 33,560
33.350 33,355
5.050 5.077
9759 9750
257,136 256,021.
33,459 33,536
33,465 33,356
5,021 5.075
10000 10000
257,116 256,061
33,467 33,550
35,370 33,364
5,032 5,072
19000 19000
253,480 253,048
23,i80 33.414
33,226 33,280
4.988 3,050
19250 19250
254,444 255,020
33,277 33,410
35.198 35,279
5,004 5,048
19500 19500
253,626 234,995
33,244 33,400
33,191 53,276
4,995 5,041
1.9750 19750
252,900 254,947
33,165 33,400
35,1.92 33,214
4,982 5,040
20000 20000
253,433 254,988
23,181 33.395
35,1.28 33,272
4,990 5,044
219
220
G. P. MUELLER PROGRAM LISTING
ABIAFOURBODY. RAPID EVALUATION OF FOUR—BODY CLUSTER CONTRIBUTIONS. 1 MUELLER, 0.?. %JOB,66H0106,2780PM,4 %FI’N,L,x PHOORAM FOURBODY
0 o 2 3 ~ 5 6
I
C C C C C C
C a C C C
COMMENTS CARDS WILL ALWAYS FOLLOW THE STATEMENTS P11 WHICH THEY REFER
11
C C
10 11
C
INITIALIZATION SECTION”*..
CALL CECIL(1) VARIABLE NAMES ARE GENERALLY SET OFF BY SINGLE DASHES
•
CECIL) 1) RETURNS TIE OVERALL FACTOR MULTIPLYING THE 1. TN INTEGRAL AND THE VALUE OF THE L TH INTEORAND EVALUATED FOE THU N—O SAMPLING POINT. THAT IS, FOR ALL OF THE INTEGRATION VARIABLES SET TO ZERO. THESE TWO QUANTITIES APPEAR IN VALUE(L+1O) AND VALUE)L+20), RESPECTIVELY, FOE L=l.NINT.
12 13 C C C DIMENSION THETA(6),T)6)
rr(6),vALUE)3o),sONU(9),STwO)9),AvE(g),
1CESAR)9),STJRE(9),RUNAV~9) CO~8)ON NINT,T,VALUE,NIN,NOUT TYPE DOUBLE TT,DEL,QUIFK,THETA ..—-IDENTIFICATION OF DIMENSIONED VARIABLES—.— —THETA— IS AN ARRAY OF SIX IRRATIONAL NUMBERS. —TO— CONTAINS THE CURRENT (6—DIM.) SAMPLING POINT. —THETA- AND -TO- ARE CARRIES IN DOUBLE PRECISION. —T— IS AN ARRAY CONTAINING THE CURRENT SAMPLING EDItOr IN SINGLE PRECISION. —DOSE— AND —STWO— ARE THE QUANTITIES DEFINED BY EQS. (14) IN THE LWP. —AND- IS THE AVERAGE VALUE OF THE INTEORAND OVER THE SAMPLING POINTS AND —CESAR— II THE CESARO MEAN OF THAT SEQUENCE. THEY ARE DERIVES FROM —lONE- AND —ITWO—, RESPECTIVELY. —VALUE-, -HUNAV— ASS —STORE— ARE DESCRIBED LATER.
C C • •
114 15 16 1~ 19 20 21 22 23 ~ 25 26 27 211 29 30 31 32 33 35
50 300 L=1,NINT RUNAV(L)—O.G DONE(L)—VALUE(L+20) 300 STWO(L).SONE)L) C C
c C C C C C
Soc
C C
35 36
37 **~I00INPUTSECTIO
35 39 SO
41 1’ c NIH.60 NOUT=61 o C —BIN- AND —HOST— LABEL THE INPUT AND OUTPUT UNITS. 1420 C C READ )HIN,160) C C THIS READS ONE CARD. COLUMNS 2—52 OF WHICH CAN RE USED AS A C HEADER FOR THE OUTPUT. 55 C C READ )NIN,17o) THETA C C HEADS DATA CARDS TWO AND THREE, WHICH CONTAIN SIX IRRATIONAL C NUMBERS IN DOUBLE PRECISION. HAIELOROVE GIVES SETS OF C IRRATIONAL NUMBERS FOR INTEGRALS OF VARIOUS DIMENSIONS, BUT C THEY POSSESS ONLY EIGHT SIGNIFICANT FIGURES. INSTEAD, WE USE C THE NON—INTEGER PORTIONS OF THE SQUARE ROOTS OF U. 3. s. C 7, 11 AND 17
C
66
O 80 CONTINUE READ )NIN,200)
MAX,JAX,NINT,NONE
C -——IDENTIFICATION OF INPUT VARIABLES—--MAX— IS THE ~IYFTALNUMBER OF SAMPLING POINTS TO BE USED. THE PROGRAM PRINTS OUT ITS CURRENT ESTIMATES, THAT IS, THE CURRENT VALUES OF —AND— AND —CESAR-, EVERY —SAX— POINTS. -HINT— IS THE NUMBER OF INTEGRALS TO NE EVALUATED. IF —NONE— IS ZERO OR H NEGATIVE INTEGER, THE PROGRAM TRANSFERS TO END AFTER SAMPLING —MAX— POINTS. IF —NONE- IS A POSITIVE INTEGER, THE PROGRAM TRANSFERS TO 80 AFTER SAMPLING —MAX— POINTS AND ATTEMPTS TO READ IN NEW DATA. AFTER EVERY MULTIPLE OF —SAX— POINTS, —RUNAV— IS SET ESUAL TO —CESAR— FOR 1.—i ,NINT. THEN —SOME— IS MODIFIED SO THAT IT ACCUMULATES THE THE DIFFERENCES BETWEEN THE CURRENT SAMPLING POINT AND —RUNAV—. THE PURPOSE OF THIS MANEUVER IS TO REDUCE ROUROOFF ERRORS. SEE THE LWP FOR A DISCUSSION. -STORE- IS USED IN THE EVALUATION OF -RUI6AV—.
C
42 43 44 45 46 4~ ~9 50 51 52 53 54 1115
57 511 59 60 61 62 63 64 65 6~ 6 69 70 71 72 73 74
77 7~ 14 So
81 82
83
1514 85 96
87
88
89 WHITE WHITE WHITE
NOIJT,160 NO1JT,220 NOUT,240
JAX,MAX,NINT THETA
THE THREE WRITE STATEMENTS IMMEDIATELY ABOVE WHITE OUT THE HEADER CARD, THE ARRAY —THETA— AND THE VARIABLES —MAX-, -SAX— AND —HINT—.
90 91 92 93 94 95 96
97
95
C C C C
C C C C C C C 460 C C C C
480 490 500 C C C C C
131 137 133 134 1 35 LOOP 300 INSERTS THE N=C CONTRIBUTIONS INTO THE SUNS —SONS- AND 1 36 —GTWO— AND SETS THE INITIAL VALUE SF THE RUNNING AVERAGE 22 13 zERO. 13 139 140 1141 N=O 1142 143 —N- IS THE INDEX THAT LABELS THE CURRENT SAMPLING POINT. 144 145 1146 147 1IAX—JAX 14i 149 AF’I’EE EACH SUCCESSIVE SET OF —SAX- POINTS HAVE BEEN EVALUATED, 150 THE PROORAM EXITS FORM ITS NORMAL ROUTINE, PRINTS OUT ITS 1 51 CURRENT ESTIMATE OF THE VALUE OF THE INTEGRALS AND RECALCULATES 1 52 THE RUNNING AVERAGE, —IRIJNAV—. —HAS— IS THE INDEX WHICH THE 153 PROGRAM TESTS TO SEE IF THE EVALUATION OF ANOTHER -JAX— POINTS 1 58 IS COMPLETE. 155 156 15 15 DO 420 K-i ,6 159 ‘117(15)—SOD 160 161 LOOP 420 SETS THE INITIAL VALUE OF THE SAMPLING POINT AT ZPEC. 162 163 1614 165 166 16 16 169 EVALUATIOH SECTIO 170 171 172 173 N-N*1 174 175 INCREMENTS THE SAMPLING POINT NUMBER, -N—. 176 177 178 179 OG500K—1 6 180 T’r)K)—Tr(lc)*THETA(iC) 181 QUIRE—TT(i5)_1.OS 182 IF (QUIRK) 490,490,580 183 l’r(K)—QUIHI( 184 T)K)’I’T)K) 185 CONTINUE 186
18
LOOT 500 CALCULATES THE COMPONENTS OF THE NEW SAMPLING POINT AND STORES THEM IH -IT—)IIJUBLE PRECISION) AND IN -T- (ROUNDED OFF TO SINGLE PRECISION).
FORMAT 15,51 H) 234567890123456789012345678901 23456789012345670901) FIRMAT 3150) FORMAT 110,315) FORMAT 18H0tHJ. OF POINTS. O(.I4,1H),I5, 15X,17HNO. OF INTEORALS—,12) 240 FORMAT (26NOSE,r OF IRRATIONAL NUMBERI/(3D30.17)) THE FIVE STATEMENTS ABOVE ARE THE FORMAT STATEMENTS FOR THE INPUT SECTION OF THE PROGRAM.
99 100 101 102 103 104 105 106 07 1001 109 110
18 1 190 191 IS’
193 1914 195 196 19 19 C 199 C ~oc SO 520 L—1,HINT 201 SONE(L SDRE)L *VALUE(L)_NUNAV(L) 202 520 STWO)L —STWO(L *SONE(L( 203 1 ‘04 C LOOP 520 ADDS THE CONTRIBUTION OF THE NEW SAMPLING POINT TO THE 205 C SUMS —GONE— AND —STWO—. NOTE THAT THE RUNNING AVERAGE IS 206 C SUBTRACTED FROM THE VALUE OF THE CURRENT SAMPLING POINT. ‘0 C
CALL 011114(2)
C C C
CECIL) 2) CALCULATES THE VALUE OF THE INTEGEANI AT TIlE HEW SAMPLING POINT AND RETURNS IT IN VALUE(L), FOR 1.—I ,NINT.
a C
160 170 200 220 O C C
C C C C C C C
75 76
iii 112 113 115 1 5 116 11 11 119 120 121 122 123 125 125 126 127 128
2111 IF )N—NAx) 460,600,600 TESTS WHETHER THE INTEORAND HAS BEER EVALUATED FOR A COMPLETE SET OF —SAX— POINTS. IF NOT, THE PROGRAM RETURNS TO 460 AND EVALUATES ANOTHER POINT. IF SO, THE PROGRAM EXITS TO 600.
C C C C , , ,
C r ~ C
Note: Columns 73 to 76 of each card image are omitted to save space,
~12 213 214 ‘15 16 21 1 219 ‘20 221 222
FOUR-BODY CLUSTER CONTRIBUTIONS ~
~
C
2 .
600 TRIO~2*N+1 TRIT...(N+i) *(~+) 1’WO(L)/TRIT STORE L)-S AVE(L (IONE(L)(TRID+O.5*RIJRAV(L)).VALUE(L+1O) 620 CESAR L)—(STORE(L)*0.5.RUNAV(L) ).vALUE(L+iO) C C C C
227
LOOP 620 CALCULATES THE ESTIMATES —AVE— AND —CESAF— FOR THE VARIOUS INTEGRALS. WRITE (HOST,660)
N, (AVE(L) ,L-i , NIMI)
WRITE (NOIFT,680) N,(CESAR)L),L1,NINT) a C
230 ‘31 232 233 234 235 ~
THESE WRITE OUT THE CURRENT ESTIMATES -AND- AND -CESAN-.
2~lO 241 242 243
c
C C
660 FORMAT (1H0,I9,5Fi5.3)
‘45 246
680 FORMAT (Ii0,5F15.3)
247
THESE TWO FORMAT STATEMENTS FOR ‘THE PRECEDING WRITE STATEMENTS MUST BE MODIFIED IF MORE THAN FIVE (5) INTEGRALS ARE EVALUATED OR IF A DIFFERENT OUTPUT FORMAT IS DESIRED.
C C C C
249
550 251 252
C C C
255 256 257
LOOP 760 CALCULATES THE NEW VALUES OF THE RUNNING MAKES TNE APPROPIATE CHANGES IN -lONE- AND -ITWO-.
~
C C C
a
C C C C C I
~TEMBFO a C
(
C C COMMON NINT,T,VALUE,NIN NOUT DIMENSION VALUE) 30) ,T( 65 ,R(6),F.(6) C C C C C C C C
———IDENTIFICATION OF SUBSCRIPTED VARIABLES—-— -VALUE- AND -T- ARE DEFINED IN THE MAIN PROGRAM. THE COMPONENTS OF —F- ARK THE RADIAL SEPARATIONS FOR WHICH TEE INTEORANDS ARE EVALUATED. THE COMPONENTS OF —R— ARE ORTAINED FROM -T- BY CHANGES OF VARIABLE. IN THIS CASE, H(K) — —FOWER~NATURALLOU(T(K)), E-1,2 AND 3. —POWER— IS A VARIABLE THAT IS READ IN BELOW. R)14) R)5) AND R)6) ARE OBTAINED FROM R)i), F(S) AND R)3) AND T(45, T)5) AND T(6) IN THE MANNER DESCRIBED IN THE LWP.
C C
)*(pOwEp*powER.pr,JER.3. *3. .TPI).(O.
C
NOTE——GENERALLY THERE IS A DIFFERENT FACTOR FOR EACH OF THE IPEDRALS.
350
(
(
9‘ C C
260 C C C C
*
358
359
*8*t
4
269
36~
VALUEçL+20~—0.5.10. 0)
365 366 367
VALUE~L+10-FACT
RETURN •
*
*#******4**•*** *
270 271 272 273 174
*8
4
*
368
369 370 371 372
373
374
375 376 377
*
*8*. 378 379
****8*******t*****•******** **~***~********~*
****..*
380
381 382 383 384
*****SECTION TWO*•***
275 276 277
*
*
386 387
A**.n*
388 520 R(I ——POWE~L4DGF(T(I)> I~R 1
391 352
281
Y—R 2
285 286 28
QYZ—SQ,1T5t~T( 5 >+. 25> AYZ——1 . 5+QYZ ATh—m*Afl.5927( Ii .—AXY.AXY)*(
396 3 3 399
1 .—AYZ*AYZ> ).C05(FIWT(6>)
XS~X*X YS—Y*Y ZSZ*Z N 4 —SQRT XS+YS—2.*X~Y.AXY R 5 —SORT YS+ZS—2.*Y.Z*AYZ R 6 SQJ1T Z5+XS—2. .Z~j.AZE
289 ‘go ‘91 292 293 ~
C
295 ‘96
C C
‘TI 29~ 299 300
C
4 401 402 403 404 405
STATEMENTS 500 THROOUN THE LAST CALCULATE THE RADIAL SEPARATIONS IN TERMS OF THE COMPONENTS OF -T- AS DESCRIBED ABOVE AND IN ‘THE LWP.
406 407 408 409
**.****.**.*44
410 1412 413
C C
411 .
V—v.V/(T(1 )wr(2)..T(3>.qxy.QYZ) C C C
414 415
—V. IS ESSENTIALLY ‘TIlE JAC1BIAN REPRESENTING THE ABOVE CHANGE OF VARIABLES. IT IS USED TM t.~P 720.
4*8*8*8*88*
fh
C ~
00 660
416 41 81 819
an
—
(R(I~—36~)620,640,640
IF
C
311 312
GO TO (101,500),M0DE TRANSFERS TO 100 OR 500 ACCORDING TO WHETHER CECIL) 1) OR CECIL(2) IS CALLED. • •
314 315 316 31 ~
C C C
),
325 326 327 32~ 330 331 332 ~
422
423 424 425
I.~F 660 CALCULATES EEP)—P(K>> POP 11—1,6. EXPONENT UN~YL0W IS INCLUDED. *
a.
A26
42 42
A CHECK TO PREVENT *
n........n...
Y—E(1 .E(2>.E(6) VALUE I .T.E( 5 VALUE I? —VALUE i).E
322 323 324
100 READ (HIM, 120) POWER WRITE (ttou’T,14o) POWER WRITE )NOUT,160) 120 FORMAT P5.2> 140 FORMAT 6HIR — —.F3.i,23H • 1.00) Y(I) 1—1,2,3) 160 FORMAT 151 18X,1HN,11X,INNINO.7X,8HDIAGOHAL,8X,7HOPENEND,6X,9HCONNECT!D,
k21
620 E(I>—EICP(—R(I)> 00 TO 660 6140 E(I)—C.0 660 CONTINUE
C SECTION ONE
360 361
00 260 L—i HINT
LOUP 260 PLACES THE OVERALL FACTOR -FACT— IN VALUE(I,f10) AND THE VALUE OF THE INTFINRAND FOR THE N—O SAMPLING POINT IN VALUE(L+20~, FOR L—i,NINT. THE FACTOR OF 0.5 IN SHE VALUE~L+20, STATEMENT ARISES BECAUSE THE N—C SAMPLING POINT ONLY RECEIVES HALF THE WEIGHT OF THE OTHERS. IN THE CASE PRESENTED HERE, THE VALUES OF THE INTEDRANDS FOR THE N—O SAMPLING POINT ARE ZERO.
•
302 303 304 305 306 307 30~
5)
IT IS MADE UP OF POUR FACTORS. (2.> ARISES BECAUSE EACH 351 SAMPLING POINT, OI~ TEAM THE N-C POIN1I’, HAS A WEIGHT 2 IN ‘THE 352 -SOME- SUM. (2. .‘I’PIrrPI) ARISES FROM VARIOUS ANGULAR INOTMS353 RATIONS (SEE EQ. 6) IN THE LW). THE NUT FACTOR ARISES FRUM 354 THE CHANGES OF VARIABLE THAT ABE DESCRIBED ABOVE MID IN THE 355 LW. TEE LAST FACTOR IS ONE—HALF RHO CUBED, WHICH APPEARS IN 356 EQUATION 2>, POE INSTANCE. WE HAVE HERE SET RHO—i .00. 357
266 267 204
301
345 346
—FACT— IS THE OVERALL FACTOR THAT MULTIPLIES EACH INFE~RAL.
,~
SUBROUTINE CECIL(MODE
34~
t4.***..
C C C C C C C
C C C C
•••
3~O 341 342
C
~
C C o C o
338
344
C C
STATEMENTS 800 THROUGH 900 INCREMENT —NAB— BY —SAX— AND RETURN It THE EVALUATION SECTION IN ORDER TO EVALUATE —SAX— MORE POINTS. IF —MAX— POINTS NAVE ALREADY BEEN EVALUATED, THE PROGRAM CHECKS TO SEE IF A NEW SET OF INTEGRALS IS TO BE CALCULATED. IF NOT, IT TRANSFERS TO END. IF SO, IT TRANSFERS
*
;~1~,~527 FACr—(2. ).(2. ~rpIwrpI
262 263 264 265
800 NAX—NAX+JAX IF (NAX—MAX) 460,460,920 820 IF (NONE) 900 900,840 640 WRITE (~ioirr810) 860 FORMAT (IHi ) GO TO 80 900 CONTINUE
337
POWER, WRITES IT OUT AND WRITES OUT THE LABELS POP ‘THE OUTP1.rT. •
AVERAGES AND 260 261
C
41 1TNG,~T IIDIAOOHAL,8X, nJrw.nW,6X,5a.tau.SJr5u,
517X,8H 255.802,8X,7H 33.062,8X,7H 33.062,24X,14H(UAC’r VALU!5>/ 617x,8H 255.802,8X,7N 33.062,8X,7H 33.062,24X,IIIH(EXAC’r VALUES)/) THE ABOVE GROUP OF STAT~IENTSREADS IN THE VARIABLE
,~
00 760 L-1,NINT RUHAV(L)—RUNAV(L)+2.*STORE(L) SOME(L (—loNE) L)-TRIO.E’TORE(L) 760 STWO(L (—0.0
~i~ifl
221
3> VALUE 3 Y*E(3 *E(4 VALUE 14 —VALUE 3>*E 5) C C C C CO C C C ~
‘
THE ABOVE FIVE STA72MENTS STORE THE VALUES OF THE POUR INTMSRAHDS IN VALUE(L), L—l 4. THERE REPRESEMr ‘tEE POUR TERMS IN THE INTHORAND OF EQ. 25 IN THE LW. •
1
*8*88*
~
*
*33 838 435 836 837 438
839
it4o
441
442 44
222
G. P. MUELLER
C
C
EACH INTBDRAND By THE JACOBIAN REPRESENTING 447 THE CHANGE OF VARIABLES DESCRIBED ABOVE. 448 SOUP 720 MULTIPLIES
449
*
a a
450 451 ‘452
RETURN END a *
* **
1453 454 455
1
1456 ‘457
1 1
SCOPE %LOAD %RUN,4,6000 TEST C OF FOURBODY .414213562373095049 .732050807568877294 .236077977499789696
458 459 460 461 462 463
.645751311064S90591 .31662479035S399849 .1231 05625617660550
1464
20000 3.0 •
250
4
0
465
466
467
1 1 1 1
1
RAPID
EVALUATION
OF
FOUR-BODY
CLUSTER
CONTRIBUTIONS
G. P. MUELLER~ Code 7660, Naval Research Laboratory, Washington, D.C. 20390, USA Received 9 February 1971
PROGRAM
SUMMARY
Title of program (32 characters maximum): FOURBODY Catalogue number: ABIA Computer for which the program is designed and others upon which it is operable Computer: CDC 3800. Installation: Naval Research Laboratory, Washington, D.C. 20390 Operating system or monitor under which the program is executed: SCOPE Programming languages used: FORTRAN High speed store required: 3789 words. No. of bits in a word: 48 Is the program overlaid? No No. of magnetic tapes required: None What other peripherals are used?Card Reader; Line Printer No, of cards in combined program and test deck: 467 Card punching code: BCD Keywords descriptive of probleni and method of solution: Nuclear, Cluster Expansion, Numerical Integration, FourBody Contributions, Haselgrove, Alpha Matter, Monte—Carlo, Equidistributed Sequence. Nature of physical problem Within the framework of the cluster expansion for— malism we present a method for rapidly evaluating the four-body contributions to the ground state energy of an extended system of bosons, specifically, alpha matter. Since the numerical results for alpha matter have been published elsewhere [1], we present here a fictitious example for which the exact results are obtamable, Method of solution These calculations involve the numerical integration of six-dimensional integrals, which we evaluate using a method due to Haselgrove [2]. For the same accuracy, this method seems to require less running time and is easier to program than more conventional methods. ~ NRC—NRL Postdoctoral Resident Research Associate.
Restrictions on the complexity of the problem The program is restricted to the evaluation of six— dimensional integrals, but can easily be modified to evaluate an integral of almost any dimension. Typical running time Approximately 15 seconds are needed on the CDC 3800 for results that are accurate to 5% for the example given here. For a more realistic problem, perhaps twice that much running time would be needed. References [1] G. P. Mueller and J. W. Clark, Nuci. Phys. .4155 (1970) 561. [2] C. B. Haselgrove, Math. Computation 15 (1961) 323.
FOUR-BODY CLUSTER CONTRIBUTIONS
215
LONG WRITE-UP 1. JINTRODUCTION The primary purpose of this paper is to bring to the attention of physicists a method for rapidly evaluating many -dimensional numerical integrals. An example to illustrate its use is drawn from the cluster expansion approach for calculating the ground-state energy per particle of an infinitely extended system of bosons. Specifically, we present a program for evaluating the four-body contributions (involving six-dimensional integrals) that arise in this approach. This integration technique should be useful whenever a manydimensional integral must be evaluated numerically.
where p is the number density of particles and each of the four integrals runs over the volume, ~2, of the system. The effective potential V*(r) is given by
v*(r)
=
V(r)
2f(r)
-
~—[f(ryc.7
-
(17f(r))2]
,
(3)
where V(r) is a phenomenological two-body potential obtained from an analysis of the relevant two-body scattering data. Terms which are negligible in the limit N, ~2—, 00~ p constant, are not included in eq. (2). The factors i~(r)are given by ~(r) =f2(r) - 1 . (4) In order to estimate the contribution E 4 we
2. CLUSTER EXPANSION FOR A BOSON SYSTEM Consider an infinitely extended system of alpha particles interacting through a two-body, central potential. The alpha particles are assumed to be structureless. The ground-state energy per alpha particle can be written [1-3] E=E2+E3+E4+...
(1)
,
where E~is the contribution to the energy per particle arising from n-body configurations. The cluster expansion method is variational in nature; E depends functionally on a trial correlation factorf(r), which, in turn, depends on several variational parameters, some of which may be fixed by auxiliary conditions. The optimal values of the free variational parameters are those which minimize E, or, more particularly, which minimize a truncated form of the series (1) representing E. Several recent applications of the cluster expansion approach [1,4] use the sum of E2 and E3 to represent E. This truncation seems justified by the apparent rapid convergence of the series; at least, E31 was genera.lly much smaller than 1E21. These calculations are obviously improved by a reasonably accurate estimate of E4, which serves both as a correction to E2 + E3 and as a further check on the rate of convergence of the cluster series. For a boson system interacting through central, two-body forces the contribution E4 can be written 3/cl) fdr 2(r E4 = ~(p 1dr2dr3dr4V*(r12)f 12) +
~
2~(r34)77(r41)fl(r24) (2)
need to evaluate integrals of the form 1 fdr 9= ci1dr2dr3dr4g1(r12g2(r23) (5) x g3(r24)g4(r34)g5(r41)g6(r13) Provided that the integrand vanishes sufficiently rapidly as any of the four radial coordinates becomes large, 9 can be written as 00 2 00 dR 00 ~= 2(2~ j dRR 2R~ dR3R~
j
x
1
1
-1
-1
f
2~
f d~12I ~23 I dcog1(R1)g2(R2)
(6)
0
X g3(R3)g4(R23)g5(R13)g6(R12) where .R12 = R~3=
+
-
2R1R2x12
+
-
~1R3x23 (7)
R~1=R~+R~ - 2R3R1[x12x23 - {(i -x~2)(1_x~3)}1/2cosq.~]
3. THE HASELGROVE METHOD In order to evaluate these many-dimensional integrals we use a modified Monte-Carlo method based on the theory of equidistributed sequences (EDS). By way of introduction, consider a MonteCarlo approach to evaluating the integral ~ 1 1
I-I
dx1
fd~...I
dxdF(xl,x2...,xd).
(8)
216
G. P. MUELLER
The simplest Monte-Carlo method involves
of bounded variation, then the sequences S
estimating Iby calculating the sum
and S2(N) will satisfy the following relations for
M
1
‘MC
M1
=
(n) (n) F(X~ , X2
(n) ~.
. . ,
X~
1(N)
N sufficiently large:
)
(9)
,
Each component of the sampling point ,x~f))is a random number lying on the sufficiently interval well [0,1]. behaved One can function show {5] F and that for for large a enough M, I- Ii~vt-’
constant/M
~
1
19
/
(10)
.
The Monte-Carlo method can be improved by the adoption of some results from the theory of equidistributed defined as follows: sequences let z [5,61.An EDS can be 1, z2, ... be a sequence of points in the interval [a,d] and let {b,c] be a subinterval of [a,dj. Of the first N members of the sequence z1, z2, Z3, .., let MNbe the number that lies in the interval [b, C]. Then the sequence z1,z2,z3, . . . is said to be ~quidistributed if MN b -c lim—=---—=-. a -d
(11)
N—~°N
A particularly useful example of an EDS is the following: let (z) represent the non-integer portion of z, i.e., —
z
-
‘12’/
1 1 ‘
‘
equal to{z]z.isThen where the greatest if 0 is aninteger irrational less number, than or the sequence y(fl)
=
(nO)
n
,
=
is equi isLri u e ~
0, 1, 2, 3,
. . .,
(13)
-
~
2N+l S ‘N’ 2” ‘2 -
N
(15)
~
(N+1)
where the set of numbers {o~}, i = 1,2,.. .,d, is any set of irrational numbers that are linearly independent over the rationals. We note that 5 1(N)/(2N+1) is just the average value of the integrand over 2the is aNmodified samplingCesàro points mean and that of the 52(N)/(N÷1) sequence F((n91),.. ~, (nO,j)), n = 0, 1, 2
4. CHANGES OF VARIABLES There are many changes of variables which will bring an integral of the form (6) into the form (8). The primary consideration is to make the integrand as smooth as possible as a function of the new variables. Given the complexity of the integrands, we often can only guess at which new variables would be most useful. For many cluster calculations, V*(r) and q(r) are dominated by exponentials or gaussians, at least for large ticularly, the change r. Forofthe variables exponential case, par1°’ ‘ = 1 2 3 ‘ 3 ‘16 ~‘i‘ — ex ~-R~ E ~‘ ‘ ‘ ‘ ‘ ~‘ in (6) seems a useful one. For the angular integrations we use 2
2
X
~.
Let the sequences S1(.N) and S2(N) be defined
=
12 X12 ~----j—~—~—-~
X23 y5
=
X23
k--T--~—~—’-~‘6 ~
by
(17)
S1(N) =F(0 ‘ 0 I
...
‘
0)
. The purpose of the transformations yielding Y4
N +
F((nB
2 n=1
)
(no ‘
)
and Y5 is to emphasize the region of integration where the magnitudes R1 -Roj and R~-RQ~
(n9..~))
...
,
N
and =
~ S1(n) n=0
.
(14)
Haselgrove [7] demonstrates that if F(x1,x2,. . . ,xd) is continuous and possesses a derivative d’F(x1 , X2,. . . , Xd) ax1ax2.. .axd
,
are small, since the integrand contains factors which depend exponentially on these quantities. These transformations effectively smooth out the
‘
integrand. 5. PROGRAM DESCRIPTION
The detailed workings of the program are described by comment cards incorporated into it. In general terms, the program is split into two parts: a main program called FOURBODY and a
FOUR-BODY CLUSTER CONTRIBUTIONS
two-part subroutine called CECIL. FOURBODY consists of the integration routine; it generates the six-dimensional sampling points, calls on CECIL(2) for the values of the integrands at the sampling points and inserts these values into the sums S 1 and S2. CECIL(1) returns the overall 2/m for factor thatand instance) multiplies also the each valueintegral of the integrands (~p~h for the sampling point (0,0,. . ., 0). FOURBODY can be used with any V(r) and f(r) or, for that matter, for any six-dimensional integral, with only the following changes to be made by the user: the variables NIN and NOUT, which define the input and output units, must be set for the user’s machine and the FORMAT statements # 660 and # 680 must be adapted to the form of the user’s output. Much of the subroutine CECIL must be modified by the user to generate the particular integrands in which. he is interested, though the logical structure of the subroutine should remain intact,
6. NUMERICAL RESULTS
217
arising from the attractive and repulsive portions of the alpha-alpha potential.
7. GENERALIZATIONS AND IMPROVEMENTS We can suggest a number of modifications which could be incorporated into the program to enhance its usefulness. In addition, the program can be modified to evaluate almost any manydimensional integral. 1) If no more than a few tens of thousands of sampling points are evaluated and if the user’s computer carries at least eight significant figures, it is unnecessary to carry the sampling points in double precision and probably unnecessary to retain the running averages (described by comment cards in the program) to reduce round-off errors. 2) At present there is no provision for restarting the program at the N + 1 sampling point after it has run and printed out the results for N points. If the printout included the values of VONE(N) and VTWO(N), then the program could easily be modified to restart at the N + 1 sampling point. 3) With little effort, FOURBODY can be generalized to perform any n-dimensional integral, where n = 4, 5, 6 For n ‘ 3 it would be more efficient (greater accuracy for given running time) to use one of the more common methods, such as gaussian quadrature. We would like to emphasize again that the present program, though it may be useful to some researchers as it stands, is meant primanly as a guide to the use of the Haselgrove method.
The program included with this paper is used to evaluate a fictitious problem. We could have supplied the CECIL subroutine that was used for alpha matter [1], for which V(r) is the difference of two gaussians and f(r) is a third order polynomial in exp (- r), but it would be quite difficult to evaluate the integrals exactly for this case. Instead, we estimate the integrals that arise for the artificial case where 2 I (r)V(r) = ~(r) = exl (- r). (18) The four integrals listed on the computer output correspond to the four terms in brackets in eq. (2); they are labeled RING, DIAGONAL, OPENEND and CONNECTED on the printout. These names arise from a graphical representation of the boson cluster expansion. The running time on a CDC 3800 is approximately 2.7 seconds per one thousand sampling points. An examination of the output indicates that 2000 points is sufficient for 10% accuracy, 4000 for 5%, 10000 for 2% and 20000 for better than 1%. For more realistic cases the accuracy is somewhat reduced, perhaps halved. In the case of alpha matter, for instance, there is con-
[1] F. Iwamoto and M.Yamada, Progr. Theoret. Phys. 17 (1957) 543. [21 3. B. Aviles, Ann. Phys. 5 (1958) 251. [3] G. P. Mueller and 3. W. Clark, Nuci. Phys. A155 (1970) 561. [4] S.-O. Backman et al., Nucl. Phys. A130 (1969) 635. [5] P. J. Davis and P. Rabinowits, Numerical integration (Blaisdell, Waltham, 1967) P. 150. [6]P.3. Davis, Interpolation and approximation (Blaisdeli, New York, 1963) p. 354.
siderable cancelation between the contributions
[7]C. B. Haselgrove, Math. Computation 15 (1961) 323.
REFERENCES
218
G.P.MUELLER TEST TSST MO.
RUN
OUTPUT
c OF tOUaROOY W
POtr.TSs 04 25o120060
CO.
N
RI40
DIAGONAL
LI
RlhG
OI~GONAL
OF
IWEGRALS=
‘
OPENFNO
OPEWNll
CGkNECTEO CGhNEClEO
2TS.@D2
33.062 33.062
33.lw2 33.062
356.RCQ 4?9.566
44.019 49.18?
43.259 46.1hn
7.253 8.208
500
295.1_14 342.347
37.256 41.479
37.962 42.?7?
5.V06 6,986
FL0
3nl.lPV
5.405
3ir.s9n
35.943 38.403
37.396
75n
39.311
6,181
loon
291.5e7
36.976
37.159n
5.748
loon
3114.040
37.371
38.415
5.035
12% 125tl
240.941
36.514
3no.29n
37.004
36.1~7 37.Rnl
5.518 5.739
15nn
273.S13
34.948
35.112
5.188
15on
294.473
36.574
37.111
5.614
1750
2C7.7~0
34.5P4
1750
2PR.271
36.118
34.37n 36.974
5.031 5.473
2non 2060
271.131 2P4.4?2
IS.572 35.957
33.*9?. 36.n~r
5.300 5.432
22% 22%
26n.541 290.476
34.396 35.742
33.40= 35.560
5.093 5.363
2soo 75nn
2C0.910 Z?'I.PPO
34.6A3 35.422
33.765 35.179
5.3Sl 5.321
2750
2r4.407 273.9!2
34.689
33.47n
5.355
35.321
34.911
5.325
3001) znotl
2Lfl.944
34.306
33.417
5.2b3
27P.ZPO
35.164
34.67.I
5.324
3250 3250
255.331 269.941
33.1593 34.991
33.7a7 34.*m
5.134 5.302
3560
2&?.9?b
5.331
2**.5z2
34.150 34.R72
33.515
3500
34.3SL
5.295
3750 3750
2cU.915 267.612
33.826 34.766
33.nn3 34.mll
5.222 5.293
4000 4onn
24~.n24 zc*..*.
33.9134 34.634
33.337 34.002
5.156 5.279
4250 4256
2=9.w3 2h5,q72
34.038 34.557
33.249 33.9c(b
5.143 5.263
4500 l5nn
2ch.7lb
33.642
2rr.790
34.495
33.+m 33.9nR
5.125 5.251
4150 l15n
2C4.614 **1.7:7
33.663 34.41s
33.122 33.130
5.061 S.234
5000
2c5.241 ZrP.wL
33.563 34.330
33.3Rb
5.039
33.774
5.216
.%S.dfi?
2-n 750 500
77so
5000
(EXACT (EXACT
VALUES) VALUESI
FOUR-BODY CLUSTER CONTRIBUTIONS 9000 9000
257,711 255,792
33,444 33.565
33,405 33,549
5.043 5.081
9250 9250
257,043 235,882
23,496 33.560
33,356 33.350
5.045 5.078
9500 9500
256,970 255,979
33.53* 33,560
33.350 33,355
5.050 5.077
9759 9750
257,136 256,021.
33,459 33,536
33,465 33,356
5,021 5.075
10000 10000
257,116 256,061
33,467 33,550
35,370 33,364
5,032 5,072
19000 19000
253,480 253,048
23,i80 33.414
33,226 33,280
4.988 3,050
19250 19250
254,444 255,020
33,277 33,410
35.198 35,279
5,004 5,048
19500 19500
253,626 234,995
33,244 33,400
33,191 53,276
4,995 5,041
1.9750 19750
252,900 254,947
33,165 33,400
35,1.92 33,214
4,982 5,040
20000 20000
253,433 254,988
23,181 33.395
35,1.28 33,272
4,990 5,044
219
220
G. P. MUELLER PROGRAM LISTING
ABIAFOURBODY. RAPID EVALUATION OF FOUR—BODY CLUSTER CONTRIBUTIONS. 1 MUELLER, 0.?. %JOB,66H0106,2780PM,4 %FI’N,L,x PHOORAM FOURBODY
0 o 2 3 ~ 5 6
I
C C C C C C
C a C C C
COMMENTS CARDS WILL ALWAYS FOLLOW THE STATEMENTS P11 WHICH THEY REFER
11
C C
10 11
C
INITIALIZATION SECTION”*..
CALL CECIL(1) VARIABLE NAMES ARE GENERALLY SET OFF BY SINGLE DASHES
•
CECIL) 1) RETURNS TIE OVERALL FACTOR MULTIPLYING THE 1. TN INTEGRAL AND THE VALUE OF THE L TH INTEORAND EVALUATED FOE THU N—O SAMPLING POINT. THAT IS, FOR ALL OF THE INTEGRATION VARIABLES SET TO ZERO. THESE TWO QUANTITIES APPEAR IN VALUE(L+1O) AND VALUE)L+20), RESPECTIVELY, FOE L=l.NINT.
12 13 C C C DIMENSION THETA(6),T)6)
rr(6),vALUE)3o),sONU(9),STwO)9),AvE(g),
1CESAR)9),STJRE(9),RUNAV~9) CO~8)ON NINT,T,VALUE,NIN,NOUT TYPE DOUBLE TT,DEL,QUIFK,THETA ..—-IDENTIFICATION OF DIMENSIONED VARIABLES—.— —THETA— IS AN ARRAY OF SIX IRRATIONAL NUMBERS. —TO— CONTAINS THE CURRENT (6—DIM.) SAMPLING POINT. —THETA- AND -TO- ARE CARRIES IN DOUBLE PRECISION. —T— IS AN ARRAY CONTAINING THE CURRENT SAMPLING EDItOr IN SINGLE PRECISION. —DOSE— AND —STWO— ARE THE QUANTITIES DEFINED BY EQS. (14) IN THE LWP. —AND- IS THE AVERAGE VALUE OF THE INTEORAND OVER THE SAMPLING POINTS AND —CESAR— II THE CESARO MEAN OF THAT SEQUENCE. THEY ARE DERIVES FROM —lONE- AND —ITWO—, RESPECTIVELY. —VALUE-, -HUNAV— ASS —STORE— ARE DESCRIBED LATER.
C C • •
114 15 16 1~ 19 20 21 22 23 ~ 25 26 27 211 29 30 31 32 33 35
50 300 L=1,NINT RUNAV(L)—O.G DONE(L)—VALUE(L+20) 300 STWO(L).SONE)L) C C
c C C C C C
Soc
C C
35 36
37 **~I00INPUTSECTIO
35 39 SO
41 1’ c NIH.60 NOUT=61 o C —BIN- AND —HOST— LABEL THE INPUT AND OUTPUT UNITS. 1420 C C READ )HIN,160) C C THIS READS ONE CARD. COLUMNS 2—52 OF WHICH CAN RE USED AS A C HEADER FOR THE OUTPUT. 55 C C READ )NIN,17o) THETA C C HEADS DATA CARDS TWO AND THREE, WHICH CONTAIN SIX IRRATIONAL C NUMBERS IN DOUBLE PRECISION. HAIELOROVE GIVES SETS OF C IRRATIONAL NUMBERS FOR INTEGRALS OF VARIOUS DIMENSIONS, BUT C THEY POSSESS ONLY EIGHT SIGNIFICANT FIGURES. INSTEAD, WE USE C THE NON—INTEGER PORTIONS OF THE SQUARE ROOTS OF U. 3. s. C 7, 11 AND 17
C
66
O 80 CONTINUE READ )NIN,200)
MAX,JAX,NINT,NONE
C -——IDENTIFICATION OF INPUT VARIABLES—--MAX— IS THE ~IYFTALNUMBER OF SAMPLING POINTS TO BE USED. THE PROGRAM PRINTS OUT ITS CURRENT ESTIMATES, THAT IS, THE CURRENT VALUES OF —AND— AND —CESAR-, EVERY —SAX— POINTS. -HINT— IS THE NUMBER OF INTEGRALS TO NE EVALUATED. IF —NONE— IS ZERO OR H NEGATIVE INTEGER, THE PROGRAM TRANSFERS TO END AFTER SAMPLING —MAX— POINTS. IF —NONE- IS A POSITIVE INTEGER, THE PROGRAM TRANSFERS TO 80 AFTER SAMPLING —MAX— POINTS AND ATTEMPTS TO READ IN NEW DATA. AFTER EVERY MULTIPLE OF —SAX— POINTS, —RUNAV— IS SET ESUAL TO —CESAR— FOR 1.—i ,NINT. THEN —SOME— IS MODIFIED SO THAT IT ACCUMULATES THE THE DIFFERENCES BETWEEN THE CURRENT SAMPLING POINT AND —RUNAV—. THE PURPOSE OF THIS MANEUVER IS TO REDUCE ROUROOFF ERRORS. SEE THE LWP FOR A DISCUSSION. -STORE- IS USED IN THE EVALUATION OF -RUI6AV—.
C
42 43 44 45 46 4~ ~9 50 51 52 53 54 1115
57 511 59 60 61 62 63 64 65 6~ 6 69 70 71 72 73 74
77 7~ 14 So
81 82
83
1514 85 96
87
88
89 WHITE WHITE WHITE
NOIJT,160 NO1JT,220 NOUT,240
JAX,MAX,NINT THETA
THE THREE WRITE STATEMENTS IMMEDIATELY ABOVE WHITE OUT THE HEADER CARD, THE ARRAY —THETA— AND THE VARIABLES —MAX-, -SAX— AND —HINT—.
90 91 92 93 94 95 96
97
95
C C C C
C C C C C C C 460 C C C C
480 490 500 C C C C C
131 137 133 134 1 35 LOOP 300 INSERTS THE N=C CONTRIBUTIONS INTO THE SUNS —SONS- AND 1 36 —GTWO— AND SETS THE INITIAL VALUE SF THE RUNNING AVERAGE 22 13 zERO. 13 139 140 1141 N=O 1142 143 —N- IS THE INDEX THAT LABELS THE CURRENT SAMPLING POINT. 144 145 1146 147 1IAX—JAX 14i 149 AF’I’EE EACH SUCCESSIVE SET OF —SAX- POINTS HAVE BEEN EVALUATED, 150 THE PROORAM EXITS FORM ITS NORMAL ROUTINE, PRINTS OUT ITS 1 51 CURRENT ESTIMATE OF THE VALUE OF THE INTEGRALS AND RECALCULATES 1 52 THE RUNNING AVERAGE, —IRIJNAV—. —HAS— IS THE INDEX WHICH THE 153 PROGRAM TESTS TO SEE IF THE EVALUATION OF ANOTHER -JAX— POINTS 1 58 IS COMPLETE. 155 156 15 15 DO 420 K-i ,6 159 ‘117(15)—SOD 160 161 LOOP 420 SETS THE INITIAL VALUE OF THE SAMPLING POINT AT ZPEC. 162 163 1614 165 166 16 16 169 EVALUATIOH SECTIO 170 171 172 173 N-N*1 174 175 INCREMENTS THE SAMPLING POINT NUMBER, -N—. 176 177 178 179 OG500K—1 6 180 T’r)K)—Tr(lc)*THETA(iC) 181 QUIRE—TT(i5)_1.OS 182 IF (QUIRK) 490,490,580 183 l’r(K)—QUIHI( 184 T)K)’I’T)K) 185 CONTINUE 186
18
LOOT 500 CALCULATES THE COMPONENTS OF THE NEW SAMPLING POINT AND STORES THEM IH -IT—)IIJUBLE PRECISION) AND IN -T- (ROUNDED OFF TO SINGLE PRECISION).
FORMAT 15,51 H) 234567890123456789012345678901 23456789012345670901) FIRMAT 3150) FORMAT 110,315) FORMAT 18H0tHJ. OF POINTS. O(.I4,1H),I5, 15X,17HNO. OF INTEORALS—,12) 240 FORMAT (26NOSE,r OF IRRATIONAL NUMBERI/(3D30.17)) THE FIVE STATEMENTS ABOVE ARE THE FORMAT STATEMENTS FOR THE INPUT SECTION OF THE PROGRAM.
99 100 101 102 103 104 105 106 07 1001 109 110
18 1 190 191 IS’
193 1914 195 196 19 19 C 199 C ~oc SO 520 L—1,HINT 201 SONE(L SDRE)L *VALUE(L)_NUNAV(L) 202 520 STWO)L —STWO(L *SONE(L( 203 1 ‘04 C LOOP 520 ADDS THE CONTRIBUTION OF THE NEW SAMPLING POINT TO THE 205 C SUMS —GONE— AND —STWO—. NOTE THAT THE RUNNING AVERAGE IS 206 C SUBTRACTED FROM THE VALUE OF THE CURRENT SAMPLING POINT. ‘0 C
CALL 011114(2)
C C C
CECIL) 2) CALCULATES THE VALUE OF THE INTEGEANI AT TIlE HEW SAMPLING POINT AND RETURNS IT IN VALUE(L), FOR 1.—I ,NINT.
a C
160 170 200 220 O C C
C C C C C C C
75 76
iii 112 113 115 1 5 116 11 11 119 120 121 122 123 125 125 126 127 128
2111 IF )N—NAx) 460,600,600 TESTS WHETHER THE INTEORAND HAS BEER EVALUATED FOR A COMPLETE SET OF —SAX— POINTS. IF NOT, THE PROGRAM RETURNS TO 460 AND EVALUATES ANOTHER POINT. IF SO, THE PROGRAM EXITS TO 600.
C C C C , , ,
C r ~ C
Note: Columns 73 to 76 of each card image are omitted to save space,
~12 213 214 ‘15 16 21 1 219 ‘20 221 222
FOUR-BODY CLUSTER CONTRIBUTIONS ~
~
C
2 .
600 TRIO~2*N+1 TRIT...(N+i) *(~+) 1’WO(L)/TRIT STORE L)-S AVE(L (IONE(L)(TRID+O.5*RIJRAV(L)).VALUE(L+1O) 620 CESAR L)—(STORE(L)*0.5.RUNAV(L) ).vALUE(L+iO) C C C C
227
LOOP 620 CALCULATES THE ESTIMATES —AVE— AND —CESAF— FOR THE VARIOUS INTEGRALS. WRITE (HOST,660)
N, (AVE(L) ,L-i , NIMI)
WRITE (NOIFT,680) N,(CESAR)L),L1,NINT) a C
230 ‘31 232 233 234 235 ~
THESE WRITE OUT THE CURRENT ESTIMATES -AND- AND -CESAN-.
2~lO 241 242 243
c
C C
660 FORMAT (1H0,I9,5Fi5.3)
‘45 246
680 FORMAT (Ii0,5F15.3)
247
THESE TWO FORMAT STATEMENTS FOR ‘THE PRECEDING WRITE STATEMENTS MUST BE MODIFIED IF MORE THAN FIVE (5) INTEGRALS ARE EVALUATED OR IF A DIFFERENT OUTPUT FORMAT IS DESIRED.
C C C C
249
550 251 252
C C C
255 256 257
LOOP 760 CALCULATES THE NEW VALUES OF THE RUNNING MAKES TNE APPROPIATE CHANGES IN -lONE- AND -ITWO-.
~
C C C
a
C C C C C I
~TEMBFO a C
(
C C COMMON NINT,T,VALUE,NIN NOUT DIMENSION VALUE) 30) ,T( 65 ,R(6),F.(6) C C C C C C C C
———IDENTIFICATION OF SUBSCRIPTED VARIABLES—-— -VALUE- AND -T- ARE DEFINED IN THE MAIN PROGRAM. THE COMPONENTS OF —F- ARK THE RADIAL SEPARATIONS FOR WHICH TEE INTEORANDS ARE EVALUATED. THE COMPONENTS OF —R— ARE ORTAINED FROM -T- BY CHANGES OF VARIABLE. IN THIS CASE, H(K) — —FOWER~NATURALLOU(T(K)), E-1,2 AND 3. —POWER— IS A VARIABLE THAT IS READ IN BELOW. R)14) R)5) AND R)6) ARE OBTAINED FROM R)i), F(S) AND R)3) AND T(45, T)5) AND T(6) IN THE MANNER DESCRIBED IN THE LWP.
C C
)*(pOwEp*powER.pr,JER.3. *3. .TPI).(O.
C
NOTE——GENERALLY THERE IS A DIFFERENT FACTOR FOR EACH OF THE IPEDRALS.
350
(
(
9‘ C C
260 C C C C
*
358
359
*8*t
4
269
36~
VALUEçL+20~—0.5.10. 0)
365 366 367
VALUE~L+10-FACT
RETURN •
*
*#******4**•*** *
270 271 272 273 174
*8
4
*
368
369 370 371 372
373
374
375 376 377
*
*8*. 378 379
****8*******t*****•******** **~***~********~*
****..*
380
381 382 383 384
*****SECTION TWO*•***
275 276 277
*
*
386 387
A**.n*
388 520 R(I ——POWE~L4DGF(T(I)> I~R 1
391 352
281
Y—R 2
285 286 28
QYZ—SQ,1T5t~T( 5 >+. 25> AYZ——1 . 5+QYZ ATh—m*Afl.5927( Ii .—AXY.AXY)*(
396 3 3 399
1 .—AYZ*AYZ> ).C05(FIWT(6>)
XS~X*X YS—Y*Y ZSZ*Z N 4 —SQRT XS+YS—2.*X~Y.AXY R 5 —SORT YS+ZS—2.*Y.Z*AYZ R 6 SQJ1T Z5+XS—2. .Z~j.AZE
289 ‘go ‘91 292 293 ~
C
295 ‘96
C C
‘TI 29~ 299 300
C
4 401 402 403 404 405
STATEMENTS 500 THROOUN THE LAST CALCULATE THE RADIAL SEPARATIONS IN TERMS OF THE COMPONENTS OF -T- AS DESCRIBED ABOVE AND IN ‘THE LWP.
406 407 408 409
**.****.**.*44
410 1412 413
C C
411 .
V—v.V/(T(1 )wr(2)..T(3>.qxy.QYZ) C C C
414 415
—V. IS ESSENTIALLY ‘TIlE JAC1BIAN REPRESENTING THE ABOVE CHANGE OF VARIABLES. IT IS USED TM t.~P 720.
4*8*8*8*88*
fh
C ~
00 660
416 41 81 819
an
—
(R(I~—36~)620,640,640
IF
C
311 312
GO TO (101,500),M0DE TRANSFERS TO 100 OR 500 ACCORDING TO WHETHER CECIL) 1) OR CECIL(2) IS CALLED. • •
314 315 316 31 ~
C C C
),
325 326 327 32~ 330 331 332 ~
422
423 424 425
I.~F 660 CALCULATES EEP)—P(K>> POP 11—1,6. EXPONENT UN~YL0W IS INCLUDED. *
a.
A26
42 42
A CHECK TO PREVENT *
n........n...
Y—E(1 .E(2>.E(6) VALUE I .T.E( 5 VALUE I? —VALUE i).E
322 323 324
100 READ (HIM, 120) POWER WRITE (ttou’T,14o) POWER WRITE )NOUT,160) 120 FORMAT P5.2> 140 FORMAT 6HIR — —.F3.i,23H • 1.00) Y(I) 1—1,2,3) 160 FORMAT 151 18X,1HN,11X,INNINO.7X,8HDIAGOHAL,8X,7HOPENEND,6X,9HCONNECT!D,
k21
620 E(I>—EICP(—R(I)> 00 TO 660 6140 E(I)—C.0 660 CONTINUE
C SECTION ONE
360 361
00 260 L—i HINT
LOUP 260 PLACES THE OVERALL FACTOR -FACT— IN VALUE(I,f10) AND THE VALUE OF THE INTFINRAND FOR THE N—O SAMPLING POINT IN VALUE(L+20~, FOR L—i,NINT. THE FACTOR OF 0.5 IN SHE VALUE~L+20, STATEMENT ARISES BECAUSE THE N—C SAMPLING POINT ONLY RECEIVES HALF THE WEIGHT OF THE OTHERS. IN THE CASE PRESENTED HERE, THE VALUES OF THE INTEDRANDS FOR THE N—O SAMPLING POINT ARE ZERO.
•
302 303 304 305 306 307 30~
5)
IT IS MADE UP OF POUR FACTORS. (2.> ARISES BECAUSE EACH 351 SAMPLING POINT, OI~ TEAM THE N-C POIN1I’, HAS A WEIGHT 2 IN ‘THE 352 -SOME- SUM. (2. .‘I’PIrrPI) ARISES FROM VARIOUS ANGULAR INOTMS353 RATIONS (SEE EQ. 6) IN THE LW). THE NUT FACTOR ARISES FRUM 354 THE CHANGES OF VARIABLE THAT ABE DESCRIBED ABOVE MID IN THE 355 LW. TEE LAST FACTOR IS ONE—HALF RHO CUBED, WHICH APPEARS IN 356 EQUATION 2>, POE INSTANCE. WE HAVE HERE SET RHO—i .00. 357
266 267 204
301
345 346
—FACT— IS THE OVERALL FACTOR THAT MULTIPLIES EACH INFE~RAL.
,~
SUBROUTINE CECIL(MODE
34~
t4.***..
C C C C C C C
C C C C
•••
3~O 341 342
C
~
C C o C o
338
344
C C
STATEMENTS 800 THROUGH 900 INCREMENT —NAB— BY —SAX— AND RETURN It THE EVALUATION SECTION IN ORDER TO EVALUATE —SAX— MORE POINTS. IF —MAX— POINTS NAVE ALREADY BEEN EVALUATED, THE PROGRAM CHECKS TO SEE IF A NEW SET OF INTEGRALS IS TO BE CALCULATED. IF NOT, IT TRANSFERS TO END. IF SO, IT TRANSFERS
*
;~1~,~527 FACr—(2. ).(2. ~rpIwrpI
262 263 264 265
800 NAX—NAX+JAX IF (NAX—MAX) 460,460,920 820 IF (NONE) 900 900,840 640 WRITE (~ioirr810) 860 FORMAT (IHi ) GO TO 80 900 CONTINUE
337
POWER, WRITES IT OUT AND WRITES OUT THE LABELS POP ‘THE OUTP1.rT. •
AVERAGES AND 260 261
C
41 1TNG,~T IIDIAOOHAL,8X, nJrw.nW,6X,5a.tau.SJr5u,
517X,8H 255.802,8X,7H 33.062,8X,7H 33.062,24X,14H(UAC’r VALU!5>/ 617x,8H 255.802,8X,7N 33.062,8X,7H 33.062,24X,IIIH(EXAC’r VALUES)/) THE ABOVE GROUP OF STAT~IENTSREADS IN THE VARIABLE
,~
00 760 L-1,NINT RUHAV(L)—RUNAV(L)+2.*STORE(L) SOME(L (—loNE) L)-TRIO.E’TORE(L) 760 STWO(L (—0.0
~i~ifl
221
3> VALUE 3 Y*E(3 *E(4 VALUE 14 —VALUE 3>*E 5) C C C C CO C C C ~
‘
THE ABOVE FIVE STA72MENTS STORE THE VALUES OF THE POUR INTMSRAHDS IN VALUE(L), L—l 4. THERE REPRESEMr ‘tEE POUR TERMS IN THE INTHORAND OF EQ. 25 IN THE LW. •
1
*8*88*
~
*
*33 838 435 836 837 438
839
it4o
441
442 44
222
G. P. MUELLER
C
C
EACH INTBDRAND By THE JACOBIAN REPRESENTING 447 THE CHANGE OF VARIABLES DESCRIBED ABOVE. 448 SOUP 720 MULTIPLIES
449
*
a a
450 451 ‘452
RETURN END a *
* **
1453 454 455
1
1456 ‘457
1 1
SCOPE %LOAD %RUN,4,6000 TEST C OF FOURBODY .414213562373095049 .732050807568877294 .236077977499789696
458 459 460 461 462 463
.645751311064S90591 .31662479035S399849 .1231 05625617660550
1464
20000 3.0 •
250
4
0
465
466
467
1 1 1 1
1