{'fraction_non_alphanumeric': 0.037372751359643006, 'fraction_numerical': 0.009099149351554875, 'mean_word_length': 4.653330705557745, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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 1/r, 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.


I. INTRODUCTION
In nearly all plane-wave density functional calculations [1], 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 plane-wave basis cutoff [2].
In the physics community, pseudopotentials are typically constructed using density functional theory (DFT) due to its success in electronic structure calculations [3,4,5]. 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, Hartree-Fock (HF) pseudopotentials or effective core potentials (ECPs) are mostly used [6,7,8,9,10]. 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 [11,12,13,14]. 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 [11]. 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 [15]) 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 [13] have recently found that normconserving 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 [16,17,18], 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 long-range behavior, while retaining desired eigenvalue matching and norm conservation of the pseudopotential [16,17]. 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 [18]. Trail and Needs fixed the long-range unphysical behavior using a similar self-consistent approach [13] that also slightly violates norm conservation. Ovcharenko et al. [12] 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 [4]. We use a selfconsistent approach to restore the Coulombic tail in a manner similar to that used by Trail and Needs [13,19]. 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 [2]. 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. 20.
The pseudopotential construction starts by choosing an electronic reference state, and then solving the Schrödinger equation of the system for the eigenstates Φ nl (r) and eigenvalues ǫ nl . The wavefunction can be written as Φ nl (r) = φ nl (r)Y lm (θ, φ)/r where Y lm (θ, φ) are the spherical harmonics. The orbitals φ nl (r) satisfy the radial Schrödinger equation: (1) whereT = −d 2 /(2dr 2 )+l(l +1)/(2r 2 ),V ion (r) = −Z/r is the bare nuclear potential, and Z is the atomic number. V HF [{φ(r)}] is the HF potential defined such that, (4) For convenience, equal and opposite self interaction terms are included in both the Hartree and exchange terms. Therefore, the Hartree potentialV H [{φ(r)}] ≡V H [ρ(r)] is a functional of the spherically averaged total electronic density ρ(r) = nl |φ nl (r)| 2 .
The second step in the construction of norm-conserving pseudopotentials typically begins with the generation of a smooth set of pseudo-orbitalsφ nl (r) to replace the AE orbitals φ nl (r), such thatφ nl (r) equals the all-electron orbital beyond some cutoff distance r c . Different normconserving pseudopotentials mainly differ in the form of φ nl (r) for r < r c , and the guiding principles used to constructφ nl (r).
In the RRKJ method, the pseudo-orbitals in the core region are expressed as a linear combination of spherical Bessel functions j l (q k r) The Bessel wave-vectors q k are chosen such that, where j ′ l (q k r) and φ ′ nl (r) are the derivatives of j l (q k r) and φ nl (r) with respect to r, respectively. N b is the number of Bessel coefficients, c nlk .
The key ingredient of the RRKJ method is the optimization of the Bessel coefficients, c nlk , to minimize the residual kinetic energy defined as follows: whereΦ nl (q) is the Fourier transform of real space pseudo-orbitalΦ nl (r) =φ nl (r)Y lm (θ, φ)/r, and q c is the target wavevector above which the contribution to the residual kinetic energy is as small as possible. For a particular q c , the planewave energy cutoff required to achieve total energy convergence of ∆T nl is given by the square of the maximum value of q c over all orbitals [4]. The optimization of ∆T nl ({c nlk }, q c ) is constrained to pseudo-orbitals that are continuous at r c , 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 N b 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 N b = 6. The screened pseudopotentialV nl scr (r) is defined as the potential that makes the desired pseudo-orbitalφ nl (r) an eigenstate of the one-electron Schrödinger equation, with the same eigenvalue ǫ nl as the corresponding all-electron state: This equation can now be inverted, thanks to the nodeless character ofφ nl (r), to obtain the screened pseudopotential: Note thatV nl scr (r) is a continuous function because of the continuity requirements imposed on the pseudo-orbital φ nl (r) 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 V nl scr (r). That is, where {φ(r)} v 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, V nl ion (r), have different behavior for large r, which in turn stems from the corresponding large-r behavior of the HF and KS orbitals. It was shown by Handeler et al. [21] that the HF orbitals decay at infinity as where α = (−2ǫ HO ) 1/2 is determined by the eigenvalue of the highest occupied orbitals ǫ HO , and β nl is an orbital dependent quantity [21]. This behavior is to be contrasted with the exponential decaying behavior of the KS orbitals, In DFT-based pseudopotentials, the exchange potential is replaced by the exchange-correlation functional of the electron density ρ(r), which decays exponentially to zero for large r. Thus, the ionic potentialV nl ion (r) = −Z v /r for large r where Z v is the valence charge density.
The effect of the special long-range behavior of the HF orbitals on the ionic pseudopotentialV nl ion (r) can be understood by writingV nl ion (r) as, where the Schrödinger equation of Eq. (1) is used. Here ρ(r) is the total electronic charge density whileρ v (r) 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 r will be r β e −αr where β = max nl β nl , while the denominator will go as r β nl +1 e −αr . Therefore, the exchange potential, as computed in Eq. (13), is not guaranteed to decay faster than 1/r. In fact, the exchange potential can even grow without bound for large r when computed in this way. This is the source of the non-local long-range non-Coulombic tail in the HF pseudopotentials [13].
In this work we follow the Trail and Needs selfconsistent approach to localize the ionic pseudopotential The upper panel shows the Ne potential V nl ion constructed using the RRKJ method before the tail correction. The inset shows ∆Vion = V nl ion + Zv/r which displays the non-physical tail of the s-potential channel for large r. The lower panel shows V nl loc−ion after the tail correction procedure. In the inset, we show the large-r behavior of ∆V loc−ion = V nl loc−ion + Zv/r. [13]. The localized ionic potentialV nl loc−ion (r), which is forced to behave as −Z v /r for large r, is defined aŝ where ρ c (r) is the core charge density, 1/ √ ξ = r loc /16 is the characteristic distance for the localizedV nl loc−ion (r) to approach −Z v /r, and r loc is the localization radius. For all pseudopotentials, the value of r loc for each l-channel is fixed to be the same as the corresponding r c . γ nl (r) is defined as p nl + q nl f (r, r loc ) where p nl and q nl are parameters which are chosen to localize the ionic pseudopoten-tialV nl ion (r) in a self-consistent manner, and f (r, r loc ) 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 p nl and q nl 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 p nl and q nl parameters are then adjusted to minimize the error in the eigenvalues and logarithmic derivatives [22]. The procedure is then repeated until p nl and q nl parameters are found that give eigenvalues and logarithmic derivatives that closely match the all-electron values to the desired tolerance.
Finally, we define f (r, r loc ) used in γ nl (r). For the case of Troullier-Martins pseudopotentials, we chose the same form as was used in Ref. [13], namely f (r, r loc ) = r 4 1 − 2 r 2 /3r 2 loc for r < r loc , and f (r, r loc ) = r 4 loc /3 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 f (r, r loc ) = 1 − r/2r loc for r < r loc , and f (r, r loc ) = r loc /2 for r ≥ r loc . To illustrate the pseudopotential localization procedure, we show in the upper panel of Fig. 1 the pseu-dopotentialsV nl ion (r) 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 II: Ionization energy study of the first-and second-row elements. The first column shows the all-electron ionization energy, while the other columns show the deviation from the AE result. Columns 2 and 3 show the pseudopotential ionization energies as obtained using the original and the localized RRKJ pseudopotentials [Eqs. (13) and (14), respectively]. Columns 4 and 5 show the same values as obtained using the TM method. For comparison, we show also the values obtained by Trail and Needs [13]. The pseudopotential parameters are summarized in Table I  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 ofV nl ion (r) + Z v /r. Note, in particular, in the upper panel, the non-decaying tail of the s-potential, which can be shown [13] to behave as a + b/r for large r where a and b 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 r < r c 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 r c 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 ≈ 0.001 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 [23]. As can be seen from both figures, both pseudopotential construction methods yield similar results.

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 [25] 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 [25] 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 I. In order to aid in comparison to previous results, the reference configurations and the construction parameters are the same as of Ref. [13]. 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 N b and q c . We set N b = 6 for all pseudopotentials, and selected the q c 's such that ∆T nl ≈ 5 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 q c value used.
In Table II 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, ∆E ion , from the all-electron energy. ∆E ion is obtained using the original pseudopotential of Eq. (13) with the unphysical tail behavior, while ∆Ẽ ion 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. [13].
For both pseudopotential construction schemes, the agreement between ∆E ion and ∆Ẽ ion 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 ∆E ion and ∆Ẽ ion . 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 II are in good agreement with the all-electron values for both the RRKJ and TM methods. On average the difference is less than 0.5 milli-Hartree (mHa) and the largest deviation is less than 1.5 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 [13]. Any small differences could be attributed to the different grids used, or the slight differences in the self-consistent HF procedure.
In Table III, 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   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 IV. 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 3 P and 5 S 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 [28], and the all-electron results are obtained using a Gaussian basis code [29]. For comparison, we have also generated LDA [30] 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 all-electron and pseudopotential results we used the Perdew-Wang [30] flavor of LDA. We did not include a non-linear core correction (NLCC) [33,34] in We show pseudopotentials obtained using the RRKJ, TM, OAL [12], and the SBKJC [27] methods. For the RRKJ and TM methods, the rc values are the same as those in Table V. We show only the s and p potentials.
our LDA pseudopotentials for LDA/HF comparison purposes. The pseudopotential DFT calculations are carried out using ABINIT [31]. 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 [32]. 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 V, we show the cutoff radii we used to generate the RRKJ pseudopotentials for this study. These values, which are different from those of Table I, are chosen to accommodate the bond lengths of these dimers without core overlap when possible. In the case of P and Cl, the d 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 [12], and SBKJC [27] pseudopotentials in real space for the two elements. We did not include a d-channel in these plots because neither the SBKJC nitrogen pseudopotential nor the OAL potentials have a d-channel. For the RRKJ and TM pseudopotentials, we used the cutoff radii as shown in Table V. 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 [4] of the pseudowavefunctions, which we show in Fig. 6. The residual kinetic energy is calculated as shown in Eq. (7) by varying q c (we plot it against q 2 c 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 VI. 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 VII. 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 2 , and to a lesser extent in P 2 and even much less in Cl 2 . The HF pseudopotential dissociation energy of N 2 shows a significant deviation from the all-electron HF result, whereas for P 2 and Cl 2 , the results are the same to within 0.01 eV. A similar difference between all-electron and pseudopotential results was also seen in Ref. [18] 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. [34], these errors are larger for first row species compared to heavier atoms and increase with more unpaired electrons, explaining why N 2 has the largest error (first row, and 3 unpaired spins) and why Cl 2 (second row, 1 unpaired spin) has the smallest. Our LDA all-electron and pseudopotential energies are in excellent agreement with previous results [18] 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 pseudopo-tentials.
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 [35,36], 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 r. 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 [13,16,18].
In this study, we present soft-core HF pseudopotentials constructed using the RRKJ procedure which optimizes the potentials to yield rapid planewave basis cut-off convergence. The long tail behavior is fixed using a self-consistent procedure following that of Trail and Needs [13], 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 [28]. 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 [25] and is available for download.