Calculation of the auxiliary functions Fm(z)

Computer Physics Communications 31(1984) 47—52 North-Holland, Amsterdam

47

CALCULATION OF THE AUXILIARY FUNCTIONS Fm(Z) L.F. ERREA, L. MENDEZ and A. RIERA Departamento de Qulmica FIsica y Quimica Cudntica, Centro Coordinado Blanco, Madrid- 34, Spain

CSIC

-

UAM, UniversidadAutdnoma de Madri4 Canto

Received 18 April 1983; in final form 11 July 1983

PROGRAM SUMMARY Title of the program: FRANKC Catalogue number: ACFL

Keywords: atomic and molecular collisions, translation factors, molecular integrals, complex Gaussian orbitals, diamagnetic susceptibilities, electron—molecule resonances, complex selfconsistent-field method, band~structurecalculations in solids

Program obtainablefrom: CPC Program Library, Queen’s University of Belfast, N. Ireland (see application form in this issue)

Nature of the physical problem

Computer: IBM 370; Installation: Centro Coordinado UAM/IBM, Canto Blanco, Madrid-34, Spain

Calculation of a set of auxiliary functions F,,,(z) which appear in the analytical calculation [1~of molecular integrals involving Gaussian type orbitals and ge~eralisedplane waves.

Operating system: OS/VS1 Programming language used: FORTRAN IV

Method ofsolution The complex z plane is dividcd into eight domains and optimized algorithms are used for each domain. For z in the lower half plane, the principle of reflection is used.

High speed storage required: 3 Kwords Restrictions on the complexity of the problem No. of bits in a word: 32

The maximum value of m in F,,, ( z) is set to be 8. The relative error is less than 0.5 )< 106 provided that Re(z)> —30.

Peripherals used: card reader, line printer No. of cards in combinedprogram and test deck: 450

Typical running time For the test run in the IBM 370: 0 mm 0.86 s.

Card punching code: EBCDIC Reference to other CPC library programs: The program FZRI, available from the CPC program library, also provides the calculation of F,,, (z), in a small region of the complex plane, with a higher precision and larger values of m than the present program. The program FZRI uses Gauss—Legendre integration, and is, therefore, much slower than the present version.

eierence [11 L.F. Errea, L. Méndez and A. Riera, J. Phys. B 12 (1979) 69.

OO1O-4655/84/$03.OO © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

L.F. Errea et al.

48

/ Auxiliary functions

F,,,(z)

LONG WRITE-UP 1. Introduction The molecular approach to atomic collisions presents two basic difficulties that cannot be eliminated without introducing substantial modifications in the theory. The first difficulty is related to the existence of constant radial couplings, and of slowly decreasing rotational couplings, between the molecular states at large internuclear separations. The second difficulty is that calculated cross sections can depend upon the choice of the origin of electronic coordinates chosen to perform the calculation. These difficulties of the molecular approach were first pointed out by Bates and McCarroll [1] and their importance in calculations has become increasingly conspicuous. They arise because of the absence in the molecular expansion of Electronic Translation Factors (ETF or TF) which account for the translational motion of the electrons with respect to the origin of electronic coordinates [2]. Analogously, translation factors should be introduced using an atomic (e.g. pseudostates) expansion. For detailed studies of the problem and its possible solutions see e.g. ref. [3,4] and references therein. When a TF of the form exp(iqu( R) (1)

)

is introduced in the molecular expansion, a well known difficulty is the evaluation of the molecular integrals; here r is an electronic coordinate, R the internuclear distance, q a constant and lim R ~u ( R) = v, the nuclear velocity. However, we have shown [5] that all molecular integrals involving TF of the form (1) and atomic basis sets of Gaussian-type orbitals (GTOs) can be evaluated analytically in terms of a set of auxiliary functions: Fm(z)=flu2mexp(_zu2)du,

tegrals. In this work we present a program which we have been using since 1978 to evaluate (2) with a precision that is sufficient in molecular calculations, and whose computational effort has been (approximately) optimized. It is worthwhile noting that the auxiliary functions Fm( z) of eq. (2) are needed in several other areas of research, and the need of efficient algorithms to perform their calculation has been pointed out. Examples of these areas are: 1) The calculation of diamagnetic susceptibilities of molecules, when basis sets of GTOs are employed [6]. In this case, u is the vector potential. 2) The use of discrete basis functions to describe scattering events [7]. In this case not only the present program but the whole formalism of ref. [5] apply directly. 3) Complex coordinate calculations of electron—atom scattering resonances, when GTOs are used to construct the atomic wavefunctions [81. The problem of evaluating the integrals, even for the electron—molecule case, is solved in ref. [5]. 4) The complex self-consistent-field method of McCurdy et a!. [9]. 5) Use of basis sets of the form sin(p.r) 2), cos(p.r) exp(—ar2) [10]. exp(—ar 6) Band-structure calculations using mixed basis sets of plane waves and Gaussians [11]. 2. Theory The entire function Fm(Z) is an extension to the complex plane of the function Fm (x), with z

=

retO

x+

=

iy,

(3)

and Fm(X) has been defined by Shavitt [12]. It is closely related to the complex error function and to the incomplete gamma function [13]:

(2) Fm(Z)

where z is a complex number and m a positive or zero integer. We have presented in ref. [5] closed form formulae in terms of finite sums, which are numerically stable and can be coded directly for the automatic computation of the molecular in-

=

‘TT1~2

—

I

1)m_~.___1Z_1/2 erf(z1~’2)]

—m-I/2

dzm

(

~y m + ~,

—

z).

(4)

By the principle of reflection [14]: Fm(Z*)

=

F,~(z),

(5)

L.F. Errea et al.

/

where * stands for complex conjugation, and only the upper half z plane needs to be considered. When a value of z in the lower half plane is given as input to the program, it calculates Fm(Z*) and then uses eq. (5). An important feature of the calculation of molecular integrals is that when .F,(z) is needed, so are Fm(Z), m = 0, 1 hence, a matrix array of those values is given as output. The parameters involved in the different algorithms we employ were chosen so as to optimize the computer time for each value of max, the maximum value taken by m, and each domain into which we have divided the z plane. We have taken an upper bound of max = 8, which allows the evaluation of all molecular integrals involving up to 3d GTOs for systems with more than one electron, and up to 5g GTOs for one-electron systems. Extending the program to max> 8, besides requiring a careful error analysis such as the one carried out in ref. [5], complicates the expressions which define the parameters and domains of the z plane. There is, however, no intrinsic difficulty. We guarantee relative errors in the real and imaginary parts of Fm that are less than 0.5 X 106; we have tested this by comparison of the results of our program with those obtained by direct numerical integration using the extended Simpson’s rule. However, the functions Fm(Z) are ~(1) for x ~ 0, but CV (e X) for x <0. Hence, as x becomes large and negative, the functions Fm increase exponentially, and oscillate for y # 0. Moreover, the region of x large and negative never appeared in ,

-

in Z

Auxiliary functions Fm(z)

49

our calculations, since it would involve unphysically large nuclear velocities [5]. Therefore, we have limited our tests to the region x> 30, which is amply sufficient for all practical purposes. For x < —30, the program still calculates Fm(Z) but we do not guarantee a relative error less than 0.5 x 10—6. It should be noticed that the asymptotic region x —s + no is essential in practical calcu2) [5]. lations, since XR = P( R In many steps of the program, double precision is used to avoid errors due to numerical cancellations. The input and output of the program are, of course, in single precision format. To calculate Fm(Z), we have divided the region x> 30, y ~ 0 into eight domains (see fig. 1). —

—

2.1. Region I; r < 3, y # 0

This corresponds to the small r region. Integration by parts of eq. (2) or use of the recursion formula of the incomplete gamma function [13] shows that the Fm functions obey the recurrence relation: (2m + 1)_I[2zFm+s(z)+e_:]. (6) When r is small, this formula can be used downwards in the same way as for the class of real Fm(X) functions [12]: for a sufficiently high order m = M, F,,,(z) is assumed to be zero, and eq. (6) is used for successively lower orders. The value of M is chosen to be a simple step function of r and 0. It is assumed that one requires Fm(Z), m = 0, 1,. ,max and max ‘s~8. Fm(Z)

. .

2.2. Region II; 3
~:

by the+ (truncated) formula e~2 by Abramowitz erf(x iy) erf(x) x[(1 —cos2xy)+isin 2xy] 2 e~2 25 given ~ + ~ ~ 2+4x2 +-~--——

—

1n Fig. 1.

X [f~(x,y) + ig~(x,y)],

(7)

50

L.F. Errea et al.

/ Auxiliary functions F,,, (z)

where

2.5. Region V; 12.25
2x(1 cosh ny cos 2xy) +n sinh fly sin 2xy, g~(x,y) = 2xcosh ny sin 2xy +n sinh ny cos 2xy. f,,(x, y)

=

<

25, 0

<

0 < 1.22

Same algorithm as for region III.

—

(8)

Then F0 is calculated through: 2z~”2 erf(z”2) (9) .F0(z) = ~“ and the recursion formula (6) is used upwards for

2.6. Region VI; r> 25, 0<0< 1.22 Same algorithm as for regions II and IV, except that F0 is not calculated using eqs. (7)—(9), but is set to 1, region thus rendering this equal asymptotic very fast.the calculation in

F 1, F’,. ~Fmax~ up to the required value of max(when max # 0): . .

Fm+i(z)=~[(2m+1)Fm(z)_e_z1.

(10)

This is the region where our program is slowest, 2.3. Region III; r> 16, 1.22 <8

.<

I~

This corresponds to the large r region. The following expression is used [5]: m~’2 ~m(Z), (11) Fm(Z) = F(m + fl/2z where F is the gamma function [13] and ~m the complementary function: —

~m(Z)

(12)

=f~u2m exp(—zu2)du.

2.7. Region FR,’ x ~ 0, y = 0 This region corresponds to the special (and important) case when no TF are used, for example when acutoff factorisintroducedin u(R)ofeq. (1). The program sets Im Fm = 0 and employs algorithms that are very similar to those used for complex z; however, the calculation of Fm(X) is much faster because it operates only with real numbers. The program calls subroutine Frank, which evaluates Fm(x); m = 0, 1,... ,max. Then: i) If x < io~, the program sets Fm(X) = (2m + 1)~. Otherwise, it uses different algorithms, depending on whether max = 0 or not. ii) When max = 0, and x < 16, a Chebyshev expansion of F~is employed [15] .E~(x)=

c

~‘

5T5(~x 1),

(14)

—

I

We evaluate ~m(Z)r~.~~L1+

+

e

~m —

z

using the asymptotic series [13]: ~

(2 m 1) 2z

coefficients calculated from those given by Clenshawhave [16]been for the error function:

—

_________________________________ (2m — 1)(2m — 3)

(2z)2

where 77, is a Chebyshev polynomial [13] and the 1/2

(13)

+

1”2

=

erf x

~~j—~‘a S 1/2

=

~

2/4) 237~3(x”

~a~T

5(~x

Our program truncates this series when either the new term to be added is greater than the previous one or it is less than 10 12 This algorithm is very fast. 2.4. Region IV; 3
<

12.25, 0

<

0

Same algorithm as for region II.

<

1.22

—

1).

(15)

Using (9), yields expansion Theevaluation values of andthis ç are given as data, (14). and the ~I/2/2

of the sum (14) is carried out as indicated in ref. [16]. iii) When max = 0 and x ~ 16, the program sets: 2x~”2. (16) F0(x)=~”

L.F. Errea et al. / Auxiliary functions F,,,(z)

51

iv) When max > 0 and x < 16.5, the program uses the same algorithm as for region I, taking the value of m = M for which FM is set to be zero to be a simple stepfunction of x and max. v) When max > 0 and x ~ 16.5, the same method as for region VI is employed.

4. Test run output

2.8. Region D; x <0, y = 0

Acknowledgements

In practice, this region is only used for small values of lxi. The program calls subroutine DAWSON, which makes use of the relationship:

We wish to thank the staff of the Centro de Cálculo CC-UAM and of the Centro Coordinado UAM-IBM for technical assistance.

The output for the test run provides values of Fm ( z) for one value of z in each of the 8 regions of fig. 1.

(17) where F(x) is Dawson’s integral [13]. We first calculate x”2 F(x1”2), using the rational approximations proposed by Cody et al. [17]; the coefficients of these expansions have been rounded off and the number of terms reduced to achieve a higher computational speed, and because a lower precision is required than in ref. [17]. Then expression (10) is used in the form: eH~Fm+i(Hxi)

=

[i

~

—

(2m + 1)

—lxl)1

e~1Fm(

(18)

which is more stable since e’~Fm(—ixi)= P(lxL’) for x + no. Multiplication by e1”1 then yields Fm(~lXl)=Fm(X). —

3. Organization of the program The submitted program is a subroutine, FRANKC. A main program is provided for the test run. Arguments of calling sequence: ZS: value of z in Fm(Z), FA: array of complex giving the output Fm ( z), M: integer. Number of values of Fm required. Called max + 1 in the preceding section.

References [1] D.R. Bates and R. McCarroll, Proc. Roy. Soc. A245 (1958) 175.Commun. At. Mol. Phys. 12 (1982) 82. A. Riera, J.B. Delos, Rev. Mod. Phys. 53 (1981) 287. A. Maciás and A. Riera, Phys. Rep. 90 (1982) 299. L.F. Errea, L. Méndez and A. Riera, J. Phys. B 12 (1979) 69. [6) G.G. Hall, Intern. J. Quantum Chem. 7 (1973) 15. [7] T.N. Rescigno, C.W. McCurdy and V. McKoy, Phys. Rev. All (1975) 825. [2) [3] [4] [5]

[8] N. 24 (1981) Moiseyev, 1254.P.R. Certain and F. Weinhold, Phys. Rev. A [9] C.W. McCurdy, T.N. Rescigno, E.A. Davidson and J.G. Laurerdale, J. Chem. Phys. 73(1980)3268. See also C.W. McCurdy and R.C. Mowry, Phys. Rev. A 25 2529. N.C. Handy and S.F. Boys, Mol. Phys. 26 [10] (1982) D.J. Allison, (1973) 715. [11] S.G. Louie, Phys. Rev. B 22 (1980) 1933, and references therein. [12] I. Shavitt, Meth. Comput. Phys. 2 (1963) 1. [13] Handbook of Mathematical Functions, eds. M. Abramowitz and A. Stegun (Dover, New York, 1965). [14] E.C. Titchmarch, The Theory of Functions (Clarendon Press, Oxford, 1932). [15] A. Riera and J.W. Linnett, Theoret. Chim. Acta 18 (1970) 265.

[16] C.W. Clenshaw, National Physical Laboratory Mathematical Tables, vol. 5 (1963). [17] W.J. Cody, K.A. Paciorek and H.C. Thacher Jr.. Math. Comput. 24 (1970) 171.

52

L.F Errea et al.

/

Auxiliary functions F,,,(z)

TEST RUN OUTPUT 2.

0.100000E+01 M •

74 M 74 7414 71 74

2’

o

•

1

• •

2 3

•• •

,~ ~

7 • 3 •

• • •

0 1 2

3 4 5

•

-‘

• • • •

F(~) •

4(0)

0.1894720.00 0.1)02191.00 0.6473220—01 49o232E—01 0. 0.3936491—01 0.3256701—01 0.2774600—01 0.2~i5531—01

•

4(2)

• •

•

6(2)

7

4(2) • 6(2) •

0.1619320.00 0.2696700—U 0.1348351—03 0.1123630—04 o.ooiosoo—os 3.1966341—05 O.3o0,900—)7 3.721074E—08 0.1952680—08 0.OC000CE+01

4 4 74 :4 1 Si

0.6967121.00

• •

• •

H • Si •

0 3

74(2) 74(2) PU)

3

,(2)

.

‘(2) FU)

6 7

74(2)

3

PU)

1

5

• • •

4(2) •

0.4000021.31 Si Si

• •

14 71 • M

74 74 H

• • •

71 •

0 1 2 3 4 5

6 7 8

4(2) ~(Z) 4(2) 4(2) 4(Z)

0.1577081+00

0.7~78101—01 0.4817661—01 0.3431)90—01 0.2635091—01 0.2126651—01 0.1775981—01 0.1521221—01

0.3775230*03 0.2414241—01 0.1372361—02 —0.1039411—02 —3.2204060—02 —3.21+1262—02

6(2)

—0.1799490—02 —0.1575460—32 —0.1334281—02

PU)

6(2)

3.1500000~02

0.8000301.01

71 • H •

0.2)85241.00 0.4519991—02 0.1891290—03

4(2)

1

PU)

74 • 2 H • .3 44 • 4

~(Z) PU)

H

•

5

H • 6 7 74 • 8

PU) 4(2) 5(2) P02) 6(2)

0.0

3.0 0.0 0.0 0.0 0.0 3.0 0.0

• • • • • •

0

•

0.1403000.02

FU) 4(0)

2 3 4

4(0) • ~() P(~) •• FU) •

o

4(2)

7 B

4(2) 4(0)

0 1

14 -I .1 .7 74

2 3 4

•

6 7

4 ‘1 4

•

4(0) ~(Z)

•

4(2)

5

• • •

l~ —J.120300(+0Z

—0.1786731.00 —0.9226491—01 —0.6039851—01 —0.4439291—01 —0.3491531—01 —0.2869781—01 —0.2432451—01 —0.2108951—01 —0.1860371—01

‘4 4 1 14 ‘2 4

74

7 8

6(Z) 6(2)

•

74 • •

—0.1260851.00 —0.3361931—00 —0.1331340—01 —0.6618241—02 —0.3114821—02 —0.2437640—02 —0.1679051—02 —0.1223780—02 —0.9317830—03

6(2)

74 II 71 •

5

•

—0.0 3)0046.31

—0.6664.730—08 —3.0141230—06

0.2894958—06 0.2157680—06

0.1144141,02

0.1021061+02 0.9182701.01 0.8319111.01

3.3~0753E.04 Q.13~9~.04 0.3324330.04 0.3226288.04 3.3125901.04 0.3324601+04 0.2323510+04 0.2203361.04 3.2725431.04

—0.3544961.33 —0.4184691.03 —0.5936011.03 —0.6908521.03 —0.7797080+03 —0.9516571.03 —0.9254401+03 —0.9838441.03 —c..~~~69E’04

0.0 •

0.1.s2551.01 0.6079151.00

0.0 0.0

•

0.41 7+186.01

0.0

0.21 55962+20 ).25~554E+00 2.2136+61.00

0.0 0.0

14 •

0

‘(2)

74 H

•

14 14 14 74 14 14

• •

1 2 3 4 5

6(0) 4(1) 5(2) ~(Z) 4(1)

6

6(3)

0.19’076E+00

0.0 0.0

7

4(5)

0.1426480.00

0.0

4(2)

2.1392330.00

0.0

•

• • •

‘8

Z~ —0.1S00030+02

—0.9997200—05 —0.4023551—05 —5.1436110—05

0.1Z1273~.02 0.1174320.02 2.1133480+02 0.1195050.02 O.1C40230+02

0.2136210.02 0.1946831+02 0.1691531+02 0.1473190.02 0.1292751.02

u.2J1030(.02

74(2) 4(0) ‘U) 4(0) PU) 6(2)

3

0.1003322.02 0.1130.82.02 2~233361.02 0.1:3~s2E.02

4(2) 4(2) 4(0) 4(Z)

0.399947~—01 —0.1951851—02 —0.8960011—04 —0.6225961—05 —0.42765SE—06 —C.14222’.e—07 0.4802991—08 0.1948955—08 0.5248331—09

0.5003002.01

‘(2)

F(Z)

14

.1443340+00 1475451—32 0.2001261—06 ).~3e.736E—05 —0.49784’.t—06 —0.8321001—07 —0.1036741—07 —0.1503230—08 —G.o2~596t—10 ~.

2~ —0.600000E•01

—0.5213101—01 0.423902E02 —0.5247600—03 —0.8117301—04 —0.135184104 —0.19235.1—05 0.7(71411—07

—0.11769.0—34

M 4 • 4 • 74 ‘7 M • A

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

0.3000001+01 •

F(Z)

0

5.30c)000*62 0.0

0.0

4(0) -(Z) ~(l) 4(2) F(Z)

0.1000001.01

14

0.

4(2) ‘(0) ~(Z) 4(Z) 4(2) F(0)

6 8

2.

o.i—so’o.oo

0.3000000.02 14 4 74 Si ‘4 1 1 4

2’

0.0

~(i)

14 74 74 71

‘4 .1 74 ‘4 14

0 I

4(2) PU)

~

4(2) PCI) 4(0)

. — • • •

2 3

•

5

• •

I

7 • 8

PCi) 4(j)

5(j) 4(3)

0.3 •-

• • • • • • •

3.1030451.06 0.1051990.06

0.0 0.0

0.3044730,05

0.0

0.3255935.05 3.0737001•05 0.UtS628+05 0.7852320+05 0.1469711+05 0.7151662+05

0.0 C.0 0.0 0.0 0.0 0.0