Toward automated screening of band gap sensitivity in 2D materials

Computational materials science relies on simple, yet efficient, measures and indicators of the modeled materials’ properties. Ideally, the desired properties should be linked to such scalar quantities that can be obtained in polynomial time and efficiently integrated within automated high-throughput screening loops for screening and sorting out the evaluated materials to the desired categories. Here, we focus on the freestanding gapped 2D materials and scalar indicator of their band gap sensitivity to the presence of additional stacked 2D layer/s. The proposed measure uses only a freestanding model of a given material, and it is based on an automated integration of the electron density of frontier orbitals extending into the vacuum within the model unit cell. The usefulness and limitations of such an approach for materials pre-screening are demonstrated on a handful of 2D materials, like, e.g. MXenes, graphane, fluorographene, or, allotropes of phosphorus.


Motivation
The success of computational materials science relies on the availability of simple, yet efficient, measures and indicators of the desired materials properties. Currently, density functional theory (DFT) combined with high-throughput screening (HTS) of large sets of materials has been used as a workhorse to sort out materials with the desired properties [1][2][3][4]. An automated HTS of the desired properties is limited by the availability of scalar measures of the inherent materials properties that could be embodied within automated procedures used for sorting out the individual materials models into the desired categories. An example of such an approach, within a given model (i.e. computational setup), is sorting out a set of materials into metals or semiconductors categories, based on a (scalar) band gap value. However, more complex materials-sorting requirements may emerge, and it is not always trivial to propose a single-valued measure of the desired properties that would enable direct automated materials classification.
Zhou et al [5] reported on a DFT based computational study of interlayer-decoupled Sc-Based MXene 2D material, Sc 2 CCl 2 , showing weak van der Waals (vdW) interactions and high carrier mobility. The interlayer decoupling of the optoelectronic properties between the individual 2D material monolayers, an important and unique property, was attributed to the electronegativity of Cl atoms and weak wave function overlap at the Fermi level. Unlike 2D materials such as phosphorene [6][7][8][9][10], showing a strong dependency of the properties with respect to the number of layers, interlayer-decoupled 2D materials properties are expected to remain robust with respect to external factors, like the presence of additional layers. An interesting question arises, how to computationally evaluate such a property on a set of 2D materials models in an automated fashion.
This work focuses on computational modeling of the interlayer (de)coupling of optoelectronic properties of 2D materials [11][12][13][14][15][16][17][18] from the HTS point of view. As a first step toward this goal, we consider (de)coupling of the band gap as one of the most important materials optoelectronic properties. Our goal is to provide a scalar measure of the band gap sensitivity (i.e. change in the band gap value caused by external factors like, e.g. presence of additional layers) of a given freestanding 2D material model to external factors that could be evaluated in a manner suitable for automated HTS approaches. To this end, we devise a measure of the band gap sensitivity based on the part of the normalized and unscaled valence band maximum (VBM) and conduction band minimum (CBM) Kohn-Sham orbitals extending into the vacuum region of the simulation cell (see equations (1) and (2) below). This choice is based on the following assumption: the more the orbitals near the Fermi level extend from the 'volume' of the incident 2D material, the more sensitive the properties, including the band gap, to the external factors, like, e.g. the presence of additional layers. This limits the proposed measure to the super-systems without charge transfer, i.e. type I/II heterojunctions [19]. The performance of the proposed indicator is evaluated on a handful of 2D materials, including Sc-based MXenes Sc 2 CT 2 (with terminal groups T = Cl, OH, or F, figure 1) [5,[20][21][22][23][24], graphane (CH) [25][26][27], fluorographene (CF) [27][28][29], black phosphorus (P) [6,8,30,31], blue P [32][33][34], MoS 2 [11,35,36], and, BN [37][38][39][40].
The proposed approach provides interesting insights into the electronic structure of 2D materials, and, it enables fully automated qualitative pre-screening of vdW materials for high sensitivity (useful, e.g. for sensors), or, low sensitivity (useful for pointing out interlayer decoupling) of the band gap to the external perturbations, using solely a freestanding 2D-material monolayer DFT slab model [12].

Models and methods
The slab models of all 2D systems were obtained from the C2DB database [41] of 2D materials, except for the black P and blue P models, obtained from Guan et al [33]. The bulk counterparts were fully optimized as trivially AA stacked systems. Atomic positions, as well as unit cell parameters, were optimized with PBE [42] xc functional and D3 dispersion correction [43], using the Quantum Espresso [44][45][46] software package. Single-point calculations with PBE [42] and HSE [47] xc functionals used the Brillouin zone sampling with 12 × 12 × 1 k-points for 2D systems, and, 12 × 12 × 4 for bulk. The nuclei were represented by the norm-conserving correlation-consistent effective core potentials (ccECPs) [48][49][50][51][52], except for black P, where ultrasoft pseudopotentials (USPP) [53] were used (for convergence problems) for HSE calculations (tests performed on other systems confirmed that the results of USPP and ccECP are well consistent to below 0.1 eV on band gaps). The kinetic energy cutoff was set to 400 Ry for both PBE as well as HSE, xc functionals (except for USPP runs using 50 Ry).
The measure of the band gap sensitivity κ is based on the contributions from the integration of the real-space-grid-projected Kohn-Sham frontier orbitals, VBM and CBM, in the unit cell of the evaluated 2D material slab model, The prime indicates that the z-integration is limited to two sub-spaces within the slab unit cell, so that it excludes the volume of the material (figure 2): 0 ⩽ z ⩽ z a and z b ⩽ z ⩽ c, z a = z min (a) − qr vdW (a) and where z min is a z-coordinate of the bottom-most surface atom in the simulation cell a, z max is the same for the top-most surface atom b, r vdW (j) is a vdW radius of an atom j, c is the size of the simulation cell along z-direction (containing vacuum region), q is a technical parameter of the method, and, ϕ i is a frontier orbital i (CBM or VBM). Note that all the considered systems possess the same surface on both (top/bottom) sides (e.g. figure 1). In non-symmetric cases, it is necessary to integrate each of the half-space separately, so that disctinct κ values would be obtained for respective surfaces.
The role of the empirical q parameter (in q · r vdW term), is to set an integration boundaries for the z-coordinate on the side of the surface (figure 2), ensuring that only the tails of the valence part of the terminal atom(s) density is included, while avoiding integration of the (residual) core density. This is particularly important in monolayer systems like boron nitride (BN), where integration would otherwise encompass the full volume of the slab model cell, resulting in a loss of resolution for the indicator (κ = 1). The results primarily rely on tails of the frontier orbitals extending into the vacuum, and, from the point of view of qualitative screening and sorting of materials, they are not overly sensitive to this choice (table 1). Clearly, for all values, κ enables sorting of materials as desired, i.e. orders of magnitude larger values are observed for exterior-sensitive systems and vice versa. The integration starting point can be chosen almost freely (e.g. atomic or covalent radii may be used as convenient choices). Based on multiple tests (using r vdW ), we pick κ ≡ κ( 1 2 ) as the most convenient value. Finally, let us comment on the choice of the far-from-surface z-coordinate integration boundaries. Since κ is based on the integration of the one-particle orbitals, and the orbital tails of the terminal atoms diminish   exponentially into the vacuum, we expect no qualitative differences in its value, if sufficiently wide vacuum region is used (sizes typical as in routine slab computations). The tests with an enlarged vacuum region in an exterior-sensitive Sc 2 C(OH) 2 system, based on modification of c lattice parameter, indeed confirm that κ converges rapidly ( figure 3). The physical band gap sensitivity δ∆ = |∆ 2D − ∆ 3D | of a given 2D material was evaluated using a difference between the band gaps of the 2D slab system model (∆ 2D ) at the k-point of interest, i.e. relevant to the selected 2D material, and its bulk counterpart (at the same k-point) (∆ 3D ) estimated using the AA stacking model.

Results and discussion
First, let us focus on comparison of Sc 2 CCl 2 and Sc 2 C(OH) 2 MXenes ( figure 1(a) and (b)). The results (table 2) confirm that 2D Sc 2 CCl 2 is an 'interlayer-decoupled' system, as long as its electronic band gap value is resistant (at both, PBE and HSE levels) with respect to the presence of additional layers by the small band gap change when going from 2D to 3D system (δ∆ = 0.05 eV/0.05 eV, at PBE/HSE level), in line with the existing literature [5]. In contrast, Sc 2 C(OH) 2 ( figure 1(b)) turns out to be a system with a  Table 2. Direct band gaps (eV) of a freestanding 2D monolayer (∆2D) and its bulk counterpart with trivial stacking (∆3D) at the same k-point, estimate of the physical band gap sensitivity δ∆ (eV), relative 2D→3D band gap change δ∆R = δ∆/∆2D, and, measure of the band gap sensitivity κ (equations (1) and (2)), together with the individual contributions (equation (1)). pronounced band gap sensitivity (again, at both, PBE and HSE levels) to the presence of additional layers (δ∆ = 1.55-2 eV). These results lead to the following question. What is the governing difference between these two MXene systems, so that their band gap sensitivity is so much different? It is not difficult to find an answer. Namely, the spatial properties of the frontier orbitals ( figure 1, VBM, CBM) show a striking difference that appears to be responsible for the difference. While in Sc 2 CCl 2 , both VBM and CBM orbitals reside mostly within the 2D material 'volume' (figure 1(a)), in Sc 2 C(OH) 2 , one of the orbitals (CBM) significantly extends into the vacuum region ( figure 1(b)). Clearly, in the case of Sc 2 CCl 2 , external 'objects' are not allowed to directly interact (by Pauli repulsion or hybridization) with the (protected) orbitals that govern ∆ 2D . In the case of Sc 2 C(OH) 2 , interactions between the 2D system and the approaching system/s necessarily deform an outward-facing part of the original CBM orbital and the band gap (of the original 2D system, at the same k-point) changes accordingly.

System k-point
From the point of materials HTS, looking for a suitable scalar parameter that would enable an automated classification of the materials, the difference between the frontier orbitals of the discussed 2D materials is numerically best represented by the normalized and unscaled part of the Kohn-Sham orbitals that reaches out from the bulk to the vacuum. The proposed parameter of the band gap sensitivity κ is nothing but the same, i.e. it is based on the integration of the parts of the frontier orbitals that extend into the vacuum (equation (2)). Comparison of κ (table 2) for Sc 2 CCl 2 and Sc 2 C(OH) 2 , leads (for both xc functionals) to the findings in line with those from comparison of δ∆. Namely, one finds that the 'insensitive' Sc 2 CCl 2 2D material shows a small value of κ (∼ 0.0002), whereas the 2D material sensitive to the presence of additional layers, Sc 2 C(OH) 2 , shows orders of magnitude larger value (κ > 0.3). Therefore, instead of δ∆ involving the 3D counterpart of the analyzed 2D model, κ (and possibly also functions based on it) provides a convenient measure of the electronic band gap sensitivity for a given 2D material, which could be readily used within automated materials classification HTS procedures. In the following, we explore the performance and limits of κ, as a quantity possibly useful for the prediction of the band-gap-(in)sensitive systems, on a handful of additional 2D systems.
The next MXene of our interest is Sc 2 CF 2 ( figure 1(c)), another halogen-based analog of Sc 2 CCl 2 . Based on its κ ∼ 1 × 10 −5 , one would attribute this system into the class of systems insensitive to the external factors, like Sc 2 CCl 2 . In fact, this is what one finds when the related ∆ 2D and ∆ 3D are compared with each other. At both considered (PBE and HSE) levels of theory, the band gap of Sc 2 CF 2 remains almost intact (δ∆ < 0.05 eV). Again, κ seems to provide a useful measure of the 2D-system band gap insensitivity to the external perturbations.
Regarding the carbonaceous gapped 2D materials, C 2 H 2 (graphane) and C 2 F 2 (fluorographene) were selected as the most prominent members of this important group of 2D materials. C 2 H 2 shows by far the largest change of the band gap, from the considered systems, when moving from the 2D system to bulk (δ∆ > 4 eV) and adequately large κ. That is, κ would correctly indicate that this material is overly sensitive to its surroundings. In addition, the band gap of C 2 H 2 is found to be more sensitive than in Sc 2 C(OH) 2 , and, κ provides the same ordering of sensitivity, i.e. δ∆(C 2 H 2 )> δ∆(Sc 2 C(OH) 2 ), as well as, κ(C 2 H 2 ) > κ(Sc 2 C(OH) 2 ). On the other hand, 2D C 2 F 2 monolayer turns out as a material with a band gap that is much less sensitive (δ∆ ∼ 0.5-0.6 eV) than in the case of C 2 H 2 (δ∆ > 4 eV, table 2), however, much more sensitive than the halogen-based MXenes like, e.g. Sc 2 CF 2 (δ∆ < 0.05 eV). At first sight, this may look surprising, because both systems are F-terminated and their surface dipoles are oriented in ways that tend to stabilize the electron density within the 2D monolayer (instead of negatively-charged surface). The key difference between the C 2 F 2 and Sc 2 CF 2 systems, from the point of view of the band gap sensitivity, is a surface F-atom coordination and the related frontier-orbital arrangements, that are correctly reflected in κ and its contributions (table 2). While in C 2 F 2 at least one of the frontier orbitals extends to the vacuum [54], in Sc 2 CF 2 , this is not the case ( figure 1(c)). Finally, δ∆ R value (at PBE level), that can be interpreted as a percentage change of the band gap when transitioning from 2D to 3D case, is more than 10 times larger for C 2 F 2 system, than the same for the Sc 2 CF 2 system. This represents 'just' a ∼17% change, compared to the corresponding ∼118% change in the C 2 H 2 system. For a well studied MoS 2 system, κ indicates a weak band gap sensitivity estimate, in line with δ∆ of the order of 0.1 eV (table 2).
Black P is known for its strongly thickness dependent properties [6][7][8][9][10], including the electronic band gap, contrary to the blue P allotrope [55]. The results (table 2) show that κ correctly predicts a difference between the band gap sensitivities of 2D monolayers of black P, that is consistent with the related δ∆.
2D hexagonal BN monolayer is an example of a system where the sensitivity of the band gap (δ∆ ∼ 0.5 eV) falls in between the 'insensitive' and 'overly sensitive' examples. The same is true for κ, i.e. it is consistent with order-of-magnitude expectations based on previous cases. The BN and blue P systems teaches us that there is a 'gray zone' and sensible materials screening thresholds based on κ must necessarily lead to either, false-positive or false-negative materials searches, as expected for simplified scalar quantities trying to describe such complex systems.

Conclusions
A promising measure for determining the sensitivity of the band gap in 2D materials to the presence of additional monolayers has been presented, based on the spatial properties of the freestanding 2D material frontier orbitals extending to the vacuum. Representative gapped 2D systems were considered for testing, including MXenes, carbonaceous materials, allotropes of phosphorus, MoS 2 or BN. In most of the studied systems, κ offered valuable insights and band gap sensitivity estimates that aligned with expectations. κ may thus find applications in band gap engineering applications [55] or HTS of 2D materials with (in)sensitive optical properties with respect to (general) external factors. HTS sorting procedures based on κ should take into account potential false-positive or false-negative results during the screening procedure design process. For screening of the band-gap-decoupled systems, threshold of κ < 0.005 seems to provide a reasonable choice. Estimates of κ based on cheap PBE xc DFT functional are very promising, since they show similar accuracy as the same based on a much more demanding HSE xc functional.

Data availability statement
All data that support the findings of this study are included within the article.