Calculation and interpretation of classical turning surfaces in solids

Classical turning surfaces of Kohn–Sham potentials separate classically allowed regions (CARs) from classically forbidden regions (CFRs). They are useful for understanding many chemical properties of molecules but need not exist in solids, where the density never decays to zero. At equilibrium geometries, we find that CFRs are absent in perfect metals, rare in covalent semiconductors at equilibrium, but common in ionic and molecular crystals. In all materials, CFRs appear or grow as the internuclear distances are uniformly expanded. They can also appear at a monovacancy in a metal. Calculations with several approximate density functionals and codes confirm these behaviors. A classical picture of conduction suggests that CARs should be connected in metals, and disconnected in wide-gap insulators, and is confirmed in the limits of extreme compression and expansion. Surprisingly, many semiconductors have no CFR at equilibrium, a key finding for density functional construction. Nonetheless, a strong correlation with insulating behavior can still be inferred. Moreover, equilibrium bond lengths for all cases can be estimated from the bond type and the sum of the classical turning radii of the free atoms or ions.


INTRODUCTION
Modern Kohn-Sham (KS) density functional theory (DFT) 1 calculations produce a KS potential, v s (r), which, while not a physical observable, has proven useful in providing physical and chemical insight. Recently, the classical turning surfaces at the energy of the highest occupied orbital of atoms and ions have provided considerable insight into the nature of chemical bonding 2 . These turning surfaces are roughly ellipsoidal for covalent bonds, have seams for ionic bonds, and are bifurcated for van der Waals dimers. But all finite systems have densities that decay to zero far from the nuclei, whereas the interiors of real solids extended over three dimensions do not. Thus, unlike molecules, solids have the possibility of having no classical turning surfaces at all. Here we augment existing codes to calculate classical turning surfaces for a variety of extended solids.
We first find that equilibrium distances between nearestneighbor atoms or atomic ions in solids can be estimated from the sum of the classical turning radii of the corresponding free atoms or atomic ions and the bond type (metallic, covalent, ionic, or van der Waals), much as found earlier for molecules. But of greater interest is the nature of the turning surface itself in solids. Physical intuition suggests that, in the limit of extreme compression, solids become metallic, while, in the less-physical limit of extreme expansion, they become insulators. In the former case, there are no classical turning surfaces, i.e., all space is classically allowed for the most energetic electrons, while in the latter, all atoms (or ions) are isolated spheres from which classical electrons could not escape. Semiclassical reasoning suggests that solids that are entirely classically allowed should be metallic, while a significant volume of forbidden regions should accompany greater inhibition of conduction, i.e., a nonzero gap. Thus we might naively expect semiconductors with covalent bonds to have classical turning surfaces.
We calculated many solids as a function of lattice parameter using the local spin density approximation (LSDA) and the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) 3 , tracking the fraction of volume that is classically forbidden. The results are shown in Fig. 1. The single most important feature is that bulk metals at equilibrium have no classically forbidden regions (CFRs). This is consistent with the original views of Kohn and Sham when developing DFT 1 , and the subsequent use of gradient expansions starting from slowly varying densities. The next most important feature is surprising: at equilibrium, covalent semiconductors also do not exhibit any CFRs. This has two important consequences. First, this shows that such systems differ fundamentally from their molecular cousins. Covalent molecules do not have forbidden regions between nuclei, but do in other directions. Covalent solids have none.
Second, this is consistent with the general behavior of density functional approximations, such as GGAs versus (global) hybrids and meta-GGAs. The exact exchange-correlation hole in most metals is short-ranged 4 , is adequately described by the LSDA, and more accurately described by GGAs 5 . In some cases (e.g., magnetic properties of transition metals), GGAs better describe metals 6,7 than meta-GGAs 8,9 and hybrids 10 , which typically do better for semiconductors [11][12][13] and insulators 12,13 presumably because metals exhibit perfect screening, whereas semiconductors and insulators do not.
But the heuristic suggested by semiclasscial physics is also not entirely without merit. If we expand the volume of a semiconductor by 40%, we found, in every case but one, the appearance of a CFR. So semiclassical reasoning is sound, but not quite quantitative. The one exception is diamond where, according to our bonding characterization, the atomic densities are strongly overlapped.
These observations are broad strokes, based on the simple solids shown in Fig. 1. A theoretical background immediately follows this overview. Our analysis begins with a generalization of the chemical bonding analysis to solids. We also examine several other interesting effects, often by studying one specific case. We show that insertion of a defect in Pt creates a CFR. An exploration of the fraction of CFR volume created as a function of lattice expansion reveals strong patterns running down columns of the periodic table. Changes of crystalline phase in Al as well as magnetic and crystalline phases in NiO are explored. As a prototype of complexity, NiO does not fit neatly into the picture of Fig. 1, regardless of phase. The Supplementary Information contains additional figures and data.
The most basic properties of an ordered solid are its conductivity [14][15][16] and bonding 17 . The standard quantum mechanical theory of conduction for ordered solids is that of Bloch bands, with insulators having filled bands below finite gaps in the eigenvalue spectrum 17 . The local bonding within a solid is often determined qualitatively from the real-space charge density.
Consider a classical electron of energy ε moving in a oneelectron effective potential v eff (r). If ε > v eff (r) everywhere, the classical electron will traverse the solid uninhibited, and the solid should be metallic. But, if the only classically allowed regions (CARs) are disjoint regions bound to atoms, the solid should be strongly insulating. Unlike the classical electron, a quantum electron can tunnel into the CFR, permitting small conductivity.
This work analyzes turning surface analogs in the framework of KS DFT 1 , an exact-in-principle quantum theory of many electrons. In KS DFT, the exact ground-state density and total energy are given by an auxiliary system of non-interacting electrons subject to an effective scalar potential v s (r). For consistency with previous work 2 , we define the turning surface of the KS potential as the set of points r c such that v s ðr c Þ ¼ ε HO ; with ε HO the highest occupied orbital energy (or Fermi energy ε F ). If ε HO > v s (r) everywhere, there is no turning surface. (One could also define a turning surface using the chemical potential μ ≥ ε HO 18 , with equality only for metals, but Eq. (1) is used in all previous work.) Highly accurate approximations of the KS potential 19,20 are computationally tractable in atoms and molecules, but more difficult in solids. The existing literature on KS turning surfaces considers finite systems exclusively 2,[21][22][23][24] . Replacing the core regions by a pseudopotential, and explicitly excluding the core regions from the analysis of the potential (as is done here), may countervail the lack of core-region structure in approximate KS potentials. Pronounced structures in the exact KS potential may feature in the low-density interstices of stretched solids, well inside CFRs.
Solids, as opposed to atoms and molecules, offer the possibility of metallic bonds, no classical turning surface (as suggested for metals by the cartoon of ref. 25 ), a sharp distinction between zerotemperature electrical conduction and insulation, and physical spontaneous symmetry-breaking (as in transition metal monoxides). Analysis of the spin-sublattices requires a further definition of the sublattice-, or spin-, CFRs v s;σ ðrÞ ¼ ε σ;HO ; σ ¼"; # : The presence of a classical turning surface, as at a monovacancy in a metal, is responsible for strong Friedel oscillations of the electron density 26 . Moreover, quantum oscillations (shell structure) in the density are observed in the CAR, while the density decays exponentially within the CFR 26 . This work presents calculations of turning surfaces for many simple solids at the LSDA and GGA levels of exchange-correlation approximations. Both usually yield close approximations to more precise KS potentials in molecules (as both KS potential and ε HO are typically too shallow by about the same amount). However, the bandgap of the exact KS potential does not match the fundamental or physical gap 18 , and typically underestimates it. The bandgaps of semilocal approximations like LSDA or GGAs are typically close to the exact KS bandgap [27][28][29] , and thus underestimate the fundamental bandgap. Hybrid functionals and meta-GGAs yield larger bandgaps when treated in a generalized KS scheme 30 . To determine the turning surface of a hybrid functional or meta-GGA, one would need to replace 20 the non-multiplicative potential operator of generalized KS theory with a multiplicative one.
Semilocal approximations to the KS potential may deteriorate under extreme expansive strains. The exact KS potential need not be analytic nor smooth, as seen in the intershell regions of atoms 31 and stretched molecular bonds 18,27,32 . However, even if a negative hydrostatic pressure could be achieved, extreme expansion of the lattice is unphysical: the work needed to stretch the lattice will eventually exceed the surface formation energy, signaling a transition to isolated clusters.

Characterizing bonding in solids
It is often useful to identify the general chemical properties of materials. Crystals are generally classified as metallic (conducting), covalently bonded, ionically bonded, hydrogen bonded, or van der Waals bonded (the latter four being insulating). However, it is difficult to quantifiably distinguish between these classifications in terms of the KS potential or density alone. The design of meta-GGAs has led to an "iso-orbital indicator" α(r) that describes the local chemical environment quantitatively 9,33 . α(r) → 0 characterizes one-and two-electron regions, α(r) → 1 uniform densities, and α(r) ≫ 1 weak bonds.
It was shown in ref. 2 that the bond type roughly determines the ratio with s AB the nearest-neighbor separation of ion sites A and B in a solid or molecule and R A and R B are their corresponding classical turning surface radii calculated for the isolated species 2 . As seen in Table 1, the value of β AB is correlated with the bond-type of a crystal. Note that, while R A and R B are spherical, s AB is determined by the lattice geometry. β AB ⪅ 0.8 generally indicates a solid with dense covalent bonding or metallic bonding, and a lack of a CFR at equilibrium. 0.8 ⪅ β AB ⪅ 1.3 indicates an ionic crystal, and 1.3 ⪅ β AB indicates a van der Waals crystal, with a CFR at equilibrium in both cases. Note that this classification differs from that in ref. 2 only in the upper bound for covalent bonding, which was 0.7 in their work, and the lower ionic bound, which they took to be 0.9. Equation (3) might be used to guess an initial geometry for a novel material of given composition.
We note a few outliers: Sn and Pb as calculated in the diamond structure, which have β AB < 0.8 and a CFR at equilibrium; and NiO and Be, which have 0.8 < β AB < 1.3, but no CFRs at equilibrium. The ground-state structures of Sn and Pb are not diamond cubic, thus the presented lattices are not true equilibrium phases. As shown later, lattice strains also induce CFRs. NiO and Be are discussed further in later sections.
Graphite and MoS 2 feature strong covalent in-plane bonds and weak out-of plane bonds. The strongly different LSDA and PBE V c / V eq values demonstrate the limitations of semilocal functionals in describing weak bonds. We have not considered hydrogen bonding-a type of weak bonding reserved for solid state hydrates-as there are numerous phases of ice that would likely exhibit a range of bonding 34 . Table 1 shows that the turning surface of a solid at its equilibrium geometry reflects its chemical bonding. CFRs arise between neighboring atoms when the atomic CARs are not overlapped (β AB > 1), making the local electron density relatively low. Metallic bonds, and dense networks of covalent bonds, lack turning surfaces and cannot be distinguished in Table 1 (although the valence electron densities of covalent solids tend to be much more inhomogeneous than those of simple metals; see Figs. 6-6 and 6-7 of ref. 35 for the densities of Si and Al). Ionic and van der Waals bonds have classical turning surfaces whose shapes could distinguish them, just as for molecules 2 .

Detailed studies
We define and report V UC as the volume per atom of the unit cell (UC), which is independent of the choice of UC (primitive, conventional, etc.). Similarly, we define the CFR fractional volume as the portion of the UC volume, which is classically forbidden. The "Methods" section describes how these quantities may be extracted from standard plane-wave DFT codes used here, the Vienna ab initio Simulation Package (VASP) and Castep.
As seen in Tables 1 and 2, the studied defect-free metals lack CFRs at equilibrium. Table 1 and Fig. 2 show that substantial expansive strains are needed to induce CFRs in metals. The fitting method and fit parameters are described in the Supplementary Information. From Tables 1 and 2 and Supplementary Table II, V c / V eq ≳ 3 for metals, except for Pt where it is 1.7.
Note also that the LSDA and PBE curves in Fig. 2 and Supplementary Fig. 6 for Al, Cu, and NaCl cross, whereas those for elemental insulators do not. For the elemental insulators, the difference between the LSDA and PBE curves is always of the same sign.
The equilibrium CFRs of insulators, as in Tables 1 and 2, depend upon the degree of insulation and the approximation used. Thus not all covalent solids have CFRs at equilibrium. Figure 3 plots the turning surface in Si at 30% expansion of the lattice parameter. In the figure, both the CAR and CFR are simultaneously connected. Supplementary Fig. 5 shows how CARs disconnect under strong expansion. Crystalline NaCl, just like its molecular form 2 , also has large PBE and LSDA CFRs. Because NaCl is a prototypical ionic solid, ionic crystals and more weakly bound crystals will likely exhibit CFRs at equilibrium.
Weakly interacting and van der Waals solids, like graphite and Ne, have PBE CFRs at equilibrium. The small (1%) PBE CFR volume in graphite (hexagonal C) at its experimental lattice constants reflects the semimetallic nature of this material. The PBE CFR in graphite lies between monolayers, as one might expect for few-  ). An asterisk indicates a Castep calculation, otherwise the result is from VASP. V c is the critical volume per atom at which a CFR appears. For those solids where fitting was performed, a numeric value of V c /V eq is given; if no fit was performed, the appropriate inequality is given. In general, V c /V eq > 1 indicates that a solid has no CFR at equilibrium, and V c /V eq < 1 indicates that it does. "I.L." corresponds to an in-layer bond, and "O.L." to an out-of-layer bond. R A and R B are the turning surface radii of the corresponding atoms or ions as calculated in ref. 2  Barring diamond C, which has much stronger covalent bonding than the other elements in the Carbon group, all narrow gap insulators studied here have V c /V eq ≲ 1.4. As semilocal functionals provide lower bounds to the bandgap-and in the case of Ge, find zero bandgap-a semilocal calculation could recognize narrow gap insulators by two properties: a nonzero bandgap and/or the presence of a CFR when the UC is expanded within 40% of the equilibrium volume. As most metals have V c /V eq ≳ 1.7 (Pt is the edge case with V c /V eq = 1.71), the expanded criterion could also identify spin-unpolarized and antiferromagnetic (AFM) NiO as non-metallic.
The classical radius of the free Ne atom is 0.87 Å, in both PBE (Ospadov et Table 2. An Ne atom in solid Ne at the equilibrium lattice constant is very similar to a free Ne atom. We can compress the Ne lattice until the CFR vanishes, as seen in Fig. 2. The Ne CFR is predicted to vanish at 0.62a for PBE. One might expect the bandgap to shrink as the CFR collapses, but the opposite is true. For the smallest lattice constant calculated here (2.85 Å), the band gap is roughly 18.57 eV, compared to a gap of about 11.51 eV (11.45 eV) at the PBE equilibrium (experimental) lattice constant, consistent with previous work that used PBE to study phases of Ne under pressure 37 . Intuition suggests that the Ne CFR should not be fully suppressed before the classical turning surfaces between adjacent atoms just touch, at a nearest-neighbor separation of 2(0.87) = 1.74 Å, using the result from ref. 2 . This is substantially smaller than the nearest-neighbor spacing in crystalline Ne for which the PBE CFR is wiped out, 2:81= ffiffi ffi 2 p % 2:00 Å. PBE and LSDA (parenthesized when different) values for the classically forbidden regions, equilibrium volume per atom, and lattice constants of select metals and insulators. For the first set of Pt monovacancy results, the cell volume and ion positions were relaxed; for the second set, the volume was fixed to the bulk value, and the ion positions were relaxed. Both sets of calculations used 31 ions in the supercell. For graphite, Ne, and NaCl, two sets of results are shown: the first at a relaxed PBE geometry, and the second at the experimental equilibrium geometry. The percent volume is taken with respect to the unit cell (percent volume per atom). Here "ds" refers to diamond structure, "hex" to simple hexagonal structure (with a four-point basis for graphite), "rs" to rock salt structure, and "zb" to Zincblende structure. The layered structure of MoS 2 is itself a prototype for dichalcogenide structure and is often referred to as the "MoS 2 structure," or by its polytype 2H b 62 , or by its space group P6/mmc 63 . The a and c parameters have the same meaning as in a simple hexagonal lattice, the z parameter (sometimes called 2z) is the spacing between neighboring sulfur layers. No LSDA calculation was performed for the NiO phases nor for Ge. Spinunpolarized (Unp.) and antiferromagnetic (AFM) calculations for NiO are listed; the AFM rs state is the correct ground state.
A.D. Kaplan et al. Thus, unexpectedly, the critical lattice constant in Ne makes the nearest-neighbor distance noticeably greater than twice the turning radius of the free atom.
As the lattice is compressed, two competing effects determine the bandgap: the bands widen, reducing the gap; and the center of the conduction band is shifted upwards with respect to the center of the valence band, widening the gap. (For an example, see the silicon density of states (DOS) at equilibrium and at a mild expansion in Supplementary Fig. 1.) This leads to a nontrivial (non-monotonic) dependence of the gap upon the lattice parameter. Note also that certain crystalline phases of Na have been experimentally observed to transition to an insulating state under strong compression, consistent with predictions made by PBE 38 .
The defect case A Pt supercell with a monovacancy defect harbors a small CFR, and as seen in Fig. 4, the CFR encapsulates the center of the vacancy perfectly. Relaxation of the supercell volume was performed by two ways: direct minimization of the stress tensor, and allowing ion positions to change within a fixed supercell volume.
The vacancy defect formation energy can be recast as the energy needed to create a curved surface within a solid 39 . The localization of the CFR to the vacancy region is a clear manifestation of this. Carling et al. 40 found that the LSDA is more accurate than GGAs for the Al monovacancy formation energy, in line with earlier results 41 for the jellium surface energy. They also found a very low electron density near the center of the vacancy, and large Friedel oscillations around it, consistent with a CFR near the center. Large voids and exterior surfaces would also give rise  Supplementary Table II. As Al, Cu, C, and Si have no CFR at their relaxed lattice parameters, each lattice must be stretched to introduce a CFR. Conversely, Ne and NaCl must be compressed to eliminate their CFRs; for completeness, the full NaCl curve is presented here. The curve for spin-unpolarized NiO is almost identical to that of NaCl (see also Supplementary   The CFR (purple) surrounds the defect, supporting the conjecture that the formation of a defect is accompanied by the formation of an internal curved surface. Regions within the PAW pseudopotential core radii are only included here to make the image clearer. For an analogous figure in Si, refer to Supplementary Fig. 2. b Logarithm of the density log 10 ½nðrÞ (blue) and ε HO -v s (r) (orange) plotted along the black line in a. The density decays exponentially as it crosses the border of the CFR, marked by gray vertical lines.
A.D. Kaplan et al. to extensive CFRs in any material. Their definition of the monovacancy volume used the liquid drop model of jellium 39 and will generally yield larger volumes than the corresponding CFR volumes.

The periodic trends
Here we consider elemental solids beyond those emphasized in Table 2. They are members of Groups 1 (alkali metals), 2 (alkaline earth metals), 14 (Group IV or Carbon group), and 18 (noble gases). The parameters of the fit functions can be found in Supplementary Tables IV-VII and the full strain curves in Supplementary Figs. 8-11. As seen in Fig. 1, the line V c /V eq ≈ 1.4 generally distinguishes metals and insulators.
Clear trends in the strain curves of elemental solids emerge as one goes down a column of the periodic table in Table 3 (see also  Supplementary Tables IV-VII). In Fig. 5a, we plot the strain curves as a function of V UC /V c , for elemental insulators. The noble gases all fall on one line, except for the lightest, He, while the Carbon group elements fall on another, except for the lightest, C. Each group has a unique, characteristic curve.
The alkali and alkaline-earth metals show similar but more complex behavior, as shown in Fig. 5b. The green line is for the heavier alkalis, the orange line is for the alkaline-earths. The lighter two alkalis, Li and Na, are shown in blue and share a shape distinct from the later alkalis. They follow the alkaline-earth curve closely, except for a dip around 1.4 V c . Moreover, Mg (in gray) is the odd one out of the alkaline-earths, rather than Be. For small strains, Mg behaves like all other alkaline-earths but, when greatly expanded, behaves more like an alkali.
Naturally, within a column of the periodic table, the critical CFR volume V c increases with atomic number, as shown in the Supplementary Information. Defining the volume of a free atom as V at ¼ 4πr 3 TS =3, with r TS the radius of the atom's classical turning surface 2 , then the ratios V c /V at are order 1 and seem to approach a column-dependent large-Z limit, with Z the nuclear charge (see Table 3 and Supplementary Tables IV-VII). The first ionization energies of the atoms exhibit similar behavior 42 .

Phases of Aluminum
As the volume of the UC changes, the equilibrium crystal phase may no longer be the ground-state phase, leading to a structural phase transition. Semilocal functionals can find states of broken symmetry of lower energy than symmetry-preserving states 43 . In addition to varying crystalline phases, we may also account for whether a spin-unpolarized phase or a symmetry-broken phase is the true ground state.
For simplicity, we have selected a few cubic phases, facecentered cubic (fcc), body-centered cubic (bcc), simple cubic (sc), and diamond structure cubic (ds), in spin-unpolarized Al to probe possible phase transitions. Using the stabilized jellium equation of state (SJEOS) 44 to fit the energy per formula unit as a function of the UC volume, PBE predicts an energy crossover for fcc and ds Al at a cubic lattice parameter a = 5.36 Å, as can be seen in Fig. 6. By following the CFR curve of lowest energy in Fig. 6, the change in inflection in the fcc Al curve is removed. A bandgap in the DOS does not accompany the crystalline phase transition. From the SJEOS parameters, we find the equilibrium volume per atom to be 16.62 Å 3 for the fcc phase (−3.74 eV/formula unit) and 28.27 Å 3 for the ds phase (−2.93 eV/formula unit). The critical transition pressure from the fcc to the ds phase is −0.070 eV/Å 3 (−11.25 GPa), as constructed from the common tangent of the equations of state. A negative hydrostatic pressure is not experimentally realizable.

Nickel monoxide
Transition metal monoxides are often poorly described by semilocal density functionals, particularly LSDA and GGAs. The spin-unpolarized state respects the full symmetry of the Hamiltonian, but calculations with approximate functionals often find a broken-symmetry AFM state of lower energy. Spin-symmetry breaking in a density-functional calculation is well known to Comparison of the fitted critical volumes per atom V c for the emergence of a CFR and the CAR volume of isolated atoms for select elemental solids.
where r TS is the turning surface radius of the corresponding neutral atom as reported in ref. 2 . Each group of the periodic table, separated by solid lines, seems to approach a column-dependent limit for V c /V at as the nuclear charge Z grows large. improve the binding energy curves of both H 2 and LiH 45 , avoiding spurious charge transfer in the dissociation limit. Spin-symmetry breaking is possible only in spin DFT, while total DFT can only get the right energies for such stretched systems with radical nonlocality. In solids, the subdivision of the lattice into spin-up and spin-down sublattices effectively doubles its size and opens a bandgap that would otherwise not be found in the spinunpolarized solution. The difficulty lower-level functionals have in describing these seemingly simple materials has garnered their title of "strongly correlated" materials. Many approximate theories outside of DFT, particularly those of Anderson 46,47 within the Mott-Hubbard model Hamiltonian framework, have had success in describing the properties of strongly correlated materials.
We present PBE-only calculations of both spin-unpolarized (Unp.) and AFM NiO in the rs and Zincblende (zb) structures. Above the Néel temperature, the Unp. rs solution gives qualitatively correct ground-state spin densities; however, no gap is present in the Unp. rs DOS. The AFM rs (AF2) state is the correct zero-temperature ground state 48  Reference 49 suggests that the AFM zb state with antiparallel spins in bilayers along the [001] conventional cubic direction (the AF5 configuration of ref. 48 ) is close in energy to the AF2-rs state near equilibrium. We have also investigated the AF5-zb state for completeness.
In Fig. 7, we plot the CFRs of AFM (AF2-rs and AF5-zb) and Unp.rs NiO. The spin-up and spin-down CFRs for AFM-rs and AFM-zb NiO exhibit identical behavior. Moreover, the similarity of the Unp. and AFM rs curves demonstrates that the crystalline structure plays a large role in determining the size of the CFR under strain.

DISCUSSION
No bulk metal that we studied had a CFR at equilibrium, but covalent semiconductors also lack CFRs at equilibrium. Thus the fractional CFR volume at equilibrium cannot be used directly as an indicator of conductivity. CFRs emerged in all narrow gap insulators studied here (excluding C) when expanded by 40% of the equilibrium volume. Thus a semilocal calculation can recognize narrow gap insulators by either a nonzero gap (C or Si) or an emergent CFR if V c /V eq ≲ 1.4 (Ge, where semilocal functionals often find zero bandgap). However, the existence of a CFR can be used as a theoretical tool in understanding the role of semilocality in describing solids.
Since standard density-gradient expansions are derived for slowly varying densities without CFRs, the absence of CFRs in metals at equilibrium, along with the short range of the exact exchange-correlation hole around an electron in a metal, suggest that the local density and its low-order derivatives suffice for an accurate approximation to the exchange-correlation energy of a metal. This argument, but based on the hole alone, was made in ref. 51 to explain why "de-orbitalization" of a meta-GGA can improve 52 the magnetic moment in solid Fe.
The lack of CFRs in narrow-gap insulators suggests that semilocal functionals may still be accurate in predicting some, but not all, of their material properties. This is evident in the accuracy of LSDA and PBE in predicting the lattice geometry of Si, diamond, and the in-plane covalent bonding of graphite. However, a degree of nonlocality is needed to predict conductivity and thermodynamic properties of insulators. Weakly bonded materials like Ne and complex materials like NiO require an even higher degree of nonlocality. A truly general purpose density functional must balance the semilocality demanded by metals and Fig. 6 Competing crystalline phases in Al and their CFRs. a CFR evolution for different spin-unpolarized cubic crystalline phases of Al, as calculated with PBE in VASP. b The fcc and ds phases become degenerate near a = 5.36 Å The fcc, bcc, and sc phases become nearly degenerate in energy as a grows large. Note that absolute energies do not have physical meaning in a pseudopotential calculation, only differences in energies. The fcc, bcc and sc phases were taken to contain one formula unit in the computational cell, and the ds phase to contain two formula units in the computational cell. the nonlocality needed to describe insulators, atoms, and molecules.
Just as in a small molecule, the classical turning surface at equilibrium in a periodic solid reflects the bonding type. There is no classical turning surface for typical metallic bonds and for dense networks of covalent bonds, but classical turning surfaces of characteristic shapes appear for ionic, hydrogen, and van der Waals bonds. Equation (3), with β AB determined by the bond type, roughly predicts the equilibrium bond length between neighboring atoms A and B, and might be used to find an initial guess for the geometry of a novel material. β AB falls in the same range for metallic as for covalent bonds: both show a strong overlap of the bonded atoms. CFRs in solids arise where neighboring atoms are weakly overlapped, making the local electron density relatively low; CFRs are absent when the atoms are strongly overlapped. β AB ⪅ 1 seems to imply no CFR for solids in equilibrium with only one bond type, as shown in Table 1.
As a predictor of electrical insulation, the existence of a gap between occupied and unoccupied states in the LSDA or GGA band structure is much more reliable than the existence of a CFR. A CFR is missing not only in metals but also in many covalently bonded semiconductors and even in some AFM insulators like NiO. The CFR in NiO is seemingly insensitive to the magnetic phase. For the solids studied here, our calculations found no CFRs for metals, large CFRs for wide-gap insulators, and the emergence of CFRs when small-gap semiconductors are mildly expanded.
A monovacancy in a metal can induce a CFR, and an expansive strain in any material can induce a CFR or increase its volume. In both cases, the CFR emerges in a region of relatively large density depletion, such as the interstice of stretched Si. In any solid, the fractional CFR volume (the volume of the CFR relative to a chosen cell volume) is sensitive to the crystalline phase, nearest-neighbor separation, and magnetic phase. Those wider-gap insulators with a CFR at equilibrium can be compressed until the CFR vanishes. Layered materials may have a CFR at equilibrium, at least when a density functional approximation overestimates the interlayer spacing sufficiently, as PBE does for MoS 2 . Our analysis of strain has been limited to homogeneous strains of cubic lattices. Applying unequal expansive strains along the high symmetry directions of the lattice would likely yield quite different CFR curves than those shown here. In materials with structures more complicated than cubic, unequal strains would likely reflect properties of the bonding along the direction of the strain, in the vein of our bonding analysis for layered materials. However, as the possibilities for this are myriad, we defer this to future study.
CFRs are also characteristic of perfect ionic and molecular crystals at equilibrium. Our analysis supports the conclusion that rare gas atoms in the crystalline phase are nearly free. Ionic crystals can have large CFRs at equilibrium (as in NaCl but not NiO). We showed that graphite and MoS 2 , where intermediaterange van der Waals interactions dominate between monolayers, have CFRs located solely between monolayers and that their corresponding monolayers have no in-plane CFR. Our work demonstrates that weakly bound solids tend to have prominent CFRs. Hydrogen-bonded crystals like ice, while not tested here, can be expected to have substantial CFR volume fractions, as suggested by Fig. 8 for the water dimer in ref. 2 .

Computational details
All calculations were performed with either the VASP 53-56 , or the Castep code 57,58 , or both. All GGA calculations used the PBE GGA 3 , and all LSDA calculations used the Perdew-Zunger parameterization of the uniform electron gas correlation energy 59 . The calculations in VASP were performed with a cutoff energy of 800 eV, a Γ-centered mesh of spacing 0.08 Å −1 , energy convergence of 10 −6 eV, and stress convergence at 10 −3 eV/Å. To determine equilibrium geometries in VASP, for metals, first-order Methfessel-Paxton smearing with parameter of 0.2 was used, and for insulators, the Blöchl tetrahedron method was used. VASP's internal methods were used to determine the relaxed cell volume.
All calculations were spin-unpolarized, except those for AFM NiO. For AFM NiO, the relaxed volume was determined by fitting to the SJEOS 44 . Convergence was aided using a linear magnetization-density mixing scheme (AMIX = 0.2, AMIX_MAG = 0.8, BMIX = BMIX_MAG = 0.0001). From ref. 48 , the ground state of AFM NiO has antiparallel spins along the [111] direction of the conventional cubic cell (AF2 ordering, a four-ion basis), which was used here. Magnetic initialization of ±2 μ B was given for the Ni atoms and 0 μ B for the O atoms.
Approximate bandgaps were calculated from the DOS, on the default resolution for VASP (300 total sampling points). Increasing the resolution of the DOS to 3000 sampling points gave changes at most of 0.05 eV only near equilibrium for Si. For Ge, no change in the DOS gap could be discerned.
In Castep, a density-mixing algorithm was used to reach self-consistency, and geometries were determined with a BFGS (Broyden-Fletcher-Goldfarb-Shanno) energy minimization scheme with the finite basis set corrected for stress 60 . After relaxation, a calculation at the equilibrium volume using the Blöchl tetrahedron method was performed to accurately determine the DOS. Accurate 61 PAW on-the-fly pseudopotentials were used throughout. Supplementary Tables VIII-LIII  present all raw data. For monolayers, a 45 × 45 × 1 k-point grid was used in conjunction with the Blöchl tetrahedron method. All other parameters remain the same from bulk calculations. The c direction was padded with 30 Å of vacuum region to reduce interactions between image monolayers.
In density functional plane-wave codes, the densities and potentials are stored on a uniform grid R, the dimensions of which are determined by the size of the UC and the plane-wave cutoff energy. Acceptable convergence of the total energy relies on suitable convergence of the potentials and densities on this grid. The values of v s (R) are obtained from this grid. In core regions, the true potential is much deeper than the pseudopotential, so these are classically allowed. Thus the PAW pseudopotential core regions were excluded from the CFR. (Frozen-core pseudopotentials were used in both VASP and Castep.) The self-consistent electronic eigenstates give ε HO (the Fermi energy ε F in a metal), and the regions where ε HO − v s (R) < 0 define the CFR. We assign equal volume to each point relative to the UC, as the real-space mesh is uniform. To find the volume of a CFR relative to a given cell volume, we need to average over a UC or repeat-unit of the periodic crystal. But such a UC must in some cases contain more than one atom. For ease of comparison, we define and report V UC as the volume per atom of the UC, which for a given material is of course independent of the choice of UC (primitive, conventional, etc.). Similarly, we define N UC as the number of mesh points per atom of the UC. Then the volume of any mesh point is V UC /N UC . If there are N CFR points at which ε HO − v s (r) < 0, the volume per atom of the CFR is The dimensionless "fractional volume" of the CFR, which was used throughout, is defined as the number of real-space mesh points within the CFR relative to the total number of mesh points in the UC.
As the fractional CFR volume v → 0, our method requires ever finer realand reciprocal-space meshes to resolve v. This need is limited by the resolution determined by the plane-wave cutoff energy. Our data for v ≪ 1 will necessarily be more noisy than for larger values of v. Despite this, we show a posteriori that reasonable fits to v(V UC ) may be found.
Each code uses differently generated pseudopotentials with different optimal basis set cutoff energies (and hence pseudopotential grid sizes, etc.), different energy minimization schemes, and different Brillouin zone integration methods. To ensure that our method is not dependent upon the numerical methods of a particular code, we have verified that the Castep and VASP results are consistent.

DATA AVAILABILITY
All data supporting the findings of this work are available in the article and its Supplementary Information. Extra data and machine readable data are publicly available at the Materials Cloud Archive, 10.24435/materialscloud:2h-zq 68 .