Recently, we developed a new method for generating effective core potentials (ECPs) using valence energy isospectrality with explicitly correlated all-electron (AE) excitations and norm-conservation criteria. We apply this methodology to the 3rd-row main group elements, creating new correlation consistent ECPs (ccECPs) and also deriving additional ECPs to complete the ccECP table for H–Kr. For K and Ca, we develop Ne-core ECPs, and for the 4p main group elements, we construct [Ar]3d10-core potentials. Scalar relativistic effects are included in their construction. Our ccECPs reproduce AE spectra with significantly better accuracy than many existing pseudopotentials and show better overall consistency across multiple properties. The transferability of ccECPs is tested on monohydride and monoxide molecules over a range of molecular geometries. For the constructed ccECPs, we also provide optimized DZ-6Z valence Gaussian basis sets.

Effective core potentials (ECPs) are widely used in electronic structure calculations of molecular and condensed systems. ECPs simplify these calculations by eliminating the core electrons, which usually do not significantly contribute to the valence electronic structure. Additionally, scalar and even spin-orbit relativistic effects can be included in the ECP in a straightforward way. For some calculations, the use of ECPs is essentially unavoidable, especially for cases involving very heavy elements. For instance, quantum Monte Carlo (QMC) methods scale poorly with nuclear charge, Z, as O(Z5.56.5);1,2 without the ability to reduce Z via an ECP, QMC simulations would quickly become computationally prohibitive for heavy elements.

In recent decades, a number of ECP tables have been developed from various construction schemes and concepts. One such strategy is to use shape-consistency/norm-conservation of the single-particle orbitals, where for a given reference state, the pseudo-orbitals are made to match the all-electron orbitals outside of a chosen cutoff radius. An alternative approach has been to fit AE Hartree-Fock/Dirac-Fock (HF/DF) atomic excitations and gaps, known as energy consistency. These published tables have led to substantial progress by enabling large-scale calculations of both molecular and condensed systems. However, many of these studies have involved methods with mixed levels of accuracy such as Density Functional Theory (DFT) or correlated wave function approaches far from convergence.

Advancements in many-body theories and progress in computational platforms have made the quality of ECPs crucial for obtaining high accuracy for larger systems and a broader range of elements. To the best of our knowledge, the use of many-body theories and incorporation of correlation effects in the development of ECPs have been rather limited.3–8 

Very recently, we have constructed correlation consistent ECPs (ccECPs) for the first two rows (periods 2 and 3)9,10 as well as the 3d transition metal series.11 In this work, we expand the ccECP library to include the 3rd-row main group elements, i.e., K, Ca, Ga, Ge, As, Se, Br, and Kr. We also derive additional ECPs, hitherto missing from the ccECP table, i.e., H, He, Li, Be, F, and Ne.

Let us briefly recall the ideas and the guidelines in the generation of ccECPs.9–11 

  1. We use a many-body construction, i.e., the incorporation of correlation effects. Fully correlated, relativistic AE calculations are taken as a reference as opposed to mean-field frameworks such as HF or DFT. Coupled cluster with singles, doubles, and perturbative triples [CCSD(T)] is used in the optimization.

  2. We use a simple and minimal parameterization for the ECP. The ccECPs are expressed in a well-established functional form12 that is captured by a few Gaussians per symmetry channel.

  3. We test the transferability of the ccECPs in molecular systems, ranging from compressed to stretched molecular geometries. This serves as a check on the ccECP outside equilibrium.

  4. We establish an open database for future developments. The ccECPs are hosted at http://pseudopotentiallibrary.org with systematic labels and updates, which is open to other contributors.

Obviously, our ccECPs do not represent the ultimate accuracy limit of pseudopotentials, although we believe that we have obtained results that are close to these accuracy limits within the given form, number of parameters, and core choices considered. Clearly, more elaborate parameterizations are possible for finer accuracy targets. For instance, core-polarization effects and explicit spin-orbit coupling operators could be added subsequently, but that is beyond the scope of this paper. We have tried to achieve chemical accuracy (1 kcal/mol) for all properties, and in many cases, we have been successful. In cases where this demand was too strict, we have constructed the ccECPs to show systematic and balanced errors that are more consistent than most of the published tables thus far.

The paper is structured as follows: We discuss the chosen functional form and parameterization for the ccECPs in Sec. II A. In Sec. II B, the objective function and optimization strategies are examined. The results for K and Ca are presented in Sec. III A and Ga–Kr are presented in Sec. III B. We provide additional optimized ccECPs for the selected elements in Sec. III C that includes H, He, Li, Be, F, and Ne. Some of these additional ECPs also include an AE version that smooths out the Coulomb divergence in the potential. Together with our previous papers, this covers all the elements from H through Kr. The optimization of valence basis sets is discussed in Sec. III D. Finally, we present conclusions and further discussions in Sec. IV.

The goal of this work is to approximate the relativistic AE Hamiltonian with a simpler valence only Hamiltonian, Hval, where we assume the Born-Oppenheimer approximation.13 The valence Hamiltonian is given by the following expression in atomic units (a.u.):

Hval=i[Tikin+ViECP]+i<j1/rij.
(1)

We employ a semilocal ECP functional form given as

ViECP=Vloc(ri)+=0maxV(ri)m|mm|,
(2)

where i is the electron index and ri is the electron-ion distance. The nonlocal term contains angular momentum projectors on the ℓmth state. max is the maximum angular momentum number of nonlocal channel projectors included in the ECP. Here, we choose max = 2 for 4p elements due to 3d10 electrons being in the core. For K and Ca, max = 1. Vloc(ri) is the local channel; where at the origin, Coulomb singularity is explicitly canceled out; and Vloc(ri)’s first derivative vanishes at the origin,12 

Vloc(r)=Zeffr(1eαr2)+αZeffreβr2+k=12γkeδkr2,
(3)

where Zeff is the effective core charge, Zeff = ZZcore, and α, β, γk, and δk are coefficients to be optimized. The choice of a nonsingular form results in smooth orbitals that reduce computational resources such as sizes of basis sets as well as higher efficiency sampling in QMC. The nonlocal term V(ri) has the following form:

V(r)=j=12βjeαjr2,
(4)

where βlj are αlj are optimized nonlocal channel coefficients. This simple and compact form of ECP is chosen so that the ECPs can be readily used in standard quantum chemistry packages and other electronic structure codes in general. The simplicity of the form also makes the optimization process easier as commented upon later.

Similarly to our previous ccECP constructions, the objective function is a weighted combination of several components such as energy consistency/isospectrality and terms related to orbital shape and norm conservation. The energy consistency part of the objective function is given as follows:

E2=sSΔEsECPΔEsAE2,
(5)

where the energy “gap” ΔEsECP is the CCSD(T) energy difference for the state s relative to the ground state of the pseudoatom. The states that we include are several neutral excitations, a number of ionized states, and electron affinity states for anions that are bounded. ΔEsAE is the corresponding quantity in the all-electron atom with 10th order scalar relativistic Douglass-Kroll-Hess (DKH) Hamiltonian.14 Therefore, relativistic effects are inherently included in our ccECPs.

It is important to note that all AE calculations in this work are fully correlated. In other words, they include valence-valence (VV), core-valence (CV), and core-core (CC) correlations. Typically, ECPs include only VV correlations; however, the inclusion of CC + CV + VV correlations in the AE reference is essential for several aspects as well as in the optimization, as will be shown in Sec. III.

Isospectrality is a crucial property of an ECP since it relies on actual observable quantities. However, as we have shown previously, the optimization of the spectrum alone can result in reduced transferability of the ECP since the atomic spectrum can be overfit at the expense of other properties. Therefore, we make use of the additional norm/shape conservation related function which is given as

N2=ϵECPϵAE2+γNECPNAE2+VECPVAE2+SECPSAE2,
(6)

where the following ECP quantities are matched to relativistic AE counterparts as explained below. ϵAE are one-particle eigenvalues, and these are matched for all angular momentum channels. To evaluate the quantities inside the brackets, a cutoff radius R is chosen as the outermost extremum of r45ϕAE. Then, the norm NAE=0Rr+1ϕAE(r)2dr within that radius, the value VAE=ϕAER and derivative/slope SAE=ddrϕAE(r)R at that radius are matched. Usually, setting γ = 0 resulted in accurate molecular properties, however, in some cases, we set γ = 1 to further improve the transferability. The total objective function is then

O2=w0E2+N2,
(7)

where w0 was chosen to be 0.05, although adjustments were made in some cases.

Calculation of E2 is the most expensive/demanding part of the objective function evaluation since we included ∼10 atomic state energies where each must be evaluated at the CCSD(T) level with a large, uncontracted basis set. In order to alleviate this cost, we proposed a recipe to significantly speed up the optimization procedure in our 3d transition metal paper.11 Here, we briefly recall that approach. It is known that the correlation energies of various ECPs, even constructed in different approaches, are essentially the same, unlike their corresponding HF energies. This allows one to map the CCSD(T) gaps to a mean-field gap by shifting the corresponding HF energy with a quantity related to differences in AE and ECP correlation energies. Therefore, the ECP correlation energies need not be calculated at every step of the optimization and the objective function is evaluated with a cheaper, appropriately shifted HF calculation. Furthermore, this enables us to use relatively small basis-sets (unc-aug-ccpwCTZ) in the optimization since the HF-level complete basis-set (CBS) values are reached quickly. This methodology allows us to use CCSD(T) in the optimization in a computationally less costly manner while keeping essentially the same level of accuracy. More detailed discussion is given in the aforementioned work.

The quantum chemistry packages Molpro15 and PySCF16 were used for the evaluation of the objective functions throughout this work. The nonlinear optimization code Donlp217 was used to optimize the ccECP parameters. Most ccECP parameters in the optimization are initialized from either BFD12 or STU18–20 pseudopotentials.

In this section, we provide the atomic spectrum results and molecular transferability tests for all the elements mentioned above. For each element, we show the quality of the spectrum by tabulating the mean absolute deviation (MAD) of the ECP gaps with respect to AE values. The MAD is given as

MAD=1Nss=1NsΔEsECPΔEsAE.
(8)

For Ga–Kr, these quantities and errors are provided at the CBS limit, and for the other lighter atoms, these are given at the unc-aug-cc-pCV5Z basis set. In order to obtain the CBS limit energy, SCF and correlation energies are extrapolated separately. The HF energy extrapolation has the following form:

EnHF=ECBSHF+aexp(bn),
(9)

where n is the basis size, ECBSHF is the HF CBS limit value, and a and b are fit parameters. The correlation energy is extrapolated as

Encorr=ECBScorr+c(n+3/8)3+d(n+3/8)5,
(10)

where c and d are other fit parameters. The basis set used in extrapolation is uncontracted aug-cc-pwCVnZ,21,22 where n = (T, Q, 5).

We also evaluated the errors of other established ECPs so as to demonstrate the degree of improvement in our ccECP pseudopotentials. The ECPs considered are BFD,12 YSTU (Y = {SDF,18 MWB,19 MDF20}), SBKJC,23,24 eCEPP,8 and CRENBL.25,26 We also present the calculations for an AE uncorrelated-core (UC), where there are only VV correlations considered and the core is of the same size as in the ECPs. The same uncontracted basis sets and extrapolation schemes are used for these core approximations for consistency. We note that only VV correlations are present in our ccECP calculations. In addition, we want to emphasize that we did not include core polarization/relaxation (CPP) terms into ccECPs, neither considered these when using the other existing tables. One reason for this choice has been our primary goal to probe for the limits of accuracy of the very simplest ECP form. Additional reason has been our aim for the widest possible use and the fact that in majority of condensed/periodic matter codes such terms are not currently implemented. Clearly, this can be addressed in future studies.

For each element, the quality/transferability of the ccECP is tested by comparing the binding energy of the ECP for the selected molecules to the fully correlated AE case which utilizes relativistic DKH Hamiltonian. The binding energy is calculated for a variety of molecular geometries, including for compressed geometries up to the dissociation bond length, having in mind potential high-pressure applications. The discrepancy Δ(r) between the ECP and AE is plotted for all calculated molecules,

Δ(r)=De(r)ECPDe(r)DKHAE.
(11)

We fit the calculated binding energies from multiple bond lengths to Morse potential,

V(r)=Dee2arre2earre,
(12)

where De is the dissociation energy, re is the equilibrium bond length, and the parameter, a, is used to calculate the vibrational frequency,

ωe=2a2Deμ,
(13)

where μ is the reduced mass of the molecule. Fitted De, re, ωe, and Ddiss (binding energy at the dissociation threshold at the short bond lengths) are compared to AE counterparts to assess the transferability of an ECP for a particular molecule.

We consider the monohydride and monoxide molecules as our test systems to probe the σ and π bonds. In some cases, the homonuclear dimer was also used if it formed the same type of bond as in the above cases. Other ECP tests are included as well for comparison. Each oxygen ECP is taken from the corresponding ECP table, and our ccECP hydrogen is used throughout all hydride tests.

All transferability tests were carried out at unc-aug-cc-pwCV5Z basis sets unless otherwise specified. From limited tests on SeO, we observed that no significant changes occur in the discrepancy of binding energies Δ(r) compared to the CBS limit (see the supplementary material).

In Table IV, we show a summary of spectrum discrepancies of all examined ECPs compared to the AE spectrum. The MAD of the entire spectrum and the MAD of some selected low-lying states (LMAD) are given, while for Li and Be, we give the MAD values only. MDFSTU and MWBSTU are combined in one column since they are mutually exclusive for each of the atoms considered in this work. Only available ECP MADs are listed with the rest left blank. The calculated AE spectrum and ECP discrepancies for each atom can be found in the supplementary material. A summary graph of MADs and LMADs for 3rd-row main group elements is shown in Fig. 1.

FIG. 1.

MAD and LMAD for the 3rd-row main group elements. For each atom, we provide MAD of the entire calculated spectrum and MAD of the selected low-lying atomic states (LMAD). The LMAD atomic states include electron affinity, first, and second ionizations only. The solid lines represent MAD, and the dashed lines represent LMAD. The shaded region indicates the band of chemical accuracy.

FIG. 1.

MAD and LMAD for the 3rd-row main group elements. For each atom, we provide MAD of the entire calculated spectrum and MAD of the selected low-lying atomic states (LMAD). The LMAD atomic states include electron affinity, first, and second ionizations only. The solid lines represent MAD, and the dashed lines represent LMAD. The shaded region indicates the band of chemical accuracy.

Close modal

All optimized ccECP parameters are provided in Tables I and II in semilocal form. These tables assume the following common format of the potential:

V(r)=kβkrnk2eαkr2,
(14)

where is the angular momentum channel. Note that the Zeffr term is not tabulated since it is added by default in most codes/packages.

TABLE I.

Parameter values for 3rd-row main group element ccECPs. For all ECPs, the highest value corresponds to the local channel loc. Note that the highest nonlocal angular momentum channel max is related to it as max = loc − 1.

AtomZeffnℓkαℓkβℓkAtomZeffnℓkαℓkβℓk
11.442 508 11.866 873 As 3.479 387 75.655 194 
  6.537 124 90.076 771   1.637 480 −3.311 453 
  9.631 219 11.534 202   3.229 364 67.961 867 
  4.508 811 27.720 235   1.666 366 −3.094 558 
  7.273 863 9.000 000   2.068 163 24.304 734 
  11.172 983 65.464 770   1.546 999 0.939 456 
  7.706 175 −10.844 336   1.285 931 5.000 000 
  5.624 917 −15.963 161   9.934 874 6.429 657 
        1.895 682 −15.012 439 
        1.728 256 2.898 814 
Ca 10 11.240 167 149.302 623 Se 4.175 363 71.379 280 
  5.353 612 23.759 329   2.144 911 0.426 199 
  13.066 548 99.204 114   4.287 722 50.948 290 
  4.027 485 13.452 161   2.095 383 5.542 881 
  7.041 332 10.000 000   1.394 037 6.204 697 
  14.014 449 70.413 317   1.696 599 0.533 957 
  13.769 362 −92.872 980   2.977 052 6.000 000 
  4.717 260 −5.753 568   7.016 674 17.862 311 
        3.960 663 −20.009 132 
        5.028 263 10.005 735 
Ga 1.857 811 21.789 310 Br 4.971 807 85.88 4347 
  0.919 506 −2.866 851   2.042 687 4.621 255 
  1.920 302 18.639 860   4.711 839 55.361 715 
  1.008 959 −1.633 697   2.384 293 11.031 410 
  0.627 509 2.035 237   3.412 863 26.410 410 
  0.326 190 −0.085 324   1.530 285 5.468 739 
  17.004 739 3.000 000   3.665 770 7.000 000 
  14.999 618 51.014 218   5.293 023 25.660 393 
  11.992 792 −39.000 626   3.176 376 13.040 262 
  14.992 823 35.446 594   2.897 544 −21.908 839 
Ge 2.894 474 43.265 429 Kr 5.490 729 92.889 552 
  1.550 340 −1.909 340   3.863 012 12.929 478 
  2.986 529 35.263 014   4.038 577 43.099 524 
  1.283 381 0.963 440   3.306 789 9.509 760 
  1.043 001 2.339 019   4.213 480 17.804 945 
  0.554 563 0.541 381   1.549 897 4.589 115 
  1.478 963 4.000 000   10.794 238 8.000 000 
  3.188 906 5.915 851   13.323 389 86.353 904 
  1.927 439 −12.033 713   9.292 050 −11.114 533 
  1.545 539 1.283 543   20.148 958 10.229 519 
AtomZeffnℓkαℓkβℓkAtomZeffnℓkαℓkβℓk
11.442 508 11.866 873 As 3.479 387 75.655 194 
  6.537 124 90.076 771   1.637 480 −3.311 453 
  9.631 219 11.534 202   3.229 364 67.961 867 
  4.508 811 27.720 235   1.666 366 −3.094 558 
  7.273 863 9.000 000   2.068 163 24.304 734 
  11.172 983 65.464 770   1.546 999 0.939 456 
  7.706 175 −10.844 336   1.285 931 5.000 000 
  5.624 917 −15.963 161   9.934 874 6.429 657 
        1.895 682 −15.012 439 
        1.728 256 2.898 814 
Ca 10 11.240 167 149.302 623 Se 4.175 363 71.379 280 
  5.353 612 23.759 329   2.144 911 0.426 199 
  13.066 548 99.204 114   4.287 722 50.948 290 
  4.027 485 13.452 161   2.095 383 5.542 881 
  7.041 332 10.000 000   1.394 037 6.204 697 
  14.014 449 70.413 317   1.696 599 0.533 957 
  13.769 362 −92.872 980   2.977 052 6.000 000 
  4.717 260 −5.753 568   7.016 674 17.862 311 
        3.960 663 −20.009 132 
        5.028 263 10.005 735 
Ga 1.857 811 21.789 310 Br 4.971 807 85.88 4347 
  0.919 506 −2.866 851   2.042 687 4.621 255 
  1.920 302 18.639 860   4.711 839 55.361 715 
  1.008 959 −1.633 697   2.384 293 11.031 410 
  0.627 509 2.035 237   3.412 863 26.410 410 
  0.326 190 −0.085 324   1.530 285 5.468 739 
  17.004 739 3.000 000   3.665 770 7.000 000 
  14.999 618 51.014 218   5.293 023 25.660 393 
  11.992 792 −39.000 626   3.176 376 13.040 262 
  14.992 823 35.446 594   2.897 544 −21.908 839 
Ge 2.894 474 43.265 429 Kr 5.490 729 92.889 552 
  1.550 340 −1.909 340   3.863 012 12.929 478 
  2.986 529 35.263 014   4.038 577 43.099 524 
  1.283 381 0.963 440   3.306 789 9.509 760 
  1.043 001 2.339 019   4.213 480 17.804 945 
  0.554 563 0.541 381   1.549 897 4.589 115 
  1.478 963 4.000 000   10.794 238 8.000 000 
  3.188 906 5.915 851   13.323 389 86.353 904 
  1.927 439 −12.033 713   9.292 050 −11.114 533 
  1.545 539 1.283 543   20.148 958 10.229 519 
TABLE II.

Parameter values for additional ccECPs. For all ECPs, the highest value corresponds to the local channel loc. Note that the highest nonlocal angular momentum channel max is related to it as max = loc − 1. ECPs indicated with “reg” do not remove any electrons and contain only local channel.

AtomZeffnℓkαℓkβℓk
1.000 000 0.000 000 
  21.243 595 1.000 000 
  21.243 595 21.243 595 
  21.776 967 −10.851 924 
He 1.000 000 0.000 000 
  32.000 000 2.000 000 
  32.000 000 64.000 000 
  33.713 355 −27.700 840 
Li(reg) 10.000 00 0.000 000 
  24.000 00 3.000 000 
  24.000 00 72.000 000 
  24.000 00 −36.443 940 
  9.136 790 0.780 200 
Be(reg) 10.000 000 0.000 000 
  18.441 179 4.000 000 
  22.941 215 73.764 716 
  30.142 700 −52.330 20 
  11.327 931 −3.360 040 
Li 1.330 248 6.752 868 
  15.000 000 1.000 000 
  15.047 997 15.000 000 
  1.806 054 −1.242 730 
Be 2.487 404 12.663 919 
  17.949 002 2.000 000 
  24.132 003 35.898 004 
  20.138 003 −12.774 998 
  4.333 171 −2.960 014 
14.780 765 78.901 772 
  12.087 585 7.000 000 
  12.838 063 84.613 094 
  12.312 346 −53.027 517 
Ne 16.554 415 81.622 057 
  14.793 512 8.000 000 
  16.582 039 118.348 096 
  16.080 735 −70.278 859 
AtomZeffnℓkαℓkβℓk
1.000 000 0.000 000 
  21.243 595 1.000 000 
  21.243 595 21.243 595 
  21.776 967 −10.851 924 
He 1.000 000 0.000 000 
  32.000 000 2.000 000 
  32.000 000 64.000 000 
  33.713 355 −27.700 840 
Li(reg) 10.000 00 0.000 000 
  24.000 00 3.000 000 
  24.000 00 72.000 000 
  24.000 00 −36.443 940 
  9.136 790 0.780 200 
Be(reg) 10.000 000 0.000 000 
  18.441 179 4.000 000 
  22.941 215 73.764 716 
  30.142 700 −52.330 20 
  11.327 931 −3.360 040 
Li 1.330 248 6.752 868 
  15.000 000 1.000 000 
  15.047 997 15.000 000 
  1.806 054 −1.242 730 
Be 2.487 404 12.663 919 
  17.949 002 2.000 000 
  24.132 003 35.898 004 
  20.138 003 −12.774 998 
  4.333 171 −2.960 014 
14.780 765 78.901 772 
  12.087 585 7.000 000 
  12.838 063 84.613 094 
  12.312 346 −53.027 517 
Ne 16.554 415 81.622 057 
  14.793 512 8.000 000 
  16.582 039 118.348 096 
  16.080 735 −70.278 859 

Detailed discussions about each element and the accuracy of the ccECPs are given in Secs. III A–III C. The core radii of ccECPs are tabulated in Table III, and valence basis sets are discussed in Sec. III D.

TABLE III.

Core radii for all ccECPs given in this work. For every channel, radial distance is taken to be the distance from the origin at which the angular momentum channel agrees with the bare coulomb potential to within 10−5 Ha. Rc is the maximum of all channels. Radii are given in angstrom.

Local channel includedNonlocal channel
AtomrsrprdrfRcrsrprdRc
0.81 0.96 0.84  0.96 0.83 0.96  0.96 
Ca 0.78 0.98 0.88  0.98 0.87 0.99  0.99 
Ga 1.95 1.82 2.78 0.59 2.78 1.95 1.82 2.78 2.78 
Ge 1.45 1.60 2.34 1.51 2.34 1.48 1.58 2.34 2.34 
As 1.61 1.61 1.62 1.61 1.62 1.47 1.46 1.46 1.47 
Se 1.18 1.33 1.64 1.08 1.64 1.18 1.33 1.64 1.64 
Br 1.33 1.27 1.55 1.18 1.58 1.33 1.28 1.55 1.55 
Kr 1.01 1.08 1.53 0.65 1.53 1.01 1.08 1.53 1.53 
0.42    0.42     
He 0.36    0.36     
Li(reg) 0.59    0.59     
Be(reg) 0.56    0.56     
Li 1.68 1.34   1.68 1.68   1.68 
Be 1.25 0.90   1.25 1.25   1.25 
0.56 0.55   0.56 0.55   0.55 
Ne 0.52 0.51   0.52 0.52   0.52 
Local channel includedNonlocal channel
AtomrsrprdrfRcrsrprdRc
0.81 0.96 0.84  0.96 0.83 0.96  0.96 
Ca 0.78 0.98 0.88  0.98 0.87 0.99  0.99 
Ga 1.95 1.82 2.78 0.59 2.78 1.95 1.82 2.78 2.78 
Ge 1.45 1.60 2.34 1.51 2.34 1.48 1.58 2.34 2.34 
As 1.61 1.61 1.62 1.61 1.62 1.47 1.46 1.46 1.47 
Se 1.18 1.33 1.64 1.08 1.64 1.18 1.33 1.64 1.64 
Br 1.33 1.27 1.55 1.18 1.58 1.33 1.28 1.55 1.55 
Kr 1.01 1.08 1.53 0.65 1.53 1.01 1.08 1.53 1.53 
0.42    0.42     
He 0.36    0.36     
Li(reg) 0.59    0.59     
Be(reg) 0.56    0.56     
Li 1.68 1.34   1.68 1.68   1.68 
Be 1.25 0.90   1.25 1.25   1.25 
0.56 0.55   0.56 0.55   0.55 
Ne 0.52 0.51   0.52 0.52   0.52 

In this subsection, we explain the results for K and Ca atoms. These atoms use Ne-core ECPs, and they involve only s and p nonlocal projectors.

Table IV shows the MAD and LMAD of the spectrum discrepancies of various core approximations for the K atom. The chosen spectrum includes neutral d excitations, electron affinity state, and all ionizations until [Ne]3s2. Here, the least MAD corresponds to MDFSTU and the least LMAD corresponds to BFD. Figure 2 shows binding energy discrepancy Δ(r) of different ECPs for KH and KO molecules. The shaded region indicates the chemical accuracy of 0.043 eV, and the vertical dashed line marks the equilibrium bond length. Although all ECPs are within chemical accuracy, the ccECP is closest to zero and flatter near the KO molecule equilibrium geometry. Note that a flat discrepancy curve will result in a better vibrational frequency from the fit. It should be noted that the CRENBL ECP also produces very accurate molecular properties for KH and KO; however, it has singularities in both the attractive/local and repulsive nonlocal channels. Our ccECP has a reasonable balance of MAD and LMAD but arguably is more accurate in transferability, especially near the equilibrium region.

TABLE IV.

A summary of the spectrum discrepancies of various ECPs for the 3rd-row main group elements. For each atom, we provide MAD of the entire calculated spectrum and MAD of the selected low-lying atomic states (LMAD). The LMAD atomic states include electron affinity, first, and second ionizations only. All values are in eV.

AtomQuantityUCBFDSDFSTUMDFSTU/MWBSTUSBKJCCRENBLccECP
MAD 0.110 808 0.174 431  0.057 621  0.123 791 0.109 417 
 LMAD 0.004 905 0.009 807  0.012 360  0.015 636 0.017 944 
Ca MAD 0.147 679 0.155 381  0.098 974  0.199 446 0.047 605 
 LMAD 0.001 868 0.004 282  0.026 807  0.002 046 0.024 634 
Ga MAD 0.400 024 0.317 379 0.302 297 0.368 861   0.283 783 
 LMAD 0.273 255 0.239 687 0.237 425 0.304 068   0.225 519 
Ge MAD 0.542 077 0.497 299 0.487 997 0.523 465 0.414 607  0.468 781 
 LMAD 0.111 700 0.096 646 0.068 051 0.092 133 0.110 203  0.071 163 
As MAD 0.599 997 0.513 494 0.452 386 1.000 536 0.490 010  0.277 827 
 LMAD 0.098 619 0.068 839 0.038 426 0.309 579 0.094 349  0.011 692 
Se MAD 0.608 203 0.448 970 0.332 396 0.466 401 0.470 802  0.307 693 
 LMAD 0.086 533 0.100 668 0.027 611 0.087 471 0.119 486  0.018 102 
Br MAD 0.594 315 0.471 870 0.572 803 0.517 565 0.532 889  0.346 153 
 LMAD 0.084 249 0.063 693 0.108 582 0.115 210 0.123 326  0.036 020 
Kr MAD 0.755 949 0.440 971  0.542 141 0.716 075  0.194 395 
 LMAD 0.099 244 0.066 169  0.121 735 0.137 783  0.046 252 
AtomQuantityUCBFDSDFSTUMDFSTU/MWBSTUSBKJCCRENBLccECP
MAD 0.110 808 0.174 431  0.057 621  0.123 791 0.109 417 
 LMAD 0.004 905 0.009 807  0.012 360  0.015 636 0.017 944 
Ca MAD 0.147 679 0.155 381  0.098 974  0.199 446 0.047 605 
 LMAD 0.001 868 0.004 282  0.026 807  0.002 046 0.024 634 
Ga MAD 0.400 024 0.317 379 0.302 297 0.368 861   0.283 783 
 LMAD 0.273 255 0.239 687 0.237 425 0.304 068   0.225 519 
Ge MAD 0.542 077 0.497 299 0.487 997 0.523 465 0.414 607  0.468 781 
 LMAD 0.111 700 0.096 646 0.068 051 0.092 133 0.110 203  0.071 163 
As MAD 0.599 997 0.513 494 0.452 386 1.000 536 0.490 010  0.277 827 
 LMAD 0.098 619 0.068 839 0.038 426 0.309 579 0.094 349  0.011 692 
Se MAD 0.608 203 0.448 970 0.332 396 0.466 401 0.470 802  0.307 693 
 LMAD 0.086 533 0.100 668 0.027 611 0.087 471 0.119 486  0.018 102 
Br MAD 0.594 315 0.471 870 0.572 803 0.517 565 0.532 889  0.346 153 
 LMAD 0.084 249 0.063 693 0.108 582 0.115 210 0.123 326  0.036 020 
Kr MAD 0.755 949 0.440 971  0.542 141 0.716 075  0.194 395 
 LMAD 0.099 244 0.066 169  0.121 735 0.137 783  0.046 252 
FIG. 2.

Binding energy discrepancies for (a) KH and (b) KO molecules with regard to CCSD(T). The shaded region indicates the band of chemical accuracy. The dashed vertical line represents the equilibrium geometry.

FIG. 2.

Binding energy discrepancies for (a) KH and (b) KO molecules with regard to CCSD(T). The shaded region indicates the band of chemical accuracy. The dashed vertical line represents the equilibrium geometry.

Close modal

Atomic data for the Ca atom are given in Table IV as well. Here, the ccECP has a significantly lower MAD in the spectrum compared to all other core approximations. The lowest LMAD corresponds to CRENBL, although its MAD is still significantly higher than the ccECP. Binding energy discrepancies for CaH and CaO are shown in Fig. 3. We see that in CaO, most core approximations underbind as the bond is compressed, whereas the ccECP has minimal errors throughout all geometries.

FIG. 3.

Binding energy discrepancies for (a) CaH and (b) CaO molecules. The rest of the notation is the same as for KH and KO systems.

FIG. 3.

Binding energy discrepancies for (a) CaH and (b) CaO molecules. The rest of the notation is the same as for KH and KO systems.

Close modal

In Table V, we give the MAD of all fitted parameters of the Morse potential in Eq. (12) for K and Ca molecules (KH, KO, CaH, and CaO) as a summary of transferability. We find that the ccECP results in almost the same transferability as CRENBL with a minor exception in ωe. It is again worthwhile to mention that CRENBL, SBKJC, and STU ECP parameterizations differ from BFD, eCEPP, and ccECP form where the Coulomb and other singularities are absent by construction. Therefore, we conjecture that the accuracy achieved in the ccECP is close to the best possible fit within the given form, and additional optimized Gaussians will likely result in only mild and diminishing improvements.

TABLE V.

Mean absolute deviations of binding parameters for various core approximations with respect to AE data for K, Ca hydride, and oxide molecules. All parameters were obtained using the Morse potential fit. The parameters shown are dissociation energy De, equilibrium bond length re, vibrational frequency ωe, and binding energy discrepancy at dissociation bond length Ddiss.

De (eV)re (Å)ωe (cm−1)Ddiss (eV)
UC 0.007(2) 0.002(1) 1.2(2.3) 0.02(3) 
BFD 0.021(2) 0.003(1) 2.8(2.3) 0.05(3) 
CRENBL 0.007(2) 0.001(1) 0.6(2.3) 0.02(3) 
MDFSTU 0.030(2) 0.002(1) 2.6(2.3) 0.03(3) 
ccECP 0.007(2) 0.001(1) 1.9(2.3) 0.02(3) 
De (eV)re (Å)ωe (cm−1)Ddiss (eV)
UC 0.007(2) 0.002(1) 1.2(2.3) 0.02(3) 
BFD 0.021(2) 0.003(1) 2.8(2.3) 0.05(3) 
CRENBL 0.007(2) 0.001(1) 0.6(2.3) 0.02(3) 
MDFSTU 0.030(2) 0.002(1) 2.6(2.3) 0.03(3) 
ccECP 0.007(2) 0.001(1) 1.9(2.3) 0.02(3) 

In summary, we see that due to the valence space being large in K and Ca, chemical accuracy can be achieved for molecules and for a large part of the atomic spectrum.

This subsection discusses the results for 4p elements: Ga–Kr. These atoms use large [Ar]3d10 core ECPs with s, p, and d nonlocal channel operators. In this set of atoms, the valence space (4s and 4p) is smaller and the energy scales are significantly higher, resulting in larger errors in both spectrum and molecular binding energies. This is especially pronounced for higher ionization states and in compressed bond regions.

Table IV shows the MADs and LMADs of the considered spectrum for Ga–Kr atoms. For each element, we choose a consistent set of excitations/gaps relative to the ground state, which represents the spectrum. The spectrum consists of electron affinity, all ionization potentials down to single electron in the valence space, 4s → 4p excited states, and various 4d occupied states.

In all cases, our ccECP has lower MAD and LMAD compared to all core approximations including the uncorrelated-core (UC). The only exception to this is Ge where SBKJC achieves lower MAD although it is higher in LMAD, while SDFSTU achieves lower LMAD but is higher in MAD. We believe the ccECP is perhaps the most reasonable compromise between these two; this becomes clearer by looking at molecular discrepancies (Fig. 5). Speaking in general terms, significant gains in accuracy are obtained for As, Se, Br, and Kr elements, while Ga and Ge show modest improvements or are on par with previously published constructions. However, we are pleased to see more consistent improvements in molecular properties in our ccECPs as follows.

We show the transferability results for these elements in Figs. 4–9 which are hydride and oxide molecules for Ga, Ge, As, Se, Br, and Kr atoms, respectively. For the GaO molecule, we omit the UC molecular binding discrepancy curve where we were unable to converge the coupled cluster molecular calculations correctly. For Kr, we calculate the KrH+ discrepancies only. Similarly to the picture seen in the MAD of the spectrum, the ccECP has the lowest errors in molecular properties when compared to all other core approximations. We observe improvements in Ga and Ge molecules, while As molecular properties show more notable gains in accuracy. Significant improvements are seen in Se, Br, and Kr where the ccECP is essentially always within chemical accuracy for both hydrides and oxides in all bond lengths up to the dissociation limit. The summary of MADs of Morse potential fit parameters is shown in Table VI. We can see that our ccECP has the smallest errors on De, re, and Ddiss when compared to other ECPs. Only BFD has a smaller error in ωe; however, this result is not fully consistent with errors for the other parameters being notably larger.

FIG. 4.

Binding energy discrepancies for (a) GaH and (b) GaO molecules. The same notation applies as for the previous cases.

FIG. 4.

Binding energy discrepancies for (a) GaH and (b) GaO molecules. The same notation applies as for the previous cases.

Close modal
FIG. 5.

Binding energy discrepancies for (a) GeH and (b) GeO molecules. The same notation applies as for the previous cases.

FIG. 5.

Binding energy discrepancies for (a) GeH and (b) GeO molecules. The same notation applies as for the previous cases.

Close modal
FIG. 6.

Binding energy discrepancies for (a) AsH and (b) AsO molecules. The same notation applies as for the previous cases.

FIG. 6.

Binding energy discrepancies for (a) AsH and (b) AsO molecules. The same notation applies as for the previous cases.

Close modal
FIG. 7.

Binding energy discrepancies for (a) SeH and (b) SeO molecules. The same notation applies as for the previous cases.

FIG. 7.

Binding energy discrepancies for (a) SeH and (b) SeO molecules. The same notation applies as for the previous cases.

Close modal
FIG. 8.

Binding energy discrepancies for (a) BrH and (b) BrO molecules. The same notation applies as for the previous cases.

FIG. 8.

Binding energy discrepancies for (a) BrH and (b) BrO molecules. The same notation applies as for the previous cases.

Close modal
FIG. 9.

Binding energy discrepancy for the KrH+ molecule. The binding curves are relative to the CCSD(T) binding curve. The shaded region indicates a discrepancy of chemical accuracy in either direction.

FIG. 9.

Binding energy discrepancy for the KrH+ molecule. The binding curves are relative to the CCSD(T) binding curve. The shaded region indicates a discrepancy of chemical accuracy in either direction.

Close modal
TABLE VI.

Mean absolute deviations of binding parameters for various core approximations with respect to AE data for Ga-Kr hydride and oxide molecules. All parameters were obtained using the Morse potential fit. The parameters shown are dissociation energy De, equilibrium bond length re, vibrational frequency ωe, and binding energy discrepancy at dissociation bond length Ddiss.

De (eV)re (Å)ωe (cm−1)Ddiss (eV)
UC 0.028(3) 0.0092(7) 11(4) 0.26(3) 
BFD 0.093(4) 0.0091(8) 9(4) 0.34(3) 
MDFSTU 0.066(3) 0.0090(8) 25(3) 0.27(3) 
SBKJC 0.034(3) 0.0052(6) 22(3) 0.24(3) 
SDFSTU 0.087(3) 0.0097(8) 25(3) 0.30(3) 
ccECP 0.022(3) 0.0030(8) 17(3) 0.11(3) 
De (eV)re (Å)ωe (cm−1)Ddiss (eV)
UC 0.028(3) 0.0092(7) 11(4) 0.26(3) 
BFD 0.093(4) 0.0091(8) 9(4) 0.34(3) 
MDFSTU 0.066(3) 0.0090(8) 25(3) 0.27(3) 
SBKJC 0.034(3) 0.0052(6) 22(3) 0.24(3) 
SDFSTU 0.087(3) 0.0097(8) 25(3) 0.30(3) 
ccECP 0.022(3) 0.0030(8) 17(3) 0.11(3) 

In these 4p main group elements, we observe that the achieved ccECP MADs are significantly higher when compared to K and Ca. In addition, we were not able to achieve chemical accuracy in some compressed oxides although all ccECP hydrides are essentially within 1 kcal/mol. This is mainly due to the smaller size of the valence space. If higher accuracy is desired, a smaller core with a larger valence space will be required, but it is beyond the scope of this work.

In this section, we provide the data for several lighter elements that were not included in our previously published sets to complete our ccECP table for elements H through Kr. The elements considered here are H, He, Li, Be, F, and Ne.

For H, He, Li, and Be, we developed local pseudopotentials, where there are no nonlocal operators and no electrons are removed. These all-electron ECPs are used to cancel the Coulomb singularity which can be useful for computational efficiency purposes, thus resulting in soft, regularized potentials. We indicate these ECPs by label “reg” whenever there is ambiguity.

Although the H ccECP was published before,11 we provide the missing molecular tests here. We show the quality of regularized H and He ECPs by transferability tests on H2 and HeH+ molecules which are shown in Fig. 10. Note that the originally published BFD ECP significantly overbinds the H2 molecule due to a rather significant radial span of the regularized Coulomb potential.

FIG. 10.

Binding energy discrepancies for (a) H2 and (b) HeH+ molecules. The same notation applies as for the previous cases.

FIG. 10.

Binding energy discrepancies for (a) H2 and (b) HeH+ molecules. The same notation applies as for the previous cases.

Close modal

For Li and Be, we also developed He-core ECPs with an s nonlocal channel apart from all-electron pseudopotentials. The spectrum and molecular tests for Li and Be are given in Table VII and Figs. 11 and 12, respectively. The ccECPs have the most accurate spectrum for both elements among all core approximations. ccECP(reg) improves it significantly further by achieving MAD = 0.000 42 eV for Li and MAD = 0.001 95 eV for Be.

TABLE VII.

A summary of the spectrum discrepancies of various ECPs for the selected elements. For each atom, we provide MAD of the entire calculated spectrum and MAD of the selected low-lying atomic states (LMAD). The LMAD atomic states include electron affinity, first, and second ionizations only. All values are in eV.

AtomQuantityUCBFDSDFSTUMWBSTUeCEPPSBKJCCRENBLccECP
Li MAD 0.026 002 0.028 279 0.033 749  0.034 373 0.027 706 0.025 479 0.024 961 
Be MAD 0.039 119 0.039 855 0.039 200  0.041 993 0.038 804 0.043 981 0.032 178 
MAD 0.087 945 0.266 950 0.187 671 0.217 170 0.021 450 0.015 619 0.018 112 0.028 159 
 LMAD 0.028 785 0.080 417 0.025 233 0.033 469 0.007 550 0.005 818 0.005 222 0.004 284 
Ne MAD 0.124 251 0.320 135  0.212 618  0.043 077 0.036 869 0.023 995 
 LMAD 0.027 683 0.055 655  0.027 137  0.014 303 0.006 541 0.002 008 
AtomQuantityUCBFDSDFSTUMWBSTUeCEPPSBKJCCRENBLccECP
Li MAD 0.026 002 0.028 279 0.033 749  0.034 373 0.027 706 0.025 479 0.024 961 
Be MAD 0.039 119 0.039 855 0.039 200  0.041 993 0.038 804 0.043 981 0.032 178 
MAD 0.087 945 0.266 950 0.187 671 0.217 170 0.021 450 0.015 619 0.018 112 0.028 159 
 LMAD 0.028 785 0.080 417 0.025 233 0.033 469 0.007 550 0.005 818 0.005 222 0.004 284 
Ne MAD 0.124 251 0.320 135  0.212 618  0.043 077 0.036 869 0.023 995 
 LMAD 0.027 683 0.055 655  0.027 137  0.014 303 0.006 541 0.002 008 
FIG. 11.

Binding energy discrepancies for (a) Li2 and (b) LiO molecules. QZ basis and ccECP oxygen were used for all cases. Otherwise, the same notation applies as for the previous cases.

FIG. 11.

Binding energy discrepancies for (a) Li2 and (b) LiO molecules. QZ basis and ccECP oxygen were used for all cases. Otherwise, the same notation applies as for the previous cases.

Close modal
FIG. 12.

Binding energy discrepancies for (a) Be2 and (b) BeO molecules. QZ basis and ccECP oxygen were used for all cases. Otherwise, the same notation applies as for the previous cases.

FIG. 12.

Binding energy discrepancies for (a) Be2 and (b) BeO molecules. QZ basis and ccECP oxygen were used for all cases. Otherwise, the same notation applies as for the previous cases.

Close modal

Note that both Li and Be dimer curves are perfectly reproduced by ccECPs although oxides show considerable errors, as for all other ECPs. SDFSTU, however, is distinctly better near the dissociation region in LiO and BeO. We estimate that this is due to the ECP parameterization and the constraints of smoothness of our ECP form. If higher accuracy is desired, our ccECP(reg) is essentially flat and is always within 0.015 eV error in Li and Be dimers and oxides.

For F and Ne, we use He-core ECPs with an s nonlocal channel. The spectrum MADs for F and Ne atoms are shown in Table VII. In F, the ccECP has the lowest LMAD, although eCEPP, SBKJC, and CRENBL have lower MAD. However, they notably underbind in the F2 molecule, whereas the ccECP is close to the AE values, as shown in Fig. 13. In Ne, the ccECP achieves the lowest MAD and LMAD in the spectrum while achieving a more constant and overall smaller error in the NeH+ molecule, as shown in Fig. 14.

FIG. 13.

Binding energy discrepancies for (a) FH and (b) F2 molecules. The same notation applies as for the previous cases.

FIG. 13.

Binding energy discrepancies for (a) FH and (b) F2 molecules. The same notation applies as for the previous cases.

Close modal
FIG. 14.

Binding energy discrepancies for the NeH+ molecule. The same notation applies as for the previous cases.

FIG. 14.

Binding energy discrepancies for the NeH+ molecule. The same notation applies as for the previous cases.

Close modal

A summary of binding parameter errors for H, He, Li, Be, F, and Ne is given in Table VIII. Here, the ccECP has the lowest errors when compared to all the other ECPs. In fact, De and re show smaller errors than the UC, reflecting the impact of taking CC and CV correlations into account.

TABLE VIII.

Mean absolute deviations of binding parameters for various core approximations with respect to AE data for the H, He, Li, Be, F, and Ne related molecules. All parameters were obtained using the Morse potential fit. The parameters shown are dissociation energy De, equilibrium bond length re, vibrational frequency ωe, and binding energy discrepancy at dissociation bond length Ddiss.

De (eV)re (Å)ωe (cm−1)Ddiss (eV)
UC 0.018(4) 0.011(2) 6.2(8.8) 0.09(5) 
BFD 0.023(3) 0.005(1) 21(12) 0.15(5) 
CRENBL 0.027(3) 0.007(2) 15.7(8.8) 0.18(4) 
SBKJC 0.016(3) 0.006(2) 15.5(8.8) 0.20(4) 
SDFSTU 0.027(4) 0.006(2) 14.9(9.1) 0.12(5) 
eCEPP 0.047(4) 0.011(2) 30(13) 0.11(6) 
ccECP 0.010(3) 0.003(1) 12(11) 0.11(5) 
ccECP(reg) 0.006(4) 0.001(3) 0.8(4.1) 0.00(5) 
De (eV)re (Å)ωe (cm−1)Ddiss (eV)
UC 0.018(4) 0.011(2) 6.2(8.8) 0.09(5) 
BFD 0.023(3) 0.005(1) 21(12) 0.15(5) 
CRENBL 0.027(3) 0.007(2) 15.7(8.8) 0.18(4) 
SBKJC 0.016(3) 0.006(2) 15.5(8.8) 0.20(4) 
SDFSTU 0.027(4) 0.006(2) 14.9(9.1) 0.12(5) 
eCEPP 0.047(4) 0.011(2) 30(13) 0.11(6) 
ccECP 0.010(3) 0.003(1) 12(11) 0.11(5) 
ccECP(reg) 0.006(4) 0.001(3) 0.8(4.1) 0.00(5) 

For each ccECP pseudoatom, we generated the corresponding Gaussian correlation consistent basis sets. The provided basis sets are cc-pVnZ and aug-cc-pVnZ with n = (D, T, Q, 5, 6). The generation scheme we followed is similar to Ref. 12 which is briefly outlined here.

First, the HF ground state energy is minimized by optimizing a set of primitive exponents. The number of primitives is chosen in a way that the exact HF energy can be achieved with mHa accuracy or less. The exponents are optimized in an “even-tempered” way, where only the lowest/highest exponent and the ratio between the exponents are optimized. The obtained HF orbital coefficients are contracted for every occupied orbital. For K and Ca, an extra d contraction is also added by exciting an electron to the 3d orbital.

Then, m = n − 1 number of correlation consistent primitive exponents are added to every symmetry channel with contraction(s). Additional (n − 1, n − 2, n − 3, …, 1) correlation consistent polarization exponents are added to the next N = n − 1 higher symmetry channels without a contraction. For instance, Ge cc-pVQZ will have 1 s and 1 p contraction and additional (3s, 3p, 3d, 2f, 1g) correlation consistent primitive exponents. To optimize the correlation consistent exponents, the configuration-interaction with singles and doubles (CISD) ground state energy of the pseudoatom is minimized. In cases where there is only 1 valence electron present, the dimer CISD ground state energy is minimized. The number and the size of contractions are kept the same for all cardinal numbers n and correlation consistent exponents are also optimized in an even-tempered way, as described in Ref. 27. Additional diffuse primitive exponents are added by dividing the lowest primitive in each channel by 2.5 to obtain aug-cc-pVnZ. For K and Ca, m = n − 1 number of core correlating exponents are added by minimizing the atomic ground state to obtain (aug-)cc-pCVnZ. For ccECP(reg) basis sets, we use optimized correlation consistent exponents from Ref. 28. All ccECP basis sets in various code formats can be found at https://pseudopotentiallibrary.org.

In this work, we extended our ccECPs to the 3rd-row main group elements and provided additional ECPs that complete the ccECP set for all elements from H through Kr. The primary departure in the generation of ECPs from previous schemes has been the use of highly accurate many-body theories and methods in the construction. Another guide in the development was to carry out thorough tests particularly in compressed hydrides and oxides. These validation tests are imperative since one needs to know the inherent biases in the ECPs in order to be able to assess the corresponding biases in subsequent calculations.

Our ccECPs have substantially improved isospectral and transferable properties compared to previously tabulated ECPs. We observed that in certain cases such as compressed GaO and GeO, the effective core approximation breaks down. While the accuracy of ccECPs in those systems can be improved with more sophisticated parameterizations, the fundamental complication is the missing core-core classical electrostatic interactions which can be improved by considering a smaller core partitioning. Previous studies have discussed the complications in ECP construction due to these effects29–31 and the significant overlap of the core tails with the valence space.10,32

On the other hand, we observed that using an ECP is in fact quite beneficial compared to using an uncorrelated core instead. Note that, for instance, in 4p element molecules, UC almost always underbinds considerably and the ccECP is significantly more accurate. This effect can be seen in spectral discrepancies too, where the ccECP has smaller errors. The main reason for this can be understood if we recall that the electrons correlate more in an ECP than in an UC, i.e., the magnitude of ECP correlation energy is larger than that of the AE valence-valence counterpart (see, for instance, Ref. 33). This is because ECP orbitals have fewer nodes which result in higher probabilities for electrons to correlate. Therefore, our ccECPs capture part of CV and CC correlations resulting in better molecular properties. For this reason, the inclusion of total AE correlation energies (CC + CV + VV) is vitally important for proper reproduction of the valence properties.

We have shown that the accuracy and transferability of ECPs can be substantially improved by the inclusion of many-body effects. In fact, we believe that electron correlation has an important contribution to the final accuracy and additional parameters such as core polarization and spin-orbit terms are of secondary importance and can be implemented subsequently.

One achievement of this work is the demonstration of systems where accurate results are not feasible with conventionally used core approximations. For example, in LiO and BeO [Figs. 11(b) and 12(b)], all [He]-core approximations have significant errors, and only by including the missing 1s2 correlations, one is able to recover chemical accuracy. The improvements represented by ccECP(reg) can further expand the applicability such as for very high pressure/short bond length situations. As stated previously, having an upfront knowledge about the ECP biases is important and one can expect these errors to enlarge in solid/bulk calculations.

Another accomplishment of this work including previous papers is that it provides independent tests of previously tabulated ECPs. For instance, we observed that CRENBL and SBKJC ECPs show very high accuracy in molecules. SDFSTU pseudopotentials also demonstrate respectable accuracy results in few valence electron atoms. Independent tests of the ccECP have also been carried out recently,34 where the authors carry out the selected molecular tests from the G2 set using CCSD(T) and fixed-node diffusion Monte Carlo methods. In both cases, they find that ccECPs have lower MAD in atomization energies from experimental values when compared to other considered ECPs. Similarly, we believe that ccECPs achieve the best overall consistency in accuracy and we expect that the obtained results for hydrides and oxides will carry over to other systems and cases.

We anticipate that ccECPs will also be used in methods such as DFT that does not build upon explicitly correlated many-body wave functions. Since ccECPs reproduce the true many-body valence part of the spectrum, we believe that they should be used as ab initio valence Hamiltonians per se, even within such approaches. In addition, our ccECPs exhibit smooth densities in the core so that one can expect at least some alleviation of DFT biases caused by the high values and large gradients of the original core densities. However, one cannot expect the same type of systematic error from cores in all-electronic calculations to automatically transfer into ccECP valence calculations (and procedures similar to nonlinear core corrections35 might then be necessary). Furthermore, we do not expect much improvement in the errors that stem from genuine DFT shortcomings in the valence space (errors in atomizations, gaps, reaction paths, etc.). Clearly, further insights and calculations in this direction could be very helpful and useful.

Future related work is the calculation of accurate atomic ground state energies of ccECPs with various highly precise methods such as CCSDT(Q). These accurate energies will serve as a benchmark for independent calculations and methods such as fixed-node diffusion Monte Carlo. This will be published in a forthcoming paper.

The additional information about ccECPs is included in the supplementary material. Therein are given calculated AE spectra for each element and also corresponding discrepancies of various core approximations. AE, UC, and various ECP molecular fit parameters for hydrides, oxides, or dimers are provided. The ccECPs in semilocal and Kleinman-Bylander projected forms as well as optimized Gaussian valence basis sets in various input formats (Molpro, GAMESS, and NWChem) can be found at websites.36 

Input and output data generated in this work are published in the Materials Data Facility.37 

We are grateful to Paul R. C. Kent for the reading of the manuscript and helpful suggestions. M. Chandler Bennett contributed to the preparation of this manuscript and helped in formulating the methodologies used for generating the ECPs.

This work was 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.

This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.

This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

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

1.
D. M.
Ceperley
,
J. Stat. Phys.
43
,
815
(
1986
).
2.
B. L.
Hammond
,
P. J.
Reynolds
, and
W. A.
Lester
,
J. Chem. Phys.
87
,
1130
(
1987
).
3.
P. H.
Acioli
and
D. M.
Ceperley
,
J. Chem. Phys.
100
,
8169
(
1994
).
4.
L.
Maron
and
C.
Teichteil
,
Chem. Phys.
237
,
105
(
1998
).
5.
E.
Fromager
,
L.
Maron
,
C.
Teichteil
,
J.-L.
Heully
,
K.
Faegri
, and
K.
Dyall
,
J. Chem. Phys.
121
,
8687
(
2004
).
6.
J. R.
Trail
and
R. J.
Needs
,
J. Chem. Phys.
139
,
014101
(
2013
).
7.
J. R.
Trail
and
R. J.
Needs
,
J. Chem. Phys.
142
,
064110
(
2015
).
8.
J. R.
Trail
and
R. J.
Needs
,
J. Chem. Phys.
146
,
204107
(
2017
).
9.
M. C.
Bennett
,
C. A.
Melton
,
A.
Annaberdiyev
,
G.
Wang
,
L.
Shulenburger
, and
L.
Mitas
,
J. Chem. Phys.
147
,
224106
(
2017
).
10.
M. C.
Bennett
,
G.
Wang
,
A.
Annaberdiyev
,
C. A.
Melton
,
L.
Shulenburger
, and
L.
Mitas
,
J. Chem. Phys.
149
,
104108
(
2018
).
11.
A.
Annaberdiyev
,
G.
Wang
,
C. A.
Melton
,
M.
Chandler Bennett
,
L.
Shulenburger
, and
L.
Mitas
,
J. Chem. Phys.
149
,
134108
(
2018
).
12.
M.
Burkatzki
,
C.
Filippi
, and
M.
Dolg
,
J. Chem. Phys.
126
,
234105
(
2007
).
13.
M.
Born
and
R.
Oppenheimer
,
Ann. Phys.
389
,
457
(
1927
).
14.
M.
Reiher
and
A.
Wolf
,
J. Chem. Phys.
121
,
2037
(
2004
).
15.
H.-J.
Werner
,
P. J.
Knowles
,
G.
Knizia
,
F. R.
Manby
,
M.
Schütz
 et al, molpro, version 2019.1, a package of ab initio programs,
2019
, see http://www.molpro.net.
16.
Q.
Sun
,
T. C.
Berkelbach
,
N. S.
Blunt
,
G. H.
Booth
,
S.
Guo
,
Z.
Li
,
J.
Liu
,
J. D.
McClain
,
E. R.
Sayfutyarova
,
S.
Sharma
 et al, PySCF: The python-based simulations of chemistry framework,
2017
, https://onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1340.
17.
P.
Spellucci
, DONLP2 nonlinear optimization code,
2009
, https://www2.mathematik.tu-darmstadt.de/fbereiche/numerik/staff/spellucci/DONLP2/.
18.
G.
Igel-Mann
,
H.
Stoll
, and
H.
Preuss
,
Mol. Phys.
65
,
1321
(
1988
).
19.
T.
Leininger
,
A.
Berning
,
A.
Nicklass
,
H.
Stoll
,
H.-J.
Werner
, and
H.-J.
Flad
,
Chem. Phys.
217
,
19
(
1997
).
20.
H.
Stoll
,
B.
Metz
, and
M.
Dolg
,
J. Comput. Chem.
23
,
767
(
2002
).
21.
J. G.
Hill
and
K. A.
Peterson
,
J. Chem. Phys.
147
,
244106
(
2017
).
22.
N. J.
DeYonker
,
K. A.
Peterson
, and
A. K.
Wilson
,
J. Phys. Chem. A
111
,
11383
(
2007
).
23.
W. J.
Stevens
,
H.
Basch
, and
M.
Krauss
,
J. Chem. Phys.
81
,
6026
(
1984
).
24.
W. J.
Stevens
,
M.
Krauss
,
H.
Basch
, and
P. G.
Jasien
,
Can. J. Chem.
70
,
612
(
1992
).
25.
L.
Fernandez Pacios
and
P.
Christiansen
,
J. Chem. Phys.
82
,
2664
(
1985
).
26.
M.
Hurley
,
L. F.
Pacios
,
P.
Christiansen
,
R.
Ross
, and
W.
Ermler
,
J. Chem. Phys.
84
,
6840
(
1986
).
27.
T. H.
Dunning
, Jr.
,
K. A.
Peterson
, and
D. E.
Woon
, “Basis sets: Correlation consistent sets,” in , edited by
P.
von Ragué Schleyer
,
N. L.
Allinger
,
T.
Clark
,
J.
Gasteiger
,
P. A.
Kollman
,
H. F.
Schaefer
, and
P. R.
Schreiner
(
Wiley
,
2002
), p. 88.
28.
B. P.
Prascher
,
D. E.
Woon
,
K. A.
Peterson
,
T. H.
Dunning
, and
A. K.
Wilson
,
Theor. Chem. Acc.
128
,
69
(
2011
).
29.
P. J.
Hay
and
W. R.
Wadt
,
J. Chem. Phys.
82
,
299
(
1985
).
30.
P.
Schwerdtfeger
,
T.
Fischer
,
M.
Dolg
,
G.
Igel-Mann
,
A.
Nicklass
,
H.
Stoll
, and
A.
Haaland
,
J. Chem. Phys.
102
,
2050
(
1995
).
31.
T.
Leininger
,
A.
Nicklass
,
H.
Stoll
,
M.
Dolg
, and
P.
Schwerdtfeger
,
J. Chem. Phys.
105
,
1052
(
1996
).
32.
A. R.
Oganov
and
P. I.
Dorogokupets
,
Phys. Rev. B
67
,
224110
(
2003
).
33.
34.
T.
Wang
,
X.
Zhou
, and
F.
Wang
,
J. Phys. Chem. A
123
,
3809
(
2019
).
35.
S. G.
Louie
,
S.
Froyen
, and
M. L.
Cohen
,
Phys. Rev. B
26
,
1738
(
1982
).
36.
See https://pseudopotentiallibrary.org/ for Pseudopotential library.
37.
G.
Wang
,
A.
Annaberdiyev
,
C. A.
Melton
,
M. C.
Bennett
,
L.
Shulenburger
, and
L.
Mitas
, Dataset for “a new generation of effective core potentials from correlated calculations: 4s and 4p main group elements and first row additions.”

Supplementary Material