A new generation of effective core potentials from correlated and spin–orbit calculations: Selected heavy elements

Most electronic structure calculations aim at valence properties such as bonding, ground and excited states, and related properties. These characteristics are determined by the valence electronic states with the energy scale of eVs and spatial ranges from covalent bonds to delocalized conduction states. On the other hand, the atomic cores of heavier elements are very strongly bonded to nuclei and spatially very localized. Therefore, the cores appear as almost rigid repulsive charge barriers around the nuclei so that, in many quantum chemical and condensed matter electronic structure calculations, the cores are routinely kept static and frozen. Building upon this understanding, it has been realized that the atomic cores can be alternatively represented by properly adjusted effective core potentials (ECPs). The effective potentials mimic the action of the core on valence states and allow for valence-only calculations with resulting gains in efficiency. On a quantitative level, this partitioning is based on significant differences between spatial and energy domains that are occupied by core vs valence states. ECPs and closely related pseudopotentials in the condensed matter context have been developed over several decades and involve a number of complete tables^1–14 as well as computational tools.^15,16

The basic construction of ECPs involves reproducing valence one-particle eigenvalues and closely related one-particle orbital norm conservation, i.e., the amount of valence charge outside an appropriate effective ion radius.¹¹ Since the number of core states and their spatial properties vary, each angular momentum symmetry channel requires a different effective potential resulting in semilocal ECPs with corresponding projectors. This simplest construction can be further generalized in several directions.^17–20 Further improvement has been introduced through the so-called energy consistency that requires reproducing energy differences such as atomic excitations within a given theory, e.g., Hartree–Fock.^7,8,11,14 Other directions for improvement have been explored within Density Functional Theory (DFT) using many-body perturbation theory.^21,22

Very recently, we have introduced correlation consistent ECPs (ccECPs)^23–26 that build upon previous constructions. Our overarching goal was to provide ccECPs that would offer the accuracy needed for many-body and highly accurate calculations. For this purpose, our construction has been based on many-body wave functions. These have enabled us to use nearly exact results to reach higher levels of accuracy and also to ascertain the robustness and transferability of the constructed potentials. The construction of ccECPs for heavier atoms involves several steps. The initial construction involves reproducing atomic excitations across a range of energies and different states including also a number of highly ionized states. We included scalar relativistic effects from the outset, and as explained below, spin–orbit effects were added as they can influence the valence accuracy for heavier elements.

We have taken into account further criteria having in mind the ultimate goal of describing with high accuracy the valence properties in real systems such as molecules or condensed matter systems. Therefore, the transferability of the ECPs has been carefully tested by probing the (effective) ion in bonded environments. In this setting, reproducing molecular bonding curves has become another important criterion both for construction and validation. Therefore, validation tests of the ccECP construction have included hydride and oxide dimers for each element. This provides insights into the simplest bonds with covalent and ionic character. Molecular calculations involved essentially the full binding curve from the stretched bond lengths to the dissociation limit at short interatomic distances so as to test the restructuring of the valence charge in high pressure bonding environments. The typical discrepancies that were observed have been within chemical accuracy (≈0.043 eV). There were a few exceptions where discrepancies were larger due to a small number of valence electrons; especially, at highly compressed bond length, we find the inclusion of norm-conservation is essential to alleviate such discrepancies. Therefore, additional cautions are needed to achieve optimal balance in several properties including spectrum and molecular bindings.

The tests and comparisons with previously generated tables have shown that ccECPs represent a significant step forward in achieving accuracy and fidelity within the ECP effective Hamiltonian model. In addition, we have found that correlated construction also provides welcome and important gains when compared with frozen core, all-electron treatments. In particular, ccECPs capture core-valence correlations that are missed in most frozen core calculations. We have found that using ccECPs provides higher valence accuracy than uncorrelated (self-consistent) cores (UC).

In this work, we extend our generation of ccECPs beyond the 3d and 4s4p elements where explicit spin–orbit interaction is an important ingredient for accurate correlated treatments. The constructed ccECPs are selected heavy elements that include 4d and 5d as well as main group atoms, namely, I, Te, Bi, Ag, Au, Pd, Ir, Mo, and W. This choice has been motivated by several considerations. One is the development of an efficient methodology for 4d and 5d transition elements that require an accurate representation of atomic spin–orbit effects. Another goal was to include elements that are prominent in a number of technologically important 2D materials.^27–30 

We mostly employ our previously developed methodology with several updates and adjustments needed due to increasing demands on the correlated treatment of large atomic cores. Intuitively, explicit spin–orbit effects are treated as further refinements and advanced features relying on the accurate spin–orbit averaged relativistic effective potential (AREP) to form the spin–orbit relativistic effective potentials (SOREP). The detailed methods for construction of these heavy element ccECPs are described in Sec. II.

The composition of this paper is as follows: Sec. II describes the form parameterizations of ccECPs. In what follows, we discuss the objective function and optimization procedure for AREP and SOREP. In Sec. III, the results are presented including the atomic properties and selected tests on molecular systems using both correlated methods based on basis set expansions such as coupled cluster with the singles, doubles, and perturbative triples [CCSD(T)] method and fixed-phase spin–orbit diffusion Monte Carlo (FPSODMC) approaches.³¹ Each element is presented in greater detail and this is followed by summaries of ccECP properties. The results are further elaborated in discussion and conclusions in Sec. IV.

The aim of this study is to accurately reproduce the properties of the relativistic all-electron (AE) Hamiltonians with a much smaller valence effective Hamiltonians H_val. Following the Born–Oppenheimer approximation, the valence Hamiltonian H_val in atomic units (a.u.) is expressed as

The full spin–orbit relativistic effective potential (SOREP) ECPs are of the form proposed by Lee,³² 

where r is the electron–ion distance and $VLJSOREP(r)$ is a local potential. This form can be split into the averaged relativistic effective core potential (AREP) and spin–orbit potential (SO) terms,³³ 

where $VℓAREP$ is the weighted J-average of V^SOREP,

The set of SO potentials V^SO are defined using V^SOREP as

where

V^AREP is defined as follows similar to our previous works:

The latter part of Eq. (7) involves the non-local |ℓm⟩⟨ℓm| spherical harmonics projectors. $VLAREP(ri)$ is again a local channel that acts on all valence electrons and parameterized as

where the Z_eff = Z − Z_core is the effective core charge and α, β, γ_k, and δ_k are optimization coefficients. With the given format of local potential, the Coulomb singularity is explicitly canceled out and first derivative at the origin vanishes. The non-local potential is expressed as

where β_ℓp and α_ℓp are parameters to be optimized. In most cases, n_ℓp were set to be n_ℓp = 2; however, we included n_ℓp = 4 terms in some cases to achieve the desired accuracy.

The ccECP SO parameters are provided as

which is also the convention adopted by codes such as dirac, molpro, and nwchem. $VlSO,ccECP$ are parameterized similar to AREP non-local potentials,

Together with parameters described in the AREP part, n_ℓp′, β_ℓp′, and α_ℓp′, are the full sets of variables to be determined and treated by the optimizer in minimizing the chosen objective functions, separately at AREP and SO levels in sequence. The optimization is in such a way that SO parameters are pursued after we have ensured the desirable AREP accuracy. AREP parameters are kept fixed during SO optimizations.

We aimed to keep the parameterization of ECP in a simple and compact form. For the fourth row elements, I, Te, Ag, Pd, and Mo, ℓ_max = 2 is used as there are 3d electrons being included in the core. For the fifth row elements, Bi, Au, Ir, and W, ℓ_max = 3 is employed because the 4f electrons are present in the core. For each channel, the chosen number of Gaussian terms varies from 2 to 4 due to different scenarios in the ECP construction of each element. Note that we did not include any core polarization potential (CPP) terms or considered these in our calculations.

The objective functions used to optimize the AREP portion of each ccECP depended on the particular element in question. The general recipe follows our previous studies.^23–25,34 We include the many-body energy-consistency and single-particle eigenvalues in the definition of the objective function,

where the $ΔEX(s)$⁠, X ∈ {ECP, AE}, denotes the atomic energy gaps referenced to the ground state for given Hamiltonians. The subset of states is chosen by picking a representative set of states for the pertinent valence space. In particular, the electron affinity (EA), neutral excitations, and various ionization levels (IPn) are included. w_s are weights corresponding to spectral states, decreasing as the energy gap increases due to the deep ionizations. $ϵX(i)$⁠, X ∈ {ECP, AE}, labels the one-particle eigenvalue in ground state for ith valence eigenvalue with a weight w_i. The ECP parameters were initialized from either MDFSTU³⁵ or BFD^7,8 ECP parameters.

All energies in the AREP case are calculated using the CCSD(T) method with large uncontracted aug-cc-p(wC)VnZ (n ∈ T, Q, 5) basis sets while adjustments are made either to add some diffuse primitives or to remove primitives that cause near-linear dependencies, as necessary. Here, AE calculations are fully correlated, which include core–core, core-valence, and valence–valence correlations and use a scalar relativistic 10th-order DKH Hamiltonian.³⁶ All ECP calculations correlate the full valence space as well (including the semi-cores if present). Molpro code³⁷ was used for all AREP calculations.

For the elements Ag, I, Mo, and W, the simple objective function given in Eq. (12) produced accurate ECPs. This objective function includes the single particle eigenvalue that is one term of the norm-conserving objective function used in previous work.²⁴ Some of the other elements in this series did not yield as accurate ECPs when this method was employed leading to two other ECP optimization schemes being used. The first of these was a modification of the objective function [Eq. (12)] by adding the full norm-conservation criteria. Refer to Ref. 24 for details. This norm-conserving method was only used for Bi.

Finally, a method similar to the single-electron fit (SEFIT), and multi-electron fit (MEFIT) optimizations developed by the Stuttgart group,¹⁰ has been employed. In this method, the ECP parameters are first optimized using single-particle energies only, and then, the exponents are kept fixed in the MEFIT optimization. This method was used for the elements Au, Ir, and Pd with the objective function given in Eq. (12) being used in the MEFIT step.

The spin–orbit coupling terms were optimized using the DIRAC³⁸ code and the Complete Open-Shell Configuration Interaction (COSCI) method. The reference AE atomic states were calculated using the exact two-component (X2C) Hamiltonian as implemented in DIRAC. The SO parameters were initialized from MDFSTU ECP values. In general, we keep the energy-consistency scheme in spin–orbit terms optimization,

where the $ΔEX(Y)$⁠, X ∈ {ECP, AE} and Y ∈ {s, m} denotes the atomic gaps using COSCI method. Here, Y labels different kinds of states included in the SO optimization. For Y = s, states with different charges referenced to the ground state were included, such as EA, IP, and IP2 similar to the AREP case. For Y = m, states with the same charges were included, particularly the ^2S+1L_J multiplet splittings with various S, L, and J values were included by referencing the lowest energy for the given electronic charge. This separation of Y ∈ {s, m} was motivated by the goal of capturing the SO effect better for small gaps. In addition, since COSCI is generally less accurate than CCSD(T), our expectation is that referencing the lowest energy within a given charge will result in a better error cancellation producing accurate multiplet gaps. This will be apparent in Sec. III where COSCI gaps are directly compared to experimental gaps.

We use three different metrics to assess the errors of the pseudoatom spectrum at the AREP level. One is the mean absolute deviation (MAD) of all considered N atomic gaps,

Another metric is the MAD of selected low-lying n gaps (LMAD),

For LMAD states, we chose electron affinity (EA), first ionization potential (IP), and second ionization potential (IP2) states only. Finally, we also consider a weighted-MAD (WMAD) of all considered N gaps as follows:

The pseudoatom spectrum errors are evaluated also for various other tabulated ECPs so as to assess the quality of our constructions. The ECPs are compared with MWBSTU,³⁹ MDFSTU,³⁵ BFD,^7,8 LANL2,⁴⁰ CRENB(S/L),^3,4 and SBKJC⁵ ECPs. In addition, to further demonstrate the improvement of ccECPs, we include uncorrelated-core (UC) calculations that are self-consistent AE calculations but with the only valence electrons correlated, and the cores are frozen with the same size of ECPs.

The summary of defined quantities, MAD, LMAD, and WMAD, for selected elements is provided in Figs. 1–3, respectively. Clearly, ccECP achieves consistent improvements overall. Detailed discussions for each element will follow in the parts below.

In addition, for each element, we show the transferability tests of all core approximations using molecular binding energies compared to the fully correlated AE case. We use the AE discrepancy of dimer binding energies as one criterion of ECP quality,

where D(r) is the binding energy at given separation r for dimer atoms.

For the SOREP pseudoatom spectrum, we adopt MAD definition defined above in Eq. (14) where the atomic gaps are chosen for the pertinent states from spin–orbit calculations. The MADs from COSCI calculations for MDFSTU and ccECPs are compared with AE values to show the optimization gains. In addition, we provide error analysis by comparing both ECP gaps to experimental values using the FPSODMC method that also provides the spin–orbit splittings.

The SOREP transferability tests are done for several mono-atomic dimers for both MDFSTU and ccECPs using AREP CCSD(T) and SOREP FPSODMC calculations, which are compared to experimental data. The binding curve is fitted to Morse potential,

where D_e labels the dissociation energy, r_e is the equilibrium bond length, and a is a fitting parameter related to vibrations,

where μ is the reduced mass of the molecule. Note that we use Morse potential for fitting all molecules both for AREP and SOREP treatments.

All optimized parameters for ccECPs are listed in Tables I and II, respectively, for the fourth row elements and the fifth row elements. Results for each element are discussed separately in Secs. III A–III J.

For iodine, the AREP atomic and molecular results are shared in Figs. 4 and 5. The construction of ccECP for iodine followed some of the adjustments used for 3s, 3p elements that are in the same group of elements. The atomic spectral properties are significantly improved when compared with previous constructions, especially in MAD and WMAD. The low-lying excitations are on par with the best available construction SBKJC while our oxide molecule shows a much improved binding curve (Fig. 5). MWBSTU shows less errors in IO dimers; however, it underbinds in short bond length IH and has larger errors in the atomic spectrum. At low IO distances, we see similar inaccuracies that we observed for phosphorus,²⁴ which are manifested by overbinding at significantly compressed bond lengths. Within the presented data, our construction represents the most consistent option with clear advantages for a broad range of valence-only calculations.

Table III gives SOREP iodine atomic excitation errors. We see that the spin–orbit splittings are reproduced very well, on par with MDFSTU results (that are very accurate to begin with) or marginally better. However, FPSODMC MAD shows an improvement over MDFSTU.

The overall accuracy is further corroborated in the calculation of the iodine dimer bonding curve that shows significant improvements already at the AREP level (Fig. 6). We see that the calculation reproduces correctly the experimental bond length unlike the previous construction that shows a bias toward smaller values. Furthermore, the two-component spinor calculations show an excellent agreement with experimental atomization energy that clearly demonstrates the quality of ccECP. Note that the inclusion of the explicit spin–orbit effect alleviates the overbinding of about 0.5 eV for both ECPs. The excellent agreement of consistent improvement between explicit spin–orbit FPSODMC and CCSD(T) results is encouraging and suggests comparable quality of the correlation description.

The atomic spectrum and molecular binding discrepancy data for Te are provided in Figs. 7 and 8. In the atomic spectrum, we observe a marginal improvement in ccECP LMAD compared with other ECPs. Regarding MAD and WMAD, ccECP shows remarkable improvement, resulting in a much better description of higher atomic excitations. The molecular data in Fig. 8 demonstrate the transferability test in TeH and TeO molecules. In TeH, we see that ccECP is well within chemical accuracy except very short bond lengths near the dissociation limit. LANL2 and CRENBS appear to be better in TeH; however, they overbind significantly in the oxide dimer TeO. In fact, we see that all other ECPs noticeably overbind in TeO through all geometries, even at the equilibrium bond length. ccECP is very accurate near the equilibrium bond length and minor overbinding results only at the shortest bond lengths. Although BFD behaves slightly better at short bond lengths, it errs significantly at the equilibrium, which is vitally important for molecular and condensed matter properties. Obviously, we achieve a better balance between these two types of bonds. A similar compromise between hydride and oxide molecules has been seen from our previous ccECP constructions in 4p elements.³⁴ 

Tellurium SOREP atomic energy gap errors are provided in Table IV. We note that the MAD of COSCI gaps for ccECP is slightly larger than for MDFSTU, and further improvements without compromising other ccECP aspects proved to be difficult. (Exploiting the greater variational freedom of additional Gaussians could provide some options in the future.) The fixed-phase DMC calculations show similar results for ccECP and MDFSTU, which compare the atomic gaps and spin–orbit splittings as referenced to experimental data. Some of the differences are larger than desired; however, this is to be expected due to small valence space, neglect of core relaxations, and for FPSODMC using a single-reference trial function. The errors are larger for higher excitations, and at this point, it is not clear how much would that affect accuracy in bonded settings. Further research might be necessary if tests in bonded systems of interest would indicate that higher ccECP accuracy might be required.

Figure 9 shows the Bi atomic spectral errors of all considered core approximations models considered. Our ccECP displays the smallest MAD and WMAD of all the approximations, while LMAD is also within chemical accuracy. Similarly, molecular errors are shown in Fig. 10 for varying bond lengths. In BiH, ccECP results in the smallest errors that are mostly within the chemical accuracy with a slight underbinding existing near the dissociation limit. For BiH, SBKJC and LANL2 ECPs show competitive errors; however, they significantly overbind in BiO molecule with up to 2 eV errors. On the other hand, ccECP errors in BiO are mostly within the chemical accuracy throughout the whole binding energy curve. Note that UC severely underbinds in both molecules and shows larger errors in the atom compared to ccECP. Overall, the data suggest that it might be possible to achieve better accuracy with proper form and optimizations even compared to AE systems with the same active space.

Another point of interest is that, in this case, the core/valence partitioning is not as clear-cut as for 5d elements that include the semi-core 5s, 5p subshells into the valence space. Specifically, this requires partitioning of n = 5 principal quantum number, where 5d¹⁰ is in the core while 5f¹⁴ could be considered semi-core or valence space. Note that this type of partitioning (i.e., the lowest one-particle eigenvalue does not correspond to ℓ = 0 channel) could result in significant errors in transition metals.¹⁴ However, the errors seen in this case are similar to what was observed in isovalent elements, such as N, P, and As,^23,24,34 with analogous core/valence definitions.

The SOREP spectral data for Bi are given in Table V. Our chosen subset of valence states is composed of EA, IP, and single d/f excitations. The inclusion of single d/f excitation states is essential to constrain the spin–orbit splitting bias for corresponding channels. Overall, we see a marginal reduction in MAD and a more noticeable improvement in the ground state and in the first IP multiplet splitting.

Figure 11 shows the bismuth dimer binding curve. The molecular calculations include AREP CCSD(T) and SOREP FPSODMC using PBE0 trial wave functions for both MDFSTU and ccECP. For both ECPs, the two-component spinor FPSODMC calculations show significant alleviation of overbinding that is present in AREP UCCSD(T). Our ccECP outperforms MDFSTU in both AREP and SOREP calculations being very close to experiments.

For silver, the averaged relativistic atomic and molecular results are shown in Figs. 12 and 13. All ECPs show quantitatively good accuracy for LMAD due to the simple closed d shell. If we consider broader states in the spectrum, MAD and WMAD reveal significant improvement achieved comparing to other ECPs. In Fig. 13, AgH is quite well described by most ECPs, and ccECP is among the best ones. Improvement is more noticeable in AgO; although most core approximations maintain the accuracy for all geometries, ccECP and MWBSTU show the best performance with almost perfect agreement with AE binding energies. Note that, in both molecules, ccECP is the closest to AE results while keeping the discrepancies flat throughout; therefore, it provides the highest accuracy of all core approximations.

The SOREP atomic and molecular results are given in Tables VI and VII. The COSCI agreement of MDFSTU and ccECP with AE gaps is remarkably good, with MAD less than 0.01 eV. The closed-shell electronic configurations lead to Ag being a special case resulting in a single determinant COSCI wave function even in spin–orbit relativistic REL-CCSD(T) calculations (Au is a similar case, too). This enables us to perform explicit spin–orbit relativistic REL-CCSD(T) calculations for the lowest charged states with results collected in Table VII. We observe the gaps errors smaller than chemical accuracy when compared to experimental data for both ECPs.