- Email: [email protected]
A general purpose code for Monte Carlo simulations
Nu*.'l¢,~e h;~lrtm~.l~, arid Method,, in I*hys=,..~,R¢,,ear,.h 22S ( 19841 132 L,t? ,~{lfih-I hdlaDd,/'.rtl'~lcr,|:tl~,
132
A GENL~AL Wmfriec
PURPOSE
~. WILCKE
CODE
FOR MONTE
CARLO
SIMULATIONS
*
Nuclear .(;t ": cture Research Lal~rato~ and Depaftmen~ of Chemis: O" and Ph)'~ic,~, C'nit~v,~i(t' . / R,*ch¢~ter. Ro~ h,'~ler. New York , t627, tJSA
Received ;:' )uly 1983
A gcn~: ,I-purpose computer code M O N T I I Y has h c c . w/httLL hJ t),t'lfl~rl|, ~'l,~l,te (~';irh, ',;mul;.,tit,lc- ~',f i*d=y :t.~l K~r'lLrll%. TO a c h _ v e a high d e g r ~ of heJtibility the c(~le i'- t)rgAnized liJ.e a Bcneral purptJ-,.¢ c o m p u t e r , oper.atin 8 t)n ,, ~.eclls! de',crih~n~ ~he time depei, lent slate o1 the syslc'm under s i m u l a : i o n , T h e instruc.ti3n ~ t <.; the " c o n t p u l e r " i', defined by it:' u~er and is therefore a d a p t a b l e . , the particular p r o b l e m studied. T h e ~ ' g a n i z a t i o n =)f M : ) N T t l Y allc,~.~, il©ralive and ¢onditio;'t;.:l ~.~ ~culi~'~n ~)1 ~)per;nit,n ~,
I. Inlrod" ~. lion
The [ : werful M o n t e Carlo methods have been applied to dye a wide range of physical problem,; [1.2]. One p a r ; ular area of success has been the modeling of systems ~y studying the time evolution of its c(.,nstitucnts A well-known example of t , i s kind is the design o' a nuclear reactor core by s i m u l a l i n g the fate of man~ hdividual neutrons as they are emitted, therrealized 'id eventually absorbed_ S o p h : ,teated special purpose codes have been puhlisbed t( ;olve such specific problems [3-5]. In general it is noZ ~articularly easy to adopt such codes to new problem, or change the models underlying the simulations [6 The code M O N T H Y discussed I'.ere i,'; an attempt t3 create a systematic approach to the simulation of ~ wide variety of physical processes with one commot approach. Since it is impossible to anticipate all possi : le models or problems a user of the MON"I'HY code m;-, want to simulate, M O N T H Y is not a readyt o - u ~ t,.de, but more akin to a computer language, which -~:cepts a c o m m a n d sequence s;.ored on disk (a "progr;.:S') describing the pl'ocess to be simulated in terms o: "instructions", which are to be defined by the u ~ r . q'L~.'se " i n s t r u c t i o n s " are introduced into the MONT'-IY code by adding user-supplied F O R T R A N subroul.~es to "he existing backbone structure of M O N T : I Y . A systematic set of simple rule,; d c ~ r i b e s how to :tdd these subrouhnes and how to prepare the c o m m a ~d files. T b e ~ rules arc described in the separate docum,;ttation of tl-~e code. In the present p;tper the more srmal concept underlying the s t t u c ' u r e ef M O N ' I ~'IY are to be discussed_ " Prose; ~ address: IBM, Thomas J. Wa:son Re.,-c.arch Lab:,ratory, ': ~)rktown Heights, New York 10598, USA. 0167-5( ;7./84/$03.00 © Elsevier Science Publishers B_V. (North i-!oiland Physics Publishing Division)
2. (hencr~l concepl.i ul M O N I ' t l Y
2.1. Introduction "Co illustrate Ihe ahstrat:t conccpt~ to be discussed below, lel us first ititr~duce ,a.,o specific example,, selected from applications of the co.de in nucle,~r phy~,ic,,. As the f~rst ez.ample we consider a highly excited nuch'u~ moving in the lai')orater3' system with a certain velc,c~ty t/ T,le nucleus :s reducing its hiL=h temperature I' by evaporating neutrons [7]. This process will continue unti; the excitation energy of the nucleus is I,:~., h~w Io .supr,ort further emissi(,n o f neutrons. Wc May ask f,:)r the ~:ngu[ar- and e n e r g / distrtbutiem (of the neutrons as obse:'vcd in the laboratory system. In a second example we want to study the prol)erlies of a large particle detector. Consider a tank filled with a .cintillating liquid, whicia is e x p o ~ d to a radiation field [8]. Several photomuh~plier tubes observe the interior of the tank through windows. W', want to calculate the probabilily th:n the photons created by ~ i n l i l l a l i o n s are actually de~ec'ed I'.y ope of the p h o t o m u h i p l i e r tubes, taking into account the absorption of the light in the liquid, multiple reflections o l the light from the walls and the efficiency of the phototubes. T h e ~ and m a n y other problems have in c,.~mm0n that it ;s possible IO de~crib¢ the Slate of Ihe system to be simulated by a suitably ch*~,,,en ~ t [,:) of N n u m - c r s x,. Ii is cc, r'=s'enient Io call this ~',et :.* "'*,late vector", even th(rugh it is not an element of Hilbert s p a t e and does not describe a system in the q u a n t u m m~h.',nical sen~,e. The emission of a neutron or the reflection of a photon from the wall of the rank can be considered as a transition betw,-.,'n two states of the ,system. We assign to each transition an operator Ox, which changes the
J¢'. H' ltgh.l-,, / M.nt¢ ('arlo itmu,"llttJn t'oth' H',: vvcl,W ,If the syslen~ from ,late ],~ ) to ,,tale ix'}:
; -'.., /6 :'-). ' l h e index K dvli.~tti~hes the kind of transition described hy an c,per;,t,,r. '[h¢ operators are not subject to any .~rmal restrictiol~s a . d m; 0 be non-linear functionals o f the s t a t e v e c t o r . COllt~liv ~toch;islic o p e r a t i o n s o r
wh;devcr else may be neces~,ary I~.) ~.intl,ki~c the system. I:ach transition may correspond to a dil!Cr,:~tkind of operation, or the same operator may he usud ~uF,c,ttedl2/ for different Iran.,,lliOllS ;In it is the case for s u b s c q u e l l l ClltiS'~ion OF severa] lleU[r~/ns f r o m a l:Ot nucleus. A gc,,eral ~lrillg o f ,'~ IransAiollS s l a v be writtcfl as
where the index KO) indicates what kind of" operation is used in the irh t~ansition. "lhe probleat is more complicated than discuss~:d so far, since the sequence of tr;m-itions is usually r:o~ predetermin.:d b y the structure ol the problem, but different for each event t t be simulated. In the second examplc cited a photon may hit a p h o t o c a t h o d e veo, soon. or the photon may undergo several reflections befo-e it is either detected or ab,,orbed_ It depends on the r a n d o m initial conditions of the first state and tile result', of subsequent stochastic ;ransiuons which sequence of operations is to be executed.
~i P~2;£ t': i(;:;
Ix (L:+2)> ~n÷2
1.13
Ihe haste concept of M O N T I I Y is to solve this control problem by structuring the code like a general purpo:.e computer. Each transition may be considered as the result of executing an instru,.:tion of a viwtual Von N e u m a n n machine, which operates on the stale vectors as data. Thus we may define conditional branchinl~ instructions, which allows the flow of operations to he altered as a function of the state vector (data) to be proce';sed. The set of all operators ( O~,. } corresponds Io the instruction set of the virtual nlachine, and each transition is the result of executing exactly one instructi,m of this set. In the physical code of M O N T I I Y each inslruction corresponds to a F O R T R A N subroutine Stq?p!i,:d I¢. the user. 2.2. Str, cim',l,, .l~erc~'lons In the following two seclions, st,>chastic and deterministic operators will be discussed_ Stochastic operators are normally us,:d to create ~,
sequence of random numbers with a spot;r:'.! distribution given by the predeternfined funct.:r, /,qq. winch describes the particular physical process to I~e simulated. The independent variable x is an element of tile stale vector. In the first example .x corresponds to the vel(~ity of a neutron and f(x) to a Ma~cwelliau distribution. The well known Von N e u m a n n method [91 is used to calculate the sequence of r a n d o m numbers with the desire:l spectral distribution. The meth~.~ is illustrated in fig. 2. If a function is bounded by the value y: in the imerval xm,,, ~< x ~< x.,~., then the Von N e u m a n n mc:thod can be applied. A random number generater is used to select uniformly distributed numbers in the intervals [xm, .. x . , . . ] and [y..,,.,y.,..], A pair (.x.,)b~' of r a n d o m numbers is used to test the validity of the inequali'.y/(:,:0) >_vo. If the inequality is true. tllen tiffs choice of x =.,c o is accepted, otherwise this trial is disc.',rded. It is plausible that this procedure ',,,'ill yield a
Ymox
• •"
~'~
[X (t~)>
iX (tt-I)>
t. o',; 0~ CALC~I.,' i lCa~;
F;g. 1. "lhe rclalion between die flow of calculations, stale vct:lor and Ihe Iran'.itkzn,, bcct~,'c~n sLalC$ iS illustrated. Each tr~ll,~it[on corre'.,p,cw,d, t,y the c.lecutiolt of ¢xacll_v one instruction of the "-irtual computer or. equivalenfl',, to the call to a
subrouline.
Ymi ~"// ~~ ///
Xmm
~,~2(x)xx~ Xmal
Fig. 2. lllustralion of the Von Neumann mezhod to gcncral¢ a sequence of random Lunll~rs with a spectral distribution given by .fix). (See t~xt for details.)
H ~. 14:
134
H/~h'k¢ / M,~nte Curb, tt,'4,1~ i,m
sr..quence .,f x-vahies with a p r o b a b i l i t y weighted by
fix). Thq Yon N e u m a n n a l g o r i t h m c a n ~ s i l y be exlended t~ ' u n c t i o n s o f several variables - a f-,ncfioa of m i n d e p r:dent variables needs Ihe g e n e r a t i o n of m + 1 r a n d o m ~"ambers f o r the test of o n e trial .sample. In general I ~:: functions te be simulated will be d e l x n d c n t o n o n e ( : m o r e p a r a m e t e r s . (The tet,.perature T of the n u c l e u s , vlitting a n e u t r o n is a p a r a m e t e r in the o r e r a lion of e : a p o r a t i n g a n e u t r o n . ) Since a p a r a m e t e r ,'nay c h a n g e ~! the result o f a transition (e.g. the nut leus reduces : ,: t e m p e r a t u r e by emitting a neutron), o n e has to be c~. :ful not to violate the c o n d i t i o n of the Von N e u m a n " a l g o r i t h m that max[f(xA]'cy,.,... This i~ illustratec ~y the d a s h e d c u r v e ]'2 in fig. 2, where for a n o t h e r , alue of the p a r a m e t e r tile m a x i m u m of the function ::xcecSs y . ~ . . The Von N e u m a n n a l g o r thin w o u l d y d a f i a t - t o p p e d d i s t r i b u t i o n for f 2 ( x ) , if y~.= were n o adjusted a c c o r d i n g l y . Since it is wasteft~l of c o m p u t e time to cl"~xtseym, unnecessarily large jr:it to be safe. is desirable to ~ t ) , ~ equal to the maxim:~m of the k ction, but this n u m b e r is not k n o w n a p r ori. "1"o ~ ¢e this p r o b l e m in an a p p r o x i m a t e w a y the f o l l o w i n , a l g o r i t h m has been a d o p t e d in M O N T H ' ( : A start val .~: for y m . is c h o s e n as v,~,~ = a - m a x ( / ' ) . ~ here a = 1.2 -: a safety f a c t o r a n d m a x ( / ) is c a l c a l a l e d for the initi , values or the p a r a m e t e r s . W h e n the p a r i m e ters c h a "~;e, the c o d e keeps track of the largest v a h e of fix) ever . n c o u n t e r e d a n d adjusts ).',~. u r ' v a r d s if n :cessary. F c " functions which are c o n t i n u o u s with r e s l ~ c l to all Ihei; p a r a m e t e r s this p r o c e d u r e will lead to very small en- )rs if a sufficiently large n u m b e r of events is sampler W h i l ' it w o u l d have been possiblc t~ store the p a r a m e . " -s as elemenls in the state "urine, it was f 3 u n d to be m ,:e c o n v e n i e n t to associate with e a c h i n s t r u c t i o n step it'. he " p r o g r a m " of the virlua', m a c h i n e a p a r a m e ter vect , which c o n t a i n s all p a r a m e t e r s for the c p e r a for c a l h , for b y this instruction. Obviously, e a c h pattie-
, ¢~d,,
u l a r o p e r a t i o n has its ov, n f~)rmat for Ihe inlerprelaticm of its ~,~,socialcd p a r a m e t e r vector_ In general the n u m e r i c t l values of the elements in the p a r a m e t e r vector will be different for re'pealed calls of the same o p e r a t o r W i t h the i n t r t ~ u c t i o n of the p a r a m e t e r vector r, the .,~quencc o f n ~ransilions c a n he '~.'rltten in it.,, final f o r m as
I., -t')> =
~^.,.,( t%.~)., b~,,,( P^-,,,) ~'~").
where e. ch of the indices 1 - - . n c o r r e s p o n d s 1o Ihe execulio:~ of o n e instruction in the virtual ntaehlnc. In the pr~.'~ding discussion a d i s t i n c t i o n has been m a d e bq:lween argumenls Of the ',lille vec[~); a n d p : l l a n w l e r s .
While tl is dlstinclion a p p e a r s q u i t e 3 a l u l a l l y in nto~,l real prol.lems, it is really s o m e w h a t : b i t r a r y T h e s a m e q u a n t i t y m a y be a i n d e p e n d e n t vat!':. ~le for one transition an(' a p a r a m e l e r in a n o t h e r t~::e, r)r a user m a y decide t')at he w o u l d r a t h e r not k ''.rp a q u a n t i t y curlslant, b~:t allow it t o fluctuate, hi o r d e r to .'qkps. for s u c h e x ( ' p t i o n s w i t h o u t the need (',f rewriting o p e r a t o r s ( s u b r o u L n e s ) , two special o p e r a t i o n s ha~,e been introd u c c d , v hich allow an e x c h a n g e of a r b i t r a r y elements in the slat( vector a n d the p a r a m e t e r vectors.
2-3. Dee rminLrttc operolor.s' Seve~ al deterministic o p e r a t i o n s arc generally ncce~,s a r y to solve a p r o b l e m . [hese o p e r a t i o m , yield ten o u t p u t .~tale vector_ svhich is u n i q u e l y dt:tcrmined b y the i n p u t st'de vector, the o p e r a t i o n executed a n d the ass()elated i: ~rameters. Typic:t] examples of these o p e r a t i o n s are kil~:malic;t] c:dculvlion~, an, I the ev:du:ttion o f n , a t h c n . a l i c a l functions. O t h e r o p c r a t i o w , i0 this clas~ tu,-'h,de b r a n c h i n g instuctions, the sorling of events into s p e c t r a a n d c o r r e l a t i o n a r r a y s a n d the g a t i n g of events. An ovcJ view of all presently defined classes of o p e r a t o r s ( i n c l u d i n g the stochastic ones) is given in table 1,
Table 1 Overvie~. of pre.s~ntly define.d class,'..s ol operators in MC N'I-HY. The clas,~ numb.crs refer Lo internal operation c~x.les and the colulrm labeled ~ gives the increment of the prol~rarn counter :3c for ol~ralo[s in each class Note thai I'~ rcfe[s to a p~,mter in the (tMe MON'I'] : f -,bhich selects the oI'~]ations !o be carried ola and not to t~.: hJrdwired prograln c.untcr r~l the .icm:d c,m,lmwr wlm'h execut~ MON'[ tlY. Class
N a m e of cla%
pc
.......................................................... "l')'pical c',amples
1000 4(}~)
Stochastic operators l')¢terminlstic o p e r a t o r s P a r a n ~ t e r ~ g,t a t e ) B r a n c h instructions E v e n t sorting Mi~ellaneous
4- 1 4- l + l vat. + 1 + 1
[;valc)oratlOll ~ F l r u [ l l , Gau~vJall ~|)¢~[ruT[] Kinematics, algebraic expre~slons Exchange elements in i-") and 9 Branch if x, - 0 Sorl events inlo a cot.clarion array Next event, gating cc..)dlti~n
7000 g0130 9000
~'~/ ~'~', t~"l/l #', I' J ,~/or//t I ('~{~/,~ .~l,llrld/'l:loll (',¢k(11'
ht~idVl,.~L.ntaiionof c,,d¢'
it,
A';. h;t> I,C~H ,li*;cussed above, the f u n d a m e n t a l ccltc c p t of M O N ' I {{I.Y i'. (o create a t'xw_~lof o p e r a t o r s (~b," w h i c h form the instrucli ;n ,~et of a virtual Von Nt.u, n a ' m m a c h i n e . In fig. J Ihc c,','crall o r g a n i z a t i o n of :his " ' m a c h i n e " is s h o w n . T h e plo~r:,m c,,tmter (pc) of the virtual "'rl{achine" is n o r m a l l y increl~.,,.mcd .by o n e after e a c h o p e r a t i o n a n d I'~)i its to the next "Jl~-,l!'=J,Jli(al. I]r:mch il131:lUCliOns c a n a,,',ign o t h e r value: to the 13rolgr;lll CtHlltter. As~,ik,:iated with e a c h inlruclioll is a paranlci._'r vector, w h i t h conluilLs the i n f o r n t a t i o n r e q u i r e d for the o p e r a t i o n t~ he cxecuied d u r i n g this i n s t r u c t i o n . T h e d , d a ,nemc~ry of Ihi.,, virtual m a c h i n e is very small a n d c o n t a i n s only the slate vector Ix.~, As has been discussed, there are special operatio,~s w h i c h allow exc h a n g e of elements of tile state vectt:r with elements of ,m a r b i t r a r y p a r a m e t e r vector. The i m p l e m e n t a t i o t | of the vir,ual m a c h i n e in a F ( ) R T R A N c~x.le is ~,traightfor~,,ard. A c o d e n u m b e r is assigned to e a c h kind of o p e r a t o r , wl,ich is stored in the " p r o g r a m memory"_ T h e p a r a m e t e r s are stored in a m a t r i x (actually two matrices, one for integer n u m b e r s a n d o n e for floating p o i n t n u m b e r s ) . E a c h row of these matrices c o r r e s p o n d s to a p a r a m e t e r vector, with the p r o g r a m c o u n t e r IX: p o i n t i n g to the next row to be a c c e s s e d U p o n execution of the p r o g r a m fo," the virtual rnachine Ihe p r o g r a m c o u n t e r is initialized to one. a n d the first i n s t r u c t i o n c o d e n u m b e r is fetched f r o m the p r o g r a m m e m o r y . In the actual F O R T R A N c o d e a c o m p u t e d G O T O s t a t e m e n t is use~J to call the first
C
CODt" Op~c~c' in= OK{=)
P};[l)
PIK ()1
3
Opce~e ~o[ OX(])
PK{3)
(::~,tr'-x ~ _'atr'.×)
..L2"._':'.~:L~_JAL'C_ ! x > ~- l x i ,
x2, . . . . . . . . . . . . .
LN}
IzxP)
Filz. 3. The or~,amzation of the MO/ITI-IY ctw.l¢ follows the ,,Iructure of a general p u r ~ cc,mpul4 r with separate program Inemo[y
a=id d a l a lUelTIOly_ "1 h e d a t a n : t m u r y
contains
only the
pre.wrnt slate vector ~.,¢>,whereas the pxogram memory_ contains the int,.lruction co,de and the pa,am¢l¢/' vector.'=.
] 35
operator, which is represeqted by a subroutine of M O N T I I Y . After the return from this s u b r o u t i n e (i.e_ execution of an instruction in the terms of tile ",'irtuai m a c h i n e ) the p r o g r a m c,..unter pc .)f the virtual m a c h i n e is i n c r e m e n t e d b y o n e (except for b r a n c h i n g operations), a n d the cycle is repeatec, for the next operation. Somew h e r e d u r i n g the sequence of instructions the d a t a in the state vector are used as addresses to i n c r e m e n t the a p p r o p r i a t e spectra a n d c o r r e l a t i o n arrays. By default every event i n c r e m e n t s the selected channel by unity. b u t an option exists to use a n a r b i t r a r y element of Ihe Mate vector as the i n c r e m e n t a l unit. This feature is useful for c a l c u l a t i n g weighted spectra. F o r n o n - b r a n c h ing ,nstrnetions the p r o g r a m c o . m t e r is i n c r e m e n t e d not b y the subrrxitine which c o r r e s p o n d s to the instruction executed, bu.. r a t h e r by the controlling routine within M O N T H Y . ~y definition b r a n c h i n g instructions do c h a n g e Ihe value of the p, o g r a m c o u n t e r in the associa t e d s u b r o u t i n e . T h e sequence described above i'-; rep e a t e d for every event to be simulated, a n d finally the a r r a y s a n d spectra are written out on dl~k to, h~; acct'ssible by a display code. O n c e the desired pool of o p e r a t o r s has Been crc:~tctl for the p r o b l e m to be solved, the desired sec,.u,:nce ~f o p e r a t i o n s has to be described by a "program" for the virtual m a c h i n e . T h e p r o g r a m also has to initialize the p a r a m e t e r vectors needed a n d define the elements of the state vector. T o achieve a r e a s o n a b l y consistent f o r m of the p r o g r a m file. several rules have been f o r m u l a t e d . All o p e r a t o r s , which describe transitions of the state vector are assigned a n a m e which starts with the syrnbol > . w h e r e a s o p e r a t i o n s which are not p a r t of the instruction s e q u e n c e of the virtual m a c h i n e are described by names s t a r t i n g with a dot. The latter o p e r a t i o n s t~pically inc l u d e definitions of labels a n d files, a n d c o r r e s p o n d to n o n - e x e c u t a b l e s t a t e m e n t s for the vir:ual machine. All p a r a m e t e r s are read using the N A M E L I S T instruction in F O R T R A N _ If the n a m e s of the p a r a n t c t e r s are c h o s e n carefully, this will lead to a r,:.asonably m n e m o n i c p r o g r a m file w i l h o u t the need for including a sophistic a t e d parser. S o m e rules not to be discussed here govern the nam.;n6 of the pointers to the elements of the state vector_ F o r illustrative purposes a very simple input p r o g r a m file is s h o w n in fig. 4. It calculates the .,-- artd the .v-components of the velocity vectors of n e u t r o n s emitted isotropically f r o m a m o v i n g source with an e n e r g y s p e c t r u m given b y a Maxwellian. This example does n o t involve a n y b r a n c h i n g . At first the title " d e m o " of the o u t p u t file c o n t a i n i n g the sorted events in a 64 x 64 matrix and the labels for the axis are defined. The sequer,ce of > - o p e r a t o r s b r a c k e t e d b y the .do loop a n d the .endkx)p statements f o r m the actu',d p r o g r a m [or the virtual machine. It c o n t a i n s the o p e r a t o r s for the following six transitio,,s: 1_ > e v a p T h e m a g n i t u d e o f the velocity of the emitted n e u t r o n in the restframe or the
136
W. 146 Wdcke / Moate (.'ash, ~,mulatio. ~ d ~ sample
C C c c c c c .read
input
liJe
I , ) t MqNTI bY
Probleml nelltrons all .' I'V~]r'~l ~kt(:'] l ~ ' ~ t r ' , l ~ : : : . l ] J y tro[,. ,~ ln>Vll,~ I so~rce w h a t i s the v e l o c i t y di~trlbutic~n ,,f th,- l~etltl,J~lm i~ll t),( laboratory systems? I a s s i g , a name t,) the file c o n t a i n i n g the r e s u l t F name2D-'demo' Se~d I label s o ~ ele~_-nts of the state v e c t o r Jx~ $'_a~-~Is labl-'vx(neutro~) ' 'vylneutron) ' Send
; el,
Stxtle
.read . ~bels c c .do io :,
Sloop >trap
I calculage £ter=iO0000 I calculate
tO000 events Send Haxwelltan ,~vapOr4tion
spectrum
$evap opener=lO,opvel.ll,ipT--2 e m a x = 2 0 . , k l n d = l , l w a t t = O Send I c a l c u l a t e an ~soteop~e a n g u l a r d i s t r i b u t i o n $ i s o t r p i p l a n a - ] o p t e t a - 1 2 o p p h i - i 3 r a n q e ~ ] 6 0 Send
>[6o~r:' >uela~d
I a d d v e l o c l t y v e c t o r s o- ne~]trons and s o u r c e Sveladd x p v l ~ l l , i p v 2 = 1 4 , J p t h l ~ 2 , i p t h 2 - 1 ~ ipphil-3,ipphI2-O opv3-20 Opth3-21 opphl~-22 opv~3=l opvy3-2 opvz3-23 $(-nd I display o n - l l n e p r o g r e s ~ of C a l c u l a t i o n on C M T $spy m.~apy-1 m y s p y - 2 d x s w - l O 0 , dyst)y~lO0$c'~d I s o r t the eventg i n t o ~ array w i t h axes v x a n d v y Ssort2D maxms=I mx(1)=l my{1)~2 dx{l)=6.6 dy{])~6.6 Send
>~py >sort; -,next end]c~ go out2[,
$out2D
mout(1)=l
Send
.end
Fig. 4. t ' , ' m p l c i n p u t tile for M O N T H Y to c a l c u l a t e i . ~ t r o p i c e m i s s i o n of n e u t r t m ~ f r o m :. m o v i n ~ ~ourcc t'-oy e a c h of the lO 5 c'-'en[-'. the six : .stations f r o m > e v a p to > next are exoeut~l. The nameli~l ~.tructures $~tring - - - S e n d suppl:,,' the ncce~,sary p a r a m e l c r s for e a ( h o1: ::'alion.
nucleus is determined by a spectrum generator. An is,tropic in-plane angular distribution is modelled by a random number &enerator yielding uniformly distributed numbers in the range of 0 to 31M) (de-
2. > i~:trp
Krees ). 3. "-
:ladd
4. > . . y
The velocity vectors of the neulron ~nd of the movin 8 nucleus are added. The resulting vdocity vectors arc displayed for on-line monitoring on the screen or a C R T termin,al.
5. > .~)rt2D
The evenl is sorted int(, a 2-dimensional array, where Ihe axes corrc,,pond t() the x- and y-coml'~mcnts of the rely,city vCt.:t Or.
6. > next
G o back to the bcginnin~ o1" the lr,.')p fur the . c x t ev(nt. "I he numl',er'~ 1 -6 i . Ihc preceding list cot "c~pond t - the value the proEram counter pc ~,..ill altain at this operatiou. The .endloop statement tells M O N T H Y that the program for the virtual machine is n o w cor~91etely r..'ad in and ready to be executed by the .go st~ ' :- ent for 10 ~ events as has been specified by the _do :~x>p ,,tatcmenl. After the calculations are finished, the .orted data are written tm a disktile to h,e displayed la er, In fig. 5 is sho~,n the result of the calculation d~scribed above, which took 5 rain CPU-time on a V A X 750_
4. Conclmlons
r'~o --.
/os/
--.,~..~"~/_%.. " ~..~.--,~p"
~c~ ,:.,.'~
j"%
%
Fig. 5. I'~ result o[ ex,','utln 8 the r,amplc input hie shown in I'i~,. 4. " , h e ~ n p h de'picls the probability lot neutron e m i ~ i o n a~
a lunc'ion of the oentroo veto,.-'ity in two dimensions.
"] he structure of a flcxihle Monic C'arh~ c(xlc suilahlc for simulating the hmc evolution of a physical syst,:m i~ d i s u s e d . The code achieves a hish dc~r,'-" of generality by operating in the same way as a standard yon Neumann computer. The instruction set of th= virtual computer formed in the code is determined by the user by writing subroutines in F O R T R A N , Experience has
14/. 14'. Wdck¢ / Me*hie Carlo .limuiarion code ~ h o w n ta~,l "~fi,, ~tru,,.lurcd de,dgl allows, e:J',y m~hfic:~*
t~,,. :.J-d d d a p t i o n of die c ~ e to actual physical I';robI, n,.. h , f~err,on~ n,~! invol',,ed in the de',ign of the code. "the atltl,,i w - u h l likc to a c k n o w l e d g e very stimulati n g di,,zur,~,ion.,, with T)r,,. I.R. Birkelund. W.U_ Schr~der, J.R- l-luizensa a n d J.l'. K,:,,,ky ab,ou~ the idea of a general p u r p o ~ Mon'.¢ C.arlo t'r:,de ~tnd the help of Mr. W. Chert a n d Mr. L . . ~ h m i c c l e r in '.vritin~ p a r t s o f the d i s p l a y syst,-m D O I . L A R . This ~ o r k was !,"ipported b y Ihe I.JS D e p , . r m | e n t of Energy.
R~ferenc~ [l! ~, Ihudrr, MonLe C~lrIo T~le:h,~,,,ls in .~l~ki~,tical physics (Spr,,~,:r. ff.¢rlin, Iglq).
137
[2] N.P. Huslcnko, D.I. Gole'nko, Y A. Shrekler, IM. Sobol and V.G Sragovich. The ) ,onte C0r'.o method (Pergamon. Oxford, 1~66). [3] J.R. ('opley, Comp. Phys. C~m~m. 9 11975)64. [4} W.|!. Veseley et at., Report UC-32, Aer~jet Nuclear Company, Idaho Falls (1969)_ [5] E. Duck. L. Kowalskl and J.M_ Alexander. Grsay rep~)rl IPNO.DRE 82.20 (1982). [61 P. Dufour and J. Schlesinger, Comp. Phys. (_ ',,ram. 9 (11,175~ 360. [71 J M. Blatt and V.F. Weisskopf, 'l'he~r¢tical nuclear physit:s (Wiley, New York, 1982). [8] J.B. Birks, 'l]~e theor5' and practice or scintillation ,:ountmg (Macmillan, New York. 1964). i:;} Von Neumann, Nat. Bur. Stand. l',la',~.l.Serie~ 12 (1951) 36.
132
A GENL~AL Wmfriec
PURPOSE
~. WILCKE
CODE
FOR MONTE
CARLO
SIMULATIONS
*
Nuclear .(;t ": cture Research Lal~rato~ and Depaftmen~ of Chemis: O" and Ph)'~ic,~, C'nit~v,~i(t' . / R,*ch¢~ter. Ro~ h,'~ler. New York , t627, tJSA
Received ;:' )uly 1983
A gcn~: ,I-purpose computer code M O N T I I Y has h c c . w/httLL hJ t),t'lfl~rl|, ~'l,~l,te (~';irh, ',;mul;.,tit,lc- ~',f i*d=y :t.~l K~r'lLrll%. TO a c h _ v e a high d e g r ~ of heJtibility the c(~le i'- t)rgAnized liJ.e a Bcneral purptJ-,.¢ c o m p u t e r , oper.atin 8 t)n ,, ~.eclls! de',crih~n~ ~he time depei, lent slate o1 the syslc'm under s i m u l a : i o n , T h e instruc.ti3n ~ t <.; the " c o n t p u l e r " i', defined by it:' u~er and is therefore a d a p t a b l e . , the particular p r o b l e m studied. T h e ~ ' g a n i z a t i o n =)f M : ) N T t l Y allc,~.~, il©ralive and ¢onditio;'t;.:l ~.~ ~culi~'~n ~)1 ~)per;nit,n ~,
I. Inlrod" ~. lion
The [ : werful M o n t e Carlo methods have been applied to dye a wide range of physical problem,; [1.2]. One p a r ; ular area of success has been the modeling of systems ~y studying the time evolution of its c(.,nstitucnts A well-known example of t , i s kind is the design o' a nuclear reactor core by s i m u l a l i n g the fate of man~ hdividual neutrons as they are emitted, therrealized 'id eventually absorbed_ S o p h : ,teated special purpose codes have been puhlisbed t( ;olve such specific problems [3-5]. In general it is noZ ~articularly easy to adopt such codes to new problem, or change the models underlying the simulations [6 The code M O N T H Y discussed I'.ere i,'; an attempt t3 create a systematic approach to the simulation of ~ wide variety of physical processes with one commot approach. Since it is impossible to anticipate all possi : le models or problems a user of the MON"I'HY code m;-, want to simulate, M O N T H Y is not a readyt o - u ~ t,.de, but more akin to a computer language, which -~:cepts a c o m m a n d sequence s;.ored on disk (a "progr;.:S') describing the pl'ocess to be simulated in terms o: "instructions", which are to be defined by the u ~ r . q'L~.'se " i n s t r u c t i o n s " are introduced into the MONT'-IY code by adding user-supplied F O R T R A N subroul.~es to "he existing backbone structure of M O N T : I Y . A systematic set of simple rule,; d c ~ r i b e s how to :tdd these subrouhnes and how to prepare the c o m m a ~d files. T b e ~ rules arc described in the separate docum,;ttation of tl-~e code. In the present p;tper the more srmal concept underlying the s t t u c ' u r e ef M O N ' I ~'IY are to be discussed_ " Prose; ~ address: IBM, Thomas J. Wa:son Re.,-c.arch Lab:,ratory, ': ~)rktown Heights, New York 10598, USA. 0167-5( ;7./84/$03.00 © Elsevier Science Publishers B_V. (North i-!oiland Physics Publishing Division)
2. (hencr~l concepl.i ul M O N I ' t l Y
2.1. Introduction "Co illustrate Ihe ahstrat:t conccpt~ to be discussed below, lel us first ititr~duce ,a.,o specific example,, selected from applications of the co.de in nucle,~r phy~,ic,,. As the f~rst ez.ample we consider a highly excited nuch'u~ moving in the lai')orater3' system with a certain velc,c~ty t/ T,le nucleus :s reducing its hiL=h temperature I' by evaporating neutrons [7]. This process will continue unti; the excitation energy of the nucleus is I,:~., h~w Io .supr,ort further emissi(,n o f neutrons. Wc May ask f,:)r the ~:ngu[ar- and e n e r g / distrtbutiem (of the neutrons as obse:'vcd in the laboratory system. In a second example we want to study the prol)erlies of a large particle detector. Consider a tank filled with a .cintillating liquid, whicia is e x p o ~ d to a radiation field [8]. Several photomuh~plier tubes observe the interior of the tank through windows. W', want to calculate the probabilily th:n the photons created by ~ i n l i l l a l i o n s are actually de~ec'ed I'.y ope of the p h o t o m u h i p l i e r tubes, taking into account the absorption of the light in the liquid, multiple reflections o l the light from the walls and the efficiency of the phototubes. T h e ~ and m a n y other problems have in c,.~mm0n that it ;s possible IO de~crib¢ the Slate of Ihe system to be simulated by a suitably ch*~,,,en ~ t [,:) of N n u m - c r s x,. Ii is cc, r'=s'enient Io call this ~',et :.* "'*,late vector", even th(rugh it is not an element of Hilbert s p a t e and does not describe a system in the q u a n t u m m~h.',nical sen~,e. The emission of a neutron or the reflection of a photon from the wall of the rank can be considered as a transition betw,-.,'n two states of the ,system. We assign to each transition an operator Ox, which changes the
J¢'. H' ltgh.l-,, / M.nt¢ ('arlo itmu,"llttJn t'oth' H',: vvcl,W ,If the syslen~ from ,late ],~ ) to ,,tale ix'}:
; -'.., /6 :'-). ' l h e index K dvli.~tti~hes the kind of transition described hy an c,per;,t,,r. '[h¢ operators are not subject to any .~rmal restrictiol~s a . d m; 0 be non-linear functionals o f the s t a t e v e c t o r . COllt~liv ~toch;islic o p e r a t i o n s o r
wh;devcr else may be neces~,ary I~.) ~.intl,ki~c the system. I:ach transition may correspond to a dil!Cr,:~tkind of operation, or the same operator may he usud ~uF,c,ttedl2/ for different Iran.,,lliOllS ;In it is the case for s u b s c q u e l l l ClltiS'~ion OF severa] lleU[r~/ns f r o m a l:Ot nucleus. A gc,,eral ~lrillg o f ,'~ IransAiollS s l a v be writtcfl as
where the index KO) indicates what kind of" operation is used in the irh t~ansition. "lhe probleat is more complicated than discuss~:d so far, since the sequence of tr;m-itions is usually r:o~ predetermin.:d b y the structure ol the problem, but different for each event t t be simulated. In the second examplc cited a photon may hit a p h o t o c a t h o d e veo, soon. or the photon may undergo several reflections befo-e it is either detected or ab,,orbed_ It depends on the r a n d o m initial conditions of the first state and tile result', of subsequent stochastic ;ransiuons which sequence of operations is to be executed.
~i P~2;£ t': i(;:;
Ix (L:+2)> ~n÷2
1.13
Ihe haste concept of M O N T I I Y is to solve this control problem by structuring the code like a general purpo:.e computer. Each transition may be considered as the result of executing an instru,.:tion of a viwtual Von N e u m a n n machine, which operates on the stale vectors as data. Thus we may define conditional branchinl~ instructions, which allows the flow of operations to he altered as a function of the state vector (data) to be proce';sed. The set of all operators ( O~,. } corresponds Io the instruction set of the virtual nlachine, and each transition is the result of executing exactly one instructi,m of this set. In the physical code of M O N T I I Y each inslruction corresponds to a F O R T R A N subroutine Stq?p!i,:d I¢. the user. 2.2. Str, cim',l,, .l~erc~'lons In the following two seclions, st,>chastic and deterministic operators will be discussed_ Stochastic operators are normally us,:d to create ~,
sequence of random numbers with a spot;r:'.! distribution given by the predeternfined funct.:r, /,qq. winch describes the particular physical process to I~e simulated. The independent variable x is an element of tile stale vector. In the first example .x corresponds to the vel(~ity of a neutron and f(x) to a Ma~cwelliau distribution. The well known Von N e u m a n n method [91 is used to calculate the sequence of r a n d o m numbers with the desire:l spectral distribution. The meth~.~ is illustrated in fig. 2. If a function is bounded by the value y: in the imerval xm,,, ~< x ~< x.,~., then the Von N e u m a n n mc:thod can be applied. A random number generater is used to select uniformly distributed numbers in the intervals [xm, .. x . , . . ] and [y..,,.,y.,..], A pair (.x.,)b~' of r a n d o m numbers is used to test the validity of the inequali'.y/(:,:0) >_vo. If the inequality is true. tllen tiffs choice of x =.,c o is accepted, otherwise this trial is disc.',rded. It is plausible that this procedure ',,,'ill yield a
Ymox
• •"
~'~
[X (t~)>
iX (tt-I)>
t. o',; 0~ CALC~I.,' i lCa~;
F;g. 1. "lhe rclalion between die flow of calculations, stale vct:lor and Ihe Iran'.itkzn,, bcct~,'c~n sLalC$ iS illustrated. Each tr~ll,~it[on corre'.,p,cw,d, t,y the c.lecutiolt of ¢xacll_v one instruction of the "-irtual computer or. equivalenfl',, to the call to a
subrouline.
Ymi ~"// ~~ ///
Xmm
~,~2(x)xx~ Xmal
Fig. 2. lllustralion of the Von Neumann mezhod to gcncral¢ a sequence of random Lunll~rs with a spectral distribution given by .fix). (See t~xt for details.)
H ~. 14:
134
H/~h'k¢ / M,~nte Curb, tt,'4,1~ i,m
sr..quence .,f x-vahies with a p r o b a b i l i t y weighted by
fix). Thq Yon N e u m a n n a l g o r i t h m c a n ~ s i l y be exlended t~ ' u n c t i o n s o f several variables - a f-,ncfioa of m i n d e p r:dent variables needs Ihe g e n e r a t i o n of m + 1 r a n d o m ~"ambers f o r the test of o n e trial .sample. In general I ~:: functions te be simulated will be d e l x n d c n t o n o n e ( : m o r e p a r a m e t e r s . (The tet,.perature T of the n u c l e u s , vlitting a n e u t r o n is a p a r a m e t e r in the o r e r a lion of e : a p o r a t i n g a n e u t r o n . ) Since a p a r a m e t e r ,'nay c h a n g e ~! the result o f a transition (e.g. the nut leus reduces : ,: t e m p e r a t u r e by emitting a neutron), o n e has to be c~. :ful not to violate the c o n d i t i o n of the Von N e u m a n " a l g o r i t h m that max[f(xA]'cy,.,... This i~ illustratec ~y the d a s h e d c u r v e ]'2 in fig. 2, where for a n o t h e r , alue of the p a r a m e t e r tile m a x i m u m of the function ::xcecSs y . ~ . . The Von N e u m a n n a l g o r thin w o u l d y d a f i a t - t o p p e d d i s t r i b u t i o n for f 2 ( x ) , if y~.= were n o adjusted a c c o r d i n g l y . Since it is wasteft~l of c o m p u t e time to cl"~xtseym, unnecessarily large jr:it to be safe. is desirable to ~ t ) , ~ equal to the maxim:~m of the k ction, but this n u m b e r is not k n o w n a p r ori. "1"o ~ ¢e this p r o b l e m in an a p p r o x i m a t e w a y the f o l l o w i n , a l g o r i t h m has been a d o p t e d in M O N T H ' ( : A start val .~: for y m . is c h o s e n as v,~,~ = a - m a x ( / ' ) . ~ here a = 1.2 -: a safety f a c t o r a n d m a x ( / ) is c a l c a l a l e d for the initi , values or the p a r a m e t e r s . W h e n the p a r i m e ters c h a "~;e, the c o d e keeps track of the largest v a h e of fix) ever . n c o u n t e r e d a n d adjusts ).',~. u r ' v a r d s if n :cessary. F c " functions which are c o n t i n u o u s with r e s l ~ c l to all Ihei; p a r a m e t e r s this p r o c e d u r e will lead to very small en- )rs if a sufficiently large n u m b e r of events is sampler W h i l ' it w o u l d have been possiblc t~ store the p a r a m e . " -s as elemenls in the state "urine, it was f 3 u n d to be m ,:e c o n v e n i e n t to associate with e a c h i n s t r u c t i o n step it'. he " p r o g r a m " of the virlua', m a c h i n e a p a r a m e ter vect , which c o n t a i n s all p a r a m e t e r s for the c p e r a for c a l h , for b y this instruction. Obviously, e a c h pattie-
, ¢~d,,
u l a r o p e r a t i o n has its ov, n f~)rmat for Ihe inlerprelaticm of its ~,~,socialcd p a r a m e t e r vector_ In general the n u m e r i c t l values of the elements in the p a r a m e t e r vector will be different for re'pealed calls of the same o p e r a t o r W i t h the i n t r t ~ u c t i o n of the p a r a m e t e r vector r, the .,~quencc o f n ~ransilions c a n he '~.'rltten in it.,, final f o r m as
I., -t')> =
~^.,.,( t%.~)., b~,,,( P^-,,,) ~'~").
where e. ch of the indices 1 - - . n c o r r e s p o n d s 1o Ihe execulio:~ of o n e instruction in the virtual ntaehlnc. In the pr~.'~ding discussion a d i s t i n c t i o n has been m a d e bq:lween argumenls Of the ',lille vec[~); a n d p : l l a n w l e r s .
While tl is dlstinclion a p p e a r s q u i t e 3 a l u l a l l y in nto~,l real prol.lems, it is really s o m e w h a t : b i t r a r y T h e s a m e q u a n t i t y m a y be a i n d e p e n d e n t vat!':. ~le for one transition an(' a p a r a m e l e r in a n o t h e r t~::e, r)r a user m a y decide t')at he w o u l d r a t h e r not k ''.rp a q u a n t i t y curlslant, b~:t allow it t o fluctuate, hi o r d e r to .'qkps. for s u c h e x ( ' p t i o n s w i t h o u t the need (',f rewriting o p e r a t o r s ( s u b r o u L n e s ) , two special o p e r a t i o n s ha~,e been introd u c c d , v hich allow an e x c h a n g e of a r b i t r a r y elements in the slat( vector a n d the p a r a m e t e r vectors.
2-3. Dee rminLrttc operolor.s' Seve~ al deterministic o p e r a t i o n s arc generally ncce~,s a r y to solve a p r o b l e m . [hese o p e r a t i o m , yield ten o u t p u t .~tale vector_ svhich is u n i q u e l y dt:tcrmined b y the i n p u t st'de vector, the o p e r a t i o n executed a n d the ass()elated i: ~rameters. Typic:t] examples of these o p e r a t i o n s are kil~:malic;t] c:dculvlion~, an, I the ev:du:ttion o f n , a t h c n . a l i c a l functions. O t h e r o p c r a t i o w , i0 this clas~ tu,-'h,de b r a n c h i n g instuctions, the sorling of events into s p e c t r a a n d c o r r e l a t i o n a r r a y s a n d the g a t i n g of events. An ovcJ view of all presently defined classes of o p e r a t o r s ( i n c l u d i n g the stochastic ones) is given in table 1,
Table 1 Overvie~. of pre.s~ntly define.d class,'..s ol operators in MC N'I-HY. The clas,~ numb.crs refer Lo internal operation c~x.les and the colulrm labeled ~ gives the increment of the prol~rarn counter :3c for ol~ralo[s in each class Note thai I'~ rcfe[s to a p~,mter in the (tMe MON'I'] : f -,bhich selects the oI'~]ations !o be carried ola and not to t~.: hJrdwired prograln c.untcr r~l the .icm:d c,m,lmwr wlm'h execut~ MON'[ tlY. Class
N a m e of cla%
pc
.......................................................... "l')'pical c',amples
1000 4(}~)
Stochastic operators l')¢terminlstic o p e r a t o r s P a r a n ~ t e r ~ g,t a t e ) B r a n c h instructions E v e n t sorting Mi~ellaneous
4- 1 4- l + l vat. + 1 + 1
[;valc)oratlOll ~ F l r u [ l l , Gau~vJall ~|)¢~[ruT[] Kinematics, algebraic expre~slons Exchange elements in i-") and 9 Branch if x, - 0 Sorl events inlo a cot.clarion array Next event, gating cc..)dlti~n
7000 g0130 9000
~'~/ ~'~', t~"l/l #', I' J ,~/or//t I ('~{~/,~ .~l,llrld/'l:loll (',¢k(11'
ht~idVl,.~L.ntaiionof c,,d¢'
it,
A';. h;t> I,C~H ,li*;cussed above, the f u n d a m e n t a l ccltc c p t of M O N ' I {{I.Y i'. (o create a t'xw_~lof o p e r a t o r s (~b," w h i c h form the instrucli ;n ,~et of a virtual Von Nt.u, n a ' m m a c h i n e . In fig. J Ihc c,','crall o r g a n i z a t i o n of :his " ' m a c h i n e " is s h o w n . T h e plo~r:,m c,,tmter (pc) of the virtual "'rl{achine" is n o r m a l l y increl~.,,.mcd .by o n e after e a c h o p e r a t i o n a n d I'~)i its to the next "Jl~-,l!'=J,Jli(al. I]r:mch il131:lUCliOns c a n a,,',ign o t h e r value: to the 13rolgr;lll CtHlltter. As~,ik,:iated with e a c h inlruclioll is a paranlci._'r vector, w h i t h conluilLs the i n f o r n t a t i o n r e q u i r e d for the o p e r a t i o n t~ he cxecuied d u r i n g this i n s t r u c t i o n . T h e d , d a ,nemc~ry of Ihi.,, virtual m a c h i n e is very small a n d c o n t a i n s only the slate vector Ix.~, As has been discussed, there are special operatio,~s w h i c h allow exc h a n g e of elements of tile state vectt:r with elements of ,m a r b i t r a r y p a r a m e t e r vector. The i m p l e m e n t a t i o t | of the vir,ual m a c h i n e in a F ( ) R T R A N c~x.le is ~,traightfor~,,ard. A c o d e n u m b e r is assigned to e a c h kind of o p e r a t o r , wl,ich is stored in the " p r o g r a m memory"_ T h e p a r a m e t e r s are stored in a m a t r i x (actually two matrices, one for integer n u m b e r s a n d o n e for floating p o i n t n u m b e r s ) . E a c h row of these matrices c o r r e s p o n d s to a p a r a m e t e r vector, with the p r o g r a m c o u n t e r IX: p o i n t i n g to the next row to be a c c e s s e d U p o n execution of the p r o g r a m fo," the virtual rnachine Ihe p r o g r a m c o u n t e r is initialized to one. a n d the first i n s t r u c t i o n c o d e n u m b e r is fetched f r o m the p r o g r a m m e m o r y . In the actual F O R T R A N c o d e a c o m p u t e d G O T O s t a t e m e n t is use~J to call the first
C
CODt" Op~c~c' in= OK{=)
P};[l)
PIK ()1
3
Opce~e ~o[ OX(])
PK{3)
(::~,tr'-x ~ _'atr'.×)
..L2"._':'.~:L~_JAL'C_ ! x > ~- l x i ,
x2, . . . . . . . . . . . . .
LN}
IzxP)
Filz. 3. The or~,amzation of the MO/ITI-IY ctw.l¢ follows the ,,Iructure of a general p u r ~ cc,mpul4 r with separate program Inemo[y
a=id d a l a lUelTIOly_ "1 h e d a t a n : t m u r y
contains
only the
pre.wrnt slate vector ~.,¢>,whereas the pxogram memory_ contains the int,.lruction co,de and the pa,am¢l¢/' vector.'=.
] 35
operator, which is represeqted by a subroutine of M O N T I I Y . After the return from this s u b r o u t i n e (i.e_ execution of an instruction in the terms of tile ",'irtuai m a c h i n e ) the p r o g r a m c,..unter pc .)f the virtual m a c h i n e is i n c r e m e n t e d b y o n e (except for b r a n c h i n g operations), a n d the cycle is repeatec, for the next operation. Somew h e r e d u r i n g the sequence of instructions the d a t a in the state vector are used as addresses to i n c r e m e n t the a p p r o p r i a t e spectra a n d c o r r e l a t i o n arrays. By default every event i n c r e m e n t s the selected channel by unity. b u t an option exists to use a n a r b i t r a r y element of Ihe Mate vector as the i n c r e m e n t a l unit. This feature is useful for c a l c u l a t i n g weighted spectra. F o r n o n - b r a n c h ing ,nstrnetions the p r o g r a m c o . m t e r is i n c r e m e n t e d not b y the subrrxitine which c o r r e s p o n d s to the instruction executed, bu.. r a t h e r by the controlling routine within M O N T H Y . ~y definition b r a n c h i n g instructions do c h a n g e Ihe value of the p, o g r a m c o u n t e r in the associa t e d s u b r o u t i n e . T h e sequence described above i'-; rep e a t e d for every event to be simulated, a n d finally the a r r a y s a n d spectra are written out on dl~k to, h~; acct'ssible by a display code. O n c e the desired pool of o p e r a t o r s has Been crc:~tctl for the p r o b l e m to be solved, the desired sec,.u,:nce ~f o p e r a t i o n s has to be described by a "program" for the virtual m a c h i n e . T h e p r o g r a m also has to initialize the p a r a m e t e r vectors needed a n d define the elements of the state vector. T o achieve a r e a s o n a b l y consistent f o r m of the p r o g r a m file. several rules have been f o r m u l a t e d . All o p e r a t o r s , which describe transitions of the state vector are assigned a n a m e which starts with the syrnbol > . w h e r e a s o p e r a t i o n s which are not p a r t of the instruction s e q u e n c e of the virtual m a c h i n e are described by names s t a r t i n g with a dot. The latter o p e r a t i o n s t~pically inc l u d e definitions of labels a n d files, a n d c o r r e s p o n d to n o n - e x e c u t a b l e s t a t e m e n t s for the vir:ual machine. All p a r a m e t e r s are read using the N A M E L I S T instruction in F O R T R A N _ If the n a m e s of the p a r a n t c t e r s are c h o s e n carefully, this will lead to a r,:.asonably m n e m o n i c p r o g r a m file w i l h o u t the need for including a sophistic a t e d parser. S o m e rules not to be discussed here govern the nam.;n6 of the pointers to the elements of the state vector_ F o r illustrative purposes a very simple input p r o g r a m file is s h o w n in fig. 4. It calculates the .,-- artd the .v-components of the velocity vectors of n e u t r o n s emitted isotropically f r o m a m o v i n g source with an e n e r g y s p e c t r u m given b y a Maxwellian. This example does n o t involve a n y b r a n c h i n g . At first the title " d e m o " of the o u t p u t file c o n t a i n i n g the sorted events in a 64 x 64 matrix and the labels for the axis are defined. The sequer,ce of > - o p e r a t o r s b r a c k e t e d b y the .do loop a n d the .endkx)p statements f o r m the actu',d p r o g r a m [or the virtual machine. It c o n t a i n s the o p e r a t o r s for the following six transitio,,s: 1_ > e v a p T h e m a g n i t u d e o f the velocity of the emitted n e u t r o n in the restframe or the
136
W. 146 Wdcke / Moate (.'ash, ~,mulatio. ~ d ~ sample
C C c c c c c .read
input
liJe
I , ) t MqNTI bY
Probleml nelltrons all .' I'V~]r'~l ~kt(:'] l ~ ' ~ t r ' , l ~ : : : . l ] J y tro[,. ,~ ln>Vll,~ I so~rce w h a t i s the v e l o c i t y di~trlbutic~n ,,f th,- l~etltl,J~lm i~ll t),( laboratory systems? I a s s i g , a name t,) the file c o n t a i n i n g the r e s u l t F name2D-'demo' Se~d I label s o ~ ele~_-nts of the state v e c t o r Jx~ $'_a~-~Is labl-'vx(neutro~) ' 'vylneutron) ' Send
; el,
Stxtle
.read . ~bels c c .do io :,
Sloop >trap
I calculage £ter=iO0000 I calculate
tO000 events Send Haxwelltan ,~vapOr4tion
spectrum
$evap opener=lO,opvel.ll,ipT--2 e m a x = 2 0 . , k l n d = l , l w a t t = O Send I c a l c u l a t e an ~soteop~e a n g u l a r d i s t r i b u t i o n $ i s o t r p i p l a n a - ] o p t e t a - 1 2 o p p h i - i 3 r a n q e ~ ] 6 0 Send
>[6o~r:' >uela~d
I a d d v e l o c l t y v e c t o r s o- ne~]trons and s o u r c e Sveladd x p v l ~ l l , i p v 2 = 1 4 , J p t h l ~ 2 , i p t h 2 - 1 ~ ipphil-3,ipphI2-O opv3-20 Opth3-21 opphl~-22 opv~3=l opvy3-2 opvz3-23 $(-nd I display o n - l l n e p r o g r e s ~ of C a l c u l a t i o n on C M T $spy m.~apy-1 m y s p y - 2 d x s w - l O 0 , dyst)y~lO0$c'~d I s o r t the eventg i n t o ~ array w i t h axes v x a n d v y Ssort2D maxms=I mx(1)=l my{1)~2 dx{l)=6.6 dy{])~6.6 Send
>~py >sort; -,next end]c~ go out2[,
$out2D
mout(1)=l
Send
.end
Fig. 4. t ' , ' m p l c i n p u t tile for M O N T H Y to c a l c u l a t e i . ~ t r o p i c e m i s s i o n of n e u t r t m ~ f r o m :. m o v i n ~ ~ourcc t'-oy e a c h of the lO 5 c'-'en[-'. the six : .stations f r o m > e v a p to > next are exoeut~l. The nameli~l ~.tructures $~tring - - - S e n d suppl:,,' the ncce~,sary p a r a m e l c r s for e a ( h o1: ::'alion.
nucleus is determined by a spectrum generator. An is,tropic in-plane angular distribution is modelled by a random number &enerator yielding uniformly distributed numbers in the range of 0 to 31M) (de-
2. > i~:trp
Krees ). 3. "-
:ladd
4. > . . y
The velocity vectors of the neulron ~nd of the movin 8 nucleus are added. The resulting vdocity vectors arc displayed for on-line monitoring on the screen or a C R T termin,al.
5. > .~)rt2D
The evenl is sorted int(, a 2-dimensional array, where Ihe axes corrc,,pond t() the x- and y-coml'~mcnts of the rely,city vCt.:t Or.
6. > next
G o back to the bcginnin~ o1" the lr,.')p fur the . c x t ev(nt. "I he numl',er'~ 1 -6 i . Ihc preceding list cot "c~pond t - the value the proEram counter pc ~,..ill altain at this operatiou. The .endloop statement tells M O N T H Y that the program for the virtual machine is n o w cor~91etely r..'ad in and ready to be executed by the .go st~ ' :- ent for 10 ~ events as has been specified by the _do :~x>p ,,tatcmenl. After the calculations are finished, the .orted data are written tm a disktile to h,e displayed la er, In fig. 5 is sho~,n the result of the calculation d~scribed above, which took 5 rain CPU-time on a V A X 750_
4. Conclmlons
r'~o --.
/os/
--.,~..~"~/_%.. " ~..~.--,~p"
~c~ ,:.,.'~
j"%
%
Fig. 5. I'~ result o[ ex,','utln 8 the r,amplc input hie shown in I'i~,. 4. " , h e ~ n p h de'picls the probability lot neutron e m i ~ i o n a~
a lunc'ion of the oentroo veto,.-'ity in two dimensions.
"] he structure of a flcxihle Monic C'arh~ c(xlc suilahlc for simulating the hmc evolution of a physical syst,:m i~ d i s u s e d . The code achieves a hish dc~r,'-" of generality by operating in the same way as a standard yon Neumann computer. The instruction set of th= virtual computer formed in the code is determined by the user by writing subroutines in F O R T R A N , Experience has
14/. 14'. Wdck¢ / Me*hie Carlo .limuiarion code ~ h o w n ta~,l "~fi,, ~tru,,.lurcd de,dgl allows, e:J',y m~hfic:~*
t~,,. :.J-d d d a p t i o n of die c ~ e to actual physical I';robI, n,.. h , f~err,on~ n,~! invol',,ed in the de',ign of the code. "the atltl,,i w - u h l likc to a c k n o w l e d g e very stimulati n g di,,zur,~,ion.,, with T)r,,. I.R. Birkelund. W.U_ Schr~der, J.R- l-luizensa a n d J.l'. K,:,,,ky ab,ou~ the idea of a general p u r p o ~ Mon'.¢ C.arlo t'r:,de ~tnd the help of Mr. W. Chert a n d Mr. L . . ~ h m i c c l e r in '.vritin~ p a r t s o f the d i s p l a y syst,-m D O I . L A R . This ~ o r k was !,"ipported b y Ihe I.JS D e p , . r m | e n t of Energy.
R~ferenc~ [l! ~, Ihudrr, MonLe C~lrIo T~le:h,~,,,ls in .~l~ki~,tical physics (Spr,,~,:r. ff.¢rlin, Iglq).
137
[2] N.P. Huslcnko, D.I. Gole'nko, Y A. Shrekler, IM. Sobol and V.G Sragovich. The ) ,onte C0r'.o method (Pergamon. Oxford, 1~66). [3] J.R. ('opley, Comp. Phys. C~m~m. 9 11975)64. [4} W.|!. Veseley et at., Report UC-32, Aer~jet Nuclear Company, Idaho Falls (1969)_ [5] E. Duck. L. Kowalskl and J.M_ Alexander. Grsay rep~)rl IPNO.DRE 82.20 (1982). [61 P. Dufour and J. Schlesinger, Comp. Phys. (_ ',,ram. 9 (11,175~ 360. [71 J M. Blatt and V.F. Weisskopf, 'l'he~r¢tical nuclear physit:s (Wiley, New York, 1982). [8] J.B. Birks, 'l]~e theor5' and practice or scintillation ,:ountmg (Macmillan, New York. 1964). i:;} Von Neumann, Nat. Bur. Stand. l',la',~.l.Serie~ 12 (1951) 36.