A new generation of effective core potentials from correlated calculations: 3d transition metal series

Effective core potentials (ECPs) and closely related pseudopotentials provide a well-known technique for simplifying electronic structure calculations to valence-only degrees of freedom. The ECP Hamiltonians replace the core states by potentials with projectors that mimic the action of the core on the valence electrons with different symmetries. For very heavy elements, ECPs become almost indispensable due to major complications that come from relativity, overwhelming energy scales of the core states, and difficulties in correlating cores in multi-atom systems. In fact, effective core calculations might even increase the overall accuracy for valence properties since ECPs capture effects that might be otherwise ignored. Clearly, one can adjust the ECP constructions to describe at least some of the core-core and core-valence correlation effects on the valence space as noted also previously.^1,2 In addition, scalar relativity impacts on valence can be incorporated into ECPs in a straightforward manner. Furthermore, when dealing with spin-orbit explicitly, ECPs conveniently conform to the transformation of the four-component Dirac formulation to the two-component formalism that again simplifies the subsequent calculations.

Over the past few decades, constructions of ECPs have evolved to a high degree of sophistication.^2–9 However, most of these approaches have been focused on reproducing the original all-electron (AE) Hamiltonian within some approximate method, typically Density Functional Theory (DFT) or Hartree-Fock (HF)/Dirac-Hartree-Fock (DHF). In particular, the norm/shape preservation outside an effective core radius has been used very extensively.^3,10–12 Many constructions were also designed to be efficient for codes with the plane wave basis. One of the recent focal points within the DFT framework has been the fidelity of ECPs for calculations of transition metal oxides.¹³ Important refinement in self-consistent wave function theories pioneered by Stoll and co-workers has been the adjustment of ECP atomic excitations to their all-electron HF or DHF values.^2,14,15

Further developments have targeted better description of the electron correlations and inclusion of the many-body effects into the ECP constructions. For example, correlated density matrices were employed as the key quantity to be matched outside the core region.¹⁶ Another step forward has been represented by refinement of the ECPs that reproduce atomic excitations from correlated calculations.¹⁷ 

Recently, we have advanced the use of correlated methods in ECP constructions in a systematic manner and we have proposed to increase the accuracy of ECP operators to a significantly higher level for first and second row elements.^1,18 This was motivated by the needs of accurate correlated methods such as Coupled Cluster (CC), selected Configurations Interaction (CI) combined with quantum Monte Carlo (QMC)¹⁹ and full CIQMC (FCIQMC)²⁰ that allow for correlated calculations of larger systems. Indeed, the accuracy of these calculations started to be hampered by the biases in the existing sets of ECPs^21–23 and therefore much more reliable, accurate, and tested ECPs became highly desirable.

Let us recall the key principles of our correlation consistent ECP (ccECP) constructions:

• reproducing the many-body spectrum of all-electron, relativistic, nearly-exact calculations for a subset of atomic states;

• adjusting the ECPs to capture the behavior of the original Hamiltonian in bonded situations for both equilibrium and non-equilibrium atomic geometries;

• using a simple semi-local form so that the potentials are given by a minimal or small set of gaussian functions; in addition, the potentials are bounded and smooth at the origin/nucleus; we also try to keep the tails in nonlocal channels short so as to diminish their actions in bonding regions; and

• finally, establishing a database and systematic labeling together with data on obtained accuracy benchmarks; in addition, we keep the table open to further updates, benchmarks, and refinements.

In this work, we apply these ideas to the 3d transition metal series. These atoms present their own set of challenges when compared with the first- and second-rows and therefore several considerations have to be taken into account. First, although the nominal valence levels are only 3d and 4s, accurate ECPs require the semicore 3s and 3p states to be included in the valence space as has been shown repeatedly in the past.^14,24 This is straightforward to understand considering that 3s, 3p, and 3d states occupy the same electronic shell and therefore bonding and hybridization of 3d levels have a significant impact also on the semicore levels. Second, the elements in the middle of the 3d series are key constituents in many magnetic materials and therefore require rather a precise description of related electronic structure properties. Often the Hund’s rule high-spin atomic ground state is significantly modified in the bonded environment. This can involve several bonding mechanisms such as charge transfer from 4s states and/or d-s, d-p, or d-d hybridizations with resulting partial or full quench of the local magnetic moment. Similarly, spin flip energies and a large number of states/multiplicities of the open d shell are important in many catalytic chemicals and materials. Particular importance for these changes in the electronic structure is related s ↔ d promotion energies as well as higher occupations of the d subshells in both neutral and ionized cases.

Several tables of ECPs for 3d transition metal atoms and beyond have been produced over the past three decades; for example, see Refs. 10, 13, 14, 17, and 25. In particular, the full periodic table developed over a period of time by the Stuttgart-Bonn-Koeln collaboration^2,14,15 has been used quite extensively.

We expand upon these advances here and construct high accuracy ccECPs for the elements Sc to Zn. We use a very compact form with a few gaussians per channel and adjust the terms so that the potentials are finite and smooth in the origin as established previously.^12,25 This form can, therefore, be used in many packages as well as with a variety of basis sets. In particular, we are interested in its usefulness also for periodic systems such as extended 2D materials or 3D solids calculations based on plane waves.

The all-electron and ECP calculations we use are based on the scalar relativistic CC method with large basis sets and basis extrapolations. Although the calculations are very accurate, obviously they have limitations. The important factors are extrapolations to the infinite basis set limits, level of correlations in the CC method, and scalar relativistic treatment. We estimate that these biases could sum up to ≈0.05 eV for almost all of the valence energy differences we consider. Therefore we deem ECPs with discrepancies within this threshold as being of comparable accuracy.

An additional point for future use is that we provide rather accurate values of atomic total energies for the ground and selected excited states. The high accuracy CC method with extrapolations provides estimations of the exact eigenvalues within systematic errors that vary from a few to about ≈10 mHa for the heaviest atoms. These data offer a useful checkpoint for many methods even for those that rely on sizable error cancellations in differences such as DFT. Since practical versions of DFT employ approximate functionals, this also provides valuable data for analysis of the performance in all-electron vs. ECP formulations since they can lead to non-negligible differences.

We point out that the presented ECPs include also core-core and core-valence correlations and therefore they represent effective Hamiltonians that have almost as broad use as the original all-electron one, i.e., for many purposes they aspire to be almost universal. That should be true essentially for any method and any valence property that does not require explicit core presence. We believe that this is important because it offers a well-defined platform to develop a systematic understanding of biases in a variety of approaches. For the transition metal series, this is particularly relevant since fully correlated relativistic all-electron calculations of such systems are out of reach for just a few atom systems. The ECPs, therefore, enable us to expand the system sizes that can be treated by many-body methods and at the same time provide a reference framework to study systematic biases.

In what follows, we first present the methods and objective functions used in the optimization process. Section III shows the atomic and molecular properties of the constructed operators compared with all-electron and other existing tables. We comment and discuss several aspects of these new constructions in Sec. IV.

The ECP electronic Born-Oppenheimer Hamiltonian (a.u.) has the following form:

For this work, we use a semi-local ECP form with a minimal number of parameters,²⁶ 

where r_i is the radial distance of the ith electron from the core’s origin and ℓ_max is the maximum angular momentum for the non-local projectors. In this work, we choose ℓ_max = 1. The non-local terms contain the projectors on the angular momentum states with ℓm quantum numbers. The local term, V_loc, is chosen to cancel the Coulomb singularity at the origin, i.e.,

where Z_eff is the effective core charge, Z_eff = Z − Z_core. To cancel the Coulomb singularity at the origin and force the first derivative to be zero at the origin, we fix β_loc,1 = Z_eff, β_loc,2 = Z_effα_loc,1 with n_loc,1 = 1 and n_loc,2 = 3 and all other n_loc,k = 2. The V_ℓ(r) potentials are expanded as follows:

where all the n_ℓ,k are fixed to n_ℓ,k = 2. All variables are treated as optimization parameters in the minimization of a chosen objective function, except for the terms we have fixed as described above. In addition, a constraint that imposes a concave shape of the potential at the origin is imposed,²⁶ 

The fully correlated calculations of transition metals and heavier elements with very large basis sets and high accuracy methods such as CCSD(T) are very costly and, eventually, impractical. Therefore, we constructed an optimization strategy that introduces correlation into the ECP without explicitly calculating the correlation energies at each step of the optimization. This approach relies on the fact that the correlation energy for a given atomic state is similar across different ECPs, and the variation for the correlation energy error is significantly smaller than the variation in the HF error, as shown in Fig. 1. For a particular excitation energy, the all-electron (AE) contribution can be written as $ΔgAE=ΔHFAE+ΔCCAE+ΔCVAE+ΔVVAE$⁠, where CC, CV, and VV represent the core-core, core-valence, and valence-valence correlation contribution to the excitation energy, respectively. When using an ECP, there are only valence electrons, so the ECP excitation energy can be written as $ΔgECP=ΔHFECP+ΔcorrECP$⁠. In order to match the AE excitation energy, the ECP can either adjust $ΔHFECP$ or $ΔcorrECP$⁠. By considering a variety of ECPs, we plot the spread of HF errors and the difference between the ECP correlation energy and the valence-valence correlation energy of the AE in Fig. 1 for a variety of ECPs, including Burkatzki-Filippi-Dolg (BFD) ECP,²⁶ Stuttgart (STU) ECP,¹⁴ and Trail-Needs ECP (eCEPP),¹⁷ and a variety of our constructed ECPs at both the HF and correlated levels. The spread of HF errors is larger than the correlation energies, which indicates that the HF contribution is more flexible across ECP parameterizations whereas the correlation energy varies much less. This is true whether the ECP form shows smooth and finite amplitudes at the nucleus or whether it contains −Z_eff/r or 1/r² divergences. This suggests that the $ΔHFECP$ can be adjusted while treating $ΔcorrECP$ as a constant during an optimization procedure, which will result in $ΔgECP=ΔgAE$ as is desired. This fact has been used in the optimization strategy below.

1. For a given atom, a set of all-electron (AE) atomic states is calculated using the CCSD(T) method with the uncontracted aug-cc-pwCVnZ-DK²⁷ basis set and extrapolated to the complete basis set (CBS) limit using data for n = T, Q, 5. We utilize the Douglass-Kroll-Hess (DKH) Hamiltonian to incorporate scalar relativistic effects for all atomic states to generate our AE reference data. For the HF, we extrapolate to the CBS limit with

where n labels the basis set size, $ECBSHF$ is the CBS limit, and a and b are other fitting parameters, and

where the correlation energy is defined as $Encorr=EnCCSD(T)−EnHF$⁠. It is known that the atomic states of some transition metals such as the ground state of V are incorrectly described by a single reference using real spherical harmonics. In order to obtain the correct atomic symmetries, we generate “symmetry equivalenced” orbitals by state-averaging the irreducible representations (irreps) of D_2h for a particular atomic state.²⁸ We use the resulting natural orbitals to calculate the correlation energies with CCSD(T). This way, we are able to obtain reference AE energy gaps from a consistent description of the atomic state.

1. The initial construction of the ECP is launched by reproducing the scalar relativistic Dirac-Hartree-Fock (DHF) spectrum. We generate a large set of random ECPs and optimize each using the nonlinear DONLP2 code by Spelluci.²⁹ The objective function, Σ, to be minimized over the spectrum S is given by

where ΔE^(s) is the excitation energy for state s relative to the neutral ground state and $ws$ is the weight assigned to that excitation. The spectrum is optimized using numerical DHF code³⁰ to avoid basis set errors and to speed up the calculations. The set of states considered includes the neutral and singly ionized s ↔ d excitations, as well as further ionizations and the anion ground state. These states are chosen due to the fact that the low-lying s ↔ d transitions and various oxidation states play a crucial role in transition metal chemistry.¹⁴ The resulting ECPs are further refined as follows.The DHF spectrum optimized ECP with the lowest objective function in step 2 is used to calculate the correlation energies using CCSD(T) for the states mentioned. This defines $ΔcorrECP$⁠, i.e., the correlation energy contribution to the gaps. For the ECPs, we utilize the same uncontracted aug-cc-pwCVnZ-DK basis set. It should be noted that for the ECP calculations, no DKH Hamiltonian is used since we intend to build the scalar relativistic effects directly into the effective potential.

1. It is clear that DHF spectrum optimized ECP atomic excitation energies $ΔEECP(s)$ are not going to match the corresponding AE excitation energies $ΔEAE(s)$ when a post-HF method such as CCSD(T) or CI is used. However, it is possible to obtain $ΔgECP=ΔgAE$ by only matching $ΔHFECP$ to shifted $ΔHFAE$⁠,

where ϵ is given as

Here AE correlation energy refers to the total correlation energy (CC+CV+VV) of the state. The new shifted gaps are inserted into the objective function in Eq. (8) to be minimized.

1. The steps 3-4 are iterated until a self-consistent ECP is obtained. For each step, ϵ is re-evaluated using the ECP from the last iteration. In every iteration, the ECP parameters are randomly perturbed within ≈1%-2% of parameter values to ensure better scanning of possible local minimas around the current values of the parameters. Usually, a set of self-consistent parameters are obtained within a few iterations.

In cases where pure energy consistency does not result in an acceptably transferable ECP (as described in Sec. III), we reduce the spectrum by removing very highly ionized states and add additional constraints to improve the transferability. Note that the full spectral fits include also core-core and core-valence correlations that become quite significant for the highly ionized states and therefore they have the tendency to steer the ECP from the optimal valence-valence description. This is especially true for our ECP form with small number of variational parameters. Although the bias is typically small (say, 0.1 eV), we opted for further refinements. This is accomplished by including one-particle eigenvalue discrepancies into the objective function. This increases the HF character of the atomic states. In these cases, we utilize a new objective function, Γ, given as

where the first term is the initial spectral objective function and ϵ_ℓ is one particle eigenvalue. The eigenvalues are weighted by γ to allow the spectrum to be minimized while keeping the one-particle eigenvalues reasonably close to the corresponding AE ones. With few iterations, one can find a compromise that reproduces a large part of the spectrum as well achieves transferability in molecular calculations. Clearly, further improvements are possible, but in this work we are mostly concerned with the demonstration of principle and with providing simple ECPs that fulfill the accuracy criteria.

We note that the objective function utilized here for the 3d transition metals is different than what was used for the first and second row elements.^1,18 As we have seen, energy consistency alone [our spectral-only ECP (ccECP.S)] does not in general result in an ECP with accurate atomic properties as well being transferable to molecular systems. In order to improve transferability, we include HF properties into our objective function as we did in the second row. Due to the large number of valence electrons for the transition metals, we found that our objective functions only needed a modest level of HF information to improve overall transferability as opposed to the second row, where the smaller number of electrons and large cores required a more complicated construction.

We present results comparing our correlation consistent ECPs, labeled ccECP, to various other core approximations relative to the all-electron calculation for both the atomic spectra and monohydrides and monoxides. The all-electron reference in all following calculations is spin-restricted CCSD(T) where we have correlated all electrons in the system, i.e., no core orbitals have been frozen. To include relativity, we utilize a 10th order Douglass-Kroll-Hess Hamiltonian.³¹ To fully correlate the core, we use an uncontracted aug-cc-pwCVnZ-dk basis²⁷ and extrapolate to the complete basis set limit to minimize the basis set error. For the ECP calculations, we perform CCSD(T) calculations without DKH (since the effects are built into the effective potential) with the same aug-cc-pwCVnZ-dk basis sets and extrapolate to the complete basis set limit. We compare our ECPs against an all-electron uncorrelated core (UC), which does not allow for any excitations from the neon core in the CCSD(T), as well as other ECPs. The uncorrelated core approximation, labeled UC, only accounts for valence-valence correlation, and any core-core or core-valence correlation is absent despite it being an all-electron calculation. The other ECPs include Burkatzki-Filippi-Dolg (BFD) ECPs²⁶ which are energy-consistent DHF ECPs designed for use in QMC, Stuttgart (STU) ECPs¹⁴ which are energy-consistent DHF ECPs, and the recently constructed Trail-Needs ECPs (eCEPP)¹⁷ which use a shape and energy consistent scheme to construct correlated ECPs.

Obviously, the methods we use have their accuracy limits. In particular, we rely on the scalar relativistic approach with averaged spin-orbit effects. This introduces an additional bias of ≈0.025 eV as a representative value that we estimated from accurate spin-orbit correlated calculations of atomic excitations.³² Further important sources of bias are the sizes of basis sets (alleviated by extrapolations) and level of correlation captured by the CCSD(T) method. On smaller systems where we were able to push the basis set limits up to n = 6, 7, i.e., 6Z, 7Z and level of correlation up to CCSDT(Q) level (triples explicitly and quadruples perturbationally). These limited calculations suggest that it is difficult to ascertain accuracies better than roughly 0.02 eV for energy differences. In addition, for the heaviest atoms, we were not able to carry out the full sequence of basis set calculations up to 5Z, so this threshold is probably mildly higher. Clearly, ECPs with discrepancies below these inherent biases of ≈0.05 eV (that include also a number of our constructions presented below) are certainly interesting for methodology reasons. However, the ECPs within such bound should be considered as being of comparable quality.

We briefly discuss the results for every element in this section. In Sec. III A, we present an overview of the atomic spectra for all transition metals. In Secs. III B–III K, for every atom we show the atomic spectrum discrepancies and transferability tests on monohydrides and monoxides. In these molecules, we test various core approximations near and out of equilibrium geometries, i.e., compressed bond lengths. In general, we show the results for our spectral-only ECP (labeled ccECP.S) as well as the ECP that compromises a small part of the spectrum for an increase in the overall transferability (this is labeled as ccECP). For oxygen ECP, we use the ECP from the corresponding table; our oxygen ccECP can be found in our previous work.¹ Although for some of the oxide molecules, we were unable to converge with 5Z basis sets due to technical difficulties. However, we note that in terms of the binding curves, discrepancies are well converged at the VQZ level. Therefore, for each molecule, we show only the largest basis set that we converged as opposed to the CBS extrapolation.

All ECP parameters are given in Tables I and II while our hydrogen ECP is given in Table III.

We first provide an overview of the atomic spectra. The optimized spectrum for ccECP.S included the neutral s²dⁿ, s¹dⁿ⁺¹, and dⁿ⁺², the first ionized s¹dⁿ and dⁿ⁺¹, all ionizations down to the Ar core, and the [Ne]3s² state. As we will see in Secs. III B–III L, the inclusion of these highly ionized states can lead to non-negligible biases and a decrease in the overall transferability. Therefore, in cases where this full spectrum is insufficient, we add in constraints and work with a reduced spectrum as described in Sec. II B.

Figure 2 shows the errors across the entire transition metal series for select atomic states ranging from some low lying excited states to a highly ionized state. In Figs. 2(a) and 2(b), we show the neutral s ↔ d transitions from the s²dⁿ to s¹dⁿ⁺¹ and dⁿ⁺², respectively. Both STU and BFD result in discrepancies on the order of 0.1 eV throughout the entire transition metal series for these low lying atomic states. UC, eCEPP, and our ccECP.S and ccECP all result in discrepancies within half of the chemical accuracy, indicated by the shaded region. Figure 2(c) shows the s²dⁿ to s¹dⁿ ionization, which is well described with all core approximations. It is important to note that STU and BFD include many of these states directly into their optimization as well. However, the neglect of directly introducing correlation into the construction can result in large errors in the atomic spectra, despite being well optimized at the HF/DHF level as illustrated by both the eCEPP and our ccECPs.

As the ionization level is increased, all previous core approximations begin to show larger discrepancies with the AE reference. In order to illustrate this, we show the ionization down to the Ar core in Fig. 2(d). The STU and BFD ECPs we see show significantly larger discrepancies, ranging from 1 to 10 eV throughout the series. The correlated eCEPPs show an improvement but still result in discrepancies of a few eV. The all-electron UC approximation results in a slightly increasing error as the atomic number increases up to roughly 1 eV for Zn. Our energy consistent ECPs (ccECP.S), however, are able to maintain a uniform accuracy within half of the chemical accuracy for each atomic state and for all ionizations.

In Fig. 3, we show the mean absolute deviation,

across all of the states including the low-lying s ↔ d transitions, anions, and ionizations down to the semicore [Ar]. Our ccECP.S ECPs are able to maintain MADs of nearly 0.01 eV for all transition atoms. In Secs. III B–III K focused on the individual atoms, we investigate how energy consistency influences the transferability of the ECPs for each element. Our adjusted ccECP constructions, which slightly compromise the atomic spectrum when compared to the ccECP.S (as illustrated in Fig. 3), still maintain higher accuracy for the MADs when compared against the other tested ECP constructions. In Table IV, we summarize the MADs for each atom and ECP. LMAD indicates we are only considering the low lying atomic spectrum, which includes the electron affinity, ionization potential, and the neutral and singly ionized s ↔ d transitions. The MAD corresponds to the LMAD, as well as additional ionizations down to the Ar core. Details about the atomic excitation energies can be found in the supplementary material.

In the case of Sc, the atomic and molecular data are given in Table IV and Fig. 4, respectively. Our spectral-only ccECP.S results in a significant improvement over the atomic spectra from previously existing tables such as STU and BFD. In fact, when comparing to the all-electron uncorrelated core UC, we see an improvement for the low-lying spectrum (indicated by LMAD) as well as the full spectrum ionized down to the Ar semi-core. Compared to the recently derived eCEPP,¹⁷ our spectral errors are only slightly larger, presumably due to our significantly smaller number of gaussian functions as well as a shorter radial range of nonlocal channels. Since for Sc the higher order effects such as spin-orbit can reach 0.010 eV,³² discrepancies below this level do not provide a genuine quality measure as we have noted also above.

For the molecular calculations, we plot discrepancies from the AE CCSD(T) binding curve for a range of bond lengths. The compressed bond regime is plotted near the dissociation threshold, i.e., where the binding energy approaches zero while the vertical line indicates the equilibrium bond length. For ScH, we see that all ECPs result in binding discrepancies that are well within chemical accuracy with regard to the AE results, indicated by the shaded region. Note that STU, BFD, and UC discrepancies vary as functions of the bond length. That might cause some shifts in the predicted vibrational frequencies, whereas a flat discrepancy should lead to better accuracy for vibrational properties. The ScO molecule probes the charge transfer and polar bond regime showing that both STU and BFD exhibit significantly larger deviations from the AE potential energy surface. Our ccECP.S results in a relatively flat discrepancy throughout the entire binding region and a very small error for the binding energy. Due to the quality of the constructed operator, we decided that no further refinements to the ccECP.S were needed. The parameters for this optimal ccECP.S are given in Table I.

In the case of Ti, the atomic and molecular data are given in Table IV and Fig. 5, respectively. For the atomic spectrum, our ccECP.S outperforms all other core approximations, for both the low lying atomic spectra (LMAD) and all ionizations down to the Ar semi-core. We see that STU and BFD result in comparable MADs of roughly 0.1 eV. For the low-lying spectra, both the eCEPP and ccECP.S have marginally smaller discrepancies than UC showing that some of the core-valence correlations have been captured. Note that our ccECP.S maintains uniform accuracy throughout the entire spectrum, resulting in a MAD of 0.013 eV; see the supplementary material for details about individual states.

When considering the hydride molecule, the BFD shows significant underbinding for the compressed geometries. On the other hand, STU maintains uniform accuracy throughout the entire binding region with mild 0.04 eV underbinding. UC, eCEPP, and our ccECP.S are all very comparable for TiH, showing deviations within 0.02 eV for all bond lengths. For TiO with a polar bond, we find that BFD is inadequate to describe the binding even at equilibrium and the error reaches up to 0.4 eV in the short bond length region. STU is also underbound near equilibrium and has the opposite behavior to BFD in that the error decreases as the bond is compressed. While UC and eCEPP are well within chemical accuracy near equilibrium, each begins to underbind as the bond length is compressed. Near dissociation, both UC and eCEPP reach the region outside the chemical accuracy. Our ccECP.S, on the other hand, is well within the desired error margin throughout the entire bonding region. Additionally, the ccECP.S has an extremely flat discrepancy near the equilibrium. Considering the accuracy of the constructed ccECP.S, we did not pursue any refinements. The parameters for the Ti ccECP.S are given in Table I.

The atomic and molecular data for the V ECPs are given in Table IV and Fig. 6, respectively. For the atomic spectra, our ccECP.S has significantly smaller discrepancies compared to all other core approximations, including the all-electron UC approximation, for both the low-lying spectra and all ionizations. If we consider VO for our ccECP.S, we see that we have quite favorable binding properties compared to other core approximations. Both BFD and STU exhibit relatively large errors throughout the entire bonding region, including near equilibrium. Both eCEPP and UC begin to underbind to nearly 0.1 eV near the dissociation threshold, whereas our ccECP.S is well within chemical accuracy. However, when considering the hydride, we found that our spectral-only ECP overbinds by a small amount of ≈0.04 eV that reaches the borderline at dissociation. In fact, this is still acceptable considering the systematic biases present although less accurate than the eCEPP binding curve. Note that ccECP.S ECPs fit a significant part of the atomic spectrum that includes ionizations past the valence electrons and down into the 3s and 3p semi-core states. However, this nudges the ECP operator away from transferability optimum as can be observed on the VH molecule. Therefore, we refined ccECP.S by reducing the considered spectrum and modifying our objective function as described in Sec. II B, resulting in a refined ccECP. The charge-transfer physics is unchanged and it is almost identical to the spectral ccECP.S. However, we see a significant improvement for the hydride molecule. In terms of the spectrum, we slightly compromise the low-lying spectrum (LMAD) while maintaining roughly the same overall MAD. The parameters for our ccECP of V are included in Table I.

The atomic and molecular data for Cr are given in Table IV and Fig. 7, respectively. Our spectral ECP, ccECP.S, has a significantly improved spectrum when compared to all other core approximations, including the all-electron UC approximation. The low-lying spectrum (LMAD) is slightly improved over eCEPP, while the MAD across the spectrum above [Ar] is significantly lower. The ccECP.S performs reasonably well for CrH, being within chemical accuracy to AE for the entire binding region. However, when we consider the oxide CrO, we see that the ccECP.S shows a clear tendency to overbind, up to ≈0.08 eV near the dissociation threshold. Our refined ccECP, optimized using Eq. (11), shows a well-balanced atomic spectrum with a mild increase in LMAD. In molecular calculations, we see a modest improvement for CrH and a very satisfactory reduction in overbinding when compared to ccECP.S. The final ccECP is well within chemical accuracy throughout the entire binding region, including near the dissociation threshold. The parameters for our ccECP are given in Table I.

The atomic and molecular results for Mn are given in Table IV and Fig. 8, respectively. For the spectral ccECP.S, we find similar uniform accuracy for both the low-lying atomic and full spectrum to [Ar]. It is interesting to see that despite this fidelity, we find the overall transferability to be lacking. While ccECP.S for MnH is within chemical accuracy for most of the binding curve, we observe a slightly larger discrepancy near the dissociation threshold. This error is more severe for MnO with overbinding by roughly 0.04 eV near equilibrium that increases to 0.11 eV at dissociation. Therefore, despite having the best atomic spectrum, the ccECP.S MnO molecule has the largest discrepancies from AE when compared to all other core approximations. We clearly observe that within the used parameterization, the atomic accuracy and molecular accuracy pull in rather opposite directions.

We improve upon our spectral only ccECP.S and reoptimize via Eq. (11) to obtain our final ccECP. In terms of the atomic spectrum, we slightly increase the LMAD to ≈0.01 eV. Some error of 0.1 eV is introduced to the highly ionized [Ar] semi-core excitation. This increase results in a MAD on the entire spectrum of only 0.024 eV, which outperforms all other core approximations including the all-electron UC. For both the hydride and oxide, we see a clear improvement. For MnH, the ccECP has a discrepancy of less than 0.01 eV for all bond lengths. The most drastic improvement comes for the oxide molecule, where the discrepancy is bounded by the chemical accuracy throughout the entire curve. The resulting ccECP has excellent atomic and molecular properties and the parameters are given in Table I.

The atomic discrepancies for Fe atoms are given in Table IV. Our spectral ccECP.S produces an atomic spectrum with uniform accuracy throughout the entire valence spectrum, with a low-lying spectrum MAD of 7 meV and 5 meV for the entire spectrum. In terms of the low-lying spectrum (LMAD), the ccECP.S is comparable to the all-electron UC approximation. In addition, the ccECP.S has a much lower MAD for the entire spectrum. While the ccECP.S is well within chemical accuracy for FeH, when we consider the polar bond, illustrated by the FeO molecule, we find significant overbinding across the entire binding region (Fig. 9). In fact, the ccECP.S overbinds well outside chemical accuracy for most of the bond lengths. It should be noted that complications related to molecules with high-spin TM atoms are clearly visible for BFD and STU ECPs with significant discrepancies in the binding energies as well as in the shapes of the potential energy surfaces.

We achieve a remarkable balance between the atomic spectrum as well as the molecular properties with our refined ccECP. The ccECP results in an increase of only 3 meV on the LMAD, and we find larger errors only for the highest ionizations. It should be noted that the absolute error of 0.43 eV for the ionization to the [Ar] semi-core corresponds to a relative error of only ≈0.08% since the excitation is sizeable, 560 eV. Additionally, despite the compromise in the spectrum, the ccECP outperforms all other ECPs on the LMAD and all the core approximations on the entire MAD. The molecular properties are significantly improved for the ccECP. For FeH, there is near perfect agreement on the binding curve, both near equilibrium and near dissociation. For FeO, we see a significant reduction in the overbinding, putting the binding discrepancy within the chemical accuracy of the AE CCSD(T) binding curve. The parameters for this ccECP are given in Table II.

The atomic and molecular data for Co are given in Table IV and Fig. 10, respectively. Our spectral construction ccECP.S shows significant improvements when compared with other core approximations. We are able to achieve for both the LMAD and MAD significantly lower discrepancies than UC, BFD, and STU. In terms of CoH, the ccECP.S construction looks very satisfactory. CoH results in a discrepancy that is well within chemical accuracy for ccECP.S across the entire binding region. However, for CoO, we see that the spectral construction is less accurate and as in previous cases it results in significant overbinding. Note that other available ECPs have a similar behavior, resulting in overbinding over 0.1 eV near the dissociation threshold.

When we introduce additional constraints into our objective function, see Sec. II B, we are able to obtain a dramatic improvement in the overall transferability. With regards to the spectrum, our ccECP is comparable to the UC approximation for the low-lying spectrum (LMAD) and even slightly better for the entire atomic spectrum (signaling a complicated landscape of the objective function with many minimas). For CoH, the binding energy discrepancy is comparable to our ccECP.S, but the overall curve is flatter. The most profound improvement comes with the CoO polar bond, where we are able to obtain a nearly flat discrepancy throughout the entire binding region and very marginal 0.01 eV overbinding. Note the overall major boost in accuracy when compared to the other ECPs we studied. The parameters for the ccECP are given in Table II.

We show the atomic and molecular data for Ni in Table IV and Fig. 11, respectively. Optimization for ccECP.S results in an accurate spectrum with a 3 meV LMAD for the low-lying spectrum and a 10 meV MAD for the entire spectrum. These discrepancies are significantly smaller than for the other core approximations of UC, BFD, and STU. For NiH, we see significant issues with BFD, whereas all other core-approximations agree with AE to within chemical accuracy. We find more sizeable errors for NiO for BFD, STU, and also for our spectral ccECP.S.

To obtain a more reasonable description of the molecular properties, we reoptimize using Eq. (11) that produces our ccECP. With refinement, we compromise on both the LMAD and MAD that mildly increase relative to ccECP.S; however, this increase is caused mainly by the very deeply lying states (>500 eV) while the low-lying spectrum is still exceptionally well described, much better than for BFD, STU, and comparable to UC. With ccECP, we see a minor improvement for NiH and a significant improvement for NiO, where we reduce the discrepancy to be well within chemical accuracy for all considered bond lengths and maintain a flat discrepancy throughout. The final parameters for our ccECP are given in Table II.

We show the atomic properties for Cu in Table IV and the molecular properties for CuH and CuO in Fig. 12. Due to the closed d-shell, LMAD includes 4 excitations, whereas previous tables for Sc–Ni included 5. Focusing first on the spectral ccECP.S, we are able to obtain an almost uniform agreement for the entire spectrum, with a LMAD of 0.018 eV and a MAD of 0.015 eV. While the LMAD is slightly higher than UC and eCEPP, the overall MAD is significantly improved since the high ionizations are in much better agreement to AE. As has been seen for the previous elements, uniform accuracy on the atomic spectrum can decrease the overall transferability; see the supplementary material for details about each state in the spectrum. While CuH is quite well described by the ccECP.S, we again find significant overbinding near dissociation for our ccECP.S. In fact, all ECPs including STU, BFD, and eCEPP significantly overbind for the CuO molecule, from roughly 0.1 eV for our ccECP.S to over 0.25 eV for BFD.

We substantially improve the overall transferability with our refined ccECP. For both CuH and CuO, we have near perfect agreement with AE along the entire binding curve. For CuO, we are able to obtain the accuracy that stands out when compared to previous ECPs. Note that behavior of biases for this element suggests that keeping the accuracy for energies >500 eV, one would need more variational freedom and perhaps a more sophisticated construction. This is due to different energy and also length scales for semicore vs valence spaces.

The atomic and molecular data for Zn are given in Table IV and Fig. 13, respectively. For our spectral ccECP.S, we are able to achieve uniform accuracy across the entire atomic spectrum, resulting in a MAD of 0.015 eV. This outperforms BFD, STU, and UC significantly, resulting in an ECP that on the spectrum is very close to the original AE Hamiltonian. Additionally, the LMAD for Zn includes only the first two states due to the fact that the anion and most of the s ↔ d transitions do not exist. For the transferability tests, we find similar behavior for both ZnH and ZnO, where the ccECP.S is within chemical accuracy near equilibrium but begins to overbind as the bond is compressed.

With the refinements from our ccECP, we are able to strike a reasonable balance between the atomic properties as well as the molecular transferability. In both the hydride and oxide, the ccECP has a relatively flat discrepancy resulting in excellent vibrational frequencies as well as energy differences well within chemical accuracy for the entire curve. In terms of the atomic spectrum, the lowest states are only 0.01 eV error and we achieve a MAD on the entire spectrum that is better than all other core approximations. The parameters for our final ccECP are given in Table II.

In Table V, we collect the results of molecular calculations for all elements and evaluate mean absolute deviations for the equilibrium molecular parameters. Note that even our ccECP.S set achieves very respectable accuracy that is comparable or better than the available tables. Clearly, the adjusted set ccECP is the most balanced overall and to the best of our knowledge shows overall the best consistency for all the calculated parameters.

In this paper, we present the correlation consistent ECPs for 3d elements. First, we constructed the ECPs using spectral optimizations for all the elements. This optimization was iterative and included states that ionize the given atom down to the [Ne]3s² ion, as well as s ↔ d transfer energies for low-lying atomic and ionic states, including bounded anions. Quite unexpectedly, for most atoms, we were able to optimize the states within 0.01 eV or even smaller discrepancies for the full span of ionization energies, sometimes with accuracies better than the known spectroscopic data and also better than the inherent systematic biases in our methodology. As discussed above, any discrepancies smaller than ≈0.05 eV are in fact comparable to estimated systematic errors in our calculations. These spectral-only results simply pointed out that such fits are indeed feasible using a rather small number of free ECP parameters.

We further probed for accuracy and transferability of the constructed operators on transition metal-hydride and oxide molecules. The hydrides appear mostly less problematic and within the error bar window of chemical accuracy for all bond lengths with a few exceptions as analyzed above. The oxide molecules and their polar bonds revealed a more complicated picture and most ECPs required refinement that has led to high fidelity and chemical accuracy along the binding curves. That set is labeled as the recommended ccECP. The spectral-only ccECP.S versions could still be useful for atomic calculations where the accuracy of highly ionized states would be of crucial interest. We note that in cases of Sc and Ti, ccECP = ccECP.S since the spectral optimization provided desired accuracy without further refinement. It is quite remarkable that such accurate fits can be constructed using just the spectral information as an input.

For high spin elements, the spectral-only ECPs, while still respectable, were less accurate than desired. For stretched and equilibrium bond length regions, the binding properties were described very well. However for shorter bond lengths, we observed overbinding in some cases on the order of $≈$ 0.1 eV. Although rather small in a relative sense, we have opted for further improvement that took into account the decreasing importance of correlations of deeply ionized states vs. the key desired accuracy of the low-lying excitations most relevant for valence properties. Adding a shift toward the HF eigenvalues into the objective function enabled us to mildly retune the 3s and 3p channels and that proved to be sufficient to get the binding curves within the chemical accuracy threshold. This caused some increase in spectral errors, however, only for very highly ionized states that are of the order of >500 eV with relative errors being still very small (≈0.08% or smaller).

For the late transition elements, namely, Co, Ni, Cu, and Zn, we have observed similar behavior with an overbinding tendency for oxide molecules at short bond lengths. Similar refinement as for the high-spin elements has enabled us to alleviate these deviations so that all the binding curves were within the chemical accuracy threshold. For these elements, we also note additional complications with the largest basis sets for Co, Cu, and Zn where we have encountered the feasibility limits of the used codes. Therefore we have restricted ourselves to the best results we could obtain with reasonable computational resources and we consider the achieved systematic comparisons as adequate for our present purposes.

Together with accuracy tests, we have also tried to estimate the exact atomic total energies for both the ground and excited states. The extrapolations from extensive basis sets aug-cc-pCVnZ with n = 2, 3, 4, 5 provide rather a consistent picture of the corresponding correlation energies; however, at present, we do not claim better systematic accuracy than ≈10 mHa. In order to provide these energies with better uncertainties, we plan further investigation in a subsequent study that will be devoted to decrease the error margins to much lower values.

An observation on the legacy ECP constructions is that our results confirm that the Stuttgart-Koeln-Bonn ECPs^2,14 established and derived by Stoll, Dolg, and co-workers over the past three decades show systematic consistency and respectable accuracy for the 3d series. The discrepancies are mostly within the 0.1–0.2 eV margins for low-lying atomic excitations and also for most of the molecular binding curves. Since the STU table was constructed with Dirac-Hartree-Fock spectral fits, this confirms that the precise self-consistent energy difference data are the dominant factor in achieving consistent behavior. However, correlation effects become important at finer resolutions that we have targeted in this work.

Clearly, our present constructions raise the bar of accuracy higher and we believe that they provide a significant step forward for correlated calculations. Further data including the basis sets, including both cc-pVnZ and aug-cc-pVnZ for n ∈ {D, T, Q, 5}, can be found in the web library http://pseudopotentiallibrary.org as well as in the supplementary material. Data for this work are also provided at the Materials Data Facility.³³ 

The presented results show that there is room for sizeable improvements of both the construction methods and practical versions of ECPs. In agreement with previous papers on first and second rows, we were aiming at offering a “minimal model” that is still more accurate than existing tables. For example, nonlocal s and p channels are described by two gaussians only (similar to the STU table). We have demonstrated that such a combination of small variational space and accuracy is indeed feasible and the resulting constructions are ready for general use. There are several directions where this work could expand further. Heavier atoms such as Co, Cu, and Zn could benefit from more variational parameters that would address differences between deeply lying 3s, 3p states vs large the valence subshell. In addition, more testing and validation is needed in a variety of chemical systems such as larger molecules, solids, and 2D materials. Indeed, we expect that such data would provide new insights and possibly point out the directions for further refinements and updates in future.

See supplementary material to find additional information about our ccECPs. We provide a detailed description for the constructed basis sets that are provided at http://pseudopotentiallibrary.org. Additionally we provide information about the atomic spectra, molecular fitting parameters, and information on the objective function utilized for each ccECP.

We are grateful to P. R. C. Kent for reading the manuscript and for comments. This work has been supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, as part of the Computational Materials Sciences Program and Center for Predictive Simulation of Functional Materials.

The calculations were performed mostly at Sandia National Laboratories and at TACC under XSEDE.

Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia LLC, a wholly owned subsidiary of Honeywell International, Inc. for the U.S. Department of Energy’s National Nuclear Security Administration under Contract No. DE-NA0003525.