Precision and efficiency in solid-state pseudopotential calculations

Despite the enormous success and popularity of density-functional theory, systematic verification and validation studies are still limited in number and scope. Here, we propose a protocol to test publicly available pseudopotential libraries, based on several independent criteria including verification against all-electron equations of state and plane-wave convergence tests for phonon frequencies, band structure, cohesive energy and pressure. Adopting these criteria we obtain curated pseudopotential libraries (named SSSP or standard solid-state pseudopotential libraries), that we target for high-throughput materials screening ("SSSP efficiency") and high-precision materials modelling ("SSSP precision"). This latter scores highest among open-source pseudopotential libraries available in the $\Delta$-factor test of equations of states of elemental solids.


INTRODUCTION
In the last three decades, atomistic electronic-structure methods have radically changed the way we think at materials theory and simulations. The 1998 Nobel prize in Chemistry given for density-functional theory (DFT) acknowledges this shift of paradigm. The ability to reduce the complexity of the many-body quantum-mechanical problem involving interacting electrons and nuclei into more tractable forms and algorithms allowed to leverage inexpensive and exponentially growing computational power, in order to provide sufficiently precise predictions for a great number of materials. What followed has been a flourishing of methods to compute more and more complex materials properties, and most notably spectroscopies (such as Raman, IR, ARPES, EELS, NMR and more). Nowadays, the (over)confidence in first-principles methods is such that they are routinely used to help interpreting experiments and guide the discovery and understanding of novel materials. In particular, systematic DFT-based computational materials screening is a fast-growing field of research, as reflected by the creation, in the last few years, of several research centres worldwide on computational materials discovery. Up to now, the most visible output has been the creation of large databases of materials properties obtained from first principles, to be compared with or to augment experimental databases such as the Pauling file [1], COD [2] or ICSD [3]. Even from a purely financial perspective, the personnel cost for plane-wave studies is of the order of 1 billion US$ per year [4], in purchasing power parity (PPP) terms, matched by substantial hardware usage. It is ever so more remarkable that in spite all of this, the efforts of verification of the precision of the underlying pseudopotentials (PSPs) or projector augmentedwave (PAW) approximations have been minimal. Only in 2016 a multi-group effort was able to establish a baseline in the calculations of the equations of state of elemental crystals [6]. In addition to the issue of precision, or of ver-ification (i.e. insuring that the removal of the core electrons from the calculations performs with the required tolerance), the issue of performance looms large -a softer, smoother PSP will allow straightforwardly faster calculations, both because the basis set is decreased as the computational cost scales with the square of the basis size, and because the minimization or iterative approaches can become more efficient or better pre-conditioned with a smaller basis.
Here, we introduce a PSP testing protocol based on extensive DFT and density-functional perturbation theory (DFPT) calculations of elemental solids, and select the optimal PSP for 85 elements of the periodic table. Our protocol, named SSSP (standard solid-state pseudopotential) testing protocol, is made of a verification part, based on the ∆−factor (i.e. the difference between allelectron and PSP equations of state) [6,7], and an extensive performance-oriented part based on plane-wave convergence tests for phonon frequencies, band structures, cohesive energies and stress tensors. We underline here that the SSSP testing protocol is a protocol based on verification and not on validation (following the nomenclature of Ref. [6]). Indeed our primary goal is to test the precision (verification) of the PSPs and thus to furnish, through our testing protocol, PSP libraries that give results as close as possible to the "exact" theoretical results for the PBE functional, as they would be obtained through a perfectly converged all-electron calculation. We do not perform any comparison with respect to experimental results, i.e. we do not test the accuracy (validation) of the PSPs. We consider up to eight (depending on the element) publicly available PSP libraries for the PBE functional [8] and test them with the PWscf and Phonon codes of the Quantum ESPRESSO (QE) distribution [9] in an automated fashion within the framework of the AiiDA [10] infrastructure for reproducible computational science. Ai-iDA also allows straightforward dissemination of results through the Materials Cloud web platform [11], a cloud service designed to enable data sharing in computational arXiv:1806.05609v2 [cond-mat.mtrl-sci] 7 Dec 2018 materials science.
In this paper, first we describe the test set of physically relevant quantities defined for our SSSP testing protocol and what are the corresponding selection criteria used in order to select the best PSPs. Then, we discuss the need for testing several properties, showing how different PSPs may yield similar results and convergence behaviours for one property but different ones for another property. Finally, we propose two optimal PSP libraries chosen according to the SSSP testing protocol and criteria.

SSSP testing protocol
In this work we investigate the precision and performance of several PSPs libraries available for the QE distribution [9]. QE is an integrated suite of open-source codes for electronic-structure calculations based on DFT which uses plane-waves as basis set and PSPs to represent the electrostatic electron-ion interactions. Nowadays QE is one of the most popular DFT codes adopted by researchers for the first-principles study of materials properties, with more than 2000 citations last year, according to Google Scholar™. All the tested PSP libraries are based on the generalized gradient approximation (GGA) for the exchangecorrelation functional of Perdew, Burke and Ernzerhof (PBE) [8] and they include the three main pseudization approaches [12]: norm-conserving (NC) [13], ultrasoft (US) [14] and projector-augmented wave (PAW) [15]. In particular, we investigate three PAW PSP libraries (pslibrary.0.3.1 [16], pslibrary.1.0.0 high accuracy [17] and the library proposed by Topsakal and Wentzcovitch for the rare-earth elements [18]), five US libraries (GBRV (versions 1.2, 1.4 and 1.5) [19], pslibrary.0.3.1 [16], and pslibrary.1.0.0 high accuracy [17]) and two NC libraries (SG15, versions 1.0 and 1.1 [20]). For a few selected elements, i.e. N, O, F and Hf, in which all the PSP libraries above perform less well in the ∆-factor test we also consider the recent NC Pseudo Dojo library [21] and, only for N, our own set (called THEOS) of US PSPs. Besides, for the elements from H to Ne in the periodic table, we test the NC library proposed by Willand et al. [22] which is tailored for systems made of light elements (see Table 1 for a list of all the PSP libraries tested).

Equation of state.
In order to assess the precision of PSPs, we compute the ∆-factor, i.e. the integral of the difference between the equations of state calculated with PSP simulations and with reference all-electron results. For this purpose we use the protocol introduced in 2014 by Lejaeghere et al. [7]. This protocol was recently exploited to compare 15 different DFT codes, including both all-electron and PSP codes, in order to verify the reproducibility of the PBE equations of state of elemental crystals across different methods and implementations [6]. The protocol consists in calculating the energy-versusvolume at 7 equidistant points centred around the reference equilibrium volume and then performing a Birch-Murnaghan fit. From the parameters of the fit some important physical quantities related to the structural and elastic properties of the system are extracted: the equilibrium volume V 0 , the bulk modulus B 0 and the first derivative of the bulk modulus B 1 . The ∆-factor, that is reported in units of meV/atom, gives an overall estimate of the discrepancy between PSPs and all-electron results in terms of these structural properties. However, as originally noticed by F. Jollet et al. [23], the ∆-factor is a stiffness dependent quantity, being proportional to B 0 . Indeed, very soft materials, as for example the noble-gas solids, are generally associated with small values of the ∆-factor, even for significant volume differences. On the other hand for very hard materials the opposite situation occurs, i.e. small volume differences give rise to large values of the ∆-factor. To solve this problem the alternative ∆ -factor was introduced, which corresponds to a ∆-factor "renormalized" to reference values of V 0 and B 0 as described in Ref. [23]. In Fig. 1 it is shown how the ∆ -factor is indeed very strongly correlated with the differences in equilibrium volume, δV 0 , obtained from the equations of state (note that the data for ∆-factor, ∆ -factor and δV 0 reported in the scatter plots of Fig. 1 are taken between pairs of PSPs for the elemental crystals tested in this work and not between PSPs and all-electron results). Instead the ∆-factor is scattered along a wide range of slopes that correspond to different values of the bulk modulus B 0 : small (large) bulk moduli correspond to small (large) slopes. In our protocol we adopt the ∆-factor (and not the ∆ -factor), as much more ∆-factor reference data has been produced in the literature [6] for several PSP libraries and codes. However, for completeness, we compute also the renormalized ∆ -factor and report this value as well. Broadly speaking two compared equations of state can be considered undistinguishable if the ∆-factor is smaller than 1 meV/atom (valid for most of the elemental solids but with the notable exception of very soft materials) or if the ∆ -factor is smaller than 3 meV/atom, where the latter corresponds to a variation in the equilibrium volume of less than 0.5% for all elemental solids (see Fig. 1). The reference all-electron results of the equation of states chosen in this work are the ones of the WIEN2k code [24] reported in Ref. [6], with the exception of the rare-earth nitrides for which we use the WIEN2k results reported in Ref. [18]. All PSP calculations needed for the ∆-factor estimation are performed at the reference wavefunction cutoff of 200 Ry using a dense Monkhorst-Pack [25] k-grid of 20 × 20 × 20 and a Marzari-Vanderbilt smearing [26] of 2 mRy. Magnetism is included for the equations of state of oxygen and chromium (antiferromagnetism), manganese (antiferrimagnetism) and iron, cobalt, nickel and the rare-earth nitrides (ferromagnetism).
Within the SSSP testing protocol we study the convergence of four different quantities as a function of the wavefunction cutoff E c , i.e. of the number of plane-waves used in the expansion of the Kohn-Sham states. The tested quantities are phonons frequencies at the zoneborder, cohesive energies, pressure, and band structures. All the calculations are performed on the ground-state structures of elemental crystals at 0 K, as provided in Ref.
[27], with the exception of fluorine for which the SiF 4 structure is used because of convergence issues of the elemental fluorine structure and of lanthanides that are not included in the test set of Ref.
[27] and for which the nitride structures of Ref. [18] are used. In total we test 85 different elements of the periodic table. In all PSP frameworks, a plane-wave representation of the charge density requires a cutoff, E ρ , higher than the wavefunction cutoff, E c [28]. Typically, convergence tests are performed by varying E c and keeping the dual, i.e. the ratio E ρ /E c , fixed. For instance, in the NC scheme the charge density is simply the modulus squared of the single-particle wavefunctions, summed over all the electrons, and in reciprocal space it reads: where sums run over the occupied bands with index n, Bloch vectors k and reciprocal lattice vectors G . Hence, the largest G-vector appearing in the charge density has modulus twice as large as the largest G-vector appearing in the wavefunction and, as plane-wave energies scale quadratically, the dual should be equal to 4 to guarantee that all Fourier components are represented. So for NC PSPs we always adopt a dual of 4, although it is known that in some cases calculations could be efficiently converged with a lower dual. In the PAW and US formalisms, wavefunctions are designed to be slowly varying in real space (i.e. to be soft) requiring substantially fewer plane-waves to be represented with respect to NC pseudo-wavefunctions. The price for working with a reduced basis set is the additional complexity in deriving expressions for observables, including the charge density which cannot be simply computed using Eq. 1. However the charge density ρ is the fundamental quantity in DFT. The physics of the problem and the energy functional determine the spatial variation of ρ and thus the cutoff E ρ , independently of the PSP scheme that is adopted. Fundamentally, that is why PAW and US PSPs require relatively high duals, where according to common knowledge a choice of a dual of 8 is usually rea-sonable to efficiently achieve good precision. However, the convergence patterns at duals equal to 12 and 16 are also checked for selected elements that show a particularly high wavefunction cutoff, namely manganese, iron, cobalt, hafnium and oxygen, as also suggested by a convergence study on Fe [29]. The convergence patterns are obtained at fixed duals as specified above and E ref c = 200 Ry as the reference wavefunction cutoff [30]. All the quantities are considered as differences with respect to the corresponding reference value calculated at E ref c . An example of the calculated convergence pattern plot is shown in Fig. 2 for the case of palladium. We choose 200 Ry as the reference wavefunction cutoff because, for each element, all the quantities tested in the SSSP testing protocol typically convergence well before that value for at least one PSP (with radon being the only exception, as discussed in the Supplementary Figure 1). Therefore, even if for some hard PSPs the convergence plots could be marginally different by using a larger value for E ref c , those modifications would not be relevant for the conclusions of our work and, in particular, the selection of the SSSP libraries would not be affected. We perform all the tests on the elemental crystals using a relatively coarse 6 × 6 × 6 Monkhorst-Pack k-grid (except for oxygen and all the lanthanides where a 10 × 10 × 10 k-grid is used instead) because in our protocol for convergence we are not directly interested in the absolute values of the tested quantities but rather on their difference with respect to the reference values computed at E ref c . We also disregard spin-polarization in all the convergence tests but we have verified for the magnetic structures that the convergence patterns are not substantially altered by the inclusion of magnetism (see Supplementary Figures 2-7).
Phonon frequencies. The convergence of vibrational properties of elemental crystals is performed by calculating, within the framework of DFPT, the phonon frequencies at the zone-border of the Brillouin zone, i.e. at the point Q = ( 1 2 , 1 2 , 1 2 ) in relative coordinates of the reciprocal lattice vectors. While the ∆-factor test is related to the structural and elastic properties of the system, by considering phonon frequencies at the border of the Brillouin zone we have access to information related to both acoustic and optical modes. The number of phonon frequencies depends on the number of atoms in the unit cell, and so on the element under investigation. In the convergence pattern plots of the SSSP testing protocol we condense the information related to the several phonon frequencies into a single number δω. It is defined as a relative average deviation (in percentage) among all the phonon frequencies ω i calculated at Q for each wavefunction cutoff E c : where N is the total number of phonon frequencies. The maximum relative deviation is similarly defined as and it is represented as an half error bar in the convergence pattern plots. If the highest phonon frequency ω max of an elemental crystal at Q is smaller than 100 cm −1 at E ref c , the absolute average deviation and the corresponding maximum deviation are computed instead of the relative ones, since a precision of a few cm −1 is often the reasonable target for a DFPT calculation.
Cohesive energies. We investigate the convergence of the energy difference between the crystalline solid and the corresponding individually isolated atoms, i.e. the cohesive energy of the elemental crystals. Since periodic boundary conditions are used in the calculations, the isolated atom is placed in a cell of lattice parameter equal to 12Å to avoid spurious interactions with the periodic images. The quantity δE coh considered in the SSSP testing protocol is defined as the absolute difference between the cohesive energy at a given cutoff E c and the one at the reference wavefunction cutoff E ref c , i.e. 200 Ry (in units of meV per atom).
Pressure. We evaluate the convergence of the stress by computing the hydrostatic pressure, which is defined as P = 1/3 Tr(σ), where σ is the stress tensor. Rather than checking convergence directly on the pressure itself (the magnitude of which depends strongly on the stiffness of the material) we evaluate it through its conversion into an equivalent volume. This allows the definition of a stiffness-agnostic and hence material's independent convergence criterion. Starting from the Birch-Murnaghan equation of state for the pressure fitted on the reference all-electron calculations we define the deviation volume V as the one closest to the equilibrium volume V 0 such that is the residual pressure of a calculation performed at the cutoff E c . With this definition, fully converged values of pressure give δP = 0 and therefore V = V 0 . Once V is known we can eventually find the relative volume deviation (in percentage) due to the residual pressure: δV press = (V − V 0 )/V 0 , which is the quantity considered in the SSSP testing protocol.
Band structure. The tests discussed so far deal with ground-state quantities only, computed either using DFT or DFPT. However, PSP calculations are often employed to study optical, transport and other properties that involve charged or neutral excitations. The majority of excited-states calculations are based on many-body perturbation theory (MBPT), e.g. G 0 W 0 and self-consistent GW [31][32][33], the Bethe-Salpeter equation (BSE) [34], or dynamical mean field theory (DMFT) [35], and are performed on-top of a DFT calculation, which provides the starting point for both self-consistent and one-shot approaches. Hence, we include band structures in our testing protocol, taking into account both the occupied bands and some of the lower lying unoccupied bands. Here, we outline a protocol for performing both convergence tests and verification of band structures by defining a bands distance (a similar idea has been proposed independently in Ref. [36]). The aim is to quantify how much two band structures "differ" by introducing a simple and computationally inexpensive metric in the band structures space. We call our bands distance η and consider two cases that are distinguished solely by the number of bands taken into account. The η v (or "eta valence") considers the occupied bands only, while in the η 10 (or "eta conduction 10") all the bands up to 10 eV above the Fermi level are considered. We always use a robust (0.3 eV) Fermi-Dirac smearing to deal with partially occupied bands, while to compute η v for insulators we use no smearing. We choose a 6 × 6 × 6 uniform k-grid, in the full Brillouin zone and with no symmetry reduction. Choosing a high-symmetry path could result in an unsatisfactory arbitrary choice, as different recipes for the standardization of paths have been introduced in the recent literature [37,38] and interesting features of the band structure may occur far from the high-symmetry lines (such as Weyl points) [39,40]. A uniform mesh is also more appropriate from the point of view of electron's nearsightedness [41]: if the energy eigenvalues are known on a sufficiently fine uniform kpoints mesh, it is possible to get an exact real-space representation of the Hamiltonian in a Wannier function basis [42] and then interpolate to an arbitrary fine mesh. Let us suppose we have two sets of bands A nk and B nk ; we define the distance between the two sets of (valence) bands as wheref , σ), where N el is the number of electrons. In order to properly align the two sets of bands, η v is defined as the minimum with respect to a rigid energy shift ω. We now consider also the low lying conduction bands by introducing η 10 , defined as in Eq. 5 but with a Fermi level up shift of 10 eV. In this way, η 10 measures the bands distance of the valence bands plus the conduction bands up to 10 eV above the Fermi energy. Finally, we also take into account the possibility that significant differences between band structures may occur only in subregions of the Brillouin zone or in small energy ranges. After computing the η, we check the slowest converging band by computing max η, defined as and request that is has to be converged with a slightly higher threshold than η itself.
In the SSSP testing protocol we use η 10 and max η 10 (in units of meV) as criteria to quantitatively study the convergence of band structures.

SSSP selection criteria
We discuss now the selection criteria used to build our optimal PSP libraries, namely the SSSP efficiency and SSSP precision libraries (version 1.1). As mentioned in the introduction, our primary goal is to define tested PSP libraries with a focus on efficiency and precision for high-throughput calculations and to suggest converged wavefunction cutoffs. The main idea behind the SSSP precision library is to provide the PSPs that are the closest to all-electron calculations in terms of ∆-factor computed at the reference wavefunction cutoff E ref c , without much consideration on the computational cost and the wavefunction cutoffs actually needed to converge all relevant quantities. On the other hand, the SSSP efficiency library is designed for practical applications that should remain affordable, and therefore PSPs are chosen such that wavefunction cutoffs are as low as possible while keeping the precision reasonable. The selection criteria are listed in Table 2. For SSSP efficiency, when possible we select PSPs with a rather small ∆-factor (below 1 meV/atom). The phonons δω should be converged within 2% (or within 2 cm −1 if the highest phonon frequency is smaller than 100 cm −1 ), the cohesive energy δE coh within 2 meV/atom, the pressure within 1% for δV press (i.e. 0.33% on the lattice parameter of a cubic crystal) and the band structure within 10 meV for η 10 and within 20 meV for max η 10 . For the SSSP precision, the criteria are slightly stricter (see Table 2) and we systematically opt for the PSP with the smallest ∆-factor. Therefore the wavefunction cutoffs of the SSSP precision are typically higher than the ones proposed for the SSSP efficiency. We underline here that in a few difficult cases the SSSP libraries are built following these criteria as general guidelines and not using these as strict rules. In practice, this means that the PSPs are chosen one-by-one through human inspection and not with an automatic procedure. This flexible approach is necessary because the convergence of some of the tested quantities is sometimes slow and/or irregular. For example it can happen that the selection criteria are not all together satisfied at a reasonable wavefunction cutoff for any of the PSPs of a given element or that the convergence patterns show outlier data points or oscillations. A clear example of this situation is given by the extremely soft noble-gas elemental solids for which the convergence of the tested quantities, in particular of the stress tensor and of the phonon frequencies, can be very noisy due to numerical instabilities. In these cases it is therefore necessary to make compromise choices that can sacrifice or increase some of the thresholds imposed by the SSSP selection criteria, if no other possibilities are available or in order to keep the computational cost reasonable.

Ghost states
We use the bands distance η 10 defined above not only for the convergence tests but also to compare the band structures of the tested PSPs for all the elemental crystals considered. However different PSPs are often generated with different combinations of semi-core states in the valence band. Hence, we compare only the bands they have in common, by taking the minimum number of electrons of all the set and cutting the exceeding low-energy bands accordingly [43]. By means of this additional criterion it is possible to automatically detect ghost states [44] in a PSP in the valence and in the conduction up to the chosen threshold (here 10 eV above the Fermi energy), as they are signalled by extremely large values (of the order of eV or more) of the bands distances when computed with respect to other ghost-free PSPs (see Fig. 3 for an example). A list of the tested PSPs having ghost states in the empty conduction bands is reported in Table 3. However, it is worth noting that standard DFT calculations for ground-state properties are unaffected by ghost states above the Fermi level. Nonetheless they could be a possible source of problems for applications related to excited-state properties (e.g. in MBPT calculations such as GW or BSE). As expected, none of the PSPs considered has ghost states in the valence as they would give unphysical results also for ground-state properties and they would be easily spotted. We stress here that the bands distance could in principle also be used for verification studies because it would allow for a quantitative comparison of PSPs band structures with reference all-electron band structures.

Correlations among tested quantities
Before giving the list of PSPs chosen by following the SSSP testing protocol, we show with an analysis of our results that an accurate selection of PSPs for generic applications in computational material science needs several independent criteria to be satisfied, based on the estimation of different physical properties. This is done also with the purpose of furnishing an a posteriori justification of the protocol we established for PSP testing. In particular we show that PSPs that give very similar results for a certain tested quantity can give, in a non-negligible number of cases, significant discrepancies in the estimation of some other quantity. We compare the PSP results for the physical properties considered in the SSSP testing protocol by calculating the discrepancies between all the available PSPs for a given element, using the data obtained at the reference wavefunction cutoff of 200 Ry. An example is shown in Fig. 4 where the differences between the equations of state (through the use of the ∆-factor) and between band structures (through the use of η 10 ) are compared for the case of sodium. We notice that, even if all the sodium PSPs we considered produce very similar equations of states-with ∆-factors always smaller than 0.3 meV-the band structure of a particular PSP shows instead substantial differences up to η 10 = 65 meV. More generally, we can study correlation between pairs of quantities by looking at scatter plots where the differences between all possible couples of PSPs for all the 85 elements are considered. In Fig. 5 we show as an example the correlation plot between the equation of state (obtained through both the ∆-factor and the ∆ -factor) and pressure (δV press ), cohesive energy (δE coh ), highest phonon frequency (δω max ) and valence band structure (η v ). Fig. 5 shows how for all such properties correlation is very weak, suggesting that the precision of a PSP is property dependent. However, we also notice (see Fig. 5) that the correlation between the ∆ -factor and δV press is higher than for the ∆-factor. This observation can be rationalized in terms of the ∆ -factor renormalization, that provides an estimate of the difference between two equations of state that is more material-independent and straightforwardly related to volume differences (see Fig. 1 and discussion therein). From our results we conclude that there is no strong correlation between pairs of tested quantities entering our selection criteria, hence the similarity between PSPs is strongly property dependent. Up to now, the ∆-factor is the only verification test present in the literature to assess the precision of DFT calculations of solids. Given the small correlations among the SSSP criteria, we stress here the importance for the electronic structure community to head for the creation of an heterogeneous set of validated all-electron reference data, which would ideally include other physical properties beyond the equation of state, such as phonons or band structures (such effort is currently coordinated by S. Cottenier). This would allow the extension of the available PSP verification tests beyond the ∆-factor for elemental crystals, potentially improving the assessment of PSPs precision. Similarly, we show that the convergence with respect to the wavefunction cutoff of a given physical property usually has a different and uncorrelated behaviour if compared to other tested quantities (see Fig. 6 for a comparison of the wavefunction cutoffs at which two tested quantities in the SSSP testing protocol reach the required precision for each PSP). Indeed, differences in the mathematical expression adopted and/or in the code-specific implementations that are needed to compute the tested quantities can result in different and independent convergence patterns, so that each quantity reaches the required precision at different wavefunction cutoffs. For example, the derivatives involved in the calculation of the stress tensor or the phonon frequencies either through direct, finite differences methods or linear response theory, can introduce different numerical noise and display a slower convergence if compared to other properties, such as equations of state or band structures, that do not require calculations of derivatives. In general, it is therefore necessary to study the convergence of each relevant quantity separately, in order to correctly estimate the optimal number of plane waves that gives results converged within the required precision for all the properties of interest.

Exchange-correlation functionals
Our study on the precision and efficiency of PSPs is restricted to the PBE functional as it is among the most popular ones in the electronic structure community and the only functional for which a verified set of reference all-electron results for solids exists [6]. So the ∆-factor test for verification can be performed, at the moment, only with the PBE functional. Still, it is worth to comment on the transferability of the convergence tests performed in the SSSP testing protocol among different functionals. For this purpose, we consider the revised PBE GGA for solids, namely the PBEsol [45] functional, and one of the most widely used functionals for the local density approximation (LDA), i.e. the PZ [46] functional. By testing some elemental crystals for the GBRV library (see Fig. 7 for the case of Ga and the Supplementary Figures 8-9 for a few more systems), we find that the convergence patterns turn out to be very similar if the PSPs are generated with the same atomic parameters (such as electronic configuration, cutoff radii, etc.), thus showing a good transferability of the convergence tests among different local and semi-local functionals. However, performing consistent tests for the transferability of more complex and nonlocal functionals, such as SCAN [47] (meta-GGA) or HSE [48] (exact-exchange), is less straightforward because, as of today, no PSP libraries for these kind of functionals exist. Indeed, in these cases, the common approach followed in all PSP DFT codes is to use PSPs generated with local or semi-local functionals, e.g. PBE, and to then "switch" to the complex functional, e.g. HSE, when performing the PSP DFT calculation including the valence electrons only. Tests performed following this approach, although useful, would not ensure the transferability of the functional under investigation and more extensive and consistent studies on the subject are therefore left to future work, in particular when HSE or SCAN PSP libraries, to only cite a few possible examples, will be available. Table 4 and Table 5 show the two SSSP libraries, efficiency and precision (version 1.1), selected according to the SSSP selection criteria specified above. The suggested wavefunction cutoffs (in Ry) and the duals are also indicated for each PSP chosen. The SSSP periodic table with all PSPs, wavefunction cutoffs and duals is also accesible interactively online on the Materials Cloud platform [11] (see Fig. 8). The average suggested wavefunction cutoffs of the SSSP efficiency and SSSP precision over all the 85 elements tested are 44 Ry and 56 Ry, respectively. A dual of 8 has been used for all PSPs except norm-conserving ones where a dual of 4 is used, and iron and manganese for which a dual of 12 is suggested. The SSSP efficiency and SSSP precision have small average ∆-factors of 0.44 meV and 0.33 meV, respectively (where the average is performed over all elements tested excluding the 15 rare-earth nitrides, following the recipe of Ref. [6]). The SSSP libraries have already proven to be a reliable tool for a number of computational studies: for instance the beta version (called version 0.7) of the SSSP libraries have enabled the high-throughput computational exfoliation of two-dimensional materials [49] and have supported the combined experimental and theoretical study of catalysts for oxygen evolution reaction [50]. On a more general level, apart from the SSSP testing protocol and libraries, our work provides a database of verification data and convergence tests that facilitates the optimal choice of PSPs and wavefunction cutoffs for custom applications. For example, some physical properties may be implemented only for some PSP types (typically only NC) or some applications may require convergence of just a subset of the quantities that we consider in the SSSP testing protocol. By a look at our plots and data, see for instance the condensed plot for palladium shown in Fig. 2, a user can quickly select the optimal PSP and wavefunction cutoff tailored for the specific application.

SSSP libraries
In summary, we propose an extensive testing protocol for PSPs to investigate precision and performance of several NC, US and PAW PSP libraries that are publicly available. We incorporate in the SSSP protocol a verification part, based on the ∆-factor, and an efficiency part, based on the plane-wave cutoff convergence tests for phonon frequencies, cohesive energies, pressures, and band structures. Leveraging the SSSP protocol, we identify two optimal PSP libraries, named SSSP efficiency and SSSP precision (version 1.1), that provide thoroughly tested and precise PSPs for 85 elements of the periodic table, selected from publicly available PSP libraries [16][17][18][19][20][21][22], for which the original authors should be acknowledged. Our effort not only is particularly relevant for high-throughput computational materials screening, where the right compromise between precision and computational cost is essential, but it substantially contributes to set high the bar of the quality of PSP calculations of solid-state materials. Building on the invaluable work behind all the PSP libraries we considered, we provide a systematic survey of PSP quality across multiple physical properties and multiple libraries and techniques (NC, US and PAW). Our work shows how the assessment of both precision and efficiency of pseudopotentials is strongly property-dependent and requires a multi-dimensional quality gauge, pointing to the need of a verification standard in the computational solid-state community. Given the importance of PSP calculations in modern materials science, nanotechnology, chemistry and physics, our findings call for more verification efforts aimed at increasing precision and efficiency of computed quantities that are routinely used to discuss novel physics, to help interpreting experiments or even to discover and design novel materials. In particular, we hope that this work will stimulate further investigations in the all-electron community, that ideally would provide more reference data for an heterogeneous set of properties elaborating on our discussion.

METHODS
All the calculations needed for this work (more than 50'000 DFT and DFPT calculations) were performed with the goal to ensure reproducibility of all the data obtained, compliant with the FAIR guiding principles for data management [51]. This is the reason why we used AiiDA [10], an open-source Python infrastructure for computational science, that is specifically designed to track the provenance of data and calculations and that allow the user to implement workflows that can run complex sequences of calculations. It is therefore particularly suited for high-throughput studies, such as the deployment of the SSSP testing protocol, where a large number of simulations are involved. In practice AiiDA can prepare and submit calculations (usually to an HPC cluster) and then retrieve and store the results inside a database, all automatically. The database can be subsequently queried by the user to extract data or other useful informations. The complete SSSP testing protocol is implemented as an AiiDA workflow, called SsspWorkflow, that can run all the convergence tests and the ∆-factor verification test. The SsspWorkflow is built on top of the PwWorkflow, a very robust lower-level workflow in charge of handling all the QE simulations and that can restart calculations in case of standard QE errors or, for example, if the user-specified wall time is too small. The SsspWorkflow allows a generic user to perform all the calculations required by the SSSP protocol in a completely automatic way. With this tool, other and new pseudopotential libraries could be easily tested in the future in order to update the subsequent versions of the SSSP libraries with more precise and efficient PSPs.  . SSSP testing protocol applied to palladium. For each pseudopotential the convergence w.r.t. the wavefunction cutoff of the zone-boundary phonons (δω), cohesive energy (δE coh ), pressure (δVpress) and bands structure (η10 and max η10) is monitored. The horizontal dashed lines correspond to the thresholds of the SSSP selection criteria (efficiency or precision); here precision is shown. On the right-hand side we report the number of valence electrons of the pseudopotential (Z), the ∆-factor and the ∆ -factor with respect to the reference all-electron results and the converged value of the highest phonon frequency (ωmax). The circle marks the pseudopotential and wavefunction cutoff chosen for the SSSP library (version 1.1). All convergence pattern plots of the 85 elements tested are available on the Materials Cloud platform [11].   and, respectively, the pressure (δVpress), cohesive energy (δE coh ), highest phonon frequency (δωmax) and valence band structure (ηv). No strong correlations among the quantities tested are observed but we note that the ∆ -factor is more correlated to δVpress than the ∆-factor. Figure 6. Comparison of the PSP wavefunction cutoffs (circles) selected strictly following the SSSP efficiency selection criteria (see Table 2) for the phonon frequencies (δω) and the cohesive energies (δE coh ). The colormap for the circles corresponds to the number of PSPs that have that pair of wavefunction cutoffs for the two quantities tested. We find that, in general, the convergence of one quantity does not imply convergence of the other.    Fig. 2) to identify the pseudopotential libraries.