- Email: [email protected]
Two-center Coulomb functions
ComputerPhy~cs Commun~cat~'~s ELSEVIER
Computer Physics Communications 103 (1997) 209-216
Two-center Coulomb functions Miyabi Hiyama a,, Hiroki Nakamura a.b a Department of Functional Molecular Science. School of Mathematical and Physical Science, The Graduate University for Advanced Studies. Myodaifi, Ok~aki 444. Japan b Division of Theoretical Studies, Institute for Molecular Science, Myodaiji. Okazaki 444. Japan
Received 19 February 1997; revised 14 April 1997
Abstract Two-center Coulomb functions are crucial to describe the ionization continua of diatomic molecules, since the latter cannot be well represented by finear combinations of one-center Coulomb functions. A code to numerically evaluate the two-center angular as well as radial regular Coulomb functions is presented here. ~ ) 1997 Elsevier Science B.V. Keywords: Schr6dinger equation in spheroidal coordinates; Two-center Coulomb function; Numerical solution
PROGRAM SUMMARY
?[tle of program: TCOULOM Catalogue identifier: ADFY Program obtainable from: CPC P r o g ~ Lib,"ary,Queen's University of Belfast, N. Ireland Licensing provisions: noae Computers: HP 735 Operating systems under which the program has been tested: UNIX Programming language used: FORTRAN77 No. of bytes in distributed program, including test data. etc.: 39946
Distribution format: ASCII Keywords: Schr0dinger equation in spheroidal coordinates, twocenter Coulomb function, numerical solution Nature of physical problem Subroutine TCOULOM is a FORTRAN77 subroutine to c a k e ] ~ the regular two-center Coulomb functions for the positive onergy, 6. This program requires six constants: R, the hUemuciear ~ ; Z! and Z2, two positive charges: m and q, ~ m numbers; and e, the energy of an outgoing electron. Atomic uni~ are employed. Method of solution There are flu~e main subroutines in this pmgrara. The first one is to estimate the separation constants. The chain equation is solved there (see Ref. J I l L The other two are to o w n wavefunctions of the angular part, ~.~(v/;k,R), and of the radial paff, II,r~( ~ ; k, R ) . The spheroidal angular variable vI ranges between - I and Io and the radial vatiabie ~ ranges from l to influX. ~'~(~,;k, R) is o ~ by solving the S c ~ g e r eqea~a numerically, where we put 7/= tanhq. The hategral range o f q is
I Present address: Department of Pure and Applied Sciences, Graduate School of Arts and Sciences, The University of Tokyo, Kom~,ba, Meguro-ku, Tokyo 153, Japan.
0010-4655/97/$17.00 ~) 1997 Elsevier Science B.V. All fights reserved. PII S0010-4655 ( 97 ) 00045-3
M. Hiyama. H. Nakomura/Computer Physics Communications 103 (1997) 209-216
210
(--co, oo), but in practice the range between -4.5 and 4.5 was found to be enough. The Numerov method [21 is used to solve the differential equation. The boundary conditions are provided by pow~ series expansions .at both ends. llmq(~;k,R) is obtained by solving the corresponding SchrOW~ger equation with use of both the forward and backward Nume~,,v procedure (see Ref. 121 ).
Restrictions on the complexity of the problem This program may be used to evaluate the two-center Coulomb functions for the positive energy e. The two positive charges are
assumed to satisfy the relation Z2 _> Z:. The wavefunction for the case of Z2 = Zi also can be calculated using this program. m and q are positive integers. They correspond to magnetic and azimuthal quantum nurabers. This program produces only regular solutions both for the angalar and radial parts.
References [ I ] L.I. Pouomamv and L.N. Somov, J. Comput. Phys. 20 (1976) 183. 12] H. Takagi and H. Nakamum, Report of The Institute for Plasma Physics, Nagoya University, IIPJ-AM-16 (1980).
LONG WRITE-UP
1. Introduction Two-center Coulomb functions represent bound as well as continuum states of a diatomic molecule with one electron. Bates et al. reported the exact wave functions of the bound states of H~" and Hell 2+ [ 1]. Helfrich et al. extended the calculations to other one-electron diatomic molecules such a LiH 3+, Hell 2+, Lille 4+, and H~[2,3]. These can he used as basis functions to evaluate bound states of diatomic molecules which have more than one electron [4]. The two-center Coulomb functions at positive energy (e > 0) is also useful to represent the ionization continuum of a diatomic molecule or a long range part of the electron scattering by a diatomic molecular ion [4,5]. In the dynamic processes involving superexcited states of diatomic molecules such as photo-ionization, autoionization, Penning ionization, and dissociative recombination, the ionization continuum naturally plays an important role. The multi-channel quantum defect theory (MQDT) presents a powerful tool to analyze these processes [6-8]. The first kind of superexcited state, whose excitation energy is higher than the lowest ionization potential, is completely embedded in the ionization continuum and autoionizes even with the internuclear distance kept fixed. The autoionization rate is represented by the electronic coupling defined by
V(R, 8) = (~¢o.tlHell~a),
(1)
where Hel is the electronic Hamiltonian of the system, and t/rcont and grd are the wave functions of the ionization continuum and of the superexcited bound state, respectively. This is one of the most important parameters to characterize the first kind of superexcited state, and is the most difficult quantity to oe properly evaluated. The two-center Coulomb functions can be used to construct the wave function l/rcont and the conventional CI code enables us to evaluate the electronic coupling. In order to facilitate this it is necessary not only to numerically compute the regular two-center Coulomb functions, but also to express them in terms of Gaussian functions. The first part, i.e. the numerical solution of the regular two-center Coulomb functions, is presented here, and the second part is discussed in an accompanying paper. Takagi and Nakamura reported the numerical solutions of the two-center Coulomb functions for the homo-nuclear case (Zi = Za) [9]. Analytical expansions of the functions for the hetero-nuclear case (Zi @ Z2) are discussed by Ponomarev and Somov [ 10] and by Komarov, Ponomarev and Slavjanov [ 1 ! ]. Here we report purely numerical solutions which can be directly applied to practical problems.
M. Hiyama.H. Nakamura/Cm,q~uterPhysicsComraunications103 (1997) 209-216
211
X
y,Zz)2" ~
Z
(O,O,R/2)
(O,O,-R/2)
Fig. I. Coordinates of the two-center system.
2. Two-center Coulomb functions The SchrOdinger equation for the two-center problem (see Fig. ! ) is written as
( - ½A,
Z,r,
~2) #b'(r;R)=e~t(r;R)"
(2)
where the two positive charges Zi and Z2 are located at ( x , y , z ) = ( 0 , 0 , - R / 2 ) and at ( 0 . 0 , R / 2 ) and
v = k2/2 (in atomic units) is the electron energy, rl and r2 ark the distances of the electron from Zl and Z2, res!w.ctively. [~. (2) can be separated in terms of the spheroidal coordinates, = (rl + r2)/R ~=(rl -r2)/R ,/,
(I _~ ~: < ~ ) ,
(-I __~] < I), (0 _< ~ < 27r).
(3)
The solution of Eq. (2) is expressed as
~bk(r; R) =~b,,,q(~,.'rl, dp;k,R) =
flmq(~;k,R)~mq(~l;k,R)
exp(:t=imq~) (2~)1/2
,
(4a)
where m is the magnetic quantum number and q represents the number of nodes of ~mq07; k. R). This function is normalized as
dr ~b*k( r; R)~t ( r') R)
= 8ram'8qq,8( - ~').
(4b)
The radial part,//mq(sc; k, R), and the angular part, ~,nq(V/; k, R), are the solutions of the following differential equations:
2[2
M. Hiyama, H. Nakamura/Computer Physics Communications 103 (1997) 209-216
~/'~
--
~-~ mq(~,k,R)
d-[-Amqq-C2(~2-1)-ka~-m2(~
2
l)-1117n,q(~;k,R)=O
for 1 ~ ~ < o~.,
(5)
and
(1--1/2)
-~nn(1/;k,R) - I - [ a m q q - c 2 ( l - 1 / 2 ) d-blrl--m2(l-1/2)-I],-~mq(TI;k,R)=O
for - 1 < , / < 1 ,
(6)
with
a=R(Z2 + Zi), b= R( Z2 - Zi) , c = kR/2 = (2e) I/2(R/2).
(7)
The behaviour of the normalized regular angular function -~mq(1/) at 1/-~ 4-1 is given ,,'~y q-rl)s
(1/~--1),
.--,,~(1/) =(1-1/~)'n~-~d+(1 _1/)s
(1/,"¢ l),
_'--mq(TI)=(l--~2)m/2Ed.~(l
s=0 and
(8)
with d_a:~.= dl + = O, and do~ = 1. The coefficients d~ are determined from the recurrence equations (9)
psds:t:+l +qsds~ +rsd~_, +tsdsi_2=O, where
ps=2(s+ 1 ) ( s + m + 1), qs=-[-,~q:b+m(m+
l) + s(s + 2m+ 1 ) ] ,
rs = 2c 2 ~ b,
(tO)
ts = --C2 •
At s¢ ~ I the regular radial function, H~(s¢; k, R), can be. expressed as
(ll)
//mq(s¢; k, R) = (s¢2 - 1 )m/2F(s¢; a, c, m ) , where F(s¢; a, c, in) satisfies the following equation:
t(t+2)F"+2(m+l)(t+l)F'+[c2t
2-t-(2c 2 + a ) t + a ÷ m ( m + l ) _ A ] F = 0 ,
(12)
with t = ~:-- 1. At t < 1, F can be expressed as a power series with respect to t (see Ref. [ 10] ). In the asymptotic region s~ ---, o~, the function/lmq(~:; k, R) behaves as 2
H~(~;k,R)~oo-ff~sin
(
a
c~¢+ ~cln2c~¢ -
~
+Araq
)
•
The constant ,~,~q appearing in Eqs. (5) and (6) is the separation constant.
(13)
M. Hiyama. H. NakamuralComputerPhysics Communications103 (1997~ 209-216
213
In order to obtain this separation constant it is convenient to use the following expansion of ~mq(T/; k, R) in terms of the associated Legendre polynomials: ~,,q (¢/; k, R) = exp( -ic( i - 71)) E d s ( k , R) P~m(¢/),
(14)
.~--0
with d_ I = 0 and d0 = 1. The coefficients oY~satisfy the following recursion relation:
(15)
Psds+l - Ksds + 8,d,~,_, = O,
where p s = (s + 2 m + l ) [ b + 2 i c ( s + m + I ) ] / [ 2 ( s + m) + 3 ] ,
Ks=--A+ ( s + m ) ( s + m + !), 8~ = s [ b - 2 i c ( s + m ) ] / [ 2 ( s + m ) - 1].
(16)
The separatiop, constant Amq is the solution of y(A)--O, where p081 :~
y(~) = Ko K|
(17)
•
-- ~¢2. . . .
The function y(A) is real, since it depends only on Ks and the product psSs+l. As for the initial guess of An,q, the following expansion is used: ( l 2 -m2)(b2+4c212)
Amq(k,R)~-l(l+l)-
2l(2l+1)(2l-I)
[(!*{- ! ) 2 - m 2 ] [ b 2 + 4 c 2 ( l + ! ) 2 ]
--I-
2(1+!)(2l+1)(2l+3)
'
(18)
which follows from Eq. (17) in the limit R ~ 0. With this separation constant the wave functions Hmq(~; k, R) and --,~(v/; k, R) are obtained from Eqs. (5) and (6). 3. Notes on the program
Ir~putfor total wave functions The subroutine TCOULOM has four steps. In this program all input data should he provided in the subroutine STEPO. The quantum numbers (rim and nq), the internuclear distance (Rwacl), the charges ( z l and z2), and the electron energy ( E e l e c ) should he specified. The range of R (Rin, Rout:), and the mesh size ( l ~ e s h ) are also to be defined here. Atomic units are used throughout.
Input for Am,~ The separa~.ion constant A,nq is evaluated in the subroutine STEPI at R = N x Panesh + Rin starting from [tin. The. initial guess for R = Rin is obtained from F,q. (18) and the values of Amz at grid points are calculated from Eq. (17) consecutively with use of the values at the smaller Rs. From these values at grid points, A~ at any given R is evaluated by using an interpolation procedure. The relative convergence criterion is given by ESPL, which is again defined in STEP0. The number of iteration for STEP1 is fixed at 30. When Amq is not converged, l ~ e s h or F.SPL should be changed.
214
M. Hiycma, 14. Nakamura/ Computer Physics Communications 103 (1997) 209-216
Table ! ~ters
used for the test runs
nm nq Rnucl k(Eelec) ,~) TomlCPO time
=-~(vl;0.2, 1.0)
~eA)(r/; !.0, 1.0)
~ ) ( r / ; 5.0, 1.0)
H0o(~;0.2, 1.0)
I1~)(~; 1.0, 1.0)
II~x~(~;5.0, 1.0)
0 0 1.0 0.2(0.02) -0.164017311895064 1.07 (s)
0 0 1.0 1.0(0.5) -0.317583235505185 2.10 (s)
0 0 1.0 5.0(12.5) -4.67992758609765 2.28 (s)
input for .=mq(7I) XOE (XFI/) is a starting (final) point for the variable q = tanh -~ r/ in the Numerov method. ESPE gives the accuracy of cenvergency. This represznts a gap between the forward and the backward solutions at the midpoint. meshE is the number of steps of the integral. XOE, XFE, ESPE, and meshE should be appropriately chosen for the desired accuracy, meshl/: has to be a even number and smaller than 2804.
Input for H~(6) The parameters XNFG, XAG, and HG should be defined in STEP0. HG is the mesh size of the variable ~. XNFG is the point to match the backward and forward solutions in the NUMEROV procedure to solve the differential equation. XAG is the end point in the asymptotic region at which the backward solution starts. When the energy is low, XAG has to be large. The conditions which should be satisfied by XNFG and HG are as follows: XNFG- I ---------g~+ i < 20000,
(19)
and
6 . 0 - I x 10 -15 HG
+ 2 < 20000.
(20)
Therefore I-!G(XNFG) has to be larger than 0.00031 (7.0). These restrictions are due to the array sizes declared in the programs. If necessary, these arrays can be changed larger. Table 1 shows an example ~ these constants. An example of the input data in STEP0 is shown at the end of this paper. H o w to use
When subroutine TCOULOM is called, the functions FETA('r/) and F G Z I ( ~ ) , which mean ~mq(r/; k, R) and Hn,a(~:; k, R), respectively, are obtained as the outputs. Actual values of the functions are evaluated at any given point by using the interpolation procedure. An example to use this subroutine is shown in the end of this paper. It is recommended to check the convergence of the results with respect to the parameters mentioned above.
4. Test run
Test runs are performed for, z2 = 2.0, z l = 1.0, and am = 0, nq = 0, with (Racle, k) = ( 1.0, 0.1 ), ( 1.0, 1.0), and ( 1.0, 5.0). Table ! shows the parameters of these test runs. Tables 2 and 3 are the outputs for ~ . n q ( ~ . R n u c l , k ) and for H,,= ~ ( 6 , P m u c l , k), respectively.
M. Hiyama, H. Nakamura/Computer Physics Communications 103 (1997) 209-216
Table 2 Output results of the angular wavefunetions in the test run r/ 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 -0. I -0.2 -0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9
~0o(7/;0.2, 1.0)
-~lx)(7/; 1.0, 1.0)
_--~(7/;5.0, 1.0)
0.9973162985840333 0.9557308639209367 0.9153733 ! 24477652 0.8762 ! 66670495677 0.8382343989174454 0.8014004220 i 82538 0.7656890876092364 0.731075 ! 787989325 0.697533905 i 538809 0.6650408973513078 0.6335722018778495 0.6031 042757735951 0.5736139814216167 0.5450785813824582 0.5174757332739121 0.4907834846956857 0.4649802681985568 0.4400448962986726 0.4159565565389642
0.973830149191.5026 0.9396652845455655 0.905701585684451 0.8719747374694967 0.8385194448865489 0.8053693919148546 0.7725572027268472 0.7401 i 44052630061 0.70807 ! 3972220394 0.6764574144994209 0.6453005021022747 0.6146274875657374 0.5844639568883672 0.5548342330011734 0.5257613567771637 0.4972670705855228 0.4693718043900561 0.4420946643843036 0.4154534241 549982
0.53¢,357658724455 0.6 | 40785188447631 0.68679 ! i 071489132 0.75126788347 t 5989 0.80552 ! 026556633 I 0.847880886346493 ! 0.877061535251989 0.89221044836736 0.892940 ! 89412654 0.8793407609438212 0.85197219002,5833 0.8118378297350333 0.760339740605765 0.6992183256723024 0.6304790920779232 0.5563099650358223 0.4789929635023357 0.4008 i 42431289343 0.3239765132291491
Table 3 Output results o ¢ the radial wavefunctions in the test run
1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 I 0.0 10.5 I 1.0
IhM)(¢; 0.2, 1.0)
IhN~(~; 1.0, 1.0)
IhM)(¢; 5.0, !.0)
0.2300408955125011 - 1.918704207939065E-02 -0.139711460560801 -0.174245628128448 -0.156423744679642 - 0 . I I 1397470304548 -5.691102047675129E-02 -4.619102365433256E-02 3.859 ! 626322 ! 5552E - 02 6.940111883059174E-02 8.702637491216445E-02 9.237634679879289E-02 8.73833603961281 I E - 0 2 7.44930749G774498E - 02 5.629092330463445E-02 3.52422345652089 ! E - 0 2 1.35243301060 i 004E-02 -7.069077014722884E-03 - 2 . 5 ]6839995540749E-02 - 3.982828068284425E - 0 2
0.2783491000464228 -7.548323046443275E-02 -0.224565763850092 -0.240695056862436 -0.181742429523908 -9.17 i 71353232.3798E-02 - 1.280009303193203E-03 7.076773033188092E-02 O. I 153993648132 ! 65 0.131036216347821 0.1212444747434783 9.26386 ! 869435803E-02 5.310757115900699E-02 i .042486081897740E-02 -2,874008266704371E-02 -5,942673150805303E-02 -7.864682434448724E-02 -8,53655 ! 321560798E-02 -8.026858591919255E-02 -6.53864 1345416336E-02
8.9484086930 ! 6534E-02 -0.502745828231001 -3.4530500 ! 8 3 1 2 ~ 9 E - 0 2 0.31608902382~204 0.1085347945461758 - 0 . i 98230 ! 24557762 - 0 . ! 582054087 ! 9254 9.321893562462223E-02 O. 16942227 ! 3027432 -6.780137922783706E-04 -0.143450952838185 -6.8750002 i 2997392E-02 9.02 ! 485762314| 9 6 E - 0 2 O. 1056548034584278 -2.538711776 i 3301 "lE-- 02 -0.107081771094568 -3.370204348708137E- 02 7.84 ! 328359862041E-02 7.265858q 11567263E-02 - 3. ! 872109706 ! 6360E- 02
2t5
216
M. Hiyama, tl. Nakamura/Computer Physics Communications 103 (1997) 209-216
The present code was tested to work well for the following ranges of the parameters: q -- 0 ,-,, 2, Im I = 0 ,,, 2, e = 0.001 ,,., 1.0 a.u., and R = 0.1 ,-~ 5.0 a.u. This should work all right even in the range slightly outside of these, if the, convergence is carefully tested.
Acknowledgements The present work is supported in part by a Grant-in-Aid for Scientific Research on Priority Area "Theory of Chemical Reactions" and "Quantum Tunneling o f Group of Atoms as System with Many Degrees of Freedom" from the Ministry o f Education, Science and Culture of Japan. One of the authors (M.H.) would like to thank the Japan Society for the Promotion o f Science for its support by the Fellowships for Japanese Junior Scientists. The authors are also indebted to Prof. Goodman of Univ. of Waterloo for useful discussions.
References [ ! I D.R. Bates, ER.S. Carson and T.R. Carson. Proc. Roy. Soc. London 234 (1956) 207. [2] K. Helfrich and H. Hartmann, Theor. Chim. Acta. 16 (1970) 263. [3] D.R. Bates and R.H.G. Reid, Adv. At. Mol. Phys. 4 (1968) 13. 14l See, for instance. S. Lee. M. Iwai and H. Nakamum, in: Molecules in Laser Fieldes, A.D. Bandrauk, ed. (1994) p. 217. 151 H. Takagl and H. Nakamura, J. Phys. B I I (1978) !..675. [6l MJ. Seaton. Rep. Prog. Phys. 46 (1983) 167. [71 A. Giusti, J. Phys. B 13 (1980) 3867. 18l H. Nakamura, int. Rev. Phys. Chem. 10 (1991) 123; J. Chinese Chem. Soc. 42 (1995) 359. 19l H. Takagi and H. Nakamura, Report of The Institute for Plasma Physics. Nagoya University,IIPJ-AM-16 (1980). [ 101 L.I. Ponomarev and L.N. Somov, J. Comput. Phys. 20 (1976) 183. [11] I.V. Komarov. L.I. Ponomarev and S.Yu. Slavjanov, Spheroidal and Coulomb spheroidal Functions (Nauka. Moscow, 1976) l in RussianI-
Computer Physics Communications 103 (1997) 209-216
Two-center Coulomb functions Miyabi Hiyama a,, Hiroki Nakamura a.b a Department of Functional Molecular Science. School of Mathematical and Physical Science, The Graduate University for Advanced Studies. Myodaifi, Ok~aki 444. Japan b Division of Theoretical Studies, Institute for Molecular Science, Myodaiji. Okazaki 444. Japan
Received 19 February 1997; revised 14 April 1997
Abstract Two-center Coulomb functions are crucial to describe the ionization continua of diatomic molecules, since the latter cannot be well represented by finear combinations of one-center Coulomb functions. A code to numerically evaluate the two-center angular as well as radial regular Coulomb functions is presented here. ~ ) 1997 Elsevier Science B.V. Keywords: Schr6dinger equation in spheroidal coordinates; Two-center Coulomb function; Numerical solution
PROGRAM SUMMARY
?[tle of program: TCOULOM Catalogue identifier: ADFY Program obtainable from: CPC P r o g ~ Lib,"ary,Queen's University of Belfast, N. Ireland Licensing provisions: noae Computers: HP 735 Operating systems under which the program has been tested: UNIX Programming language used: FORTRAN77 No. of bytes in distributed program, including test data. etc.: 39946
Distribution format: ASCII Keywords: Schr0dinger equation in spheroidal coordinates, twocenter Coulomb function, numerical solution Nature of physical problem Subroutine TCOULOM is a FORTRAN77 subroutine to c a k e ] ~ the regular two-center Coulomb functions for the positive onergy, 6. This program requires six constants: R, the hUemuciear ~ ; Z! and Z2, two positive charges: m and q, ~ m numbers; and e, the energy of an outgoing electron. Atomic uni~ are employed. Method of solution There are flu~e main subroutines in this pmgrara. The first one is to estimate the separation constants. The chain equation is solved there (see Ref. J I l L The other two are to o w n wavefunctions of the angular part, ~.~(v/;k,R), and of the radial paff, II,r~( ~ ; k, R ) . The spheroidal angular variable vI ranges between - I and Io and the radial vatiabie ~ ranges from l to influX. ~'~(~,;k, R) is o ~ by solving the S c ~ g e r eqea~a numerically, where we put 7/= tanhq. The hategral range o f q is
I Present address: Department of Pure and Applied Sciences, Graduate School of Arts and Sciences, The University of Tokyo, Kom~,ba, Meguro-ku, Tokyo 153, Japan.
0010-4655/97/$17.00 ~) 1997 Elsevier Science B.V. All fights reserved. PII S0010-4655 ( 97 ) 00045-3
M. Hiyama. H. Nakomura/Computer Physics Communications 103 (1997) 209-216
210
(--co, oo), but in practice the range between -4.5 and 4.5 was found to be enough. The Numerov method [21 is used to solve the differential equation. The boundary conditions are provided by pow~ series expansions .at both ends. llmq(~;k,R) is obtained by solving the corresponding SchrOW~ger equation with use of both the forward and backward Nume~,,v procedure (see Ref. 121 ).
Restrictions on the complexity of the problem This program may be used to evaluate the two-center Coulomb functions for the positive energy e. The two positive charges are
assumed to satisfy the relation Z2 _> Z:. The wavefunction for the case of Z2 = Zi also can be calculated using this program. m and q are positive integers. They correspond to magnetic and azimuthal quantum nurabers. This program produces only regular solutions both for the angalar and radial parts.
References [ I ] L.I. Pouomamv and L.N. Somov, J. Comput. Phys. 20 (1976) 183. 12] H. Takagi and H. Nakamum, Report of The Institute for Plasma Physics, Nagoya University, IIPJ-AM-16 (1980).
LONG WRITE-UP
1. Introduction Two-center Coulomb functions represent bound as well as continuum states of a diatomic molecule with one electron. Bates et al. reported the exact wave functions of the bound states of H~" and Hell 2+ [ 1]. Helfrich et al. extended the calculations to other one-electron diatomic molecules such a LiH 3+, Hell 2+, Lille 4+, and H~[2,3]. These can he used as basis functions to evaluate bound states of diatomic molecules which have more than one electron [4]. The two-center Coulomb functions at positive energy (e > 0) is also useful to represent the ionization continuum of a diatomic molecule or a long range part of the electron scattering by a diatomic molecular ion [4,5]. In the dynamic processes involving superexcited states of diatomic molecules such as photo-ionization, autoionization, Penning ionization, and dissociative recombination, the ionization continuum naturally plays an important role. The multi-channel quantum defect theory (MQDT) presents a powerful tool to analyze these processes [6-8]. The first kind of superexcited state, whose excitation energy is higher than the lowest ionization potential, is completely embedded in the ionization continuum and autoionizes even with the internuclear distance kept fixed. The autoionization rate is represented by the electronic coupling defined by
V(R, 8) = (~¢o.tlHell~a),
(1)
where Hel is the electronic Hamiltonian of the system, and t/rcont and grd are the wave functions of the ionization continuum and of the superexcited bound state, respectively. This is one of the most important parameters to characterize the first kind of superexcited state, and is the most difficult quantity to oe properly evaluated. The two-center Coulomb functions can be used to construct the wave function l/rcont and the conventional CI code enables us to evaluate the electronic coupling. In order to facilitate this it is necessary not only to numerically compute the regular two-center Coulomb functions, but also to express them in terms of Gaussian functions. The first part, i.e. the numerical solution of the regular two-center Coulomb functions, is presented here, and the second part is discussed in an accompanying paper. Takagi and Nakamura reported the numerical solutions of the two-center Coulomb functions for the homo-nuclear case (Zi = Za) [9]. Analytical expansions of the functions for the hetero-nuclear case (Zi @ Z2) are discussed by Ponomarev and Somov [ 10] and by Komarov, Ponomarev and Slavjanov [ 1 ! ]. Here we report purely numerical solutions which can be directly applied to practical problems.
M. Hiyama.H. Nakamura/Cm,q~uterPhysicsComraunications103 (1997) 209-216
211
X
y,Zz)2" ~
Z
(O,O,R/2)
(O,O,-R/2)
Fig. I. Coordinates of the two-center system.
2. Two-center Coulomb functions The SchrOdinger equation for the two-center problem (see Fig. ! ) is written as
( - ½A,
Z,r,
~2) #b'(r;R)=e~t(r;R)"
(2)
where the two positive charges Zi and Z2 are located at ( x , y , z ) = ( 0 , 0 , - R / 2 ) and at ( 0 . 0 , R / 2 ) and
v = k2/2 (in atomic units) is the electron energy, rl and r2 ark the distances of the electron from Zl and Z2, res!w.ctively. [~. (2) can be separated in terms of the spheroidal coordinates, = (rl + r2)/R ~=(rl -r2)/R ,/,
(I _~ ~: < ~ ) ,
(-I __~] < I), (0 _< ~ < 27r).
(3)
The solution of Eq. (2) is expressed as
~bk(r; R) =~b,,,q(~,.'rl, dp;k,R) =
flmq(~;k,R)~mq(~l;k,R)
exp(:t=imq~) (2~)1/2
,
(4a)
where m is the magnetic quantum number and q represents the number of nodes of ~mq07; k. R). This function is normalized as
dr ~b*k( r; R)~t ( r') R)
= 8ram'8qq,8( - ~').
(4b)
The radial part,//mq(sc; k, R), and the angular part, ~,nq(V/; k, R), are the solutions of the following differential equations:
2[2
M. Hiyama, H. Nakamura/Computer Physics Communications 103 (1997) 209-216
~/'~
--
~-~ mq(~,k,R)
d-[-Amqq-C2(~2-1)-ka~-m2(~
2
l)-1117n,q(~;k,R)=O
for 1 ~ ~ < o~.,
(5)
and
(1--1/2)
-~nn(1/;k,R) - I - [ a m q q - c 2 ( l - 1 / 2 ) d-blrl--m2(l-1/2)-I],-~mq(TI;k,R)=O
for - 1 < , / < 1 ,
(6)
with
a=R(Z2 + Zi), b= R( Z2 - Zi) , c = kR/2 = (2e) I/2(R/2).
(7)
The behaviour of the normalized regular angular function -~mq(1/) at 1/-~ 4-1 is given ,,'~y q-rl)s
(1/~--1),
.--,,~(1/) =(1-1/~)'n~-~d+(1 _1/)s
(1/,"¢ l),
_'--mq(TI)=(l--~2)m/2Ed.~(l
s=0 and
(8)
with d_a:~.= dl + = O, and do~ = 1. The coefficients d~ are determined from the recurrence equations (9)
psds:t:+l +qsds~ +rsd~_, +tsdsi_2=O, where
ps=2(s+ 1 ) ( s + m + 1), qs=-[-,~q:b+m(m+
l) + s(s + 2m+ 1 ) ] ,
rs = 2c 2 ~ b,
(tO)
ts = --C2 •
At s¢ ~ I the regular radial function, H~(s¢; k, R), can be. expressed as
(ll)
//mq(s¢; k, R) = (s¢2 - 1 )m/2F(s¢; a, c, m ) , where F(s¢; a, c, in) satisfies the following equation:
t(t+2)F"+2(m+l)(t+l)F'+[c2t
2-t-(2c 2 + a ) t + a ÷ m ( m + l ) _ A ] F = 0 ,
(12)
with t = ~:-- 1. At t < 1, F can be expressed as a power series with respect to t (see Ref. [ 10] ). In the asymptotic region s~ ---, o~, the function/lmq(~:; k, R) behaves as 2
H~(~;k,R)~oo-ff~sin
(
a
c~¢+ ~cln2c~¢ -
~
+Araq
)
•
The constant ,~,~q appearing in Eqs. (5) and (6) is the separation constant.
(13)
M. Hiyama. H. NakamuralComputerPhysics Communications103 (1997~ 209-216
213
In order to obtain this separation constant it is convenient to use the following expansion of ~mq(T/; k, R) in terms of the associated Legendre polynomials: ~,,q (¢/; k, R) = exp( -ic( i - 71)) E d s ( k , R) P~m(¢/),
(14)
.~--0
with d_ I = 0 and d0 = 1. The coefficients oY~satisfy the following recursion relation:
(15)
Psds+l - Ksds + 8,d,~,_, = O,
where p s = (s + 2 m + l ) [ b + 2 i c ( s + m + I ) ] / [ 2 ( s + m) + 3 ] ,
Ks=--A+ ( s + m ) ( s + m + !), 8~ = s [ b - 2 i c ( s + m ) ] / [ 2 ( s + m ) - 1].
(16)
The separatiop, constant Amq is the solution of y(A)--O, where p081 :~
y(~) = Ko K|
(17)
•
-- ~¢2. . . .
The function y(A) is real, since it depends only on Ks and the product psSs+l. As for the initial guess of An,q, the following expansion is used: ( l 2 -m2)(b2+4c212)
Amq(k,R)~-l(l+l)-
2l(2l+1)(2l-I)
[(!*{- ! ) 2 - m 2 ] [ b 2 + 4 c 2 ( l + ! ) 2 ]
--I-
2(1+!)(2l+1)(2l+3)
'
(18)
which follows from Eq. (17) in the limit R ~ 0. With this separation constant the wave functions Hmq(~; k, R) and --,~(v/; k, R) are obtained from Eqs. (5) and (6). 3. Notes on the program
Ir~putfor total wave functions The subroutine TCOULOM has four steps. In this program all input data should he provided in the subroutine STEPO. The quantum numbers (rim and nq), the internuclear distance (Rwacl), the charges ( z l and z2), and the electron energy ( E e l e c ) should he specified. The range of R (Rin, Rout:), and the mesh size ( l ~ e s h ) are also to be defined here. Atomic units are used throughout.
Input for Am,~ The separa~.ion constant A,nq is evaluated in the subroutine STEPI at R = N x Panesh + Rin starting from [tin. The. initial guess for R = Rin is obtained from F,q. (18) and the values of Amz at grid points are calculated from Eq. (17) consecutively with use of the values at the smaller Rs. From these values at grid points, A~ at any given R is evaluated by using an interpolation procedure. The relative convergence criterion is given by ESPL, which is again defined in STEP0. The number of iteration for STEP1 is fixed at 30. When Amq is not converged, l ~ e s h or F.SPL should be changed.
214
M. Hiycma, 14. Nakamura/ Computer Physics Communications 103 (1997) 209-216
Table ! ~ters
used for the test runs
nm nq Rnucl k(Eelec) ,~) TomlCPO time
=-~(vl;0.2, 1.0)
~eA)(r/; !.0, 1.0)
~ ) ( r / ; 5.0, 1.0)
H0o(~;0.2, 1.0)
I1~)(~; 1.0, 1.0)
II~x~(~;5.0, 1.0)
0 0 1.0 0.2(0.02) -0.164017311895064 1.07 (s)
0 0 1.0 1.0(0.5) -0.317583235505185 2.10 (s)
0 0 1.0 5.0(12.5) -4.67992758609765 2.28 (s)
input for .=mq(7I) XOE (XFI/) is a starting (final) point for the variable q = tanh -~ r/ in the Numerov method. ESPE gives the accuracy of cenvergency. This represznts a gap between the forward and the backward solutions at the midpoint. meshE is the number of steps of the integral. XOE, XFE, ESPE, and meshE should be appropriately chosen for the desired accuracy, meshl/: has to be a even number and smaller than 2804.
Input for H~(6) The parameters XNFG, XAG, and HG should be defined in STEP0. HG is the mesh size of the variable ~. XNFG is the point to match the backward and forward solutions in the NUMEROV procedure to solve the differential equation. XAG is the end point in the asymptotic region at which the backward solution starts. When the energy is low, XAG has to be large. The conditions which should be satisfied by XNFG and HG are as follows: XNFG- I ---------g~+ i < 20000,
(19)
and
6 . 0 - I x 10 -15 HG
+ 2 < 20000.
(20)
Therefore I-!G(XNFG) has to be larger than 0.00031 (7.0). These restrictions are due to the array sizes declared in the programs. If necessary, these arrays can be changed larger. Table 1 shows an example ~ these constants. An example of the input data in STEP0 is shown at the end of this paper. H o w to use
When subroutine TCOULOM is called, the functions FETA('r/) and F G Z I ( ~ ) , which mean ~mq(r/; k, R) and Hn,a(~:; k, R), respectively, are obtained as the outputs. Actual values of the functions are evaluated at any given point by using the interpolation procedure. An example to use this subroutine is shown in the end of this paper. It is recommended to check the convergence of the results with respect to the parameters mentioned above.
4. Test run
Test runs are performed for, z2 = 2.0, z l = 1.0, and am = 0, nq = 0, with (Racle, k) = ( 1.0, 0.1 ), ( 1.0, 1.0), and ( 1.0, 5.0). Table ! shows the parameters of these test runs. Tables 2 and 3 are the outputs for ~ . n q ( ~ . R n u c l , k ) and for H,,= ~ ( 6 , P m u c l , k), respectively.
M. Hiyama, H. Nakamura/Computer Physics Communications 103 (1997) 209-216
Table 2 Output results of the angular wavefunetions in the test run r/ 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 -0. I -0.2 -0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9
~0o(7/;0.2, 1.0)
-~lx)(7/; 1.0, 1.0)
_--~(7/;5.0, 1.0)
0.9973162985840333 0.9557308639209367 0.9153733 ! 24477652 0.8762 ! 66670495677 0.8382343989174454 0.8014004220 i 82538 0.7656890876092364 0.731075 ! 787989325 0.697533905 i 538809 0.6650408973513078 0.6335722018778495 0.6031 042757735951 0.5736139814216167 0.5450785813824582 0.5174757332739121 0.4907834846956857 0.4649802681985568 0.4400448962986726 0.4159565565389642
0.973830149191.5026 0.9396652845455655 0.905701585684451 0.8719747374694967 0.8385194448865489 0.8053693919148546 0.7725572027268472 0.7401 i 44052630061 0.70807 ! 3972220394 0.6764574144994209 0.6453005021022747 0.6146274875657374 0.5844639568883672 0.5548342330011734 0.5257613567771637 0.4972670705855228 0.4693718043900561 0.4420946643843036 0.4154534241 549982
0.53¢,357658724455 0.6 | 40785188447631 0.68679 ! i 071489132 0.75126788347 t 5989 0.80552 ! 026556633 I 0.847880886346493 ! 0.877061535251989 0.89221044836736 0.892940 ! 89412654 0.8793407609438212 0.85197219002,5833 0.8118378297350333 0.760339740605765 0.6992183256723024 0.6304790920779232 0.5563099650358223 0.4789929635023357 0.4008 i 42431289343 0.3239765132291491
Table 3 Output results o ¢ the radial wavefunctions in the test run
1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 I 0.0 10.5 I 1.0
IhM)(¢; 0.2, 1.0)
IhN~(~; 1.0, 1.0)
IhM)(¢; 5.0, !.0)
0.2300408955125011 - 1.918704207939065E-02 -0.139711460560801 -0.174245628128448 -0.156423744679642 - 0 . I I 1397470304548 -5.691102047675129E-02 -4.619102365433256E-02 3.859 ! 626322 ! 5552E - 02 6.940111883059174E-02 8.702637491216445E-02 9.237634679879289E-02 8.73833603961281 I E - 0 2 7.44930749G774498E - 02 5.629092330463445E-02 3.52422345652089 ! E - 0 2 1.35243301060 i 004E-02 -7.069077014722884E-03 - 2 . 5 ]6839995540749E-02 - 3.982828068284425E - 0 2
0.2783491000464228 -7.548323046443275E-02 -0.224565763850092 -0.240695056862436 -0.181742429523908 -9.17 i 71353232.3798E-02 - 1.280009303193203E-03 7.076773033188092E-02 O. I 153993648132 ! 65 0.131036216347821 0.1212444747434783 9.26386 ! 869435803E-02 5.310757115900699E-02 i .042486081897740E-02 -2,874008266704371E-02 -5,942673150805303E-02 -7.864682434448724E-02 -8,53655 ! 321560798E-02 -8.026858591919255E-02 -6.53864 1345416336E-02
8.9484086930 ! 6534E-02 -0.502745828231001 -3.4530500 ! 8 3 1 2 ~ 9 E - 0 2 0.31608902382~204 0.1085347945461758 - 0 . i 98230 ! 24557762 - 0 . ! 582054087 ! 9254 9.321893562462223E-02 O. 16942227 ! 3027432 -6.780137922783706E-04 -0.143450952838185 -6.8750002 i 2997392E-02 9.02 ! 485762314| 9 6 E - 0 2 O. 1056548034584278 -2.538711776 i 3301 "lE-- 02 -0.107081771094568 -3.370204348708137E- 02 7.84 ! 328359862041E-02 7.265858q 11567263E-02 - 3. ! 872109706 ! 6360E- 02
2t5
216
M. Hiyama, tl. Nakamura/Computer Physics Communications 103 (1997) 209-216
The present code was tested to work well for the following ranges of the parameters: q -- 0 ,-,, 2, Im I = 0 ,,, 2, e = 0.001 ,,., 1.0 a.u., and R = 0.1 ,-~ 5.0 a.u. This should work all right even in the range slightly outside of these, if the, convergence is carefully tested.
Acknowledgements The present work is supported in part by a Grant-in-Aid for Scientific Research on Priority Area "Theory of Chemical Reactions" and "Quantum Tunneling o f Group of Atoms as System with Many Degrees of Freedom" from the Ministry o f Education, Science and Culture of Japan. One of the authors (M.H.) would like to thank the Japan Society for the Promotion o f Science for its support by the Fellowships for Japanese Junior Scientists. The authors are also indebted to Prof. Goodman of Univ. of Waterloo for useful discussions.
References [ ! I D.R. Bates, ER.S. Carson and T.R. Carson. Proc. Roy. Soc. London 234 (1956) 207. [2] K. Helfrich and H. Hartmann, Theor. Chim. Acta. 16 (1970) 263. [3] D.R. Bates and R.H.G. Reid, Adv. At. Mol. Phys. 4 (1968) 13. 14l See, for instance. S. Lee. M. Iwai and H. Nakamum, in: Molecules in Laser Fieldes, A.D. Bandrauk, ed. (1994) p. 217. 151 H. Takagl and H. Nakamura, J. Phys. B I I (1978) !..675. [6l MJ. Seaton. Rep. Prog. Phys. 46 (1983) 167. [71 A. Giusti, J. Phys. B 13 (1980) 3867. 18l H. Nakamura, int. Rev. Phys. Chem. 10 (1991) 123; J. Chinese Chem. Soc. 42 (1995) 359. 19l H. Takagi and H. Nakamura, Report of The Institute for Plasma Physics. Nagoya University,IIPJ-AM-16 (1980). [ 101 L.I. Ponomarev and L.N. Somov, J. Comput. Phys. 20 (1976) 183. [11] I.V. Komarov. L.I. Ponomarev and S.Yu. Slavjanov, Spheroidal and Coulomb spheroidal Functions (Nauka. Moscow, 1976) l in RussianI-