Combining Slater-type orbitals and effective core potentials

TL;DR

This paper introduces a new methodology for evaluating matrix elements of effective core potentials using Slater-type orbitals, enabling more accurate quantum chemical calculations involving core electrons.

Contribution

It develops a general recursive and numerical scheme for matrix element evaluation of ECPs with STOs, including core polarization and spin-orbit effects, and validates it with atomic and molecular data.

Findings

Accurately calculated excitation energies and polarizabilities for alkaline earth metals.

Predicted molecular properties of barium and strontium hydrides consistent with experimental data.

Demonstrated the effectiveness of the method for complex quantum chemical calculations.

Abstract

We present a general methodology to evaluate matrix elements of the effective core potentials (ECPs) within one-electron basis set of Slater-type orbitals (STOs). The scheme is based on translation of individual STO distributions in the framework of Barnett-Coulson method. We discuss different types of integrals which naturally appear and reduce them to few basic quantities which can be calculated recursively or purely numerically. Additionally, we consider evaluation of the STOs matrix elements involving the core polarisation potentials (CPP) and effective spin-orbit potentials. Construction of the STOs basis sets designed specifically for use with ECPs is discussed and differences in comparison with all-electron basis sets are briefly summarised. We verify the validity of the present approach by calculating excitation energies, static dipole polarisabilities and valence orbital energies for the alkaline earth metals (Ca, Sr, Ba). Finally, we evaluate interaction energies, permanent dipole moments and ionisation energies for barium and strontium hydrides, and compare them with the best available experimental and theoretical data.

Tables

Table 1. Table 1: Results of the calculations for the calcium atom (see the main text for technical details). The abbreviation IP stands for first ionisation potential of the system. Small-core PP subtracts 10 electrons from the system (ECP10MDF) whilst large-core PP − - 18 electrons (ECP18SDF). All values are given in wavenumbers, cm -1 .

state	large-core PP	small-core PP		exp${}^{\text{a}}$
	CCSD	CCSD	CC3
3P	15097.0	15173.2	15195.3	15263.1
3D	20941.1	20856.1	21299.6	20356.6
1D	22216.8	22878.6	22859.0	21849.6
1P	23429.8	24845.8	23879.6	23652.3
3S	31651.2	31828.7	31545.5	31539.5
1S	33411.0	33890.9	33336.9	33317.3
IP	49405.2	49821.9	—	49305.9

Table 2. Table 2: Results of the calculations for the strontium atom (see the main text for technical details). The abbreviation IP stands for first ionisation potential of the system. Small-core PP subtracts 28 electrons from the system (ECP28MDF) whilst large-core PP − - 36 electrons (ECP36SDF). All values are given in wavenumbers, cm -1 .

state	large-core PP	small-core PP		exp${}^{\text{a}}$
	CCSD	CCSD	CC3
3P	14579.9	14546.3	14597.2	14702.9
3D	18442.2	18155.0	18393.7	18253.8
1D	20380.4	20584.7	20411.1	20149.7
1P	21451.1	22701.9	21797.5	21698.5
3S	29201.4	29189.7	28939.3	29038.8
1S	30634.4	31063.1	30508.6	30591.8
IP	46006.2	46284.4	—	45932.2

Table 3. Table 3: Results of the calculations for the barium atom (see the main text for technical details). The abbreviation IP stands for first ionisation potential of the system. Small-core PP subtracts 46 electrons from the system (ECP46MDF) whilst large-core PP − - 54 electrons (ECP54SDF). All values are given in wavenumbers, cm -1 .

state	large-core PP	small-core PP		exp${}^{\text{a}}$
	CCSD	CCSD	CC3
3D	9419.4	8923.7	9178.1	9357.8
1D	11609.6	11653.5	11391.4	11395.5
3P	12986.2	12823.6	12925.9	13085.5
1P	17578.9	19527.3	18284.6	18060.3
3S	26281.3	26269.3	26141.9	26160.3
1S	27275.0	—${}^{\text{b}}$	—${}^{\text{b}}$	26757.3
IP	42156.4	42245.8	—	42034.9

Table 4. Table 4: Dipole polarisabilities of the ground state of the calcium, strontium and barium atoms. All values are given in the atomic units.

atom	theory	large-core PP	small-core PP	exp
Ca	HF	164.50	181.60	169$\pm$17${}^{\text{a}}$
	CCSD	170.38	159.14
	CCSD(T)	—	156.12
Sr	HF	205.14	231.94	186$\pm$15${}^{\text{b}}$
	CCSD	221.48	203.16
	CCSD(T)	—	198.52
Ba	HF	280.36	327.48	268$\pm$22${}^{\text{c}}$
	CCSD	323.48	284.70
	CCSD(T)	—	276.62

Table 5. Table 5: Outermost valence orbital energies calculated with the pseudopotentials compared with the reference all-electron Dirac-Hartree-Fock (DHF). All values have their signs reversed and are given in the atomic units.

atom	shell	large-core PP	small-core PP	all-electron DHF${}^{\text{a}}$
Ca	4$s$	0.2064	0.1967	0.1963
Sr	5$s$	0.1930	0.1813	0.1813
Ba	6$s$	0.1760	0.1630	0.1632

Table 6. Table 6: Dissociation energy of the barium hydride (see the main text for technical details) calculated with small-core pseudopotential (ECP46MDF). The abbreviations “ae” and “fc” stand for all-electron and frozen-core, respectively. The quantity in the last column ( Δ ​ fci Δ fci \Delta\mbox{fci} ) is the difference between the dissociation energies calculated at the frozen-core FCI and CCSD(T) levels. The row denoted ∞ \infty lists values extrapolated to the complete basis set. All values are given in wavenumbers, cm -1 .

basis	ae-CCSD(T)	fc-CCSD(T)	$\Delta \text{fci}$	total
wtcc-2	13249.6	14239.9	$+$4.6	13254.2
wtcc-3	15701.2	15975.0	$+$0.3	15701.5
wtcc-4	16393.9	16355.8	$-$1.2	16392.7
wtcc-5	16645.7	16411.9	$-$1.9	16643.8
$\infty$	16903.9	16563.8	$-$2.4	16901.5
exp${}^{\text{a}}$	—	—	—	16350.0
rev exp${}^{\text{b}}$	—	—	—	16910.6

Table 7. Table 7: Dissociation energy of the strontium hydride (see the main text for technical details) calculated with small-core pseudopotential (ECP28MDF). The abbreviations “ae” and “fc” stand for all-electron and frozen-core, respectively. The quantity in the last column ( Δ ​ fci Δ fci \Delta\mbox{fci} ) is the difference between the dissociation energies calculated at the frozen-core FCI and CCSD(T) levels. The row denoted ∞ \infty lists values extrapolated to the complete basis set. All values are given in wavenumbers, cm -1 .

basis	ae-CCSD(T)	fc-CCSD(T)	$\Delta \text{fci}$	total
wtcc-2	12157.2	13221.7	$+$28.1	12185.3
wtcc-3	13561.4	14280.3	$+$31.0	13592.4
wtcc-4	13881.9	14428.6	$+$29.2	13911.1
wtcc-5	13982.5	14474.3	$+$28.7	14011.2
$\infty$	14103.9	14530.6	$+$28.0	14131.9
other theor.	—	—	—	14259.8${}^{\text{a}}$
	—	—	—	14114.6${}^{\text{b}}$

Table 8. Table 8: Molecular properties of strontium and barium hydrides calculated with the small-core pseudopotentials. The abbreviations μ 𝜇 \mu and IP stand for the absolute values of the permanent electronic dipole moment and the (vertical) ionisation potential. IPs and dipole moments are given in units of wavenumbers (cm -1 ) and Debyes (D), respectively.

	SrH		BaH
	ae-CCSD	ae-CCSD(T)	ae-CCSD	ae-CCSD(T)
IP	42707.5	42917.6	38453.6	38791.6
$\mu$	13.49	13.53	14.30	14.38

Equations

$${\chi }_{nlm}(\text{r};\beta )=r^{n-1}{e}^{-\beta r}{Y}_{lm}(\theta ,\varphi ),$$

$$r_b^{n-1}{e}^{-\beta r_b}=\sum _{m=0}^{\infty }\frac{2k+1}{2}{P}_m(\cos {\theta }_a){\zeta }_{nm}(\beta ,r_a;R),$$

$${\zeta }_{nm}(\beta ,r_a;R)={\int }_0^{\pi }d{\theta }_a\sin {\theta }_a{P}_m(\cos {\theta }_a)r_b^{n-1}{e}^{-\beta r_b}.$$

$$\displaystyle\begin{split}\zeta_{n+2,m}&=\big{(}r_{a}^{2}+R^{2}\big{)}\,\zeta_{nm}-\frac{2\,r_{a}\,R}{2m+1}\\ &\times\Big{[}m\,\zeta_{n,m-1}+(m+1)\,\zeta_{n,m+1}\Big{]},\end{split}$$

$$\displaystyle\begin{split}\zeta_{1m}=\frac{\beta\,r_{a}\,R}{2m+1}\Big{[}\zeta_{0,m-1}-\zeta_{0,m+1}\Big{]}.\end{split}$$

$$\displaystyle\zeta_{10}=\beta\,r_{a}\,R\Big{[}\zeta_{00}-\zeta_{01}\Big{]}+e^{-\beta(r_{a}+R)}.$$

$${\zeta }_{0m}=\frac{2\beta }{\pi }{i}_m(\beta r_{<})k_m(\beta r_{>}),$$

$$\displaystyle\begin{split}\frac{1}{\sqrt{\pi}}\,r_{b}^{n-1}\,e^{-\beta r_{b}}\,Y_{lm}(\theta_{b},\phi_{b})&=(-1)^{l-m}\,\big{(}2l+1\big{)}\sum_{\lambda=0}^{l}\frac{r_{a}^{\lambda}}{\sqrt{2\lambda+1}}\,Y_{\lambda m}(\theta_{a},\phi_{a})\,(-R)^{l-\lambda}\begin{pmatrix}\lambda&l-\lambda&\lambda\\ m&0&-m\end{pmatrix}\\ &\times\binom{2l}{2\lambda}^{1/2}\sum_{k=0}^{\infty}\sqrt{2k+1}\,Y_{k0}(\theta_{a},\phi_{a})\,\zeta_{n-l,k}(\beta,r_{a};R),\end{split}$$

$$\begin{array}{rl}{\hat{H}}_v& =\sum _i^{n_v}\left[-\frac{1}{2}{\nabla }_i^2+\sum _a\left[-\frac{Q_a}{r_{ia}}+{\hat{U}}_{\text{mbox}PP}^a(r_{ia})\right]\right]\\ & +\sum _{i<j}^{n_v}\frac{1}{r_{ij}}+\sum _{ab}\frac{Q_a{Q}_b}{r_{ab}}+{\hat{U}}_{\text{mbox}CPP},\end{array}$$

$$\begin{array}{rl}{\hat{U}}_{\text{mbox}PP,av}^a(r_{ia})& =U_L^a(r_{ia})+\sum _{l=0}^{L-1}\sum _{m=-l}^{+l}|lm{⟩}_a\\ & \times \left[U_l^a(r_{ia})-U_L^a(r_{ia})\right]{}_a⟨lm|,\end{array}$$

$$\begin{array}{r}{\hat{U}}_{\text{mbox}PP,so}^a(r_{ia})=\sum _{l=0}^{L-1}\sum _{m=-l}^{+l}\frac{2\Delta U_l^a(r_{ia})}{2l+1}|lm{⟩}_a{\mathbf{l}}_{ia}\cdot {\mathbf{s}}_i{}_a⟨lm|,\end{array}$$

$$r^2{U}_l^a(r)=\sum _k{A}_{kl}^a{r}^{n_{kl}}{e}^{-B_{kl}^a{r}^2},$$

$${\hat{U}}_{\text{mbox}CPP}^a={\hat{U}}_{\text{mbox}CPP}^{[0],a}+{\hat{U}}_{\text{mbox}CPP}^{[1],a}+{\hat{U}}_{\text{mbox}CPP}^{[2],a},$$

$${\hat{U}}_{cpp}^{[0],a}=-\frac{1}{2}{\alpha }_a\sum _{b,c\ne a}{Q}_b{Q}_c\frac{{\mathbf{R}}_{ab}\cdot {\mathbf{R}}_{ac}}{R_{ab}^3{R}_{ac}^3},$$

$$\displaystyle\begin{split}\hat{U}_{cpp}^{[1],a}&=-\frac{1}{2}\alpha_{a}\sum_{i}\frac{1}{r_{ia}^{4}}C^{2}\big{(}r_{ia},\delta_{a}\big{)}\\ &+\alpha_{a}\sum_{i}\sum_{b\neq a}Q_{b}\,\frac{\mathbf{r}_{ia}\cdot\mathbf{R}_{ab}}{r_{ia}^{3}R_{ab}^{3}}\,C(r_{ia},\delta_{a}),\end{split}$$

$${\hat{U}}_{cpp}^{[2],a}=-{\alpha }_a\sum _{j<i}\frac{{\mathbf{r}}_{ia}\cdot {\mathbf{r}}_{ja}}{r_{ia}^3{r}_{ja}^3}C(r_{ia},{\delta }_a)C(r_{ja},{\delta }_a),$$

$$\displaystyle C(r,\delta)=\Big{(}1-e^{-\delta r^{2}}\Big{)}^{\bar{n}},$$

$$\displaystyle\mathcal{F}^{\,0}_{n}\big{(}x,y\big{)}=\int_{0}^{\infty}dr\;r^{n}\,e^{-xr-yr^{2}},$$

$$\displaystyle\mathcal{F}^{\,>}_{n}\big{(}x,y\big{)}=\int_{1}^{\infty}dr\;r^{n}\,e^{-xr-yr^{2}},$$

$$\displaystyle\mathcal{F}^{\,<}_{n}\big{(}x,y\big{)}=\int_{0}^{1}dr\;r^{n}\,e^{-xr-yr^{2}}.$$

$$\displaystyle\mathcal{G}^{\,0}_{n}\big{(}x,y\big{)}=\int_{0}^{\infty}dr\;\ln r\,r^{n}\,e^{-xr-yr^{2}},$$

$$I_{\text{mbox}bca}=⟨{\chi }_{n_b{l}_b{m}_b}({\text{r}}_b;{\beta }_b)|{\hat{U}}_{\text{mbox}PP,av}^c(r_c)|{\chi }_{n_a{l}_a{m}_a}({\text{r}}_a;{\beta }_a)⟩.$$

$$I_{\text{mbox}bca}=I_{\text{mbox}bca}^{\text{mbox}loc}+I_{\text{mbox}bca}^{\text{mbox}nloc}$$

$$\begin{array}{rl}& I_{\text{mbox}bca}^{\text{mbox}loc}=⟨{\chi }_{n_b{l}_b{m}_b}({\text{r}}_b;{\beta }_b)|U_L^c(r_c)|{\chi }_{n_a{l}_a{m}_a}({\text{r}}_a;{\beta }_a)⟩,\end{array}$$

$$\begin{array}{rl}& I_{\text{mbox}bca}^{\text{mbox}nloc}=\sum _{KM}^{L-1}{\int }_0^{\infty }d{r}_c{r}_c^2⟨{\chi }_{n_b{l}_b{m}_b}({\text{r}}_b;{\beta }_b)|KM{⟩}_c\\ & \times \left[U_K^c(r_c)-U_L^c(r_c)\right]{}_c⟨KM|{\chi }_{n_a{l}_a{m}_a}({\text{r}}_a;{\beta }_a)⟩.\end{array}$$

$$\begin{array}{rl}& I_{\text{mbox}aaa}=I_{\text{mbox}aaa}^{\text{mbox}loc}+I_{\text{mbox}aaa}^{\text{mbox}nloc}=\\ & \{\begin{array}{ll}⟨{\chi }_{n_a^{\prime }{l}_a^{\prime }m}({\text{r}}_a;{\beta }_a^{\prime })|U_L^a(r_a)|{\chi }_{n_a{l}_{a}m}({\text{r}}_a;{\beta }_a)⟩,& l_a\ge L\\ ⟨{\chi }_{n_a^{\prime }{l}_a^{\prime }m}({\text{r}}_a;{\beta }_a^{\prime })|U_{l_a}^a(r_a)|{\chi }_{n_a{l}_{a}m}({\text{r}}_a;{\beta }_a)⟩,& l_a<L,\end{array}\end{array}$$

$$\displaystyle\begin{split}&I_{\mbox{\scriptsize aaa}}=\sum_{k}A_{kl}^{a}\,\mathcal{F}^{\,0}_{n_{a}+n_{a}^{\prime}+n_{kl}-2}\Big{(}\beta_{a}+\beta_{a}^{\prime},B_{kl}^{a}\Big{)},\end{split}$$

$$\begin{array}{rl}& \frac{(-1{)}^{l_b}}{\sqrt{2l_a+1}}⟨{\chi }_{n_b{l}_{b}m}({\text{r}}_b;{\beta }_b)|U_L^a(r_a)|{\chi }_{n_a{l}_{a}m}({\text{r}}_a;{\beta }_a)⟩=\\ & \frac{2l_b+1}{2}\sum _{\lambda =0}^{l_b}{\left(\genfrac{}{}{0ex}{}{2l_b}{2\lambda }\right)}^{1/2}\left(\begin{array}{ccc}\lambda & l_b-\lambda & l_b\\ m& 0& -m\end{array}\right){\left(-R\right)}^{l_b-\lambda }{\mathcal{A}}_{\lambda },\end{array}$$

$$\displaystyle\begin{split}&\mathcal{A}_{\lambda}=\sum_{k}\big{(}2k+1\big{)}\begin{pmatrix}\lambda&k&l_{a}\\ 0&0&0\end{pmatrix}\begin{pmatrix}\lambda&k&l_{a}\\ m&0&-m\end{pmatrix}\\ &\times\sum_{m}A_{mL}^{a}\,W_{\lambda+n_{a}+n_{mL}-1,n_{b}-l_{b},k}\big{(}\beta_{a},\beta_{b},B_{mL};R\big{)}.\end{split}$$

$$\displaystyle\begin{split}W_{lmn}\big{(}\alpha,\beta,\gamma;R\big{)}=\int_{0}^{\infty}dr\;r^{l}\,\zeta_{mn}\big{(}\beta,r;R)\,e^{-\alpha r-\gamma r^{2}}.\end{split}$$

$$\displaystyle\begin{split}W_{l,m+2,n}&=W_{l+2,mn}+R^{2}\,W_{lmn}-\frac{2R}{2n+1}\\ &\times\Big{[}n\,W_{l+1,m,n-1}+\big{(}n+1\big{)}W_{l+1,m,n+1}\Big{]},\end{split}$$

Full text

I Introduction

In the last 40 years, Gaussian-type orbitals boys50 ; boys56 (GTOs) have undeniably taken the role of the routine one-electron basis set for ab initio calculations in molecular physics and quantum chemistry. Nonetheless, a considerable interest has remained in the field of Slater-type orbitals (STOs) slater30 ; slater32 or more general exponential-type orbitals (ETOs) homeier91 ; shull59 . This is motivated mainly by the superior analytical properties of STOs (i.e., fulfilment of the nuclear cusp condition kato57 and correct long-range decay agmon82 ; cycon87 ) and their formal simplicity.

The biggest obstacle connected with use of STOs is calculation of many-centre two-electron integrals which are unavoidable in any molecular study. Interestingly enough, there has been a significant progress on this issue in recent years. In fact, looking at only the past 15 years, there are many notable works of Bouferguene et al. bouferguene05a ; bouferguene05b ; bouferguene07 ; bouferguene09 , Rico et al. steinborn00 ; rico00 ; rico01a ; rico01b ; rico04a ; rico04b ; rico05 ; rico06 ; rico08a ; rico08b ; rico13 , Hoggan et al. hoggan01 ; hoggan02 ; hoggan03 ; hoggan04 ; hoggan09 ; hoggan10 , Pachucki pachucki09 ; pachucki12 ; pachucki13a ; pachucki13b ; pachucki16 , and others guidotti03 ; budzinski04 ; weniger00 ; weniger02 ; weniger12 ; harris02 ; barnett00 ; avery04 ; avery14 ; winkler13 ; gebremedhin13 ; vukovic10 ; lesiuk12 ; lesiuk14a ; lesiuk14b ; lesiuk15 ; lesiuk16 . In particular, for the diatomic systems STOs can be now used routinely lesiuk15 .

State-of-the-art ab initio electronic structure calculations are important for the new field at the border of chemistry and physics - studies of ultracold molecules. Experimental advances in laser cooling and trapping of neutral atoms have opened a door for the formation of ultracold diatomic molecules by photoassociation julienne06a , magnetoassociation julienne06b , and vibrational cooling demille10 techniques. To interpret the experimental observations, ab initio calculations of the potential energy curves and coupling matrix elements between the electronic states are crucial. Somewhat surprisingly, spectroscopic and collisional studies of ultracold molecules mostly involve molecules with heavy atoms. See, for instance, Refs. skomorowski12a ; skomorowski12b ; mcguyer13 ; mcguyer15a ; mcguyer15b ; mcguyer16 for joint experimental and theoretical studies of new spectroscopic features of the strontium molecule. Electronic structure calculations can also be used to predict new schemes for the formation of ultracold diatomic molecules moszynski03 ; tomza11 ; tomza13a ; krych11 ; tomza13b ; tomza14 . Accurate interatomic interaction potentials are also of great importance in search for a new physics, see e.g., the work on the YbF molecule which is used in measurements of the permanent electric dipole moment (EDM) of the electron hudson11 , and determination of the proton-electron mass ratio time variation zelevinsky08 . One can also point out the work of Schwertweger et al. schwertweger11 on the Sr2 molecule where time dependence of the fine structure constant is considered. Other examples of physically important diatoms include RbYb molecule nemitz09 ; micheli06 ; zuchowski13 (a promising candidate for quantum computing), BaH lane15 ; tarallo16 ; moore16 (preparation of ultracold hydrogen atoms) and others.

It must be stressed that in a majority of the studies mentioned above, accurate first principles calculations were fundamental in understanding and interpretation of the experimental data. In particular, computations of the potential energy surfaces and the corresponding coupling matrix elements appear to be of prime importance. This is also the area where the Slater-type orbitals are the most advantageous.

Unfortunately, in accurate ab initio calculations for heavy elements one typically encounters two additional major problems. First, the number of occupied orbitals becomes fairly substantial. This, by necessity, calls for extended basis sets with high angular momenta, increasing the overall cost of correlated electronic structure calculations. The second obstacle is the relative importance of the relativistic effects; for heavier elements they are of a similar magnitude (or larger) as the electron correlation contribution visscher96a ; visscher96b . Moreover, additivity of the latter two effects for heavy atoms is at best questionable. visscher96a ; visscher96b

There are several approaches available in the literature to handle the aforementioned problems and most of them are based on the Dirac-Coulomb(-Breit) equations dirac28 ; bethe75 . This is done, e.g., by constructing an approximate four-component spinor expanded in a kinetically balanced basis set dyall91 ; visscher94 ; jong98 ; styszynski00 , or by decoupling the small and large components of the spinor, so that the equations take a familiar two-component form hess85 ; wolf02 ; reiher04 ; reiher06 ; liu05 ; barysz02 ; kedziera04a ; kedziera04b . Another idea developed independently relies on the so-called regular approximations lenthe93 ; lenthe96 . In this paper we consider the effective core potential (ECP) approach dolg12 which may be viewed as a slightly less rigorous method than the former ones. However, little accuracy is typically sacrificed (especially for weakly bound systems) and the calculations simplify to a great extent.

The fundamental idea behind ECPs is that the inner core orbitals of heavy elements are inert and do not change significantly in different chemical environments. Therefore, their influence on the valence space can be modelled with a proper pseudopotential (PP) foot1 which is, by definition, universal for a given element. This leads to an approximate two-component valence relativistic wavefunction, obtained as an eigenfunction of the valence only Hamiltonian. This approach has two unquestionable advantages. First, the inner core orbitals are removed from explicit consideration, so that the size of the one-electron basis set is considerably reduced. Second, the scalar relativistic effects can be straightforwardly included in the pseudopotential (by a proper parametrisation).

The main goal of this work is to combine the methodology of effective core potentials with the one-electron basis set of Slater-type orbitals. We propose a general method to evaluate all necessary matrix elements by using analytical or seminumerical techniques. Efficiency of the proposed algorithm is sufficient to perform general large-scale calculations. Further in the paper, we consider the so-called core polarisation potentials fuentealba82 ; muller84a ; muller84b which rely on the assumption that the core is additionally polarisable. This captures the first-order effects of the implicit core-valence correlations and significantly improves the description when the large-core pseudopotentials are used. We also briefly consider effective spin-orbit pseudopotentials cohen79 ; ermler81 ; pitzer88 which allow for an approximate computation of the spin-orbit splittings and couplings. Finally, we present results of test calculations for alkaline earth metals (Ca, Sr, Ba) and predict properties of the corresponding hydrides.

In the paper we rely on the known special functions to simplify the derivations and the final formulae. Our convention for all special functions appearing below is the same as in Ref. stegun72 .

II Theory

In this section we introduce some important formulae which are used further in the paper. This is necessary to introduce the notation and specify precise meaning of several quantities. This short mathematical preface may be useful for readers who are not entirely familiar with employed concepts.

II.1 Slater-type orbitals and spatial translations

Slater-type orbitals (STOs) slater30 ; slater32 have the following general form

[TABLE]

where $Y_{lm}$ are the spherical harmonics defined according to the Condon-Shortley phase, $n$, $l$ are nonnegative integers satisfying $n>l$, and $\beta>0$ is a real parameter. The orbitals defined above are not normalised; we find this convention to be more robust for the purposes of the paper. In order to restore the proper unity normalisation Eq. (1) has to be multiplied by $S_{n}(\beta)=(2\beta)^{n+1/2}/\sqrt{(2n)!}$.

Throughout the paper we shall extensively use the translation method for STOs in order to shift them to a convenient point in space. This is achieved with help of the famous Barnett-Coulson $\zeta$-function method coulson37 ; barnett51 ; barnett63 . Translation of the $ns$ STOs is given by the following two-range formula

[TABLE]

where $P_{m}$ are the ordinary Legendre polynomials, $R$ is the distance between centres $a$ and $b$, and $\zeta_{nm}$ are given by the integral representation

[TABLE]

From now on, we drop the parentheses from the definition of the $\zeta$-function, i.e., it is assumed that $\zeta_{nm}:=\zeta_{nm}(\beta,r_{a};R)$ unless explicitly stated otherwise.

The above formal definitions are not particularly useful in actual applications. Instead, the following recursive relations provide a starting point for further developments

[TABLE]

and

[TABLE]

The last formula is not valid for $m=0$ and the explicit expression should be used instead

[TABLE]

To initiate the recursive process one requires the following starting values

[TABLE]

where $i_{m}$ and $k_{m}$ are the modified spherical Bessel functions of the first and second kind stegun72 , respectively, and $r_{<}=\min(r_{a},R)$, $r_{>}=\max(r_{a},R)$. For convenience of the reader, we gathered all properties of the modified spherical Bessel functions which are important here in the Supplemental Material supplement . Equations in the Supplemental Material are referenced with prefix “S”, e.g., the sixth equation in the Supplemental Material is referenced as Eq. (S6).

In order to spatially shift STOs of the form (1) one needs to combine Eq. (2) with the well-known translation formula for the regular solid harmonics, Ref. arfken , pp. 797. This leads to

[TABLE]

with the usual notation for the Wigner $3J$ symbols, Ref. arfken , pp. 270.

A short comment on the nature of the above expression is necessary. It is well-known that the strongest drawback of the Barnett-Coulson method is that it leads, in general, to infinite series. Worse than that, these series tend to converge extremely slowly; in some cases a logarithmic convergence pattern has been observed lowdin56 ; flygare66 . This limits the applicability of the method significantly and forced some authors to apply convergence acceleration techniques petersson67 ; bouferguene98 . Fortunately, this issue is absent in all final formulae derived in this paper. In most cases, the sum over $k$ truncates as a result of the triangle conditions for the Wigner $3J$ symbols, see Ref. arfken , pp. 803.

Before the end of the present section we would like to point out that there exist some other methods for translation of STOs, including one-range and two-range formulae, Refs. danos65 ; filter80 ; homeier92 ; novosadov ; steinborn80 ; weniger85 ; guseinov76 ; guseinov78 ; guseinov80 ; guseinov85 ; guseinov88 , yet we have not found these alternative formulations to be particularly advantageous in the present case compared to the standard Barnett-Coulson scheme, Eq. (8). General theory of addition theorems was given in a pedagogical way by Weniger weniger00 ; weniger02 . Note that mathematical correctness (and usefulness) of some formulations of the addition theorems is still subject to a debate weniger12 .

II.2 Pseudopotentials parametrisation

As already mentioned in the introduction, in calculations involving ECPs one considers the following valence-only Hamiltonian dolg12

[TABLE]

where $i,j,...$ denote the electrons, $a,b,...$ denote the nuclei, $\hat{U}_{\mbox{\scriptsize PP}}^{a}$ is the pseudopotential of the core $a$ with charge $Q_{a}$, and $n_{v}$ is the number of valence electrons. The term arising from the core polarisation potential ($\hat{U}_{\mbox{\tiny CPP}}$) will be specified further in the text.

Let us briefly discuss the construction of the effective core potentials. They are divided into the spin-averaged and spin-dependent terms, $U_{\mbox{\scriptsize PP}}^{a}=U_{\mbox{\scriptsize PP,av}}^{a}+U_{\mbox{\scriptsize PP,so}}^{a}$. Typically, the first term is included explicitly in the electronic structure calculations whilst the second is treated perturbatively. Both of these potentials are represented in a semi-local form

[TABLE]

and

[TABLE]

where $L$ is the highest angular momentum of the orbitals in the core $a$, $\mathbf{l}_{ia}$ is the (orbital) angular momentum operator corresponding to the centre $a$ and the electron $i$, $\mathbf{s}_{i}$ is the spin operator of the electron $i$, and $\langle lm|_{a}$ are projection operators on the spherical harmonics $Y_{lm}$ placed at the centre $a$. Presence of the projection operators assures that orbital components of different angular momenta connect with proper radial functions. Parenthetically, it is observed that $U_{l}^{a}(r_{ia})$ are nearly identical for $l>L$ which justifies the rearrangements in Eqs. (10) and (11).

To specify a family of pseudopotentials a precise form of the radial components, $U_{l}^{a}(r)$, must be given. It is very common to use a short linear combination of the radial Gaussian functions dolg12

[TABLE]

where $n_{kl}$, $A_{kl}^{a}$ and $B_{kl}^{a}$ are adjustable parameters. Their determination for a given atom is far from trivial and strategies of the so-called energy-consistent preuss55 ; topp73 ; flad79 ; preuss81 , shape-consistent phillips59 ; kahn76 ; baran88 ; wang00 , and other dolg12 pseudopotentials were developed.

II.3 Effects of the core polarisation

The so-called core polarisation potentials muller84a ; muller84b ; fuentealba82 (CPP) constitute a method to improve upon the approximations underlying the ordinary ECPs. The core is allowed to be polarisable, i.e., reorientation of valence electrons in a molecular environment creates an induced dipole moment of the core. By simple electrostatic arguments, the value of this dipole moment is assumed to be proportional to the strength of the electric field at the core. This gives rise to the total potential $\hat{U}_{\mbox{\tiny CPP}}=\sum_{a}\hat{U}_{\mbox{\tiny CPP}}^{a}$ in the form

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the consecutive terms are the scalar, one- and two-electron components. In the above expression $\alpha_{a}$ is the polarisability of the core $a$, determined from separate theoretical calculations or by semiempirical adjustment; $C(r_{ia},\delta_{a})$ is the cutoff function assuring that the potential is regular when electron $i$ is at the core $a$. The form of the cutoff function as well as the value of the cutoff parameter $\delta_{a}$ are arbitrary. The following expression is frequently used

[TABLE]

where $\bar{n}$ is either $1$ (Stoll-Fuentealba fuentealba82 ) or $2$ (Müller-Meyer muller84a ; muller84b ). The optimal values of $\delta_{a}$ are determined by numerical experimentation for each atom separately.

II.4 Basic integrals

Further in the text we show that all matrix elements involving averaged and spin-orbit pseudopotentials can be expressed through the following family of one-dimensional integrals

[TABLE]

For evaluation of the core polarisation potentials matrix elements one additionally requires integrals with a logarithmic (albeit integrable) singularity, e.g.

[TABLE]

and similarly for $\mathcal{G}^{\,>}_{n}$ and $\mathcal{G}^{\,<}_{n}$.

The issue of calculation of the above integrals is fairly technical and marred with numerical problems. For completeness, in the Supplemental Material we present an exhaustive description of the methods we recommend for calculation of these basic quantities. Additionally, a special case of these integrals were considered in Ref. silkowski15 . Note that Eqs. (18) and (20) are well-defined only for $n\geq 0$; this restriction does not hold for Eq. (19).

III Spin-averaged and spin-orbit matrix elements

In the remainder of the paper we consider calculation of ECPs matrix elements for diatomic systems only. This is mainly because the issue of exact calculation of the two-electron integrals for polyatomic molecules has not been fully resolved yet. Consequently, we adopt a coordinate system where both atoms are located on the $z$ axis and separated by a distance $R$.

For calculations on general polyatomic systems one needs the following matrix elements involving the spin-averaged potentials

[TABLE]

By using Eq. (10) and after simple manipulations one can rewrite the above expression as

[TABLE]

where

[TABLE]

and

[TABLE]

The approach adopted here depends on the relative location of the centres.

The main difficulty connected with the calculation of the pseudopotentials matrix elements is the presence of the projection operators if the orbitals are placed on different centres. Additionally, pseudopotentials are typically parametrised in terms of the Gaussian-type expansions, Eq. (12), which leads to mixed Slater-Gaussian type basic integrals. The latter are usually not easily expressible through the standard elementary and special functions, and new techniques need to be developed to handle them.

III.1 Spin-averaged potentials, $I_{\mbox{\scriptsize aaa}}$ type

Let us first consider the atomic case, $a=b=c$. Due to orthogonality of the spherical harmonics the matrix element simplifies to

[TABLE]

provided that $l_{a}=l_{a}^{\prime}$ and $m_{a}=m_{a}^{\prime}=m$. Otherwise, the result vanishes due to the spherical symmetry of the integrand. Evaluation of the remaining integrals is now elementary; making use of Eqs. (1) and (12)

[TABLE]

where \mathcal{F}^{\,0}_{n}\big{(}x,y\big{)} is defined by Eq. (18).

III.2 Spin-averaged potentials, $I_{\mbox{\scriptsize baa}}$ type

Let us now consider the first of the two-centre matrix elements, $I_{\mbox{\scriptsize baa}}$. One can easily see that they obey formally the same expression as (26), but $l_{a}$ and $l_{b}$ do not need to be equal. However, the requirement $m_{a}=m_{b}$ still holds as a consequence of the axial symmetry. Translating STO from the centre $b$ to the point $a$

[TABLE]

where

[TABLE]

The remaining one-dimensional integration is confined to the following formula

[TABLE]

A straightforward approach to Eq. (30) is to use a quadrature of some sort and treat the integrals in a purely numerical fashion. However, the integrand possesses a derivative discontinuity (i.e., a cusp) at $r=R$. This makes integration with standard Gaussian quadratures difficult. For a reasonable performance one would need to divide the integration range into two subintervals, \big{[}0,R\big{]} and \big{[}R,\infty\big{]}, and treat each of them separately, possibly with different integration rules. This, in turn, increases the computational costs as the integration needs to be performed for a large set of parameters $l$, $m$, $n$.

An alternative approach which we propose here relies on the recursive relations (4) and (5). By inserting them into the definition (30) one arrives at

[TABLE]

and

[TABLE]

To initiate the above recursions one needs $W_{l0n}$ which can be expressed with help of Eq. (7) as

[TABLE]

where

[TABLE]

The latter two integrals can be integrated numerically to a very good precision. This approach is reasonable if one does not care about the timings of the calculations (e.g., for benchmark purposes). However, to reach a computational cost comparable with GTOs and use ECPs efficiently for large systems better procedures are required, preferably recursive. They are described in detail in next paragraphs.

Let us begin with the first class of integrals, $\mathcal{I}_{ln}$ (we drop the parentheses from now on). By using the relation (S5) one arrives at

[TABLE]

This recursion needs to be carried out in the direction of decreasing $n$ in order to maintain the numerical stability. To start the process (37) one requires $\mathcal{I}_{lN}$ for two neighbouring (large) $N$ and $\mathcal{I}_{0n}$. We propose to evaluate both of them by inserting the power series expansion of $i_{n}(x)$ around the origin into the definition (35)

[TABLE]

Since the above summation is infinite and for practical reasons needs to be truncated, it is helpful to estimate in advance how many terms are required to achieve convergence.

We first note that the rate of convergence of Eq. (38) is not significantly affected by a change of values $l$, $\alpha$, and $\gamma$, the only important variables being $\beta$, $n$. The sum (38) converges faster when $\beta$ decreases or $n$ increases. Therefore, we can consider the worst-case scenario of $\mathcal{I}_{00}$ as a function of $\beta$. Making use of the relationship \mathcal{F}^{\,<}_{n}\big{(}\alpha,\gamma\big{)}\leq\frac{1}{n+1} one arrives at the formal upper bound

[TABLE]

One can assume that the convergence pattern of the above series is very similar to the original $\mathcal{I}_{00}$. The number of terms necessary to achieve convergence for a given $\beta$ can be estimated by solving the equality $\beta^{2k}=\epsilon\,(2k)!\,(2k+1)^{2}$ and rounding up to the closest integer value ($\epsilon$ is the prescribed accuracy goal). We obtained numerical solutions of Eq. (39) for a finite set of $\beta$ and fitted them with a linear function, giving $n_{\mbox{\scriptsize terms}}=0.68\beta+29.5$. This estimation is reliable for all $\beta$, but it tends to overshoot $n_{\mbox{\scriptsize terms}}$ slightly, especially for smaller $\beta$.

The method based on the infinite summation is quite successful for small and moderate $\beta$ but becomes tedious when the values of $\beta$ get large. It typically occurs for stretched molecules or for extended basis sets with high exponents. To avoid laborious summations in such situations we present a large $\beta$ asymptotic expansion of the functions $\mathcal{I}_{ln}$. The derivation begins by rewriting Eq. (35) as a difference of two integrals over the intervals $[0,\infty]$ and $[1,\infty]$. In the first integral one needs to exchange the variables to $\beta x$ and subsequently expand the Gaussian function under the integral sign in power series. The remaining integral can be recognised as the Legendre function of the second kind $Q_{n}$ by means of the analytic continuation. This finally leads to the asymptotic formula for the first part

[TABLE]

where the subscript in $Q_{n}$ denotes differentiation with respect to the main argument. Calculation of the Legendre functions and their derivatives is a standard task as has been discussed many times in the literature olver83 ; gil98 ; schneider10 . Let us pass to the second part, i.e., the integral over $[1,\infty]$. Note that in this integral (contrary to the former) the argument of the Bessel function is always large for large $\beta$. Therefore, we can use the large-argument expansion of the Bessel function given by Eq. (S7) in Supplemental Material. This straightforwardly leads to the formula

[TABLE]

By combining Eqs. (40) and (41) one obtains the final large-$\beta$ asymptotic expansion of the integrals (35). Explicit expressions for the coefficients in the above expression are given in Supplemental Material.

Passing to the second class of integrals, $\mathcal{K}_{ln}$, and inserting the explicit formula for $k_{n}(z)$, Eq. (S3), leads to

[TABLE]

which makes the evaluation elementary. Note that all terms in the above sum are strictly positive, so that no cancellations are possible and the final result acquires the same accuracy as the supplied values of ${F}^{\,>}_{l}$.

III.3 Spin-averaged potentials, $I_{\mbox{\scriptsize bab}}$ type

In the case of $I_{\mbox{\scriptsize bab}}$ configuration there are no simplifications analogous as in Eq. (26) and we must use Eqs. (24) and (25) as they stand. Therefore, the local and nonlocal parts need to be treated separately in this case. Considering the local part, note that in Eq. (24) both STOs occupy the same centre. Therefore, one can expand the product of two STOs into a linear combination of STOs by using standard relations for coupling of the angular momenta. As a result, the integrals $I_{\mbox{\scriptsize bab}}^{\mbox{\scriptsize loc}}$ can easily be expressed in terms of

[TABLE]

where we have made use of Eq. (8) and integrated over the angles. In the above expression a handful of quantities has been introduced, i.e., $n_{ab}=n_{a}+n_{b}-1$, $\beta_{ab}=\beta_{a}+\beta_{b}$, and $|l_{a}-l_{b}|\leq l_{ab}\leq l_{a}+l_{b}$ (note that $n_{ab}>l_{ab}$).

The remaining one-dimensional integrals in Eq. (43) are of the same class as defined by Eq. (30) but with $\alpha=0$. Theoretically, this brings a degree of simplification and allows for a more robust scheme. However, we found that it is not worth increasing the size and complexity of the code by including separate routines for the case $\alpha=0$. Therefore, we recommend that the case $\alpha=0$ is treated with general techniques described above. There are no singularities or numerical instabilities in these expressions as $\alpha$ approaches zero, so that the codes can be reused with no changes.

Let us now consider the calculation of the nonlocal term, $I_{\mbox{\scriptsize bab}}^{\mbox{\scriptsize nloc}}$. This case is much more troublesome due to the fact that the coupling of the angular momenta cannot be used before the translation of the orbitals. Therefore, both STOs need to be shifted independently from the centre $b$ to the centre $a$. After some algebra one finds

[TABLE]

where

[TABLE]

and analogously

[TABLE]

Finally, the innermost integral can be expressed as

[TABLE]

reducing all matrix elements to a definite one-dimensional integration. Let us note at this point that the integrals Eq. (47) are invariant with respect to permutation $1\leftrightarrow 2$ of all indices (including the nonlinear parameters), i.e., $U_{l\,n_{1}n_{2}}^{m_{1}m_{2}}(\beta_{1},\beta_{2},\gamma;R)=U_{l\,n_{2}n_{1}}^{m_{2}m_{1}}(\beta_{2},\beta_{1},\gamma;R)$.

Clearly, the integrals $U_{l\,n_{1}n_{2}}^{m_{1}m_{2}}$ are the most complicated quantities appearing in the theory. Since they are five-index objects, any numerical integration is expected to be prohibitively expensive. Therefore, the recursive techniques are preferred despite the necessity to operate in many dimensions. Derivation of the recursive formulae for the basic integrals $U_{l\,n_{1}n_{2}}^{m_{1}m_{2}}$ follows along a line similar as in the previous subsection. Let us insert Eqs. (4) and (5) into the definition of $U_{l\,n_{1}n_{2}}^{m_{1}m_{2}}$. After some rearrangements one obtains

[TABLE]

and

[TABLE]

The exceptions from the above relation are the integrals with $n_{1}=0$ which have to be calculated according to Eq. (6) instead

[TABLE]

The recursion relations which allow to increase the second pair of indices can be obtained by using the aforementioned symmetry property.

The above relations allow to calculate $U_{l\,n_{1}n_{2}}^{m_{1}m_{2}}$ with nonzero $m_{1}$, $m_{2}$ starting solely with the integrals $U_{l\,n_{1}n_{2}}^{00}$. The latter obey the relationship

[TABLE]

where

[TABLE]

which results directly from Eqs. (4)-(7). Let us note that some of the indices of $U_{l\,n_{1}n_{2}}^{m_{1}m_{2}}$ must be increased only few times at most. In fact, the maximal value of $6$ for the indices $n_{1}$, $n_{2}$ is sufficient to cover the whole known periodic table. Moreover, in accurate calculations with Slater-type orbitals for light systems lesiuk15 one typically uses even-tempered sequences of functions with $n=l+1$. This reduces the necessary values of $n_{1}$, $n_{2}$ to $1$. A similar observation is valid for the $W_{lmn}$ integrals defined in the previous section, Eq. (30).

Evaluation of the integrals (52) and (53) follows a very similar strategy as adopted previously. By using the power series expansion of $i_{n}(z)$ one easily arrives at

[TABLE]

The corresponding expression involving the second pair of indices is obtained by using the symmetry relation $\mathcal{J}_{n_{1}n_{2}}^{l}(\beta_{1},\beta_{2},\gamma)=\mathcal{J}_{n_{2}n_{1}}^{l}(\beta_{2},\beta_{1},\gamma)$. Both these formulae are useful for small or moderate $\beta_{1}$ or $\beta_{2}$, but fail otherwise due to slow convergence of the infinite series. In this case one needs the large $\beta_{1}$ or $\beta_{2}$ asymptotic expansion which can be derived analogously as Eqs. (40) and (41).

Finally, evaluation of the second class of integrals $\mathcal{K}_{n_{1}n_{2}}^{l}$ relies on the explicit expression for the modified Bessel functions, Eq. (S3). By inserting it twice into the definition (53) and rearranging one obtains

[TABLE]

III.4 Spin-orbit potentials

The effective spin-orbit potentials are of very similar form as the scalar pseudopotentials. In fact, they differ only due to presence of the angular momentum and spin operators, Eq. (11). Additionally, there is no local part in the spin-orbit pseudopotentials. After some manipulations one can show that the necessary matrix elements

[TABLE]

can be rewritten without a loss of generality as

[TABLE]

To derive this expression one uses the fact that the projection operators ${}_{c}\langle lm|$ are idempotent and that they commute with the spin-orbit operator. The imaginary unit has been added to make all matrix elements real as the orbital angular momentum operator is, in general, complex valued. The only new objects present in Eq. (57) are matrix element of the angular momentum operator, $\langle lm|\,\mathbf{l}\cdot\mathbf{s}\,|lm^{\prime}\rangle$. Explicit expressions for these integrals can be derived with standard algebra of the angular momentum (see Ref. arfken , pp. 793).

Standard quantum chemistry packages compute all basic matrix elements over spatial orbitals and the spin component is added later by proper construction of an approximate wavefunction. This is the approach we adopt here. The integrals (56) and (57) are evaluated for all Cartesian components separately and stored for further manipulations.

IV Core polarisation matrix elements

In order to evaluate the core polarisation correction to the Hamiltonian, Eqs. (13)-(16), only two distinct matrix elements are necessary. They read

[TABLE]

with $i=1,2$, and

[TABLE]

The Gaussian factors in these definitions come from the adopted cutoff function, Eq. (17). Note that instead of Cartesian coordinates in Eq. (59) we use pure spherical components corresponding to $M=-1,0,+1$. The total contribution to the Hamiltonian can be assembled by combining these matrix elements with geometric and molecular data according to Eqs. (13)-(16).

Starting with the atomic-type integrals, one can straightforwardly integrate over the angles in the spherical coordinate system, giving after some rearrangements

[TABLE]

provided that $m_{a}^{\prime}=m_{a}+M$ (otherwise the result vanishes). The form of the expression for the matrix element involving $\hat{V}_{\mbox{\scriptsize CPP}}^{(2)}(r)$ depends on the value of $n_{a}+n_{a}^{\prime}$. It reads

[TABLE]

for $n_{a}+n_{a}^{\prime}\geq 4$,

[TABLE]

[TABLE]

for $n_{a}+n_{a}^{\prime}=3$ and $n_{a}+n_{a}^{\prime}=2$, respectively. Let us remind that the above matrix elements are nonzero if and only if $l_{a}=l_{a}^{\prime}$ and $m_{a}=m_{a}^{\prime}$.

Passing to the two-centre matrix elements, we first note that calculation of $I_{\mbox{\scriptsize baa}}^{\mbox{\scriptsize CPP}(1)}$ and $I_{\mbox{\scriptsize bab}}^{\mbox{\scriptsize CPP}(1)}$ is almost exactly the same as for the local components of the spin-averaged potentials described in Sections III.2 and III.3. Thus, there is no need to repeat the details of the derivation and we provide only a short sketch for convenience of the reader. Considering $I_{\mbox{\scriptsize baa}}^{\mbox{\scriptsize CPP}(1)}$, the major difference as compared with the derivation given in III.2 is that two spherical harmonics placed on the centre $a$ need to be coupled first. Next, translation of the STO from the centre $b$ to the centre $a$ enables to integrate over the angles and the Jacobian cancels the apparent $1/r^{2}$ singularity introduced by the potential (59). This allows to expand the Gaussian damping function with help of the binomial theorem and the final result is written as a linear combination of the $W_{lmn}$ integrals defined by Eq. (30). A similar conclusion is found for the $I_{\mbox{\scriptsize bab}}^{\mbox{\scriptsize CPP}(1)}$ integrals class. Two STOs present on the centre $a$ need to be expanded into a linear combination of STOs giving an analogue of Eq. (43). Once we translate the distribution from centre $b$ to centre $a$ and integrate over the angles, the singularity vanishes and rest of the derivation is straightforward. The final result can also be written in terms of the integrals (30).

Unfortunately, calculation of the matrix elements involving the potential $\hat{V}_{\mbox{\scriptsize CPP}}^{(2)}(r)$ is more involved. This is due to the fact that the apparent singularity is not automatically cancelled by the Jacobian and thus the damping factor in Eq. (60) cannot be expanded that easily. As a result, in addition to the ordinary integrals introduced in Sec. III.2 and III.3 one requires

[TABLE]

where the analogy with Eqs. (35) and (36) is obvious. However, the values of $l$ are not restricted to nonnegative integers here since $l=-1,-2$ are also necessary. For the integrals $\widetilde{\mathcal{K}}_{ln}^{p}$ this is not problematic because of the integration range. Only the integrals $\widetilde{\mathcal{I}}_{-1,n}^{p}$, and $\widetilde{\mathcal{I}}_{-2,n}^{p}$ are troublesome. To bring them into a closed form we introduce the following quantities

[TABLE]

so that $\widetilde{\mathcal{I}}_{-1,n}^{p}$ and $\widetilde{\mathcal{I}}_{-2,n}^{p}$ can now be simplified by integration by parts. This gives

[TABLE]

and

[TABLE]

where the notation for the nonlinear parameters \big{(}\alpha,\beta,\gamma\big{)} has been suppressed when it is clear from the context. Finally, calculation of the integrals $\mathcal{M}_{ln}$ is reminiscent of the methods introduced in Sec. III.2. For example, for small and moderate $\beta$

[TABLE]

This finalises the present section of the paper.

V Numerical examples

Throughout the paper we use atomic units for calculated quantities unless explicitly stated otherwise. The approximate conversion factors are $1\,a_{0}=0.52\,917\,$Å for lengths (Bohr radius), $1\,\mbox{a.u.}=219\,474.63$ cm*-1* for energies, and $1\,\mbox{a.u.}=2.54\,158$ Debye (D) for dipole moments.

V.1 Basis set optimisation

While there are many families of pseudopotentials available in the literature, the same cannot be said about the relevant Slater-type basis sets. Therefore, we performed optimisation of the valence STOs basis sets for three elements - calcium, strontium and barium (Ca, Sr, Ba). The last known element of the rare-earth metals (radium, Ra) is not considered here because it is highly radioactive and thus not enough confirmed experimental data is available to constitute a comprehensive test case. For all elements we adopted the Stuttgart-Dresden family of energy consistent pseudopotentials. The so-called small-core pseudopotentials (10 valence electrons) are described in Ref. lim06 whilst the large-core counterparts (2 valence electrons) are given in Ref. fuentealba85 .

In general, construction of the STOs basis sets for pseudopotential calculations is similar as in the recent paper concerning the beryllium dimer lesiuk15 . Therefore, we shall not repeat the minutiae of the procedure and illuminate only the most important differences. First, instead of the conventional even-tempered stencil for the nonlinear parameters (exponents) of each angular momentum we use the following extended scheme (well tempering)

[TABLE]

where $i=0,\ldots,n_{l}$, $l$ is the angular momentum and $\alpha_{l}$, $\beta_{l}$, $\gamma_{l}$, $\delta_{l}$ are variational parameters optimised for each $l$. For $l>2$ we set $\delta_{l}=0$ to reduce the number of parameters. The second difference is the choice of the target function - total atomic valence correlation energy, i.e., we do not freeze any additional orbitals in the valence space. Let us mention that there are many similarities between the basis set optimisation procedures in the all-electron systems and for the valence-only pseudopotential. However, the latter case is much more technically challenging. This is mainly due to occurrence of numerous local minima and problematic behaviour of the pseudoorbitals near the nucleus causing linear dependencies problem.

The basis sets optimised in the course of the present work are constructed according to the correlation consistency principle dunning89 . They are abbreviated wtcc-$l$ (well-tempered correlation-consistent), where $l$ is the highest angular momentum present in the basis set. For example, for the valence-only ten electron systems (small-core pseudpotentials) the smallest basis set (wtcc-2) has composition $10s8p3d$ and the largest (wtcc-5) - $13s11p7d5f4g2h$. This includes two sets of additional diffuse functions which were trained to maximise the atomic polarisability calculated at the closed-shell Hartree-Fock level. All basis sets used in this work can be obtained from the authors upon request.

V.2 Test results

In order to check the accuracy of the new basis sets and correctness of the procedures given in this work we performed extensive numerical tests. For each atom (Ca, Sr, Ba) we evaluated the first three excitation energies and the first ionisation potential (IP). The results are given in Tables 1-3. Additionally, in Tables 4 and 5 we provide ground-state dipole polarisabilities (static) and outermost $ns$ valence orbitals Hartree-Fock energies, respectively. All calculations were performed both with large- and small-core pseudopotentials (2 and 10 valence electrons, respectively). In the case of the large-core pseudopotentials the corresponding core polarisation potential was included by default. All valence two-electron calculations were performed with the CCSD method purvis82 and its variants for the excited and ionised states (EOM, IP-EOM geertsen89 ; stanton93 ; piecuch02 ).

For the 10 electron systems (small-core pseudopotentials) the calculations are slightly more involved. For the excited states we used the EOM-CC3 method koch97 as implemented in the code for excited state properties of Tucholska et al tucholska14 ; tucholska17 with all orbitals active. For the ionised states we used the IP-EOM2 method gour05 , and the polarisabilities were evaluated at the CCSD and CCSD(T) raghavachari89 levels by using two-point finite difference method with displacement of 10*-4* a.u. All calculations were performed with help of a locally modified versions of the Gamess gamess1 ; gamess2 and AcesII aces2 program packages, with an exception of the computations at the CC3 level of theory where we used a program written by one of us (AMT). In all calculations presented in this section the largest basis sets available in each case are used - wtcc-5 for the small-core pseudopotentials and wtcc-3 for the large-core counterparts.

Let us begin the analysis with the atomic excitation spectra and consider the strontium atom as an example. The overall picture is more-or-less the same for the remaining atoms and we shall comment on the differences further in the text. One can see that both the small-core and large-core pseudopotentials give a very good agreement with the experimental data. However, the small-core pseudopotential combined with the CC3 method performs better, as could have been expected. The average deviation from the experimental data is around 0.6% for the small-core, and 0.9% for the large-core potentials. One can safely say that the ECP-MDF/CC3 level of theory is very reliable. On average, excitation energies are expected to be less than 1% away from the experimental data. Additionally, no significant increase of the error is observed for any particular spatial symmetry or spin state. This suggests that the new basis sets have no inherent bias, which is a desirable feature in a molecular work.

Excitation energies for barium are in only slightly worse agreement with the experiment than in the case of strontium. The average error is around 0.9% for the small-core and 1.7% for the large-core pseudopotentials. Unfortunately, we observe a significant error for the 3D and 1D states of calcium with both pseudopotentials. This behaviour is surprising because the remaining excitation energies are in a good agreement with the experiment. Therefore, our first suspicion was that 3D and 1D states are highly diffused and the basis set is not saturated well enough. However, we found that further extension of the basis set changed the results by less that 100 cm*-1* which is not enough to explain the discrepancy. As a result, we presume that this increase in the error is an inherent problem of the given pseudopotentials. We note that in the original papers describing the pseudopotentials fuentealba85 ; lim06 errors obtained for Ca were in fact significantly larger than for the other elements.

Let us also compare our results for the strontium atom with the values obtained by Skomorowski et al. skomorowski12b In this work the same pseudopotential (ECP28MDF) was used in combination with a custom-made GTOs basis set and the EOM-CC3 method. Both basis sets are roughly of the same size, so a fair comparison is possible. Skomorowski et al. skomorowski12b give 14570.8 cm*-1* and 21764.3 cm*-1* for the nonrelativistic 3P and 1P states, respectively. These results are very similar to the values given in Table 2; any differences are probably accidental, suggesting that both basis sets are of a similar quality for the P states. However, the situation is different for the D states. The authors of Ref. skomorowski12b report 18668.8 cm*-1* for the 3D state and 20650.3 cm*-1* for the 1D state. Clearly, errors with respect to the experimental values are much larger than for the P states, and also by few hundreds cm*-1* larger than calculated with our basis sets (cf. Table 2).

Next, we would like to check the quality of the basis sets for properties different than the atomic spectra. First, let us consider the static dipole polarisabilities calculated with both families of pseudopotentials. The results are given in Table 4. The large-core pseudopotentials underperform considerably - the calculated values differ by more than 10% from the experimentally determined ones (and lie outside the corresponding error bars). The only exception is the calcium atom, but this agreement is probably accidental. A completely different picture is found for the small-core pseudopotentials. Here, calculated values are reasonably close to the experiment and lie within the given error bars. We estimate the basis set error to be smaller than $1$ a.u. by observing the effect of additional diffuse functions. Omission of the higher cluster operators brings an uncertainty of $1-2$ a.u. assuming that the results converge geometrically with the excitation level. Therefore, one can expect that the theoretical limits are $2-3$ au. below the values given in Table 4. This is still slightly above the experiment for Sr and Ba and somewhat below for Ca. The remaining discrepancy might be a result of an inherent pseudopotential error or a systematic error in the experimental data.

Lastly, we would like to consider the outermost valence $ns$ orbital energies calculated with the pseudopotentials and compare them with all-electron Dirac-Hartree-Fock (DHF) values which we treat as a reference. Note that this quantity is very important for chemical bonding phenomena and it is connected with some important descriptors such as the electronegativity, etc. The results are given in Table 5. Remarkably, the small core pseudopotentials reproduce 3-4 significant digits for all atoms. The large-core counterparts are not that accurate and overestimate the energy by 5%-10%. This alone allows that to predict that small-core pseudopotentials are expected to be much more reliable in molecular studies.

V.3 Results for diatomic systems

To keep length of the paper within reasonable limits we concentrate here on two molecules - strontium hydride and barium hydride (SrH and BaH). Both of them have attracted a significant attention recently lane15 ; abe10 ; aymar12 ; gao14 ; liu16 . We present results obtained with the more reliable small-core pseudopotentials only. Analogous results for the large-core effective potentials can be obtained from the authors upon request.

For each of the molecules we evaluate the interaction energy ($D_{e}$) of the ground $X^{2}\Sigma^{+}$ state at the experimentally determined geometry. We set the interatomic distance to $R=2.1461$ and $R=2.2319$ for SrH and BaH, respectively, in accordance with the most recent experimental studies shayesteh04 ; ram13 . Additionally, we evaluate the permanent dipole moment of both molecules and their vertical ionisation energy.

The procedure for evaluation of the aforementioned quantities is as follows. The interaction energy (i.e., the well depth) is evaluated at the all-electron CCSD(T) level of theory by using the new basis sets, wtcc-$l$, with $l=2,3,4,5$. Next, valence full triples correction is added, obtained at a difference between the frozen-core full CI (FCI) and frozen-core CCSD(T) values. All results are extrapolated towards the complete basis set by using the ordinary $L^{-3}$ formula. The ionisation energy is evaluated as a difference between the extrapolated CCSD(T) energies of the molecule and the corresponding ion at a fixed geometry. Permanent dipole moments of the molecules are evaluated with the finite field method by using displaced CCSD(T) energies. In contrast to the other quantities, basis set convergence pattern for the dipole moment is not entirely regular. Thus, the extrapolation to the complete basis set has not been attempted and we simply give values calculated with the largest available basis sets.

Let us begin the analysis with the interaction energy of the barium hydride. This quantity is interesting because of a controversy connected with interpretation of the experimental data. The original experimental work of Kopp et al. kopp66 gives the value $D_{e}<16350.0$ cm*-1*. However, in a recent paper of Moore et al. lane15 a significantly larger value has been obtained from ab initio calculations, $D_{e}=16895.12$ cm*-1*. The discrepancy can be explained by assuming that the asymptote of one of the electronic states has been incorrectly identified. By selecting the correct Ba$(^{3}\mbox{D}_{3})$ asymptote instead of Ba$(^{3}\mbox{D}_{1})$, a revised experimental value is obtained, $D_{e}<16910.6$ cm*-1*. Our ab initio results are given in Table 6 and the final value, $D_{e}=16901.5$ cm*-1*, supports the revision of the experimental data. The difference between the theoretical value and the original experimental result ($\approx 550$ cm*-1*) is too large to be explained by the basis set error or the pseudopotential error. Moreover, the agreement between our result and the value of Moore et al. lane15 is striking. Note that the $\Delta\mbox{fci}$ correction is very small for this molecule, of the order of few wavenumbers, indicating that the CCSD(T) method works exceptionally well for this molecule.

Passing to the second molecule, strontium hydride, the corresponding results are given in Table 7. Unfortunately, for this system we have no direct experimental results at our disposal. However, we can compare our results with values reported in other theoretical papers. The most recent result of Liu et al. liu16 gives $D_{e}=14114.6$ cm*-1*, i.e., differing merely by 17 cm*-1* or about 0.1%. Somewhat older paper of Gao et al. gao14 gives $D_{e}=14259.8$ cm*-1* - a slightly larger deviation from our value. However, let us note that a significantly smaller basis set was used in this work. Overall, it appears that the newest theoretical values converge towards the most probable result around $D_{e}=14100$ cm*-1*. Parenthetically, the values of the $\Delta\mbox{fci}$ correction are by an order of magnitude larger for SrH than for BaH, indicating that the former possesses a much more pronounced multireference character.

Finally, in Table 8 we present vertical ionisation potentials and permanent electronic dipole moments calculated for both molecules. Unfortunately, these values are not directly comparable with any experimental data available. Nonetheless, they can be used for comparison with other theoretical results, e.g., note that the permanent dipole moments of SrH reported here are substantially larger than the values given by Gao et al. gao14

To conclude this section we would like to comment on the computational efficiency of the procedures for calculation of the pseudopotentials matrix elements. In all applications reported here we found these quantities to be much more computationally expensive than the standard one-electron integrals, both in the atomic and diatomic systems. However, this cost is still insignificant compared to the two-electron matrix elements of the electron-electron repulsion operator. Therefore, calculations of the effective core potentials matrix elements do not constitute any significant bottleneck within the present approach.

VI Conclusions

We have presented a general theory to evaluate matrix elements of effective core potentials in one-electron basis set of Slater-type orbitals. As a rule, we have used the Barnett-Coulson translation method for STOs whenever possible. It generates transparent formulae and all infinite summations truncate. As a result, the matrix elements are reduced to relatively simple one-dimensional integrals. We have presented a scheme to evaluate them to a very good precision.

Next, we have shown that the matrix elements of the spin-orbit pseudopotentials are reduced to the same basis quantities as averaged effective potentials and only minor modifications are necessary to accomplish the calculations. Somewhat larger changes are necessary to facilitate computations with the core polarisation potentials due to the apparent singularities in the potential. Additional one-dimensional integrals with logarithmic singularities appear and we have discussed their evaluation in details.

Finally, various numerical examples have been provided to verify the validity of the present approach. First, we have shown a set of test results for the calcium, strontium and barium atoms, and compared the excitation energies, dipole polarisabilities and valence orbital energies with reliable reference (exact or experimental) data. In all cases we have found a very good agreement. Lastly, we have considered two molecular systems (strontium and barium hydrides) and evaluated interaction energies, permanent dipole moments and ionisation energies; deviations from the available experimental values have been found surprisingly small.

In this paper we have concentrated mainly on the diatomic molecules. However, the present approach can probably be extended to an arbitrary polyatomic case with relative ease. This may be important for calculations in the spirit of density functional theory adf01 , but also for general quantum chemical calculations for polyatomic systems in the STOs basis in the face of recent improvements in many-centre STOs integrals technology jablczynska17 .

The code for evaluation of matrix elements of the effective core potentials in the STOs basis described in this paper has been incorporated in the Kołos program: general purpose ab initio program for electronic structure calculations with Slater-type orbitals, geminals, and Kołos-Wolniewicz functions.

Acknowledgements.

ML acknowledges the Polish Ministry of Science and Higher Education for the support through the project “Diamentowy Grant”, number DI2011 012041. AMT was supported by the Polish National Science Centre, project number 2016/21/N/ST4/03734.