BINARY COLLISION APPROXIMATION CASCADE PROGRAM CONSTRUCTION

11

11.1

BINARY COLLISION APPROXIMATION CASCADE PROGRAM CONSTRUCTION

Introduction

This chapter describes the construction of a binary collision approxima­ tion c o m p u t e r p r o g r a m for the simulation of collision cascades. A collision cascade is a branching sequence of atomic collisions that is induced by the ejection of a P K A . T h e P K A is ejected by a primary radiation particle ( P R P ) . T h e radiation-induced defect states of greatest technological i m p o r t a n c e in iron and nickel alloys, for example, a p p e a r to be those p r o d u c e d by collision cascades initiated by P K A with energies above 3 k e V . T h e s e defect states (displacement spikes) are sufficiently complex that any technologically useful description of their structure and gross defect content probably can be obtained only from c o m p u t e r experiments. A t least that has b e e n the case so far. T h e only m e t h o d currently available for simulating large cascades is the binary collision approximation m e t h o d . It is c o m m o n practice to refer to a collision cascade in terms of the energy of the P K A which initiated the cascade. H e n c e , the term " 2 0 k e V cascade" signifies a cascade initiated by a 2 0 k e V P K A . Binary collision approximation ( B C A ) computer experiments were first r e p o r t e d by Yoshida [1961]. H e treated cascades in an a m o r p h o u s model of 4 g e r m a n i u m . A m a p of a 1 0 e V cascade in g e r m a n i u m , given by Yoshida's simulation, appears in fig. 11.1. This m a p clearly illustrates the sequential branching aspect of cascade evolution. B C A c o m p u t e r simulation of collision cascades in crystals was first r e p o r t e d by Beeler and Besco [1962,1963]. Their initial work concerned the p e n e t r a t i o n of iodine and krypton ions in B e O , which has the wurtzite structure. A cascade initiated by a 9 k e V Be P K A in B e O is shown in fig. 11.2, and a d a m a g e region m a p for this cascade appears in fig. 11.3. T h e Κ A trajectory m a p in fig. 11.2 shows only the trajectories directly associated with the 9 k e V B e P K A and the two most energetic S K A it p r o d u c e d . These trajectories reveal the central structure of the cascade, unobscured by peripheral trajectory lines. T h e d a m a g e m a p in fig. 11.3 indicates the production of vacancy strings as the initial vacancy defect structure. This work on B e O was followed by a B C A simulation study on cascades in bcc 649

650

Binary Collision

Approximation

[§11.1

Fig. 11.1 Vacancies (dots) and interstitials (open circles) produced by a lOkeV collision cascade in an amorphous model of germanium (Yoshida [1961]).

iron(m2), c o p p e r ( m ) and tungsten(m) (Beeler [1966]). A K A trajectory m a p for a 5 k e V cascade in bcc iron(m2), from this work, appears in fig. 11.4. This figure shows trajectories for all K A whose energies exceeded 8 e V and thereby illustrates b o t h the central and peripheral trajectory structures for a cascade in a crystalline solid. T h e objective of the work by Beeler and Besco was to describe the defect state produced by an isolated cascade at the atomistic level. They used the C A S C A D E Program (Besco and Beeler [1963]) to generate a collision cascade, and the C L U S T E R P r o g r a m , written by B a u m g a r d t , to analyze the defect cluster populations and deployment in the defect state (displacement spike) produced by the cascade. These two programs subsequently were combined into the C A S C A D E - C L U S T E R P r o g r a m by Besco and B a u m g a r d t [1965]. C A S C A D E - C L U S T E R was designed to accept an arbitrary initial defect state in which a new cascade was generated. This feature m a d e it possible to examine the defect state p r o d u c e d by a succession of different cascades in a given local region (Beeler

§11.1]

11.1

Introduction

651

Fig. 11.2. Trajectory map for a 9 k e V collision cascade in a crystalline model of B e O . The initiating P K A was a B e atom (Beeler [1971]).

|(X)

Τ

40 &

>

DAMAGE REGION (Z=47-50Ä) 9 keV Be PKA IN BeO

Fig. 11.3. D e p l o y m e n t of vacancies ( • ] in the damaged region produced by the 9 k e V cascade shown in fig. 11.2 (Beeler [1971]).

[1967]). L e k k e r k e r et al. [1971] developed a B C A cascade p r o g r a m m o d e l e d after the C A S C A D E P r o g r a m . Their p r o g r a m was n a m e d SIC. Both C A S C A D E and SIC assume elastic atomic collisions; they do not account for the effects of inelastic atomic collisions.

652

Binary Collision

Approximation

[§11.1

Fig. 11.4. Trajectory map for a 5 k e V collision cascade in bcc iron(m2). The short, thick line represents the P K A trajectory. The three main S K A trajectories are represented by heavy dotted lines. For higher order Κ A the trajectory lines are alternately solid and dashed lines (Beeler [1966]).

R o b i n s o n and T o r r e n s [1974] wrote a B C A collision cascade simulation p r o g r a m called M A R L O W E . They used the Firsov [1959] model for simulating the effect of inelastic atomic collisions and w o r k e d on an improved m e t h o d for describing multiple collision events which are characteristic of either focussed collision chains or channeling. A B C A p r o g r a m , called C O L L I D E , was written by Beeler and Beeler [1973]. This p r o g r a m uses an a m o r p h o u s solid model and was written to describe d a m a g e energy transport in M e V energy range self-ion b o m b a r d m e n t of austenitic steels and nickel alloys. C O L L I D E also uses the Firsov inelastic collision model. T h e r e are two main types of B C A collision models for cascade simulation: (1) the simple B C A m o d e l , applicable to a m o r p h o u s solids; and (2) the modified B C A m o d e l which should be used for a crystalline solid. In the simple B C A m o d e l , only o n e target a t o m is considered p e r collision. In the modified B C A m o d e l , simultaneous collision events between the projectile a t o m and several target atoms are considered when special collision geometries arise, such as those in focussed collision chains and channeling. Besco and Beeler [1963] accounted for multiple collision events in the

§11.2]

11.2 General features

of BCA computer

program

construction

653

C A S C A D E P r o g r a m using the forbidden target algorithm. In this scheme, a selected target a t o m for a particular K A cannot be any one of the four previous target atoms for that K A . This implicit m e t h o d treats multiple target instances on a sequential basis. Experience shows that the forbidden target m e t h o d is sufficiently general to simulate focussed replacement collision chains and long-range channeling in cubic and hexagonal struc­ tures. R o b i n s o n and T o r r e n s [1974] used an explicit scheme to model multiple collisions. In their algorithm, all target atoms whose distances from the projectile were less than a prescribed increment D R were considered as simultaneous targets. E a c h of the B C A schemes used thus far to model multiple target collision events has serious problems at low K A energies. T h e most reliable practical a p p r o a c h for describing low-energy, multiple target collision events, at this t i m e , appears to be the use of cataloged dynamical m e t h o d results.

11.2

General features of BCA computer program construction

A d v a n c e s in cascade simulation technique are n e e d e d to reduce the cost of cascade simulation for cascade energies above the equivalent of 2 k e V in tungsten. T h e dynamical m e t h o d can be used up to that level of cascade energy at an acceptable cost. But beyond that level of cascade energy, the cost of the dynamical m e t h o d is not acceptable. Because of this circum­ stance, it is worthwhile to examine the construction and working action of a cascade simulation p r o g r a m in detail, in order to clearly point out the specific technique i m p r o v e m e n t s that are n e e d e d . A B C A p r o g r a m which uses the forbidden target algorithm to simulate multiple collision events will b e described in detail. All of the major points of interest are illustrated by this example. A simple B C A cascade simulation program generates a continually branching sequence of two-atom collisions in each of which a moving projectile a t o m collides with an initially stationary target a t o m . T h e cascade is initiated by the ejection of a P K A from a normal atom site by a P R P . This P K A is the first projectile atom in the simulation. It subsequently collides with the first target atom and, in general, both the projectile and target a t o m e m e r g e from the collision with sufficient energy that each of t h e m induces a new collision, i.e., each of them serves as the projectile in a subsequent twoa t o m collision. A branching collision chain process, diagrammed schemati­ cally in fig. 11.5, develops in which the projectile population doubles at each generation. In fig. 11.5, " 0 " represents the initial P K A ejection event and " 1 " represents the first collision in the cascade. T h e first collision constitutes proliferation m o d e l , i.e., two projectiles per collision. T h e m a x i m u m

654

Binary Collision

g= I

g=2

Approximation

g=3

[§11.2

g=4

<2>

CASCADE

BRANCHING

Fig. 11.5. Schematic diagram of a branching collision chain for a binary collision model. Notation is explained in the text.

generation g = 1. T h e second generation, g = 2, consists of collisions 2 and 3. G e n e r a t i o n s g = 3 and g = 4 also are depicted in t h e figure. Given a collision n u m b e r , c, E(c,l) is the largest and E(c,2) is the smallest emergent a t o m energy for that collision. Em(c) is the inelastic energy loss due to electronic excitation during the collision. Only o n e collision event can be handled at a time in a B C A cascade simulation. Because of this, the energy, ejection point vector and ejection velocity vector for each additional projectile created must be held in storage until its turn occurs in t h e cascade collision sequence. This projectile

§11.2]

11.2 General features ofBCA

computer program

construction

655

candidate storage block is called the Q U E U E table. If two projectiles always e m e r g e d from each collision, it would be necessary to store two projectile candidates in the Q U E U E table after each collision. It can be seen from fig. 11.5 that the total n u m b e r of collisions, N C O L ( g ) , completed at the e n d of t h e gth generation would b e N C O L ( g ) = (2* - 1) and that the total n u m b e r of projectiles N Q ( g ) stored in the Q U E U E table at the end of the gth generation would be N Q ( g ) = V. This proliferation of projectile candidates is o n e of the critical problems in a B C A cascade simulation. Projectile candidate proliferation in a cascade does not continue indefi­ nitely, of course. A t each collision, the projectile kinetic energy is subdivided into t h r e e parts: (1) target a t o m kinetic energy; (2) a reduced projectile a t o m kinetic energy; and (3) electronic excitation energy. A s a consequence of this progressive energy subdivision, either o n e or both collision partners eventually will fail to induce a new collision and, as a c o n s e q u e n c e , the cascade eventually dies out. T h e failure to induce a new collision occurs because a target atom must receive a kinetic energy greater than the displacement energy threshold in order to join in the cascade (see ch. 10), and the projectile a t o m must maintain a kinetic energy that at least exceeds the sublimation energy in o r d e r to remain in the cascade. T h e s e effects are illustrated in fig. 11.6, which is a collision sequence diagram for a 2 0 0 e V cascade in fee iron(m2). In the particular simulation concerned, both t h e projectile and target a t o m cutoff energies were assigned the value 25 e V . N o t e that the energy of o n e of the collision p a r t n e r s , after collision, often is less than the cutoff energy. T h e underlined n u m b e r at each collision symbol is the collision n u m b e r . T h e Κ A with the largest energy was selected as the projectile a t o m for each collision. A n y ΚA with energy ^ 25 e V was used as a projectile a t o m . A n y K A with energy < 2 5 e V was d r o p p e d from the cascade. A vacancy ( • ) was formed when the post-collision energy of each collision p a r t n e r was ^ 2 5 e V . A n interstitial ( · - · ) was formed when the post-collision energy of each collision p a r t n e r was < 25 e V . T h e rise and fall of the projectile population in the Q U E U E table of the C O L L I D E P r o g r a m is illustrated in fig. 11.7 for 2 , 4 , 6 and 8 k e V cascades in fee iron(m2). In these particular instances, the projectile selection rule was to select the projectile that had the largest energy. A sharp displacement energy threshold Ed = 2 0 e V was assumed. H e n c e , each K A history was t e r m i n a t e d when the K A energy fell below 2 0 e V . T h e Q U E U E table population N Q initially rises along the line

•

Vacancy φ-φ

Split

Interstitial

Fig. 11.6. Schematic diagram of a 2 0 0 e V collision cascade in fee iron(m2). Underscored numbers 0 , 1 , 2 , etc. are collision numbers. Numbers not underscored are K A energies in e V . In the model selected for this simulation, an interstitial was formed w h e n the energies of the projectile and target atoms after collision both were < 25 e V .

Primary R a d i a t i o n Particle

§11.2]

11.2 General features ofBCA

computer program

construction

657

140

Fig. 11.7. Projectile populations ( N Q ) in the Q U E U E table of the C O L L I D E Program for 2, 4, 6 and 8 k e V cascade simulations in fee iron(m2). N C O L is the collision number.

NQ = 0.62*NCOL+1. This initial rise is less than the N Q = N C O L 4- 1 behavior of the unchecked proliferation m o d e l , i.e., two projectiles p e r collision. T h e m a x i m u m Q U E U E table population occurred at NCOL = NCOLT/2, w h e r e N C O L T is t h e total n u m b e r of collisions in the cascade. Table 11.1 gives some additional information. It shows that the n u m b e r of collisions in these cascades increased at the rate of about 64 collisions per k e V of P K A energy and that N Q m a x increased by 17 entries p e r k e V of P K A energy. B o t h N C O L T and N Q m a x increase markedly with decreasing K A history termination energy a n d , therefore, so does the computation time. For this

Binary Collision

658

[§11.2

Approximation

Table 11.1. Q U E U E table characteristics as a function of cascade energy E. These characteristics pertain to the C O L L I D E Program and cascades in fee i r o n ( m l ) .

E,keV

NC0LT

a

NQ

b r

r

/E

NCOLT/NQ

17.0

3.79

16.0

4.02

97

16.2

3.97

134

63.9

16.8

3.81

129

3^

257

64

6

385

8

511

645

NQ

64.5 64.2 64.2

2 4

10

NCOLT/E

I63

64.5

16.2

3.98

20

1272

316

63.6

15.8

4.02

30

1876

451

62.5

15.0

4.16

^ o t a l number of c o l l i s i o n s i n t h e cascade Maximum number of e n t r i e s o c c u r r i n g i n t h e QUEUE t a b l e reason, it is economically advantageous to use a ΚA history termination scheme that allows history termination at a relatively high energy. This m a t t e r is discussed in §11.11. A flow chart of the computational activity in a B C A cascade simulation appears in fig. 11.8. A cascade simulation begins with a reading of the input data which defines what is to be simulated and the particular way in which the simulation is to proceed. T h e input data normally submitted gives the following information: (1) crystal structure (bcc, fee, h e p ) , (2) length unit (hlc for bcc and fee structure; the basal plane lattice constant for hep structure), (3) box size (see §11.3), (4) crystal size (see §11.3), (5) a t o m mass and atomic n u m b e r , (6) P K A energy, E P , (7) P K A ejection direction, ( C A P , C B P , C G P ) , (8) P K A ejection point, ( Χ Ρ , Υ Ρ , Ζ Ρ ) , (9) potential function p a r a m e t e r s , (10) defect production criteria. T h e subroutine that performs this j o b will be called R E A D . After the input data has b e e n read, the program utilizes it to set u p the conditions u n d e r which the simulation is to be run. This involves setting u p t h e specified crystal structure and size, and initializing all table values and computational p a r a m e t e r s . It includes, as well, assigning the P K A energy,

READ (1)

659

Read i n p u t data

SETUP (2)

1 . Set up c r y s t a l geometry 2 . I n i t i a l i z e a l l t a b l e values 3 . I n i t i a l i z e a l l c o m p u t a t i o n a l parameters 4 . Assign PKA ( f i r s t p r o j e c t i l e ) e j e c t i o n p o i n t 5 . A s s i g n PKA e j e c t i o n d i r e c t i o n 6 . Assign PKA energy

TARSEL (3)

1 . Make a l i s t o f p o s s i b l e t a r g e t atoms 2 . S e l e c t t h e a p p r o p r i a t e t a r g e t atom(s) from t h i s l i s t

KOLIDE (4)

Compute the f o l l o w i n g c o l l i s i o n r e s u l t s : 1 . P r o j e c t i l e atom e x i t e n e r g y , EP 2 . P r o j e c t i l e atom e x i t d i r e c t i o n (CAP,CBP,CGP) 3 . P r o j e c t i l e atom e x i t p o i n t (ΧΡ,ΥΡ,ΖΡ) 4 . T a r g e t atom e x i t e n e r g y , ET 5 . T a r g e t atom e x i t d i r e c t i o n (CAT,CBT,CGT) 6 . T a r g e t atom e x i t p o i n t (ΧΤ,ΥΤ,ΖΤ) 7. I n e l a s t i c energy l o s s , EINE 8 . C o l l i s i o n t i m e , TC

BUKEEP (5)

Assay c o l l i s i o n r e s u l t s : 1 . Score vacancy p r o d u c t i o n 2 . Score i n t e r s t i t i a l p r o d u c t i o n 3 . S t o r e vacancy c o o r d i n a t e s i n d e f e c t t a b l e 4 . Store i n t e r s t i t i a l coordinates i n d e f e c t t a b l e 5. Enter a p p r o p r i a t e KA atom(s) i n t h e QUEUE table

STOP TEST

Are t h e r e any s u f f i c i e n t l y e n e r g e t i c a v a i l a b l e t o c o n t i n u e t h e cascade?

(7)

projectiles

(6)

ASSIGN (7)

A s s i g n c o l l i s i o n parameters f o r t h e n e x t Projectile

PRINT (8)

1 . P r i n t cascade r e s u l t s 2 . P l o t t r a j e c t o r y map 3 . P l o t d e f e c t p o s i t i o n map

STOP

End o f s i m u l a t i o n

•*(3)

(6)

Fig. 1 1 . 8 . Flow chart for a B C A cascade simulation program.

660

Binary Collision

Approximation

[§11.12

ejection direction and ejection point coordinates values to the first projectile a t o m . T h e subroutine that does this job will be called S E T U P . T h e target a t o m which the projectile atom collides with is selected from a list of possible target candidates. T h e particular target atom selected d e p e n d s u p o n the projectile takeoff (ejection) point, the projectile takeoff (ejection) direction, and the center-to-center distance b e t w e e n the p r o ­ jectile atom and the target atom candidate. T h e r e are two principal target a t o m selection m e t h o d s : (1) the cell m e t h o d , and (2) the neighbor shell m e t h o d . In the cell m e t h o d , a cubic cell is constructed a r o u n d the projectile a t o m position such that each corner of the cell corresponds to either a n o r m a l a t o m site or an octahedral interstitial site. A list of possible target atoms can then be m a d e on the basis of the cell structure. In the neighbor shell m e t h o d , the list of possible target atoms is t a k e n to be all first and second neighbors of the a t o m site nearest to the projectile a t o m . T h e target selection subroutine will be called T A R S E L . Given the projectile atom position, energy and direction, and the target a t o m position, then the collision angle, collision time and energy assign­ ments for each of the collision partners can be calculated by numerical integration of the classical equations of motion for the collision. These quantities are then used to calculate the following collision exit items: (1) (2) (3) (4) (5) (6) (7) (8)

projectile atom exit energy, E P , projectile atom exit direction, ( C A P , C B P , C G P ) , projectile a t o m exit point, ( Χ Ρ , Υ Ρ , Ζ Ρ ) , target a t o m exit energy, E T , target a t o m exit direction, ( C A T , C B T , C G T ) , target a t o m exit point, ( Χ Τ , Υ Τ , Ζ Τ ) , inelastic energy loss, E I N E , collision time, T C .

T h e subroutine that does these calculations will be called K O L I D E . After a collision is completed, it is necessary to d e t e r m i n e w h e t h e r a vacancy, an interstitial or both were created by the collision, and at what location each of these elemental defects was created. In addition, it is necessary to record the location of each defect produced in the defect table and to add projectiles to the Q U E U E table. T h e subroutine that does these b o o k k e e p i n g chores will be n a m e d B U K E E P . A cascade termination test must be m a d e after each collision. This is a simple test: if there are no projectiles available to continue the cascade, the simulation is stopped. T h e termination test normally follows the b o o k k e e p ­ ing operations. A s long as projectiles are available, a subroutine called A S S I G N selects the projectile for the next collision event.

§11.3]

11.3 Box method for locating defects and C

CASCADE SIMULATION PROGRAM

atoms

661

Text

Section

C CALL READ CALL SETUP

11.3

1 0 CALL TARSEL

11.4

CALL KOLIDE

11.5

CALL BUKEEP

11.10,

15

IF(NQ.EQ.O)

-

11.9 11.11

GO TO 2 0

CALL ASSIGN

11.10

GO TO 1 0 2 0 CALL PRINT CALL EXIT END

Fig. 11.9. Cascade simulation program. The text section number in the right-hand column gives the section or sections in which the associated subroutine is described.

W h e n the cascade is completed, the P R I N T subroutine prints the cascade collision history and the lists of vacancy and interstitial positions for the primary defect state p r o d u c e d . A workable cascade simulation p r o g r a m can be outlined succinctly in t e r m s of the seven subroutines we have just discussed. This p r o g r a m is listed in fig. 11.9. In statement 15, N Q is the Q U E U E table population, i.e. the n u m b e r of available projectile atoms. W h e n N Q = 0, no projectiles remain to p e r p e t u a t e the cascade. A detailed account of each subroutine in this p r o g r a m is given in the succeeding sections of this chapter. T h e pertinent section n u m b e r s are cited in parentheses at the right of each C A L L s t a t e m e n t in fig. 11.9.

11.3

11.3.1

Box method for locating defects and atoms

Introduction

A scheme for specifying the locations of normal atoms and elemental defects in a crystal is the foundation for any cascade simulation p r o g r a m . T h e so-called Box M e t h o d is a general scheme for locating normal atom sites and octahedral interstitial sites. Two important elemental defects that position themselves at octahedral interstitial sites are carbon and helium impurity a t o m s . In the Box M e t h o d , the crystal volume is partitioned into a set of equal, space-filling subvolumes called boxes. This type of partitioning is illustrated in fig. 11.10 in which a cubic volume is partitioned into 125 cubic

662

Binary Collision Approximation

[§11.3

6

1 1

16

21

2

7

12

17

22

3

8

13

18

23

4

9

14

19

24

10

15

20

25

5

Fig. 11.10. (a) A cubic volume partitioned into 125 cubic boxes of equal volume, (b) B o x numbering scheme for the first 25 boxes in the bottom layer KB = 1. The box indices ( I B , J B , K B ) are defined by eq. (11.5).

§11.3]

11.3 Box method for locating defects and

atoms

663

boxes of equal v o l u m e . Partitioning into cubic boxes is a natural p r o c e d u r e for either bcc or fee crystals. Partitioning into rectangular parallepiped boxes is a natural p r o c e d u r e for hexagonal crystals. E a c h box contains a small n u m b e r of mesh points. In the case of bcc or fee crystals, the distance between successive mesh points along either the x-, yor z-axis is a half-lattice constant (hlc). One-fourth of the mesh points represent n o r m a l a t o m sites in a bcc crystal. In a fee crystal, one-half of the mesh points represent normal a t o m sites, In either instance, the remaining m e s h points represent octahedral interstitial sites. A cubic box containing 64 mesh points is described in fig. 11.11. E a c h mesh point is assigned the indices ( I , J , K ) which are the absolute coordinates of the mesh point in hlc units.

11.3.2

Indices of the nearest mesh

point

Given the absolute coordinates of a space point ( Χ Ι , Υ Ι , Ζ Ι ) in Ä , the indices ( I , J , K ) of the nearest mesh point are given by I = I N T ( X I / A N G + 0.5) J = I N T ( Y I / A N G + 0.5) K = I N T ( Z I / A N G + 0.5),

(11.1)

w h e r e A N G is the n u m b e r of A per hlc. T h e basis for this algorithm is illustrated in fig. 11.12 for two dimensions. T h e extension to three dimensions is straightforward. In t h r e e dimensions, the algorithm defines a cubic region with side A N G about each mesh point.

11.3.3

Box indices from mesh point

indices

T h e particular box which contains a given mesh point (I,J,K) is identified by the Box Indices ( I B , J B , K B ) which are given by IB = (I 4- N B D - 1)/NBD J B = (J + N B D - 1)/NBD K B = (K + N B D - 1 ) / N B D ,

(11.2)

w h e r e N B D is the box edge length in hlc. In the case of the box illustrated in fig. 11.11, N B D = 4 and we h a v e , for this particular case, IB = (I + 3)/4;

J B = (J + 3)/4;

K B = (K + 3)/4.

(11.3)

664

[§11.3

Binary Collision Approximation

JR=1

JR=2

IR=1

Ο 1

Ο 5

Ο 9

Ο 13

IR=2

Ο 2

Ο 6

Ο 10

Ο 14

IR=3

Ο 3

Ο 7

Ο η

Ο 15

IR=4

Ο 4

Ο 8

ο 12

Ο 16

JR=3

JR=4

KR=1

JR=1

JR=2

JR=3

JR=4

IR=1

Ο 17

Ο 21

Ο 25

Ο 29

IR=2

Ο 18

Ο 22

Ο 26

Ο 30

IR=3

Ο 19

Ο 23

Ο 27

Ο 31

IR=4

Ο 20

Ο 24

Ο 28

Ο 32

KR=2

§11.3]

Ii j Box method for locating defects and atoms

665

JR=1

JR=2

IR=1

Ο 33

Ο 37

Ο 41

Ο 45

IR=2

Ο 34

Ο 38

Ο 42

Ο 46

IR=3

Ο 35

Ο 39

Ο 43

Ο 47

IR=4

Ο 36

Ο 40

Ο 44

Ο 48

JR=1

JR=2

JR=3

JR=4

IR=1

Ο 49

Ο 53

Ο 57

Ο 61

IR=2

Ο 50

Ο 54

Ο 58

Ο 62

IR=3

Ο 51

Ο 55

Ο 59

Ο 63

IR=4

Ο 52

Ο 56

Ο 60

Ο 64

JR=3

JR=4

KR=3

KR=4

Fig. 11.11. Relative mesh point indices ( I R , J R , K R ) and mesh point position number ( N P O S ) for each mesh point in a cubic box containing 64 mesh points. ( I R , J R , K R ) and N P O S are defined by eqs. (11.6) and (11.7), respectively.

Binary Collision

666

Approximation

J-l

J

0+1

1-1

Ο

Ο

Ο

1+1

ο

Ο

[§11.3 γ

>

ο

Fig. 11.12. The region associated with a mesh point in an jcy-plane.

11.3.4

Box number

(NBOX)

T h e box n u m b e r N B O X is defined in terms of the n u m b e r of boxes N B E along each edge of t h e crystal (see fig. 11.10) and t h e box indices (IB,JB,KB). N B O X = IB + (JB - 1)*NBE + (KB - 1)*(NBE**2).

(11.4)

Given N B O X , we can obtain the box indices using the following recipe: KB NT JB IB

= ( N B O X - 1)/(NBE**2) + 1 = N B O X - ( K B - 1)*(NBE**2) =(NT-1)/NBE + 1 = N T - (JB - 1 ) * N B E .

(11.5)

Five boxes p e r edge are shown in fig. 11.10. This crystal size and the N B D = box size shown in fig. 11.11 will be used in the examples given h e r e . In practice, a larger crystal size usually is needed. For example, a crystal having a b o u t 100 boxes p e r edge with a box size N B D = 4 is n e e d e d to contain a 5 0 k e V cascade in iron.

11.3.5

Relative mesh point indices (IR,JR,

KR)

T h e indices ( I R , J R , K R ) of a mesh point, relative to the box in which it resides, are defined as

§11.3]

11.3 Box method for locating defects and

atoms

IR I - N B D * ( I B - 1) J-NBD*(JB-1) JR K R = Κ - N B D * ( K B - 1).

667

(11.6)

T h e s e indices are used to c o m p u t e the position n u m b e r N P O S for a given mesh point in its box.

11.3.6

Mesh point position number

(NPOS)

T h e position n u m b e r N P O S for a mesh point in its box is defined as N P O S = I R + ( J R - 1 ) * N B D + ( K R - 1)*(NBD**2).

(11.7)

Given N P O S , we can obtain the relative mesh point indices using the following recipe: KR NT JR IR

11.3.7

= ( N P O S - 1)/(NBD**2) + 1 = N P O S - ( K R - 1)*(NBD**2) =(NT-1)/NBD + 1 = N T - ( J R - 1)*NBD.

Utility ofNBOXand

(11.8)

NPOS

N B O X and N P O S are used to search for defects in a given local region during the selection of target atoms and to record the locations of newly p r o d u c e d defects. T h e periodicity of the crystal m a k e s the location of n o r m a l a t o m sites completely predictable; hence it is necessary to store only t h e locations of defects. Defect positions are stored and retrieved using the arrays K B O X ( N ) and L I S T ( N P O S , N ) . K B O X ( N ) is a table of box n u m b e r s for boxes which contain o n e or m o r e defects. L I S T ( N P O S , N ) tells what type of defect is located at position n u m b e r N P O S in the box with box n u m b e r K B O X ( N ) . For e x a m p l e , L I S T = 0 could be used to d e n o t e a vacancy and L I S T = 4 could be used to d e n o t e an interstitial. In o r d e r to illustrate the evaluation and use of the indices and location n u m b e r s used in the Box M e t h o d , consider the bcc crystal with N B E = 5 (fig. 11.10) and N B D = 4 (fig. 11.11). Suppose it contains a vacancy at mesh point (4,6,18) and an interstitial at mesh point (9,11,19). It will be assumed that the vacancy was p r o d u c e d first and that the box n u m b e r associated with the vacancy will be the first entry in the K B O X table. F o r the case of the vacancy we have the following box indices, relative

Binary Collision

668

Approximation

[§11.3

indices, box n u m b e r and position n u m b e r : IB = ( 4 + 4 - l ) / 4 = l JB = ( 6 + 4 - l ) / 4 = 2 K B = (18 + 4 - l ) / 4 = 5

I R = 4 - 4*(1 - 1 ) = 4 JR= 6-4*(2 - 1 ) = 2 K R = 18 - 4*(5 - 1 ) = 2.

T h e box n u m b e r is N B O X = 1 + (2 - 1)*5 + (5 - 1)*25 = 106 and t h e position n u m b e r in the box is N P O S = 4 + (2 - 1)*4 + (2 - 1)*16 = 24. F o r the case of the interstitial, the indices, box n u m b e r and position n u m b e r are IB = ( 9 + 4 - l ) / 4 = 3 J B = (11 + 4 - l ) / 4 = 3 K B = (19 + 4 - l ) / 4 = 5

I R = 9 - 4*(3 - 1) = 1 J R = 11 - 4*(3 - 1) = 3 K R = 19 - 4*(5 - 1) = 3.

T h e box n u m b e r is N B O X = 3 + (3 - 1)*5 + (5 - 1)*25 = 113 and the position n u m b e r is N P O S = 1 + (3 - 1)*4 + (3 - 1)*16 = 4 1 . O n the basis of these results, the arrays K B O X and L I S T would a p p e a r as follows: Ν KBOX(N) 1 106 (vacancy) 2 113 (interstitial) LIST(24,1) = 0 (vacancy) LIST(41,2) = 4 (interstitial). L I S T ( N P O S , N ) is used routinely during cascade simulation in the T A R S E L routine which selects target atom sites. After a potential target a t o m site has b e e n identified, the first test m a d e is the defect occupancy test. If this test shows that the site either is vacant or is occupied by a single impurity a t o m , t h e subsequent description of the collision process usually is simple. If the site either is occupied by a split interstitial or is part of a multiple target complex, the subsequent description of t h e collision process is not simple a n d , usually, must resort to dynamical m e t h o d result catalogs for a physically reasonable solution. A split interstitial target configuration seldom occurs in t h e simulation of an isolated cascade in a virgin crystal pertinent to P K A energies associated with n e u t r o n irradiation in fission

§11.3]

11.3 Box method for locating defects and

atoms

669

reactors. A split interstitial target configuration tends to be c o m m o n p l a c e , however, in the simulation of either high-fluence fission reactor n e u t r o n irradiation or self-ion irradiation by ions with energies above 2 M e V in iron and nickel. W h e n the target is a split interstitial, its axial orientation, relative to the projectile direction, determines the o u t c o m e of the collision. In this regard, L I S T ( N P O S , N ) can b e set u p to identify interstitial orientation as well as interstitial position. F o r e x a m p l e , in a bcc crystal t h e r e are six distinct axial orientations for a (110) split interstitial. H e n c e , the assignments L I S T ( N P O S , N ) = 4,5,6,7,8,9, for example, could be used to identify each of these distinct axial orientations. Additional discussion of split interstitial target configurations is given in § 11.9.

11.3.8

Normal atom site mesh

points

F o r a box whose mesh points are indexed in the way shown in fig. 11.11, t h e indices ( I , J , K ) for a normal atom site mesh point in a bcc crystal are either all even integers or all odd integers. All other index combinations represent octahedral interstitial sites. In the case of a fee crystal, the indices for a normal a t o m site mesh point are either all odd or o n e is odd and the o t h e r two are even. A g a i n , any o t h e r combintion of indices represents an octahedral interstitial site. A simple numerical test can be m a d e to identify a n o r m a l a t o m site mesh point. Consider the sum MSUM = MOD(I,2) + MOD(J,2) + MOD(K,2).

(11.9)

Since M O D ( N , 2 ) is zero w h e n Ν is an even integer and is o n e when Ν is o d d , M S U M is equal to the n u m b e r of mesh point indices which are o d d integers. H e n c e , for a normal a t o m site mesh point in a bcc crystal M S U M = 0 or 3

(bcc crystal)

(11.10)

and for a normal site mesh point in a fee crystal M S U M = 1 or 3

(fee crystal).

(11.11)

In the case of a fee crystal, we also can test N S U M = (I + J + K ) . If M O D ( N S U M , 2 ) = 1, the mesh point is a normal a t o m site; otherwise it is an octahedral interstitial site.

11.3.9

Rectangular

boxes

In general, the crystal studied will be a rectangular parellelepiped as will b e the boxes contained within it. Let N D X , N D Y and N D Z be the box

Binary Collision

670

Approximation

[§11.4

dimensions along the x-, y- and z-directions, respectively. Given the mesh point index ( I , J , K ) , the associated box indices are given by IB = (I + N D X - 1)/NDX J B = (J + N D Y - 1)/NDY K B = (K + N D Z - 1 ) / N D Z .

(11.12)

T h e box n u m b e r is N B O X = IB + (JB - 1)*NXE + ( K B - 1 ) * N E X Y ,

(11.13)

w h e r e N X E and N Y E are the n u m b e r of boxes along the x- and y-edges of t h e crystal, respectively, and N E X Y = N X E * N Y E . Given N B O X , the associated box indices are given by KB NT JB IB

= ( ( N B O X - 1)/NXY) + 1 = N B O X - (KB - 1)*NEXY = ((NT-1)/NXE) + 1 =NT-(JB-1)*NXE.

(11-14)

T h e mesh point relative indices are I R = I - N D X * ( I B - 1) JR = J - N D Y * ( J B - 1 ) K R = Κ - N D Z * ( K B - 1)

(11.15)

and the mesh point position n u m b e r in the box is NPOS = IR + (JR - 1)*NDX + (KR - 1)*NDXY,

(11.16)

where N D X Y = N D X * N D Y . Given N P O S , the associated relative mesh point indices are KR NT JR IR

11.4

11.4.1

= = = =

( ( N P O S - 1)/NDXY + 1 NPOS - (KR - 1)*NDXY ( ( N T - 1)/NDX) + 1 N T - (JR - 1)*NDX.

(11-17)

Target atom selection

Cell

method

For simplicity, we will consider target a t o m selection in a p u r e metal in which the only elemental defects are vacancies and split interstitials. E a c h of these defect types is located at a normal atom site. T h e projectile atom is started from the point ( Χ Ι , Υ Ι , Ζ Ι ) along the direction D whose direction

§11.4]

11.4 Target atom

selection

671

> PX

PY

/

(XI, ΥΙ,ΖΙ)

X *K

Xf

Fig. 11.13. Cell method target selection geometry in two dimensions. Notation is explained in the text.

XFR

^YBK

Fig. 11.14. Cell method target selection geometry in three dimensions. Notation is explained in the text.

cosines are C A P , C B P , C G P ) . T h e geometry is shown in fig. 11.13 for two dimensions and in fig. 11.14 for t h r e e dimensions. T h e p a r a m e t e r s ( X B K , X F R ) , ( Y B K , Y F R ) and ( Z B K , Z F R ) define three pairs of planes

672

Binary Collision

Approximation

[§11.4

which b o u n d the cubic cell containing ( Χ Ι , Υ Ι , Ζ Ι ) . W e are interested in the points P X , P Y and P Z w h e r e the ray D intersects the x-, y- and z-plane faces of the cell. R X , R Y and R Z represent the distances from the point ( Χ Ι , Υ Ι , Ζ Ι ) to the points P X , PY and P Z , respectively. Assuming that all distances and coordinates are expressed in hlc, we have X B K = A I N T ( X l ) and X F R = X B K + 1.0; Y B K = A I N T ( Y l ) and Y F R = Y B K + 1.0; and Z B K = A I N T ( Z l ) and Z F R = Z B K + 1.0,

(11.18)

where the function A I N T ( B ) gives the integer part of the real n u m b e r Β in decimal form. T h e plane intersected depends upon the direction cosines for the projectile a t o m direction. Let X N E X T , Y N E X T and Z N E X T be the coordinates of the first x-, y- and z-planes intersected. In particular, X N E X T = X F R if C A P . G T . O X N E X T = X B K if CAP.LT.O YNEXT = YFRifCBP.GT.O Y N E X T = Y B K if CBP.LT.O

(11.19)

Z N E X T = Z F R if C G P . G T . O Z N E X T = Z B K if CGP.LT.O. T h e distances from the point ( X I , Υ Ι , Ζ Ι ) to the planes designated above are R X = ABS((XNEXT - X1)/CAP) R Y = A B S ( ( Y N E X T - Y1)/CBP) R Z = ABS((ZNEXT - Z1)/CGP).

(11.20)

T h e plane first intersected is that with the minimum distance of the t h r e e stated above. Let that minimum distance be R M . T h r e e cases arise: (1) If R M = R X , then set N V = 1 and

(2)

XI = X N E X T YI = Y l + R X * C A P ZI = Z l + RX*CGP If R M = R Y , then set N V = 2 and

(3)

XI = X I + R Y * C A P YI = Y N E X T Z l = ZI + RY*CGP If R M = R Z , then set N V = 3 and XI = XI + R Z * C A P YI = Y l + R X * C B P ZI = Z N E X T .

(11-21)

(11.22)

(11.23)

§11.4]

11.4

Target atom

selection

673

Having the intersection point ( Χ Ι , Υ Ι , Ζ Ι ) in Ä we now find the mesh point ( I , J , K ) nearest to ( Χ Ι , Υ Ι , Ζ Ι ) and construct the indices of the nine possible target a t o m mesh points that this index specifies. These mesh points are as follows: IfNV = l(X-plane) I I I

J-l

K-l

J

K-l

J+l

K-l

I I I

J-l Κ J Κ J +1 Κ

I I I

J-l

K+l

J

K+l

J+l

K+l.

(11.24)

IfNV = 2(Y-plane)

1-1 I

1+1

J J J

K-l

1-1

K-l

I

K-l

1+1

1-1

Κ Κ Κ

I

J J J

Κ Κ Κ

1-1

J+l

I

J+l

1+1

J+l

1+1

J J J

K+l

J J J

K+l

(11.25)

K+l.

I f N V = 3 (Z-plane)

1-1

J-l

I

J-l

1+1

J-l

Κ Κ Κ

1-1 I

1+1

Κ Κ K.

(11.26)

T h e appropriate set of nine mesh point indices defines mesh point coordinates and occupancy. Specifically, (I,J,K) implies N B O X and N P O S . Search of K B O X reveals if any defect exists in the box concerned. If a defect does exist, L I S T ( N P O S , N ) tells what type of defect resides at the site. If no defect exists, the o d d - e v e n rules for normal atom sites reveal if (I,J,K) is a n o r m a l a t o m site. E a c h mesh point which contains an a t o m is tested to see if that a t o m will be the target a t o m . This is d o n e by computing the a t o m position vector R P relative to ( Χ Ι , Υ Ι , Ζ Ι ) . T h e c o m p o n e n t s of R P are (XP2=FLOAT(I)-Xl,

YP2=FLOAT(J)-Yl,

ZP2=FLOAT(K)-Zl).

(11.27)

If R P exceeds the potential function cutoff R C A , the atom u n d e r consideration cannot be a target a t o m and we proceed to the next target a t o m possibility. If R P is less than R C A , we then examine the scalar product PROJ = XP2*CAP + YP2*CBP + ZP2*CGP.

(11.28)

P R O J is the projection of R P on the direction in which the projectile a t o m is moving. H e n c e , if P R O J is negative, the target candidate lies behind the projectile a t o m relative to the projectile arm direction of motion, and cannot be a target a t o m . Only a target candidate for which P R O J is positive can be the target a t o m . T h e target candidate selected is that o n e with the smallest positive P R O J of the lot considered. If two or m o r e candidates are

674

Binary Collision

Approximation

[§11.4

nearly equidistant from the projectile and also have nearly equal P R O J values at the minimum level, a multiple target collision event has b e e n e n c o u n t e r e d . Multiple target collisions are treated in §11.4.2 and §11.9. If no target candidate satisfies the conditions for a valid target a t o m , the p r o g r a m must then trace the projectile along the direction ( C A P , C B P , C G P ) from the intersection point ( Χ Ι , Υ Ι , Ζ Ι ) . This brings the trace into a new mesh point cell and sets u p a new intersection point for the trace ray. It may be necessary to trace through several cells before finding a target if the local region concerned has a high vacancy content. T h e target selection p r o c e d u r e described in this section is awkward and tedious. Unfortunately, all known target selection procedures are awkward and tedious, and also are expensive because of the necessity to test each possible target site for vacancy or interstitial occupancy. E a c h defect occupancy test requires a new search of the K B O X table for a specific box n u m b e r . T h e time n e e d e d to m a k e this search can be minimized by using a hashing technique. If a straightforward, linear search technique is used, the cost of simulation can become prohibitive. Hashing techniques are described in A p p e n d i x A 11.1.

11.4.2

Target selection: neighbor

method

If we need not be concerned with defects that can be located at octahedral interstitial sites, the neighbor search m e t h o d can be used to select a target a t o m . T w o important defects which can be located at an octahedral interstitial site are carbon and helium impurity atoms, for example. T h e neighbor search m e t h o d consists of determining the normal a t o m site (IS,JS,KS) nearest to the projectile atom and then regarding each of the first and second neighbors of (IS,JS,KS) as being a possible target atom site. Let N A Y B O R be the n a m e of the subroutine which m a k e s u p a list of coordinates for possible target a t o m sites. A straightforward example of N A Y B O R for a bcc crystal is as follows: C C

NEIGHBOR TABLE SUBROUTINE SUBROUTINE NAYBOR NC = 0 I F ( N N A B . G T . 14) G O T O 20 D O 10 Ν = 1,14 I T = IS + J X ( N ) J T = JS + J Y ( N ) K T = KS + J Z ( N ) CALL VACTEST I F ( V A C ) G O T O 10

§11.4]

11.4 Target atom

selection

675

NC = NC + 1 IN(NC) = IT JN(NC) = JT KN(NC) = KT 10 C O N T I N U E MTX = NC RETURN C 20 IN(1) = IS J N ( 1 ) = JS KN(1) = KS D O 30N=2,15 I T = IS + J X ( N - l ) J T = JS + J Y ( N - l ) K T = KS + J Z ( N - l ) CALL VACTEST I F ( V A C ) G O T O 30 NC = NC + 1 IN(NC) = IT JN(NC) = JT KN(NC) = KT 30 C O N T I N U E MTX = NC RETURN END T h e arrays Ν 1 2 3 4 5 6 7

JX, JY JX 1 -1 -1 1 1 -1 -1

and J Z a r e : JY JZ 1 1 1 1 -1 1 -1 1 1 -1 1 -1 -1 -1

Ν 8 9 10 11 12 13 14

JX 1 2 0 -2 0 0 0

JY -1 0 2 0 -2 0 0

JZ -1 0 0 0 0 2 -2

with ( J X , J Y , J Z ) being the relative coordinates of the eight first-neighbor sites in a bcc crystal for Ν = 1 , . . . , 8, and the relative coordinates for the six second-neighbor sites in a bcc crystal for Ν = 9, ..., 14. T h e index N N A B determines how many target sites need to be listed. If the reference site (IS,JS,KS) should be listed in the target table, then N N A B is set equal to 15; otherwise it is set equal to 14. In this regard, it is convenient to include the coordinates of the reference site that is associated with each

676

Binary Collision

Approximation

[§11.4

projectile a t o m candidate in the Q U E U E table along with the candidate's kinetic energy, position and direction. If the reference site was the target a t o m site for the projectile's last collision, the projectile energy can then be e n t e r e d in the Q U E U E table as a negative n u m b e r . This negative sign then serves as a signal that the reference site should not b e included in the table of possible targets for the projectile concerned. H e n c e , when a projectile is selected from the Q U E U E table and its energy Q E is negative, the p r o g r a m responds by first setting E P = - Q E and N N A B = 14 before calling N A Y B O R for a list of possible target atoms. T h e neighbor table must be edited to exclude each neighbor site for which either a vacancy is present, or the vector from the projectile a t o m position to the neighbor site has a negative projection on the projectile direction ( C A P , C B P , C G P ) , or the impact p a r a m e t e r exceeds the impact p a r a m e t e r cutoff SM A X , or the distance from the projectile position to the neighbor site exceeds the potential function cutoff R C A . SM A X normally is t a k e n to b e one-half of the interatomic distance. T h e result of this editing is the set of arrays, XTP(M), YTP(M), ZTP(M), P R O ( M ) , STP2(M), w h e r e ( Χ Τ Ρ , Υ Τ Ρ , Ζ Τ Ρ ) are neighbor site coordinates, P R O is the projec­ tion of the distance between the projectile position and the site position on the projectile direction, and STP2 is the square of the impact p a r a m e t e r . T h e s e arrays either represent a single target atom site or a collection of M T X target atom sites with P R O values within an interval D E L P of the minimum projection P R O M . A straightforward routine for building these arrays is as follows: Subroutine V A C T E S T , called in the D O 10 and D O 30 loops, is a routine which identifies vacancy atom sites using the K B O X scheme described in §11.3. T h e mesh point indices used are ( I T , J T , K T ) . If a vacancy is present, the routine returns the signal V A C = . T R U E . . If there is no vacancy at the site, it returns the signal V A C = . F A L S E . . A test for split interstitial occupancy is automatically included in the K B O X scheme but we have omitted it h e r e for sake of simplicity. A split interstitial target configuration is a special variation of the multiple target configuration collision p r o c e d u r e and is discussed in §11.9. If M T X = 1, the routine has selected a single target a t o m site. This is the so-called normal result for a binary collision simulation. If M T X > 1, the routine has selected a set of multiple targets. T h e multiple target case is discussed further in §11.9. If M T X = 0, the routine has been unable to find a target atom site. This occurs only w h e n the projectile is passing through a region which contains a multitude of vacancies. In such a case, the projectile is stepped forward

§11.5]

11.5 Collision

geometry

671

I YL

Fig. 11.15. Binary collision geometry in the local laboratory system.

along the direction ( C A P , C B P , C G P ) until a new reference site is found which is a first neighbor of the previous reference site. 11.5

Collision geometry

This section deals with the collision geometry for an elastic collision of two equal masses. This special case illustrates a procedure for determining the coordinates of the takeoff point for each collision p a r t n e r , after collision, in t h e fixed laboratory coordinate system. A procedure for determining the takeoff points w h e n the masses are not equal and when inelastic loss occurs is described in §11.8. Let the coordinate axes for the local laboratory coordinate system be designated as X L , Y L and Z L , w h e r e the direction of Z L is that of the projectile a t o m velocity as it enters the collision. Fig. 11.15 describes the collision geometry in the local laboratory system. A s the projectile a t o m

678

Binary Collision

Approximation

[§11.5

approaches the stationary target atom from a distance beyond the potential function cutoff radius RC, the two atoms begin their collision interaction w h e n their separation distance r first falls below RC. In fig. 11.15, t h e beginning of the collision corresponds to the projectile atom being at point A and the stationary target atom being at point C. T h e coordinates of point A are (0,0,—D) and those of point C are (0,5,0), w h e r e 5 is the impact p a r a m e t e r . In the fixed laboratory coordinate system, the coordinates of A and C are ( Χ Ι , Υ Ι , Ζ Ι ) and ( X 2 , Y 2 , Z 2 ) , respectively. T h e center of mass is located at point Β at the start of the collision. T h e coordinates of point Β are (0,£>/2,5/2). During the first half of the collision, the a t o m pair separation distance decreases monotonically from RC to RO, the distance of closest approach. A t this juncture, the center of mass is located at point Μ whose coordinates are [ 0 , ( - D / 2 + v c i 0 ) , 5 / 2 ] , w h e r e v c is the center of mass velocity and t0 is the time at mid-collision. Thereafter, the a t o m pair separation distance increases monotonically from RO to RC. T h e collision ends when r = RC because the potential function is zero for r ^ RC. Calculation of the center of mass scattering angle 0 c m, the projectile scattering angle 0, the target scattering angle φ and the collision time 2t0 is discussed in §11.7. These four quantities are assumed to be known in the following discussion. A t the end of the collision, the center of mass is at point I and the projectile and target atoms are at points G and K, respectively. T h e coordinates of point I are [0,(—D/2 4- 2 v c i 0 ) , 5 / 2 ] . T h e coordinates of point G are 0 XG YG 5/2 - (flC/2)sin(0 c 'cm Z G = -D/2 4- 2vct0 + ( f l C / 2 ) c o s ( 0 c m + y),

(11.29)

where γ = arc sin(5/i?C). Point F is used as the takeoff point for the projectile atom when its next collision is calculated. Similarly, point J is used as the takeoff point for the target a t o m w h e n , as a projectile a t o m , its next collision is calculated. Let the coordinates of point F be (Ο,Ο,ΖΡΙ). T h e n , from the figure, it follows that tan0= YG/(ZG-ZP1) Z P 1 = Z G + YG/tanO.

(11.30)

For the special case of equal masses, Z P 1 turns out to be Z P 1 = (2vct0 - D). Let the coordinates of point J be ( 0 , 5 , Z P 2 ) . T h e n from the figure we have Z P 2 = Z P 1 4- 5/tan<£.

(11.31)

§11.5]

11.5 Collision

geometry

679

Let ux, uy and uz be unit vectors along the X L , Y L and Z L axes, respectively, and let the c o m p o n e n t s of these unit vectors in the fixed laboratory system be designated as follows: ux = ( C X X , C X Y , C X Z ) , uz = ( C Z X , C Z Y , C Z Z ) .

uy = ( C Y X , C Y Y , C Y Z ) ,

Sine uz is along Z L , and hence along the projectile velocity at the beginning of the collision, we have C Z X = C A P , C Z Y = C B P and C Z Z = C G P .

(11.32)

T h e local laboratory coordinate system Y L and Z L axes lie in the collision plane which is defined by the vectors uz and RC = [(X2 — X I ) , ( Y 2 — Y 1 ) , ( Z 2 — Z l ) ] . T h e X L axis is normal to the collision plane and is directed outward toward the reader in fig. 11.15. T h e direction n u m b e r s for the X L axis in the fixed coordinate system are the components of the vector product RC X M 2 , which are XN = RCY*CZZ - RCZ*CZY YN = RCZ*CZX - RCX*CZZ ZN = RCX*CZY - RCY*CZX. F r o m these direction n u m b e r s , the direction cosines for the X L axis are CXX = XN/BN, CXY = YN/BN, CXZ = ZN/BN.

(11.33)

where B N = S Q R T ( X N * * 2 + YN**2) is the magnitude of the vector product. Since uz is a unit vector, it is also true that B N = ÄCsiny = 5 ; hence we can use C X X = X N / S , C X Y = Y N / S , and C X Z = Z N / S . T h e direction cosines for the Y L axis are the components of the vector product uz x ux. These direction cosines are CYX = CXZ*CZY - CXY*CZZ CYY = c x x * c z z -

cxz*czx

(11.34)

CYZ = CXY*CZX - CXX*CZY. In the local laboratory coordinate system, the direction of the projectile a t o m exit asymptote is defined by the unit vector (0,—sin#,cosf9) and the direction of the target atom exit asymptote is defined by the unit vector (O,sin0,cos0). T h e direction cosines of the projectile atom exit asymptote in the fixed laboratory coordinate system are then CAP = -SINT*CYX + COST*CZX CBP = -SINT*CYY + COST*CZY CGP = -SINT*CYZ + COST*CZZ

(11.35)

680

Binary Collision Approximation

[§11.6

and those for the target a t o m exit asymptote are C A T = SINF*CYX + COSF*CZX CBT = SINF*CYY + COSF*CZY C G T = SINF*CYZ + COSF*CZZ.

(11.36)

In these expressions, SINT = sinö, C O S T = cosö, S I N F = sin0 and C O S F = cos0. T h e projectile a t o m takeoff point coordinates (Ο,Ο,ΖΡΙ) transform to XP = X I + (D + ZP1)*CZX YP = Y l + (D + ZP1)*CZY ZP = Z l + (D + ZP1)*CZZ,

(11.37)

in the fixed laboratory coordinate system. Having the takeoff point, we call Subroutine S I T E to find the nearest normal atom site to ( Χ Ρ , Υ Ρ , Ζ Ρ ) and store these reference site coordinates in the Q U E U E table along with projectile a t o m takeoff point, direction and energy E P . If the reference point is the a t o m site for the last target hit by the projectile, E P is multiplied by — 1 to signal this fact. In most instances, the reference site will be the a t o m site for the last target. T h e target a t o m takeoff point coordinates (0,S,ZP2) transform to X T = X 2 + S*CYX + Z P 2 * C Z X Y T = Y2 4- S*CYY + Z P 2 * C Z Y Z T = Z 2 + S * C Y Z + ZP2*C£>Z,

(11.38)

in the fixed laboratory coordinate system. T h e target a t o m reference site is the target a t o m itself.

11.6

Collision time and angle

T h e collision calculation will be described for a finite potential cutoff m o d e l . Gaussian integration will be used. T h e integrals for the collision angles and the collision time arise from the equation

§11.6]

11.6 Collision time and angle

2

2

681

112

dr/dt = [ v i ( l - s /? ) - (Alm^Vir)} ,

(11.39)

w h e r e vx is t h e projectile initial velocity in t h e laboratory coordinate system, πΐχ is t h e projectile mass, s is t h e impact p a r a m e t e r , r is t h e distance between projectile a n d target a n d V(r) is t h e central force potential energy for t h e projectile-target pair. T h e time integral from t h e time at which t h e atoms begin t o interact when r = Rc to t h e time of closest a p p r o a c h , r = R0, is

2

dt= ( - 1 / V l ) | * °

1/2

[(1 -s^r )

- V(r)/Er]- dr.

(11.40)

2

In e q . (11.40), Er = ra0v /2, w h e r e ra0 = m 2 / ( m 1 + m2) is t h e reduced 2 mass a n d m2 is t h e target a t o m . Substituting r = Ro/(l — u ) a n d d r = 2u2 2 du/(l — w ) , this becomes r ίο

L J

rο a t

—2RoH(u)udu

=( ^

L

y

ίΛ

ο

(1 —

l

2 w

(- ) n 41a

)

m

2

where 2

2

2

H(u) = [1 = 5 ( 1 - w ) / / ^ - y ( r ) / £ r ] ~

1 /2

(11.41b)

and m

UL=(l-Ro/Rc) .

(11.41c) 2

Bringing t h e quotient (u /Ro) 2

V(R0)/Er

through t h e radical and using t h e relation

2

= (i-s /R 0)

(11-42)

we obtain

'»- "'ί» ΤΓ^Ϋ

'

<1 ι)


where 2

2

2

F(u) = [s {s - u ) + (Rl/u Er){V(R0)

- V(Ro/(l

2

- " ))}Γ

1 /2

(11.43b)

Binary Collision

682

Approximation

[§11.7

tQ is the time for one-half the collision and dcm = t0Vi is the total distance m o v e d by the center of mass during the collision. T h e scattering angle in the center-of-mass coordinate system is given by

(11.44)

T h e integrals appearing in the expression for t0 and 0 c m can be integrated numerically at each collision using G a u s s - L e g e n d r e integration. Let the integrand be G ( u ) ; then according to the G a u s s - L e g e n d r e integration procedure

(11.45)

where ut = (uJ2)(l + *,·), the {A,-} are the G a u s s - L e g e n d r e coefficients and the are the G a u s s - L e g e n d r e arguments for nth order Gauss L e g e n d r e integration. G a u s s - L e g e n d r e coefficients and arguments are listed on p . 135 of Scheid [1968], for example, for η = 2 , 4 , 6 , . . . , 16.

11.7

Energy transfer

Following A p p e n d i x Β in Evans [1955], the scattering angle φ for the target always is φ = 0.5(180 - 0 c m ) degrees

(11.46)

independent of the projectile and target a t o m masses. F r o m this relation we obtain sin0 = c o s ( 0 c m/ 2 ) , cos0 = s i n ( 0 c m/ 2 ) and t a n 0 = s i n 0 c m/ ( l - c o s 0 c m) . (11.47) T h e velocity of the target a t o m after collision is 2

2

v 2 = 2 [ v / ( l + > l ) ] ( l - c o s Ö c m) ,

(11.48)

§11.7

11.7 Energy

transfer

683

(l+Acose

)/A

(a)

(1+fAcose

)/A

(b)

Fig. 11.16. Relation b e t w e e n the projectile atom scattering angle θ in the local laboratory coordinate system and the center of mass scattering angle Qcm: (a) elastic collision; (b) Leibfried's model for an inelastic collision (see §11.8).

w h e r e A = m2lmi. m1 is the projectile a t o m mass and m2 is the target a t o m mass, ν is the velocity magnitude of the projectile as it begins the collision. 2 2 W h e n A = l,v 2 = ( v / 2 ) ( l - c o s 0 c m) . Fig. 11.16 shows how the relation between the projectile a t o m scattering angle θ in the local laboratory system and the center of mass scattering angle 0 c m. F r o m this figure we obtain

Binary Collision

684

t a n 0 = Asinecml(l

Approximation

+ ^ c o s 0 c m) .

[§11.8

(11.49)

T h e energy of t h e target a t o m after collision is E T A = 0 . 5 m 2 v | = 2EQA(1

2

- cosf3 c m)/(l +A)

(11.50)

and that for the projectile a t o m is EPA = £ 0 - E T A

(11.51)

2

where E0 = O.Sm\V is the kinetic energy of the projectile at the beginning of t h e collision. T h e relations listed above pertain to a purely elastic collision. W h e n an inelastic loss occurs during a collision, the exit energies for the projectile and target are lessened by this dissipation. Let the inelastic loss be Ein. T h e elastic kinetic energy available is then £ a v l = E0 - Ein. O n e approxima­ tion is to c o m p u t e E T A as E T A = 2 £ a v l^ ( l - c o s 0 c m) / ( l + Λ ) by merely substituting £ a EPA = £ a

vl

2

(11.52)

for E0. T h e n E P A would be given by

vl

- ETA.

(11.53)

Robinson and T o r r e n s [1974] use Firsov's model to c o m p u t e an inelastic loss t e r m for each collision in a cascade simulation. Their recipe is 5

Exn = A/(l + B * R 0 ) e V ,

(11.54)

where 5/3

1/2

A = 0.05941(Zi + Z2) (E0/m1)

(11.55)

and £ = 0.3042(Ζ! + Ζ 2 )

1 / 3

.

(11.56)

Ζ ι and Z 2 are the atomic n u m b e r s for the projectile and target atoms, respectively, and R O is the distance of closest approach in A .

11.8

Collision asymmetry due to inelastic loss

T h e previous two sections emphasize the fact that an elastic collision is symmetric in space and time. Specifically, the center of mass travels the same distance in the local laboratory system during the first half of the collision, as the a t o m separation distance decreases from R C to RO, as it does during t h e second half of the collision, as the a t o m separation distance increases from RO to R C . In addition, the atom velocities are constant in the center of mass system and each a t o m trajectory is symmetric about ( 0 c m/ 2 ) .

11.8 Collision

§11.8]

asymmetry

due to inelastic

loss

685

lost

This symmetry is when inelastic loss occurs. This loss causes a m o n o t o n i c decrease in the relative energy Er during the collision. A s a consequence, the center of mass moves a different distance during the first half of the collision, in the local laboratory system, than it does during the second half of the collision. In addition, the atom velocities are not constant in the center of mass system and their trajectories are not symmetric about ( 0 c m/ 2 ) . N o general solution for the inelastic loss effect has been developed; however, Leibfried [1965] developed an approximation model for this process. In his m o d e l , the first half of the collision is considered to be purely elastic and the center of mass scattering angle is taken as being the center of mass scattering angle for a purely elastic collision. T h e relative energy during the second half of the collision is reduced relative to that during the first half such that we have

EP=fE?\ where 5

/ = (1 -

EJE^f .

O n this basis he obtained 2

2

E T A = [ / s i n ( 0 c m/ 2 ) + 0.25(1 -f) (4AE0/(l

2

+ A) ]

(11.57)

and E P A = E0-

E T A - Ein

(11.58)

for the target and projectile energies after collision. T h e effect that the inelastic loss has on the scattering angles is given by t a n 0 = ( £ 4 s i n 0 c m) / ( l + / , 4 c o s 0 c m)

(11.59)

tan0 = / s i n 0 c m/ ( l - / c o s 0 c m ) .

(11.60)

In this m o d e l , the distance traveled by the center of mass in the local laboratory system during the first half of the collision is x0 = vct0 and that traveled during the second half of the collision is x0/f. W e can incorporate Leibfried's model into the takeoff point expressions eqs. (11.29)-(11.31) by replacing the term 2vct0 in eq. (11.29) with the term v c i 0 ( l + 1//), by using t a n 0 i n eq. (11.30) as it is expressed in eq. (11.59) and by using eq. (11.31) as it is expressed in eq. (11.60).

tannin

11.9

Multiple targets

Knock-on a t o m trajectories routinely pass near the center of nearly symmetrical arrangements of two or m o r e target atoms in either a

686

Binary Collision

Approximation

[§11.9

replacement collision chain or during channeling. In these instances, a scheme for simulating nearly simultaneous collisions with several target atoms is required. T h e frequency with which these multiple target events occur in a cascade is sufficiently large that any cascade simulation p r o g r a m must contain a multiple target routine. Besco used the forbidden target scheme to carry a K A through a multiple collision event in the C A S C A D E Program. A s illustrated in §11.4.2, the target selection routine can detect the occurrence of a multiple target event by sensing the n u m b e r of targets whose projections on the ΚA path direction lie within an interval of width D E L P of the minimum positive projection P R O M . W h e n e v e r the presence of two or m o r e targets was detected, Besco's routine would flag this occurrence. This flag signaled the p r o g r a m to retain the projectile atom and start it from the same position through a sequence of binary collisions. T h e exit direction from the nth collision was used as the takeoff direction for the (n + l ) t h collision. T h e flag also signaled the p r o g r a m to reject as a possible target any target used after the flag was established. This forbidden target provision m a d e it impossible for the projectile a t o m to collide twice with the same target during the multiple target sequence. W h e n the set of targets with projections within the interval D E L P of the minimum projection P R O M had been exhausted, the multiple target flag was removed and the program proceeded in its normal binary collision m o d e . Besco's t r e a t m e n t of a multiple target event is an implicit t r e a t m e n t . In contrast, R o b i n s o n and Torrens use an explicit scheme in their M A R L O W E P r o g r a m . M A R L O W E selects target atom candidates using the neighbor m e t h o d . T h e specific tests used to identify nearly simultaneous targets are illustrated in fig. 11.17. These tests involve both the projection P R O and the impact p a r a m e t e r 5. A s before, P R O is the projection of the distance between the projectile atom and the target candidate o n t o the projectile direction. A n a t o m is considered to lie in the correct relative position if P R O is greater than P R O m i n > 0, where P R O m i n is the projection for the previous target. If, in addition, the impact p a r a m e t e r is less than 5 m a x, the a t o m is considered to be a viable target atom candidate. A multiple target event is identified in the following way: Let the target a t o m candidate with the smallest projection P R O * be designated as 1. For each other target atom candidate /, the quantities D E L i , = P R O , - PROx 2 Ol = S? + (DELU) 2 2 O\ = S + ( D E L ! , ) , are c o m p u t e d (see fig. 11.17b). Collision with target candidate / is considered to b e simultaneous with collision with candidate 1 if the three

§11.9]

11.9 Multiple

targets

687

PRO

Ρ

Φ-

4 4

(b)

PR0]

PRO« 2

I

Fig. 11.17. Geometrical quantities in the simultaneous collision scheme used by Robinson and Torrens [1974]. Ρ is the projectile atom and Τ is a target atom candidate, (a) P R O m ni is the projection associated with either the previous target or the reference site R. (b) T w o simultaneous target atoms.

inequalities DELU < D P m i n, Du < 5 m a x and Dn < 5 m a x are satisfied. D P m i n is an empirically determined p a r a m e t e r . A collision with each of the

Binary Collision

688

Approximation

[§11.10

"simultaneous" targets is treated as if it were occurring alone. T h e several deflections of the projectile atom then are added vectorially. This t r e a t m e n t of a multiple target collision event underestimates the projectile atom energy loss since the free recoil of the projectile, assumed in each of the individual collision calculations, does not, in fact, occur. This disadvantage is offset, to some extent, by the fact that spurious defocussing of a focussed collision sequence is not induced by this m e t h o d . T h e sequential approach can induce spurious defocussing. T h e results of the sequential (Besco) approach and the simultaneous collision (Robinson and Torrens) approach for a 3 2 8 e V projectile and two simultaneous targets are c o m p a r e d in fig. 11.18. T h e defocussing effect induced by the sequential m e t h o d is illustrated clearly in this example. N o satisfactory strategy for dealing with multiple collision events, within the framework of the B C A , has been devised. T h e root p r o b l e m , of course, is a many-body p r o b l e m . T h e difficulty is trying to do many-body work with a two-body tool. T h e interaction of a projectile atom with symmetrical rings of target atoms has been studied by Weijsenfeld [1967] and by A n d e r s o n and Sigmund [1965, 1966]. Unfortunately, their treatments are not applicable when the target a t o m arrangement is asymmetrical.

11.10

Q U E U E table

A s shown in fig. 11.7, a plot of the n u m b e r of projectile atoms in a cascade versus the n u m b e r of collisions in the cascade has a p e a k at about N C O L T / 2 , where N C O L T is the total n u m b e r of collisions in the cascade. Since only one collision can be computed at a time, a Q U E U E table is required which lists the energies, positions and directions of projectile atom candidates for subsequent collision calculations. Q U E U E table design can affect the quality of a cascade simulation and certainly affects the cost of a cascade simulation. In this context, quality of simulation m e a n s the degree of realism of the cascade structure developed by the simulation. Projectile atoms have been selected from the Q U E U E table according to t h r e e different prescriptions: (1) select the last entry in the Q U E U E table, (2) select the entry with the largest velocity, (3) select the entry most recently produced. Case 1: After each collision, the energies of the emerging projectile and target atoms are examined. If each energy exceeds a prescribed cutoff energy, t h e n that a t o m with the largest velocity is retained and used as the projectile for the next collision. T h e slower atom is stored in what becomes

§11.10]

11.10 QUEUE

table

689

34eV

S 1 = 1.25Ä 260eV

>

S2 =

(b)

1.258

Simultaneous collisions

Fig. 11.18. Sequential and simultaneous methods compared. The scatterings are based on a Moliere potential with al2 = 0 . 0 9 6 0 Ä . Ρ is the projectile atom; Tx and T 2 are target atoms.

the last position in the Q U E U E table. If the energy of only o n e a t o m exceeds the cutoff energy, that a t o m is retained and used as the projectile for the next collision. If the energy of neither atom exceeds the cutoff energy, the last entry in the Q U E U E table is t a k e n as the projectile for the next collision. This prescription provides a simulation in which the principal trajectory line of the cascade and its branches are developed first. T h e secondary trajectory line and its branches are developed second, and so on. It does not necessarily develop the cascade branches in chronological order.

690

Binary Collision

Approximation

[§11.11

Case 2: T h e Q U E U E table entry with the largest velocity is selected as the projectile atom for the next collision. This p r o c e d u r e traces out the trajectory of the fastest current K A . It develops the cascade branches in very nearly the true chronological order. Case 3: If the p r o g r a m keeps a running account of the elapsed time as the cascade develops, then each projectile can be associated with a definite collision time. In this case, a Q U E U E table entry with the earliest collision time is selected as the next projectile. Because the projectile and target atoms for a given collision have the same collision time, the fastest of these two entries is selected as the projectile for the next collision. In principle, this approach gives the best description of cascade chronology. A cascade simulation which uses the "last entry" Q U E U E table look-up p r o c e d u r e (Case 1) costs less than o n e which uses either of the other two procedures. A s things stand, it is not clear that the quality of cascade simulation is significantly improved by using either the "largest velocity" or "earliest t i m e " p r o c e d u r e . Either procedure requires a binary tree Q U E U E table structure; otherwise the simulation cost is prohibitive. T h e binary tree table structure is described in A p p e n d i x A l l . 2 .

11.11

Vacancy and interstitial production criteria

Dynamical, many-body interaction computer experiments indicate that a well-defined displacement energy threshold, Ed(D), exists for each P K A ejection direction, D , in a crystal when SCIC are a d o p t e d . In addition, a single displacement event is characterized by a separation distance vector S(E,D) and an interstitial axis orientation vector A(E,D). T h e dynamical c o m p u t e r experiments also indicate that only a single displacement event is produced whenever the P K A energy lies between Ed(D) and an u p p e r U U b o u n d energy E (D). In iron, E (D) is less than AEd(D). A s the P K A energy sweeps from the lower b o u n d , Ed(D), to the upper b o u n d of this energy range, the position and orientation of the interstitial may change, but no m o r e than o n e stable displacement event is produced. This behavior also is observed when T C I C are adopted. Taking into account the influence of thermal vibrations, no m o r e than o n e displacement event is p r o d u c e d in fee iron(m2) c o m p u t e r experiments, for example, when the P K A energy lies in the range 20 ^ Ε ^ 75 e V . T h e basic mechanism for single displacement production is a replacement collision chain. H e n c e , given the Κ A energy and direction in a B C A model cascade simulation, it is possible to stipulate whether or not a vacancy will be p r o d u c e d at the K A ejection site, provided Ed(D) is used as the production criterion. T h e site at which the associated interstitial is formed also can be

§11.11]

11.11 Vacancy and interstitial production

criteria

691

d e t e r m i n e d provided the B C A m o d e l cascade simulation p r o g r a m can accurately simulate replacement collision chains for K A energies below E"(D). Unfortunately, the B C A does not always properly describe replace­ m e n t collision chains in this low energy range. T h e two capabilities stated h e r e are the criteria for judging w h e t h e r or not a B C A cascade simulation p r o g r a m can properly describe vacancy and interstitial production. A n o t h e r way of providing an accurate description of vacancy and interstitial production in a B C A cascade simulation p r o g r a m is to utilize a stored catalog of S(E,D) and A(E,D) results from dynamical m e t h o d c o m p u t e r experiments in assigning the interstitial position. T h e vacancy production criterion is, of course, Ed(D). In this approach, the B C A cascade simulation p r o g r a m is used to develop ΚA trajectories for energies above U U E (D). W h e n e v e r the energy of a K A falls below E (D), the catalog is used to provide accurate, dynamical m e t h o d results for S(E,D) and A(£,£>). T h e p r o c e d u r e takes less time than would a continued B C A simulation of the K A trajectory. In addition, it significantly cuts down on the size of the Q U E U E table. T h e C O L L I D E P r o g r a m contains this type of vacancy and interstitial production criteria option. Cascade simulation p r o g r a m s which use a single-valued displacement energy criterion for establishing vacancy production and interstitial produc­ tion cannot provide a useful defect distribution. Their utility is limited to estimating the total n u m b e r of displacements. W h e n e v e r the defect distributions they provide are realistic, it is a fortuitous occurrence. T h e vacancy and interstitial production criteria used in C A S C A D E , SIC and in a recent revision of M A R L O W E allow a Κ A to be ejected when its energy is comparable with the cohesive energy per atom. W h e t h e r or not it travels afar or is pushed back o n t o its initial site then d e p e n d s u p o n the a r r a n g e m e n t of atoms in the local environment. This approach rejects the use of a single "displacement energy" and substitutes, as far as we can in a B C A , many-body influences in determining vacancy production. T h e ability of C A S C A D E and SIC to handle replacement chains is surprisingly apt, even though they use the sequential approach for multiple target collisions. Several three-dimensional models of C A S C A D E cascades w e r e constructed. In viewing these models it is immediately obvious that all of the C A S C A D E cascades were dominated by replacement chains. P h o t o g r a p h s of low energy cascade models appear in figs. 11.19 and 11.20. It is pertinent to r e m a r k here that the shapes of the C A S C A D E P r o g r a m displacement spikes (depleted zones) simulated during 1963-65 are striking­ ly similar to the displacement spike shapes observed several years later by Seidman's group in field ion misroscope (FIM) experiments. These F I M experiments were performed during 1969-70 and the results a p p e a r in B e a v a n et al. [1971], Scanlan et al. [1971] and Seidman [1973,1975]. A comparison of displacement spike shapes given by F I M observations and c o m p u t e r experiments is given in fig. 11.21.

Binary Collision

692

Approximation

[§A11.1

Fig. 11.19. Collided atom region for a 5 0 0 e V collision cascade in bcc iron(m2) as given by the C A S C A D E Program. T h e cascade consisted mainly of two [001] replacement chains and a [Ϊ01] replacement chain. Black elements represent atoms that were permanently moved from their initial sites. White elements represent struck atoms which returned to their initial sites. Gray elements represent interstitials.

Appendix A l l . l

Hashing technique

T h e defect occupancy test portion of the target selection routine requires a search of the table K B O X ( N ) for a specific box n u m b e r N B O X . If N B O X is listed in K B O X , its index is required; if N B O X is not listed in K B O X , a quick reply to this effect is required. T h e task at hand is to show how to enter a new box n u m b e r in K B O X and how to search K B O X for a specific box number. T h e straightforward, linear search procedure operates by making a new box n u m b e r the last entry in the K B O X table. T h e search for a particular box n u m b e r is performed as follows: D O 10 Κ = 1, K M X I F ( N B O X . E Q . K B O X ( K ) ) G O T O 20 10 C O N T I N U E D E F E C T = .FALSE. G O T O 30 20N = K DEFECT = .TRUE. 30 C O N T I N U E RETURN

§A11.1]

AI 1.1 Hashing

technique

693

Fig. 11.20. Collided atom region for a 1 k e V collision cascade in bcc iron(m2) as given by the C A S C A D E Program. The cascade was dominated by a [111] replacement collision chain. Black elements represent atoms that were moved permanently from their initial sites. White elements represent struck atoms which returned to their initial sites.

If K B O X contains m a n y entries, the linear search m e t h o d is prohibitively expensive. A hashing technique can be used to search for a specific entry in a table at a fraction of the cost of a linear search. A hashing technique uses some property of the box n u m b e r N B O X to define its probable index Ν in the K B O X table. T h e recipe used to specify the probable index is called a hashing function. T h e probable index also is called the " h o m e index". T w o simple hashing functions suitable to searching K B O X will be described. T h e y are the m o d u l a r hashing function and the r a n d o m hashing function. T h e a u t h o r is indebted to A l t o n Coulter for teaching him about hashing functions for collision cascade simulation programs. (a) M o d u l a r hashing function: in this case the hashing function is Ν = M O D (NB O X , N M A X ) + 1,

(All.l)

where N M A X is the m a x i m u m n u m b e r of entries for K B O X . This function is used to obtain the p r o b a b l e index N . If another box n u m b e r already is stored in K B O X ( N ) , a linear m e t h o d is then used to find an unused location in the table. A s an example of how a box n u m b e r is entered in K B O X , consider the

Fig. 11.21. Displacement spike (depleted zone) shapes in bcc metals given by computer experiments and by field ion microscopy, (a) A 2 0 k e V displacement spike observed by Beavan et al. [1971] in tungsten, (b) A 1 5 k e V displacement spike simulated in bcc iron(m2) by Beeler and B e s c o (Beeler [1966]). B o t h the FIM observations and the computer experiments indicate: (1) a spongy vacancy-rich region at the spike center and (2) elongation of the spike along (110) directions.

§A11.1]

AI 1.1 Hashing

I s NSUM . L T . 1 5

-®-Λ

technique

Table i s

695

Full

Stop

Calculation

Ν = M0D(NB0X,15) + 1

λ Is

KBOX(N) . L T . 0

[

Ν = Ν + 1

(γ)

»

Assiqn KBOX(N) = NBOX NSUM = NSUM + 1

- > RETURN

A R o u t i n e f o r S t o r i n g NBOX

IF(N.GT.15)N=1

Fig. A l l . l 125-box system of fig. 11.10 a n d suppose that box n u m b e r s 14, 97, 57, 4 3 , 106, 69 a n d 118 w e r e e n t e r e d in K B O X in t h e order of the stated sequence. A s s u m e that K B O X is initialized by setting K B O X ( N ) = - 1 for all Ν a n d that N M A X = 15. Using t h e hashing function ( A l l . l ) , as d i a g r a m m e d in fig. A l l . l , t h e following K B O X table results:

Ν 1 2 3 4 5

KBOX(N) 118 106 -1 -1 -1

Ν 6 7 8 9 10

KBOX(N) -1 -1 97 -1 69

Ν 11 12 13 14 15

KBOX(N) -1 -1 57 43 14

N o t e that the hashing function assigns the same p r o b a b l e index, Ν = 14, for 43 a n d 118. Because of this conflict a n d because N B O X = 14 previously h a d b e e n placed in K B O X ( 1 5 ) , N B O X = 118 was placed in K B O X ( l ) . T h e hashing function a p p r o a c h is efficient provided less t h a n 6 6 % of t h e table positions a r e occupied. T h e p r o c e d u r e for using e q . ( A l l . l ) t o search for a specific N B O X is d i a g r a m m e d in fig. A l l . 2 . This search p r o c e d u r e b e c o m e s a linear search w h e n K B O X ( N ) for t h e p r o b a b l e index Ν is already occupied.

696

Binary Collision

NT = 0

Μ = M0D(NB0X,15)

Approximation

A Routine

for

Searching

f o r NBOX

[§A.11.1

+ 1

NT = NT+1

DEFECT = . F A L S E . Is

NT.GT.15

NBOX n o t i n

RETURN

Table

τ IF

KBOX(M).LT.O

DEFECT = . T R U E . Is

KBOX(M) = NBOX

Μ

RETURN

Μ = M+l I F ( M . G T . 1 5 ) M=l

Fig.

A11.2

(b) R a n d o m hashing function: in this scheme, N B O X is used as t h e seed integer for a r a n d o m n u m b e r generator routine and the first r a n d o m n u m b e r given by the routine is used t o obtain the probable index Ν according t o t h e recipe N = RN*15.0 + 1.0,

(A11.2)

w h e r e R N is a r a n d o m n u m b e r uniformly distributed o n (0,1). If the location K B O X ( N ) is used, t h e second r a n d o m n u m b e r generated by t h e routine is used t o obtain a n o t h e r N . This process is continued until an unfilled location is found. A n example of placing box n u m b e r s in K B O X o n the basis of eq. (A11.2) was w o r k e d o u t for the same sequence of box n u m b e r s used in the m o d u l a r hashing function example. T h e r a n d o m n u m b e r generator routine a d o p t e d is

§A11.2]

A11.2

Heap table

technique

697

KX = NBOX 10 K B I G = 24298*KX + 99991 KY = MOD(KBIG,199017) R N = KY/199017 KX = KY N = RN*15.0 + 1.0 T h e key n u m b e r s which describe the progress of this example are:

NBOX

KY

RN

Ν

KBOX(N)

14 97 57 43 106

42129 68693 91858 149720 88358 27279 184417 180917 133761

0.21168 0.34516 0.46156 0.75230 0.44397 0.13707 0.92664 0.90905 0.67211

4 6 7 12 7 3 14 14 11

free free free free used free free used free

69 118

T w o sampling attempts w e r e required in each instance to find free locations for N B O X = 106 and N B O X = 1 1 8 . T h e K B O X table produced is

Ν

KBOX(N)

1 2 3 4 5

Appendix A I 1.2

-1 -1 106 14 -1

Ν 6 7 8 9 10

KBOX(N) 97 57 -1 -1 -1

Ν 11 12 13 14 15

KBOX(N) 118 43 -1 69 -1

H e a p table technique

T h e Q U E U E table is a group of one-dimensional arrays which define projectile a t o m candidates. Specifically, it consists of the projectile a t o m energy array Q E ( K ) , the projectile atom velocity array Q V ( K ) , the

698

Binary Collision Approximation

[§A11.2

Fig. A l l . 3

projectile a t o m takeoff point coordinates arrays ( Q X ( K ) , Q Y ( K ) , Q Z ( K ) ) , the projectile atom direction cosines arrays ( Q A ( K ) , Q B ( K ) , Q G ( K ) ) and the projectile a t o m type array Q T ( K ) . Most workers feel that it is best to use the atom with the largest velocity as the projectile a t o m in the next collision calculation. T h e most efficient scheme yet devised for this purpose requires the Q U E U E table to be structured as a binary tree array which also is called a h e a p table. In such an array, the projectile a t o m candidate with the largest velocity always is stored in the Q U E U E table arrays with the index Κ = 1. A binary tree array is illustrated in fig. A l l . 3 for the velocity table Q V ( K ) . A n y given entry Q V ( K ) is less than or equal to the entry Q V ( K / 2 ) . For example, each of the entries QV(10) and Q V ( l l ) is less than or equal to the entry Q V ( 5 ) since, in integer arithmetic, each of the quotients 10/2 and 11/2 is equal to 5. Similarly, each of the entries QV(4) and Q V ( 5 ) is less than orequaltoQV(2). T h e properties of the binary tree table defined above d e t e r m i n e how we add a new entry to a binary tree array. Suppose Q V contained 10 entries (N = 10) and a new entry V P was to be added to the array. T h e entry routine, in effect, sets Q V ( l l ) = V P and then tests to see if V P ^ Q V ( 5 ) . I f this relationship is satisfied, the entry task is finished. If V P exceeds Q V ( 5 ) , the routine m a k e s the switch Q V ( l l ) = QV(5) and then tests Q V ( 5 ) to see if Q V ( 5 ) ^ Q V ( 2 ) , and so forth. T h e n u m b e r of entries in a filled h e a p table with Μ rows (see fig. A l l . 3 ) is Ν = (2

M

— 1)

§A11.2]

A II. 2 Heap table

technique

699

If a h e a p contains Μ rows, t h e n the addition of a new entry requires n o more than (Μ - 1) relative magnitude tests and entry location switches. H e n c e , if a h e a p contained 1000 entries (10 rows), the maximum n u m b e r of m a g n i t u d e tests required in making a new entry would be 9. A m a x i m u m of 1000 magnitude tests would b e required if a linear structure, o r d e r e d m a g n i t u d e array system were used. T h e following routine can be used to m a k e a new entry in a binary tree array. Let Ν be the n u m b e r of entries in the Q V ( K ) table and let V P be the new entry to b e a d d e d to the array. C

A D D VP T O QV(K) W H I C H CONTAINS Ν ENTRIES I F ( N . L T . l ) G O T O 30 N=N+1 I=N 10 J = I / 2 I F ( ( I . L E . l ) . O R . ( Q V ( J ) . G T . V P ) ) G O T O 20 QV(I)=QV(J) I=J G O T O 10 20 Q V ( I ) = V P RETURN 30 Q V ( 1 ) = V P N=l RETURN

This routine constructs a binary tree array as illustrated in fig. A l l . 3 . T h e following routine removes Q V ( 1 ) from the array Q V and t h e n restructures the r e m a i n d e r of the table so that it retains the fundamental properties of the binary tree array illustrated in fig. A l l . 3 . C C C

R E M O V E QV(1) F R O M QV A N D R E S T R U C T U R E T H E A R R A Y . T H E R E A R E Ν E N T R I E S IN Q V AT THE BEGINNING. IF(N.LT.l)GOTO40 VP=QV(1) QLAST=QV(N) N=N-1 1=1 J=2 10 I F ( J . G T . N ) G O T O 30 J=2*I QL=QLAST IF(J.LE.N) Q L = Q V ( J )

Binary Collision

700

Approximation

[§A11.2

QR=QLAST IF((J+1).LE.N) Q R = Q V ( J + 1 ) I F ( ( Q L A S T . G E . Q L ) . A N D . ( Q L A S T . G E . Q R ) ) G O T O 30 I F ( Q L . L E . Q R ) G O T O 20 QV(I)=QL I=J G O T O 10 20 Q V ( I ) = Q R I=J+1 G O T O 10 30 Q V ( I ) = Q L A S T RETURN 40N=0 RETURN END Fig. A l l . 4 illustrates the development of a binary tree array Q V ( K ) as new entries V P = 14, 97, 57, 43, 106, 69 and 118 are m a d e in the order of the stated sequence. It is a simple matter to follow the steps defined in the routine for adding V P to Q V ( K ) and to understand how a binary tree array is produced in which the element in Q V ( 1 ) is always the largest element. If we take the final array in fig. A l l . 4 , for which Ν = 7, and remove Q V ( 1 ) = 118 according to the removal routine, the tactics of the removal routine b e c o m e apparent. In particular, the routine first finds Q L = 97 and Q R = 106 which are in Q V ( 2 ) and Q V ( 3 ) , respectively, which tie directly to Q V ( 1 ) . It then elevates the larger of these two elements (106) into Q V ( 1 ) which, by definition of the binary t r e e , must contain the largest element in the t r e e . T h e element elevated into Q V ( 1 ) initially was in Q V ( 3 ) . N o w that Q V ( 3 ) is "vacant", it must be filled by switching the contents of either Q V ( 6 ) or Q V ( 7 ) into Q V ( 3 ) . Again, is is the larger of t h e two elements which must be switched, and the routine proceeds on this basis. In the present example, this is the end of the restructuring part of the binary tree array which always follows removal of Q V ( 1 ) , the largest element. A s in this example, the restructuring process always continues down to the row containing Q L A S T .

References A n d e r s e n , Η. H. and P. Sigmund [1965] Phys. Lett. 15 237 A n d e r s e n , Η. H. and P. Sigmund [1966] Kgl. D a n . Vid. Sels. Mat.-Fys. Medd. 3 4 N o . 15 B e a v a n , L. Α . , R. M. Scanlan and D . N . Seidman [1971] Acta Met. 19 1339 B e e l e r , J. R., Jr. [1966] Phys. Rev. 1 5 0 4 7 0

§Ref]

References

4

5

6

7

701

4

5

6

Fig. A l l . 4

Beeler, J. R., Jr. [1967] Neutron irradiation effects saturation in alpha-iron, in: Hasiguti, R. R. e d . , Lattice Defects and their Interactions, Gordon and Breach, pp. 619-652 Beeler, J. R., Jr. [1971] The use of computer experiments to predict radiation effects in solid materials, in: Kriegel, W. W. and H. Palmour III, eds., Ceramics in Severe Environments, Plenum Press, pp. 553-574 Beeler, J. R., Jr. and M. F. Beeler [1973] A t o m bombardment simulation program ( C O L L I D E ) , in: U S A - A E C Report O R O - 3 9 1 2 - 2 1 , April, pp. 2-21 Beeler, J. R., Jr. and D . G. B e s c o [1962] Knock-on cascades and point defect configurations in binary materials, in: Radiation D a m a g e in Solids, Vol. 1, I A E A , Vienna, pp. 4 3 - 6 3

702

Binary Collision

Approximation

[§Ref

Beeler, J. R., Jr. and D . G. Besco [1963] J. Appl. Phys. 34 2873 B e s c o , D . G. and N . R. Baumgardt [1965] C A S C A D E and C L U S T E R , U S A - A E C Report G E M P - 3 5 6 , April B e s c o , D . G. and J. R. Beeler, Jr. [1963] Computer programs describing collision cascades in binary materials III. Body-centered cubic and face-centered cubic structures, U S A - A E C Report G E M P - 2 4 3 , August Evans, R. D . [1955] The A t o m i c Nucleus, McGraw-Hill Firsov, Ο. B. [1959] Sov. Phys.-JETP 36 1076 Leibfried, G. [1965] Bestrahlungseffekte in Festkörpern, Teubner, Stuttgart Lekkerker et al. [1971] The SIC and R E L A X programs, Delft Reactor Group, informal report Robinson, Μ. T. and I. M. Torrens [1974] Phys. Rev. B 9 5008 Scanlan, R. M., D . L. Styris and D . N . Seidman [1971] Phil. Mag. 23 1459 Scheid, F. F. [1968] Theory and Problems of Numerical Analysis, Schaum Outline Series, McGraw-Hill Seidman, D . N . [1973] J. Phys. F:3 393 Seidman, D . N . [1975] Field-ion microscope studies of the defect structure of the primary state of damage of irradiated metals, in: A S M Materials Science Seminar: Radiation D a m a g e in Metals, N o v e m b e r , Cincinnati, Ohio Weijsenfeld, Ο. H. [1967] Phillips R e s . R e p . Suppl. N o . 2 Yoshida, M. [1961] J. Phys. Soc. Japan 16 44