Optimized norm-conserving Hartree-Fock pseudopotentials for plane-wave calculations
We report Hartree-Fock (HF) based pseudopotentials suitable for plane-wave calculations. Unlike typical effective core potentials, the present pseudopotentials are finite at the origin and exhibit rapid convergence in a plane-wave basis; the optimized pseudopotential method [A. M. Rappe et al. , Phys. Rev. B 41 1227–30 (1990)] improves plane-wave convergence. Norm-conserving HF pseudopotentials are found to develop long-range non-Coulombic behavior which does not decay faster than , and is non-local. This behavior, which stems from the nonlocality of the exchange potential, is remedied using a recently developed self-consistent 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 all-electron 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 all-electron results.
pacs:02.70.Ss, 71.15.-m, 31.25.-v
Present address: ]Physics Department, Clark Hall, Cornell University, Ithaca NY 14850
In nearly all plane-wave 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 electron-nuclear Coulomb interaction by a weaker (often angular-momentum dependent) potential. The advantage is the reduced computational cost, due to reduction in the number of electrons and in the plane-wave 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 plane-wave basis. On the other hand, in the chemistry community, Hartree-Fock (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 plane-wave calculations, mainly because of their singular behavior at the origin and slow convergence in a plane-wave basis.
Recently, there has been a growing interest in developing soft-core 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 DFT-based 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 exchange-correlation 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 norm-conserving HF pseudopotentials constructed using standard norm-conserving pseudopotential methods develop a non-decaying and non-local 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 long-range 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 long-range 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 self-consistent manner by altering the pseudo-orbitals 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 solid-state calculations engel_exact_xc . Trail and Needs fixed the long-range unphysical behavior using a similar self-consistent approach trail_needs1 that also slightly violates norm conservation. Ovcharenko et al. lester_psp have recently developed HF-based pseudopotentials which are finite at the origin and are expressed in terms of Gaussian basis, but are still relatively hard for typical plane-wave calculations, as we show in Sec. IV.
In this paper, we construct soft-core HF based pseudopotentials following the Rappe-Rabe-Kaxiras-Joannopoulos (RRKJ) method RRKJ . We use a self-consistent approach to restore the Coulombic tail in a manner similar to that used by Trail and Needs trail_needs1 ; casino_psp . Pseudopotentials are highly non-unique 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 plane-wave 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 plane-wave convergence of the atoms, and we investigate the transferability by looking at the ionization energy, electron affinity, and excitation energy of the first- and second-row elements. In Sec. IV, we study several dimers at the HF level, and compare both the all-electron 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 norm-conserving pseudopotentials starting from an all-electron (AE) atomic calculation pickett . However, in this case, we solve the Hartree-Fock instead of the Kohn-Sham 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:
where , is the bare nuclear potential, and is the atomic number. is the HF potential defined such that,
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 norm-conserving pseudopotentials typically begins with the generation of a smooth set of pseudo-orbitals to replace the AE orbitals , such that equals the all-electron orbital beyond some cutoff distance . Different norm-conserving pseudopotentials mainly differ in the form of for , and the guiding principles used to construct .
In the RRKJ method, the pseudo-orbitals in the core region are expressed as a linear combination of spherical Bessel functions
The Bessel wave-vectors are chosen such that,
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:
where is the Fourier transform of real space pseudo-orbital , 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 pseudo-orbitals 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:
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 non-linearity of Eq. (7). Throughout this study, we used .
The screened pseudopotential is defined as the potential that makes the desired pseudo-orbital an eigenstate of the one-electron Schrödinger equation, with the same eigenvalue as the corresponding all-electron state:
This equation can now be inverted, thanks to the nodeless character of , to obtain the screened pseudopotential:
Note that is a continuous function because of the continuity requirements imposed on the pseudo-orbital 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 exchange-potential contributions of the valence electrons from the screened potential . That is,
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 pseudo-orbitals are nodeless.
Up to this point, the pseudopotential construction recipe mirrors exactly the procedure done in regular norm-conserving 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
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 DFT-based pseudopotentials, the exchange potential is replaced by the exchange-correlation 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 long-range behavior of the HF orbitals on the ionic pseudopotential can be understood by writing as,
where the Schrödinger equation of Eq. (1) is used. Here is the total electronic charge density while is the pseudo-valence charge density. In the HF case, the exchange potential is a non-local 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 non-local long-range non-Coulombic tail in the HF pseudopotentials trail_needs1 .
In this work we follow the Trail and Needs self-consistent approach to localize the ionic pseudopotential trail_needs1 . The localized ionic potential , which is forced to behave as for large , is defined as
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 self-consistent manner, and is a polynomial function whose explicit form depends on the pseudopotential construction method (see below).
The self-consistent 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 all-electron values to the desired tolerance.
Finally, we define used in . For the case of Troullier-Martins 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||TM111Ref. trail_needs1 .|
|Average abs. error|
|Maximum abs. error|
|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 self-consistent 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 long-range behavior of . Note, in particular, in the upper panel, the non-decaying 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 all-electron and pseudo-orbitals after the localization procedure. The lower panel shows the difference in the pseudo-orbitals 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 core-norm will be different from the all-electron 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 non-zero value.
As mentioned before, the localization procedure conserves the all-electron 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 norm-conservation and the change in the logarithmic derivative at due to the localization procedure for the first- and second-row 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||TM111Ref. trail_needs1 .|
|Ground state||Excited state|
|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 second-row 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 Troullier-Martins 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 pseudo-orbital, 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 all-electron value (shown in the first column), and we only report the difference, , from the all-electron 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 self-consistent 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 non-Coulombic 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 all-electron values for both the RRKJ and TM methods. On average the difference is less than milli-Hartree (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 self-consistent HF procedure.
In Table 3, we show a similar study for the electron affinity of the first- and second-row 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 DFT-based 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 all-electron values. The absolute average deviations from the all-electron 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).
Iv Study of several dimers
In this section, we apply our HF pseudopotentials in a study of three dimers and compare them with all-electron calculations. The pseudopotential calculations are done using a HF planewave basis code cpmd , and the all-electron 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 all-electron counterparts on an equal footing with the HF results. In our all-electron and pseudopotential results we used the Perdew-Wang lda flavor of LDA. We did not include a non-linear 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 all-electron calculations are performed using the correlation-consistent cc-pV5Z basis set for N and cc-pV(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 6-311++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 second-row 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 all-electron 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 all-electron HF result, whereas for P and Cl, the results are the same to within eV. A similar difference between all-electron 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 all-electron and pseudopotential energies are in excellent agreement with previous results engel_exact_xc without a NLCC. This study also showed that the all-electron 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 core-valence interactions in the atomic reference state are frozen into the core and only valence-valence exchange terms are subtracted during the descreening step. As shown by our results and others shirley ; hock-engel , 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 core-valence contributions.
V Summary and conclusion
Constructing HF pseudopotentials suitable for planewave calculations is a non-trivial task considering the extreme nonlocality of the exchange potential. Pseudopotentials constructed using a typical norm-conserving procedure develop a non-local and a non-Coulombic 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 soft-core 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 self-consistent procedure following that of Trail and Needs trail_needs1 , which leads to a negligibly small violation of norm-conservation. These pseudopotentials are applied in HF calculations of several atomic properties yielding results in good agreement with the all-electron 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 all-electron and pseudopotential results are in agreement with each other, and the values are consistent with a similar comparison done using LDA pseudopotential and LDA all-electron calculations.
Generation of HF based pseudopotentials has been released in version 3.0 of the GPL package OPIUM opium and is available for download.
This work is supported by ONR grants N000140110365 and N000140510055, by the Department of Energy Office of Basic Energy Sciences grants DE-FG02-07ER15920 and DE-FG02-07ER46365, by the Air Force Office of Scientific Research, Air Force Materiel Command, USAF, Grant FA9550-07-1-0397 and by the NSF grant EAR-0530813. We are grateful to H. Krakauer and S. Zhang for many valuable discussions.
- (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. Aspuru-Guzik, 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).
- (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 1990-2005, Copyright MPI für Festkörperforschung Stuttgart 1995-2001, 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. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, J. L. Andres, C. Gonzalez, M. Head-Gordon, 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. Cara-cas, 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).