Normconserving pseudopotentials with chemical accuracy compared to allelectron calculations
Abstract
By adding a nonlinear core correction to the well established Dual Space Gaussian type pseudopotentials for the chemical elements up to the third period, we construct improved pseudopotentials for the Perdew Burke Ernzerhof (PBE) [J. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996)] functional and demonstrate that they exhibit excellent accuracy. Our benchmarks for the G21 test set show average atomization energy errors of only half a kcal/mol. The pseudopotentials also remain highly reliable for high pressure phases of crystalline solids. When supplemented by empirical dispersion corrections [S. Grimme, J. Comput. Chem. 27, 1787 (2006); S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010)] the average error in the interaction energy between molecules is also about half a kcal/mol. The accuracy that can be obtained by these pseudopotentials in combination with a systematic basis set is well superior to the accuracy that can be obtained by commonly used medium size Gaussian basis sets in allelectron calculations.
I Introduction
During the last decades, density functional theory (DFT) has proven its pivotal role for computational studies in the fields of condensed matter physics and quantum chemistry. Particularly the KohnSham (KS) formalism of DFT has gained enormous popularity as an ab initio method applicable to relatively large systems. An essential ingredient for many large scale implementations of KSDFT are pseudopotentials which are also frequently denoted as effective core potentials. By eliminating the strongly bound core electrons pseudopotentials reduce the number of occupied electronic orbitals that have to be treated in an electronic structure calculation. In addition the valence wavefunctions arising from a pseudopotential are much smoother than the allelectron valence wavefunction in the core region, since the orthogonality constraints to the rapidlyvarying wavefunctions carrying core electrons are missing. Since it is not necessary to describe rapidly varying wavefunctions the size of the basis set used for their representation can be reduced. These two factors lead to a significant reduction of the computational effort of a pseudopotential calculation compared to an allelectron calculation.
Even though it is well known that the valence electrons are responsible for the majority of chemical and physical properties of atoms, pseudopotentials have to be constructed very carefully in order to reproduce the properties of the allelectron atom accurately. If a pseudoatom, i.e an atom described by a pseudopotential, reproduces the allelectron behavior accurately for any chemical environment the pseudopotential is said to be transferable.
Pseudopotentials (PSP’s) are an essential ingredient of most electronic structure codes and different solutions are implemented in presentday DFT codes. Traditional normconserving (NC) approaches, e.g. [4] are formally the simplest approach, since they give rise to pseudowavefunctions which lead to a valid charge density. By introducing atomic like orbitals as additional basis functions any atomic Hamiltonian arising either from an allelectron potential or from a norm conserving pseudopotential [5] can be transformed into an Linearized Augmented Plane Wave (LAPW) like Hamiltonian [6]; [7]; [8]. The widespread ProjectorAugmented Wave (PAW) methods [9] and the ultrasoft pseudopotentials [10] are derived by such a transformation from an allelectron atom Hamiltonian. The number of required basis functions is reduced by this transformation, but the calculation of the charge density is more complicated and a generalized eigenvalue problem has to be solved even for the case of an orthogonal basis set. For applications in quantum chemistry, effective core potentials [11]; [12]; [13] are often optimized for a certain basis set and usually employed for heavier elements only.
In this paper, the Dual Space Gaussian pseudopotentials of GoedeckerTeterHutter (GTH) and HartwigsenGoedeckerHutter (HGH) [14] PSP are generalized by the inclusion of a NonLinear Core Correction (NLCC) term. These new pseudopotentials are able to provide an accuracy that is comparable to that of the very best allelectron calculations.
The starting point for understanding why pseudopotentials work is the subdivision of space into muffintin spheres centered on the atom in a molecule or solid and the remaining interstitial region [15]. A nonselfconsistent Schrödinger equation can be solved exactly in the interstitial region if one knows the scattering properties on the surface of the muffintin spheres [16]. The scattering properties are typically specified by the logarithmic derivative as a function of energy. This function is the quotient of the radial outward derivative and the functional value of the wavefunction on the surface of the muffintin sphere. In this way the boundary conditions are specified which are necessary to integrate the Schrödinger equation, a second order partial differential equation where the amplitude of the solution is fixed by a normalization constraint. A necessary but not sufficient condition for a pseudopotential to be transferable is therefore that the logarithmic derivatives of the allelectron and pseudoatom agree over the relevant energy interval. The construction of pseudopotentials is typically done using as the reference state a neutral isolated atom which has been spherically symmetrized. This symmetrization can be achieved by using identical and generally fractional occupation numbers for all the orbitals with the same and quantum numbers, e.g. for the set of , and orbitals. The norm conservation condition [17] ensures that the logarithmic derivative function describes well the scattering properties of a muffintin sphere containing the charge distribution of this reference configuration. In a selfconsistent calculation the charge distribution changes however when the free atom is inserted in a molecule or solid and the potential in the muffintin region will in general differ from the potential within a muffintin sphere of the same radius around the reference atom. Hence the scattering properties change and the pseudopotential constructed using the charge distribution in the muffintin sphere of the isolated atom might not well reproduce these modified scattering properties of a new chemical environment. Due to screening effects there exists however an invariant muffintin sphere within which the total electronic charge distribution is nearly independent of the chemical environment [18]. The radius of this invariant muffintin sphere is a fraction of the covalent bondlength and thus considerably smaller than the muffintin radii used in other methods such as the LAPW method. The scattering properties of this invariant muffintin sphere hardly vary as a function of the chemical environment of the atom. If the separable terms of a pseudopotential as well as the difference between the local part of the pseudopotentials and the pure coulombic potential remain localized within this invariant sphere the pseudopotential is expected to be highly transferable. This recipe was followed in the construction of the GTH [19] and HGH [14] pseudopotentials which are indeed well transferable for nonspinpolarized systems.
Despite the fact that the total charge in the invariant muffintin sphere is nearly identical in different chemical environments, the spin polarization is not, as illustrated in Fig. 1. Shown are the changes in the radial charge and spin densities if one adds half an electron to the unoccupied spin channel of the orbital. For spin polarized calculations the concept of an invariant muffintin sphere is therefore not applicable. One possibility to overcome this problem is to construct pseudopotentials which have an explicit spin dependence [20]. The other possibility is to include nonlinear core corrections (NLCC) [21] in the pseudopotential.
In the NLCC schemes the spin and charge densities in the muffintin sphere are not just the ones from the valence electrons treated explicitly by the pseudopotential but they are both the respective sums of the valence charge and the core charge given by the nonlinear core correction. Since the core electrons can be considered to be frozen, i.e to be invariant with respect to different chemical environments, this core charge is fixed once and for all. It is obvious that the spin polarization, i.e. the quantity is very poorly represented without a core charge. If for instance all valence orbitals are spin up then the spin down charge density would be zero and the spin polarization would be equal to one. In the real atom is not equal to one in the core region since the core electrons are never spin polarized. Since the density of the core electrons is much larger than the density of the valence electrons in the core region the spin polarization actually has typically small values. These correct small values of are reestablished by the core charge of a NLCC pseudopotential and exchange correlation functionals can provide reliable total energies. Nonlinear core corrections have therefore the potential to substantially improve the description of spin polarized states.
Whereas previous implementations of NLCC pseudopotentials [22] tried to faithfully represent the core charge, we follow here another approach. In the spirit of the GTH pseudopotentials where all terms have simple parametrized analytical forms we also represent the core charge density just as one single Gaussian. The amplitude and width of this Gaussian core charge distribution are then optimized by a fitting procedure in the same way as the other parameters of the pseudopotential.
Ii Methodology
The procedure for constructing the NLCC pseudopotentials is very similar to the one used for the construction of the GTH and HGH pseudopotential. In contrast to the original GTH and HGH pseudopotentials which were fitted to a single atomic configuration, the new NLCC pseudopotentials are fitted not only to the ground state but also to several excited and ionized electronic configurations where half an electron is added or removed possibly to or from different valence orbitals. The atoms are always considered to be spherically symmetric, but some of the configurations used for the fit have a spin polarization. The parameters of the dual space Gaussian ansatz [19] are fitted such that:

The occupied and first few unoccupied valence eigenvalues of the allelectron and pseudo atom match for all configurations used in the fit.

The charge inside the inert region of the pseudo atom matches the allelectron valence charge in the same region for all the orbitals for which the eigenvalues are matched for all configurations used in the fit. This means that the pseudopotential is norm conserving for all the configurations used in the fit.

A high precision of a.u. is achieved for valence eigenvalues and charge integrals of the nonpolarized ground state, whereas only a moderate precision of a.u. is enforced for all other orbitals and configurations considered.

The total energy differences of the allelectron atom are reproduced for all configurations used in the fit.

The spin polarization energies of the allelectron atom are reproduced for all spin polarized configurations used in the fit.
Since the considered quantities are fitted for several configurations atomic transferability is already built in to these new pseudopotentials by construction.
The core charge is represented by a single Gaussian which is optimal for numerical efficiency. It is initialized such that it approximates well the physical core charge density and it is held fixed during an initial stage of the fit. Then both the amplitude and the width of the Gaussian are released and considered as fitting parameters. As a consequence the total amount of core charge and the width can differ from the physical value.
The parameters of the core charge constitute thus a small set of only two additional degrees of freedom. Yet this allows to optimally reproduce atomic polarization energies without degrading the transferability and accuracy of other atomic properties. It was found that the inclusion of a more complicated core charge is not beneficial. Furthermore, it should be emphasized that the novel NLCC pseudopotentials are not harder than their HGH counterparts. The smoothness of the core charge seems to play an important role for the fact that the hardness is not affected, and roughly the same grid spacings or energy cutoffs can be employed as for conventional HGH pseudopotentials.
In particular, pseudopotentials with NLCC were generated for boron, carbon, nitrogen, oxygen, fluorine, aluminium, silicon, phosphorous, sulfur and chlorine. Very weak spin dependences are expected for the rare gasses, and for all remaining chemical elements up to the third row, NLCC are found to be unnecessary, as HGH pseudopotentials are available that either include semicores (sodium and magnesium) or leave no core states at all (hydrogen, lithium and beryllium). For the special case of hydrogen, it was found that the multi configuration fit gave slightly improved results even though obviously no core charge was added. Since the focus of this paper is on systems made out of light elements, no relativistic effects such as spinorbit coupling were included in the pseudopotentials.
Iii Results
To assess the accuracy of the new pseudopotentials extensive calculations were performed for different test sets. The accuracy of covalent bond formation energies was examined for the standard G21 test set [23]; [24]; [25]; [26]. For the assessment of the accuracy of nonbonded interactions the S22 [27] test set was used. To check the performance for materials under high pressure we chose carbon, silicon, silicon carbide and boron nitride as test systems.
All pseudopotential calculations were done with the BigDFT package [28]. The BigDFT code uses a systematic wavelet basis set which allows to obtain the exact density functional solution with arbitrarily small error bounds. The parameter were set such that an accuracy of at least Hartree was obtained. The LibXC library [29] is used within BigDFT for the evaluation of the exchange correlation functional. Semiempirical vander Waals corrections were added in BigDFT according to the DFTD2[2] and DFTD3[3] methods for the calculations of the S22 test set.
To obtain reliable allelectron reference values for the atomization energies of the G21 test set, we performed allelectron calculations with the NWChem software package [30] using one of the largest available Gaussian type basis sets, namely an augmented correlation consistent polarized valence quintuple zeta Gaussian type basis set (augccpV5Z). Care was taken to disable symmetry detection and to check for the lowest energy spin multiplicity. For the chemical elements Li, Be Na and Mg, the augccpV5Z set was not available, so the corresponding quadruple zeta set (augccpVQZ) was used to compute the atomization energies of , , , , , and . To obtain the atomization energies of the relaxed molecules, geometry optimizations were carried out using the very same basis set.
Atomic allelectron calculations of the spin polarization energies were done with our nonspherical atomic code, which expresses the wavefunctions as a product of spherical harmonics and radial functions. The radial function are given numerically on a logarithmic grid. The settings were chosen such that a precision of at least Hartree can be obtained for the total energy. This required angular integration grids of 232 points and multipole representations up to . The atomic LSDA reference energies from the National centre of science and technology (NIST)[31] where reproduced within the given precision for all elements considered.
We calculated the atomization energies also with the three different sets of PAW [9] potentials available in VASP [32]. Those PAW potentials are derived from the allelectron atomic Hamiltonian and aim at allelectron accuracy. In order to obtain the required high precision some parameters had to be set to tighter values than the default values. We had to use for the general accuracy (prec = High Accurate and lasph = true) to activate nonspherical gradient corrections inside the PAW spheres. It was carefully checked that the calculations were converged with respect to the size of the periodic simulation cell. Furthermore, care was taken that the correct spin multiplicity and nonfractional occupations were produced. Hard PAW potentials were available for all required elements except for Li, Be, Na and Mg, for which semicore potentials were used instead. For comparison, all energies were recomputed with a set of default potentials. The third set consists of soft potentials for the elements B, C, N, O and F and default potentials otherwise. For the periodic solids, allelectron calculations have been performed using the fullpotential linearized augmented plane wave (FLAPW) and augmented plane wave plus local orbitals (APW + lo) methods as implemented in the WIEN2k[33] software package. We used a reduced muffintin radii for all atomic sorts in order to avoid their overlap up to the highest studied pressures. The sphere radii were kept fixed throughout the whole set of lattice parameters to obtain the best possible error cancellation. Semicore states were treated as valence, because high compressions can lead to an overlap of their wavefunctions, which will give a contribution to the energy. Inside the spheres, the partial waves were expanded up to lmax = 10. The number of plane waves was limited by a cutoff parameter RMTKmax = 9.0 for all the compounds under consideration. The charge density was Fourier expanded with Gmax = 14 a.u. For the majority of the systems we used a very dense kpoints grid (151515) to ensure total energy convergence.
All the calculations were done at zero electronic temperature, i.e. no Fermi smearing was used. Zero point energies were not included in any of our results.
iii.1 Atomization energies of the G21 test
Atomization energies are frequently used to assess the quality of various exchange correlation functionals as well as other approximations used in electronic structure calculations. The Gaussian G21 test set [23]; [24]; [25]; [26] is a standard benchmark set of 55 molecules in this context. Since this test set does not contain molecules with the chemical elements B, Al and Mg we added the molecules BH, BH2, AlH, AlH2, Mg2 und MgH. We used this augmented test set to compare our pseudopotential results with allelectron calculations. Because of Hund’s rule most isolated atoms are strongly spin polarized. When an atom is inserted into a molecule or solid, its spin polarization is typically strongly reduced. Since standard pseudopotentials are based on a nonspin polarized reference configuration they can typically better describe atoms in molecules or solids than isolated atoms themselves. Since the atomization energy is the difference between the total energy of the molecule and the sum of the total energies of its constituent isolated atoms, the largest contribution to the error in the atomization energy of a pseudopotential calculation comes actually from the atomic energies.
The atomization energies of the molecules in the G21 test set were first computed using conventional HGH pseudopotentials [34] for the PBE exchange correlation functional. A comparison with allelectron data is shown in figure 2. The spin multiplicity of systems with a net magnetic moment are indicated in brackets and omitted for closed shell systems. Deviations of kcal/mol are indicated with a (green) shading to relate the errors to the requirements for chemical accuracy.
It is found that the direct computation of the electronic atomization energies with the conventional pseudopotentials leads to significant disagreement with the results obtained in an allelectron calculation. An rather high mean absolute deviation (MAD) of 6.83 kcal/mol to the electron reference values for all 55 molecules in the G21 set is found. However, the main contribution to the error in the atomization energies comes from the estimation of the energy of the isolated atoms. Therefore, the atomization energies can be improved significantly by a two step procedure where the atomization energies are calculated as a sum of two terms. The first term is the energy difference between the molecule and the sum of the total energies of isolated, spherical and nonspinpolarized atoms. It thus can be considered as the atomization energy with respect to a set of nonphysical atoms. This energy difference is calculated with the HGH pseudopotentials and is fairly accurate since no strong spin polarizations are involved. The second term is the difference in total energy between the real, i.e nonspherical and spin polarized, atom and the previously defined nonphysical atom. This second term is calculated with our allelectron program for nonspherically symmetric atoms and is therefore exact. Since the atomic spin polarization energies and energy terms for breaking the spherical symmetry are only a property of the atoms they can be considered as a set of atomic correction terms for the accurate calculation of atomization energies. The atomic correction terms for the chemical elements considered in this study are listed in Table 1.
H  Li  B  C  N  O  F  Na  Al  Si  P  S  Cl 
25.4  6.8  10.4  31.7  72.0  43.7  16.0  5.1  7.1  19.7  43.1  23.8  7.6 
It has to be stressed that these atomic spin polarization energies drop out in most instances such as in the calculation of energy differences in a chemical reaction where only molecules are involved. Using this two step scheme, the errors in the atomization energies are decreased considerably to a MAD of 1.56 kcal/mol. Because of the cancellation effect, this is the accuracy that can be expected in the majority of energy differences calculated with the standard HGH type pseudopotentials. The above MAD value was obtained with the bond lengths and angles fixed as given on the computational chemistry comparison and benchmark database[35]. Using instead the equilibrium geometries of each method, i.e. the HGH pseudopotential and the allelectron calculation, the MAD value is slightly decreased to 1.52 kcal/mol. This last value is actually more relevant in practice since atomization energies for an unknown system necessarily have be calculated with theoretically determined geometries.
The new Gaussian type pseudopotential with a NLCC can however still considerably improve the accuracy without the need of using a two step procedure. The error of a direct computation of the atomization energies decreases to a MAD of only 0.52 kcal/mol. Using the equilibrium geometries obtained with the new NLCC pseudopotential the error drops again slightly to 0.51 kcal/mol. More important than this small improvement of the MAD is the fact that the result could be improved for the few molecules where the error was well above the average.
In figure 3 the accuracy of the HGH pseudopotentials with NLCC is shown for relaxed molecular geometries and compared with the results of PAW calculations published by Paier et. al.[36]. In this work hard PAW potentials were used and we were able to reproduce their results. In essence, the absolute errors of the new NLCC pseudopotentials are comparable with those using hard PAW potentials. As shown in Figure 4 the accuracy however goes down significantly when one uses the default or even the soft PAW potentials of the VASP package [32]. Furthermore, the same figures shows the discrepancies between allelectron results obtained in two large Gaussian basis sets while keeping the molecular geometries fixed. Even at this size the differences between the two basis sets are not negligible compared to the deviations to other methods and the accuracy of the pseudopotential method is indeed close to the discrepancies between different choices of allelectron reference values. This is quite surprising given the fact that these simple chemical compounds show only straightforward covalent type bonding properties which are certainly easier to describe with a Gaussian basis set than other more complex bonding patterns. It has also to be stressed that the computational cost rises very steeply when one goes from a medium size basis set to these very large basis sets. This is in contrast to the wavelet method where a modest decrease of about 15 percent in the grid spacing results in an gain of a factor of ten in accuracy because of the high order convergence rate of .
The accuracy problems of Gaussian basis sets become even more evident if one employs medium size or small standard basis sets in an allelectron calculation. The 631G, 631++G*, 631+G** and 6311++G(3df,3pd) basis sets were employed to compare the relative accuracy of the pseudopotential method with the incompleteness of and disagreement between standard Gaussian basis of various sizes. Figure 4 clearly shows that the accuracy obtained with these basis set is considerably lower than the accuracy with the NLCC pseudopotentials or also with the standard HGH pseudopotential within the two step procedure described above.
A summary of the deviations in the atomization energies averaged over the molecules of the G21 test set is given for fixed and relaxed geometries in tables 2 and 3, respectively. Indicated are the MAD, RMSD, mean signed deviation (MSD), maximum absolute deviation (maxAD) and minimum absolute deviation (minAD).
The last row of table 3 describes the change in the allelectron reference values when going from the fixed, experimental (CCCBDB) to relaxed geometries in the augccpV5Z basis set. This gain in energy upon geometry relaxation is significant compared to the assessed accuracy of the pseudopotential based methods, which are found to be very reliable for geometry optimizations.
MAD  RMSD  MSD  maxAD  minAD  
NLCCHGH  0.52  0.65  0.15  1.52  0.01 
HGH Krack  6.82  9.13  6.74  23.98  0.05 
Twostep  1.56  2.09  0.91  5.86  0.04 
PAW medium  1.20  1.89  1.14  7.15  0.01 
PAW soft  4.84  8.53  3.77  30.23  0.05 
6311++G(3df,3pd)  0.43  0.53  0.36  1.48  0.02 
631++G*  6.53  7.76  6.53  27.78  0.27 
631G  22.13  31.16  22.13  151.88  0.37 
MAD  RMSD  MSD  maxAD  minAD  
NLCC  0.51  0.63  0.16  1.50  0.03 
HGH Krack  6.85  9.13  6.76  23.94  0.10 
Twostep  1.52  2.05  0.88  5.73  0.01 
PAW Paier  0.46  0.56  0.43  1.13  0.01 
allelectron geopt  0.29  0.70  0.29  4.21  0.00 
iii.2 Accuracy of the equilibrium geometries
In order to compare the accuracy of the equilibrium geometries of the pseudopotential and allelectron calculations, the optimized geometry of each molecule is aligned with its allelectron counterpart, such that the RMSD is minimized [37]. The resulting RMSD values are shown in figure 5. It is observed that conventional HGH pseudopotentials yield already very good agreement with the allelectron data. Nevertheless, the inclusion of NLCC leads to a systematic improvement of the equilibrium geometries.
It can be noted that in a previously mentioned work [36], a similar test was carried out for the bond lengths of some dimers in order to verify the accuracy of the PAW method. It is found that our pseudopotential approach yields geometry data of at least the same or even better accuracy, and that the high precision is maintained when moving to more complicated geometries.
iii.3 Evaluation of pressure of extended systems
Next we benchmark pseudopotential (PSP) calculations for extended systems. A few crystalline systems made of light elements (diamond carbon, Silicon Carbide, Bulk Silicon and Boron Nitride) were selected and the pressure at a given lattice parameter was then compared between different approaches.
The details on how the stress energy tensor can be calculated in GGA for NLCC terms are given in the appendix. Figure 6 shows the difference between the LAPW results and the PSP results at the same lattice parameter. In addition we also show the results for the hard PAW potentials. In order to show the relative accuracy of the pressure we specify the lattice constant along the x axis of the figure in terms of the pressure. At the highest pressures the lattice constant is reduced by about 10 percent compared to the lattice constant with zero pressure. With NLCC PSP, an excellent relative accuracy of about is found. In this case, it can be noticed that the inclusion of a NLCC term improves further the results even though the systems under pressure are not spinpolarized. Results of similar quality can be obtained within the hard PAW scheme described above.
iii.4 Dispersioncorrected functionals
Long range van der Waals interactions are missing in all standard LDA and GGA density functionals. Adding semiempirical classical van der Waals interactions has however recently been demonstrated to give a rather accurate description of weakly bonded systems and is now frequently used. We will therefore examine the accuracy of the semiempirical models in the context of our pseudopotential calculations with a systematic wavelet basis set.
In BigDFT, we implemented two semiempirical models to correct dispersion energies and energy gradients DFTD2[2] and DFTD3[3]. The parameters of these models were separately fitted for each exchange correlation functional based on thermochemical data for weakly interacting systems. Since BigDFT uses a wavelet basis and pseudopotentials Figure 7 and Table 4 show the comparison of interaction energies of the benchmark database S22[27], with a reference calculation using Coupled Cluster CCSD(T) in the complete basis limit (CBS)[38]. The inclusion of dispersion correction D2 into BigDFT clearly improves the description of weak interactions within PBE, even though the S22 data set was not used as the fitting data set. The rootmeansquaredeviation (RMSD) between the CCSD(T) reference values and the NLCCDFT interaction energies is 0.58 kcal/mol. The absolute maximum difference corresponds to acetic acid dimer (COOH), where the overestimation is 1.57 kcal/mol, that means an 8% of the total interaction energy. The largest relative error of 35 % is found for the methane dimer whose interaction energy is only 0.2 kcal/mol. The errors for these systems are comparable to those that are obtained when PBED2 and PBED3 are used with a large basis set (augccpVDZ and augccpVTZ). On average the PBED2 scheme performs better with BigDFT than with any Gaussian basis set, while the PBED3/BigDFT results are comparable to the results obtained with PBED3/augccpVTZ electron calculations.
MAD  RMSD  MSD  maxAD  minAD  
NLCCno corr.  2.59  3.61  2.59  10.08  0.05 
HGH Krackno corr.  2.64  3.66  2.64  10.18  0.01 
NLCCD2  0.51  0.64  0.39  1.57  0.05 
HGH KrackD2  0.50  0.58  0.34  1.25  0.13 
NLCCD3  0.48  0.64  0.03  1.14  0.05 
HGH KrackD3  0.47  0.64  0.02  1.95  0.03 
augccpVDZD2  1.05  1.15  1.05  2.33  0.49 
augccpVTZD2  0.53  0.68  0.51  1.60  0.02 
augccpVTZD3  0.44  0.57  0.14  1.43  0.07 
Iv Discussion and conclusions
We have shown that our new NLCC PSP’s give very high accuracy for a wide range of applications. In particular they give atomization energies with chemical accuracy compared to allelectron calculations for the G21 test set. This accuracy can easily be obtained with a systematic basis set such as wavelets where one has to change only a single parameter to obtain arbitrarily high accuracy. Obtaining such a high accuracy with Gaussian basis sets requires using the largest available basis sets and is therefore frequently not feasible in practice. Contrary to a widespread belief, PAW calculations do not necessarily give allelectron accuracy. Soft PAW potentials can actually lead to appreciable errors. Well constructed hard PAW potentials on the other hand give very high accuracy and are together with our new normconserving pseudopotentials in practice the only feasible way to highest quality results for large systems.
V Appendix: NLCC HGH pseudopotentials in Kohnsham DFT formalism
The PSP format is based on HGHKrack form [34]:
(1) 
where the first part is a local potential
(2) 
and the second part the nonlocal term which is separated into different channels , each one defined in terms of separable projectors
(3)  
(4) 
The core charge of the new PSP’s is given by
(5) 
H  1  1  
0.20000  4.07312  0.68070  
B  5  3  
0.43250  4.26853  0.59951  
0.37147  6.30164  
0.33352  0.43364  
C  6  4  
0.31479  6.92377  0.96360  
0.30228  9.57595  
0.36878  0.00996  
0.27440  0.76008  
N  7  5  
0.24180  10.04328  1.39719  
0.25697  12.96802  
0.15686  0.73453  
0.24612  0.66037  
O  8  6  
0.26100  14.15181  1.97830  
0.22308  18.37181  
0.26844  0.10004  
0.25234  0.44314  
F  9  7  
0.20610  19.86716  2.79309  
0.19518  23.47047  
0.17154  0.61254  
Al  13  3  
0.35000  1.20404  2.14849  
0.46846  2.69262  0.00000  
2.15425  
0.54697  2.13804  
0.48775  0.38780  
Si  14  4  
0.33000  0.07846  0.79378  
0.42179  2.87392  0.02559  
0.00000  2.59458  
0.48800  2.47963  
0.44279  0.41540  
P  15  5  
0.34000  1.62258  0.72412  
0.38209  3.47754  0.01267  
3.47461  
0.43411  3.37859  
0.39868  0.45667  
S  16  6  
0.33000  1.49043  0.73314  
0.37046  6.18605  0.00000  
2.57761  
0.39772  3.89113  
0.38622  0.57500  
Cl  17  7  
0.32000  0.27448  
0.32659  4.20336  0.00000  
4.55652  
0.36757  4.22908  
0.42148  0.29324 
The core density is then used in the KohnSham total energy expression as follows:
(6) 
where and are the XC energy and potential respectively, is the Hartree potential and the ’s are KS wavefunctions, whose summed squares give the valence density . Eq.6 ensures HellmannFeynman condition at selfconsistency, .
The contribution to the stress tensor coming from the XC term with NLCC can be shown to be given by
(7) 
where is the supercell volume and . The formula shows that the gradient of is needed to evaluate , even for a LDA computation. A detailed derivation of DFT of stress (without NLCC) was shown by Dal Corso and Cresta [39].
Vi Acknowledgements
We acknowledge support from the Swiss National Science Foundation (SNF). Computer calculations were also performed at the Swiss National Supercomputing Center (CSCS) in Lugano. Partially, this research used resources of the Argonne Leadership Computing Facility at Argonne National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under contract DEAC0206CH11357.
References
 J. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
 S. Grimme, Journal of Computational Chemistry 27, 1787 (2006).
 S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
 D.R. Hamann, M. Schlüter, C. Chiang, Phys. Rev. Lett. 43, 1494 (1979),
 G. Lippert, J.Hutter and M. Parrinello, Theor. Chem. Acc. 103, 124 (1999).
 David J. Singh und Lars Nordstrom, Planewaves, Pseudopotentials, and the LAPW Method, Springer, New York, 2006.
 S. Goedecker and K. Maschke, Phys. Rev. B 45, 1597 (1992).
 S. Goedecker and K. Maschke, Phys. Rev. B 42, 8858 (1990).
 P. Blöchl, Phys. Rev. B 50 (1994) 17953, G. Kresse, D. Joubert, Phys. Rev. B 59 (1999) 1758, M. Torrent et al., Comput. Mat. Sci. 42 (2008) 337351.
 D. Vanderbilt, Phys. rev. B 41, 7892 (1990).
 A. Bergner, M. Dolg, W. Kuechle, H. Stoll, H. Preuss, Mol. Phys. 80, 1431 (1993).
 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)
 P. J. Hay and W. R. Wadt, J. Chem. Phys. 82, 284 (1985).
 C. Hartwigsen, S. Goedecker, and J. Hutter, Phys. Rev. B 58, 3641 (1998).
 Richard P Martin (2004). Electronic Structure: Basic Theory and Applications. Cambridge University Press. ISBN 0521782856.
 B. J. Austin, V. Heine, and L. J. Sham Phys. Rev. 127, 276 (1962).
 G. B. Bachelet, D. R. Hamann, and M. Schlüter, Phys. Rev B 26 (4199) (1982).
 S. Goedecker and K. Maschke, Phys. Rev. A 45, 88 (1992).
 S. Goedecker, M. Teter, and J. Hutter, Phys. Rev. B 54, 1703 (1996).
 S. C. Watson and E. A. Carter Phys. Rev. B 58, R13309 (1998).
 S. G. Louie, S. Froyen, M. L. Cohen, Phys. Rev. B 26, 1738 (1982).
 D. Porezag, M. R. Pederson, A. Y. Liu, Phys. Rev. B 60, 14132 (1999)
 L. A. Curtiss, K. Raghavachari, P. C. Redfern, and J. A. Pople, J. Chem. Phys. 106, 1063 (1997).
 J. A. Pople, M. HeadGordon, D. J. Fox, K. Raghavachari, and L. A. Curtiss, J. Chem. Phys. 90, 5622 (1989).
 L. A. Curtiss, C. Jones, G. W. Trucks, K. Raghavachari, and J. A. Pople, J. Chem. Phys. 93, 2537 (1989).
 L. A. Curtiss, P. C. Redfern, K. Raghavachari, and J. A. Pople, J. Chem. Phys. 109, 42 (1998).
 P. Jurecka, J. Sponer, J. Cerný, and P. Hobza, PCCP 8, 1985 (2006).
 L. Genovese et al., J. Chem. Phys. 129, 014109 (2008).
 Miguel A. L. Marques, Micael J. T. Oliveira, and Tobias Burnus, Comput. Phys. Commun. 183, 2272 (2012).
 M. Valiev, E.J. Bylaska, N. Govind, K. Kowalski, T.P. Straatsma, H.J.J. van Dam, D. Wang, J. Nieplocha, E. Apra, T.L. Windus, W.A. de Jong, Comput. Phys. Commun. 181, 1477 (2010).
 See http://physics.nist.gov/PhysRefData/DFTdata/Tables/ptable.html for atomic total energies and eigenvalues.
 G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
 P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka, and J. Luitz, WIEN2k, An Augmented Plane WaveLocal Orbitals Program for Calculating Crystal Properties (Karlheinz Schwarz, TU Vienna, Austria, 2001).
 M. Krack, Theor. Chem. Acc. 114, 145 (2005).
 NIST Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 101, Release 15b, August 2011, Editor: Russell D. Johnson III. http://cccbdb.nist.gov/
 J. Paier, R. Hirschl, M. Marsman and G. Kresse, J. Chem. Phys. 122, 234102 (2005).
 W. Kabsch, Acta Crystallographica 32, 922 (1976).
 T. Takatani, E.G. Hohenstein, M. Malagoli, M.S. Marshall, and C.D. Sherrill, J. Chem. Phys. 132, 144104 (2010).
 A. Dal Corso and R. Resta, Phys. Rev. B 50, 43274331 (1994)