A new generation of effective core potentials for correlated calculations

In recent decades, many-body electronic structure methods have enjoyed steady progression in accuracy and efficiency; this has largely been driven by algorithmic improvements coupled with the increasing power of computational resources. These advances can be readily seen, for instance, in methods that rely on wave function expansions, such as coupled cluster (CC) methods, as well as stochastic approaches that sample either particle configurations or the Fock space of antisymmetric wave functions. For example, quantum Monte Carlo (QMC) methods that utilize explicitly correlated many-body wave functions often reach accuracies of 1-2 kcal/mol for systems containing elements from several rows of the periodic table.^1,2 For practical reasons, these methods often approximate properties of atomic cores, especially when calculating systems with heavier elements. In particular, the cores are either left uncorrelated, fully frozen, or even eliminated entirely using pseudopotentials or effective core potentials (ECPs).³ Note, however, that ECPs are generated within approximate theories as they are often tuned to match Density Functional Theory (DFT) or Hartree-Fock (HF) all-electron (AE) atomic properties in the valence space.

In correlated treatments of very heavy elements, the use of ECPs is essentially unavoidable since explicit correlation of atomic cores would become another bottleneck for such calculations. For instance, in QMC, the corresponding computational cost grows as $O(Z5–6)$⁠, where Z is the nuclear charge. Although ECPs provide an important efficiency boost, they inevitably introduce errors that can compromise the overall quality of results. Furthermore, the increasing accuracy of the many-body methods mentioned earlier reveals that ECP imperfections can become a dominating source of systematic bias.^4,5

In recent decades, a number of ECP tables as well as on-demand computational packages have been established using several types of electronic structure approaches. This includes DFT constructions such as the comprehensive table of Bachelet et al.⁶ based on relativistic DFT and ECPs generated by the OPIUM package.⁷ We also mention recent advances in DFT constructions with multiple projectors.⁸ Several sets of tables suitable for many-body calculations have been developed using Hartree-Fock/Dirac-Fock (HF/DF) methods.^9–16 

Most of the existing ECP constructions have relied on reproducing one-particle properties such as norm/charge-conservation⁶ or closely related shape consistency,^11,12,17 HF/Dirac-Fock (DF) excitation energies⁹ or their extensions⁸ including a recent composite scheme designed to accurately build-in many-body effects.^18–20 Comprehensive reviews^10,21 contain details as well as many relevant references which we skip here. To the best of our knowledge, the incorporation of many-body approaches into ECP constructions have been rather limited so far, first by the QMC work of Acioli and Ceperley,²² also the studies of Maron and Teichteil²³ and Fromager et al.,²⁴ and recent studies of Trail and Needs.^18–20 Overall, ECPs have arguably been some of the most successful ways to eliminate unnecessary degrees of freedom through effective Hamiltonians. They have helped to save huge amounts of computational time and effort over the past few decades and are largely responsible for enabling well-converged many-body solid state calculations.

As indicated above, the current generation of ECPs is becoming a significant hurdle in achieving a level of accuracy consistent with high-accuracy all-electron correlated calculations. Presently, valence-only correlated approaches often require additional ECP testing, redeveloping, adapting, modifying, and revalidation in order to get the systematic errors originating in ECPs below a desired threshold. Such practice is highly undesirable and indeed hampers a sizable fraction of QMC and other correlated wave function calculations. This situation is similar to that of past DFT calculations, where pseudopotentials were largely constructed by hand. Optimization approaches²⁵ and systematic testing^26,27 versus all-electron approaches and solid-state calculations have recently been very successful and show that there is a significant room for improvements. However, the substantially higher cost of many-body calculations makes a direct use of these approaches for ECP development quite challenging.

In this work, we envision and assess practical approaches for constructing a new generation of ECPs. We have several goals in mind. First, an effective Hamiltonian that includes an ECP is much more useful if it is reasonably simple and constructed for broad use in basically all theories, i.e., its universality is similar to the original all-electron Hamiltonian. One goal, therefore, is to probe for effective and simple forms that are appropriate for use in both mean-field and/or explicitly correlated methods. In this way, the effective ECP Hamiltonian provides a well-defined reference so that systematic errors from subsequent approximations and approaches are on full display and can be clearly benchmarked.

Another key goal is to increase the accuracy of the new ECPs beyond previous constructions, using measures that we define later. We explore a few strategies for constructing ECPs using correlated wave function methods from the outset. Ideally, these effective Hamiltonians should reproduce, as closely as possible, the behavior of the valence electrons with the original Hamiltonian in a vast range of chemical environments, regardless of whether one calculates a molecule or a condensed system, in an equilibrium or non-equilibrium conformation, and within weak or strong bonding settings. Therefore, our goal is to construct an effective Hamiltonian that mimics the many-body valence spectrum, the spatial structure of eigenstates, and overall scattering properties of the original, relativistic, all-electron atom’s Hamiltonian. Clearly, some compromises will have to be made, and in this work, we investigate the accuracy limits of ECPs of a simple semi-local form with almost a minimal number of parameters and we derive such ECPs for a small set of testing elements. As a guiding principle, we have in mind isospectrality for a subset of states, i.e., we demand that the all-electron and ECP spectra are as close as possible to a set of valence states. The isospectrality is a very general property and applies even in cases when the Hilbert spaces for the two isospectral operators are different, e.g., due to different boundary conditions and different spatial domains.

Note that when we refer to the spectrum, we have in mind the spectrum for not only the neutral atom but also states of cations, anions, and binding curves for molecular systems containing this species. In addition to spectra, we explore the use of many-body spatial information such as correlated single-body density matrices and natural orbitals as well as constructions that are matched to a combination of spatial and spectral information and further iterated constructions. These explorations help us to identify trade-offs between the accuracy for spectra and spatial many-body properties versus demands on the constructions useful for a variety of electronic structure packages. Additional demands such as systematic record of benchmarks from correlated calculations and appropriate updates are further specified in Sec. II.

In Sec. II, we outline a list of properties that we argue are necessary for any ECP that leads to optimal results within high-accuracy valence-only many-body calculations. Section III outlines the inverse-problem approach that we have utilized throughout this work as well as the form of our ECPs. Section IV covers general optimization methods and constructions for generating ECPs consistent with our desired set of properties. Results for the first- and second-row atoms, namely, boron, carbon, nitrogen, oxygen, and sulfur are shared in Sec. V. An analysis of each pseudoatom’s dimer properties as predicted from coupled-cluster with singles, doubles and perturbative triples (CCSD(T)) compared to all-electron calculations is also given in Sec. V. An extensive test of transferability is applied to our most accurate ECP constructions in Sec. VI in order to further ascertain their performance within various molecular environments. We conclude in Sec. VII.

For this new generation of ECPs, we envision including the following set of properties:

• (a)

  Many-body construction is the key and perhaps the most salient point in this list. Some ideas in this direction have been tried before, for example, Acioli and Ceperley²² explored correlated wave functions to generate ECPs and required that the pseudoatom’s density matrix matched that of the full atom beyond a given core radius. Extending on this methodology, Trail and Needs have generated tables of ECPs from MCHF wave functions^18,19 and CCSD(T) excitation energies²⁰ to improve their previously generated table.¹⁶ Additionally, attempts were made to build core effects into ECPs in the studies of Maron and Teichteil²³ and Fromager et al.²⁴ We believe that a more systematic and accurate approach here can still be achieved, especially when considering properties outside of equilibrium (see below). In our construction, we use Coupled Cluster (CC) and Configuration Interaction (CI) methods that are very effective for atomic and small molecular systems, and moreover, we analyze how our generated ECPs perform in non-equilibrium geometries.

• (b)

  The simplicity of the ECP operator has significant benefits since it allows its use in many methods (QMC, quantum chemistry based on expansions in basis sets, and DFT, using established codes). We, therefore, build upon a simple semi-local form that has been in use for some time. Per need, a number of extensions could prove to be useful and could be incorporated in future. This might involve different choices of core-valence partitions, several types of representations such as numerical on radial meshes besides usual Gaussian expansions, fixed radial cutoffs, etc. More desired properties could be added as the project develops, e.g., how to optimally deal with “fat” cores of heavy elements, possible inclusion of core polarization and core relaxation terms, minimizing the fixed-node errors generated by the core part of the wave function without compromising the accuracy, and minimizing the locality approximation in QMC applications following.²⁸ For plane wave applications, the suitability of recasting the ECP into the Kleinman-Bylander form²⁹ could be included as well.

• (c)

  Testing on a set of systems in order to delineate the accuracy limits for energy differences, equilibrium bond lengths, and other properties. Potentially, some of these systems could be included into the retuning set, if important or necessary. For use and further improvements, it will be very useful to have documented systematic errors of generated ECPs documented upfront.

• (d)

  Systematic labeling and updates, i.e., keeping the data, developments, and history on a website³⁰ that can be eventually updated by interested contributors at large.

One way to approach the ECP construction is to formulate it as an inverse problem, i.e., finding an effective operator that produces a subset of valence atomic states such that outside the core they reproduce the all-electron properties as closely as possible. We assume that for this subset both the all-electron and pseudized atoms can be solved in the same systematic and consistent framework, ideally exactly or, in practice, nearly exactly, using the state-of-the-art many-body methods.

What makes the two operators, the all-electron Hamiltonian and ECP Hamiltonian, close? Assuming we are interested in the valence subset of states and properties, there are the following two key aspects:

• (i)

  the two spectra should be the same/very close and

• (ii)

  the spatial characteristics (one- and multi-particle many-body density matrices) of the corresponding two sets of eigenstates outside the core region should be the same/very close.

Note that, in general, the two Hamiltonians differ substantially, in the number of particles/degrees of freedom, presence/absence of ionic Coulomb singularities, and relativity. Since in the reference calculations all correlations (including cores) are assumed to be present, (i) and (ii) implicitly demand that the ECPs capture, as much as possible, also the impact of core-core and core-valence correlations on the valence properties. Our results in Sec. V show that to a certain extent, this is indeed the case.

Finding the desired solution in this setting essentially defines an inverse problem. We can expect that this problem will be ill-conditioned with many nearly optimal or non-unique solutions. Why is this the case? The inclusion of many states into the optimization set often leads to a frustrated optimization problem, e.g., improvement for one state increases penalties for other states. With sufficient number of such frustrating couplings, one ends up with a problem that is qualitatively similar, say, to finding a ground state of a spin glass (a well-known problem in statistical mechanics).

One way how to deal with such ill-conditioned problems is to impose appropriate constraints that limit the space of possible solutions that makes the problem solvable but could result in compromises on the resulting accuracy or incurred biases. Often, it is quite difficult to find the right set of constraints. If the objective function is over-constrained, the search for a minimum is fast, but the solution might be too biased by the constraints. On the other hand, in an under-constrained formulation, the optimization could be very inefficient. Therefore, the goal is to find the best trade-off(s) between these two limits.

The valence-only electronic Born-Oppenheimer Hamiltonian we consider 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 electron i from the core’s origin. The non-local terms contain the projectors on the lmth angular momentum state. The local term, V_loc, is chosen as

where Z_eff is the effective core charge, Z_eff = Z − Z_core. The V_l potential was chosen to consist of a single gaussian term,

All variables labeled by Greek letters are treated as optimization parameters in the minimization of a chosen objective function. Additionally, a constraint that forces a concave shape at the origin is imposed,¹⁵ 

In the case(s) of combined constructions in Sec. IV, the number of Gaussians in each channel is extended.

Note that the chosen form imposes a very significant restriction on the variational freedom of the effective operator. We have used this simple form not only to simplify the optimization problem but also to investigate the true many-body accuracy limit for this “minimal model” version. The results are actually very encouraging. As presented in Secs. V–VI, even with such restricted variational freedom, we were able to construct more accurate effective potentials than the existing ones as well as to illustrate the presented concepts above.

In this section, we give a summary of a few atomic and molecular properties predicted from our constructed ECPs for several atoms from the first and second rows. For comparison, along with our results, we juxtapose the predictions of these properties from correlated all-electron and uncorrelated core/correlated valence for a handful of other ECPs that have been used in many-body calculations in recent years. Furthermore, we present results that utilize different strategies to construct the ECPs, namely, ECPs built from all-electron spectral data only, from spatial data only as well as combined and iterated constructions.

For correlated calculations, a number of codes can be used in a non-relativistic setting. For relativistic calculations, however, the choices are more limited. Nevertheless, it is interesting to see how far one could push the accuracy limits with currently existing methods and codes.

For the all-electron spectral-only references, we used the Molpro quantum chemistry package³¹ to calculate a subset of states from each atom’s spectrum using the CCSD(T) method. To account for scalar relativistic effects, in the all-electron reference, we used the Douglas-Kroll-Hess Hamiltonian throughout. For each atom, the uncontracted aug-cc-pCV5Z basis set (which includes core state correlation functions) was chosen in order to minimize finite basis set errors as much as possible. We used this same basis for both pseudo-atom and all-electron cases so that the basis set errors would largely cancel when gap discrepancies between the two cases were calculated. We carried out limited tests with even larger basis (6Z), and we found out that the gained accuracy had a marginal effect on the spectral differences when compared with other discrepancies such as binding curves for molecular systems presented in Sec. V. In the supplementary material, we show the marginal level of errors from the extrapolation to complete basis set limit for both Burkatzki-Filippi-Dolg (BFD) and all-electron spectral gaps for the first and second row atoms N and S, respectively.

To generate the ECPs, we followed an energy consistent scheme, as described in Ref. 15 and references therein, and minimized the differences in the all-electron and ECP excitation energies in a least-squares way, with our objective function defined as

where s labels a given excited state and ΔE_s is the energy gap between the excited state and the neutral atomic ground state.

We chose to include as many ionizations as possible in the atomic references from the most deeply ionized single-valence state up to (at least) the Mth anionic state where M is the predicted number of bound anions from the all-electron CCSD(T) spectrum. The spectrum here, therefore, includes not only neutral atomic states but also numerous states for the anions and cations. In our investigation, we observed that when the most deeply ionized states were not included in the reference, their corresponding excitations tended to have significant discrepancies from the all-electron excitations and consequently the pseudoatom’s transferability could be negatively affected. Additionally, for each ionization, we expanded the reference to include the bound ground states of all possible total spin channels of that ion, e.g., for an ion with four electrons in the valence space, we included the ground states of the quintet, triplet, and singlet spin channels into the reference, provided that each state was predicted to be bound from the all-electron calculation. The latter choice was motivated by a desire to minimize any possible contamination from energetically lower states within the same symmetry channel.

In order to incorporate spatial information into the pseudopotential, we explored the utilization of the single-body density matrix following the previous work of Acioli and Ceperley. For the all-electron reference and pseudoatom, we generated the single-body density matrices from a singles and doubles configuration interaction (CISD) wave function. We expect that if the pseudoatom density matrix and the all-electron density matrix agree beyond a chosen core region defined by cutoff radius r_c, then the pseudoatom valence will mimic that of the all-electron atom. Given an AE atom and a pseudoatom, we can measure the difference between the density matrices outside of r_c as

If we express each density matrix through the natural orbitals {ϕ_i} and occupation numbers n_i of a correlated calculation such as CISD, we can rewrite the discrepancy between the density matrices as

where the expression $ϕi|ϕjrc$ is the overlap of two natural orbitals evaluated from r_c to infinity.

An important point worth noting here is that one should not expect to rigorously match the all-electron and pseudoatom density matrices. The reason is that the density matrix of the original atom generally contains a contribution from the tails of the core states beyond the cutoff radius. Since it is difficult to disentangle the tails in the many-body setting, this “contamination” will be present. Note that in traditional one-particle formulations this is a non-issue since HF or Kohn-Sham self-consistent orbitals decompose the core and valence states into distinct sets. The closest analogy to the one-particle decomposition comes from the natural orbitals; however, even there the occupation numbers and the natural orbitals will slightly differ (reflecting thus the presence or the absence of the core states). This becomes particularly relevant whenever smaller core radii are imposed. The fits can be further elaborated to approximately take this into account; however, here we opted for a simple agreement assuming that this effect will not dominate the achieved discrepancies.

We further constructed ECPs where we include both the spectral and spatial information into the objective function. Using the integration of the density matrix beyond a cutoff radius as defined earlier, we define our spatial objective function as

where we sum over all states considered. We construct a new objective function where we include both spatial information and spectral information as

where f₀ and g₀ are the optimal values of the respective objective functions. In the case of sulfur in Sec. V E, the two objective functions were left unweighted, given the similarity between f₀ and g₀ in that case.

Throughout, we have utilized the nonlinear DONLP2 optimization code of Spellucci³² for generating the final parameter values of our ECPs with respect to the various objective functions which we have considered. The specifics of the method implemented in DONLP2 are outlined in Refs. 33 and 34, where generally, the solver extends the sequential quadratic programming method, an iterative method that relies on the second derivatives of both the objective function and the constraints such that it can be applied to nonlinear problems. We considered this choice for the chosen solver to be appropriate, given the nonlinear “smoothness” constraint that has been imposed on the parameter sets of our ECPs.

The final procedure that we pursued took molecular data into account. Having generated molecular binding curve predictions from an assortment of ECP constructions, we noticed that there was a potential to improve the ECPs further. In particular, the optimization could be guided to fulfill additional criteria such as reproducing the dimer binding curve in a few iterations. This type of optimization loop produced very small discrepancies from the all-electron CCSD(T) curves along the desired range of geometries even at the steep repulsive side of the binding curves. In particular, for the corresponding dimers, we found that this strategy is able to produce discrepancies as small as 0.05 eV or lower along the binding curve all the way up to the dissociation limit at the repulsive side. Due to time demands of correlated calculations and rather slow coupling between different codes, we guided this level of optimization by interventions to speed up the search. This was also very useful to understand the qualitative relationships between shapes/amplitudes of the ECP versus its properties, and therefore we leave the automation of this part for future. We comment on these constructions at each presented atomic case.

In all of the systems presented below, we compare our ECPs to all-electron CCSD(T), all-electron CCSD(T) with an uncorrelated core (UC), i.e., no excitations from the core, and various ECPs often used in QMC calculations, namely, Burkatzki-Filippi-Dolg (BFD),¹⁵ Trail-Needs DF (TN-DF)¹⁶ ECPs, and where applicable, the shape-consistent correlated Trail-Needs (TN-CEPP) ECPs¹⁸ as well as the recent shape and energy consistent correlated Trail-Needs (TN-eCEPP) ECPs.²⁰ Note that when utilizing the TN-CEPP and TN-eCEPP, we only test the effective core potential and do not include the semi-empirical core polarization potentials utilized by Trail and Needs, given that we were interested in comparing all effective cores within a consistent level of approximation. A key question that we chose to pursue was how well an ECP alone would capture core-valence effects without the need for additional adjustments or approximations. By not including CPPs along with the ECPs allowed us to isolate this effect. Note that Trail and Needs have studied the effect of the core polarization potential (CPP) in Fig. 6 of Ref. 20, illustrating that CPPs do have only minor impact on the dissociation energies over a wide variety of molecules. As explained in results, we found that much larger discrepancies almost invariably appear at shorter molecular bond lengths, and therefore we have focused on addressing this aspect in our constructions.

Here we present our initial results for a number of explicitly correlated ECPs using the various objective functions described earlier. In the first row, we present optimized ECPs for boron, carbon, nitrogen, and oxygen. For the second row elements, we show our results for sulfur.

Boron is perhaps the lightest element where time savings from an appropriate ECP become significant, while the valence space is sufficiently large to represent the atom in chemical settings. Table I shows several constructions of the ECP obtained by minimizing the atomic spectral error only (Spectral) and another that has been adjusted to reproduce the binding curve of the ground state of B₂ in the range of bond lengths (Constructed). Note the very high accuracy that was obtained with minimizing only the spectral discrepancy. A significantly improved dimer solution has been found by minimizing a compromise between reproducing the spectrum and the binding curve. Note that the shortest bond length corresponds approximately to the dissociation point due to nucleus-nucleus repulsion. In solids, such distance between the atoms would correspond to very high pressures roughly beyond 500 GPa. Figure 1 shows the B₂ potential energy surface for the ³Σ_g state, and Fig. 2 provides a set of discrepancies from the all-electron CCSD(T) for our constructed ECPs compared against previous ECPs: BFD, TN-DF, the recently correlated constructions by Trail and Needs (TN-CEPP and TN-eCEPP), and the all-electron uncorrelated core result.

Right away we can observe the following:

• (a)

  The constructed ECPs provide a significantly more accurate picture of the molecular system, both for the valence spectrum optimization (Spectral) and more accurate binding curve optimization (Constructed), than the existing examples. Remarkably, our ECPs show smaller errors for both the spectrum and binding curve than the all-electron correlated valence-valence calculation with an uncorrelated core (UC). In UC, all one-particle states are solved self-consistently including the cores (i.e., it is not the atomic frozen core in that sense); however, only valence-valence correlations are invoked in the CCSD(T) method while core-valence and core-core correlations are neglected. The fact that our effective Hamiltonians provide a better description than UC for the dimer and atomic spectra is quite surprising. We, therefore, observe that the inverse problem formulation enables to mimic, within a certain level of accuracy, effects of core excitations on the valence space.

• (b)

  The simple form of the ECP described earlier is able to accommodate higher accuracy demands.

Note that overall accuracy in reproducing correlated properties is roughly at the level of 0.02 eV, significantly higher than that in previous constructions. Both of these observations are quite unexpected and will be further elaborated in the cases of other elements. Since these solutions were satisfactory, we have not pursued further improvements by including spatial correlations explicitly into the objective function. Note that the spectrum and refinement alone were sufficient to determine the high accuracy solution within the minimal and significantly restricted representation. Our best ECP parameters (Constructed) are given in Table II.

For the sake of clear identification and recording of progress, we introduce here the label for this and subsequent ECPs as “correlation consistent ECP,” version 0.1 or ccECP.0.1, in short. The same label is assigned to all tabulated ECPs in this paper.

Note that the notation can be further expanded, for example, as ttECP.xx.yy.ENSCO. Here the theory would be labeled as tt = cc,hf,df,…, while major.minor release is represented as above by xx.yy, respectively. In addition, fitted quantities could include ENSCO = excitations-norms/shapes-spatialDM-combined/iterated-other approaches. As far as systems used for fit is concerned, one could possibly affix a label ADHOCG = atom-dimer-hydride-oxide-cluster-general/other system. As an example, the full label for our boron case would be ccECP.0.1.EC.

Since the spectral only optimization worked quite well for B, only needing a slight readjustment with our constructed ECP, we test the same for the C atom. The atomic spectrum for various ECPs is given in Table III, and dimer discrepancies in the ¹Σ_g state are given in Fig. 3. The spectral only construction does remarkably well on minimizing the energy differences for the atomic spectra and produces an ECP that is improved over other ECPs. While there is slight overbinding, the spectral curve is relatively flat indicating accurate vibrational frequencies compared to the all-electron. The electron affinity (EA) for the spectral differs with the all electron CCSD(T) value by only −1 meV. The ionization potential (IP) is an order of magnitude better and only differs by +0.1 meV. In addition, a minor compromise in reproducing the atomic spectra produces a much better dimer binding that is within 0.05 eV across the entire curve and has a flat discrepancy, i.e., the same vibrational frequencies as the all-electron curve.

The parameters for our best carbon ECP (Constructed) is given in Table IV.

In constructions of nitrogen, we explored several objective functions and refinement strategies, as shown in Fig. 4 for the molecule in the ¹Σ_g state and Table V for the atomic properties. The spectrum only optimization performs rather well and is a significant improvement over many of the previous constructions. The results show that using the spatial information only such as one-particle density matrices appears less favorable overall. We conjecture that this is caused due to very weak correlation “signal” in the corresponding density matrices. The difference between the mean-field and correlated density matrices is very tiny so that the correlation effect is overwhelmed by the one-particle character that is dominant. On the other hand, the spectrum alone provides a much stronger signal since the eigenvalues basically determine the wave function tails essentially exactly. When combining the spatial and spectral information, the spatial signal was not sufficient to significantly change the parameters, and we obtain essentially the same ECP for both constructions. Although the spectral ECP has excellent atomic properties and a reasonable dimer, we sought to construct a better ECP overall. The constructed ECP was generated by adding additional constraints that were not included in the original objective function in order to alter systematic trends we observed in our various constructions to minimize the dimer discrepancy. In particular, we observed that pairs of ECPs with dimer discrepancies of opposite trends could be produced based on different constraints and subsequently the ECPs could be combined as a linear sum leading to increased agreement with the molecular properties from all-electron CCSD(T) while also preserving a high level of accuracy on the atomic spectrum. This procedure could certainly be used to define a new objective function where the dimer fit is directly included; however, this proved to be computationally inefficient, and we therefore chose to improve the molecular properties by hand as described. The constructed ECP has a discrepancy of below 0.05 eV across the entire curve, although we compromise the spectrum when compared to our spectral construction, particularly for high-spin atomic states. However, if the high-spin atomic states are not considered, the atomic spectrum is better than all previous generations of ECPs. Note that in molecular systems the nitrogen almost invariably appears in low/lowest spin configurations.

Another observation is that ECPs with a combination of minimizing the dimer binding and spectrum discrepancies could be constructed without generating an overly large negative impact to the atomic spectra, provided an appropriate constraint is utilized. Though some accuracy of the atomic properties would diminish as the dimer binding is improved, as can be seen in Table V, we found that a reasonable balance could be obtained between the two when compared to other ECP constructions.

The parameters for our best nitrogen ECP (Constructed) is given in Table VI.

As in the previous case, we test a variety of objective functions when constructing an ECP for O. We present the results in Fig. 5 for the molecule in the ³Σ_g state and in Table VII for the atomic properties. Using a spectral-only construction, we are able to construct an ECP that provides a binding curve discrepancy to the all-electron CCSD(T) curve within 0.05 eV along the entire curve. Since the spectral ECPs are reasonably flat compared to other ECP constructions and the all-electron uncorrelated core (UC), the vibrational frequencies agree quite well to the all-electron CCSD(T) calculations. As in N, the spatial information alone appears insufficient in generating a high-quality ECP. When combining both spatial information and spectral information (Spec/Space), we actually obtain an ECP that is essentially flat across a wide range of bond lengths and shows the best dissociation energy. This also results in a slight compromise on the atomic properties. However, this ECP begins to overbind significantly in the shorter bond length regime which would correspond to high pressure in solids. Based on this and the quality of our spectral only ECP, we believe that the spectral only optimization for oxygen produces the best ECP overall and we do not pursue iterated constructions. The parameters for spectral oxygen are given in Table VIII.

As a last illustration, we include our progress from the second-row atom, sulfur. Construction of the ECP for sulfur has followed similar steps as in previous cases. In Fig. 6, we show the impact on the overall accuracy of the S₂ dimer’s ground state (⁠$Σg(−)3$⁠) potential energy surface from leaving the core uncorrelated though we see that near equilibrium its agreement with the all-electron prediction is at the level of ≈0.03 eV.

For this case, we utilized a Ne core ECP and attempted a set of optimization strategies similar to the investigation of N and O. In this case, we again consider constructing ECPs to mimic all-electron many-body spectral, spatial, and combined properties, in turn. Figures 7 and 8 show CCSD(T) binding energy discrepancies from the all-electron S₂ molecule (⁠$Σg(−)3$⁠) for various approximations to the sulfur cores including our generated ECPs.

For the spectral case, we again have used an all-electron reference formed from bound excitations and included all possible valence ground states of total spin and charge whereby we minimized the spectral discrepancies exclusively with respect to ECP’s parameters. In Table IX, it is shown that the MAD from the all-electron excitation energies is no more than 0.03 eV in this case. This agreement is nearly an order of magnitude improvement over the other approximations to the core. Additionally, with the spectral objective function, we observe that the agreement with the all-electron binding energy at equilibrium is very good with an error of no more than 0.01 eV. For shorter separations of the dimer, however, the spectral ECP undershoots the all-electron binding energy by tenths of eV, and as a result of this large change from equilibrium, we see a non-negligible slope and curvature in the discrepancy at equilibrium which negatively impacts the agreement with the all-electron prediction of the ground state vibrational frequency.

Using the spatial information for the same set of states, we generated the single-body density matrices from all-electron CISD wave functions and subsequently imposed that the single-body density matrix of the pseudoatom’s CISD wave function matches beyond a core radius. The resulting ECPs constructed from this procedure were not adequately transferable which we attribute again to marginal electronic correlations signal in the natural orbitals as we argued for the cases of N and O. It can also be observed from Table IX that the atomic properties are generally negatively affected by optimizing spatial data alone.

For sulfur, we also attempted to match both the all-electron excitation energies and single-body density matrices, simultaneously. For this case, the binding energy discrepancies from the all-electron atom are shown in Fig. 7. Here we see that the error is quite uniform over the entire region plotted, and moreover, it is no more than about 0.06 eV.

Lastly, we again considered additional molecular constraints on the objective function in order to further improve the ECP’s dimer properties as described for the nitrogen atom. The constraints were constructed in such a way as to reach a balance between the spectral and molecular properties. We show its discrepancies in Fig. 7. For this ECP, we see that the dimer’s error is mostly within 0.05 eV throughout the plotted region and the spectral properties are an improvement over both UC and previously generated ECPs. The parameters for sulfur’s constructed ECP are shared in Table X.

One of the desired properties for our ECPs is transferability, i.e., high-quality performance in systems that were not used in our optimization procedure. In the case of our spectral ECPs, the optimization of the parameters only involved atomic properties; for these ECPs, the dimer discrepancies compared to the all-electron results illustrate the transferability of our ECPs to a degree. However, for our constructed ECPs, the parameters were tuned in order to produce improved dimer properties while not sacrificing the atomic properties, and thus the transferability should be verified also on independent examples. Therefore, it is desirable to test those cases in other bonding environments in order to illustrate their transferability beyond the dimers. For this purpose, we have calculated potential energy surface discrepancies for hydrides, oxides, and a handful of additional molecules BH₃, BN, BS, and CN and in order to ascertain their level of transferability. All molecular calculations were performed with either the Molpro quantum chemistry package³¹ or the gaussian09 code.³⁵ 

We show binding energy discrepancies with respect to all-electron CCSD(T) for the NH, OH, NO, SH, and SO molecules in Figs. 9–13, respectively. Furthermore, we summarize the discrepancies in the binding parameters, D_e, r_e, and ω_e of all molecules considered in this work with respect to all-electron CCSD(T) values in Fig. 14, where D_e is the dissociation energy, r_e is the equilibrium bond length, and ω_e is the vibrational frequency. For each case, the parameters and their errors were obtained from fitting the potential energy surface to the Morse potential near equilibrium, where the potential is given as

where a is related to the vibrational frequency by

and μ is the reduced mass of the molecule.

To make the comparison of the pseudoatom transferabilities easier, we also share the MADs of these parameters, both at equilibrium and at the dissociation limit at short bond lengths corresponding to high pressures, in Table XI. From the table, it is shown that our ECPs perform the best overall, where the mean-absolute deviation from the all-electron dissociation energy is smaller than all considered core approximations; furthermore, this improvement and better overall balance are achieved with very limited variational freedom as given by the choice of the ECP form. The only exception is a marginally better MAD for the vibrational frequencies of TN-CEPP with 5(3) versus our 7(2) cm⁻¹.

Our paper presents advancement in the construction of effective core potentials for accurate, correlated valence-only calculations. A key difference from previous constructions being the consistent use of nearly exact many-body approaches to build the ECPs and balancing this with refinement by molecular data that improves the transferability. We introduce isospectrality of all-electron and ECP Hamiltonians on the subspace of valence states as a foundation in formulating the objective functions. Additional criteria that were explored included many-body spatial information from density matrices. This was followed by analyzing the accuracy on dimers, and if significant differences were observed, we included the dimer information to boost the accuracy and transferability. This was done in an iterative manner, and combined constructions were used for efficiency reasons. We were especially careful in obtaining consistent binding curves within a large portion of the bond, i.e., from separations near the dissociation limit to at least equilibrium. We further ascertained the transferability in those cases on hydride, oxide, and selected other molecules that involved the considered element.

The calculations are done at the relativistic level; therefore, energy differences are very close to the actual experimental values for the relevant quantities.

We were able to decrease the atomic spectral errors for B, C, N, O, and S atoms to significantly smaller values when compared to previous constructions.

Interestingly, the spectral discrepancies are also smaller than those in uncorrelated core calculations. Therefore the constructed operators are able to effectively accommodate at least some of the effects from core-core and core-valence correlations as it is clear from comparison with uncorrelated core calculations. Similar argument applies to the improved description of simple molecular systems for a range of bond lengths including very short bonds that are relevant for calculations of materials at high pressures. Rather surprisingly, this implies that for a range of valence properties these ECPs will provide the same or even more accurate results than all-electron treatments with uncorrelated cores. Note the caveat “for a range of valence properties” in the previous statement since this is true only for states/properties that are not too far from the one used in optimization. Clearly, inclusion of larger range of energies and states that can be described would require either additional terms in the effective ECP Hamiltonian or, in general, including more subshells into the valence space (as is common for heavier elements).

Our experience from these constructions shows the following:

• (a)

  the accuracy of even the simplest forms of semi-local ECPs can be boosted very significantly by appropriate many-body constructions;

• (b)

  optimization is quite involved since minor changes in the inputs and the objective function [e.g., increasing the size of the basis set and boosting the spectral accuracy from CCSD to CCSD(T)] complicate the optimizations with many local minima solutions, i.e., we see a clear hallmark of an ill-conditioned problem;

• (c)

  there are differences from element to element influenced by the extent and choice of excitations, occupation, and other details showing thus further complications from the ill-conditioned character of the task.

The obtained results and experience offer a route for constructing comparable quality ECPs for more involved elements such as 3d transition series and beyond. Some adjustments in the constructions and expanding the parameter sets can be expected. Overall prospects for heavier elements would be limited basically only by the availability of codes that can do accurate correlated all-electron relativistic calculations for atomic and small molecular systems.

In this work, we have pursued the construction of ECPs from a many-body framework as opposed to one-particle schemes. In addition, we have included further considerations into the construction such as molecular systems at and away of the equilibrium for better overall transferability.

One-particle constructions have generally involved norm/shape and charge conservations, one-particle eigenvalues, logarithmic derivatives, differences between approximate total energies, etc. Unfortunately, all of these were subject to the biases of the (approximate) method used to solve the underlying atomic problem. The expectation and hope were that in subsequent calculations the systematic cancellations of the underlying biases would be the same as in all-electron setting. In many cases, this worked reasonably well. However, this was not true systematically, and for some elements, states, and bonding environments, the errors are significant. For example, in DFT, all-electron and ECP calculations for just a single atom could lead to large differences. In addition, special adjustments might be necessary such as the nonlinear core corrections for transition metal atoms.³⁶ These complications, however, defy the universal use of such ECPs in other methods, what is one of our stated goals. In addition, there has always been a lingering (and valid) question that why an ECP constructed in DFT or HF methods should be any good in a correlated many-body approach. Clearly, this was one of our motivations to use a many-body framework consistently throughout the construction, in testing and also in further developments.

The validity and limits of the constructed ECPs also become much more transparent if the information about systematic errors built into ECPs is provided upfront so as to clear the stage for subsequent calculations. The simplicity of the basic ECP form is also highly desirable since that enables wide use in many settings and approaches.

Note that presented ECPs are of the simplest semilocal type with minimal form and size, with only one or two Gaussians per channel. It is realistic to expect that more careful and more elaborate fits could further tune the properties beyond accuracies we have currently obtained. The corresponding optimizations are still challenging and are hampered by complicated couplings between the atomic solvers and optimization methods and require further refinement, higher efficiency, and more robustness. At the same time, the presented results offer encouraging examples of the accuracy that can be obtained with present day capabilities and show that there is significant room for further improvements and expansions.

Input and output files for the calculations performed in this study are available³⁷ via the Materials Data Facility.

See supplementary material for basis set extrapolation analysis and a listing of correlation consistent basis sets for each pseudoatom presented in this work.

We are grateful for stimulating discussions with Jerry Whitten during this study. We also would like to thank Paul Kent for a careful reading of the manuscript while in preparation which lead to a number of vital improvements.

The majority of this work (development of the methods, calculations, tests, and writing of the paper) 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 initial theoretical and conceptual considerations were supported by ORNL/UT Batelle, LLC, Subcontract No. 4000144475.

The calculations for this work were performed mostly at Sandia National Laboratories, while some of the calculations have been carried out 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 National Nuclear Security Administration under Contract No. DE-NA0003525.