Surface green's function for a rumpled crystal surface

Computer Physics Communications 51(1988) 381—390 North-Holland, Amsterdam

381

SURFACE GREEN’S FUNCTION FOR A RUMPLED CRYSTAL SURFACE Franti~ekMACA

and Matthias SCHEFFLER

Fritz-Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 4—6, D-I000 Berlin 33, Germany Received 18 March 1988

The computer program allows to calculate the Green’s function of a three-dimensional system with two-dimensional translation symmetry. Examples for its application are the crystal surface and the interface between two crystals. The crystal can be composed from rumpled atomic layers, that is, the atoms in the layer unit cell may have small differences in heights. The layer-by-layer KKR scheme is used. An overlayer with a coincidence structure of a substrate layer may be involved; the unit cell of any layer may be composed of up to four different atoms. The Green’s function is evaluated in a spherical wave expansion, with basis functions centred at the atoms. It directly allows to calculate the surface electron charge density and the local and projected density of states.

PROGRAM SUMMARY Title of program: RUMPGF Catalogue number: ABFF Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland (see application form in this issue) Computer: CRAY X-MP/24; Installation: IPP der MaxPlanck-Gesellschaft, Garching Operating system: COS Programming language used: FORTRAN High speed storage required: 315 Kwords No. of bits in a word: 64

CPC: 38 (1985) 403; version 2: cat, no.: AAXZ, title: SURFACE GREEN’S Function vers. 2, ref in CPC: 47 (1987) 349 Both papers are referred to below as paper or computer code I. Nature of the physical problem The computer program (as that of paper I) allows to calculate, in any given range of energy E, the Green’s function for the one-electron Hamilton operator with a potential for a three-dimensional system with two-dimensional translation symmetry, i.e. a crystal surface or interface. The Green’s function satisfies the Bloch-function-type periodic boundary conditions parallel to the surface, for a given value of k~ 1and the outgoing-wave boundary conditions normal to the surface. In generalization to paper I, the new code allows to handle rumpled atomic layers, that is, the atoms composing a unit cell may have small differences in heights. Therefore, it is now possible to treat more complex systems, as for example, reconstructed clean surfaces or systems with small adsorbate—substrate distances.

No. of lines in combined program and test deck: 4787 Keywords: Green’s function, surface, density of states, charge density, layer KKR method, multiple scattering, solid state, condensed matter Reference to other published versions of this program: cat, no.: AADF, title: SURFACE GREEN’S FUNCTION, ref in Permanent address: Institute of Physics, Czechoslovak Academy of Sciences, Na Slovance 2, 180 40 Prague 8, Czechoslovakia.

Method of solution The layer-by-layer KKR scheme is used. The crystal (with or without an adsorbate layer) is divided into layers parallel to the surface. All the layers may be rumpled and their two-dimensional unit cell may contain up to four atoms. The potential of each layer is treated in the muffin-tin approximation. The intra-layer scattering is calculated with the method of Kambe [1], and the inter-layer scattering is calculated using the doubling scheme proposed by Pendry 121. The Green’s function is evaluated in the region of one rumpled layer in a sphericalwave expansion with basis functions centered at the atoms in

OO1O-4655/88/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

F. Mdca, M. Scheffler / Surface Green’s function

382

that layer. Beginning with the top (surface) layer this procedure is repeated layer-by-layer for the specified number of layers. When needed, the Green’s function for the hulk layers can also be evaluated. Restriction on the complexity of the problem All substrate layers are assumed to be identical. They may differ from the top layer. The structure of the top layer must he a superstructure of the substrate, including 1 >< 1. The maximum number of atoms per two-dimensional unit cell in one layer is restricted to four. The maximum array dimensions

are set for 45 plane waves and 3 phase shifts (9 spherical waves) per atom. Running time The time for the test run (one energy and one k~point) is S s on a CRAY-MP, using 45 plane waves and 3 phase shifts. References [1] K. Kambe. Z. Naturforsch. 23a (1968) 1280. [2] J.B. Pendry. Low Energy Electron Diffraction (Academic Press. London and New York, 1974).

LONG WRITE-UP 1. Introduction In this paper we describe a method to calculate the Green’s function of a surface with two-dimensional periodicity but rumpled geometry. The Green’s function is related directly to the local density of states and the electronic charge density. These quantities can be obtained for any layer of the semi-infinite crystal, and they are given in terms of a partial wave (angular quantum number) expansion. The program described below is a generalization of the computer code of paper I [1] to systems with rumpled atomic layers. It can handle a clean rumpled crystal surface, possibly formed by reconstruction, as well as an adsorbate system where the adsorbate—substrate distance is small. For such cases the multiple-scattering processes in the surface region cannot be accurateley described as a scattering between planar atomic layers in a plane-wave basis (bad convergence with respect to the number of plane waves). It is necessary to use a spherical-wave basis and the method of Ewald summation given by Kambe [2]. For example, for the adsorbate system S on Pd(100) it is more accurate to treat the S—Pd double layer as one “rumpled layer” than as two layers. Energy levels or resonances may then shift by some tenths of an electron volt (see also section 5) as a consequence of the higher accuracy of the present approach. Below, the term “layer” is always used for “rumpled layer”; such a layer may contain up to four different atoms per two-dimensional unit cell. It is to be noted that the bulk layers may also be rumpled. Time as well as memory are now increased compared to computer code I. For systems composed of planar layers with larger spacings the results naturally agree with those of computer code I. In the following section we summarize the method. The new subroutines, not presented in paper I, and the other alterations due to the generalization of the computer program to rumpled layers are described in section 3. The information about input parameters is given in section 4, and the results of test runs are presented in section 5.

2. Layer Green’s function The calculation of the layer- and k11resolved Green’s function proceeds in the same way as described in paper I (k1 is the two-dimensional wave vector). The semi-infinite system of muffin-tin type potentials is divided into rumpled layers parallel to the surface and the layer KKR method of Kambe and Scheffler [3] is applied. There are two successive steps: 1) The layer for which the Green’s function should be obtained is at first left out and the Green’s function, Gempty, for the system in which this layer is missing, is evaluated.

F. Maca, M. Scheffler / Surface Green’s function

383

2) The scattering at the left-out layer is then included, and gives the complete Green’s function valid in the region of the layer of interest. The empty-layer Green’s function Gempty ( r, r’; E, k1~)is evaluated in a spherical-wave basis centred at positions ‘7, of all the missing atoms inside the two-dimensional unit cell. The energy E refers to the muffin-tin zero. Following our previous notation [1] this reads: ~

r’; E,

kii)Gempty(r+T,,

=

r+’5 E, k11)

~

~

Im ,i’rn’

+ Ci’,~t’m’ jj(/~iKr)ji’(~5Jkr’)}Ytm(Q(r))Y~m’(Q(r’)),

(1)

where r’~i

—

I,n.l’m’

—

A’]

.L

tm,l’m’

~] /pn.l’m’

1 if i =1’ 1= 1’, m m’ and ~i/m,jt’m’ 0 otherwise. are KKR structure constants and B~,j’m’ are the scattering corrections due to the potential of the remaining layers. The definition spherical Bessel functions[4]. j~and functions h ~ as well as 2(r))agreeofwith Pendry’s convention The Hankel layer KKR structure constants sphericalare harmonics evaluatedY,m(~ from and

~j/m jt’m’

=

=

=

A~,i’m’

A7mi’m’=~E’(1m,~ where m”

m’

=

—

(3)

m and we introduced the Clebsch—Gordan-type constants

E’(lm, i’m’)

2±m+mfY,,

=

4~(_1)(t_I

~~(s2)1’,(~2)1’~~ ~(s2)dQ.

(4)

_I)/

Mathematical details and the derivation of the formulae for the calculation of the lattice sums D~can be found in Kambe’s original paper [2] (compare refs. [4,5] too). Our phase and spherical harmonics conventions differ from Kambe’s and we obtain:

=

~

—

[(2!

+ 1)(i +

m)!(i

—

m)!]1~2~exp( +i(k 11

t-lniI n=O

1

K~ ~

2n-1

+

g)(T,11

—

~)} exp{ —im~(g)}

Ik11+gI/%/~K) (i-m_s)!(i+m-s)I

M

[~(~_~)]2n_s(

Z~ng(a)~(~fl_s)

~ (5)

where A is the two-dimensional unit-cell area and M min(2n, i m I). p (g) are a.zimutal angles of (k11 + g), where the vectors g span the two-dimensional reciprocal net corresponding to the layer Bravais lattice given by translation vectors R~.~ and ~ are the parallel and normal components of the =

—

F. Maca, M. Scheffler / Surface Green ‘sfunction

384

three-dimensional vector ‘r,. For the definition of the integrals ~,g(a) see eq. (3.17) in ref. [2]. calculated by a recursion method. The expression for D~2>reads with our phase convention: D~2>=

L\,,.,~ are

~ xfu_3/2~ exp{u_

(~IR,

+ ~_~l/2)2/u}

du.

(6)

cs is the parameter that determines the division of the series between reciprocal (D~~)and real (D~2>) space. Kambe’s value a EA/4’rr is used in the code. The scattering corrections ~ in eq. (2) are obtained as described in paper I by transformation of the plane-wave expansion (compare eq. (11) in paper I and version 2). =

G° W~+

Genipty —

+

=

(7)

+

where G>~is the free-electron (vacuum) Green’s function and the scattering operators (right-hand side of eq. (7)) are defined in paper I. These operators represent the multiple scattering processes between the regions “below” (semi-infinite bulk) and “above” (surface area) the left-out layer. Using the definition of phase factors 4,(g) =exp(iK~~)=exp(—iK;.i~) we

(8)

obtain ~

=

(4~)2i’~’( 1)’”~ {~‘~~(g)/~,(g’) + (-1)’~’/~,(g)/~,(g’) -

+

( —1) ~

+

(— 1)~ “~“~‘/~,(g)~ 1(g’)}

xexp{i(K~ - K~.T,)}Y>~,(~(K;))~>,,,( Q(K;.)).

(9)

Kg denotes the normal component of three dimensional vector K which has been introduced by eq. (3) of paper I. There are no changes compared to code I in the part of the program which calculates the plane-wave scattering between layers.

3. Program structure The crystal is considered to be composed of identical (rumpled) atomic layers parallel to the surface. Only the top layer, possibly including adsorbate atoms, may be different from the other layers. The unit cell of any layer can include up to four different atoms. Only small alterations have been made in the main program. The subroutine CELMG which generates Clebsch—Gordan-type coefficients (compare eq. (4)) has been generalized (removing the restriction / + 1’ + m + m’ even). The subroutine XMATKC for evaluating the layer KKR structure constants as well as the subroutines where the layer transmission and reflection coefficients are computed (MSMOTB, MSMATB) have been rewritten completely. The generalization to rumpled substrate layers required also changes in RBULKA. =

F Maca, M. Scheffler / Surface Green’s function

385

Table 1 The input data for unit 2 Format

Variable

Comment

Control parameters and definition of the coordinate system 413 IPOT Different (= 0) or all identical (= 1) muffin-tin potentials in the total problem. ISTR Overlayer and substrate layers are identical (= 1) or different (= 0). IBULK Controls the calculation of the scattering matrices for the substrate layer > 2: The substrate geometry is read in only once; its unit cell must be identical to that of the adlayer. 2: The substrate geometry is read in twice: At first referring to the bulk unit cell and then referring to the adlayer unit cell. IARES If IARES = 1, the density of states is calculated for all different atoms in the rumpled-layer unit cell. If IARES = 0, the density of states is averaged over all atoms in the unit cell. F12.5 SPA Unit of length in atomic units used for some of the following input data.

3F12.5 3F12.5 3F12.5

SURF(3) ASM1(3) ASM2(3)

3F12.5

DIS(3, 1)

1A1

ICONTR

The layer geometry data 213 NATL NAT1

(3F12.5)

S(i, 3)

3F12.5 3F12.5

ARI1 ARI2

3F12.5

DIS(3, j)

The atomic potential data 3F12.5 Z(i) RMT(i) VSHIFT 416 NMESH KRMT NF ISHIFT

The following vectors SURF, ASM1, ASM2 are referred to the Cartesian coordinates z, x, y. SURF, ASM1, ASM2 must form a right-handed orthogonal system. The length of these vectors is not important, as they will be normalized in the program. The inwards pointing surface normal. A vector perpendicular to SURF, i.e. parallel to the surface. A vector perpendicular to SURF and ASM1 which completes the orthogonal right-handed coordinate system used in the program. The displacement of the first layer from the surface barrier. Because the barrier has no structure parallel to the surface, only DIS(1, 1) is relevant; referred to the direction SURF, ASMI, ASM2. The length is in units of SPA. A character which is read in twice at the beginning of the layer and atom data, and either twice or only once everytime at the end of the atom data. We call the first of the two ICONTR LAYER and the second ICONTR ATOM (see TEST RUN, INPUT FILE). 1) If both * E, then the layer and atom data follow. 2) If ICONTR LAYER * E but INCONTR ATOM = E, then only the atom data follow. 3) If ICONTR LAYER E, then no more layer or atom data follow. In this case the second ICONTR is no more read in.

Number of atoms of the layer unit cell. Number of different atoms in the layer unit cell. If NUMBER OF ATOMS specifications are NATL = 3 and NAT1 = 2 (or NATL = 4 and NAT1 = 3), the convention used in the program is such that the second atom is identical to the third. Position of the ith atom in the unit cell, referring to the directions SURF. ASMI, ASM2. The length is in units of SPA. A vector forming one side of the unit cell of the layer. A vector forming the other side of the unit cell of the layer. ARI1 and ARI2 are in units of SPA, referred to the coordinate axes z, x, y. The vector describing the displacement of the next layer deeper into the crystal. The three components refer to the directions SURF, ASM1, ASM2. The lengths are in units of SPA.

Atomic number of the ith atom. Muffin-tin radius. A constant shift of the muffin-tin potential. Number of points of the potential (Herman—Skillman-type grid). The first point beyond the muffin-tin radius. The point beyond which the grid interval stays constant. If > 0 then the potential is shifted by VSHIFT if =0 then VSHIFT is ignored.

386

F. Maca, M. Scheffler

/ Surface Green’s function

All arrays originally separated according to odd and even values of / + m have been reorganized. The important changes are in the subroutine MSBB where the empty Green’s function, Gempty, is calculated. The common block WS is used as work space to save computer memory. Most array dimensions depend on the input data (number of atoms in the unit cell, maximum number of beams, maximum number of phase shifts).

Table 1 (continued) Format

Variable

Comment

The atomic potential data (continued)

(5E14.6)

(RATP(J),ATP(J).J = 1,NMESH) The grid points and radial potential times r in hartrees. The last two radial point must he beyond RMT( 1). The muffin-tin potential must be tabulated on a Herman—Skiliman-type grid. The first radial point, RATP(1), is zero. Every successive group of ten points (e.g. RATP(J). J = 2 11) has equidistant spacing with an interval twice of that of the previous group. The doubling of the interval continues until J is greater or equal to NF. after which the interval is kept constant.

Control parameters

416

ISP IC1 lC2

IC4

Parameter controlling the integration of the radial Schrodinger equation. ISP is an input parameter for subroutine CALPFI (predictor—corrector method). Controls the amount of output on unit 6. Determines the type of the overlayer calculation: = 0: overlayer on substrate; = 1: overlayer in jellium; = 2: overlayer on jellium; = 3: isolated monolayer. Determines if the radial Schrodinger equation is solved, or if data are read in: = 0: calculate density of states; = 1: solve Schrodinger equation, write results on unit 4 and calculate density of states: = 2: the phaseshifts and radial itegrals are read in from unit 1: = 3: solve Schrodinger equation and write results on unit 4.

Calculation parameter

SF12.5

THETA,FI AK2,AK3

SF12.5

416

EOMEGA EVPR EVPI EEMIN, EESTEP, EEMAX NL LLFIXA LMAX N

The polar and azimuth angles (in degrees) specifying together with parameter EOMEGA (photon energy) the wave vector k 11. If THETA <0, then k1~is input directly by AK2, AK3 below. The components of the wave vector k~1referred to the adlayer reciprocal lattice vectors (used only if THETA <0). An energy (in eV) used to calculate the wave vector k11. It is only used if THETA 0. Real part of the inner potential (in eV). Imaginary part of the inner potential (in eV). Define the energies in eV for which the Green’s function should be calculated (these energies refer to the muffin-tin zero). Gives the number of layers for which the Green’s function should be calculated; if NL < 0 the values for the bulk are computed too. Number of layers used in RBULKA ( = 2, 4, 8, 16,...). Largest /-value in the spherical-wave expansion. Number of beams used in the plane-wave representation.

F. Maca, M. Scheffler / Surface Green’sfunction

387

4. Input data Meaning and sequence of the input parameters are practically the same as in paper I [1]. They are listed in table 1, together with short explanations. Few typing errors, present in table 1 of I, are now corrected. Concerning the sulphur and palladium potentials we give in the print out of the input file only the first and the last line (lines 190 to 610, 680 to 1220, 1380 to 1940 are deleted). The complete list of these potentials is, of course, present on the tape and it was listed in paper I.

5. Test runs The input data for the test runs are for the system (~/~ X I~)Son Pd(100). The relatively small distance between the top sulphur adlayer and the first palladium substrate layer requires more plane waves for an accurate description of the multiple scattering between these atomic planes. If the concept of a rumpled surface layer is introduced with three atoms in the unit cell (the sulphur atom and two palladium atoms of first substrate layer) the scattering problem can be handled properly in spherical waves. The results, i.e. the trace of the imaginary part of the Green’s function integrated over the muffin-tin spheres (cf. I, formula on p. 408) of the sulphur atom and of the palladium atoms in different atomic layers are shown in fig. 1. For comparison also bulk Pd is shown. The angular momentum resolved components and their sum are

I

100

~ 100

d

p

s I

I

I

I

I

I

x20

I

I~~’

I

TOTAL I

I

I

K100

x20

~ ,150

100

~

ENERGY BELOW E

1 (eV)

Fig. 1. The trace of the imaginary part of the Green’s function integrated over the muffin-tin spheres of a sulphur (top), a first palladium (second row) a second palladium (third row) and a bulk palladium layer (bottom) for = (0, 0). From the left to right: s-, p-, d-contributions and their sum are displayed.

388

F. Mdca, M. Scheffler / Surface Green’s function

presented. For the first Pd layer the two Pd atoms in the superstructure unit cell are equivalent. However, in the second Pd layer they are different. Therefore, we show a full and a broken line for these two atoms. 45 plane waves and 3 phase shifts have been used. We note that the sharp p-like states (induced by sulphur) are shifted by 0.2 eV toward lower energies in comparison to the results of code I, if both computations have been performed with 45 plane waves. This is a direct consequence of the better treatment of the substrate—adsorbate multiple scattering, now attained by the new code. In the ‘~/~J x adsorbate system, the two-dimensional primitive cell of the bulk layer is smaller than the adlayer primitive cell. As a consequence, the surface Brillouin zone of the adsorbate system (adlayer + substrate) represents a folded-back bulk layer zone. Therefore, the point k1~ (0, 0) in our adsorbate adlayer corresponds to two =

k11-points ((0, 0) and (2rr/a) (0.5, 0.5)) of the clean Pd(100) surface Brillouin zone, a is the nearest neighbour distance in the Pd crystal.

Acknowledgements We thank Ch. Droste for discussion. In particular we are gratefull to K. Kambe for discussions and helpful comments on the manuscript.

References [11 F. Maca and M. Scheffler, Comput. Phys. Commun. 38 (1985) 403, second version: 47 (1987) 349. [2] K. Kambe, Z. Naturforschg. 23a (1968) 1280. [31 K. Kambe and M. Scheffler, Surface Sci. 89 (1979) 262. [4) J.B. Pendry, Low Energy Electron Diffraction (Academic Press, London, New York, 1974). [5) A.P. Shen, Phys. Rev. B 3 (1971) 4200.

F. Maca, M. Scheffler / Surface Green’s function

TEST RUN INPUT FILE UNIT 1N2 10= 20= 30= 40= 50= 60= 70= 80= 90= 100= 110= 120=

0

0 2 1 7.42 1.0 0.0 0.0 0.355

3

2 0.0 0.331 0.331

130=

0.0

140= 150= 160=

111

0.0 0.0 1.0 0.0

0.0 —0.5 0.0

0.0 0.0 0.5 0.0 1.0 0.0 0.0

1.0

0.0 0.831 16.00000

170=

0.0 1.0 0.0 0.0

101

0.0 —1.0 1.7970

101

IPOT, ISTR,IBULK, IARES SPA SURF ASM1 ASM2 SURF BARRIER DISPLACEMENT ICONTR LAYER ICONTR ATOM NUMBER OF ATOMS 1. ATOM POSITION 2. ATOM POSITION 3. ATOM POSITION ARI1 ARI2 OhS Z,RMT,VSHIFT

0

180= 0.0 —0.160000E+02 0.175660E—03 190= 620= O.269641E+01 0.000000E+OO 630= 640=E 650= 46.0 2.62337 0.0 660= 141 131 131 0

670= 680=

0.0

—0.460000E+02

1230= 1240= 1250=

O.393523E+01

1260= 1270= 1280= 1290=

1

1300=

1310= 1320=

2

0.0

0.5

—0.5 —0.5

0.5 0.5 0.0

1 0.0

0.0

0.0

1330= 0.0 0.5 1340= 0.0 1.0 1350= 0.0 0.0 1360= 0.5 —0.5 1370= 46.0 2.62337 1380= 141 131 131 0 1390= 0.0 —O.460000E+02 1400= 1950= 0.3935235+01 0.000000E+00 1960=5 1970= ISP Id IC2 IC4 1980= 16 0 0 0

0.5 0.0 1.0 0.0 0.0

0.00 3.10000

1990=

—30.00000 12.00000

2010=

—1

2 02 0 = 202 1= XE OR

128

0.6405485—04

0.000000E+0O

0.0 0.5

2000=

O.351320E—03

ICONTR LAYER ICONTR ATOM Z,RMT,VSHIFT

0.320274E—04 —0.459931E+02

1

0.0 0.0 0.0

—0.159920E+02

45. —0.10000

2

45

ICONTR LAYER ICONTR ATOM NUMBER OF ATOMS 1. ATOM POSITION ARIl ARI2 DIS NUMBER OF ATOMS 1. ATOM POSITION 2. ATOM POSITION ARI1 ARI2 DIS Z,RMT,VSHIFT

O.320274E—04 —0.459931E+02

0.6405485—04

ICONTR

0.00 21.20000 0.05000 3.10000 NL,LLFIXA,LMAX,N

TRETA,FI.

389

F. Maca, M. Scheffler / Surface Green’sfunction

390

TEST RUN OUTPUT FILE UNIT 0U3 10= 20= 30= 40= 50= 60= 70= 80= 90= 100= 110= 120=

130= 140= 150= 160=

170= 180=

190= 200= 210=

220= 230= 240=

SPA

7.42000

SURF 1.00000 0.00000 0.00000 ASM1 0.00000 1.00000 0.00000 ASM2 0.00000 0.00000 1.00000 SURF. BARR. DISPL. 0.35500 0.00000 0.00000 LAYER 1 AOl 0.00000 1.00000 0.00000 AR2 0.00000 0.00000 1.00000 DIS 0.83100 —1.00000 0.00000 Z 16.0 RMT 1.797 ISHIFT 0 VSRIFT 0.00 Z 46.0 RMT 2.623 ISHIFT 0 VSHIFT 0.00 LAYER 2 AR1 0.00000 1.00000 0.00000 AR2 0.00000 0.00000 1.00000 DIS 0.50000 —0.50000 0.00000 Z 46.0 RMT 2.623 IsHIFT 0 VSHIFT 0.00 LAYER 1 NATMAX = 3 NAT1 = 2 0.00 0.00 0.00 0.33—0.50 0.00 0.33 0.00 0.50 LAYER 2 NATMAX = 1 NAT1 = 1 0.00 0.00 0.00 LAYER 3 NATMAX = 2 NAT1 = I 0.00 0.00 0.00 0.00 0.50 0.50 ISP Idi 1C2 1C4 IPOT ISTR IBULK IARES 16 0 0 0 0 0 2 1 NL,LLFIXA,LMAX,NX —1 128 2 45 VPR,VPI,OMEGA 12.00 —0.10 21.20 THrTA,FI,AK2,AK3 —30.00 45.00 0.00 0.00

250= 260= 1 45 O.310E+01—0.100E+OO 0.0005+00 O.000E+00 LAYER,NBEAMS.ENEItGY,K// 270= O.216E+OO—0.224E—O1 O.127E+O1—O.123E+01 0,448E—09—0.264E—08—0.280E—13 0.318E—13 280= 0.264E—08 O.448E—09—O.838E—O1—O.119E+00—0.448E—09 O.264E—08—O.171E—13—0.471E—15 290=—O.264E—08—O.448E—09 O.127E+O1—O.123E÷O1—O.804E+00—O.5895—O1—0.450E—11—0.423E—10 300= O.266E—15 0.5895—15 O.312E—1O—O.200E—10—O.377E—01 0.5265—01 0.423E—10—0.450E—11 310=—O.806E+OO—O.448E—02 O.127E—1O—O.103E—1O—O.117E—15 0.2465—16—O.312E—10 0.200E—10 320= 0.827E—16 0.9165—15 0.103E—10 O.127E—10—0.819E+OO—O.641E—02—O.1275—10 0.1035—10 330= O.309E—15—O.583E—15—O.200E—1O—O.3125—10 0.9545—16 0.729E—17—O.103E—10—0.127E—1O 340=—O. 806E+0O—O.448E—02 0.4505—11 O.423E—10—O. 3775—01 0. 526E—O1 0. 200E—1O 0.3125—10 350=—O.183E—15—O.531E—15—0.423E—10 O.450E—11—O.804E+00—0.589E—01 0.0005+00 O.000E+O0 360= O.000E+OO O.000E+OO

370= 380= 1 45 O.310E+O1—O.100E+O0 0.0005+00 0.0005+00 LAYER,NBEAMS,ENERGY,X// 390= 0.464E+00—O.1S1E+0O—0.541E+00—O.333E+O0 0.144E—09—O.593E—09—0.403E—O1 O.944E—01 400=—O.278E—10—0.286E—09—O.520E+00—O.250E—01—0.144E—09 0.5935—09—O.403E—01 0.9445—01 410= O.278E—10 O.286E—09—O.541E+00—0.333E+0O—0.115E+00—O.534E+O1 0.153E—07—0.322E—08 420= O.114E+01—O.135E+00—O.959E—09—O.133E—07—0.433E+00 0.494E+O1—0.720E—08—0.909E—08 430=—O.469E+O1—0.141E+O1 O.392E—08 0.292E—08—O.236E+O1—O.103E+01 0.959E—09 0.1335—07 440= O.114E+O1—O.135E+OO O.441E—08—0.155E—08—O.257E+01—O.441E÷0O—0.392E—O8—0.2915—08 450= O.114E+O1—O.135E+O0 O.128E—07 0.290E—08—O.236E+O1—0.103E+O1—0.441E—08 O.155E—08 460=—O.469E+O1—O.141E+O1—O.153E—07 O.322E—08—O.433E+OO 0.494E+01—O.128E—07—0.2905—08 470= 0.1145+01—0.1355+00 0.7205—08 O.909E—08—O.115E+00—O.534E+01 O.000E÷00 0.000E+00 480= 0.000E+OO 0.0005+00 490= 500= 1 45 0.3105+01—0.1005+00 0.000E+00 0.0005+00 LAYER,NBEAMS,ENERGY,K// 510= O.464E+0O—0.151E+OO—O.541E+OO—O.3335+00—O.286E—09 0.278E—1O 0.403E—O1—O.944E—01 520= 0.5935—09 0.144E—09—O.520E+O0—O.250E—01 O.286E—09—0.278E—10 0.403E—O1—O.944E—01 530=—O.593E—09—O.1445—09—0.541E+00—O.333E+0O—0.115E+OO—0.534E+01—0.909E—08 0.720E—08 540=—O.114E+01 O.135E-i.O0—O.290E—08 0.1285—07—0.4335+00 0.494E+01 0.3225—08 0.1535—07 550=—O.469E+fJl—O.1415+Ol—O.1555—o8—O.4415—08 O.236E+O1 0.1035+01 0.2905—08—0.1285—07 560=—O.114E+O1 0.1355+00—0.2925—08 O.392E—08—0.257E+O1—0.441E+0O 0.155E—08 O.441E—08 570=—0.114E+Ol 0.1355+00—0.1335—07 0.9595—09 0.236E+O1 0.103E+01 0.2925—08—0.3925—08 580=—O.469E+Ol—0.141E+01 0.9095—08—0.7205—08—0.4335+00 0.4945+01 0.133E—07—O.959E—09 590=—O.114E+O1 0.135E+0O—O.322E—08—O.153E—07—0.115E+OO—0.5345+01 0.000E+00 O.000E+00 600= 0.0005+00 0.0005+00 610=*EOR