{'fraction_non_alphanumeric': 0.02942672382474751, 'fraction_numerical': 0.009663590786892393, 'mean_word_length': 5.026269702276707, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reproducibility of Hybrid Density Functional Calculations for Equation-of-State Properties and Band Gaps

Hybrid density functional (HDF) approximations usually deliver higher accuracy than local and semilocal approximations to the exchange-correlation functional, but this comes with drastically increased computational cost. Practical implementations of HDFs inevitably involve numerical approximations -- even more so than their local and semilocal counterparts due to the additional numerical complexity arising from treating the exact-exchange component. This raises the question regarding the reproducibility of the HDF results yielded by different implementations. In this work, we benchmark the numerical precision of four independent implementations of the popular Heyd-Scuseria-Ernzerhof (HSE) range-separated HDF on describing key materials' properties, including both properties derived from equations of states (EOS) and band gaps of 20 crystalline solids. We find that the energy band gaps obtained by the four codes agree with each other rather satisfactorily. However, for lattice constants and bulk moduli, the deviations between the results computed by different codes are of the same order of magnitude as the deviations between the computational and experimental results. On the one hand, this means that the HSE functional is rather accurate for describing the cohesive properties of simple insulating solids. On the other hand, this also suggests that the numerical precision achieved with current major HSE implementation is not sufficiently high to unambiguously assess the physical accuracy of HDFs. It is found that the pseudopotential treatment of the core electrons is a major factor that contributes to this uncertainty.


Introduction
The reproducibility of the results, both experimental and computational ones, is essential for scientific researches. In computational materials science which heavily relies on computer simulations, this is particularly a concern. Ideally, one would like to see that the same results be obtained regardless which computer code is used. As one of the most widely In recent years, the hybrid density functionals (HDFs) are getting more and more popular in computational materials science. By incorporating a fraction of Hartree-Fock exchange (HFX), HDFs successfully addressed one big issue of local-density approximation (LDA) and GGA-type functionals -the severe underestimation of the band gaps of semiconductors and insulators. Popular HDFs include B3LYP, 3,4 PBE0 5 and HSE. [6][7][8] Among these, the HSE functional belongs to the screened HDFs in which the Coulomb interaction is split into short-range and long-range components and only the HFX based on the short-range Coulomb interaction is mixed into the exchange-correlation (XC) functionals. Incorporating the HFX only for the short-range part of the Coulomb interaction is arguably better suited for describing narrow-and middle-gap semiconductors than global HDFs. 6 In practice, the HSE functional shows remarkable performance not only for band gaps, but also for structural and cohesive properties. 9 The price to pay, however, is the significantly increased computational cost, arising from the evaluation of the short-range HFX, compared to GGAs.
In particular, using the localized variant of the RI approximation, one can achieve a linear scaling of the computational cost for building the HFX matrix with respect to the system size N . [23][24][25] The localized RI (LRI) based HDFs have been implemented within the Fritz Haber Institute ab initio molecular simulations (FHI-aims) code package 26 and recently also in the atomic-orbital-based ab initio computation at UStc (ABACUS) code package. 27,28 Besides these efforts, Lin developed an adaptively compressed exchange operator (ACE) 29 scheme, which, together with a two-level self-consistent field iterator procedure, can significantly accelerate the HDF calculations using PW basis sets. Such a scheme has been implemented in the DGDFT 30,31 and the Quantum ESPRESSO (QE) 32,33 code packages. Other techniques to speed up the HFX calculations include the exploration of the localization offered by Wannier representation, [34][35][36] intepolative separable density fitting approaches, 37,38 and representing NAOs in terms of Gaussian-type orbitals, 39 etc.
As alluded to above, the numerical precision and the reproducibility of PBE implementations in mainstream DFT codes have been established via the ∆-value project. 2,40 However, this conclusion may not be automatically extended to HDFs as implemented in these codes.
In fact, the numerical precision of HDF calculations does not only rely on the treatment of electron-ion interactions (pseudopotentials vs. all-electron) and basis set types (e.g., plane waves vs. atomic orbitals), but also on the numerical approximations involved in the evaluation of the HFX potential and energy. To our knowledge, a systematic assessment of the numerical precision of different HDF implementations has not yet been reported. In this work, utilizing the HSE functional available in four computer codes, including VASP, FHIaims, QE, and ABACUS, whose HDF implementations differ considerably in algorithmic and numerical details, we assess the reproducibility of HSE results by computing the pairwise ∆ values among these four codes, based on the energy-volume curves of 20 binary compounds. These 20 compounds have been widely used as a test set for evaluating the performance of HDFs. 9,24,41-44 Such a ∆-value based assessment will be corroborated by additional properties including the lattice constants, bulk moduli, and band gaps. Naturally, benchmarking the numerical precision of HDF implementations over all mainstream DFT codes and across different types of materials requires community-wide efforts and goes beyond the scope of the present paper. Nevertheless our work provides a solid basis that will hopefully stimulate more benchmark efforts in this direction.

Methods
In this section, we recapitulate the essentials of the numerical frameworks and computational setups employed in the four codes considered in this work. First of all, different strategies are used to treat the effects arising from the nuclei and core electrons: The VASP code employs the PAW method whereas QE and ABACUS make use of the SG15 optimized normconserving Vanderbilt-type (ONCV) pseudopotentials. 45,46 In contrast, FHI-aims treats all electrons on an equal footing. Scalar relativistic effects are included either in the PAW or norm-conversing pseudopotentials (NCPP), or via the "atomic ZORA" scheme 47,48 in the all-electron case. The spin-orbit coupling effect is not included in the HSE calculations here.
Regarding the basis sets, VASP and QE both use PWs to expand the Bloch wavefunctions, whereas numerical atomic orbital (NAO) basis sets are used in the FHI-aims and ABACUS codes, although different NAO construction strategies are adopted in these two codes, in order to best represent their respective all-electron (AE) or pseudo-wavefunctions. Moreover, kinetic energy cutoffs of 550 eV and 100 Rydberg (Ry) are used for the PW expansion of wavefunctions in VASP and QE calculations, respectively. For FHI-aims, the default "tight" NAO basis sets 26 are adopted, whereas the so-called "DZP-DPSI" 49 basis sets are employed in ABACUS. Additionally, 100 Ry energy cutoff is used in ABACUS for the Hartree potential calculation and for determining the uniform quadrature grid. In Table 1, we summarize the computational settings for the four codes, together with the convergence thresholds in the self-consistent calculations. Please note that these settings are roughly the default ones that are typically used in production calculations. Our intention here is not to push towards the limit of high-precision settings of each code, but rather to check the numerical precision that is achieved in routine calculations. An investigation of the influence of difficult levels of basis sets and numerical settings on the achieved precision of the results has recently been reported by Carbogno et al. 50 for three different computer codes, including FHI-aims.
As mentioned above, additional complexity is involved in the evaluation of the HFX component of the HDFs. For VASP, the standard HSE implementation as available in VASP 5.4.1 is used. As for QE calculations, we use its PWSCF module v.6.5 where the ACE technique 29 is used as the default to calculate the HFX contributions. 33 Finally, for NAObased codes, the LRI approximation is used in both FHI-aims and ABACUS for evaluating the HFX term. To render the numerical errors associated with LRI sufficiently small, an enhanced auxiliary basis set has to constructed, compared to conventional global Coulombmetric based RI scheme. 22 In FHI-aims, this is achieved by adding a hydrogen-like 5 g function with an effective nuclear charge Z = 6 for each element that is solely used to generate extra auxiliary basis functions (ABFs). 22 In ABACUS, additional optimized ABFs can be included to complement the "on-site" ones to increase the numerical accuracy. 24 Furthermore, an efficient pre-screening procedure is invoked to filter out negligibly small contributions to achieve a linear-scaling build of the HFX. The settings of the ABF basis sets, and the pre-screening parameters for ABACUS HSE calculations are given in the Supporting Information (SI) (cf. Table S1 and discussions therein). Finally, the band gaps are estimated on a discrete 8×8×8 k-point mesh, and the experimental lattice parameters are used for all the band structure calculations. Table 1: Numerical frameworks and computational settings for the four codes included in the present benchmark studies. "Method" describes how the effects from the nuclei and core electrons are treated. "Basis sets" tell the type of basis functions to expand the wavefunctions. E cutoff : the kinetic energy cutoffs for the PW expansions (in case of VASP and QE), and for the Hartree potential calculation and for determining the uniform quadrature grid (in case of ABACUS). η etot , η ρ denotes the convergence threshold for the self-consistency cycle, based on the criteria of total energy, charge density, respectively. (n/a: this metric is not used in the particular code or not set in this benchmark test).

Equation of states
We start by examining the pairwise ∆ values evaluated in between the four codes, given by an average of the individual ∆ values for the 20 materials. Namely, where measures the root-mean-square difference between the energy-volume (E-V ) curves (with aligned energy minima) generated by two codes i and j for a material α in the test set. In Eq. 2, V 0,α is the averaged equilibrium volume of the α-th material, as predicted by the two codes. A change of the volume within a range of −6% to 6% around the equilibrium point is included in the evaluation of the ∆ value. The thus determined 4 × 4 ∆ matrix is shown in Figure 1 for both the PBE and HSE functionals. From Figure 1 Figure 1 suggests that the approach to describe the valence-ion interactions, has a bigger influence than the choice of basis sets. Indeed, a ∆ value of 3.7 meV/atom is obtained between the two NCPP-based codes -QE and ABACUS, which is smaller than the discrepancy between the two PW-based codes (5.3 meV/atom) and between the two NAO-based codes (4.9 meV/atom). Comparing the two panels in Figure 1 Table S2 of SI. Figure 2b shows  Despite these discrepancy, some general features can be drawn from this benchmark study.
For example, all four codes predict that the HSE bulk moduli of C and BN are appreciably larger than the experimental value, pointing to an intrinsic accuracy limitation of HSE for these hard materials.

Band gaps
One important success of the HSE functional is that it significantly improves over local and semi-local functionals for describing the band gaps of semiconductors and insulators, thanks to the HFX component and the generalized KS framework. Thus, the reproducibility of the calculated HSE band gaps is also an important aspect to check, in addition to the EOSbased properties. In Figure 4a, we present the MAD matrix for the band gaps, including the experimental and four sets of HSE results. The actual values of the HSE band gaps, as given by the four codes, as well as the experimental results are presented in  properties; on the other hand, this also means that the numerical precision of HSE implementation still needs to be improved to unambiguously assess its intrinsic accuracy with respect to experimental reference data. Among the computational results, it seems that the prescription for describing the ion-cores plays a significant role in determining the numerical precision of the HSE implementations. Specifically, the all-electron (FHI-aims) and PAW- solids. 52 We expect that this will also holds for the HSE functional and the benchmark set presented in this work.
As HDF calculations get increasingly popular in computational materials sciences, the reproducibility of this type of calculations becomes a prominent issue. Our benchmark calculations in this work only cover a limited set of codes and materials, yet interesting features and trends have been revealed. We hope that this work will stimulate more comprehensive investigations along this line.

Supporting Information Available
The Supporting Information is available free of charge .
• Description of the computational parameters used in ABACUS calculations, and tables containing the computational and experimental values of lattice constants, bulk moduli, and band gaps for 20 solids.
The input and output files of ABACUS calculations are available from the authors upon reasonable request.