Optimized normconserving HartreeFock pseudopotentials for planewave calculations
Abstract
We report HartreeFock (HF) based pseudopotentials suitable for planewave calculations. Unlike typical effective core potentials, the present pseudopotentials are finite at the origin and exhibit rapid convergence in a planewave basis; the optimized pseudopotential method [A. M. Rappe et al. , Phys. Rev. B 41 1227–30 (1990)] improves planewave convergence. Normconserving HF pseudopotentials are found to develop longrange nonCoulombic behavior which does not decay faster than , and is nonlocal. This behavior, which stems from the nonlocality of the exchange potential, is remedied using a recently developed selfconsistent procedure [J. R. Trail and R. J. Needs, J. Chem. Phys. 122, 014112 (2005)]. The resulting pseudopotentials slightly violate the norm conservation of the core charge. We calculated several atomic properties using these pseudopotentials, and the results are in good agreement with allelectron HF values. The dissociation energies, equilibrium bond lengths, and frequency of vibrations of several dimers obtained with these HF pseudopotentials and plane waves are also in good agreement with allelectron results.
pacs:
02.70.Ss, 71.15.m, 31.25.vPresent address: ]Physics Department, Clark Hall, Cornell University, Ithaca NY 14850
I Introduction
In nearly all planewave density functional calculations kohn_nobel , the use of pseudopotentials is essential in order to eliminate the atomic core states and the strong potentials that bind them. This results in replacing the electronnuclear Coulomb interaction by a weaker (often angularmomentum dependent) potential. The advantage is the reduced computational cost, due to reduction in the number of electrons and in the planewave basis cutoff pickett .
In the physics community, pseudopotentials are typically constructed using density functional theory (DFT) due to its success in electronic structure calculations BHS ; RRKJ ; TM . Pseudopotentials based on methods other than DFT are less widely available. One of the key ingredients of most of these pseudopotentials is their “softness,” meaning that they provide rapid convergence in a planewave basis. On the other hand, in the chemistry community, HartreeFock (HF) pseudopotentials or effective core potentials (ECPs) are mostly used dolg_rev ; SBKJC ; hay_wadt_psp ; CLP_psp ; dolg_psp . Most of the available ECPs, which are typically expressed in a Gaussian basis, are of limited use in planewave calculations, mainly because of their singular behavior at the origin and slow convergence in a planewave basis.
Recently, there has been a growing interest in developing softcore HF pseudopotentials greeff_lester ; lester_psp ; trail_needs1 ; trail_needs2 . HF pseudopotentials have been applied in diffusion Monte Carlo calculations, and the results suggest that they are better suited for correlated calculations than DFTbased pseudopotentials greeff_lester . In addition, HF pseudopotentials could be useful in certain calculations where a negatively charged reference state is needed; HF tends to bind electrons more strongly than DFT. Also, hybrid exchangecorrelation functionals, which include some proportion of HF exact exchange (notably B3LYP B3LYP ) have proved successful in electronic structure calculations, and whether HF or B3LYP based pseudopotentials would be useful in these type of calculations remains to be explored.
Trail and Needs trail_needs1 have recently found that normconserving HF pseudopotentials constructed using standard normconserving pseudopotential methods develop a nondecaying and nonlocal tail. These pseudopotentials would generally lead to erroneous results, and even infinite energies in solid state calculations. Previously, pseudopotentials based on exact exchange within the optimized potential method (OPM) were also found to develop a similar, but local and eventually decaying long range structure bk_exact_xc ; stadele_exact_xc ; engel_exact_xc , which led to erroneous results in applications.
Different schemes have been proposed to cure the longrange behavior in the HF and OPM type pseudopotentials mentioned above. Originally with the OPM method, different groups developed procedures for reducing or eliminating the unphysical longrange behavior, while retaining desired eigenvalue matching and norm conservation of the pseudopotential bk_exact_xc ; stadele_exact_xc . This procedure was later implemented in a selfconsistent manner by altering the pseudoorbitals over all space, which led to a small violation of the norm conservation of the core charge. Nevertheless, the pseudopotentials proved transferable in both atomic and solidstate calculations engel_exact_xc . Trail and Needs fixed the longrange unphysical behavior using a similar selfconsistent approach trail_needs1 that also slightly violates norm conservation. Ovcharenko et al. lester_psp have recently developed HFbased pseudopotentials which are finite at the origin and are expressed in terms of Gaussian basis, but are still relatively hard for typical planewave calculations, as we show in Sec. IV.
In this paper, we construct softcore HF based pseudopotentials following the RappeRabeKaxirasJoannopoulos (RRKJ) method RRKJ . We use a selfconsistent approach to restore the Coulombic tail in a manner similar to that used by Trail and Needs trail_needs1 ; casino_psp . Pseudopotentials are highly nonunique entities, and it is useful to have several alternative forms available with different properties. For example, one of the advantages of the RRKJ construction scheme is the softness of the resulting pseudopotentials.
The rest of the paper is organized as follows. We first review the general construction scheme of the RRKJ pseudopotentials, and the procedure we used to “localize” the pseudopotentials (remove the long tail behavior). The crucial steps in the validation of the pseudopotentials are the study of their convergence and transferability, i. e. how well the pseudopotentials converge in a planewave basis set, and how well they perform in environments different than the reference configuration. These steps are explored in Sections III and IV. In Sec. III, we study planewave convergence of the atoms, and we investigate the transferability by looking at the ionization energy, electron affinity, and excitation energy of the first and secondrow elements. In Sec. IV, we study several dimers at the HF level, and compare both the allelectron and pseudopotential results. Finally, in Sec. V, we conclude by summarizing the main points.
Ii Pseudopotential construction
We follow the usual density functional theory approach for the construction of the normconserving pseudopotentials starting from an allelectron (AE) atomic calculation pickett . However, in this case, we solve the HartreeFock instead of the KohnSham equations. Our HF solver is adapted from the code of Ref. f_fisher, .
The pseudopotential construction starts by choosing an electronic reference state, and then solving the Schrödinger equation of the system for the eigenstates and eigenvalues . The wavefunction can be written as where are the spherical harmonics. The orbitals satisfy the radial Schrödinger equation:
(1) 
where , is the bare nuclear potential, and is the atomic number. is the HF potential defined such that,
(2) 
(3) 
(4) 
For convenience, equal and opposite self interaction terms are included in both the Hartree and exchange terms. Therefore, the Hartree potential is a functional of the spherically averaged total electronic density .
The second step in the construction of normconserving pseudopotentials typically begins with the generation of a smooth set of pseudoorbitals to replace the AE orbitals , such that equals the allelectron orbital beyond some cutoff distance . Different normconserving pseudopotentials mainly differ in the form of for , and the guiding principles used to construct .
In the RRKJ method, the pseudoorbitals in the core region are expressed as a linear combination of spherical Bessel functions
(5) 
The Bessel wavevectors are chosen such that,
(6) 
where and are the derivatives of and with respect to , respectively. is the number of Bessel coefficients, .
The key ingredient of the RRKJ method is the optimization of the Bessel coefficients, , to minimize the residual kinetic energy defined as follows:
(7)  
where is the Fourier transform of real space pseudoorbital , and is the target wavevector above which the contribution to the residual kinetic energy is as small as possible. For a particular , the planewave energy cutoff required to achieve total energy convergence of is given by the square of the maximum value of over all orbitals RRKJ .
The optimization of is constrained to pseudoorbitals that are continuous at , with continuous first and second derivatives [in fact, one of these local constraints would be already satisfied due to Eq.(6)], and that satisfy norm conservation of the core charge:
(8) 
This procedure defines all the Bessel coefficients in Eq. (5), and must be at least 3 to allow the constraints to be satisfied. Any additional parameters will be optimized by minimizing the residual kinetic energy due to the nonlinearity of Eq. (7). Throughout this study, we used .
The screened pseudopotential is defined as the potential that makes the desired pseudoorbital an eigenstate of the oneelectron Schrödinger equation, with the same eigenvalue as the corresponding allelectron state:
(9) 
This equation can now be inverted, thanks to the nodeless character of , to obtain the screened pseudopotential:
(10) 
Note that is a continuous function because of the continuity requirements imposed on the pseudoorbital and its first two derivatives.
The last step in the construction of pseudopotentials is to remove the screening effects of the valence electrons, and to obtain the ionic pseudopotential. This is done by subtracting the Hartree and exchangepotential contributions of the valence electrons from the screened potential . That is,
(11) 
where includes only the valence orbitals. Note that the last term in Eq. (11), the exchange potential, is orbital dependent, and its explicit construction is feasible because the pseudoorbitals are nodeless.
Up to this point, the pseudopotential construction recipe mirrors exactly the procedure done in regular normconserving pseudopotentials based on DFT, except that the original Hamiltonian is the HF and not the KS Hamiltonian. However, the resulting ionic potentials, , have different behavior for large , which in turn stems from the corresponding large behavior of the HF and KS orbitals. It was shown by Handeler et al. handler_limits_hf that the HF orbitals decay at infinity as
(12) 
where is determined by the eigenvalue of the highest occupied orbitals , and is an orbital dependent quantity handler_limits_hf . This behavior is to be contrasted with the exponential decaying behavior of the KS orbitals, for large .
In DFTbased pseudopotentials, the exchange potential is replaced by the exchangecorrelation functional of the electron density , which decays exponentially to zero for large . Thus, the ionic potential for large where is the valence charge density.
The effect of the special longrange behavior of the HF orbitals on the ionic pseudopotential can be understood by writing as,
(13)  
where the Schrödinger equation of Eq. (1) is used. Here is the total electronic charge density while is the pseudovalence charge density. In the HF case, the exchange potential is a nonlocal functional of the HF orbitals as shown in Eq. (4). In particular, the leading behavior of the numerator for large will be where , while the denominator will go as . Therefore, the exchange potential, as computed in Eq. (13), is not guaranteed to decay faster than . In fact, the exchange potential can even grow without bound for large when computed in this way. This is the source of the nonlocal longrange nonCoulombic tail in the HF pseudopotentials trail_needs1 .
Atom  

H  0.50  0.50  0.50  10.8 
He  0.60  0.60  0.60  12.8 
Li  2.19  2.37  2.37  2.7 
Be  1.41  1.41  1.41  6.2 
B  1.88  1.96  1.96  4.0 
C  1.10  1.10  1.10  8.3 
N  0.94  0.88  0.84  10.6 
O  0.80  0.75  0.99  12.6 
F  0.70  0.64  0.89  15.0 
Ne  0.63  0.57  0.63  17.1 
Na  2.70  2.85  2.85  2.5 
Mg  2.38  2.38  2.38  3.2 
Al  1.94  2.28  2.28  3.9 
Si  1.67  2.01  2.06  4.6 
P  1.48  1.71  1.71  5.8 
S  1.33  1.50  1.50  5.8 
Cl  1.19  1.34  1.34  6.5 
Ar  1.09  1.20  1.31  7.3 
In this work we follow the Trail and Needs selfconsistent approach to localize the ionic pseudopotential trail_needs1 . The localized ionic potential , which is forced to behave as for large , is defined as
(14) 
where is the core charge density, is the characteristic distance for the localized to approach , and is the localization radius. For all pseudopotentials, the value of for each channel is fixed to be the same as the corresponding . is defined as where and are parameters which are chosen to localize the ionic pseudopotential in a selfconsistent manner, and is a polynomial function whose explicit form depends on the pseudopotential construction method (see below).
The selfconsistent procedure is performed as follows. First, starting with an initial guess for the and parameters, the HF equations with a new ionic potential, as defined in Eq. (14), are solved for a new set of orbitals and eigenvalues. The and parameters are then adjusted to minimize the error in the eigenvalues and logarithmic derivativeslogdernote1 . The procedure is then repeated until and parameters are found that give eigenvalues and logarithmic derivatives that closely match the allelectron values to the desired tolerance.
Finally, we define used in . For the case of TroullierMartins pseudopotentials, we chose the same form as was used in Ref. trail_needs1 , namely for , and otherwise. For RRKJ pseudopotentials, we investigated a few different forms and found the results to be insensitive to the choices. We will report our results using for , and for .
AE  This work/RRKJ  This work/TM  TM^{1}^{1}1Ref. trail_needs1 .  

Atom  
H  0.5000  
He  0.8617  
Li  0.1963  
Be  0.2956  
B  0.2915  
C  0.3964  
N  0.5129  
O  0.4368  
F  0.5776  
Ne  0.7293  
Na  0.1819  
Mg  0.2428  
Al  0.2020  
Si  0.2812  
P  0.3690  
S  0.3317  
Cl  0.4335  
Ar  0.5430  
Average abs. error  
Maximum abs. error 
AE  RRKJ  TM  

Atom  
B  
C  
N  
O  
F  
Al  
Si  
P  
S  
Cl  
Average abs. error  
Maximum abs. error 
To illustrate the pseudopotential localization procedure, we show in the upper panel of Fig. 1 the pseudopotentials of Ne as constructed using the RRKJ method with the parameters given in Table I. In the lower panel, we show the same Ne pseudopotentials after the selfconsistent localization procedure. As can be seen, the differences between the upper and lower panels are rather small, and are in the core as well as in the valence regions. Most of the changes are within the core region, however. The insets show the longrange behavior of . Note, in particular, in the upper panel, the nondecaying tail of the potential, which can be shown trail_needs1 to behave as for large where and are real numbers.
In Fig. 2, we show in the upper panel the allelectron and pseudoorbitals after the localization procedure. The lower panel shows the difference in the pseudoorbitals before and after the localization procedure. Again, changes are found in both the core and valence regions. One obvious drawback of the localization procedure is that the corenorm will be different from the allelectron value, i.e. Eq. (8) will be violated slightly. This can be seen in Fig. 2 because the changes in the orbitals for are all positive, and will integrate to a nonzero value.
As mentioned before, the localization procedure conserves the allelectron eigenvalues, and the logarithmic derivative. In practice, the logarithmic derivative is only conserved approximately. Figures 3 and 4 show a systematic study of the violations of the normconservation and the change in the logarithmic derivative at due to the localization procedure for the first and secondrow elements. The study was done using both the RRKJ and TM pseudopotentials. From Fig. 3, we see that there is less than electron shift between the core and the valence. The relative error between the logarithmic derivatives before and after applying the localization procedure is shown in Fig. 4 logdernote2 . As can be seen from both figures, both pseudopotential construction methods yield similar results.
Atom  AE  RRKJ  TM  TM^{1}^{1}1Ref. trail_needs1 .  

Ground state  Excited state  
H  0.3750  
0.4444  
He  0.7302  
0.8061  
Li  0.0677  
0.1408  
Be  0.0615  
0.2389  
B  0.0784  
0.2353  
C  0.0894  
0.3402  
N  0.4127  
0.4565  
O  0.3809  
0.6255  
F  0.5220  
0.8781  
Ne  0.6734  
1.7565  
Na  0.0725  
0.1263  
Mg  0.0679  
0.1843  
Al  0.0858  
0.1441  
Si  0.0913  
0.2146  
P  0.3006  
0.3023  
S  0.2672  
0.4260  
Cl  0.3733  
0.5653  
Ar  0.4824  
1.1597  
Average abs. error  
Maximum abs. error 
Iii Pseudopotential Tests on atomic properties
In this section, we study the transferability of the pseudopotentials by looking at several atomic properties namely ionization energies, electron affinities, and excitation energies of the first and secondrow elements. All of these calculations are performed by solving the HF equations in real space with the same code opium which is used to generate the pseudopotentials.
We show results obtained using HF pseudopotentials constructed using both the RRKJ and TM methods. For the TM pseudopotentials, we used the same code opium which is used to generate the RRKJ pseudopotentials but following the TroullierMartins construction scheme.
The parameters for all pseudopotentials used in the study of the atomic properties are summarized in Table 1. In order to aid in comparison to previous results, the reference configurations and the construction parameters are the same as of Ref. trail_needs1 . In general, however, one has to choose these parameters according to the targeted applications. Moreover, unlike the TM method, the RRKJ method has two additional adjustable parameters for each pseudoorbital, namely and . We set for all pseudopotentials, and selected the ’s such that meV/electron for each orbital. As mentioned before, the energy cutoff required to achieve this level of energy convergence in the target calculations is approximately the square of the largest value used.
In Table 2 we present the ionization energies. For each element, we calculated the ionization energy using our RRKJ and TM pseudopotentials. These values are compared with the allelectron value (shown in the first column), and we only report the difference, , from the allelectron energy. is obtained using the original pseudopotential of Eq. (13) with the unphysical tail behavior, while is calculated using the localized HF pseudopotential of Eq. (14). We show also for comparison the same results as obtained using the HF TM pseudopotentials as reported in Ref. trail_needs1 .
For both pseudopotential construction schemes, the agreement between and indicates that the selfconsistent procedure that is used to localize the pseudopotential did not change significantly the original potential. It is important to stress that the effects of the nonCoulombic tail in atomic calculations are minimal, and this is why there is a good agreement between and . However, in any solid calculation with periodic boundary conditions, the results obtained using pseudopotentials with the tail problem would be erroneous.
Also, the pseudopotential results in Table 2 are in good agreement with the allelectron values for both the RRKJ and TM methods. On average the difference is less than milliHartree (mHa) and the largest deviation is less than mHa. This is a clear indication of the good quality of the pseudopotentials. Finally, the results obtained using the TM construction scheme are in excellent agreement with the equivalent results obtained by Trail and Needs trail_needs1 . Any small differences could be attributed to the different grids used, or the slight differences in the selfconsistent HF procedure.
In Table 3, we show a similar study for the electron affinity of the first and secondrow elements. The electron affinity is the difference in energy between the neutral atom and the negatively charged ion in their ground state configurations. The results agree to within mHa with previously published HF electron affinities EA .
Again, these results show that the pseudopotentials have good transferability properties. Also, we note here that the HF pseudopotentials generally bind an extra electron more strongly than the DFTbased pseudopotentials. This is why the study of the electron affinity is feasible.
We report the results of excitation energies in Table 4. This is the difference in energy between the ground state configuration shown in the first column and the two excited state configurations shown in the third column. As before, the pseudopotential results are in good agreement with the allelectron values. The absolute average deviations from the allelectron values are consistent across the different pseudopotential schemes, and is less than 1.2 mHa. The largest deviation is with carbon and is 5–6 mHa (the excitation energy between the and states).
Atom  

N  0.91  0.91  0.91  10.6 
P  1.58  1.58  1.90  5.5 
Cl  1.68  1.68  1.90  6.2 
E(Ry)  

N  P  
RRKJ  112  30 
TM  145  50 
OAL  210  55 
SBKJC  375  90 
LDA  HF  
AE  PSP  AE  PSP  
N  11.61  10.96  5.10  4.92  
2.07  2.06  2.01  2.02  
2388  2380  2720  2722  
P  6.25  6.02  1.70  1.70  
3.57  3.57  3.50  3.50  
794  788  912  905  
Cl  3.62  3.55  0.84  0.84  
3.74  3.74  3.73  3.72  
561  559  610  613 
Iv Study of several dimers
In this section, we apply our HF pseudopotentials in a study of three dimers and compare them with allelectron calculations. The pseudopotential calculations are done using a HF planewave basis code cpmd , and the allelectron results are obtained using a Gaussian basis code gaussian . For comparison, we have also generated LDA lda pseudopotentials using exactly the same parameters as those of the HF pseudopotential, and we compare them with their allelectron counterparts on an equal footing with the HF results. In our allelectron and pseudopotential results we used the PerdewWang lda flavor of LDA. We did not include a nonlinear core correction (NLCC)nl_core ; nl_core_porezag in our LDA pseudopotentials for LDA/HF comparison purposes. The pseudopotential DFT calculations are carried out using ABINIT Abinit . The allelectron calculations are performed using the correlationconsistent ccpV5Z basis set for N and ccpV(5+d)Z for P and Cl dunning . The differences between the quadrupole and quintuple basis sets are negligible; e.g., in the binding energies the differences are less than 0.02 eV. We also verified these values using 6311++G(3df,3pd) basis sets.
In Table 5, we show the cutoff radii we used to generate the RRKJ pseudopotentials for this study. These values, which are different from those of Table 1, are chosen to accommodate the bond lengths of these dimers without core overlap when possible. In the case of P and Cl, the core cutoff radius is allowed to extend beyond half of the dimer bond length to make softer pseudopotentials. The atomic reference configurations are the same as before.
Before presenting the molecular results, we comment on the planewave convergence of the RRKJ pseudopotentials in comparison to other pseudopotentials or ECPs. We examined the pseudopotentials of nitrogen and phosphorus as representatives of the first and secondrow elements. Figure 5 shows the RRKJ, TM, OAL lester_psp , and SBKJC EMSL pseudopotentials in real space for the two elements. We did not include a channel in these plots because neither the SBKJC nitrogen pseudopotential nor the OAL potentials have a channel. For the RRKJ and TM pseudopotentials, we used the cutoff radii as shown in Table 5. Note that the SBKJC ECP has a divergence at the origin, and is thus not suitable for planewave calculations. More generally, direct examination of pseudopotentials in real space gives little information about the planewave basis needed for converged eigenstates and eigenvalues.
One way to monitor the size of the planewave basis needed is to look at the residual kinetic energy convergence RRKJ of the pseudowavefunctions, which we show in Fig. 6. The residual kinetic energy is calculated as shown in Eq. (7) by varying (we plot it against which is the cutoff energy, for convenience). The planewave basis cutoff for each pseudopotential is determined by the potential channel requiring the largest cutoff energy. The estimated planewave basis cutoffs based on Fig. 6 are summarized in Table 6. As can be seen, RRKJ pseudopotentials give the smallest planewave basis cutoffs.
We summarize the results for the spectroscopic properties of the dimers in Table 7. The LDA and HF results for the equilibrium bond length and the harmonic vibrational frequencies of the dimers are shown to be well reproduced by both types of pseudopotentials. In the LDA dissociation energies there is a large discrepancy between allelectron and pseudopotential results for N, and to a lesser extent in P and even much less in Cl. The HF pseudopotential dissociation energy of N shows a significant deviation from the allelectron HF result, whereas for P and Cl, the results are the same to within eV. A similar difference between allelectron and pseudopotential results was also seen in Ref. engel_exact_xc using the OPM method.
Most of errors in the LDA results are due to linear descreening of the pseudopotential, which is typically repaired by including a NLCC in the target calculation. Furthermore, because all three dimers are closed shell, the linear descreening effects are largest in the isolated atom calculations and therefore mostly affect the dissociation energies. As shown by Porezag et al. nl_core_porezag , these errors are larger for first row species compared to heavier atoms and increase with more unpaired electrons, explaining why N has the largest error (first row, and 3 unpaired spins) and why Cl (second row, 1 unpaired spin) has the smallest. Our LDA allelectron and pseudopotential energies are in excellent agreement with previous results engel_exact_xc without a NLCC. This study also showed that the allelectron results are in good agreement with the pseudopotentials values after including a NLCC, indicating that this is the dominant error for LDA pseudopotentials.
At the HF level there is also an error in the exchange contribution to the pseudopotential since the terms arising from corevalence interactions in the atomic reference state are frozen into the core and only valencevalence exchange terms are subtracted during the descreening step. As shown by our results and others shirley ; hockengel , this type of descreening error is much smaller in magnitude in HF pseudopotentials compared to LDA. This would suggest that HF derived pseudopotentials are advantageous over those from LDA, at least in better descreening of the corevalence contributions.
V Summary and conclusion
Constructing HF pseudopotentials suitable for planewave calculations is a nontrivial task considering the extreme nonlocality of the exchange potential. Pseudopotentials constructed using a typical normconserving procedure develop a nonlocal and a nonCoulombic tail for large . This would generally lead to erroneous results, especially in solid calculations. Several schemes have been introduced to cure the long tail behavior which lead to a small violation of the conservation of the core charge bk_exact_xc ; engel_exact_xc ; trail_needs1 .
In this study, we present softcore HF pseudopotentials constructed using the RRKJ procedure which optimizes the potentials to yield rapid planewave basis cutoff convergence. The long tail behavior is fixed using a selfconsistent procedure following that of Trail and Needs trail_needs1 , which leads to a negligibly small violation of normconservation. These pseudopotentials are applied in HF calculations of several atomic properties yielding results in good agreement with the allelectron values. We also apply them in the study of the dissociation energies, equilibrium bond lengths, and frequency of vibrations of several dimers, using a HF planewave code cpmd . The allelectron and pseudopotential results are in agreement with each other, and the values are consistent with a similar comparison done using LDA pseudopotential and LDA allelectron calculations.
Generation of HF based pseudopotentials has been released in version 3.0 of the GPL package OPIUM opium and is available for download.
Vi Acknowledgments:
This work is supported by ONR grants N000140110365 and N000140510055, by the Department of Energy Office of Basic Energy Sciences grants DEFG0207ER15920 and DEFG0207ER46365, by the Air Force Office of Scientific Research, Air Force Materiel Command, USAF, Grant FA95500710397 and by the NSF grant EAR0530813. We are grateful to H. Krakauer and S. Zhang for many valuable discussions.
References
 (1) W. Kohn, Rev. Mod. Phys. 71, 1253 (1999).
 (2) For an excellent review see W. E. Pickett, Comput. Phys. Rep. 9, 115 (1989).
 (3) G. B. Bachelet, D. R. Hamann, and M. Schlüter, Phys. Rev. B 26, 4199 (1982).
 (4) A. M. Rappe, K. M. Rabe, E. Kaxiras, and J. D. Joannopoulos, Phys. Rev. B 41, R1227 (1990).
 (5) N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
 (6) For a review of ECPs, see for example, M. Dolg, in Modern Methods and Algorithms of Quantum Chemistry, Proceedings, Second Edition, ed. J. Grotendorst, John von Newmann Institute for computing, Jülich, NIC Series, Vol 3 (2000).
 (7) J. S. Binkley, J. A. Pople, W. J. Hehre, J. Am. Chem. Soc. 102, 939 (1980); W. J. Stevens, H. Basch, M. Krauss, J. Chem. Phys. 81, 6026 (1984); W. J. Stevens, M. Krauss, H. Basch, P. G. Jasien, Can. J. Chem. 70, 612 (1992); T. T. Cundari and W. J. Stevens, J. Chem. Phys. 98, 5555 (1993).
 (8) P. J. Hay and W. R. Wadt, J. Chem. Phys. 82, 284 (1985); 82, 270 (1985); 82, 299 (1985).
 (9) P. A. Christiansen, Y. S. Lee, and K. S. Pitzer, J. Chem. Phys. 71, 4445 (1979).
 (10) M. Dolg, U. Wedig, H. Stoll, and H. Preuss, J. Chem. Phys. 86, 1866 (1987); J. Sheu, S. Lee, and M. Dolg, J. Chin. Chem. Soc. (Taipei) 50, 583 (2003).
 (11) C. W. Greeff and W. A. Lester, Jr., J. Chem. Phys. 109, 1607 (1998).
 (12) I. Ovcharenko, A. AspuruGuzik, and W. A. Lester, Jr., J. Chem. Phys. 114, 7790 (2001).
 (13) J. R. Trail and R. J. Needs, J. Chem. Phys. 122, 014112 (2005).
 (14) J. R. Trail and R. J. Needs, J. Chem. Phys. 122, 174109 (2005).
 (15) A. D. Becke, J. Chem. Phys. 98, 5648 (1993).
 (16) D. M. Bylander and L. Kleinman, Phys. Rev. Lett. 74, 3660 (1995); Phys. Rev. B 52, 14566 (1995); 54, 7891 (1996); 55, 9432 (1997).
 (17) M. Städele, J. A. Majewski, P. Vogl, and A. Görling, Phys. Rev. Lett. 79, 2089 (1997).
 (18) E. Engel, A. Höck, R. N. Schmid, R. M. Dreizler, and N. Chetty, Phys. Rev. B 64, 125111 (2001).

(19)
http://www.tcm.phy.cam.ac.uk/ ~mdt26/casino2_pseudopotentials.html
 (20) C. F. Fischer, Comp. Phys. Comm. 98, 255 (1996).
 (21) G. S. Handler, D. W. Smith, and H. J. Silverstone, J. Chem. Phys. 73, 3936 (1980).
 (22) As in Ref. trail_needs1, , the deviation in the logarithmic derivative is computed neglecting changes in the Hartree and exchange potentials.
 (23) The logarithmic derivative is the ratio between the wavefunction and its slope, thus comparisons between logarithmic derivatives can be misleading if either value is close to zero. In all cases shown in Fig. 4, the magnitude of the wavefunctions and slopes ranged between and , respectively.
 (24) X. Gonze, R. Stumpf, and M. Scheffler, Phys. Rev. B 44, 8503 (1991).
 (25) OPIUM pseudopotential package http://opium.sourceforge.net.
 (26) E. Clementi and A. D. Mclean, Phys. Rev. 133, A419 (1963).; E. Clementi, A. D. Mclean, D. L. Raimondi, and M. Yoshimine, Phys. Rev. 133, A1274 (1963).
 (27) The ECPs are obtained from the Extensible Computational Chemistry Environment Basis Set Database (http://www.emsl.pnl.gov/forms/basisform.html).
 (28) CPMD version 3.9.2, Copyright IBM Corp 19902005, Copyright MPI für Festkörperforschung Stuttgart 19952001, http://www.cpmd.org.
 (29) Gaussian 98, Revision A.11.4, M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, V. G. Zakrzewski, J. A. Montgomery, Jr., R. E. Stratmann, J. C. Burant, S. Dapprich, J. M. Millam, A. D. Daniels, K. N. Kudin, M. C. Strain, O. Farkas, J. Tomasi, V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, S. Clifford, J. Ochterski, G. A. Petersson, P. Y. Ayala, Q. Cui, K. Morokuma, N. Rega, P. Salvador, J. J. Dannenberg, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. Cioslowski, J. V. Ortiz, A. G. Baboul, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. Gomperts, R. L. Martin, D. J. Fox, T. Keith, M. A. AlLaham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, J. L. Andres, C. Gonzalez, M. HeadGordon, E. S. Replogle, and J. A. Pople, Gaussian, Inc., Pittsburgh PA, 2002.
 (30) J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
 (31) X. Gonze, J. M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G. M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jollet, M. Torrent, A. Roy, M. Mikami, Ph. Ghosez, J. Y. Raty, and D. C. Allan, Comp. Mat. Sci. 25, 478 (2002).
 (32) T. H. Dunning, J. Chem. Phys. 90, 1007 (1989); T. H. Dunning, Jr., K. A. Peterson, and A. K. Wilson, J. Chem. Phys. 114, 9244 (2001).
 (33) S. G. Louie, S. Froyen, and M. L. Cohen Phys. Rev. B 26, 1738 (1982).
 (34) D. Porezag, M. R. Pederson, and A. Y. Liu , Phys. Rev. B 60, 14132 (1999).
 (35) E. L. Shirley, R. M. Martin, G. B. Bachelet, and D. M. Ceperley, Phys. Rev. B 42, 5057 (1990).
 (36) A. Höck and E. Engel, Phys. Rev. A 58, 3578 (1998).