- Email: [email protected]
Basis set limitations on the ab initio calculation of stopping cross-sections via generalized oscillator strengths
Chemical Physics 309 (2005) 89–94 www.elsevier.com/locate/chemphys
Basis set limitations on the ab initio calculation of stopping cross-sections via generalized oscillator strengths J.A. Nobel a, S.B. Trickey a
a,*
, John R. Sabin a, Jens Oddershede
b
Quantum Theory Project, Departments of Physics and of Chemistry, University of Florida, P.O. Box 118435, Gainesville, FL 32611-8435, USA b Kemisk Institut, Syddansk Universitet, Odense, Denmark Received 7 May 2004; accepted 26 May 2004 Available online 3 August 2004
Abstract Ab initio treatment of the first Born electronic stopping cross-section for a bare proton incident upon atomic Be is attempted via a large-basis (15s13p6d3f contracted to 9s13p6d3f) polarization propagator calculation of the generalized oscillator strength distribution (GOSD). The calculated GOSD fails to satisfy the Bethe sum rule for moderately high values of collisional momentum transfer in spite of asymptotic augmentation of the calculated GOSD by a hydrogen-like contribution. A simple model shows that the BSR violation introduces unsystematic spurious contributions to the stopping cross-section. The failure illustrates a shortcoming in the conventional approach to optimizing finite basis sets and shows that predictive calculation of electronic stopping in real systems requires continued development and use of judiciously selected approximate methods. Ó 2004 Elsevier B.V. All rights reserved.
1. Salutation Our colleague Notker R€ osch has worked in many aspects of the theory and computation of molecules and clusters. On the occasion of his 60th birthday, we salute his intensity and creativity in extending and applying quantum chemical methods to challenging systems. This note focuses on an aspect of a problem that he also has studied, namely the selection and optimization of finite basis sets or the use of unusual functions in order to obtain high-fidelity treatment of important molecular properties [1,2]. The problem we treat illustrates that, for a relatively simple response property, more ingenious modes of basis augmentation than addition of functions and optimization of their exponents must be found. 2. Motivation Scattering of charged particles from atoms, clusters, surfaces, or solids is a ubiquitous mode of probing *
Corresponding author. Tel.: +352-392-1597; fax: +352-392-8722. E-mail address: [email protected]fl.edu (S.B. Trickey).
0301-0104/$ - see front matter Ó 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.chemphys.2004.05.027
and/or altering the system properties. At sufficiently low energies, such scattering is the essence of many chemical reaction steps. At higher energies, such scattering is basic to ion implantation, radiation damage, various radiation therapies, and particle detection in nuclear and high-energy physics. As in all many-body problems, the question arises as to which of these scattering phenomena can be treated accurately from well-defined approximations to first-principles descriptions. The predominant approach to electronic structure problems in both quantum chemistry and materials physics is to select a specific approximation and convert the problem into a set of linear-algebra tasks by introduction of a one-particle basis. In quantum chemistry, the basis of choice is Gaussian-type orbitals. Customarily the consequences of practical restriction to a finite basis are controlled by exploiting the variation principle: the size of the basis, the number of different symmetries, and the orbital exponents are adjusted to drive some fiduciary energy (e.g., the isolated atom energy) down to a specified accuracy. Non-variational quantities such as response properties (susceptibilities, cross-sections, etc.) are another matter. A typical approach is to refine the basis in order to satisfy a set of relevant sum rules to
90
J.A. Nobel et al. / Chemical Physics 309 (2005) 89–94
satisfactory precision. Here, we revisit a particular response problem and show that there still is no systematic way of refining a basis to satisfy the relevant sum rules. We view this as a major unsolved challenge in computational molecular physics.
3. Summary of physical problem and methods In the first Born approximation, the stopping crosssection SðvÞ of a structureless point particle of charge Z1 and velocity ðvÞ incident on a target of atomic number Z2 gives the well-known dependence upon generalized oscillator strengths Fn0 ðqÞ first considered by Bethe and reviewed at length by Inokuti [3,5,6]. By definition, the electronic stopping cross-section is Z 1 dE dr ¼ dE E: ð1Þ SðvÞ ¼ N dx dE Here N and dr=dE are the material density and the total inelastic differential cross-section for electronic excitations, respectively. In Hartree atomic units (energies are in Hartree1 EH ¼ 27:2116 eV and lengths are in Bohr Bethe’s expression is radii: 1 a.u. ¼ 0.529177 A), ( Ip Z 4pZ12 X qmax dq SðvÞ ¼ Fn 0ðqÞ q v2 En0 =v n6¼0 ) Z qmax Z 1 dF ðq; EÞ dq dE þ : ð2Þ dE q E=v Ip The two contributions are from the discrete and continuum states, respectively, En0 is the usual energy difference between the ground and nth state, and Ip , qmax are the system ionization potential and the maximum momentum transfer [determined from the collision kinematics; see [7], Eqs. (17) and (18), respectively]. The upper limit of Ip on the sum is shorthand to indicate inclusion of terms up to a state energy difference of that amount. A general review of molecular stopping crosscalculations from oscillator strength distributions is in [4]. The exact result, Eq. (2) is relatively simple. In view of that simplicity and early calculations of atomic He generalized oscillator strengths (hereafter GOS),[8–10], there might seem to be little need for the continued development and use of various approximate schemes (examples may be found in [11,12] and references therein) for the calculation of SðvÞ in the first Born approximation. Such optimism would be tempered a bit by the need to include all excited states, both discrete and continuum, in Eq. (2). However, the issue was addressed in [7] at least for He. Those authors used a hydrogenic large-q GOS distribution to approximate the large-q contributions from high-lying states which are represented poorly or not at all by the finite basis set.
Another concern would be whether the ab-initio propagator techniques and large-q approximations used in [7] (or alternatively, the explicitly correlated wave functions used in [8–10]) could be utilized to treat the more complicated systems (solids, thin films, polymers, etc.) typically treated with more approximate methods. The issue is simple to state: are modern ab initio techniques for calculating stopping feasible for many-electron systems? This question is more-or-less unaddressed in the literature even for atoms beyond He, let alone more complicated systems. A simple, relevant test therefore is the modest step of treating the Be atom by direct application of the techniques of [7]. Note that though the Li atom might seem the logical next step, it is in fact rather less interesting. As an alkali its dominant high energy excitations will be hydrogenic. That fact makes Li well-suited to the methodology of [7] a priori, hence does not provide as serious a test as does Be. A brief summary of the methodology in question helps both in framing the issue and in defining notation. For transitions between initial and final states j0i and jni with energy difference En0 ¼ En E0 , the GOS is written as Fn0 ¼
2En0 2 jh0jeiqr jni j q2
with the second-quantized position operator g X ri;j ayi aj : r¼
ð3Þ
ð4Þ
i;j¼1
Here g is the total number of basis functions, and the matrix representation of the position operator in that basis is rij ¼
Z2 X hijr‘ jji:
ð5Þ
‘¼1
For the choice of operators A ¼ eiqr ; B ¼ eiqr
ð6Þ
the poles and residues, respectively, of the polarization propagator X h0jAjnihnjBj0i h0jBjnihnjAj0i hhA; BiiE ¼ ð7Þ E En þ E 0 E En þ E0 n6¼0 (with the sum over all excited states) correspond to the excitation energies and generalized transition moments necessary to calculate the GOS distribution (GOSD). All the relevant details can be found in [13]. Actual implementation for this study (except for the choice of basis set of course; see below) was the same as in [7], i.e., the calculations were done using the MUNICH [14] program package and within the random phase approximation (RPA).
J.A. Nobel et al. / Chemical Physics 309 (2005) 89–94
4. Results and discussion As customary in modern molecular calculations, a contracted Gaussian–type orbital (CGTO) basis was used. It was constructed starting with the CGTO (15s,9p,5d) ! [9s,9p,5d] basis set of Graham et al. [15] With that basis, those authors found that a full-configuration-interaction treatment gave a very good description of the low-lying excitations of the Be atom. Since higher-lying excitations are relevant to our problem, we added four primitive p-functions, one primitive d-function, and a new f-function sub-manifold of three primitives. From these, the basis was contracted as (15s13p6d3f) ! [9s13p6d3f]. Exponents and contraction coefficients are provided in Table 1. The additional diffuse orbitals were chosen in an attempt to improve fulfillment of the Bethe sum rule; see discussion below. The total closed-shell HF energy from both the original and augmented bases was )14.572608 EH compared to )14.572543 EH quoted in [15] and a HF limit Table 1 Basis exponents and contraction coefficients ‘-type exponent coefficients s
s
s s
p
p
d
f
3630.0 532.31 117.80 32.660 10.480 3.6680 3.6680 1.3540 0.38900 0.15020 0.05241 0.02000 0.00800 0.00300 0.00100 154.000 31.2234 6.71000 1.44200 0.41030 0.13970 0.04922 0.01600 0.00700 0.00300 0.00100 0.00050 0.00010 0.65000 0.20000 0.07000 0.01500 0.00320 0.00060 0.65000 0.07000 0.00320
0.000839 0.006735 0.035726 0.138635 0.385399 0.547688 0.213406 0.814692 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
91
[16] of )14.573023 EH . For a complete basis set, the Thomas–Reiche–Kuhn (TRK) sum rule for the dipole oscillator strengths fn0 X fn0 ¼ Z2 ð8Þ n6¼0
is exact in the RPA. Since this is the small q limit of the GOSD, the calculated value of the TRK sum is a relevant indicator of the quality of a finite basis set. A detailed discussion of the significance of the TRK sum rule in evaluating RPA stopping calculations is given by Cabrera-Trujillo et al. [4]. The essential point is that we obtain 4.0077865 (3.9923177) for atomic Be in the length (velocity) formulation, so that the augmented basis set appears to be quite reasonable in quality, for this particular sum rule. A basis of very high quality by usual quantum chemistry standards is, nevertheless, severely challenged in the calculation of the stopping cross-section. Details are in [4,7]. Briefly, any finite basis has, explicitly or implicitly, an upper limit on the atomic orbital angular momentum which it describes. Whatever that limit, a sufficiently large but finite value of q always can be found such that the expansion of Eq. (6) in powers of q has significant multipoles of higher angular momentum than that limit. The problem manifests itself quite strongly in the failure of the calculated Bethe sum rule (BSR) to behave as it should. The exact result [5] is X Fn0 ðqÞ ¼ Z2 ð9Þ n6¼0
for all values of q. The calculations on He [7] showed, however, that the calculated value of the sum in Eq. (9) is acceptably close to Z2 only for q 6 3 a.u. qL . For larger q the calculated sum falls rapidly to zero, a dramatic sign of the finite basis set having been overrun. The approximation used in [7] to eliminate this largeq misbehavior was to append an analytic hydrogenic GOSD to the RPA-propagator GOSD, with a suitably chosen value of q for the cut-over. Thus, Holm-Mortensen et al. [7] approximated the BSR (9) as Z 1 EL X dF ðq; EÞ Fn0 ðqÞ þ dE Z2 ð10Þ dE hyd EL n6¼¼0 with EL q2L =2 adjusted to make the approximation as good as possible. The hydrogenic form for dF j chosen dE hyd in [7] was from Ketkar and Bonham [17]. Fig. 2 of [7] shows that, for atomic He, the value EL ¼ 137 eV ) qL 3:2 a.u. provides a calculated GOSD which satisfies Eq. (10) to within a few percent for all q. Holm-Mortensen et al. [7] argued that the augmentation of the GOSD with a high-energy hydrogenic part is a reasonable procedure because the critical quantity Eq. (6) is a one-electron operator and the
92
J.A. Nobel et al. / Chemical Physics 309 (2005) 89–94
hydrogenic GOSD for large q therefore is physically sensible. The appending of the hydrogenic GOSD is equivalent logically (though not operationally) with selecting an approximate complement to the finite Gaussian basis and then restricting the propagator calculation solely to hydrogenic excitations within that complementary space. Here we test whether that argument is in any sense general, since in the Be atom there are manyelectron correlations not found in the He atom studied in [7]. The merit of choosing Be rather than Li should be clear, but perhaps it is worth adding that atomic Be is notorious (among quantum chemists) as a test in its own right because of the near-degeneracy of its ground state: so-called static correlation. For the test calculations the augmented BSR Eq. (10) was generalized to read Z 1( EL X dF ðq; EÞ Z2 Fn0 ðqÞ þ hðE E2s Þ dE hyd;2s EL n6¼¼0 ) dF ðq; EÞ þ hðE E1s Þ dE; ð11Þ dE hyd;1s where the Heaviside step functions turn on the correspondingly parameterized hydrogenic GOSDs for energies above some suitably selected excitation thresholds Enl . The most straightforward test case, in the sense of building in the least additional arbitrariness, is the one actually used, namely Enl ¼ EL , so that the entire hydrogenic augmentation came into play at once. The parameters of the Ketkar–Bonham model were chosen (using the same notation as in [7]) to be best-available atomic values; see Table 2. As before, EL was deter-
Table 2 Parameters used for the Ketkar-Bonham hydrogen-like GOSD Parameter (units) a
Ip (eV) n (a.u., inverse length)b k (a.u., charge)b a b
1s value
2s value
123.33 3.6848 0.956
9.16 3.6848 0.956
Ref. [23]. Ref. [24].
mined empirically to make Eq. (11) as close to an equality as possible. The result for the best value we could obtain EL ¼ 193:2 eV is in Fig. 1 (which is to be compared with Fig. 2 of [7]). Even the best value for EL results in a calculated BSR which falls (at q 2:45 a.u.) as much as 45% below the proper constant, Z2 . We were unable to find a combination of basis and EL value which does any better than this clearly unsatisfactory outcome. Clearly the appending procedure is unsuccessful for the Be atom, hence the procedure is not general. A brief calibration of the quality of the underlying propagator calculation may serve to reinforce the severity of this finding. Consider the excitation energies for the low-lying singlet transitions. Table 3 shows the present results, the RPA results from a smaller (9s9p5d) basis [18], and the experimental results.[19] The present results extend the earlier one in only one crucial way, namely the transition to the 4f state which could not be treated in the basis of [18], which had no f-type functions. Transitions to this state affect the GOSD and BSR only for q 6¼ 0. However, this is the domain in which improvement is needed. The differences between RPA and experimental excitation energies listed in Table 3
4
3.8
3.6
Calculated Z2 (q)
3.4
3.2
3
2.8
2.6
2.4
2.2 0
2
4
6
8
10
12
q (au)
Fig. 1. Calculated Z2 ðqÞ from Eq. (11).
14
16
18
20
J.A. Nobel et al. / Chemical Physics 309 (2005) 89–94 Table 3 Low-lying singlet excitation energies (eV) for the Be atom: comparison of the current calculation with a previous, small-basis RPA calculation [18] and with experiment [19] Final state
This calc.
Smaller basis RPA
Expt.
2s2p; 1 P0 2s3s; 1 S 2s3p; 1 P0 2s3d; 1 D 2s4s; 1 S 2s4p; 1 P0 2s4d; 1 D 2s4f; 1 F 2s5s; 1 S
4.80 6.12 6.72 6.94 7.26 7.48 7.59 7.68 7.72
4.80 6.12 6.72 6.94 7.26 7.48 7.59 – 7.72
5.28 6.78 7.46 7.99 8.09 8.31 8.53 Forbidden 8.60
For the forbidden 2s4f; 1 F our calculation gives zero dipole oscillator strength.
could be eliminated by use of higher-order polarization propagator approximations.[18] Such improvements in excitation energies would be of little consequence for the issue addressed here, namely the behavior of the calculated GOSD (which is determined primarily by the limitations of the finite basis set). It might be argued that a calculated GOSD that failed the BSR still could give a decent stopping cross-section. The BSR after all is not an integrated quantity whereas SðvÞ is. A qualitative argument shows that this expectation cannot be true in general. First, the polarization propagator produces a spectrum corresponding to the bound transitions of the system plus a set of pseudo-transitions to states that are the best representation of continuum states possible in a L2 basis. Within a finite L2 basis therefore, the stopping cross-section Eq. (2) has the approximate form nmax Z qmax 4pZ12 X dq ð12Þ SðvÞ ¼ Fn0 ðqÞ : q v2 n6¼0 En0 =v The limit nmax is determined by the basis set size. If it were the case that the calculated generalized oscillator strength sum were to saturate at or before nmax and fulfill the BSR, then that sum could be moved outside the integral (essentially as happens in the dipole approximation) and the result would be a modified mean excitation and Bethe logarithm. Failure to fulfill the BSR, on the other hand, introduces a spurious q-dependence that cannot be pulled out of the integral. We can model that dependence roughly (see Fig. 1) as nmax X
Fn0 ðqÞ ¼ Z2 ;
q 6 ðq0 D=2Þ;
n6¼0
4 2 ðq q Þ 1 ; ¼ Z2 1 þ b 0 D2 ðq0 D=2Þ 6 q 6 ðq0 þ D=2Þ; ¼ Z2 ; q P ðq0 þ D=2Þ:
ð13Þ
93
Note that the location q0 , width D, and depth b of the deficiency region in the calculated BSR all are velocitydependent, depending on the details of the basis. Insertion of this form into Eq. (12) and integration gives terms quadratic, linear, and logarithmic in q. Since the kinematic limits of q have velocity dependence of v1 and v, the net effect of the BSR violation is to introduce spurious v-dependence into SðvÞ that can range from v2 to v2 in addition to shifting the logarithmic term. Therefore the gross qualitative failure to fulfill the BSR such as we find for atomic Be is certain to have highly basis-set-specific effects on SðvÞ. In view of that specificity, we did not bother to compute SðvÞ for the particular basis used here since nothing general would have been learned.
5. Conclusions Though the pseudo-excitations from the polarization propagator treated in a finite L2 basis have no physical significance [20], the spectrum calculated with a highquality basis does satisfy the basic sum rule, namely TRK. There is, however, a severe finite basis problem for excitations at high momentum transfer. Striking evidence is in the failure of the calculated BSR, Eq. (9), to be independent of the value of momentum transferred. Augmentation of the calculated GOSD with simple hydrogenic distributions above some critical q, a successful technique for He, fails completely to resolve the problem for Be. The challenge involved in resolving the problem is easy to understand. The finite basis produces a kind of trench in the Bethe surface defined by the GOSD. The trench is roughly parallel with the E axis, hence over some finite range of q the BSR is not satisfied. A successful augmentation procedure would have to fill that trench so as to satisfy the BSR without building in spurious traits in the augmented GOSD. It is this challenge which suggests that attempts to tune various parameters, e.g., the Enl , in the model tested here would not be a good idea. Such tuning could well lead to a wholly unphysical augmented Bethe surface which nonetheless satisfied the BSR. In particular, the fulfillment of Eq. (9) is no guarantee of a correct GOSD, especially if that fulfillment is achieved by grafting on hydrogenic behavior for relatively low Enl . The calculations presented here are a state-of-the-art treatment of the GOSD. Since it is fundamental to firstprinciples calculation of electronic stopping in the first Born approximation, these calculations present a sobering challenge. Occasionally it has been suggested to the authors that approximate treatments of electronic stopping (such as those studied by us and others) no longer merit attention. Rather, the argument goes, one should proceed straightway with a proper GOSD
94
J.A. Nobel et al. / Chemical Physics 309 (2005) 89–94
calculation. The present calculation shows otherwise. If the state-of-the-art cannot treat a free Be atom adequately, it is clearly unrealistic to expect predictive GOSD calculations of electronic stopping for polymers, thin films, etc. In particular, the current status of timedependent density functional theory calculations for extended systems (see references in [21]) or dielectric treatments (RPA plus DFT orbitals and energies; see Mathar et al. [22]) makes clear that it is quite certain that even the level of refinement and generality reached in the present work is unattainable for those systems. Continued development and use of approximate stopping theories for systems more complicated than isolated atoms therefore seems not only warranted but inescapable. Acknowledgements We thank Stephan P.A. Sauer for useful conversations. J.R.S. and S.B.T. were supported in part by the US Army Research Office under Contract DAA L03-91G-0119. S.B.T. also was supported in part by US NSF Grant DMR-0218957. J.R.S. also was supported in part by US NSF Grant INT-9016299. J.O. was supported in part by the Danish Natural Science Research Council (Grant No. 11-0924). References [1] V. Nasluzov, N. R€ osch, unpublished. [2] I. Wilhelmy, N. R€ osch, Chem. Phys. 185 (1994) 317. [3] H. Bethe, Ann. Phys. (Leipzig) 5 (1930) 325.
[4] R. Cabrera-Trujillo, J.R. Sabin, J. Oddershede, Adv. Quantum Chem. 46 (2004) 121. [5] M. Inokuti, Rev. Mod. Phys. 43 (1971) 297. [6] M. Inokuti, Y. Itikawa, J.E. Turner, Rev. Mod. Phys. 50 (1978) 23. [7] E. Holm-Mortensen, J. Oddershede, J.R. Sabin, Nucl. Instrum. Meth. B 69 (1992) 24. [8] Y.K. Kim, M. Inokuti, Phys. Rev. 175 (1968) 176. [9] Y.K. Kim, M. Inokuti, Phys. Rev. 181 (1969) 205. [10] Y.K. Kim, M. Inokuti, Phys. Rev. 184 (1969) 38. [11] J. Wang, R.J. Mathar, S.B. Trickey, J.R. Sabin, J. Phys.: Condens. Matter 11 (1999) 3973. [12] D.E. Meltzer, J.R. Sabin, S.B. Trickey, J.Z. Wu, Nucl. Instrum. Meth. B 82 (1993) 493. [13] J. Oddershede, Adv. Chem. Phys. 69 (1987) 201. [14] G.H.F. Diercksen, W.P. Kraemer, Reference Manual to the MUNICH System of Programs, Max Planck Inst. f€ ur Astrophysik, Special Report, 1981. [15] R.L. Graham, D.L. Yeater, J. Olsen, P. Jorgensen, R. Harrison, S. Zarrabian, R.J. Bartlett, J. Chem. Phys. 85 (1986) 6544. [16] C. Froese-Fischer, The Hartree–Fock Method for Atoms, Wiley, New York, 1977. [17] N. Ketkar, R.A. Bonham, J. Chem. Phys. 84 (1986) 6091. [18] J. Geertsen, J. Oddershede, G.E. Scuseria, Int. J. Quantum Chem. S21 (1987) 475. [19] S. Bashkin, J.A. Stoner, Atomic Energy Levels Grotrian Diagrams, North-Holland, Amsterdam, 1975. [20] J. Geertsen, J. Oddershede, J.R. Sabin, Phys. Rev. A 34 (1986) 1104. [21] H. Appel, E.K.U. Gross, K. Burke, Phys. Rev. Lett. 90 (2003) 043005; M.E. Casida, D.R. Salahub, J. Chem. Phys. 113 (2000) 8918. [22] R.J. Mathar, S.B. Trickey, J.R. Sabin, Adv. Quantum Chem. 45 (2004) 277; Nucl. Instrum. Meth. B 155 (1999) 249. [23] D.A. Shirley, R.L. Martin, S.P. Kowalczyk, F. McFeely, L. Ley, Phys. Rev. B 15 (1977) 544. [24] E. Clementi, D.L. Raimundi, J. Chem. Phys. 38 (1963) 2686.
Basis set limitations on the ab initio calculation of stopping cross-sections via generalized oscillator strengths J.A. Nobel a, S.B. Trickey a
a,*
, John R. Sabin a, Jens Oddershede
b
Quantum Theory Project, Departments of Physics and of Chemistry, University of Florida, P.O. Box 118435, Gainesville, FL 32611-8435, USA b Kemisk Institut, Syddansk Universitet, Odense, Denmark Received 7 May 2004; accepted 26 May 2004 Available online 3 August 2004
Abstract Ab initio treatment of the first Born electronic stopping cross-section for a bare proton incident upon atomic Be is attempted via a large-basis (15s13p6d3f contracted to 9s13p6d3f) polarization propagator calculation of the generalized oscillator strength distribution (GOSD). The calculated GOSD fails to satisfy the Bethe sum rule for moderately high values of collisional momentum transfer in spite of asymptotic augmentation of the calculated GOSD by a hydrogen-like contribution. A simple model shows that the BSR violation introduces unsystematic spurious contributions to the stopping cross-section. The failure illustrates a shortcoming in the conventional approach to optimizing finite basis sets and shows that predictive calculation of electronic stopping in real systems requires continued development and use of judiciously selected approximate methods. Ó 2004 Elsevier B.V. All rights reserved.
1. Salutation Our colleague Notker R€ osch has worked in many aspects of the theory and computation of molecules and clusters. On the occasion of his 60th birthday, we salute his intensity and creativity in extending and applying quantum chemical methods to challenging systems. This note focuses on an aspect of a problem that he also has studied, namely the selection and optimization of finite basis sets or the use of unusual functions in order to obtain high-fidelity treatment of important molecular properties [1,2]. The problem we treat illustrates that, for a relatively simple response property, more ingenious modes of basis augmentation than addition of functions and optimization of their exponents must be found. 2. Motivation Scattering of charged particles from atoms, clusters, surfaces, or solids is a ubiquitous mode of probing *
Corresponding author. Tel.: +352-392-1597; fax: +352-392-8722. E-mail address: [email protected]fl.edu (S.B. Trickey).
0301-0104/$ - see front matter Ó 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.chemphys.2004.05.027
and/or altering the system properties. At sufficiently low energies, such scattering is the essence of many chemical reaction steps. At higher energies, such scattering is basic to ion implantation, radiation damage, various radiation therapies, and particle detection in nuclear and high-energy physics. As in all many-body problems, the question arises as to which of these scattering phenomena can be treated accurately from well-defined approximations to first-principles descriptions. The predominant approach to electronic structure problems in both quantum chemistry and materials physics is to select a specific approximation and convert the problem into a set of linear-algebra tasks by introduction of a one-particle basis. In quantum chemistry, the basis of choice is Gaussian-type orbitals. Customarily the consequences of practical restriction to a finite basis are controlled by exploiting the variation principle: the size of the basis, the number of different symmetries, and the orbital exponents are adjusted to drive some fiduciary energy (e.g., the isolated atom energy) down to a specified accuracy. Non-variational quantities such as response properties (susceptibilities, cross-sections, etc.) are another matter. A typical approach is to refine the basis in order to satisfy a set of relevant sum rules to
90
J.A. Nobel et al. / Chemical Physics 309 (2005) 89–94
satisfactory precision. Here, we revisit a particular response problem and show that there still is no systematic way of refining a basis to satisfy the relevant sum rules. We view this as a major unsolved challenge in computational molecular physics.
3. Summary of physical problem and methods In the first Born approximation, the stopping crosssection SðvÞ of a structureless point particle of charge Z1 and velocity ðvÞ incident on a target of atomic number Z2 gives the well-known dependence upon generalized oscillator strengths Fn0 ðqÞ first considered by Bethe and reviewed at length by Inokuti [3,5,6]. By definition, the electronic stopping cross-section is Z 1 dE dr ¼ dE E: ð1Þ SðvÞ ¼ N dx dE Here N and dr=dE are the material density and the total inelastic differential cross-section for electronic excitations, respectively. In Hartree atomic units (energies are in Hartree1 EH ¼ 27:2116 eV and lengths are in Bohr Bethe’s expression is radii: 1 a.u. ¼ 0.529177 A), ( Ip Z 4pZ12 X qmax dq SðvÞ ¼ Fn 0ðqÞ q v2 En0 =v n6¼0 ) Z qmax Z 1 dF ðq; EÞ dq dE þ : ð2Þ dE q E=v Ip The two contributions are from the discrete and continuum states, respectively, En0 is the usual energy difference between the ground and nth state, and Ip , qmax are the system ionization potential and the maximum momentum transfer [determined from the collision kinematics; see [7], Eqs. (17) and (18), respectively]. The upper limit of Ip on the sum is shorthand to indicate inclusion of terms up to a state energy difference of that amount. A general review of molecular stopping crosscalculations from oscillator strength distributions is in [4]. The exact result, Eq. (2) is relatively simple. In view of that simplicity and early calculations of atomic He generalized oscillator strengths (hereafter GOS),[8–10], there might seem to be little need for the continued development and use of various approximate schemes (examples may be found in [11,12] and references therein) for the calculation of SðvÞ in the first Born approximation. Such optimism would be tempered a bit by the need to include all excited states, both discrete and continuum, in Eq. (2). However, the issue was addressed in [7] at least for He. Those authors used a hydrogenic large-q GOS distribution to approximate the large-q contributions from high-lying states which are represented poorly or not at all by the finite basis set.
Another concern would be whether the ab-initio propagator techniques and large-q approximations used in [7] (or alternatively, the explicitly correlated wave functions used in [8–10]) could be utilized to treat the more complicated systems (solids, thin films, polymers, etc.) typically treated with more approximate methods. The issue is simple to state: are modern ab initio techniques for calculating stopping feasible for many-electron systems? This question is more-or-less unaddressed in the literature even for atoms beyond He, let alone more complicated systems. A simple, relevant test therefore is the modest step of treating the Be atom by direct application of the techniques of [7]. Note that though the Li atom might seem the logical next step, it is in fact rather less interesting. As an alkali its dominant high energy excitations will be hydrogenic. That fact makes Li well-suited to the methodology of [7] a priori, hence does not provide as serious a test as does Be. A brief summary of the methodology in question helps both in framing the issue and in defining notation. For transitions between initial and final states j0i and jni with energy difference En0 ¼ En E0 , the GOS is written as Fn0 ¼
2En0 2 jh0jeiqr jni j q2
with the second-quantized position operator g X ri;j ayi aj : r¼
ð3Þ
ð4Þ
i;j¼1
Here g is the total number of basis functions, and the matrix representation of the position operator in that basis is rij ¼
Z2 X hijr‘ jji:
ð5Þ
‘¼1
For the choice of operators A ¼ eiqr ; B ¼ eiqr
ð6Þ
the poles and residues, respectively, of the polarization propagator X h0jAjnihnjBj0i h0jBjnihnjAj0i hhA; BiiE ¼ ð7Þ E En þ E 0 E En þ E0 n6¼0 (with the sum over all excited states) correspond to the excitation energies and generalized transition moments necessary to calculate the GOS distribution (GOSD). All the relevant details can be found in [13]. Actual implementation for this study (except for the choice of basis set of course; see below) was the same as in [7], i.e., the calculations were done using the MUNICH [14] program package and within the random phase approximation (RPA).
J.A. Nobel et al. / Chemical Physics 309 (2005) 89–94
4. Results and discussion As customary in modern molecular calculations, a contracted Gaussian–type orbital (CGTO) basis was used. It was constructed starting with the CGTO (15s,9p,5d) ! [9s,9p,5d] basis set of Graham et al. [15] With that basis, those authors found that a full-configuration-interaction treatment gave a very good description of the low-lying excitations of the Be atom. Since higher-lying excitations are relevant to our problem, we added four primitive p-functions, one primitive d-function, and a new f-function sub-manifold of three primitives. From these, the basis was contracted as (15s13p6d3f) ! [9s13p6d3f]. Exponents and contraction coefficients are provided in Table 1. The additional diffuse orbitals were chosen in an attempt to improve fulfillment of the Bethe sum rule; see discussion below. The total closed-shell HF energy from both the original and augmented bases was )14.572608 EH compared to )14.572543 EH quoted in [15] and a HF limit Table 1 Basis exponents and contraction coefficients ‘-type exponent coefficients s
s
s s
p
p
d
f
3630.0 532.31 117.80 32.660 10.480 3.6680 3.6680 1.3540 0.38900 0.15020 0.05241 0.02000 0.00800 0.00300 0.00100 154.000 31.2234 6.71000 1.44200 0.41030 0.13970 0.04922 0.01600 0.00700 0.00300 0.00100 0.00050 0.00010 0.65000 0.20000 0.07000 0.01500 0.00320 0.00060 0.65000 0.07000 0.00320
0.000839 0.006735 0.035726 0.138635 0.385399 0.547688 0.213406 0.814692 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
91
[16] of )14.573023 EH . For a complete basis set, the Thomas–Reiche–Kuhn (TRK) sum rule for the dipole oscillator strengths fn0 X fn0 ¼ Z2 ð8Þ n6¼0
is exact in the RPA. Since this is the small q limit of the GOSD, the calculated value of the TRK sum is a relevant indicator of the quality of a finite basis set. A detailed discussion of the significance of the TRK sum rule in evaluating RPA stopping calculations is given by Cabrera-Trujillo et al. [4]. The essential point is that we obtain 4.0077865 (3.9923177) for atomic Be in the length (velocity) formulation, so that the augmented basis set appears to be quite reasonable in quality, for this particular sum rule. A basis of very high quality by usual quantum chemistry standards is, nevertheless, severely challenged in the calculation of the stopping cross-section. Details are in [4,7]. Briefly, any finite basis has, explicitly or implicitly, an upper limit on the atomic orbital angular momentum which it describes. Whatever that limit, a sufficiently large but finite value of q always can be found such that the expansion of Eq. (6) in powers of q has significant multipoles of higher angular momentum than that limit. The problem manifests itself quite strongly in the failure of the calculated Bethe sum rule (BSR) to behave as it should. The exact result [5] is X Fn0 ðqÞ ¼ Z2 ð9Þ n6¼0
for all values of q. The calculations on He [7] showed, however, that the calculated value of the sum in Eq. (9) is acceptably close to Z2 only for q 6 3 a.u. qL . For larger q the calculated sum falls rapidly to zero, a dramatic sign of the finite basis set having been overrun. The approximation used in [7] to eliminate this largeq misbehavior was to append an analytic hydrogenic GOSD to the RPA-propagator GOSD, with a suitably chosen value of q for the cut-over. Thus, Holm-Mortensen et al. [7] approximated the BSR (9) as Z 1 EL X dF ðq; EÞ Fn0 ðqÞ þ dE Z2 ð10Þ dE hyd EL n6¼¼0 with EL q2L =2 adjusted to make the approximation as good as possible. The hydrogenic form for dF j chosen dE hyd in [7] was from Ketkar and Bonham [17]. Fig. 2 of [7] shows that, for atomic He, the value EL ¼ 137 eV ) qL 3:2 a.u. provides a calculated GOSD which satisfies Eq. (10) to within a few percent for all q. Holm-Mortensen et al. [7] argued that the augmentation of the GOSD with a high-energy hydrogenic part is a reasonable procedure because the critical quantity Eq. (6) is a one-electron operator and the
92
J.A. Nobel et al. / Chemical Physics 309 (2005) 89–94
hydrogenic GOSD for large q therefore is physically sensible. The appending of the hydrogenic GOSD is equivalent logically (though not operationally) with selecting an approximate complement to the finite Gaussian basis and then restricting the propagator calculation solely to hydrogenic excitations within that complementary space. Here we test whether that argument is in any sense general, since in the Be atom there are manyelectron correlations not found in the He atom studied in [7]. The merit of choosing Be rather than Li should be clear, but perhaps it is worth adding that atomic Be is notorious (among quantum chemists) as a test in its own right because of the near-degeneracy of its ground state: so-called static correlation. For the test calculations the augmented BSR Eq. (10) was generalized to read Z 1( EL X dF ðq; EÞ Z2 Fn0 ðqÞ þ hðE E2s Þ dE hyd;2s EL n6¼¼0 ) dF ðq; EÞ þ hðE E1s Þ dE; ð11Þ dE hyd;1s where the Heaviside step functions turn on the correspondingly parameterized hydrogenic GOSDs for energies above some suitably selected excitation thresholds Enl . The most straightforward test case, in the sense of building in the least additional arbitrariness, is the one actually used, namely Enl ¼ EL , so that the entire hydrogenic augmentation came into play at once. The parameters of the Ketkar–Bonham model were chosen (using the same notation as in [7]) to be best-available atomic values; see Table 2. As before, EL was deter-
Table 2 Parameters used for the Ketkar-Bonham hydrogen-like GOSD Parameter (units) a
Ip (eV) n (a.u., inverse length)b k (a.u., charge)b a b
1s value
2s value
123.33 3.6848 0.956
9.16 3.6848 0.956
Ref. [23]. Ref. [24].
mined empirically to make Eq. (11) as close to an equality as possible. The result for the best value we could obtain EL ¼ 193:2 eV is in Fig. 1 (which is to be compared with Fig. 2 of [7]). Even the best value for EL results in a calculated BSR which falls (at q 2:45 a.u.) as much as 45% below the proper constant, Z2 . We were unable to find a combination of basis and EL value which does any better than this clearly unsatisfactory outcome. Clearly the appending procedure is unsuccessful for the Be atom, hence the procedure is not general. A brief calibration of the quality of the underlying propagator calculation may serve to reinforce the severity of this finding. Consider the excitation energies for the low-lying singlet transitions. Table 3 shows the present results, the RPA results from a smaller (9s9p5d) basis [18], and the experimental results.[19] The present results extend the earlier one in only one crucial way, namely the transition to the 4f state which could not be treated in the basis of [18], which had no f-type functions. Transitions to this state affect the GOSD and BSR only for q 6¼ 0. However, this is the domain in which improvement is needed. The differences between RPA and experimental excitation energies listed in Table 3
4
3.8
3.6
Calculated Z2 (q)
3.4
3.2
3
2.8
2.6
2.4
2.2 0
2
4
6
8
10
12
q (au)
Fig. 1. Calculated Z2 ðqÞ from Eq. (11).
14
16
18
20
J.A. Nobel et al. / Chemical Physics 309 (2005) 89–94 Table 3 Low-lying singlet excitation energies (eV) for the Be atom: comparison of the current calculation with a previous, small-basis RPA calculation [18] and with experiment [19] Final state
This calc.
Smaller basis RPA
Expt.
2s2p; 1 P0 2s3s; 1 S 2s3p; 1 P0 2s3d; 1 D 2s4s; 1 S 2s4p; 1 P0 2s4d; 1 D 2s4f; 1 F 2s5s; 1 S
4.80 6.12 6.72 6.94 7.26 7.48 7.59 7.68 7.72
4.80 6.12 6.72 6.94 7.26 7.48 7.59 – 7.72
5.28 6.78 7.46 7.99 8.09 8.31 8.53 Forbidden 8.60
For the forbidden 2s4f; 1 F our calculation gives zero dipole oscillator strength.
could be eliminated by use of higher-order polarization propagator approximations.[18] Such improvements in excitation energies would be of little consequence for the issue addressed here, namely the behavior of the calculated GOSD (which is determined primarily by the limitations of the finite basis set). It might be argued that a calculated GOSD that failed the BSR still could give a decent stopping cross-section. The BSR after all is not an integrated quantity whereas SðvÞ is. A qualitative argument shows that this expectation cannot be true in general. First, the polarization propagator produces a spectrum corresponding to the bound transitions of the system plus a set of pseudo-transitions to states that are the best representation of continuum states possible in a L2 basis. Within a finite L2 basis therefore, the stopping cross-section Eq. (2) has the approximate form nmax Z qmax 4pZ12 X dq ð12Þ SðvÞ ¼ Fn0 ðqÞ : q v2 n6¼0 En0 =v The limit nmax is determined by the basis set size. If it were the case that the calculated generalized oscillator strength sum were to saturate at or before nmax and fulfill the BSR, then that sum could be moved outside the integral (essentially as happens in the dipole approximation) and the result would be a modified mean excitation and Bethe logarithm. Failure to fulfill the BSR, on the other hand, introduces a spurious q-dependence that cannot be pulled out of the integral. We can model that dependence roughly (see Fig. 1) as nmax X
Fn0 ðqÞ ¼ Z2 ;
q 6 ðq0 D=2Þ;
n6¼0
4 2 ðq q Þ 1 ; ¼ Z2 1 þ b 0 D2 ðq0 D=2Þ 6 q 6 ðq0 þ D=2Þ; ¼ Z2 ; q P ðq0 þ D=2Þ:
ð13Þ
93
Note that the location q0 , width D, and depth b of the deficiency region in the calculated BSR all are velocitydependent, depending on the details of the basis. Insertion of this form into Eq. (12) and integration gives terms quadratic, linear, and logarithmic in q. Since the kinematic limits of q have velocity dependence of v1 and v, the net effect of the BSR violation is to introduce spurious v-dependence into SðvÞ that can range from v2 to v2 in addition to shifting the logarithmic term. Therefore the gross qualitative failure to fulfill the BSR such as we find for atomic Be is certain to have highly basis-set-specific effects on SðvÞ. In view of that specificity, we did not bother to compute SðvÞ for the particular basis used here since nothing general would have been learned.
5. Conclusions Though the pseudo-excitations from the polarization propagator treated in a finite L2 basis have no physical significance [20], the spectrum calculated with a highquality basis does satisfy the basic sum rule, namely TRK. There is, however, a severe finite basis problem for excitations at high momentum transfer. Striking evidence is in the failure of the calculated BSR, Eq. (9), to be independent of the value of momentum transferred. Augmentation of the calculated GOSD with simple hydrogenic distributions above some critical q, a successful technique for He, fails completely to resolve the problem for Be. The challenge involved in resolving the problem is easy to understand. The finite basis produces a kind of trench in the Bethe surface defined by the GOSD. The trench is roughly parallel with the E axis, hence over some finite range of q the BSR is not satisfied. A successful augmentation procedure would have to fill that trench so as to satisfy the BSR without building in spurious traits in the augmented GOSD. It is this challenge which suggests that attempts to tune various parameters, e.g., the Enl , in the model tested here would not be a good idea. Such tuning could well lead to a wholly unphysical augmented Bethe surface which nonetheless satisfied the BSR. In particular, the fulfillment of Eq. (9) is no guarantee of a correct GOSD, especially if that fulfillment is achieved by grafting on hydrogenic behavior for relatively low Enl . The calculations presented here are a state-of-the-art treatment of the GOSD. Since it is fundamental to firstprinciples calculation of electronic stopping in the first Born approximation, these calculations present a sobering challenge. Occasionally it has been suggested to the authors that approximate treatments of electronic stopping (such as those studied by us and others) no longer merit attention. Rather, the argument goes, one should proceed straightway with a proper GOSD
94
J.A. Nobel et al. / Chemical Physics 309 (2005) 89–94
calculation. The present calculation shows otherwise. If the state-of-the-art cannot treat a free Be atom adequately, it is clearly unrealistic to expect predictive GOSD calculations of electronic stopping for polymers, thin films, etc. In particular, the current status of timedependent density functional theory calculations for extended systems (see references in [21]) or dielectric treatments (RPA plus DFT orbitals and energies; see Mathar et al. [22]) makes clear that it is quite certain that even the level of refinement and generality reached in the present work is unattainable for those systems. Continued development and use of approximate stopping theories for systems more complicated than isolated atoms therefore seems not only warranted but inescapable. Acknowledgements We thank Stephan P.A. Sauer for useful conversations. J.R.S. and S.B.T. were supported in part by the US Army Research Office under Contract DAA L03-91G-0119. S.B.T. also was supported in part by US NSF Grant DMR-0218957. J.R.S. also was supported in part by US NSF Grant INT-9016299. J.O. was supported in part by the Danish Natural Science Research Council (Grant No. 11-0924). References [1] V. Nasluzov, N. R€ osch, unpublished. [2] I. Wilhelmy, N. R€ osch, Chem. Phys. 185 (1994) 317. [3] H. Bethe, Ann. Phys. (Leipzig) 5 (1930) 325.
[4] R. Cabrera-Trujillo, J.R. Sabin, J. Oddershede, Adv. Quantum Chem. 46 (2004) 121. [5] M. Inokuti, Rev. Mod. Phys. 43 (1971) 297. [6] M. Inokuti, Y. Itikawa, J.E. Turner, Rev. Mod. Phys. 50 (1978) 23. [7] E. Holm-Mortensen, J. Oddershede, J.R. Sabin, Nucl. Instrum. Meth. B 69 (1992) 24. [8] Y.K. Kim, M. Inokuti, Phys. Rev. 175 (1968) 176. [9] Y.K. Kim, M. Inokuti, Phys. Rev. 181 (1969) 205. [10] Y.K. Kim, M. Inokuti, Phys. Rev. 184 (1969) 38. [11] J. Wang, R.J. Mathar, S.B. Trickey, J.R. Sabin, J. Phys.: Condens. Matter 11 (1999) 3973. [12] D.E. Meltzer, J.R. Sabin, S.B. Trickey, J.Z. Wu, Nucl. Instrum. Meth. B 82 (1993) 493. [13] J. Oddershede, Adv. Chem. Phys. 69 (1987) 201. [14] G.H.F. Diercksen, W.P. Kraemer, Reference Manual to the MUNICH System of Programs, Max Planck Inst. f€ ur Astrophysik, Special Report, 1981. [15] R.L. Graham, D.L. Yeater, J. Olsen, P. Jorgensen, R. Harrison, S. Zarrabian, R.J. Bartlett, J. Chem. Phys. 85 (1986) 6544. [16] C. Froese-Fischer, The Hartree–Fock Method for Atoms, Wiley, New York, 1977. [17] N. Ketkar, R.A. Bonham, J. Chem. Phys. 84 (1986) 6091. [18] J. Geertsen, J. Oddershede, G.E. Scuseria, Int. J. Quantum Chem. S21 (1987) 475. [19] S. Bashkin, J.A. Stoner, Atomic Energy Levels Grotrian Diagrams, North-Holland, Amsterdam, 1975. [20] J. Geertsen, J. Oddershede, J.R. Sabin, Phys. Rev. A 34 (1986) 1104. [21] H. Appel, E.K.U. Gross, K. Burke, Phys. Rev. Lett. 90 (2003) 043005; M.E. Casida, D.R. Salahub, J. Chem. Phys. 113 (2000) 8918. [22] R.J. Mathar, S.B. Trickey, J.R. Sabin, Adv. Quantum Chem. 45 (2004) 277; Nucl. Instrum. Meth. B 155 (1999) 249. [23] D.A. Shirley, R.L. Martin, S.P. Kowalczyk, F. McFeely, L. Ley, Phys. Rev. B 15 (1977) 544. [24] E. Clementi, D.L. Raimundi, J. Chem. Phys. 38 (1963) 2686.