Plane-wave pseudopotential implementation and performance of SCAN meta-GGA exchange-correlation functional for extended systems

Development of new exchange-correlation (XC) functionals in density functional theory (DFT) remains one of the most active research areas in theoretical chemistry and theoretical condensed matter physics.^1,2 Perdew and co-workers have advocated for a systematic approach in which more accurate XC functionals are constructed by satisfying larger numbers of exact physical behaviors.^2(c),5 This is done through introducing increasingly more complicated dependence on the density in the XC functional as often depicted in the “Jacob’s Ladder” of DFT.^2(c),2(g) One of the most recent developments in this non-empirical construction of XC functionals is the strongly constrained and appropriately normed, SCAN functional which is a meta-GGA functional that has an additional dependence on the kinetic energy density (KED) beyond GGA functionals.^8,21 The SCAN functional was designed so that all known constraints are satisfied at the semi-local level, and its performance appears to be quite promising for a wide range of problems in molecular and condensed matter sciences.^8,9

In modern first-principles calculations based on DFT, in addition to the physical approximation of using a particular form of XC functional, implementation of DFT has become quite important for its practical applications. Indeed, the numerical aspect has become more widely appreciated nowadays as discussed in the recent work¹⁰ that addresses the reproducibility of different DFT implementations for solid-state systems. The numerical consideration is of great importance especially with advanced XC functionals because novel techniques and further approximations are often introduced in the implementation and also because a greater accuracy is expected at the same time. In addition to the usual issue regarding how Kohn-Sham (KS) wavefunctions and electron density are represented in terms of basis functions (e.g., planewaves and Gaussians),¹¹ further approximations are often introduced in the practical implementation of the XC functionals with explicit KS orbital dependence.¹ For instance, techniques such as unitary Wannier transformation,¹² resolution-of-identity,¹³ and recursive subspace bisection approach¹⁴ are often used to efficiently deal with the Hartree-Fock exchange in hybrid XC functionals. For meta-GGA functionals,¹⁵ further dependence on the kinetic energy density (KED) often makes calculations more sensitive to numerical parameters such as the integration grid,¹⁶ and the KED dependence of the XC potential requires more attention in its numerical implementation.¹⁷ 

The implementation based on planewaves and pseudopotentials has been widely used in condensed matter sciences¹⁸ especially because the periodic boundary conditions are naturally incorporated while other implementations such as those based on Gaussians¹⁹ and numerical atomic orbitals²⁰ have rapidly emerged in the last few decades. In this work, we discuss an implementation of the SCAN meta-GGA functional²¹ in the traditional planewave-pseudopotential (PW-PP) formalism and present its performance for some extended systems of practical importance.

The SCAN functional belongs to a class of meta-GGA functionals that have an explicit dependence on the kinetic energy density, $τσ$⁠,¹⁵ 

where $𝜀xc$ is the exchange-correlation energy density. The kinetic energy density is defined as

where $ψσi$ are the KS single-particle wavefunctions, and the summation is taken over all occupied states. In the SCAN meta-GGA functional,²¹ the kinetic energy density (KED) is used to define a new dimensionless variable,

where $τW(r)$ is given by the Weizsacker KED, the single-orbital limit $τW(r)=∇ρ(r)2∕8ρ(r)$⁠, and the uniform-density limit $τunif(r)=310(3π2)2∕3ρ(r)5∕3$⁠. The spin dependence is omitted from equations for conciseness in the rest of the paper. The exchange-correlation energy density is generally expressed concisely using the density, the spin polarization, the reduced density gradient, and now the new variable α. This reduced KED variable allows the SCAN XC functional to distinguish different types of bonds.²¹ 

As in the Car-Parrinello extended-Lagrangian formalism,^18(b),22 one can find the KS wavefunctions that minimize the energy functional by taking the functional derivative with respect to the complex conjugate of the wavefunctions. By taking a unitary transform of the resulting wavefunctions, one can express the resulting equation in the diagonal form

where $𝜀i$ are KS eigenvalues, and $υ$ and $Vh$ are external potential and Hartree potential, respectively. If the XC functional is dependent on only the electron density and its gradient, we find

and the proper KS equation with the local multiplicative XC potential, $δExc[ρ]δρ(r)$⁠, is readily recognizable. However, in the context of meta-GGA functionals (and also for hybrid XC functionals), the direct minimization with respect to the single-particle wavefunctions does not yield such a local multiplicative XC potential. For the SCAN functional, a meta-GGA functional with the dependence on the kinetic energy density, we have

in which the last term is unique to the meta-GGA, and it is not present in LDA or GGA functionals. For obtaining the “proper” KS approach to DFT, the optimized effective potential (OEP) approach²³ should be employed to obtain the local multiplicative XC potential in principle. However, such a numerical approach would be computationally too demanding and not practical for the PW-PP implementation even with some approximations such as the Krieger-Li-Iafrate scheme.²⁴ Thus, we here follow the so-called generalized Kohn-Sham (gKS) scheme^9,25 in the present implementation as often done for hybrid XC functionals with the Hartree-Fock exchange by minimizing the functional with respect to the KS wavefunctions. This gKS approach has been discussed in the literature under different names for meta-GGA implementation.²⁶ In our earlier work with the TPSS meta-GGA functional,²⁷ we also followed this gKS approach. Note that this gKS scheme is generally no less accurate than the proper KS approach with the local multiplicative XC potential as obtained in the OEP scheme.⁹ 

We present here the details of our meta-GGA PW-PP implementation in the context of the SCAN functional since it was not discussed in our earlier work.^27(a) Implementation of the meta-GGA functional within the PW-PP framework requires two key considerations beyond the typical GGA functional: (1) evaluation of the τ-dependent part of XC potential with the PW basis functions and (2) PP generation by inverting the atomic KS equation with the τ-dependent XC potential. The total energy for the meta-GGA functional can be conveniently written as

where the τ-dependent part of XC potential in the last term is defined as

The mathematical form of the exchange-correlation energy density is given in the original SCAN publication.²¹ The exchange part is given in terms of the exchange enhancement to the uniform electron gas limit, and the correlation part involves a weighting function with the reduced kinetic energy density. The last term in Eq. (7) is the additional ingredient that is unique to meta-GGA XC functionals, and it can be conveniently expressed also as

The PW-PP implementation of meta-GGA XC functionals remains the same as that of GGA XC functionals, except for the evaluation of this last term that depends on the kinetic energy density. We focus here on the evaluation of this last term since other aspects are already well documented in the literature.^18(a),18(b) KS single-particle wavefunctions are represented in real space and reciprocal space as

where $Ω$ is the simulation cell volume and the summation over G vectors is for those planewaves whose modulus is less than the pre-defined cutoff G_cut. The wavefunction gradient is simply given by

The wavefunctions are stored in the reciprocal space as $ci(G)$⁠, and the kinetic energy density, $τ(r)=∑i12∇ψi(r)2$⁠, is evaluated in real space once the wavefunction gradient is calculated. A generalized KS equation can be solved using a number of standard approaches already implemented for DFT calculation, and the meta-GGA KS equation is

In the PW-PP implementation, essentially all terms except the XC potential are conveniently evaluated in the reciprocal space. For the SCAN meta-GGA functional, the XC potential part is evaluated in the real space as follows. To calculate the last τ-dependent term in the left hand side in Eq. (12), we first evaluate $∇ψi(r)$ in real space from the $ci(G)$ as in Eq. (11). $Vτ(r)∇ψi(r)$ is obtained in real space first, and then it is Fourier transformed to the reciprocal space. The last term in the generalized KS equation in the reciprocal space is

For calculating the total energy, the contribution of the kinetic energy density, $∫drVτ(r)τ(r)$⁠, is directly calculated in real space. In a typical PW-PP calculation, fast Fourier Transform (FFT) is used with the smallest FFT grid size, which is determined by G_cut (dependent of the simulation cell) and also by a particular FFT code employed in the calculation for performing the discrete Fourier transform (e.g., FFTW²⁸). Note that many FFT codes are machine/hardware-dependent. At the same time, the calculation of XC potential depends explicitly on this grid as it is calculated numerically on the FFT grid in real space. Generally, GGA functionals usually show very minor dependence on the FFT grid in practice. In the case of meta-GGA functionals like the SCAN, an extra term with the kinetic energy density is calculated in real space, and therefore we also need to examine the extent to which the FFT grid influences the calculated properties as discussed in Sec. III. The above formalism can be extended straightforwardly with the k-point sampling in the Brillouin zone in the usual manner.¹¹ 

An important practical consideration for the planewave-based implementation is the use of pseudopotentials for atoms so that “core” electrons need not be explicitly represented by planewaves in calculations.²⁹ In principle, these pseudopotentials should be generated for the specific XC functional used in the calculations. However, for DFT calculations with advanced XC functionals, it has become quite common to use pseudopotentials that are generated for a different and more convenient XC approximation in practical calculations (e.g., using PBE³⁰ pseudopotentials in PBE0³¹ calculations). At the same time, we need to examine the extent to which such numerical approximation influences the calculated physical properties for the present case of the SCAN meta-GGA functional.

For generating the pseudopotentials, we need to solve the atomic KS equation with spherically symmetric potential, and the KS wavefunction is written as

in the spherical polar coordinate. The angular part,

is the same as in the LDA/GGA functionals. The radial part for the meta-GGA functional is

Noting that

we can write the kinetic energy density in the spherical coordinate as

where O_nl are the occupations. We assumed here that we have spherically symmetric kinetic energy density, which is true only for closed-shell atoms in principle. To solve the radial KS equation on a radial grid, the usual logarithmic grid is employed because it is necessary to have a denser grid close to the nucleus.^29(b) 

Here, the Troullier-Martins (TM) scheme³² is followed for the pseudopotential generation. This requires matching the pseudo- and the all-electron energies and wavefunctions for the valence electrons in a reference configuration, which is typically taken to be the ground state configuration of the isolated atom. In the TM scheme, we write the radial KS wavefunction inside the cutoff radius, r ≤ r_c, as

where the polynomial function is given by

The coefficients in this polynomial function are uniquely determined by satisfying a set of constraints that are given by the TM scheme for the chosen r_c. In addition to the norm conservation of the pseudo-wavefunctions and ensuring the smoothness of the pseudopotential at the origin r = 0, the continuity of the wavefunction and its derivatives up to the 4th order at r = r_c are required in the pseudopotential generation. Ensuring the continuity of the wavefunction and the 1st derivative yields same equations as in the case of GGA functionals, then the parameterized polynomial function must have;

and

at r = r_c. For the 2nd derivative at r = r_c, we have

where $V(r)$ is defined to be the part of the potential that does not have the kinetic energy density dependence

For convenience, let us define

and its first and second derivatives

and

Then, the 3rd and 4th derivatives of the polynomial function must satisfy

and

at r = r_c. This set of constraints on the polynomial function is used to uniquely determine its coefficients. The potential and other meta-GGA terms are that of the all-electron calculation in the above equations. Once the coefficients are determined, we obtain pseudo-electronic density and pseudo-kinetic energy density from $Rps(r)$⁠. Then, the screened potential for r ≤ r_c is obtained with $Vτ(r)$ from valence electrons by inverting the radial KS equation^29(b) as

Note that $dVτ(r)dr$ needs to go to zero as fast as r does; otherwise, the last term diverges for l > 0. The final form of the psedopotentials is obtained by subtracting the contributions from the Hartree and the XC potential from the screened potential in the usual manner.^29(b) When one needs to take into account the overlap effects of density and the kinetic energy density between core and valence electrons in some elements, non-linear core correction³³ and/or projector augmented-wave (PAW) approach³⁴ might be necessary.

To study the performance of the SCAN meta-GGA functional using the planewave-pseudopotential (PW-PP) formalism, we considered three representative cases that are particularly of interest in condensed matter sciences. These are (a) crystalline silicon in the semiconducting diamond phase and the metallic beta-tin phase^8,35 and crystalline germanium, (b) martensitic phase transition energetics of phosphorene,³⁶ and (c) physisorption of a single water molecule on a graphene sheet.³⁷ We discuss how key numerical approximations influence the calculations in the PW-PP formalism for the SCAN functional. Given that the SCAN depends on the kinetic energy density in addition to the electron density and its gradient, we examine the influence of planewave (PW) cutoff as well as the FFT grid, which is used for evaluating the XC potential in real space. We also assess the influence of using pseudopotentials that are generated at a different level of XC potential. Especially with hybrid functionals like the PBE0, it has become quite common to use pseudopotentials that are generated at a different lower level of the XC approximation (e.g., PBE) in practical calculations. The SCAN functional was implemented in the Quantum-Espresso code^18(c) based on our earlier TPSS meta-GGA implementation,^27(a) and ld1 code by Paolo Giannozzi was modified as described in Sec. II to include the SCAN meta-GGA XC functional. Details of the TM pseudopotentials used for this work are summarized in Table I. The ground states are used as reference states, and the cutoff radii are taken to be the same as those found for the PBE functional in pslibrary.³⁸ For the Si, P, and Ge pseudopotentials, 3d or 4d projector channels are included as usually done.^29(b) 

The periodically repeated unit cell of crystalline silicon was used to perform the energy convergence test first as an example here. The Monkhorst-Pack k-point grids of 16 × 16 × 16 and 12 × 12 × 14 were used for sampling the Brillouin zone using the unit cell for the diamond and beta-tin silicon phases, respectively. As shown in Fig. 1, the total energy using the SCAN functional converges similarly to the case of the PBE GGA functional although the SCAN has an additional dependence on the kinetic energy density. The FFT grid dependence of the total energy is shown in Fig. 2. The smallest possible FFT mesh is automatically generated by the particular FFT library (FFTW²⁸ in this work) for a given PW cutoff energy and a given simulation cell size. The PW cutoff of 100 Ryd (for wavefunctions) yielded 50 grid points per axis for FFT. As discussed in Sec. II, the numerical evaluation of the XC potential is performed in real space, and the total energy can therefore depend on the FFT grid, especially for meta-GGA functionals like the SCAN. Figure 2 shows how the total energy changes with the increasing (denser grids) grids up to 500 points per axis for both the SCAN and PBE functionals. Note that increasing the FFT grid dramatically increases computational cost in the PW-PP implementation. Although the total energy convergence with respect to the FFT grid is more problematic with the SCAN than with the PBE functional, the absolute magnitude of the observed fluctuations in the energy can be considered negligible for most applications in this case.

Figure 3 shows the band structure of the crystalline silicon in the diamond lattice structure using the PW cutoff of 100 Ryd and the FFT grid of (50,50,50). The convergence with respect to the FFT grid was checked with the FFT grid of (300, 300, 300). The SCAN bandgap was determined to be 0.93 eV, which is significantly improved from the PBE value of 0.59 eV (see Table II). Our value using the PW-PP formalism is very close to the SCAN energy gap of 0.97 eV using the PAW^34,39 potential (at the PBE level) as reported in Ref. 8. We note, however, that Zeng-hui Yang et al. have recently reported that the SCAN bandgap using the OEP scheme⁹ (0.78 eV) is noticeably smaller than the SCAN bandgap in the gKS scheme as employed here. The comparison to the earlier meta-GGA TPSS^27(a),40 is given in Table II, and noticeable improvement is observed with the SCAN approximation over the TPSS approximation of the meta-GGA. Also, the SCAN functional does not suffer from the overestimation of the bandgap by PBE0 hybrid approximation.³¹ We tested the pseudopotential dependence by performing the SCAN calculation with the PBE GGA pseudopotential. The SCAN bandgap is somewhat smaller, yielding 0.83 eV when the PBE pseudopotential is used.

Figure 4 shows the total energy difference between the diamond and beta-tin phases of the crystalline silicon. The beta-tin phase is metallic and the diamond phase is insulating, and LDA and most GGA functionals are known to perform rather poorly in describing the energy difference between the two phases.^35(b) Recent Diffusion Monte Carlo (DMC) calculations predict the energy difference between these two phases to be 475 ± 10 meV/atom⁴¹ or 424 ± 20 meV/atom,^35(a) depending on computational details, such as the extrapolation scheme for finite-size errors. Note that these reported DMC values include a separate core polarization correction of −30 meV.^35(a) LDA and the PBE GGA significantly underestimate the energy difference, yielding 200 meV/atom and 290 meV/atom, respectively. The energy difference using the SCAN functional in this work is 480 meV/atom, which is 65 meV/atom larger than the corresponding value reported by Sun et al. based on the implementation of the SCAN functional using the PAW potential formalism. Note that the PBE functional was used for the PAW in their work.^8,39 We believe that this difference derives from the treatment of the core electrons. As shown in Fig. 4, our SCAN implementation indeed yields 440 meV/atom when the PBE GGA pseudopotential is used, instead of 480 meV/atom. While the diamond phase is not very sensitive to the pseudopotential, the beta-tin phase shows appreciable dependence. We also tested the dependence on the pseudopotential cutoff radius with r_c = 1.5 a.u., 1.8 a.u., 2.0 a.u., and we found negligible dependence as shown in the  Appendix.

The SCAN functional yields the bulk modulus of 99.30 GPa for silicon in the diamond structure (Table III), which is close to the reported value of 99.44 GPa in Ref. 8. Experimental value is 99.2 GPa.⁴ This is notable especially because the PBE GGA functional performs worse than LDA functional for calculating the bulk modulus of silicon and many other semiconducting solids,⁴ and the SCAN functional appears to remedy this shortcoming of the PBE GGA functional and earlier meta-GGA functional TPSS.^27(a),40 Using the PBE GGA pseudopotential in the SCAN calculation yields the bulk modulus of 95.03 GPa. For obtaining the bulk modulus, calculations with different lattice parameters are needed to numerically take the derivative of the total energy. We found it particularly important to use the same FFT grid for all lattice parameters because the FFT grid changes could introduce discontinuities in the total energy (Fig. 2). We checked the convergence of the bulk modulus by comparing to the calculation with the larger FFT grid of (300, 300, 300), which yields the same value of 99.30 GPa.

Given the great performance of the SCAN functional for the silicon, we also considered the notorious case of the crystalline germanium. GGA-based DFT calculations predict that the material is metallic, and also the GGA-PBE predicts the bulk modulus that is very different from the experimental value and worse than the LDA prediction (see Table III). The Monkhorst-Pack k-point grid of 16 × 16 × 16 was used for this calculation. Figures 5 and 6 show the convergence behavior with respect to plane-wave cutoff energy and FFT grid, and we report only the converged values in the followings. As summarized in Table III, the SCAN yields the bulk modulus of 73.38 GPa, compared to the PBE value of 58.73 GPa, TPSS value of 57.15 GPa, and experimental value of 75.8 GPa. The use of the PBE functional pseudopotential within the SCAN calculation yielded a value of 62.76 GPa, which is noticeably different from the SCAN result. Figure 7 shows the band structure of the crystalline germanium in the diamond lattice structure. The SCAN bandgap was determined to be 0.57 eV, which represents a great improvement toward the experimental value of 0.74 eV (see Table II). We note, however, that Zeng-hui Yang et al. have recently reported that the SCAN bandgap using the OEP scheme⁹ for the germanium is vanishingly small, unlike the SCAN bandgap in the gKS scheme as employed here. We tested the pseudopotential dependence by performing the SCAN calculation with the PBE GGA pseudopotential. The SCAN bandgap is noticeably smaller, yielding 0.19 eV when the PBE pseudopotential is used.

Phosphorene is a single layer allotrope of phosphorus. Among various phases of phospherene, black and blue phospherene phases are predicted to have a similar energy, and the martensitic phase transitions involve two transition states with modest barriers.³⁶ Here we use this particular phase transition as an example to study the performance of the SCAN functional in our PW-PP implementation. The Monkhorst-Pack k-point grid of 8 × 8 × 1 was used for sampling the Brillouin zone using the unit cell with four atoms. Figure 8(a) shows the total energy difference between the black and the blue phases as a function of the PW cutoff using the default FFT grid given by the FFT library. The energy difference is converged satisfactorily to the desired accuracy already at the PW cutoff energy of ∼40 Ry. Figure 8(b) shows the total energies of the black and blue phospherene phases and the energy difference between the two phases as a function of the FFT grid density. The PW cutoff energy of 60 Ry was used here and the default FFT grid was used as the reference for FFT grid density (FFT grid density = 1). As can be seen, the calculated energetics was not sensitive to the FFT grid in this particular case.

For the calculation of the martensitic phase transition energetics from the black phase to the blue phase (Table IV), the PW cutoff of 100 Ry was used. The SCAN functional predicts the overall reaction barrier of 0.406 eV/atom at the second transition state (TS2), and this represents great improvement from the PBE functional’s 0.335 eV/atom and TPSS functional’s 0.370 eV/atom. The SCAN functional value is comparable to the value of 0.397 eV/atom we obtained using the PBE0 hybrid XC functional. Although the energy barrier is still underestimated with the SCAN functional, it improves over other XC functionals toward the value of 0.475(3) eV/atom we obtained using DMC calculation.³⁶ As in the case with using the PBE, PBE0, and TPSS approximations, the SCAN predicts that the blue and black phosphorenes are essentially degenerate as predicted by the DMC calculation. We also examined the influence of using the PBE GGA pseudopotential in the SCAN calculation as often done for convenience when advanced orbital-dependent XC functionals are used. We found that this approximation introduced a rather small error of at most 0.03 eV/atom in this case (Table IV).

Physisorption of a single water molecule on an extended material surface like graphene is believed to arise from the interplay between electrostatic and dispersion interactions. It remains rather difficult to apply DFT calculations reliably to such a problem because of its significant sensitivity to the XC approximation. For example, the PBE GGA functional yields the adsorption energy of only 27 meV while LDA functional predicts the adsorption energy of 151 meV.³⁷ These are both quite different from the DMC prediction of 70 ± 10 meV and the random-phase approximation (RPA) of 98 meV for the adsorption energy.³⁷ Also, we note that the PBE0 hybrid and TPSS meta-GGA functionals do not improve much over the PBE functional. In our calculation using the SCAN functional, we used the simulation cell containing 60 carbon atoms with 12.28 Å × 12.762 Å in the planar directions and a 20 Å in the orthogonal direction. Γ point was used for the k-point sampling of the Brillouin zone. Figure 9(a) shows how the adsorption energy converges as a function of the PW cutoff energy using two different FFT grids: (1) the default FFT grid which is determined by the PW cutoff and the FFT library and (2) a dense converged FFT grid of (300, 320, 486). Unlike for the previous two cases of the crystalline silicon and phase transition energetics of phosphorene, the adsorption energy was found to be quite sensitive to the FFT grid. Figure 9(b) shows the total energies of the adsorbed and isolated structure of water molecule on graphene, along with the adsorption energy as a function of the FFT grid density. The PW cutoff energy of 60 Ry was used here, and its default FFT grid was used as the reference for FFT grid density (FFT grid density = 1). The figure shows that FFT grid needs to be much denser (as much as 1.5 times) than the default value in order to obtain the accurate SCAN adsorption energy.

Using the PW cutoff of 100 Ry, the adsorption energy was calculated to be 82 meV as shown in Figs. 9 and 10. This value is comparable to the adsorption energy from the DMC calculation and the RPA calculation by Ma et al.³⁷ Although the SCAN functional yields an accurate adsorption energy, the overall interaction energy profile differs noticeably from that of the RPA calculation (Fig. 10). For instance, the interaction energy at the separation distance of ∼5 Å is only about 10 meV in the SCAN calculation while the RPA value is approximately 30 meV. This is likely due to the missing long-range dispersion interaction in the SCAN approximation as discussed in the paper.⁴² When we use the recently developed SCAN + rVV10 by Peng et al.,⁴² the adsorption profile is noticeably improved for large separation distances, becoming closer to the RPA profile. However, the rVV10 addition also appears to make the adsorption energy overestimated as seen in Fig. 10. For this problem, the dependence of the adsorption energy on the pseudopotential was negligible.

We presented implementation and performance of the strongly constrained and appropriately normed, SCAN, meta-GGA exchange-correlation approximation in the planewave-pseudopotential (PW-PP) formalism using the Troullier-Martins pseudopotential scheme. We studied its performance by applying the implementation to (a) crystalline silicon in the semiconducting diamond phase and the metallic beta-tin phase, and crystalline germanium (b) martensitic phase transition energetics of phosphorene, and (c) a single water molecule physisorption on a graphene sheet. Given the much-improved accuracy over the GGA functionals and its relatively low computational cost compared to hybrid XC functionals, the SCAN functional is very promising for various practical applications of density functional theory calculations for condensed matter systems. Importantly, the SCAN formulation represents a significant improvement from the earlier TPSS formulation of the meta-GGA XC functional by Perdew and co-workers for all cases we studied here. At same time, the SCAN meta-GGA functional requires more careful attention to numerical details. The cutoff energy for converging the planewave representation of the Kohn-Sham wavefunctions for the SCAN functional was found similar to that of other types of XC functionals even though the SCAN has an additional dependence on the kinetic energy density. At the same time, we found that the SCAN meta-GGA functional shows more significant dependence on the fast Fourier transform grid, which is used for evaluating XC potential in real space in the PW-PP formalism, than other more conventional functionals like the PBE GGA. For example, in the cases of calculating the bulk modulus of silicon/germanium and the water physisorption, the FFT-grid dependence had to be checked rather carefully. We also examined the influence of using pseudopotentials that are generated at a different/lower level of XC approximation. This approximation has become quite common in practice when using advanced XC functionals with explicit orbital dependence in recent years. The accuracy of this approximation in the PW-PP formalism depends on specific properties one is interested in. For example, the energetic difference between the metallic and semiconductor phases of silicon as well as the bulk modulus of silicon and germanium is rather sensitive to this approximation while the water physisorption energy is not.

This work is supported by the National Science Foundation under Grant No. CHE-1565714. We thank the National Energy Research Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02- 05CH11231 for computational resources. An award of computer time was also provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC02-06CH11357. We thank Dr. X. Wang for helpful discussions.

Figure 11 shows how the pseudopotential cutoff radius affects the total energies of the metallic beta-tin phase (high-pressure phase) and semiconducting diamond phase of the crystalline silicon. See Table I for the pseudopotential parameters. As can be seen, the dependence is very minor.