A new generation of effective core potentials from correlated calculations: 2nd row elements

The effective core potential (ECP) approximation has been vital in electronic structure calculations for several decades. An ECP is used to replace an atom’s core (nucleus and tightly bound electrons) and is constructed to mimic the original core’s effects on the valence electrons. This approximation provides important benefits such as the removal of large energy scales associated with the core electrons while simultaneously enabling the implicit incorporation of the core’s relativistic effects. It generates smoother charge densities around the nucleus which, for instance, allows the use of much smaller plane wave basis cutoff energies or less extensive gaussian basis sets. Consequently, in many cases, the availability of accurate ECPs can be the deciding factor for whether a calculation is even feasible or not. Given that the scope and accuracy of many-body electronic structure methods are continually increasing (true for both stochastic and basis set approaches; see Refs. 1–5), it is important that the errors of the ECP approximation be kept as small as possible or at least kept comparable to other systematic errors that may be present.

Historically, ECPs have generally been formulated within effective mean-field (single-particle) theories where the concept of core-valence partitioning follows rather naturally. A number of approaches have been established in the Density Functional Theory (DFT);^6–12 see also review Ref. 13. For the review of Dirac-Fock energy consistent constructions, see Ref. 14. Although generating valence-only Hamiltonians is more straightforward in an effective one-particle setting, the performance of such ECPs within many-body theories is not guaranteed and can lead to results of mixed quality¹⁵ which can require adjustments to make them reliable within these theories.^16,17 In order to achieve more consistent and systematic accuracy within many-body theories, we have recently proposed a general framework and provided examples of a new generation of correlation consistent ECPs (ccECPs)¹⁸ that aim to reach beyond the current status. The key principles included the use of many-body constructions and benchmarking to ascertain the quality of the ECPs from the outset. The first component of our strategy was based on near-isospectrality, within a subspace of valence states, between the all-electron (AE) and ECP Hamiltonians. Furthermore, the ECP properties were studied not only for atoms but also for small molecules to probe for transferability in bonded situations in order to validate their overall quality.

Remarkably, for the first row, we found that it was possible to derive ECPs that show significant improvements in various valence properties at the CCSD(T) (coupled-cluster with up to double excitations and triples treated perturbationally) level of theory when compared with existing ECP tables or the all-electron uncorrelated-core (UC) approximation where excitations out of the core are suppressed. Quite surprisingly this was achieved with simple, “minimal,” non-divergent ECP parameterizations, i.e., with very few gaussians per nonlocal channel.

In order to maintain the simplicity and direct connection of the ECP to the original Hamiltonian, the solution was formulated as an inverse problem. Namely, we searched for an effective Hamiltonian operator that reproduces a given set of many-body eigenenergies and eigenstates, outside the core region, as closely as possible. Not surprisingly, this can lead to an ill-conditioned problem (high sensitivity to inputs) and/or a complicated optimization landscape with many near-degenerate minima. Consequently, such optimizations can be rather difficult. Such difficulties can be understood conceptually, for instance, if one considers the well-known example from the spectral theory that two operators with different domains can have the same spectrum, i.e., they can be exactly isospectral. This non-uniqueness suggests that the problem requires adequate and proper constraining in order to arrive at a desirable solution with a reasonable amount of effort. For the B, C, N, O, and S elements, we were able to accomplish just that and produced ECPs for these atoms that exceed the accuracy of existing tabulated ECPs and that have led to significant improvements in transferability.

We also aimed to construct ECPs that are simple and can be employed in a variety of electronic structure packages and methods, i.e., as widely applicable as possible. To this end, we chose a well-established form^19,20 that is non-singular and behaves quadratically at the origin and is parameterized by a short sum of radial powers multiplied by gaussians. This form leads to one-particle orbitals that also behave quadratically at the origin which makes their expansion in gaussians more efficient. Additionally, the smooth form provides computational savings for quantum Monte Carlo by reducing fluctuations.

Another aspect of our recently constructed ECPs was to probe for the transferability in bonds both near and out of the equilibrium configurations, having in mind, for example, high-pressure applications in future. Therefore we monitored discrepancies of dimer binding curves including short bond lengths up to the dissociation limit, where the binding energy goes to zero for the compressed bond. Note that this limit is particularly challenging since it directly probes whether the ECP correctly mimics the (missing) core charge at short lengths. In cases where larger errors were present, we included the dimer into the construction as an additional constraint which improved the quality of the resulting ECPs. Indeed, we found that the key molecular properties were more accurately reproduced in this way. As an independent probe of transferability, we then tested hydride and oxide molecules. Further benchmarks and testing will come from future applications for systems where all-electron correlated calculations can be carried out.

In the present work, we extend our previous collection of correlation-consistent ECPs to include the entire 2nd row of elements Na through Ar (following the convention that H and He comprise the zeroth row). For these elements, we have further improved and modified the methodology used in generating the ECPs which we describe in more detail in what follows. The constructed ECPs that we present here are not meant to be definitive. It is clear that with more elaborate forms, say, with a larger number of free parameters, the accuracy might still improve. Another consideration is that different applications might require further refinement or even new constructions. For example, alkaline elements have nominally large cores and one or two valence electrons. Clearly, there are applications where this would not provide an accurate replacement of the all-electron Hamiltonian (high pressures) and the inclusion of the outermost shell from the core to the valence space might be needed. This implies that one needs more options to guarantee the accuracy in various physical or chemical applications. Therefore keeping track of updates and any new constructions, as well as maintaining benchmark sets is an important part of providing clear choices and validations for future use of ECPs.

The paper is organized as follows. In Sec. II, we discuss the different components of a number of atomic energy gaps of a variety of tabulated 2nd row ECPs and how this affects their agreement with all-electron gaps at the many-body level as well as how this motivated the particular construction we use for the ECPs presented in this work. We outline this construction in Sec. III. In Sec. IV, we describe the particular parameterization we use for our ECPs. Section V outlines the scheme we followed to construct the ECPs’ corresponding valence basis sets. In Sec. VI, we share our results and present our ECPs together with several existing ones for comparison and highlight the shortcomings of Ne-core ECPs in particular settings (short polar bonds) and therefore, in addition, present He-core ECPs for the entire 2nd row that significantly improve the accuracy of ECPs with Ne-core partitioning and serve also as benchmarks for comparisons. In Sec. VII, we give our concluding remarks.

We have looked in a more quantitative manner at the contributions from valence-valence correlation energies as well as Hartree-Fock (HF) energies for a variety of atomic excitations among a test bed of tabulated ECPs. This provides additional insights into how various components of the energy are modified and respond to ECP differences (form, parameterization, etc.) which ultimately enabled us to improve our ECP construction scheme and helped to improve their overall performance. Specifically, for a number of tabulated ECPs, we have analyzed the errors (relative to relativistic all-electron results) of the HF and correlation contributions to the CCSD(T) spectral energies for the Si and P atoms and we have looked at how the errors for these energy components change between different ECPs. Relativity in all-electron calculations was included through the 10th-order Douglas-Kroll-Hess Hamiltonian for both HF and CCSD(T) approaches.

In Figs. 1 and 2, for Si and P, we show examples of the spread of HF and CCSD(T) valence-valence correlation errors across a number of previously tabulated ECPs, which employ a wide variety of parameterizations, for a number of excitation energies. In particular, we plot the spread of the following quantities over the set of ECPs

where ΔE_s represents the energy gap between a given state and the ground state, and

where VV_corr is the valence-valence correlation energy. In these figures, the spread (represented by the shaded areas) are bounded by the maximum discrepancy (upper bound) and minimum discrepancy (lower bound) among the various ECPs. For the ECP cases, the valence-valence correlation is taken to be the fully correlated CCSD(T) energy of a particular state. For the all-electron case, the valence-valence correlation is taken to be the UC CCSD(T) correlation energy of a particular state.

There are several observations that can be gleaned from these plots. The first is that the ECPs’ correlation energies are on average larger than the corresponding all-electron valence-only correlated values with the increase varying between 0 and 0.15 eV. This has been known for some time; see the analysis by Dolg.²⁵ The reason is the absence of the radial node(s) and increased smoothness of ECP (pseudo)orbitals since this increases the probability for unlike spin pairs of electrons to encounter each other, thus increasing the correlation energy. However, the next striking feature is the rigidity of the correlated energy differences across various ECPs, regardless of construction, whether it shows singular or bounded/smooth behavior of local and nonlocal terms at the origin or other details. Essentially the ECP correlation energies are almost invariant to the form or the construction of ECP, typically within 0.02 eV or less. Consequently, the largest differences between ECPs for atomic spectra come from fluctuations that originate in the HF component of energy differences. Finally, the third point is that aiming for atomic spectra accuracy better than 0.1 eV requires the HF component to “compensate” for the rigid correlation contribution component for a number of states. Given this rigidity of ECP correlations, we incorporate mean-field data into the construction of the 2nd row ECPs along with many-body data in order to avoid large mean-field biases which we found helps to improve transferability (shown below in Sec. V). In Sec. III we give the specifics of the procedure we use to construct the ECPs.

The strategy that we employ to generate our ECPs is a combination of many-body energy consistency and single-body norm conservation. We provide the details of each and how they are combined to form our full objective function in this section.

In order to account for correlation when generating our ECPs, we have employed a spin-unrestricted CCSD(T) energy consistency scheme as part of our full objective function. The details of this scheme are as follows: For a given atom, we first generate a large number of all-electron reference atomic states which includes various ionizations and spin-multiplicities. The states are calculated with the fully correlated (both core and valence electrons) spin-unrestricted CCSD(T) [UCCSD(T)] method. We account for relativity by using the 10th order Douglas-Kroll-Hess Hamiltonian.²⁶ We formulate this part of the objective function by calculating the same atomic states with a parameterized ECP at the same level of theory, and then vary the parameters until the gaps between these states and the neutral ground state agree with that from the all-electron results to a satisfying level. The discrepancies of the ECP gaps relative to the all-electron gaps are minimized in the least squares way which defines the first component of our objective function, $E2$⁠, given by

where s labels a given atomic state, ΔE_s is the energy gap between the state and the neutral ground state, and “ECP” and “AE” label the ECP and all-electron gaps, respectively.

The energetics from all-electron and ECP were calculated with the Molpro quantum chemistry package.²⁷ For all atoms, the uncontracted aug-cc-pCV5Z basis set,²⁸ where core state correlation functions are present, was utilized. If linear dependencies developed as a result of uncontracting the basis, the problematic terms with nearly identical exponents were removed. This same basis was used for calculating both the all-electron and ECP gaps in order to construct similar basis set errors between the two which would largely cancel when calculating their difference in Eq. (3). Quantitative estimates of the residual finite basis set errors were given in our previous work¹⁸ for a test set of atoms from the 1st and 2nd row.

We chose not to use the $E2$ objective function exclusively when training our ECPs given that in doing so we observed an overall large increase in error of the HF gaps and subsequently poorer transferability when calculating potential energy surfaces of diatomic molecules (dimers and oxides) at the CCSD(T) level of theory. We suspect this is due to the rigidity in the ECP correlation energies as outlined in Sec. II, namely, in order to reduce the total error of the CCSD(T) gap, the magnitude of the HF error increases to compensate the error from the more rigid correlations. We observed that reaching small gap errors at the many-body level at the full expense of the HF gap errors showed tendencies to increase transferability errors. In fact, the best results came from a balance between the errors in correlation and HF components. In addition, significant variational freedom inside the core is proved to be helpful and we have therefore included a single-body norm conservation into our objective function, as outlined in Subsection III B, in order to strike this balance.

In order to maintain accuracy at the HF level of theory, we employed a norm conservation scheme as another component of our full objective function. We have taken a strategy similar to what has been used in previous studies,^20,22,24,29 however, the specifics of our approach are as follows. For the Ne-core ECPs, we first generated all-electron scalar relativistic HF radial orbitals, $ϕlAE$⁠, where the 10th order DKH (Douglas-Kroll-Hess) hamiltonian was employed. For each angular momentum channel, we chose a cutoff radius, R_l, corresponding to the outer most extremum of $r45ϕlAE$⁠. At this radius, we then obtained the orbital’s norm inside this point, $NlAE=∫0Rlrl+1ϕlAE(r)2dr$⁠, its value, $VlAE=ϕlAE(Rl)$⁠, and its slope, $SlAE=ddrϕlAE(r)Rl$⁠. The norm-conserving objective function is then formed as

where we also include the errors in the orbital eigenvalues ϵ_l. This procedure is similar to what is generally done in norm-conserving schemes; however, we allow the core more freedom by not matching it to a particular functional form but rather focus on matching the quantities mentioned. For the elements Al–Ar, we took reference s and p orbitals from the atom’s neutral ground state and took the reference d orbital from the p → d excitation. For Na and Mg, we took the single-valence doublet states as the reference in each channel.

To form our full objective function, O², we take $E2$ and $N2$ as a linear combination

where $w$₀ is a chosen weight, typically taken as 5 × 10⁻². To construct our ECP, we vary its parameters until O² is minimized. The value of $w$₀ is chosen so as to strike a balance between the accuracy of the many-body and single-body properties of the ECP—or viewed in another way, we are allowing a violation of single-body norm conservation in order to achieve higher accuracies on the many-body spectral properties.

When using the conventional Ne-core and 3s3p valence, the transferability of the ECPs in particular environments was observed to be less than satisfying. This is particularly true for compressed polar bonds (as will be shown). In order to reach high many-body and single-body accuracies in these regimes as well, we have also generated He-core ECPs for all of the second row atoms. Though more expensive than their Ne-core counterparts, the He-core ECPs still provide a non-negligible computational savings when considering their use in quantum Monte Carlo (QMC) which scales with the atomic number as $O(Z5.5−6.5)$ and within planewave codes where calculations with the full Coulomb potential are not feasible. For these small-core ECPs, the optimization procedure was nearly identical, except that we took the neutral ground states for all references in the norm-conserving component and we generated the cutoff radii from the semi-core orbitals. In this case, for each angular momentum channel, the cutoff radius was taken as the inner most finite extremum of $rϕlAE$⁠. We then used these radii as the match points of the norm, value, and slope of the valence orbitals while also matching their eigenvalues. Furthermore, we found that a much larger weight could be used on the UCCSD(T) spectral component (two or three orders of magnitude) than the weight used for the Ne-core ECPs without introducing a significant negative impact in transferability. This is due to a much smaller discrepancy in the correlation energies (valence and semi-core electrons) between the He-core and all-electron cases.

We have applied the construction outlined in this section to all the second row atoms, Na–Ar. For all optimizations, we have used the DONLP2 solver of Spellucci³⁰ to generate the optimal parameter sets for all ECPs. The methodologies used in DONLP2 are outlined in Refs. 31 and 32. Though finding the global minimum in these types of non-linear optimizations can be difficult, our incorporation of the many-body spectral component and mean-field spatial component (i.e., additional constraints) should alleviate issues related to the ill-conditioned nature of this problem. Further justification of our chosen objective function came after analyzing the accuracy of our ECPs in settings that were not included as part of their training, which we cover in Sec. VI.

The parameter values of our ccECPs for all atoms and all core sizes are shared in Tables I and II and are also available at the website in Ref. 33.

The ECPs we have constructed have the following form:

where r_i is the radial distance of electron i from the core’s origin and l_max is chosen to be the maximum occupied angular momentum channel from the ground state core. The non-local terms include lmth angular momentum projectors, |lm⟩⟨lm|.

The local potential, V_loc, is taken to be

where Z_eff = Z − Z_core is the pseudo-atom’s total core charge. We fix β_loc,1 = Z_eff, β_loc,2 = Z_effα_loc,1, while n_loc,1 = 1, n_loc,2 = 3 so that V_loc(r) is finite everywhere and its first derivative vanishes at the origin. Similarly, the set of non-local potentials, {V_l}, are parameterized by

where the powers {n_l,k} are all equal to 2. Except for β_loc,1 and β_loc,2, which are fixed as given above, the Greek letter parameters are all varied in the minimization of our objective function. For an initial round of optimization, for V_loc and all V_l’s, we set the number of gaussian terms to one. For a given Ne-core ECP, we included 1-2 more gaussians in these channels for additional flexibility if further optimization rounds were required to reach agreement within roughly 0.05-0.1 eV of the low-lying all-electron many-body spectrum. Pushing the accuracies further, generally led to significant mean-field errors for reasons discussed in Sec. II. In the case of the He-core ECPs, we again only included one gaussian term in these functions initially. However, for these ECPs, further terms were added until we observed agreement with all-electron data within roughly the chemical accuracy window. This simple form keeps the volume of the parameter space as small as practically possible and therefore helps to simplify the optimization problem.

For each atom and each core size, we have generated a corresponding basis set. The general recipe that we followed to generate the basis sets (similar to Ref. 19) was the same for both core sizes, and we outline it in this section.

To generate the basis sets for the He-core ECPs, we have minimized the atomic ground-state CCSD(T) energy for all atoms. The basis sets for the Ne-core ECPs minimize the atomic ground-state CCSD(T) energy for the atoms Al through Ar and the dimer ground-state CCSD(T) energy for the atoms Na and Mg. For a given basis set, we start by generating contractions, one for each occupied HF ground-state orbital, by starting from a set of even-tempered primitives and optimizing the largest exponent until the ground-state CCSD(T) energy is minimized. The resulting exponents are then taken as the contraction primitives whose coefficients are taken to be the ground-state HF expansion coefficients. We take the same contractions for each VDZ, VTZ, VQZ and V5Z basis. For a given basis, we append additional even-tempered uncontracted primitives following the correlation-consistent scheme³⁴ and determine the largest exponent among these primitives by minimizing the ground-state CCSD(T) energy.

In addition to these basis sets, we have also generated augmented varieties which follow the same scheme, but include an additional even-tempered diffuse primitive into each symmetry channel. All of our basis sets are shared at the website given in Ref. 33.

For all 2nd row atoms, we share in Tables III–X the all-electron CCSD(T) valence ionization potentials and (if bound) the electron affinity along with the discrepancies of the same quantities from various core approximations which include all-electron UC approximation, a number of previously tabulated ECPs,^19,22–24 and our Ne-core and He-core correlation consistent ECPs, ccECP[Ne] and ccECP[He], respectively. We also share the discrepancies from all-electron UCCSD(T) potential energy surfaces for one or two diatomic molecules for all the atoms in Figs. 3–10. For the cases of oxides, we utilized our previously published oxygen ECP¹⁸ when calculating potential energy surfaces for our ECPs, and for all other ECPs we used that table’s corresponding oxygen ECP. For the molecular discrepancies, we span a large range of geometries in all cases, from near the dissociation threshold (⪅0.05 Å) on the left to just past equilibrium. For argonium, we used a well-tuned unpublished hydrogen ECP for all cases other than the uncorrelated core. Our ECP parameter sets are shared in Table II.

In Subsections VI A–VI H, we provide a discussion for each pseudoatom’s atomic and molecular accuracies, in turn, and summarize the results. In Subsection VI I, we share a comparison of the mean absolute deviations of the binding parameters of the molecules presented in this work between our ECPs and the all-electron UC approximation as well as the other tabulated ECPs we have benchmarked against in this work to provide a more global picture of the accuracy of our potentials. In this section, we also present total atomic energy components of our ECPs along with their core radii for future reference.

In the case of sodium, the atomic and molecular results of our Ne-core and He-core ECPs are shared in Table III and Fig. 3, respectively. For the case of the Ne-core, our construction led to a significant improvement of the atom’s first ionization potential, when compared to the all-electron UCCSD(T) result, relative to the previously tabulated ECPs as well as the uncorrelated core (UC) result—the error being roughly a factor of 3 smaller than the other ECPs and roughly a factor of 2 smaller than UC. The electron affinity of our Ne-core ECP is also in good agreement with the all-electron result with an error that is less than 0.01 eV. In the case of our He-core ECP, the agreement with the all-electron atom is better still and for each quantity, the errors are at sub meV scales.

For the Na₂ and NaO molecules, the accuracies of both of our core partitions are also quite high with our Ne-core and He-core ECPs being within 0.05 eV and 0.01 eV of the all-electron UCCSD(T) result over the full range of geometries, respectively.

For magnesium, the atomic and molecular results of our Ne-core and He-core ECPs are shared in Table IV and Fig. 4, respectively. The Ne-core partition shows an increased accuracy over UC and the previously tabulated ECPs with respect to the first and second ionization potentials. The He-core partition further increases the accuracy by more than an order of magnitude, with errors less than a meV for these quantities.

In the case of the MgO molecule, our Ne-core ECP shows better agreement with all-electron UCCSD(T) at shorter atomic separations than both UC and previously tabulated ECPs and at equilibrium, the agreement is within chemical accuracy. Our He-core ECP improves upon this further, with errors less than 0.01 eV across all geometries plotted.

For aluminum, the atomic and molecular results of our Ne-core and He-core ECPs are shared in Table V and Fig. 5, respectively. In the Ne-core case, the errors over the ionization potentials and the electron affinity were not improved over the other core approximations but instead remained comparable. We note the lower variational freedom for our case due to bounded forms and smaller radial range when compared with some of the constructions. However, in the case of the He-core, we observed errors less than a hundredth of an eV.

The Al₂ molecule is described well by both core partitions, both being within chemical accuracy scales. For the AlO molecule, similar to previously tabulated ECPs, our Ne-core ECP is overbound by tenths of eV as the molecule is compressed toward dissociation. This potentially signifies a limitation of this particular choice of core for this type of environment, namely, compressed polar bonds. We suspect that key contributors to the error in this regime are a breakdown of the assumed point-charge interactions of the cores and also the neglect core-valence overlaps.^14,37 For our smaller He-core ECP, we see that these difficulties evaporate and we observe errors well within chemical accuracy.

For silicon, the atomic and molecular results of our Ne-core and He-core ECPs are shared in Table VI and Fig. 6, respectively. In the Ne-core case, like Al, the accuracies among the ionization potentials and electron affinity were not improved over the other core approximations, though the He-core ECP shows errors smaller than chemical accuracy for all of these quantities. The results from the Ne-core ECP potentially reveal the limitations of this particular core partition.

In the case of the silicon dimer, both core partitions perform very well and their errors at the UCCSD(T) level are no more than about 0.02 eV. The SiO molecule shows similar and profound discrepancies for the Ne-core ECP, as was seen in AlO. At short atomic separations, the pseudo-molecule is severely overbound for all tested ECPs including our Ne-core partition. Interestingly, at equilibrium, the agreement for our ECP and some others is much better and is within chemical accuracy. For our He-core ECP, the errors are significantly reduced and we observe chemical accuracy agreement across all geometries plotted. On this element, it is clearly visible that an atomic spectrum is not sufficient to assess the overall quality of the ECPs. Seemingly, BFD, TN, and STU have somewhat smaller energy discrepancies than our Ne-core ECP. However, molecular calculations reveal lower or even significantly lower transferability. In fact, our Ne-core ECP is on par with SBKJC that however is unbounded (1/r² terms) so that our variational freedom is then significantly lower in comparison.

For phosphorous, the atomic and molecular results of our Ne-core and He-core ECPs are shared in Table VII and Fig. 7, respectively. For the case of our Ne-core partition, we observed that comparable accuracies of the atomic properties were only achievable with respect to the other core approximations. However, significant improvement was again achieved for the case of our He-core partition.

For the case of P₂, both core partitions perform well with the Ne-core ECP obtaining maximum errors that are only marginally larger than chemical accuracy and the He-core ECP’s errors being a factor of 2-3 smaller still. The results from PO are also quite good for both partitions. At equilibrium, both He- and Ne-core ECPs show errors no more than about 0.03 eV. For our He-core ECP, this level of accuracy is maintained as the bond is compressed; however, the Ne-core ECP at these shorter separations begins to overbind. Note that our constructions compete in overall optimality with SBKJC while having smaller variational freedom in bounded potentials. In addition, the molecular results show rather nonsystematic behavior for STU ECPs that are different (underbinding) from the rest of the row (typically, overbinding).

For sulfur, the atomic and molecular results of our Ne-core and He-core ECPs are shared in Table VIII and Fig. 8, respectively. In this case, the accuracies of the Ne-core ECP at the UCCSD(T) level are rather good with errors within chemical accuracies throughout all bond lengths plotted for both S₂ and SO molecules. And for our He-core ECP, we again see the benefit of this particular choice of core-valence partitioning where the atomic accuracies show large improvements over the other Ne-core approximations and ≈0.02 − 0.03 eV accuracy is achieved for both the sulfur dimer and oxide across all geometries plotted.

For chlorine, the atomic and molecular results of our Ne-core and He-core ECPs are shared in Table IX and Fig. 9, respectively. For this case, we only observed a significant improvement of the ionization potentials and electron affinity for our He-core ECP—its errors reaching no more than about 0.1 eV even for very deeply ionized cases.

Our construction, however, led to He- and Ne-core ECPs that reproduce the all-electron Cl₂ and ClO potential energy surfaces quite well—where both large and small core ECPs remain well within chemical accuracy for the full range of geometries. To provide another option with lower spectral errors, given the larger errors of our Ne-core ECPs in deeply ionized cases, we also generated a Ne-core ECP that weighted the energy component of the objective function more heavily to reduce the spectral errors further at the cost of slightly larger errors in the accuracy of Cl₂ and ClO. We share this alternative ECP in the supplementary material.

For argon, the atomic and molecular results of our Ne-core and He-core ECPs are shared in Table X and Fig. 10, respectively. With respect to the atomic properties, the Ne-core ECP performs well in this case with a total mean absolute deviation that is comparable to the all-electron uncorrelated core results. The He-core again performs even better where the total mean absolute deviation is improved over our Ne-core ECP by another 0.4 eV. To test the ECPs’ performances in a molecular setting, we calculated the errors from the all-electron UCCSD(T) ArH + binding curve. Both core partitions show high accuracies for this case where the Ne-core and He-core errors are no larger than ≈0.01 eV in magnitude. Similar to chlorine, we also generated a Ne-core ECP that weighted the energy component of the objective function more heavily to reduce the spectral errors further at the cost of larger errors in the accuracy in the molecule. We share this alternative ECP in the supplementary material.

We compiled all the results for molecular equilibrium parameters into the Table XI. The statistics provide useful information for an overall performance of the presented constructions as well as comparison with previously published tables. We have achieved rather significant improvements in the molecular properties due to the more flat discrepancy curves than other ECPs. In general, the ECPs have a tendency to overbind for short bonds and the cause is the missing core charge and the related core-valence repulsion. Interestingly enough, except for Al and Si, we were able to keep these biases typically below 0.05 eV with mildly increased deviations (0.1 eV) at very small bond length for Mg. Clearly, minor improvements might still be possible although we expect that one would encounter diminishing returns unless the ECP form is further generalized. On the other hand, He-core framework clearly eliminates these deficiencies and provides very high accuracy across the whole row. Finally, we list total HF and correlation energies from UCCSD(T) as extrapolated to the complete basis set limit in Table XII along with the core radii of our ccECPs. While HF energies are saturated to about 0.1 mHa or better in all cases, for Ne-core, we estimate extrapolation errors on correlation energies to be within about 1 mHa. For He-core correlation energies, we guess underestimation by about 3 mHa for Na and up to about 6 mHa for Ar, showing limitations of 5Z basis set and excitations up to perturbational triples to fully capture the many-body effects of 16 electrons. This estimate comes from an approximate analysis of core-core, core-valence, and valence-valence correlation energy components in all-electron calculations³⁸ with an added correction from Dolg²⁵ that reflects increase in correlations in the ECP setting. The core radii, r_l, of each ccECP is taken as the radial distance from the origin at which the channel’s full potential agrees to within 10⁻⁵ Ha of the bare Coulomb potential, −Z_eff/r. Additionally, we share the cutoffs, r_l,nl, of each channel’s non-local potential which we define as the point where the non-local potential drops below 10⁻⁵ Ha. The non-local radii have implications primarily related to efficiency of QMC calculations as further explained elsewhere.³⁹ Additional data can be found also in the supplementary material.

In this work, we offer several new advances and insights into the construction of correlation consistent ECPs and apply this to the construction of potentials for the 2nd row elements.

In certain aspects, the 2nd row elements appear to be somewhat more complicated than the 1st row even though the corresponding valence spaces are very similar and the total and kinetic energies are lower (by about 20% when compared with the 1st row) indicating thus smoother electronic densities. Intuitively, therefore, accurately solving systems of 2nd row elements might be considered “easier” than, say, describing the deeper lying and more localized p-subshells in systems containing 1st row elements. However, the small number of valence electrons combined with larger core sizes generates larger biases than in the 1st row. In particular, the core-valence effects in correlation (important for spectra) compete with the absence of electrostatic core-valence interactions (important for short bonds) and therefore these tendencies are more difficult to reconcile. We found that the fits with many-body spectra only had limited accuracy, but when we included more of the mean-field spatial information, we were able to reach the most acceptable compromises.

We were lead to the particular construction scheme introduced in this work after analyzing the behavior of the one-particle HF energies together with the correlation energies. Not too surprisingly, we found that the correlation energies for particular states and for a broad variety of existing ECPs with different parameterizations were essentially constant—the changes were on the order of 0.01-0.02 eV. This appears consistent with the key tenet of effective core potential theory that implies that the correlation will faithfully “follow” the constructed effective valence operators. However, we found this to be correct essentially only up to about 0.1 eV accuracy for atomic spectra. As soon as higher accuracy was sought and molecules were considered, the interplay between various energy components became more complicated and required refinements.

The testing of transferability proved to be very revealing. In particular, for the aluminum oxide and silicon oxide molecules at short bond lengths, we see a significant departure from the assumed accuracy standards of ≈0.05 eV for both ours and existing ECPs. We traced the deficiency to the missing electron core tails and likely core-valence effects which are absent in the ECP and would need to be included to reach higher accuracies. Similar problems have been noticed before; see for example, Ref. 40. Here, we present a systematic view on this issue for Ne-core ECPs within polar bonds at short separations and, in addition, we point out that elimination of this issue might require more elaborated solutions. We have tested the inclusion of classical effective repulsive charges that would correspond to the electronic core-core interactions that could reasonably restore the desired accuracy. In order to introduce this in a more thought-through fashion, we leave this aspect for future studies. This is a notable point that can affect many oxide calculations⁴⁰ such as Al and Si perovskites at high pressures (for example, these are important components of the Earth’s mantle). Due to the importance of oxides, it is highly desirable to probe for these particular aspects in calculations of solids and other systems.

The atomic spectra and Table XI show that apart from the inherent bias for a few short polar bonds, our constructions show consistent accuracy in both atomic and molecular calculations and offer close to an optimal balance overall. Our results, somewhat unexpectedly, are on par with (or in a few cases significantly improve upon) the ECPs from the SBKJC table. Except for the Na atom (which is a bit special having only one valence electron), SBKJC ECPs appear to represent perhaps the most systematic consistency among the tested previous constructions. Note that the SBKJC ECPs are not bound at the nucleus, diverging as 1/r² in the repulsive channels and as −Z_eff/r in the local channel. Due to our choice of the ccECP form with bounded potentials and minimized radial extent of the nonlocal terms (although not explicitly enforced) for several elements, we hit a comparably optimal point as has been achieved in the SBKJC table before. In some cases, we were able to actually find significantly better compromises. As mentioned before, we note that Ne-core ECPs for atoms such as Na and Mg with one or two valence electrons might prove to have only limited applicability in some bonding environments where shallow core states contribute to the chemical effects (for example, in compounds at high pressures).

It is also clear that core polarization and relaxation potentials^41–44 can alleviate some of the mentioned deficiencies. In this study, however, we wanted to see the limits of the simplest ECP form before employing new terms. In addition, oxide molecules for elements in the middle might be difficult to “fix” since the root of the problem appears to be the missing repulsive (essentially classical) charge from core electrons and this might require a more general and more extensive modification of the ECP form.

On the other hand, since we report here on the existing limits of the Ne-core ECPs in their current forms, we have decided to provide another option for high accuracy calculations that offer a significant boost in accuracy and eliminate many compromises mentioned above. To this end, we constructed high accuracy He-core ECPs that significantly decrease the discrepancies with all-electron relativistic calculations. The gain is by an order of magnitude or more, with typical errors of 0.001-0.005 eV on the atomic side and ≈0.01 eV for molecular systems. Obviously, these operators imply additional cost that is most pronounced in stochastic methods due to the increase of total energies and fluctuations. Note that some of that is compensated from a much smaller effective core radius of nonlocal terms.³⁹ 

Overall, we believe that the constructed ECPs provide both new insights as well as practical new options for high accuracy valence-only calculations. The supplementary material contains further information including very accurate CCSD(T) total energies extrapolated to the complete basis set limit (better than 1 mHa accuracy) for atomic states and relativistic CCSD(T) molecular potential energy surfaces and binding parameters from all-electron calculations as well as the discrepancies from all core approximations considered in this work—experimental binding parameters are also included for comparison. All ccECPs and basis sets are accessible at the website.³³ Besides this data, the presented analysis offers a clear comparison between existing ECPs and should open new perspectives for further advances in this important research area.

Input and output files for the calculations performed in this study are available and maintained at the Materials Data Facility.⁴⁵ 

See supplementary material for additional calculated properties of our ccECPs including total energy components of various atomic states and the potential energy surfaces of the molecules considered in this work along with experimental binding properties.

We would like to thank P. R. C. Kent for reading the manuscript and for helpful suggestions. We are also grateful for the thorough reading of the manuscript and insightful suggestions from one of the referees which lead to significant revisions and improvements to the manuscript. 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 calculation has been carried out at TACC under XSEDE.

Sandia National Laboratories is a multimission laboratory managed and operated by the 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.