- Email: [email protected]
Ion-Atom-Wave — Calculation of single ionization cross sections in ion—atom collisions
Computer Physics Communications ELSEVIER
Computer Physics Communications 114 (1998) 385--400
Ion-Atom-Wave
-
Calculation of single ionization cross sections in ion-atom collisions
B.S. Nesbitt 1, S.EC. O'Rourke, D.S.E Crothers The Department of Applied Mathematics & Theoretical Physics, The Queen's University of Belfast, Belfast BT7 1NN, N. Ireland Received 4 May 1998
Abstract A Fortran 90 program is presented, namely Ion-Atom-Wave ( l A W ) , which calculates total, single and double differential cross sections for the single ionization of atomic and molecular hydrogen-like targets and helium-like targets by both light-and-heavy-ion impact. The cross sections are evaluated within the continuum distorted-wave final-state ( C D W ) and the continuum distorted-wave eikonal-initial-state (CDW-EIS) approximations. IAW calculates cross sections within the wave treatment. This program extends a previously published program ION (D.S.E Crothers and M. McCartney, Comput. Phys. Commun. 72 (1992) 288) which calculated the total, single and double differential cross sections for the ionization of hydrogen-like targets by fully stripped ions. © 1998 Elsevier Science B.V. PACS: 34.50 Fa; 34.70+e Keywords: Single ionization; Total differential cross section; Single differential cross section; Double differential cross section; Continuum distorted wave eikonal-initial-state model; Hydrogen-like target; Helium-like target; Wave treatment
PROGRAM SUMMARY ~tle of program: Ion*Atom-Wave Catalogue identifier: ADJI Program Summary URL: http: //www.cpc.cs.qub.ac.uk/cpc/summaries / ADJI Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland Licensing provisions: none Computer for which the program is designed and others on which it is operable:
Computers: PC with four 200 MHz Pro Pentiums, DEC Alpha 4100, DEC Alpha 2000-500; Installations: PC and Dec Alpha 2000-500 at the Theoretical & Computational Research Division, Department of Applied Mathematics & Theoretical Physics, Queen's University of Belfast; DEC Alpha 4100 at the Computer Centre, Queen's University of Belfast Operating systems under which the program has been tested: Redhat Linux 4.2, Digital Open VMS 6.2, Digital UNIX 4.0a Programming language used: Fortran 90 with High Performance Fortran directives Memory required to execute with typical data: For test data, 30 Kb x four processors or 100 Kb for one processor
1 E-mail: [email protected] 0010-4655/98/$ - see front matter (~) 1998 Elsevier Science B.V. All fights reserved. PII S 0 0 1 0 - 4 6 5 5 ( 9 8 ) 0 0 0 7 2 - 1
386
B.S. Nesbitt et al./ Computer Physics Communications 114 (1998) 385-400
No. of bits in a word: 32 No. of processors used: Four No. of bytes in distributed program, including test data, etc.: 121376 Distribution format: uuencoded compressed tar file Keywords: Single ionization,total differentialcross section, single differential cross section, double differential cross section, continuum distorted wave eikonal-initial-statemodel, atomic/molecular hydrogen-like target, helium-liketarget, wave treatment Nature of physical problem The program calculates total, single and double differential cross sections for the single ionization of hydrogen-like atoms/molecules and helium-likeatoms by light-or-heavy-ionimpact, where the projectile ions are assumed to be structureless.
Method of solution The program allows the user to calculate cross sections using the CDW or CDW-EIS [ 1] approximationswithin the wave treatment. The cross sections are evaluated using the QUADPACKsubroutines QAGE, QKI5 and QK61 [2]. Typical running time Times vary according to input data. The average time taken to evaluate the double differential cross section for eight points of energy is 23 seconds, while using only one processor takes approximately 65 seconds. References 11] D.S.E Crothers, J.E McCann, J. Phys. B 16 (1983) 3229. [2] R. Piessens, E. de Doncker-Kapenga, C.W. Uberhuber, D.K. Kahane, QUADPACK, A Subroutine Package for Automatic Integration, Springer Series in Comput. Math., Vol. 1 (Springer, New York, 1983); QUADPACK is freely available at w ~ . n e t l i b . o r g or at the UK mirror, www.hensa.ac.uk/cgi-bin/cftp/unix. hensa, ac. uk/mirrors/netlib.
LONG WRITE-UP
1. Introduction The study of i o n - a t o m / m o l e c u l e ionization has provided a substantial part of our knowledge on the atomic structure and on the dynamics of such systems. Ion-matter interactions have many practical and important applications, namely radiation therapy, material modification, ion implementation and heavy ion pumped fusion. This area of research has involved many theoretical advances, details and further references can be found in Refs. [ 1,2]. We adopt the continuum-distorted-wave ( C D W ) and the continuum-distorted-wave final-state eikonal initialstate ( C D W - E I S ) approximations. The CDW-EIS model has been widely applied and has been shown to be versatile and successful in predicting many features of ionization at intermediate to high impact energies. The plan of this paper is as follows. In Section 2.1 we present the total, single and double differential cross sections for the single ionization of hydrogen-like atoms/molecules and helium-like atoms. In Section 2.2 we present the essential formulae for the CDW and CDW-EIS transition amplitude (a fuller description can be found in O ' R o u r k e and Crothers [ 3 ] ). Computational and numerical details are discussed in Section 3. Atomic units have been adopted throughout unless stated otherwise.
2. Theory In this paper we present a computer program which calculates the CDW and CDW-EIS total, single and double cross sections for the ionization process p Z ; + T ( Z r - 1 ) + ( I s ) ___~p Z ; + T Z 4 + e - ,
(1)
where P I T is the projectile/target of charge Z + / Z +. In particular, we consider atomic helium-like and atomic/molecular hydrogen-like targets.
B.S. Nesbitt et al./Computer Physics Communications 114 (1998) 385-400
387
2.1. Wave treatment of cross sections Within the wave treatment, the total cross section Q is commonly defined as
a = t kdk 0
ksinOdO 0
(2)
o'(k) do~, 0
where k is the momentum of the ejected electron (in the lab frame), 0 is the polar angle of k with respect to the polar axis v and ~b is the azimuthal angle of k with respect to the collision plane. The single differential cross section (SDCS) with respect to electron emission energy is defined as
dQ =
k sin OdO
dE~
0
o'(k) d~b,
(3)
0
where Ek is the kinetic energy of the ejected electron. The double differential cross section (DDCS) with respect to emitted electron energy and polar angle of emission is given as 2~
dQ
dEk d(cos 0)
= k/o(k)
d~b,
(4)
o
where ~r (k) denotes the triply differential cross section
cr(k) = 2~raZo T qdq Igif(q) l2 ,
(5)
~e
~lv
where ao is the Bohr radius and Ae = ½(k 2 + Z ~ ) ,
q=-
7/+---v-v
(6)
,
(7)
where v is the projectile velocity (in the lab frame), q is the transverse component of the change in relative momentum of the heavy particle and Rif denotes the scattering amplitude which we examine in the next section.
2.2. The scattering amplitudes The equations in Section 2.1 are quite general and are independent of the theoretical approximations used to calculate the scattering amplitude Rif. We now consider the CDW/CDW-EIS approximation to the scattering a m p l i t u d e R i f in the post form. For helium the scattering amplitude is defined as
[Rif(q)12= Z2Z21N(~)N(()N(P)I2A IA(I)2F~(iv;i(;
1;~-) -iuA(2)2Fl (iu + 1;i( + 1;2;r)l 2 ,
2~2u2a2T2
(8)
where Zr is the effective charge of the target. The three Sommerfeld parameters are defined by
= _L'r k'
(9)
B.S. Nesbitt et al./ Computer Physics Communications 114 (1998) 385-400
388
Zp
v = - - u,
(10)
(= Zv P
(11)
where p is the momentum of the ejected electron relative to the projectile and N(a) is the Coulomb densityof-the=continuum-states factor which is given by
N(a) = exp ( T r - ~ ) F ( 1 - i a ) .
(12)
The other physical quantities common to both CDW and CDW-EIS models are 5
A(I) = Z
baz3al2G( A )Ba '
(13)
baZ3/2G(h)BaI'23.'
(14)
.~=1 5
A(2) = Z 3.=1
a(a) =
(15)
o,,,T/33.
B3. = [q2 + (1 + ira)q" k ] , a~
1
2
(16)
2
~[Z3. + ( q + k)] ,
(17)
/33. = - [ q . k + k2(1 + i~3.) ] ,
(18)
(~
=
Za
(19)
k In particular, for the CDW model,
~CDW) = a(6B~ + ~,C~) yB,~(a + fl) A(cDw) = { 1 exp(--2~'v)
(20) (q2
_
(q2
_
2Ae > 0), 2Ae < 0),
(21)
and
r(CDW) _ fly -- at~
(22)
~,(a + / 3 ) '
where a , / 3 , y, 6 and Ca are common to both models and are defined as o~= lq2,
(23)
fl=q.v,
(24)
y=q.p+ot,
(25)
~3= p . v - pv + /3 ,
(26)
Ca= ( l + p )
[Ae-k'v(l
For the CDW-EIS model we have
+i('O] + v [ q ' k + k2(l
(27)
B.S. Nesbin et al./ Computer Physics Communications 114 (1998) 385--400
389
A (cDw-EIs) = e x p ( - 2 ~ p ) ,
(28)
f~(CDW-EIS) _
(29)
flyBa (8Ba + y C a ) ,
and 7"(CDw-EIS) = 1 -- Ct~
(30)
The coefficients Za and ba in the equations above come from the description of the Roothan-Hartree-Fock (RHF) wave functions which are given by ~bi(rr) = ~
baRa(rr)Ytm(O, (b) ,
(31)
where Ytm are spherical harmonics and R~(x) = ~
1
(2Z)n"+l/2Xn"-I e -z~x
(32)
The coefficients Za and ba are given in Clementi and Roetti [4].
3. P r o g r a m m i n g
and numerical details
3.1. Language and precision Ion-Atom-Wave is written in free-format Fortran 90 and incorporates a number of High Performance Fortran (HPF) directives. The program takes full advantage of Fortran 90's capabilities and as such the precision is not fixed; instead the programming module Precision defines the compile-time constants which determine the precision for real and complex variables (details are discussed below). This feature greatly enhances the portability of the code, enabling it to be compiled on anything from a supercomputer to a PC, with minimal fuss. Data is made globally accessible through module declarations, therefore there has been no need to use common blocks.
3.2. Input and output The program reads data from two input files and prints the results to two output files. The first input file is info_total.dat. This file may be set up by running the driver program Ion-Atom-Driver which asks the user interactively for all the necessary input. This data is not formatted and the structure of the file is as follows: target 0 = 1= 2 =
- Integer. Target used for collision. atomic hydrogen-like target. atomic helium-like target. molecular hydrogen-like target.
n o - Integer. 1 or 2 depending on chosen theory.
1 = CDW. 2 = CDW-EIS. my_section - Integer. Informs the program of the type of cross section required. 0 = total cross section.
390
B.S. Nesbitt et al./Computer Physics Communications 114 (1998) 385-400
1 = single differential cross section. 2 = double differential cross section. Z P - Real. The charge of the incoming projectile. Projectile_Energy_Start - Real. Starting value of the projectile energy in keV for the total cross section. The single and double differential cross sections are evaluated with respect to ejected electron energy and as such the projectile energy remains constant throughout. Projectile_Energy_Finish - Real. The for the total cross section since the to ejected electron energy, therefore program simply supplies a value of
finishing point of the projectile energy (keV). This value is only required single and double differential cross sections are evaluated with respect if my_section = 1,2, a dummy value must be specified; i.e., the driver 1.0.
- Real. Increment value used to determine the number of points of projectile energy. As stated above this only applies to total cross sections and a dummy value must be specified if
Projectile_Energy_Increment
my_section = 1,2.
Electron_Energy_Start - Real. Starting value of the ejected electron energy in eV. This value is not required for total cross sections (my_section = 0) although a dummy value is required. Electron_Energy_Finish - Real. The finishing point of ejected electron energy in eV. Again this is not required for total cross section although a dummy value must be given. Electron_Energy_Step - Real. Increment value used to determine the number of points of ejected electron energy. Again this is not required for total cross section although a dummy value must be given. a n g l e - Real. The polar angle of emission for the ejected electron, measured in radians. This value is only
required for the double differential cross section therefore if my_section = 0, 1 then a dummy value must be given. tol - Real. Global tolerance used throughout the code. Typical values may be given as 1.e-8, although one may run into problems if tol ~ e p m a c h (see below for definition of epmach). a c e - Real. Tolerance used for the nested integrals. This value is scaled by a factor of 0.1 as each new quadrature
routine is called. Therefore, for double differential cross sections, a typical value would be a c c = 0.0001. For single and total differential cross sections, reliable values are a c c = 0.001 and a c c = 0.01, respectively. The value of acc for total cross sections may appear to a little large for a tolerance, but anything lower may lead to computational difficulties, bearing in mind that we are dealing with four-fold quadrature. In the description of the input data above, it was stated that for certain cross sections dummy values were required. The program needs to read in a certain amount of data, the dummy values allow it to complete this task. As stated earlier, the program may evaluate cross sections for H, H2 and He-like systems. The user must provide an input file helium.dat which provides the RHF coefficients Za and ba (Eqs. (13) and ( 1 4 ) ) ; for He, the file contains the following values [4]: 5 1.41714
7.6837724e-1
2.37682 4.39628 6.52699 7.94252
2.2345803e-1 4.0818690e-2 -9.940800e-3 2.2963900e-3
B.S. Nesbittet al./ ComputerPhysicsCommunications114 (1998)385-400
391
The first value, 5, lets the program know how many Za and ba RHF coefficients need to be read in (this value is given by the upper limit in the summation of Eqs. (13) and (14)). For H and H2, the user must create a simple file hydrogen.dat which contains the following values:
5 1.0 0.0 0.0 0.0 0.0
1.0 0.0 0.0 0.0 0.0
For H, He and H2-1ike targets, the effective charges are defined in the code as ZT = 1.0_rap ZT = 2.0_rap* ( 2 4 . 5 8 / R y d b e r g ) * * 0 . 5 _ m p ZT = 1. 064_rap respectively. For H and H2-1ike targets, the final cross section is divided by two. Note that the factor of 2 arises from indistinguishability of the target's two electrons. Also, since our model is based on the independent electron model, we approximate the H2 target as an effective one-electron hydrogenic system, therefore Zr = 1.064. On completion, the program creates two output files; he_data.dat for detailed cross section information and total_output.dat which contains the raw data which is very useful for plotting graphs. Examples are given at the end of the paper.
3.3. Structure of program The code consists of five main programming units; the main program and four Fortran 90 modules organised as follows: Module
Precision ]
] Module Global-Data I Module
Complex_Functions I
I Module
Cross_Section I
Main Program
Ion_Atom ]
3.4. Module precision This module is the first stage in Ion-Atom-Wave and defines the compile time constants are given by mp = s e l e c t e d _ r e a l _ k i n d ( l O , 3 0 ) c_mp = m p
This defines the precision for real and complex variables using the extension real(kind--mp) :: x, y c o m p l e x ( k i n d = c mp) :: z
rnp and c_mp which
B.S. Nesbitt et al. / Computer Physics Communications 114 (1998) 385-400
392
also real constants are defined as 1.0rap instead of the usual 1.0 or 1.d0, for complex numbers the notation is z = cmplx(1.0_rap, t . 0rap, c_mp). These features greatly enhance the portability of the code. Module Precision contains one subroutine Machine_Accuracy which defines the following real constants:
epmach, uflow, tolerance. epmach - Real. Machine accuracy. uflow - Real. Smallest positive spacing between two real numbers. tolerance - Real. Machine dependent tolerance; i.e. t o l e r a n c e lower bound on tolerances required for numerical routines.
= 1.e+6 rap*epmach. This helps to set a
The constants mp, c_mp, epmach, uflow and tolerance are globally accessible.
3.5. Module Global_Data This module contains one subroutine namely Get_Data. Subroutine Get_Data is called once by the main program and reads the input data from the files info_totaLdat and hartree.dat associated with units getl = 8 and get2 = 9, respectively. The data discussed in Section 3.2 may then be accessed by other programming units by using the use Global_Data statement.
3.6. Module Complex_Functions This module contains five Fortran 90 real and complex programming functions. The module accesses the constants mp, c_mp, epmach and tolerance from module Precision. Sqmodn(x) - This function evaluates IN(a)l 2 using Eq. (6.1.31) of Ref. [6], where N(a) is given by Eq. (12).
function
Cgamma(z) - To evaluate the analytic continuations of the Gauss hypergeometric function in Eq. (8) we need to evaluate the complex gamma function F ( x + iy). This is best achieved using the Lanczos approximation [ 8 ].
function
Cdigam(z) - This function calculates the digamma function ~ ( z ) for complex argument z using an asymptotic expansion (Eq. (6.3.18) of Ref. [6])
function
1
g,(z) ~ lnz
l
2z
1
1
12z 2 + -120z ---q
252z 6 + . . .
(z --~ oo in larg(z)l < ~') •
(33)
If IRe(z)l < 15, the recurrence property of the digamma function (Eq. (6.3.6) of Ref. [ 6 ] ) , 1
~ b ( n + z ) - ( n - 1) + z
+
1
(n-2)
+z
+..-+
1
~
+
1
~
+~p(1 + z )
'
(34)
is used to retain accuracy when using the asymptotic expansion. If R e ( z ) < 0, the analytic continuation is used (Eq. (6.3.7) of Ref. [6]), ~O(1 - z) = ~/'(z) + ~ c o t ~ z .
(35
3.6.1. Gauss hypergeometric functions The most difficult special function to evaluate in the program is the Gauss hypergeometric function which is required for the evaluation of the scattering amplitude (Eq. ( 8 ) ) . Four programming units are used to evaluate 2F1 (a, b; c; z ); Feq, F_n, V1 and W1. The functions are required to evaluate the analytic continuations given below.
B.S. Nesbitt et al./ ComputerPhysics Communications114 (1998) 385--400
393
The Gauss hypergeometric series is defined as -~(a)n(b)nz
2Fl(a,b;c;z) =
n=O
n
(36)
(c)nnt
where the Pochhammer symbol is given by (a)n=a(a+l)(a+2)...(a÷n-1)
and
(a)0=l.
(37)
All the important properties of the Gauss hypergeometric function can be found in Ref. [6]. In particular, it is necessary to employ the following analytic conditions to enable numerical evaluation over the whole range of z ~ [-c~, +oo]: 2F1 (a, b; c; z ) = (1 - z ) - a 2 F l ( a , c - a ; c ; z / ( 1 - z ) )
([arg(1
2Ft(a,b;a+ b_c 2Fl(a,b;c;z) = F(c)F(c-a-b) F ( c - a) F ( c - b)
-z)l
<~r),
(38)
+ l ; l _ z ) + ( l _ z)C_a_b F ( C ) F ( a + b - c ) F ( a ) F ( b)
X2Fl(C-a,c-b;c-a-b+l;1-z) (I arg(1 - z)[ < cr) , F ( c ) F ( b - a) 2Fl(a,b;c;z) = ( - - z ) - a 2 F l ( a , 1 -- c + a; 1 -- b + a; 1 / z ) F ( b ) F ( c - a) l " ( c ) F ( a - b) 2Fl(b, 1 - c + b ; 1 -a+b;l/z) (larg(-z)l q-F(a)F(c-b)
(39)
<~),
(40)
and 2F1 (a, a; c; z ) =
r(c) F ( a ) F ( c - a)
oo
( _ z )_a ~-~ (a)n(1 -- c + a ) n z _ n
n--0
(n!)2
x [ l n ( - - Z ) + 2~b(n+ 1) -- ~ , ( a + n) -- ~ b ( c - a - n)]
(larg(--Z) < 7r),
(41)
as required in Eq. (8). The program unit V1 evaluates 2Fl(a,b; 1;z), unit W1 evaluates 2Ft(1 + ia, 1 + i b ; 2 ; z ) and the unit Feq evaluates 2 F l ( a - 1 + ib, a - 1 + i b ; b ; z ) , where a, b and z are real. More generally, the program unit F_n evaluates 2 F l ( a , b ; c ; z ) for complex a, b and c and real z provided that Re(c - a - b) > 0 and - 1 < z < 1. It is possible to remove the restrictions on Eqs. ( 3 8 ) - ( 4 1 ) provided it can be determined whether 1 - z and - z approach the negative real axis from above or below; i.e., whether we should write - 1 as e iTr or e -i~r. To determine this we note that the Gauss hypergeometric functions which appear in the CDW and CDW-EIS models occur as solutions of the standard Nordsieck integral [7], /
d~3rrexp(
~ ~r
+ iq- r ) I F I ( - i a l ; 1; i ( p l r + Pl" r) ) 1FI (ia2; 1; i(p2r + P2" r) )
271" ( O l 8 =
tr----~ ots--+fl8
)ia2 (OtS~ tal \Ys/
2FL(--ial,ia2; 1;7") ,
(42)
where
7"----and
flsV8 - ~ 8 8 8
(43)
394
B.S. Nesbitt et al. / Computer Physics Communications 114 (1998) 385~100
c~ = ½(q2 + h2) , Y8 = P l " q -- i a p j
/3 = P2" q -- i ap2, + c~,~,
68
=
Pl
"
P2 +/33-
(44) (45)
In general it is necessary to assume, for convergence, that Re(A) > 0. Clearly for a 4: 0, Eq. (43) has an imaginary part and the sign of this imaginary part determines whether the real axis is being approached from above or below as a --~ 0 +. For C D W theory as A ~ 0 + it is easy to show that the sign of I m ( r ) is given by the sign of - a [ ~ 3 p ( a + / 3 ) + y v ( ~ + Y - ce - / 3 ) ] .
(46)
Note that setting Pl = P, P2 = v and taking ,~ --~ 0 + in Eqs (44) and (45) yield the corresponding quantities in Eqs. ( 2 3 ) - ( 2 6 ) for C D W theory. Similarly, for CDW-EIS theory it can be shown that as A ~ 0 + then the sign of Im(~-) is given by the sign of - a [ 6 p f l + v y ( 6 + y) ] .
(47)
It is of great importance to have a test which is robust computationally and avoids division by zero. Eqs. (46) and (47) provide the definition of the variable TEST in function CS(x) in module Cross_Section. 3. 7. Module Cross_Section
Module Cross_Section uses modules Precision, Global_Data and Complex_Functions and contains subroutines and functions used to evaluate the cross sections. The following real variables are declared public and are therefore globally accessible: Z T - effective charge of the target. velocity - velocity v of the projectile in a.u. g l o b a l & - momentum of ejected electron k. The following real variables are declared private and can only be accessed by subroutines and functions contained with module Cross_Section. q_min - lower limit of the integral in Eq. (5). delta_epsilon
Ae given by Eq.
(6).
( given by Eq. (9).
xi nu
-
- u given by Eq. (10).
zeta - ~" given by Eq. (11). smnnu smnxi
-
N(v) given by Eq.
(12).
- N(().
snt
- sin 0.
cst
-
cosO.
k_dot_v - the dot product k • v. p_dot_v - the dot product p - v. AB - ] N ( u ) N ( ( ) N ( ( ) I .
B.S. Nesbitt et aL/ Computer Physics Communications 114 (1998) 385-400
395
csp - cos ~b.
The total cross sections are obtained by evaluating four fold nested integrals. Obviously we require extremely robust and accurate quadrature rules; additionally they must be efficient to evaluate cross sections within a reasonable time. The nested integrals are evaluated using a set of Gauss-Kronrod rules. These subroutines originate from the highly acclaimed set of automatic integrators known as QUADPACK [ 5 ]. For our purposes we have modified the QAGE, QK15 and QK61 subroutines and converted them to Fortran 90. In general the 61-point Gauss-Kronrod rule is used throughout the calculations except for peak problems (such as the ECC peak, defined below) where the 15-point rule is preferred. The subroutines and functions used to evaluate the cross sections and set up the quadrature are as follows: subroutine Totai_Differential(Modei_Energy, TDCS) - This subroutine is called by the main program and evaluates the outer integral over k where T D C S = total cross section if my_section = O. = single differential cross section if my_section = 1. = double differential cross section if my_section = 2. The Model_Energy can be either the projectile energy or the energy of the ejected electron, depending on the type of cross section to be evaluated. For single and double differential cross sections, Model_Energy :-- Electron_Energy then by default the projectile energy is given as Projectile_Energy_Start. Alternatively, for total cross sections Model_Energy -- Projectile_Energy. The function DDl(dummy_k) is then called. For single and double differential cross sections this function is called for specific values of k (the momentum of the ejected electron). For total cross sections subroutine QUAD4 is called to integrate the function DD1 over a range of k. In theory the outer integral should be evaluated from 0 to cx~, in practice, however, this is never necessary. Instead we truncate the range of integration to k_max, where k_max is related to the energy of the incoming ion Ek by k_max =
,-~ 3.5 x v,
(48)
which ensures that contributions to the total cross section from ionized electrons up to and including the Bethe ridge are included. Also, the point in momentum space where k = v has a simple pole in the DDCS, which is commonly referred to as the ECC peak (electron capture to the continuum'. The pole is due to go = Z J p , where p2 = k 2 + v2 _ 2 k - v ,
(49)
which only occurs in the forward direction (0 = 0) and where p is defined in function DD2. We avoid this point in the following manner; if I v - k J < • t h e n if v < k then v=v-~v else v=v+6v endif endif
where • = tolerance and ~v > •. This allows the singularity to be jumped over without losing vital information.
396
B.S. Nesbitt et al. /Computer Physics Communications 114 (1998) 385-400
function DDl(dummy_k) - If my_section = 0 or 1 (i.e., evaluation of total or single differential cross sections), DD1 calls subroutine QUAD3 to integrate the function DD2 over a range 0 to ¢r for a given value of k. If my_section = 2 (double differential cross section) then DD1 evaluates DD2 at a particular value angle which is read from module Global_Data. function DD2(polar) - This function calls subroutine QUAD2 to integrate function DD3 over ~b for specific values of k and 0. The integrand of the integral over ~b, the azimuthal angle of k with respect to the collision plane only depends on cos ~b, therefore the transformation f f ( c o s O) 0
f(-
d~b =
-,r
f ( - cosO) d~
(50)
0
can be made which halves the range of integration which in turn halves the time taken to evaluate this integral.
function DD3(polar) - Function DD3 calls QUAD1 to integrate CS over the range - 1 to +1 for a given value of ~b, O and k. function CS(x) - Function CS is the central integrand defined by Eq. (5). The range of the integration is changed from 7/ E [0, o~] to x E [ - 1 , +1] by the transformation r/2= 1 - x 1+ x
(51)
This function also evaluates the variable TEST which determines which of the four analytic continuations to use (Eqs. ( 3 8 ) - ( 4 1 ) ) as described in Section 3.6. For the CDW model TEST is defined by Eq. (46), for the CDW-EIS model it is defined by Eq. (47). Module Cross_Section also contains the functions V1, W1 (which where discussed in Section 3.6) and the QUADPACK routines QUAD1-QUAD4.
3.8. Main program Ion_Atom The main program uses modules Precision, Global.Data and Cross_Section. The first executable statement is to call subroutine Machine_Accuracy which determines the machine dependent constants. The next statement calls subroutine Global_Data which reads the input data from file. This data can then be accessed by the other programming units. The output files he_data.dat and total_info.dat, associated with units 6 and 4, respectively, are then initialised.
3.8.1. Parallelization of code The main program incorporates a number of High Performance Fortran (HPF) directives allowing large sections of the code to be run concurrently in a SIMD (single instruction multiple data) fashion. This is achieved as follows. Consider the total cross section, the user initially supplies starting and ending values for the kinetic energy of the projectile, Projectile_Energy_Start and Projectile_Energy_Finish, respectively. An increment value is also supplied, Projectile_Energy_Increment, enabling the program to create an allocatable array Energy of energy points which contains N elements where the number of elements is determined from N = Etinish - - E s t a r t + 1. Eincrement
The program is then parallelised as follows:
(52)
B.S. Nesbitt et al. / Computer Physics Cormnunications 114 (1998) 385--400
397
!HPF£ PROCESSORS procs(number_of_processors) !HPF£ DISTRIBUTE (BLOCK) ONTO proes :: Energy, Cross !HPF£ INDEPENDENT do loop = 1,N call Total_Differential(Energy(loop), Differential_Cross_Section) Cross(loop) = Differential_Cross_Section end do where the do loop evaluates the cross section N times for N different values of energy. For single and double differential cross sections, projectile energy is replaced with ejected electron energy. Note that each iteration of the do loop is completely independent of all others enabling the HPF directive to allocate the N calls to the subroutine Total_Differential to an array of processors.
Acknowledgements Brian Nesbitt acknowledges financial support from the Department of Education, Northern Ireland (DENI). Francesca O'Rourke acknowledges current financial support by a Royal Society University Research Fellowship and Derrick Crothers acknowledges EPSRC support under grant GR/L21891.
References [1] [2] [3] [4] [5]
D.S.E Crothers, J. E McCann, J. Phys. B 16 (1983) 3229. P.D. Fainstein, V.H. Ponce, R.D. Rivarola, J. Phys. B 24 (1991) 3091. S.EC. O'Rourke, D.S.E Crothers, J. Phys. B 30 (1997) 2443. E. Clementi, C. Roetti, At. Dam Nucl. Tables 14 (1974) 177. R. Piessens, E. de Doncker-Kapenga, C.W. Uberhuber, D.K. gahane, QUADPACK, A Subroutine Package for Automatic Integration, Springer Series in Comp. Math., Vol. i (Springer, New York, 1983); QUADPACK is freely available at ~ n ~ . n e t l i b . o r g or at the UK mirror, ~ r m ~ . h e n s a . a c . u k / c g i - b i n / c f t p / u n i x . hensa, ae. uk/mirrors/net lib~
[6] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions (Dover, New York, 1970). [71 A. Nordsieck, Phys. Rev. 93 (1954) 785. [8] C. Lanczos, J.S.I.A.M. Num. Analy. B 1 (1964) 86.
B.S. Nesbitt et al./Computer Physics Communications 114 (1998)385-400
398
TEST RUN OUTPUT Input Data Test data for file info_total.dat generated by the program Ion-Atom-Driver: 1 2 2 I. 000000000000000 75. 00000000000000 1. 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 i. 000000000000000 O. 2000000000000000 90. 00000000000000 2. 000000000000000 O. i000000000000000 1. O000000000000000E-O i0 1. O000000000000000E-O04
Test data for file helium.dat: 5 1.41714D0 7. 6837724D-01 2. 37682D0 2. 2 3 4 5 8 0 3 D - 0 1 4 . 3 9 6 2 8 D 0 4. 0 8 1 8 6 9 0 D - 0 2 6.52699D0 -9.940800D-03 ?.94252D0 2.2963900D-03
Test data for file hydrogen.dat: 5 1.0 0.0 0.0 0.0 0.0
1.0 0.0 0.0 0.0 0.0
Output Data Ion-Atom-Wave
Ion-Atom-Wave calculates cross sections for the single ionization of atomic and molecular hydrogen-like targets and helium-like targets by both light and heavy ion impact.
Input .
.
.
.
.
.
data .
.
.
.
.
used for .
.
.
.
.
.
.
.
.
.
this .
.
.
.
.
run .
.
.
:
.
The target is a helium atom The CDW-EIS model was used. double differential cross section evaluated for the following values :
B.S. Nesbitt et at~Computer Physics Communications 114 (1998) 385-400 Program started with the following electron energy (eV) : Program ended with the following electron energy (eV) : The electron energy was incremented by Energy of the incoming ion (keY) Angle of electron emission (in radians)
0.20000000E+O0 0.90000000E+02 0.20000000E+01 0.75000000E+02 0.10000000E+O0
Numerical details
Tolerance used throughout the code : 0.22204460E-09 Tolerance used fo r the quadrature : 0.10000000E-03
PROGRAM OUTPUT
The first column contains the electron energy (eV) and the second column contains the cross sections
0.200000000000000 2.200000000000000 4.200000000000000 6.200000000000000 8.199999999999999 10.20000000000000 12.20000000000000 14.20000000000000 16.20000000000000 18.20000000000000 20.20000000000000 22.20000000000000 24.20000000000000 26.20000000000000 28.20000000000000 30.20000000000000 32.20000000000000 34.20000000000000 36.20000000000000 38.20000000000000
40.20000000000000 42.20000000000000 44.20000000000000 46.20000000000000 48.20000000000000 50.20000000000000 52.20000000000000 54.20000000000000 56.20000000000000
2.3394624948705482E-OO2*le-20sqr.metres/eV/sr 3.0481267769151714E-OO2*le-20sqr.metres/eV/sr 3.2186809295606460E-OO2*le-20sqr.metres/eV/sr 8.2579413664121801E-OO2*le-20sqr.metres/eV/sr 3.2480520852528200E-002*le-20sqr.metres/eV/sr 3.2224803008163452E-OO2*le-20sqr.metres/eV/sr 3.1975245277550073E-O02*le-20sqr.metres/eV/sr
3.1824782971607055E-OO2*le-20sqr.metres/eV/sr 3.1830796483122169E-OO2*1e-20sqr.metres/eV/sr 3.2038672844994195E-OO2*le-20sqr.metres/eV/sr 3.2481768424916561E-002*le-20sqr.metres/eV/sr 3.3207862743912936E-O02*le-20sqr.metres/eV/sr
3.4263373599848164E-OO2*1e-20sqr.metres/eV/sr 3.5705354273153703E-OO2*le-20sqr.metres/eV/sr 3.7597027469323460E-OO2*le-20sqr.metres/eV/sr
3.9994014059866678E-OO2*le-20sqr.metres/eV/sr 4.2902732491859850E-OO2*le-20sqr.metres/eV/sr
4.6179049743754010E-OO2*le-20sqr.metres/eV/sr 4.9337403545042326E-OO2*le-20sqr.metres/eV/sr 5.1366182908066736E-OO2*le-20sqr.metres/eV/sr 5.0990543575597416E-OO2*le-20sqr.metres/eV/sr
4.7703676460025074E-OO2*1e-20sqr.metres/eV/sr 4.2414824071991944E-OO2*1e-20sqr.metres/eV/sr 3.6591704956531815E-002*le-20sqr.metres/eV/sr 3.1219229408178922E-O02*le-20sqr.metres/eV/sr
2.6653835839971593E-OO2*1e-20sqr.metres/eV/sr 2.2902354699969769E-OO2*le-20sqr.metres/eV/sr 1.9848356642432175E-OO2*le-20sqr.metres/eV/sr
1.7356716559748971E-OO2*le-20sqr.metres/eV/sr
58.20000000000000
1.5808863641622825E-OO2*le-20sqr.metres/eV/sr
60.20000000000000 62.20000000000000 64.20000000000000 66.20000000000000 68.20000000000000
1.3609851200881467E-002*le-20sqr.metres/eV/sr
70.20000000000000 72.20000000000000 74.20000000000000
1.2186184246292624E-OO2*le-20sqr.metres/eV/sr 1.0981554216192674E-OO2*le-20sqr.metres/eV/sr 9.9527930646778721E-OO3*le-20sqr.metres/eV/sr 9.0666130655147580E-OO3*le-20sqr.metres/eV/sr 8.2971324892396683E-OO3*ie-20sqr.metres/eV/sr 7.6240400385084737E-OO3*1e-20sqr.metres/eV/sr 7.0312447236584117E-OO3*le-20sqr.metres/eV/sr
399
B.S. Nesbitt et aL / Computer Physics Communications 114 (1998) 385-400
400
76.20000000000000 78.20000000000000 80.20000000000000 82.20000000000000 84.20000000000000 86.20000000000000 88.20000000000000
6.5058772606735949E-OO3*le-20sqr.metres/eV/sr 6.0375535399216493E-003*le-20sqr.metres/eV/sr
5.6178225246002788E-OO3*le-20sqr.metres/eV/sr 5.2397512887970713E-OO3*ie-20sqr.metres/eV/sr 4.8976089630747618E-003*le-20sqr.metres/eV/sr 4.5866258592475269E-003*le-20sqr.metres/eV/sr 4.3028056253776375E-003*le-20sqr.metres/eV/sr
Computer Physics Communications 114 (1998) 385--400
Ion-Atom-Wave
-
Calculation of single ionization cross sections in ion-atom collisions
B.S. Nesbitt 1, S.EC. O'Rourke, D.S.E Crothers The Department of Applied Mathematics & Theoretical Physics, The Queen's University of Belfast, Belfast BT7 1NN, N. Ireland Received 4 May 1998
Abstract A Fortran 90 program is presented, namely Ion-Atom-Wave ( l A W ) , which calculates total, single and double differential cross sections for the single ionization of atomic and molecular hydrogen-like targets and helium-like targets by both light-and-heavy-ion impact. The cross sections are evaluated within the continuum distorted-wave final-state ( C D W ) and the continuum distorted-wave eikonal-initial-state (CDW-EIS) approximations. IAW calculates cross sections within the wave treatment. This program extends a previously published program ION (D.S.E Crothers and M. McCartney, Comput. Phys. Commun. 72 (1992) 288) which calculated the total, single and double differential cross sections for the ionization of hydrogen-like targets by fully stripped ions. © 1998 Elsevier Science B.V. PACS: 34.50 Fa; 34.70+e Keywords: Single ionization; Total differential cross section; Single differential cross section; Double differential cross section; Continuum distorted wave eikonal-initial-state model; Hydrogen-like target; Helium-like target; Wave treatment
PROGRAM SUMMARY ~tle of program: Ion*Atom-Wave Catalogue identifier: ADJI Program Summary URL: http: //www.cpc.cs.qub.ac.uk/cpc/summaries / ADJI Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland Licensing provisions: none Computer for which the program is designed and others on which it is operable:
Computers: PC with four 200 MHz Pro Pentiums, DEC Alpha 4100, DEC Alpha 2000-500; Installations: PC and Dec Alpha 2000-500 at the Theoretical & Computational Research Division, Department of Applied Mathematics & Theoretical Physics, Queen's University of Belfast; DEC Alpha 4100 at the Computer Centre, Queen's University of Belfast Operating systems under which the program has been tested: Redhat Linux 4.2, Digital Open VMS 6.2, Digital UNIX 4.0a Programming language used: Fortran 90 with High Performance Fortran directives Memory required to execute with typical data: For test data, 30 Kb x four processors or 100 Kb for one processor
1 E-mail: [email protected] 0010-4655/98/$ - see front matter (~) 1998 Elsevier Science B.V. All fights reserved. PII S 0 0 1 0 - 4 6 5 5 ( 9 8 ) 0 0 0 7 2 - 1
386
B.S. Nesbitt et al./ Computer Physics Communications 114 (1998) 385-400
No. of bits in a word: 32 No. of processors used: Four No. of bytes in distributed program, including test data, etc.: 121376 Distribution format: uuencoded compressed tar file Keywords: Single ionization,total differentialcross section, single differential cross section, double differential cross section, continuum distorted wave eikonal-initial-statemodel, atomic/molecular hydrogen-like target, helium-liketarget, wave treatment Nature of physical problem The program calculates total, single and double differential cross sections for the single ionization of hydrogen-like atoms/molecules and helium-likeatoms by light-or-heavy-ionimpact, where the projectile ions are assumed to be structureless.
Method of solution The program allows the user to calculate cross sections using the CDW or CDW-EIS [ 1] approximationswithin the wave treatment. The cross sections are evaluated using the QUADPACKsubroutines QAGE, QKI5 and QK61 [2]. Typical running time Times vary according to input data. The average time taken to evaluate the double differential cross section for eight points of energy is 23 seconds, while using only one processor takes approximately 65 seconds. References 11] D.S.E Crothers, J.E McCann, J. Phys. B 16 (1983) 3229. [2] R. Piessens, E. de Doncker-Kapenga, C.W. Uberhuber, D.K. Kahane, QUADPACK, A Subroutine Package for Automatic Integration, Springer Series in Comput. Math., Vol. 1 (Springer, New York, 1983); QUADPACK is freely available at w ~ . n e t l i b . o r g or at the UK mirror, www.hensa.ac.uk/cgi-bin/cftp/unix. hensa, ac. uk/mirrors/netlib.
LONG WRITE-UP
1. Introduction The study of i o n - a t o m / m o l e c u l e ionization has provided a substantial part of our knowledge on the atomic structure and on the dynamics of such systems. Ion-matter interactions have many practical and important applications, namely radiation therapy, material modification, ion implementation and heavy ion pumped fusion. This area of research has involved many theoretical advances, details and further references can be found in Refs. [ 1,2]. We adopt the continuum-distorted-wave ( C D W ) and the continuum-distorted-wave final-state eikonal initialstate ( C D W - E I S ) approximations. The CDW-EIS model has been widely applied and has been shown to be versatile and successful in predicting many features of ionization at intermediate to high impact energies. The plan of this paper is as follows. In Section 2.1 we present the total, single and double differential cross sections for the single ionization of hydrogen-like atoms/molecules and helium-like atoms. In Section 2.2 we present the essential formulae for the CDW and CDW-EIS transition amplitude (a fuller description can be found in O ' R o u r k e and Crothers [ 3 ] ). Computational and numerical details are discussed in Section 3. Atomic units have been adopted throughout unless stated otherwise.
2. Theory In this paper we present a computer program which calculates the CDW and CDW-EIS total, single and double cross sections for the ionization process p Z ; + T ( Z r - 1 ) + ( I s ) ___~p Z ; + T Z 4 + e - ,
(1)
where P I T is the projectile/target of charge Z + / Z +. In particular, we consider atomic helium-like and atomic/molecular hydrogen-like targets.
B.S. Nesbitt et al./Computer Physics Communications 114 (1998) 385-400
387
2.1. Wave treatment of cross sections Within the wave treatment, the total cross section Q is commonly defined as
a = t kdk 0
ksinOdO 0
(2)
o'(k) do~, 0
where k is the momentum of the ejected electron (in the lab frame), 0 is the polar angle of k with respect to the polar axis v and ~b is the azimuthal angle of k with respect to the collision plane. The single differential cross section (SDCS) with respect to electron emission energy is defined as
dQ =
k sin OdO
dE~
0
o'(k) d~b,
(3)
0
where Ek is the kinetic energy of the ejected electron. The double differential cross section (DDCS) with respect to emitted electron energy and polar angle of emission is given as 2~
dQ
dEk d(cos 0)
= k/o(k)
d~b,
(4)
o
where ~r (k) denotes the triply differential cross section
cr(k) = 2~raZo T qdq Igif(q) l2 ,
(5)
~e
~lv
where ao is the Bohr radius and Ae = ½(k 2 + Z ~ ) ,
q=-
7/+---v-v
(6)
,
(7)
where v is the projectile velocity (in the lab frame), q is the transverse component of the change in relative momentum of the heavy particle and Rif denotes the scattering amplitude which we examine in the next section.
2.2. The scattering amplitudes The equations in Section 2.1 are quite general and are independent of the theoretical approximations used to calculate the scattering amplitude Rif. We now consider the CDW/CDW-EIS approximation to the scattering a m p l i t u d e R i f in the post form. For helium the scattering amplitude is defined as
[Rif(q)12= Z2Z21N(~)N(()N(P)I2A IA(I)2F~(iv;i(;
1;~-) -iuA(2)2Fl (iu + 1;i( + 1;2;r)l 2 ,
2~2u2a2T2
(8)
where Zr is the effective charge of the target. The three Sommerfeld parameters are defined by
= _L'r k'
(9)
B.S. Nesbitt et al./ Computer Physics Communications 114 (1998) 385-400
388
Zp
v = - - u,
(10)
(= Zv P
(11)
where p is the momentum of the ejected electron relative to the projectile and N(a) is the Coulomb densityof-the=continuum-states factor which is given by
N(a) = exp ( T r - ~ ) F ( 1 - i a ) .
(12)
The other physical quantities common to both CDW and CDW-EIS models are 5
A(I) = Z
baz3al2G( A )Ba '
(13)
baZ3/2G(h)BaI'23.'
(14)
.~=1 5
A(2) = Z 3.=1
a(a) =
(15)
o,,,T/33.
B3. = [q2 + (1 + ira)q" k ] , a~
1
2
(16)
2
~[Z3. + ( q + k)] ,
(17)
/33. = - [ q . k + k2(1 + i~3.) ] ,
(18)
(~
=
Za
(19)
k In particular, for the CDW model,
~CDW) = a(6B~ + ~,C~) yB,~(a + fl) A(cDw) = { 1 exp(--2~'v)
(20) (q2
_
(q2
_
2Ae > 0), 2Ae < 0),
(21)
and
r(CDW) _ fly -- at~
(22)
~,(a + / 3 ) '
where a , / 3 , y, 6 and Ca are common to both models and are defined as o~= lq2,
(23)
fl=q.v,
(24)
y=q.p+ot,
(25)
~3= p . v - pv + /3 ,
(26)
Ca= ( l + p )
[Ae-k'v(l
For the CDW-EIS model we have
+i('O] + v [ q ' k + k2(l
(27)
B.S. Nesbin et al./ Computer Physics Communications 114 (1998) 385--400
389
A (cDw-EIs) = e x p ( - 2 ~ p ) ,
(28)
f~(CDW-EIS) _
(29)
flyBa (8Ba + y C a ) ,
and 7"(CDw-EIS) = 1 -- Ct~
(30)
The coefficients Za and ba in the equations above come from the description of the Roothan-Hartree-Fock (RHF) wave functions which are given by ~bi(rr) = ~
baRa(rr)Ytm(O, (b) ,
(31)
where Ytm are spherical harmonics and R~(x) = ~
1
(2Z)n"+l/2Xn"-I e -z~x
(32)
The coefficients Za and ba are given in Clementi and Roetti [4].
3. P r o g r a m m i n g
and numerical details
3.1. Language and precision Ion-Atom-Wave is written in free-format Fortran 90 and incorporates a number of High Performance Fortran (HPF) directives. The program takes full advantage of Fortran 90's capabilities and as such the precision is not fixed; instead the programming module Precision defines the compile-time constants which determine the precision for real and complex variables (details are discussed below). This feature greatly enhances the portability of the code, enabling it to be compiled on anything from a supercomputer to a PC, with minimal fuss. Data is made globally accessible through module declarations, therefore there has been no need to use common blocks.
3.2. Input and output The program reads data from two input files and prints the results to two output files. The first input file is info_total.dat. This file may be set up by running the driver program Ion-Atom-Driver which asks the user interactively for all the necessary input. This data is not formatted and the structure of the file is as follows: target 0 = 1= 2 =
- Integer. Target used for collision. atomic hydrogen-like target. atomic helium-like target. molecular hydrogen-like target.
n o - Integer. 1 or 2 depending on chosen theory.
1 = CDW. 2 = CDW-EIS. my_section - Integer. Informs the program of the type of cross section required. 0 = total cross section.
390
B.S. Nesbitt et al./Computer Physics Communications 114 (1998) 385-400
1 = single differential cross section. 2 = double differential cross section. Z P - Real. The charge of the incoming projectile. Projectile_Energy_Start - Real. Starting value of the projectile energy in keV for the total cross section. The single and double differential cross sections are evaluated with respect to ejected electron energy and as such the projectile energy remains constant throughout. Projectile_Energy_Finish - Real. The for the total cross section since the to ejected electron energy, therefore program simply supplies a value of
finishing point of the projectile energy (keV). This value is only required single and double differential cross sections are evaluated with respect if my_section = 1,2, a dummy value must be specified; i.e., the driver 1.0.
- Real. Increment value used to determine the number of points of projectile energy. As stated above this only applies to total cross sections and a dummy value must be specified if
Projectile_Energy_Increment
my_section = 1,2.
Electron_Energy_Start - Real. Starting value of the ejected electron energy in eV. This value is not required for total cross sections (my_section = 0) although a dummy value is required. Electron_Energy_Finish - Real. The finishing point of ejected electron energy in eV. Again this is not required for total cross section although a dummy value must be given. Electron_Energy_Step - Real. Increment value used to determine the number of points of ejected electron energy. Again this is not required for total cross section although a dummy value must be given. a n g l e - Real. The polar angle of emission for the ejected electron, measured in radians. This value is only
required for the double differential cross section therefore if my_section = 0, 1 then a dummy value must be given. tol - Real. Global tolerance used throughout the code. Typical values may be given as 1.e-8, although one may run into problems if tol ~ e p m a c h (see below for definition of epmach). a c e - Real. Tolerance used for the nested integrals. This value is scaled by a factor of 0.1 as each new quadrature
routine is called. Therefore, for double differential cross sections, a typical value would be a c c = 0.0001. For single and total differential cross sections, reliable values are a c c = 0.001 and a c c = 0.01, respectively. The value of acc for total cross sections may appear to a little large for a tolerance, but anything lower may lead to computational difficulties, bearing in mind that we are dealing with four-fold quadrature. In the description of the input data above, it was stated that for certain cross sections dummy values were required. The program needs to read in a certain amount of data, the dummy values allow it to complete this task. As stated earlier, the program may evaluate cross sections for H, H2 and He-like systems. The user must provide an input file helium.dat which provides the RHF coefficients Za and ba (Eqs. (13) and ( 1 4 ) ) ; for He, the file contains the following values [4]: 5 1.41714
7.6837724e-1
2.37682 4.39628 6.52699 7.94252
2.2345803e-1 4.0818690e-2 -9.940800e-3 2.2963900e-3
B.S. Nesbittet al./ ComputerPhysicsCommunications114 (1998)385-400
391
The first value, 5, lets the program know how many Za and ba RHF coefficients need to be read in (this value is given by the upper limit in the summation of Eqs. (13) and (14)). For H and H2, the user must create a simple file hydrogen.dat which contains the following values:
5 1.0 0.0 0.0 0.0 0.0
1.0 0.0 0.0 0.0 0.0
For H, He and H2-1ike targets, the effective charges are defined in the code as ZT = 1.0_rap ZT = 2.0_rap* ( 2 4 . 5 8 / R y d b e r g ) * * 0 . 5 _ m p ZT = 1. 064_rap respectively. For H and H2-1ike targets, the final cross section is divided by two. Note that the factor of 2 arises from indistinguishability of the target's two electrons. Also, since our model is based on the independent electron model, we approximate the H2 target as an effective one-electron hydrogenic system, therefore Zr = 1.064. On completion, the program creates two output files; he_data.dat for detailed cross section information and total_output.dat which contains the raw data which is very useful for plotting graphs. Examples are given at the end of the paper.
3.3. Structure of program The code consists of five main programming units; the main program and four Fortran 90 modules organised as follows: Module
Precision ]
] Module Global-Data I Module
Complex_Functions I
I Module
Cross_Section I
Main Program
Ion_Atom ]
3.4. Module precision This module is the first stage in Ion-Atom-Wave and defines the compile time constants are given by mp = s e l e c t e d _ r e a l _ k i n d ( l O , 3 0 ) c_mp = m p
This defines the precision for real and complex variables using the extension real(kind--mp) :: x, y c o m p l e x ( k i n d = c mp) :: z
rnp and c_mp which
B.S. Nesbitt et al. / Computer Physics Communications 114 (1998) 385-400
392
also real constants are defined as 1.0rap instead of the usual 1.0 or 1.d0, for complex numbers the notation is z = cmplx(1.0_rap, t . 0rap, c_mp). These features greatly enhance the portability of the code. Module Precision contains one subroutine Machine_Accuracy which defines the following real constants:
epmach, uflow, tolerance. epmach - Real. Machine accuracy. uflow - Real. Smallest positive spacing between two real numbers. tolerance - Real. Machine dependent tolerance; i.e. t o l e r a n c e lower bound on tolerances required for numerical routines.
= 1.e+6 rap*epmach. This helps to set a
The constants mp, c_mp, epmach, uflow and tolerance are globally accessible.
3.5. Module Global_Data This module contains one subroutine namely Get_Data. Subroutine Get_Data is called once by the main program and reads the input data from the files info_totaLdat and hartree.dat associated with units getl = 8 and get2 = 9, respectively. The data discussed in Section 3.2 may then be accessed by other programming units by using the use Global_Data statement.
3.6. Module Complex_Functions This module contains five Fortran 90 real and complex programming functions. The module accesses the constants mp, c_mp, epmach and tolerance from module Precision. Sqmodn(x) - This function evaluates IN(a)l 2 using Eq. (6.1.31) of Ref. [6], where N(a) is given by Eq. (12).
function
Cgamma(z) - To evaluate the analytic continuations of the Gauss hypergeometric function in Eq. (8) we need to evaluate the complex gamma function F ( x + iy). This is best achieved using the Lanczos approximation [ 8 ].
function
Cdigam(z) - This function calculates the digamma function ~ ( z ) for complex argument z using an asymptotic expansion (Eq. (6.3.18) of Ref. [6])
function
1
g,(z) ~ lnz
l
2z
1
1
12z 2 + -120z ---q
252z 6 + . . .
(z --~ oo in larg(z)l < ~') •
(33)
If IRe(z)l < 15, the recurrence property of the digamma function (Eq. (6.3.6) of Ref. [ 6 ] ) , 1
~ b ( n + z ) - ( n - 1) + z
+
1
(n-2)
+z
+..-+
1
~
+
1
~
+~p(1 + z )
'
(34)
is used to retain accuracy when using the asymptotic expansion. If R e ( z ) < 0, the analytic continuation is used (Eq. (6.3.7) of Ref. [6]), ~O(1 - z) = ~/'(z) + ~ c o t ~ z .
(35
3.6.1. Gauss hypergeometric functions The most difficult special function to evaluate in the program is the Gauss hypergeometric function which is required for the evaluation of the scattering amplitude (Eq. ( 8 ) ) . Four programming units are used to evaluate 2F1 (a, b; c; z ); Feq, F_n, V1 and W1. The functions are required to evaluate the analytic continuations given below.
B.S. Nesbitt et al./ ComputerPhysics Communications114 (1998) 385--400
393
The Gauss hypergeometric series is defined as -~(a)n(b)nz
2Fl(a,b;c;z) =
n=O
n
(36)
(c)nnt
where the Pochhammer symbol is given by (a)n=a(a+l)(a+2)...(a÷n-1)
and
(a)0=l.
(37)
All the important properties of the Gauss hypergeometric function can be found in Ref. [6]. In particular, it is necessary to employ the following analytic conditions to enable numerical evaluation over the whole range of z ~ [-c~, +oo]: 2F1 (a, b; c; z ) = (1 - z ) - a 2 F l ( a , c - a ; c ; z / ( 1 - z ) )
([arg(1
2Ft(a,b;a+ b_c 2Fl(a,b;c;z) = F(c)F(c-a-b) F ( c - a) F ( c - b)
-z)l
<~r),
(38)
+ l ; l _ z ) + ( l _ z)C_a_b F ( C ) F ( a + b - c ) F ( a ) F ( b)
X2Fl(C-a,c-b;c-a-b+l;1-z) (I arg(1 - z)[ < cr) , F ( c ) F ( b - a) 2Fl(a,b;c;z) = ( - - z ) - a 2 F l ( a , 1 -- c + a; 1 -- b + a; 1 / z ) F ( b ) F ( c - a) l " ( c ) F ( a - b) 2Fl(b, 1 - c + b ; 1 -a+b;l/z) (larg(-z)l q-F(a)F(c-b)
(39)
<~),
(40)
and 2F1 (a, a; c; z ) =
r(c) F ( a ) F ( c - a)
oo
( _ z )_a ~-~ (a)n(1 -- c + a ) n z _ n
n--0
(n!)2
x [ l n ( - - Z ) + 2~b(n+ 1) -- ~ , ( a + n) -- ~ b ( c - a - n)]
(larg(--Z) < 7r),
(41)
as required in Eq. (8). The program unit V1 evaluates 2Fl(a,b; 1;z), unit W1 evaluates 2Ft(1 + ia, 1 + i b ; 2 ; z ) and the unit Feq evaluates 2 F l ( a - 1 + ib, a - 1 + i b ; b ; z ) , where a, b and z are real. More generally, the program unit F_n evaluates 2 F l ( a , b ; c ; z ) for complex a, b and c and real z provided that Re(c - a - b) > 0 and - 1 < z < 1. It is possible to remove the restrictions on Eqs. ( 3 8 ) - ( 4 1 ) provided it can be determined whether 1 - z and - z approach the negative real axis from above or below; i.e., whether we should write - 1 as e iTr or e -i~r. To determine this we note that the Gauss hypergeometric functions which appear in the CDW and CDW-EIS models occur as solutions of the standard Nordsieck integral [7], /
d~3rrexp(
~ ~r
+ iq- r ) I F I ( - i a l ; 1; i ( p l r + Pl" r) ) 1FI (ia2; 1; i(p2r + P2" r) )
271" ( O l 8 =
tr----~ ots--+fl8
)ia2 (OtS~ tal \Ys/
2FL(--ial,ia2; 1;7") ,
(42)
where
7"----and
flsV8 - ~ 8 8 8
(43)
394
B.S. Nesbitt et al. / Computer Physics Communications 114 (1998) 385~100
c~ = ½(q2 + h2) , Y8 = P l " q -- i a p j
/3 = P2" q -- i ap2, + c~,~,
68
=
Pl
"
P2 +/33-
(44) (45)
In general it is necessary to assume, for convergence, that Re(A) > 0. Clearly for a 4: 0, Eq. (43) has an imaginary part and the sign of this imaginary part determines whether the real axis is being approached from above or below as a --~ 0 +. For C D W theory as A ~ 0 + it is easy to show that the sign of I m ( r ) is given by the sign of - a [ ~ 3 p ( a + / 3 ) + y v ( ~ + Y - ce - / 3 ) ] .
(46)
Note that setting Pl = P, P2 = v and taking ,~ --~ 0 + in Eqs (44) and (45) yield the corresponding quantities in Eqs. ( 2 3 ) - ( 2 6 ) for C D W theory. Similarly, for CDW-EIS theory it can be shown that as A ~ 0 + then the sign of Im(~-) is given by the sign of - a [ 6 p f l + v y ( 6 + y) ] .
(47)
It is of great importance to have a test which is robust computationally and avoids division by zero. Eqs. (46) and (47) provide the definition of the variable TEST in function CS(x) in module Cross_Section. 3. 7. Module Cross_Section
Module Cross_Section uses modules Precision, Global_Data and Complex_Functions and contains subroutines and functions used to evaluate the cross sections. The following real variables are declared public and are therefore globally accessible: Z T - effective charge of the target. velocity - velocity v of the projectile in a.u. g l o b a l & - momentum of ejected electron k. The following real variables are declared private and can only be accessed by subroutines and functions contained with module Cross_Section. q_min - lower limit of the integral in Eq. (5). delta_epsilon
Ae given by Eq.
(6).
( given by Eq. (9).
xi nu
-
- u given by Eq. (10).
zeta - ~" given by Eq. (11). smnnu smnxi
-
N(v) given by Eq.
(12).
- N(().
snt
- sin 0.
cst
-
cosO.
k_dot_v - the dot product k • v. p_dot_v - the dot product p - v. AB - ] N ( u ) N ( ( ) N ( ( ) I .
B.S. Nesbitt et aL/ Computer Physics Communications 114 (1998) 385-400
395
csp - cos ~b.
The total cross sections are obtained by evaluating four fold nested integrals. Obviously we require extremely robust and accurate quadrature rules; additionally they must be efficient to evaluate cross sections within a reasonable time. The nested integrals are evaluated using a set of Gauss-Kronrod rules. These subroutines originate from the highly acclaimed set of automatic integrators known as QUADPACK [ 5 ]. For our purposes we have modified the QAGE, QK15 and QK61 subroutines and converted them to Fortran 90. In general the 61-point Gauss-Kronrod rule is used throughout the calculations except for peak problems (such as the ECC peak, defined below) where the 15-point rule is preferred. The subroutines and functions used to evaluate the cross sections and set up the quadrature are as follows: subroutine Totai_Differential(Modei_Energy, TDCS) - This subroutine is called by the main program and evaluates the outer integral over k where T D C S = total cross section if my_section = O. = single differential cross section if my_section = 1. = double differential cross section if my_section = 2. The Model_Energy can be either the projectile energy or the energy of the ejected electron, depending on the type of cross section to be evaluated. For single and double differential cross sections, Model_Energy :-- Electron_Energy then by default the projectile energy is given as Projectile_Energy_Start. Alternatively, for total cross sections Model_Energy -- Projectile_Energy. The function DDl(dummy_k) is then called. For single and double differential cross sections this function is called for specific values of k (the momentum of the ejected electron). For total cross sections subroutine QUAD4 is called to integrate the function DD1 over a range of k. In theory the outer integral should be evaluated from 0 to cx~, in practice, however, this is never necessary. Instead we truncate the range of integration to k_max, where k_max is related to the energy of the incoming ion Ek by k_max =
,-~ 3.5 x v,
(48)
which ensures that contributions to the total cross section from ionized electrons up to and including the Bethe ridge are included. Also, the point in momentum space where k = v has a simple pole in the DDCS, which is commonly referred to as the ECC peak (electron capture to the continuum'. The pole is due to go = Z J p , where p2 = k 2 + v2 _ 2 k - v ,
(49)
which only occurs in the forward direction (0 = 0) and where p is defined in function DD2. We avoid this point in the following manner; if I v - k J < • t h e n if v < k then v=v-~v else v=v+6v endif endif
where • = tolerance and ~v > •. This allows the singularity to be jumped over without losing vital information.
396
B.S. Nesbitt et al. /Computer Physics Communications 114 (1998) 385-400
function DDl(dummy_k) - If my_section = 0 or 1 (i.e., evaluation of total or single differential cross sections), DD1 calls subroutine QUAD3 to integrate the function DD2 over a range 0 to ¢r for a given value of k. If my_section = 2 (double differential cross section) then DD1 evaluates DD2 at a particular value angle which is read from module Global_Data. function DD2(polar) - This function calls subroutine QUAD2 to integrate function DD3 over ~b for specific values of k and 0. The integrand of the integral over ~b, the azimuthal angle of k with respect to the collision plane only depends on cos ~b, therefore the transformation f f ( c o s O) 0
f(-
d~b =
-,r
f ( - cosO) d~
(50)
0
can be made which halves the range of integration which in turn halves the time taken to evaluate this integral.
function DD3(polar) - Function DD3 calls QUAD1 to integrate CS over the range - 1 to +1 for a given value of ~b, O and k. function CS(x) - Function CS is the central integrand defined by Eq. (5). The range of the integration is changed from 7/ E [0, o~] to x E [ - 1 , +1] by the transformation r/2= 1 - x 1+ x
(51)
This function also evaluates the variable TEST which determines which of the four analytic continuations to use (Eqs. ( 3 8 ) - ( 4 1 ) ) as described in Section 3.6. For the CDW model TEST is defined by Eq. (46), for the CDW-EIS model it is defined by Eq. (47). Module Cross_Section also contains the functions V1, W1 (which where discussed in Section 3.6) and the QUADPACK routines QUAD1-QUAD4.
3.8. Main program Ion_Atom The main program uses modules Precision, Global.Data and Cross_Section. The first executable statement is to call subroutine Machine_Accuracy which determines the machine dependent constants. The next statement calls subroutine Global_Data which reads the input data from file. This data can then be accessed by the other programming units. The output files he_data.dat and total_info.dat, associated with units 6 and 4, respectively, are then initialised.
3.8.1. Parallelization of code The main program incorporates a number of High Performance Fortran (HPF) directives allowing large sections of the code to be run concurrently in a SIMD (single instruction multiple data) fashion. This is achieved as follows. Consider the total cross section, the user initially supplies starting and ending values for the kinetic energy of the projectile, Projectile_Energy_Start and Projectile_Energy_Finish, respectively. An increment value is also supplied, Projectile_Energy_Increment, enabling the program to create an allocatable array Energy of energy points which contains N elements where the number of elements is determined from N = Etinish - - E s t a r t + 1. Eincrement
The program is then parallelised as follows:
(52)
B.S. Nesbitt et al. / Computer Physics Cormnunications 114 (1998) 385--400
397
!HPF£ PROCESSORS procs(number_of_processors) !HPF£ DISTRIBUTE (BLOCK) ONTO proes :: Energy, Cross !HPF£ INDEPENDENT do loop = 1,N call Total_Differential(Energy(loop), Differential_Cross_Section) Cross(loop) = Differential_Cross_Section end do where the do loop evaluates the cross section N times for N different values of energy. For single and double differential cross sections, projectile energy is replaced with ejected electron energy. Note that each iteration of the do loop is completely independent of all others enabling the HPF directive to allocate the N calls to the subroutine Total_Differential to an array of processors.
Acknowledgements Brian Nesbitt acknowledges financial support from the Department of Education, Northern Ireland (DENI). Francesca O'Rourke acknowledges current financial support by a Royal Society University Research Fellowship and Derrick Crothers acknowledges EPSRC support under grant GR/L21891.
References [1] [2] [3] [4] [5]
D.S.E Crothers, J. E McCann, J. Phys. B 16 (1983) 3229. P.D. Fainstein, V.H. Ponce, R.D. Rivarola, J. Phys. B 24 (1991) 3091. S.EC. O'Rourke, D.S.E Crothers, J. Phys. B 30 (1997) 2443. E. Clementi, C. Roetti, At. Dam Nucl. Tables 14 (1974) 177. R. Piessens, E. de Doncker-Kapenga, C.W. Uberhuber, D.K. gahane, QUADPACK, A Subroutine Package for Automatic Integration, Springer Series in Comp. Math., Vol. i (Springer, New York, 1983); QUADPACK is freely available at ~ n ~ . n e t l i b . o r g or at the UK mirror, ~ r m ~ . h e n s a . a c . u k / c g i - b i n / c f t p / u n i x . hensa, ae. uk/mirrors/net lib~
[6] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions (Dover, New York, 1970). [71 A. Nordsieck, Phys. Rev. 93 (1954) 785. [8] C. Lanczos, J.S.I.A.M. Num. Analy. B 1 (1964) 86.
B.S. Nesbitt et al./Computer Physics Communications 114 (1998)385-400
398
TEST RUN OUTPUT Input Data Test data for file info_total.dat generated by the program Ion-Atom-Driver: 1 2 2 I. 000000000000000 75. 00000000000000 1. 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 i. 000000000000000 O. 2000000000000000 90. 00000000000000 2. 000000000000000 O. i000000000000000 1. O000000000000000E-O i0 1. O000000000000000E-O04
Test data for file helium.dat: 5 1.41714D0 7. 6837724D-01 2. 37682D0 2. 2 3 4 5 8 0 3 D - 0 1 4 . 3 9 6 2 8 D 0 4. 0 8 1 8 6 9 0 D - 0 2 6.52699D0 -9.940800D-03 ?.94252D0 2.2963900D-03
Test data for file hydrogen.dat: 5 1.0 0.0 0.0 0.0 0.0
1.0 0.0 0.0 0.0 0.0
Output Data Ion-Atom-Wave
Ion-Atom-Wave calculates cross sections for the single ionization of atomic and molecular hydrogen-like targets and helium-like targets by both light and heavy ion impact.
Input .
.
.
.
.
.
data .
.
.
.
.
used for .
.
.
.
.
.
.
.
.
.
this .
.
.
.
.
run .
.
.
:
.
The target is a helium atom The CDW-EIS model was used. double differential cross section evaluated for the following values :
B.S. Nesbitt et at~Computer Physics Communications 114 (1998) 385-400 Program started with the following electron energy (eV) : Program ended with the following electron energy (eV) : The electron energy was incremented by Energy of the incoming ion (keY) Angle of electron emission (in radians)
0.20000000E+O0 0.90000000E+02 0.20000000E+01 0.75000000E+02 0.10000000E+O0
Numerical details
Tolerance used throughout the code : 0.22204460E-09 Tolerance used fo r the quadrature : 0.10000000E-03
PROGRAM OUTPUT
The first column contains the electron energy (eV) and the second column contains the cross sections
0.200000000000000 2.200000000000000 4.200000000000000 6.200000000000000 8.199999999999999 10.20000000000000 12.20000000000000 14.20000000000000 16.20000000000000 18.20000000000000 20.20000000000000 22.20000000000000 24.20000000000000 26.20000000000000 28.20000000000000 30.20000000000000 32.20000000000000 34.20000000000000 36.20000000000000 38.20000000000000
40.20000000000000 42.20000000000000 44.20000000000000 46.20000000000000 48.20000000000000 50.20000000000000 52.20000000000000 54.20000000000000 56.20000000000000
2.3394624948705482E-OO2*le-20sqr.metres/eV/sr 3.0481267769151714E-OO2*le-20sqr.metres/eV/sr 3.2186809295606460E-OO2*le-20sqr.metres/eV/sr 8.2579413664121801E-OO2*le-20sqr.metres/eV/sr 3.2480520852528200E-002*le-20sqr.metres/eV/sr 3.2224803008163452E-OO2*le-20sqr.metres/eV/sr 3.1975245277550073E-O02*le-20sqr.metres/eV/sr
3.1824782971607055E-OO2*le-20sqr.metres/eV/sr 3.1830796483122169E-OO2*1e-20sqr.metres/eV/sr 3.2038672844994195E-OO2*le-20sqr.metres/eV/sr 3.2481768424916561E-002*le-20sqr.metres/eV/sr 3.3207862743912936E-O02*le-20sqr.metres/eV/sr
3.4263373599848164E-OO2*1e-20sqr.metres/eV/sr 3.5705354273153703E-OO2*le-20sqr.metres/eV/sr 3.7597027469323460E-OO2*le-20sqr.metres/eV/sr
3.9994014059866678E-OO2*le-20sqr.metres/eV/sr 4.2902732491859850E-OO2*le-20sqr.metres/eV/sr
4.6179049743754010E-OO2*le-20sqr.metres/eV/sr 4.9337403545042326E-OO2*le-20sqr.metres/eV/sr 5.1366182908066736E-OO2*le-20sqr.metres/eV/sr 5.0990543575597416E-OO2*le-20sqr.metres/eV/sr
4.7703676460025074E-OO2*1e-20sqr.metres/eV/sr 4.2414824071991944E-OO2*1e-20sqr.metres/eV/sr 3.6591704956531815E-002*le-20sqr.metres/eV/sr 3.1219229408178922E-O02*le-20sqr.metres/eV/sr
2.6653835839971593E-OO2*1e-20sqr.metres/eV/sr 2.2902354699969769E-OO2*le-20sqr.metres/eV/sr 1.9848356642432175E-OO2*le-20sqr.metres/eV/sr
1.7356716559748971E-OO2*le-20sqr.metres/eV/sr
58.20000000000000
1.5808863641622825E-OO2*le-20sqr.metres/eV/sr
60.20000000000000 62.20000000000000 64.20000000000000 66.20000000000000 68.20000000000000
1.3609851200881467E-002*le-20sqr.metres/eV/sr
70.20000000000000 72.20000000000000 74.20000000000000
1.2186184246292624E-OO2*le-20sqr.metres/eV/sr 1.0981554216192674E-OO2*le-20sqr.metres/eV/sr 9.9527930646778721E-OO3*le-20sqr.metres/eV/sr 9.0666130655147580E-OO3*le-20sqr.metres/eV/sr 8.2971324892396683E-OO3*ie-20sqr.metres/eV/sr 7.6240400385084737E-OO3*1e-20sqr.metres/eV/sr 7.0312447236584117E-OO3*le-20sqr.metres/eV/sr
399
B.S. Nesbitt et aL / Computer Physics Communications 114 (1998) 385-400
400
76.20000000000000 78.20000000000000 80.20000000000000 82.20000000000000 84.20000000000000 86.20000000000000 88.20000000000000
6.5058772606735949E-OO3*le-20sqr.metres/eV/sr 6.0375535399216493E-003*le-20sqr.metres/eV/sr
5.6178225246002788E-OO3*le-20sqr.metres/eV/sr 5.2397512887970713E-OO3*ie-20sqr.metres/eV/sr 4.8976089630747618E-003*le-20sqr.metres/eV/sr 4.5866258592475269E-003*le-20sqr.metres/eV/sr 4.3028056253776375E-003*le-20sqr.metres/eV/sr