Effective Core Potentials for Calculations of Continuum Spectra of Molecules Using the Molecular R‑Matrix Method

TL;DR

This paper introduces a method to use effective core potentials in molecular calculations for plasma modeling, enabling accurate simulations of electron collisions and photoionization.

Contribution

The paper presents a novel implementation of ECPs in the UKRmol+ suite for continuum spectral calculations.

Findings

ECP integrals over B-spline orbitals were derived and implemented for molecular calculations.

Sample calculations for electron collisions and photoionization of various molecules were successfully performed.

The method is applicable to large targets and high-energy electrons in plasma modeling.

Abstract

Implementation of effective core potentials (ECPs) into
the molecular
scattering suite UKRmol+ is presented together with a set of calculations
for a range of targets relevant for plasma modeling. Continuum description
in scattering and photoionization calculations for large targets or
high-energy electrons often requires the use of numerical continuum
functions and the associated molecular integrals. We derive expressions
for ECP integrals over B-spline-type orbitals using their momentum-space
representation and describe their implementation. Sample calculations
are presented for electron collision with ethylene (C2H4), bromine (Br2), silicon tetrabromide (SiBr4), and tungsten hydride (WH), as well as photoionization of
methyl iodide (CH3I).

Figures

Elastic scattering cross-section for C2H4.

Electronic excitation cross-sections for C2H4. Solid lines are the cross-sections calculated using an ECP, dashed lines, are the cross-sections calculated using the all-electron 6–31G basis set.

Eigenphases for C2H4. Solid lines are the results using the 2-electron ECP for each C atom, dashed lines are the results from the all-electron calculation.

Elastic scattering and momentum transfer cross-sections for Br2.

Electronic excitation cross-sections for Br2 Π states (left) and Σ and Δ states (right). Dashed lines are the cross-sections obtained using the all-electron cc-pVTZ basis set and solid lines are the cross-sections obtained using an ECP.

Eigenphases for Br2. Solid lines are the results using the 28-electron ECP for each Br atom and dashed lines are the results from the all-electron calculation.

Elastic scattering and momentum transfer cross-sections for SiBr4, compared to the measured grand total cross-section given in ref

Vibrational excitation cross-sections for SiBr4.

Total ionization cross-section for SiBr4.

Elastic scattering cross-section for WH using different inner-region sphere radii. The calculation using a sphere radius of 24 Bohr used B-spline continuum functions (see text for details).

Electronic excitation cross-sections for the lowest sextet and quartet excited states of WH.

Total ionization cross-section for WH.

Photoionization of CH3I into the ground state of the ion. Left: photoionization cross-section. Right: photoelectron angular distributions. Calculations with ECPs are compared to all-electron calculations and experimental measurements.

Tables

1. Equilibrium Geometry of C2H4 in the Centre-of-Mass Frame

atom	x (Å)	y (Å)	z (Å)
C	0	0	0.657
C	0	0	–0.657
H	0	0.9145	1.2208
H	0	–0.9145	1.2208
H	0	0.9145	–1.2208
H	0	–0.9145	–1.2208

2. Excitation Energies for Electronic States of C2H4 Calculated with and without ECPs

	excitation energy (eV)
state	all-electron	ECP
1 3 B 1u	4.805	4.790
1 3 B 1g	9.402	9.338
1 1 B 1g	9.767	9.705
1 1 B 1u	10.15	9.93

3. Comparison of Vertical Excitation Energies for Br2 Calculated Using All-Electron Basis Sets with Those Obtained Using a 28-Electron ECP for Each Br

	excitation energy (eV)
state	cc-pVTZ	ECP	ref
1 3Π u	2.67	2.54	2.64
1 1Π u	3.61	3.50	3.56
1 3Π g	5.26	5.14	5.22
1 Σ3g−	5.59	5.38	5.71
1 1Δ g	5.95	5.75	6.07
1 1Π g	6.17	6.07	6.10
2 Σ1g+	6.26	6.07	6.39
1 Σ1u−	6.82	6.62	6.91
1 3Δ u	6.92	6.72	6.99
1 Σ3u+	7.00	6.80	7.10

4. Equilibrium Geometry of SiBr4 in the Centre-of-mass Frame

atom	x (Å)	y (Å)	z (Å)
Si	0	0	0
Br	0	–1.8087	1.2790
Br	0	1.8087	1.2790
Br	1.8087	0	–1.2790
Br	–1.8087	0	–1.2790

5. Comparison of Vibrational Frequencies for SiBr4 Calculated Using QEC with Data from NIST; Γ Is the Symmetry Group of the Mode in the T d Point Group

			frequency (cm–1)
label	mode	Γ	NIST	QEC
ν1	symmetric stretch	A 1	249	253.49
ν2	bending deformation	E	90	86.41
ν3	asymmetric stretching	T 2	487	505.78
ν4	bend	T 2	137	135.91

6. Excitation Energies for Electronic States of WH Calculated Using QEC Compared to Those of Ref

	excitation energy (eV)
state	this work	ref
1 4Δ	0.855	0.749
1 4Π	1.046	0.954
1 4Σ+	1.875	1.829
1 2Σ+	2.173	-
1 2Δ	2.207	1.947
2 2Δ	2.221	2.070
1 6Π	2.233	1.515
1 2Π	2.243	-
1 2Σ–	2.315	-
2 4Π	2.321	2.555
2 4Δ	2.402	2.244
2 2Π	2.409	-
3 2Π	2.452	-
3 4Δ	2.457	-
4 2Π	2.464	-
3 2Δ	2.464	2.072
1 6Δ	2.495	2.099
1 4Σ–	2.514	-
3 4Π	2.557	-
2 4Σ+	2.836	-
4 4Π	3.110	-
2 6Σ+	3.354	2.772
5 4Π	3.402	-
2 4Σ–	3.424	-
6 4Π	3.500	-
7 4Π	3.599	-
2 6Π	3.759	4.223

Full text

Introduction

1

Electron dynamics in atomic and molecular systems is challenging not only because of the correlated motion of electrons but also due to the possible presence of relativistic effects in systems containing heavy atoms. In addition, a practical complication arises in the calculation of the related molecular two-electron integrals, which generally scale with the fifth power of the number of the molecular orbitals.

One way of alleviating these challenges is to take advantage of the orbital occupation model to reduce the number of electrons explicitly included in the molecular Hamiltonian, replacing them with effective potentials, representing an average effect of those electrons on the remaining ones. ?,? Naturally, the core- and inner-valence electrons can be satisfactorily represented by the mean field when dynamics occurring in the valence space are considered. Methods of this kind have been developed for bound state quantum chemistry calculations starting in the 1960s. Since then, they have become ubiquitous and a standard part of bound state calculations. Two major strands of effective potentials have been developed: pseudopotentials and model potentials. The latter aims to preserve the correct nodal structure of the valence orbitals, while the former leads to so-called pseudovalence orbitals whose representation may benefit from using specialized Gaussian-type atomic bases.? For that reason, the choice of the pseudopotential often comes with a recommended optimal Gaussian-type valence basis set. The pseudopotentials, more commonly called effective core potentials (ECPs), have become the dominant effective representations of the core. Since the late 1970s, ECPs have been extended to incorporate scalar and spin–orbit relativistic effects of the core on the valence shells. The history and properties of the various forms of the effective and model potentials used to date have been described in comprehensive reviews. ?,?

The use of ECPs in calculations of continuum spectra of molecules is much less common. To the best of our knowledge, the working implementations of those are in the Schwinger multichannel codes of the Brazilian group ?,? and in the complex Kohn codes. ?,? While ab initio multi-electron approaches for continuum states are able to represent electron correlation at a high level, their extension to the inclusion of relativistic effects remains a challenge. Relativistic effects have been included in all-electron atomic scattering and photoionization calculations starting with the R-matrix implementation in the early 1980s.? Fully relativistic atomic codes based on the Dirac equation also exist; see e.g. ?,? More recently, relativistic effects have been included in time-dependent simulations of atoms in external fields,? but their extension to molecular calculations remains an open and challenging problem and is not the focus of this work.

ECPs thus form a middle ground enabling the study of relativistic effects of the core electrons on the valence ones by making relatively minor modifications of the existing scattering codes. At the same time, the use of ECPs leads to a valuable saving of computational time for the two-electron integral problem since the calculations of continuum spectra necessarily include large bases of continuum functions, which makes the integral calculations relatively more demanding compared to the bound-state problems of quantum chemistry.

In this work, we describe the implementation of ECPs in the molecular R-matrix codes UKRmol+? and in the related interface QEC? and the first applications to electron scattering from molecules containing heavy atoms. Section 2 describes the implementation in detail. In Section 3, we provide examples and discuss the performance of the calculations, including ECPs. Finally, in Section 4, we discuss the limitations of the present approach and the future prospects for its extension.

Methods

2

In atomic units, the molecular fixed-nuclei Hamiltonian for n _ v _ valence electrons with n c core electrons represented by ECPs has the following form

$$[eqn]$$ $$[eqn]$$

where r _ i _ represent the spatial coordinates of the valence electrons, R _ i _ represents the positions of the atomic nuclei, and Z _ i _ represents their charges. The one-electron contribution h(i) consists of the electron’s kinetic energy and the nuclear attraction energy between the electron and each nucleus whose bare charge is effectively screened by the n c(j) core electrons represented by the ECP centered on that atom.

The ECP comprises a local term and a semilocal term, which is further split between scalar and spin–orbit components

$$[eqn]$$

where the semilocal terms contain the projectors on the real spherical harmonics $[eqn]$ . Here, r and the angular-momentum projectors are taken with respect to the atom where the ECP is centered. The radial potentials are all parametrized by a linear combination of Gaussians multiplied by powers of the radial distance

$$[eqn]$$

For reasons of integrability, the coefficients m _ j _ must not be smaller than zero.

Implementing ECPs in UKRMol+

2.1

ECPs can be used within the standard workflow of the R-matrix method.? The implementation of ECPs in UKRmol+ is therefore reduced to evaluating the matrix elements of the ECP potential (3) for the basis of one-electron atom- and center-of-mass-centered functions representing the bound and the continuum electrons, respectively. The ECP integrals have been implemented in the GBTOlib library,? which calculates all molecular integrals required by UKRmol+. In GBTOlib, the continuum can be represented using an arbitrary combination of GTOs and B-spline-type orbitals (BTOs) centered at the center-of-mass,? while the bound electrons are represented using the standard atom-centered GTOs of quantum chemistry. Orthogonalization of the continuum orbitals against the set of orbitals representing the bound electrons then leads to the mixing of all one-electron basis functions in the expansion of the molecular continuum orbitals. This, in turn, requires calculation of the ECP matrix elements for all combinations of the one-electron basis functions.

Evaluation of ECP matrix elements between GTOs can be performed using known analytic expressions? or using highly optimized numerical routines employed for the calculation of 2-electron hybrid integrals.? Therefore, the pure GTO integrals will not be discussed here further. However, besides integrals over GTOs, we require also integrals involving the numerical B-spline orbitals if those are used to represent the continuum. These include 1-, 2-, and 3-center integrals between atom-centered GTOs and BTOs centered on the center-of-mass. For that reason, the formulas for the matrix elements will be presented first in the general form suitable for a numerical quadrature, regardless of the orbital type used. The general expressions are then specialized to particular cases.

For completeness, we mention that both the semilocal and local terms involving BTOs and ECPs centered on the center-of-mass are calculated directly in real space trivially by means of radial quadratures.

One-Electron Basis Functions

2.1.1

The general form of the one-electron basis functions is

$$[eqn]$$

where R(r) is the radial part of the orbital.

For a contracted GTOs, the expression specializes to?

$$[eqn]$$ $$[eqn]$$

where N G normalizes the orbital to unit charge density, S _ lm _ is the real solid harmonic, c _ j _ are the contraction coefficients, and α_ j _ are the contraction exponents.

In the case of BTOs, the radial part of the orbital is the numerical B-spline function. The resulting BTO as employed in GBTOlib is?

$$[eqn]$$ $$[eqn]$$

where N _ B _ i _ _ is again the normalization factor and B _ i _(r) is the ith B-spline function drawn from a basis of B-splines? covering a selected radial range. B-splines are piecewise polynomials with compact support, which allow the construction of a highly accurate and linearly independent basis for the oscillating wave function of the unbound particle.

For clarity, the normalization factors for the GTOs and BTOs will be omitted from the following expressions but are assumed to be included as overall multiplication factors.

Translation of the BTOs

2.1.2

The atomic character of the ECPs and the angular-momentum projectors employed in the semilocal terms requires that the integration over the semilocal component be performed most conveniently with respect to the atomic nucleus where the ECP is centered. This is also the approach used in the analytic evaluation of ECP integrals over GTOs, and we use it for the evaluation of the mixed and pure BTO integrals too. When computing ECP integrals, we encounter 1-, 2-, and 3-center integrals.

Computing the GTO-only integrals makes use of the analytic form of the translation of the GTO with respect to an arbitrary center. ECP integrals involving numerical BTOs must be handled differently. One option is to compute the integrals in real space by rotating the coordinate system along the line connecting the center-of-mass with the atomic nucleus, allowing the analytic computation of the integral over the azimuthal angle and the subsequent numerical integration over the remaining two coordinates. However, the finite support and the piecewise character of the radial B-splines complicate the angular integrals. While this is certainly a solvable problem, we have implemented an alternative approach, making use of the translation property of Fourier transforms to perform arbitrary translations easily in momentum space.

Fourier transform of the BTO is simple

$$[eqn]$$ $$[eqn]$$

and follows trivially from the partial wave expansion of the plane-wave in the basis of real spherical harmonics

$$[eqn]$$

We note that the Fourier transformation of the BTO preserves the angular momentum character of the orbital.

The factor τ_ i,l _(k) contains the Riccati-Bessel function $[eqn]$ and is analogical to the spherical Bessel transform investigated earlier by Talman in the context of multicenter integrals over numerical orbitals. ?−? ? An efficient computation of this factor is crucial, in particular, for B-splines localized further away from the origin, where the integrand becomes highly oscillatory. The methods of Talman make use of logarithmic radial grids, which are applicable to orbitals with an infinite range, such as Slater orbitals, but not to B-splines extending over the finite radial range [a _ i _, b _ i ]. Instead, to calculate τ i,l _(k), we implemented Levin quadrature,? which is tailored to integrands of this type for arbitrary values of k.

The real-space representation of the BTO is obtained from the inverse Fourier transform

$$[eqn]$$

It follows that the translation of the BTO to another center A is performed by applying a phase shift to the momentum space coefficients

$$[eqn]$$

An equivalent relation can be obtained for any function equipped with a Fourier transform. In particular, we will use it for the BTO-only integrals to translate the local component of the ECP to the center-of-mass.

To perform the ECP integrals over BTOs, we need to evaluate, at a given distance r _ A _ from the atom, a projection of the BTO onto spherical harmonics centered on the atom, A. This projection is straightforwardly derived using the Fourier representation (13) and the formulas listed above

$$[eqn]$$

The partial wave expansion and expression (10) were used to simplify the derivation. Here, the symbol $[eqn]$ stands for the Gaunt coefficient for real spherical harmonics.?

Semilocal Term

2.2

The general expression for the matrix element of the semilocal term can be written as a radial integral with respect to the atomic center A

$$[eqn]$$ $$[eqn]$$

where the two functions ϕ_ i _ and ϕ_ j _ stand for any type of basis function. Computing the integral involving the semilocal potential requires calculation of the projections of the basis functions on the atom-centered spherical harmonics. When one of the functions is a BTO, we employ Formula (16). For hybrid integrals, the GTO partial waves are obtained from a plane-wave expansion provided by GTOlib. Here, we provide the working formulas for the particular classes of integrals involving BTOs.

BTO–BTO Class

2.2.1

The expression for this integral is arrived at by employing the angular projections for the BTOs in the momentum-space representation (16) and integrating the result over the radial coordinate r _ A _ centered on the ECP

$$[eqn]$$ $$[eqn]$$

The momentum-space integrals are currently computed by using sufficiently large momentum-space quadratures. The BTOs located farther away from the center-of-mass lead to a highly oscillatory momentum-space representation τ_ i,l _(k), but in many cases, the heavy atoms containing ECPs are localized close to the center so that they do not overlap with those BTOs in real space, and the corresponding ECP integral is zero. Although dense and wide momentum-space grids are required, the computation is still manageable and provides orders of magnitude computational savings compared with the all-electron calculation. Nevertheless, the development of more optimal quadratures suitable for these oscillatory integrals is highly desirable and will be the subject of follow-up work.

GTO–BTO Class

2.2.2

To evaluate this class, we need to combine the numerically evaluated angular projection of the GTO in real space with the angular projection of the BTO evaluated using the Fourier representation. The result is

$$[eqn]$$ $$[eqn]$$

Local Term

2.3

$$[eqn]$$ $$[eqn]$$

This type of integral can be formally recast in the form of a projection of the product of the two basis functions on the spherical harmonic X 00 followed by integration with respect to the ECP center

$$[eqn]$$

BTO–BTO Class

2.3.1

The pure BTO case can be computed fairly straightforwardly by representing the local potential in the momentum space, followed by its translation to the center-of-mass while keeping the BTOs in their real-space representation

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

GTO–BTO Class

2.3.2

In this case, one is tempted to combine the terms from the local potential (r ^ n ^ exp­[−γr ^2^]) with the atom-centered Gaussian to obtain overlap-type integral between BTOs and GTOs. Those integrals can be handled very efficiently using GBTOlib. However, the radial parts of the local potential may contain odd powers of the radial coordinate, preventing the polynomial-type translation of the prefactor to the center-of-mass. This makes the BTO-GTO class for general local potentials less amenable to analytic techniques. Instead, this type of integral can be handled by a direct 3D quadrature with respect to the center of the BTO and will be implemented in future releases of GBTOlib.

QEC Implementation

2.4

QEC? is an expert system that performs electron-molecule scattering calculations, using both the UKRMol+ codes? and the MOLPRO quantum chemistry package. ?,? Although QEC has mainly been applied to systems containing fairly light elements, ?−? ? it has made some use of ECPs in the past. Originally, ECPs were implemented in QEC to calculate electron-impact ionization cross-sections? using the binary encounter Bethe (BEB) method of Kim and Rudd.? This was straightforward since BEB calculations require input from MOLPRO only and are independent of the R-matrix calculation. The use of ECPs within the R-matrix method now enables QEC to calculate the full set of cross-sections for molecules containing heavy elements.

QEC uses the ECPs provided by the Stuttgart/Köln group.? These potentials are named ECPnXYZ, where n is the number of core electrons represented by the ECP, and the XYZ label indicates the level of theory used when constructing the pseudovalence orbitals. The X label indicates the reference system used to generate the ECP, X = S, if a single-valence-electron ion was used, and X = M, if a neutral atom was used. The YZ label denotes the method used, YZ = HF for the nonrelativistic Hartree–Fock method, YZ = WB for the quasi-relativistic Wood-Boring method, and YZ = DF for the relativistic Dirac–Fock method. The calculations reported in this work use semirelativistic ECPs labeled ECPnMWB.

A range of valence-electron basis sets can be used with these potentials. Associated basis sets named ECPnMWB? are provided in MOLPRO. In some cases, these are double-ζ bases, though for some elements, a range of triple and quadruple-ζ basis sets are available, as well as basis sets including diffuse functions. The Karlsruhe (named def2-XYZ) basis sets? are also compatible with the Stuttgart/Köln ECPs for some heavy elements.

Results and Discussion

3

In this section, we provide a set of electron scattering cross-sections for a number of molecular targets. In Section, we provide an initial proof-of-principle example for a light system, C_2_H_4_, to demonstrate that calculations that use ECPs for light elements can reproduce cross-sections obtained using all-electron basis sets. In Section, we provide a further proof-of-principle test case using the bromine molecule, for which calculations using ECPs may be benchmarked against an all-electron calculation. Electron scattering cross-sections are then presented in Sections and ? for the heavier SiBr_4_ and WH systems, which are beyond the reach of all-electron calculations.

The ECP capability has been implemented for photoionization calculations, too. In Section, we present one example for valence photoionization of the CH_3_I molecule.

C2H4

3.1

As an initial proof-of-principle test case, we calculate cross-sections for C_2_H_4_, a system containing only light elements for which cross-sections can be calculated using ECPs and valence basis sets as well as all-electron basis sets. Cross-sections calculated using these two methods ought to show a high degree of agreement given the absence of any heavy elements in C_2_H_4_. We apply the two-electron ECP, ECP2MWB, to each C atom, with the associated basis set for the four valence electrons. Since this basis set contains s and p orbitals, we choose the 6-31G basis set as the all-electron basis set, since it also contains only s and p orbitals for C. In both calculations, the 6-31G basis set is also applied to the H atoms. In both calculations, the unbound electron was treated using a GTO basis with maximum angular momentum l max = 4, and the R-matrix radius was set to 14 Bohr. The equilibrium geometry is given in Table and is optimized using the MOLPRO quantum chemistry package? using Hartree–Fock orbitals and the cc-pVTZ basis set. The R-matrix calculations are carried out using the irreducible representations of the D _2h _ point group and the close-coupling method, which uses a complete active space (CAS) description of the excited target states and CASSCF orbitals. The ground state configuration is $[eqn]$ . The calculation using the 6-31G basis set used an active space in which 4 electrons are frozen (1a _ g _ ^2^,1b _1u _ ^2^), and the remaining 12 electrons are distributed among 10 orbitals, (2 – 3a _ g _, 1b _3u _, 1 – 2b _2u _, 2 – 3b _1u _, 1b _2g _, 1 – 2b _3g _). The calculation using an ECP for each C atom is consistent with this setting, with the ECP handling a total of four electrons, and the active space is composed of the 12 remaining electrons distributed among the orbitals (1 – 2a _ g _, 1b _3u _, 1 – 2b _2u _, 1 – 2b _1u _, 1b _2g _, and 1 – 2b _3g _). Six additional virtual orbitals were included in both calculations.

1: Equilibrium Geometry of C2H4 in the Centre-of-Mass Frame

Figure shows the elastic scattering cross-section for C_2_H_4_, calculated with and without ECPs. As expected, the two calculations produce almost identical cross-sections at all energies, including the character of the resonance at around 3 eV.

Elastic scattering cross-section for C2H4.

Figure shows the electronic excitation cross-sections for C_2_H_4_. Four excited states are found to lie below the ionization potential of 10.5 eV. The vertical excitation energies of these states are given in Table. Once again, the two methods produce nearly identical cross-sections, with the main difference appearing to be in the vertical excitation energy of the ^1^ B _1u _ state, which differs by 0.2 eV in the two calculations. Overall, excitation of the lowest ^3^ B _1u _ state dominates, and both calculations capture the profile of this cross-section equally well.

Electronic excitation cross-sections for C2H4. Solid lines are the cross-sections calculated using an ECP, dashed lines, are the cross-sections calculated using the all-electron 6–31G basis set.

2: Excitation Energies for Electronic States of C2H4 Calculated with and without ECPs

In addition to the cross-section data, perhaps a more detailed and fundamental test of the similarity in the two approaches is made by comparing the calculated eigenphases in each case. Figure shows the eigenphases produced by the two calculations. Excellent agreement is observed between the all-electron and ECP-based calculations for all irreducible representations. Discontinuities in the derivatives of the eigenphases are expected to occur at the vertical excitation energies of the excited states. Both calculations show these discontinuities at almost identical energies, which confirms that both can capture electronic excitation processes to an almost identical level of accuracy.

Eigenphases for C2H4. Solid lines are the results using the 2-electron ECP for each C atom, dashed lines are the results from the all-electron calculation.

Br2

3.2

The Br_2_ molecule serves as a further proof-of-principle test case in which cross-sections can be calculated using ECPs and valence basis sets, as well as all-electron basis sets. We compare cross-sections for Br_2_ calculated using the all-electron cc-pVTZ basis set with those obtained using a 28-electron scalar-relativistic core potential (ECP28MWB) for each Br atom, with the remaining 7 valence electrons of each Br treated using the ECP28MWB_VTZ basis set. In both cases, the unbound electron was treated using a basis of GTOs with angular momentum l ≤ 4. The R-matrix radius was set to 14 Bohr. All cross-sections were calculated by using the close-coupling method. The calculations used the irreducible representations of the D 2h _ point group. Using these irreducible representations, the ground state configuration of Br_2 may be expressed as

$$[eqn]$$

The calculation using the cc-pVTZ basis set used a target active space in which 62 electrons were frozen

$$[eqn]$$

leaving 8 valence active electrons distributed among 8 active orbitals, $[eqn]$ . The active space settings for the calculations using ECPs mirror this calculation: 56 electrons are handled using the ECPs, leaving a residual target with 14-electron, 6 of which are kept frozen, occupying the lowest a _ g _ and b 1u _ orbitals of the valence-electron basis set, $[eqn]$ , and the remaining 8 electrons are allowed to occupy the active orbitals $[eqn]$ . Six additional virtual orbitals were also included in both calculations. The equilibrium geometry used for Br_2 assumed an internuclear separation of 2.2756 Å.

Figure shows the elastic scattering and momentum transfer cross-sections for Br_2_, calculated with and without the use of ECPs. A strong level of agreement is observed in both cross-sections, including the positions of some small resonance features. Small differences in the cross-sectional values are apparent at low energies. It is tempting to ascribe these differences to the semirelativistic nature of the ECP and valence basis set. However, it is not clear if this is the case, since a nonrelativistic ECP for Br, that would provide definitive proof of such an effect, is not available.

Elastic scattering and momentum transfer cross-sections for Br2.

Figure shows the calculated electronic excitation cross-sections for Br_2_. Again, very good agreement is observed between the two results, including the resonant features of the 1^3^Π_ u _ state above 4 eV, as well as a similar feature seen for excitation to the 1^1^Π_ u _ state. This also extends to the vertical excitation energies given in Table. The largest relative error in the excitation energies for the first few excited states was found to be 4%. The excitation energies are also in good agreement with recent nonrelativistic calculations.? The main difference found is that the ^1^Π_ g _ and $[eqn]$ states are almost degenerate in the calculation using ECPs, while these states are distinctly separated in the nonrelativistic all-electron calculations given here and in ref ?.

Electronic excitation cross-sections for Br2 Π states (left) and Σ and Δ states (right). Dashed lines are the cross-sections obtained using the all-electron cc-pVTZ basis set and solid lines are the cross-sections obtained using an ECP.

3: Comparison of Vertical Excitation Energies for Br2 Calculated Using All-Electron Basis Sets with Those Obtained Using a 28-Electron ECP for Each Br

Figure shows the calculated eigenphases for Br_2_. Excellent agreement is observed between the eigenphases obtained using the all-electron cc-pVTZ basis set (dashed lines) and those obtained using ECPs (solid lines). The good agreement seen here confirms that both calculations capture resonant features and electronic excitation processes at almost identical energies.

Eigenphases for Br2. Solid lines are the results using the 28-electron ECP for each Br atom and dashed lines are the results from the all-electron calculation.

SiBr4

3.3

Following our proof-of-principle calculations of Br_2_ cross-section data, we calculated cross-sections for silicon tetrabromide, a compound used as a precursor in the deposition of silicon nitride. Experimental data on SiBr_4_ are relatively scarce, though measurements of resonance energies and dissociative electron attachment ion yields are available. ?,? Previous calculations of SiBr_4_ cross-sections may be found in ref ? and ? , where cross-sections were calculated using the Schwinger multichannel method with pseudopotentials at both the static exchange (SE) and static exchange plus polarization (SEP) levels.

SiBr_4_ is a 154-electron system, and although all-electron basis sets could be used for both Si and Br, an all-electron calculation would be impractical on computational grounds. Both the calculation runtime and memory requirements of the R-matrix method scale strongly with the number of electrons in the target molecule. Our calculations used a 28-electron scalar-relativistic core potential (ECP28MWB) for each Br atom, and the remaining 7 valence electrons of each Br were treated using the ECP28MWB_VTZ basis set. A cc-pVTZ basis set was used for the Si atom. These settings reduce the number of electrons to be treated using basis sets from 154 to 42. A pure GTO continuum basis with a maximum angular momentum of l max = 4 was used to represent the unbound electron, and the R-matrix radius was set to 14 Bohr. The equilibrium geometry used for SiBr_4_ is given in Table. Each of the Si–Br bond lengths is 2.2152 Å. The tetrahedral SiBr_4_ molecule belongs to the T _ d _ point group in its equilibrium geometry, and our calculations are carried out using irreducible representations of the C _2v _ point group, the highest Abelian subgroup of T _ d _.

4: Equilibrium Geometry of SiBr4 in the Centre-of-mass Frame

Figure shows the calculated elastic scattering and momentum transfer cross-sections for SiBr_4_, calculated using the SEP method. Prominent resonances appear at around 1.2 and 6 eV, in good agreement with electron transmission measurements, ?,? as well as calculations using the Schwinger multichannel method.? A symmetry analysis of these resonances shows that they are of T 2 and E symmetry, as also shown in ref ?. In the SEP calculations, pseudoresonances begin to appear just below 8 eV due to the neglect of electronic excitation processes in the SEP method.

Elastic scattering and momentum transfer cross-sections for SiBr4, compared to the measured grand total cross-section given in ref

Table gives vibrational frequencies calculated by MOLPRO; these agree with values from NIST? to within 5% for all modes. Figure shows the vibrational excitation cross-sections for v = 0 → 1, calculated using the method of Kokoouline et al.? The dipole-allowed excitation of the T 2 asymmetric stretching mode is strongest at all energies, though significant contributions appear from all modes, including the symmetric stretch and bending deformation modes that are not dipole-allowed.

5: Comparison of Vibrational Frequencies for SiBr4 Calculated Using QEC with Data from NIST; Γ Is the Symmetry Group of the Mode in the T d Point Group

Vibrational excitation cross-sections for SiBr4.

Figure shows the total ionization cross-section for SiBr_4_, calculated using the Binary-Encounter Bethe method.? The ionization potential obtained using Koopman’s theorem was found to be 11.59 eV, in reasonable agreement with the measured value of 10.62 eV.?

Total ionization cross-section for SiBr4.

WH

3.4

Progressing to even heavier systems, a set of scattering cross-sections for WH was calculated. Tungsten (W) is a preferred plasma-facing material in experimental fusion plasma, including the ITER. Emission from the WH (in practice WD) electronic band ^6^Π → ^6^Σ^+^ have been observed in such plasma environments.? These emissions provide the easiest way to monitor W erosion and molecule formation, with electron collisions providing the main mechanism of populating the ^6^Π state. Electronic excitation cross-sections for this species are therefore important.

For WH (WD), the cc-pVTZ basis set was used for the H atom, and a 60-electron scalar-relativistic core potential (ECP60MWB) was used for the W atom, leaving the remaining 14 valence electrons of W to be treated with a basis set. In this case, the choice of basis set required careful consideration of the angular momentum components retained in the basis. The basis set that accompanies the ECP60MWB core potential for W contains g functions and could only be adequately contained with an inner region sphere of around 20 Bohr. This in turn affected the choice of the continuum basis set. Gaussian-type orbitals centered on the center-of-mass are suitably accurate for inner-region radii of up to 14 Bohr. For larger inner regions, B-spline-type orbitals (BTOs) are available in the R-matrix codes,? although this feature is not available in QEC.

In order to perform a calculation using the standard parameters of the available QEC package (a purely GTO continuum basis and inner region radius of up to 14 Bohr), alternative target basis sets were considered. In particular, the Karlsruhe family of basis sets was considered, as these can be used in conjunction with the Stuttgart ECPs for some elements.? The def2-TZVP basis set was chosen, as it contains up to and including f basis functions. Using this basis, calculations could be completed using an inner region radius of 14 Bohr. Given the numerical issues encountered when using the Stuttgart basis set, it was thought prudent to test the sensitivity of cross-sections obtained using the def2-TZVP basis set to the inner region radius.

An additional calculation was carried out using an inner region sphere of 24 Bohr. These calculations used BTOs of order 4 starting at a radius of 11 Bohr. One of the main computational demands of such calculations is the evaluation of mixed GTO/BTO integrals. These integrals were calculated using Gauss-Legendre quadrature with stepsize Δr = 0.25, in which the Coulomb potential was represented by Legendre expansion with a maximum angular momentum of 8.

The equilibrium bond length calculated by MOLPRO was found to be 1.8377 Å. This calculation used the cc-pVTZ basis set for H, while W was treated by using the 60-electron core potential ECP60WMB and the def2-TZVP basis set. All cross-sections were calculated using the close-coupling method since the first excitation threshold was found to be close to 1 eV. SEP calculations are likely to display pseudoresonances starting at energies close to this threshold, leaving a physically meaningful cross-section only at very low energies. The active space for the close-coupling calculations allowed 7 active electrons to occupy 10 orbitals. The active space consisted of 8 frozen electrons in the orbitals $[eqn]$ , and the 7 active electrons occupying the orbitals (3–5a 1, 2–3b 1, 2–3b 2, 1a 2). Two additional virtual orbitals were also included.

Figure shows the elastic scattering cross-section for WH, calculated using inner-region sphere sizes of 14 Bohr and a GTO continuum basis, as well as a calculation using an inner region radius of 24 Bohr and a mixed GTO/BTO continuum basis. In both cases, the continuum basis retained a maximum angular momentum of l = 4. Since WH is a highly polar molecule, the dipole Born correction was added to the cross-section to account for higher angular momenta. This correction increases the cross-section significantly at low energies. A series of prominent resonances are visible between energies of 2 and 4 eV, indicating the likely presence of excited target states in this energy range. As can be seen in Figure, very little sensitivity to the sphere size was observed in the cross-sections.

Elastic scattering cross-section for WH using different inner-region sphere radii. The calculation using a sphere radius of 24 Bohr used B-spline continuum functions (see text for details).

The vertical excitation energies calculated using QEC may be compared with those given in ref ?. As shown in Table, we find a significant number of additional excited states that are not listed in ref ?, particularly doublet and quartet states. The level of agreement between our calculations and those given in ref ? is generally very good, with perhaps the exception of the 1^6^Π state, which appears in our calculation at a significantly higher energy than is found in ref ?. Additionally, the energetic order of the 2^4^Π, 2^4^Δ, and 1^6^Δ states differs in the two calculations, although these states lie within a narrow energy range of less than 0.5 eV in both cases. The most likely sources of these differences are the spin–orbit interaction, which is neglected in our calculations, and the correlation description. The calculations in ref ? use a large configuration-interaction model to obtain the set of excited states, and excitation energies are generally rather sensitive to the level of detail included in such models. Their calculations also assess the influence of the spin–orbit interaction on the vertical excitation energies. They find that the excitation energy of the 1^6^Π state is altered by around 0.3 eV once the spin–orbit interaction is included.

6: Excitation Energies for Electronic States of WH Calculated Using QEC Compared to Those of Ref

Figure shows the calculated excitation cross-sections for the lowest few excited states of WH. The Born correction is applied to each cross-section to account for high angular momenta. Given the high density of target states present in WH, which can increase the required computational resources significantly, the calculation presented here retains all states with vertical excitation energies below 4 eV. Higher lying states are likely to have negligible excitation cross-sections over the energy range below 5 eV that is explored here.

Electronic excitation cross-sections for the lowest sextet and quartet excited states of WH.

Figure shows the total ionization cross-section for WH, calculated using the Binary-Encounter-Bethe (BEB) method.? The ionization energy given by Koopman’s theorem is 6 eV.

Total ionization cross-section for WH.

Photoionization of CH3I

3.5

To complement the scattering calculations above, we include here as an illustration of the new functionality a photoionization calculation for the CH_3_I molecule using the simple SE model.

The all-electron calculation used R-matrix radius of 30 au and 45 B-splines of order 6 up to l max = 6 and the 6-311G** GTO basis set.? The hybrid two-electron integrals were calculated using a combination of a semianalytic approach for the (BG|GG) class and the Legendre expansion with maximum angular momentum 20 for the mixed exchange integrals. The length of the elementary interval for application of the 21-point radial Gauss-Legendre quadrature was Δr = 0.1 a.u.

The calculations with ECPs replaced the innermost 28 electrons with the pseudopotential optimized by Peterson et al.? alongside the cc-pVTZ-PP atomic basis set. The rest of the computational setup was kept the same as that in the all-electron case.

Figure shows the photoionization cross-sections and the angular distributions (β-parameters) for the ionic ground state (X ^2^ E) from the two models in comparison with the experimental measurements of Holland et al.? The results with and without ECPs agree in general with each other, except for the angular distributions at higher energies, where the all-electron data appear to be slightly closer to the experiment. We see only a qualitative agreement of the theory results with the experimental data, which is to be expected given the 1-electron nature of the SE model. The case of CH_3_I is complicated due to the Cooper minimum around 45 eV and, at higher energies, the multichannel effects of coupling to the 4d shell of iodine? and the spin–orbit splitting of the ionic ground state. Accurate description of those effects requires sophisticated modeling of polarization and electron correlation together with inclusion of relativistic effects for the active electrons, which goes beyond the scope of this work.

Photoionization of CH3I into the ground state of the ion. Left: photoionization cross-section. Right: photoelectron angular distributions. Calculations with ECPs are compared to all-electron calculations and experimental measurements.

Conclusions

4

We have implemented ECPs into the molecular scattering and photoionization code UKRmol+ and the QEC interface. The implementation has been thoroughly tested for a range of molecules relevant for plasma modeling and demonstrated the ability of the codes to perform accurate calculations for targets containing heavy atoms. The new functionality is an important extension of the code’s applicability since reliable all-electron bases for heavy atoms are typically not available, thus making all-electron calculations for those elements inaccurate and missing the important (scalar) relativistic effects.

Our implementation allows the use of ECPs with both GTOs and BTOs for the continuum description, the latter with a restriction to the semilocal ECPs. However, the semilocal ECP is often the only and dominant component of the ECP. Using ECPs, scattering calculations for heavy atoms often become computationally comparable to calculations for light elements. Additionally, the implementation of ECP integrals over BTOs allows us to perform scattering and photoionization calculations for extended energy ranges, up to 100 eV, as often required in various applications.

Future work will focus on the inclusion of the spin–orbit ECPs, optimization of the momentum-space integrals over BTOs, and completion of the implementation of the local-type ECP integrals involving BTOs.

Implementation of ECPs also opens the route to representation of molecules embedded in extended environments with the effective electrostatic interaction represented by potentials of Gaussian type.

The approach based on ECPs naturally does not account for the relativistic effects in the valence shell. Instead, those effects must be explicitly included using appropriate terms in the Hamiltonian.? A promising route is the inclusion of the dominant relativistic effects in the nonrelativistic Hamiltonian using perturbation theory, which can be conveniently combined with the R-matrix approach.?