Adding HF to the MP2/SLATM-L2 pair produces a marginal gain for the latter methods, whereas HF presents a strong gain in reliability: the final *Q*95 values are 2.5 (HF), 1.8 (MP2) and 1.7 kcal/mol (SLATM-L2). However, this comes at the price of a large number of system rejections: for ∼2.0 kcal/mol, only 1/4th of the 6211 systems are selected by their similarity. For comparison, this number is about 2/3rd for the MP2/SLATM-L2 comparison.

In this context of uncorrelated error sets with different accuracy levels, one sees that similarity selection has a notable positive impact on the reliability of predictions by any of the methods, even the most accurate ones. It is striking that MP2 or SLATM-L2 might benefit from comparison with HF, but, as already discussed for band gaps (Figure 9) , this proceeds mainly by elimination of systems with large errors.

#### **4. Conclusions**

We asked whether picking only results that are similar to different methods would improve the accuracy of their predictions (in spite of possibly eliminating a significant part of the calculations done). The use of probabilities to treat reliability and similarity was illustrated on two benchmark data sets, one of band gap calculations with different density functional approximations, the other of effective atomization energies with two *ab initio* methods and one machine-learning method.

For the properties and methods studied, the thresholds for reliability and similarity were chosen quite generously. For the band gap data set, we found that similarity of the density functional results had only a marginal impact on improving the prediction accuracy. This is consistent with previous findings that the differences between density functional approximations are less important when considering the error distributions [1,2], or taking into account experimental uncertainty [19].

For the effective atomization energies data set, in which the error sets are uncorrelated, notable improvements of reliability after similarity selection were observed for all methods, even the most accurate ones. Roughly, we observed two categories of results:

Note that the size of the data sets might have an impact on the uncertainty of all the statistics. For the smaller datasets, this uncertainty might be comparable with the observed differences between statistics. Bootstrapping approaches, such as the ones used in our previous works [1,3], could be used to this effect. This was not the focus of the present study, and uncertainty management will be considered in forthcoming research.

**Author Contributions:** Conceptualization, A.S. and P.P.; methodology, A.S. and P.P.; software, A.S. and P.P.; validation, A.S. and P.P.; formal analysis, A.S. and P.P.; investigation, A.S. and P.P.; resources, A.S. and P.P.; data curation, A.S. and P.P.; writing—original draft preparation, A.S. and P.P.; writing review and editing, A.S. and P.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data is contained within the article.

**Acknowledgments:** This paper is dedicated to Karlheinz Schwarz, who not only was a pioneer in the development of density functional theory, but also of computer programs using it for large systems. This led to a large number of applications, enabling us to perform statistical analyses as presented in this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Accuracy and Precision in Electronic Structure Computation: Wien2k and FPLO**

**Manuel Richter 1,\*, Seo-Jin Kim 2, Klaus Koepernik 1, Helge Rosner 2 and Arnulf Möbius 3**

**Abstract:** Electronic structure calculations in the framework of density functional theory are based on complex numerical codes which are used in a multitude of applications. Frequently, existing experimental information is used as a gauge for the reliability of such codes. However, their results depend both on the chosen exchange-correlation energy functional and on the specific numerical implementation of the Kohn-Sham equations. The only way to disentangle these two items is a direct comparison of two or more electronic structure codes. Here, we address the achievable numerical accuracy and numerical precision in the total energy computation of the two all-electron density-functional codes Wien2k and FPLO. Both codes are based on almost independent numerical implementations and largely differ in the representation of the Bloch wave function. Thus, it is a highly encouraging result that the total energy data obtained with both codes agree within less than 10−6. We here relate the term numerical accuracy to the value of the total energy *E*, while the term numerical precision is related to the numerical noise of *E* as observed in total energy derivatives. We find that Wien2k achieves a slightly higher accuracy than FPLO at the price of a larger numerical effort. Further, we demonstrate that the FPLO code shows somewhat higher precision, i.e., less numerical noise in *E* than Wien2k, which is useful for the evaluation of physical properties based on derivatives of *E*.

**Keywords:** density functional theory; DFT codes; electronic structure calculation; numerical accuracy and precision

#### **1. Introduction**

Computations of electronic structure provide an important input to research in physics, chemistry, and materials science, and density functional theory (DFT) is behind a major part of such computations. While DFT is a formally exact theory, its practical use relies on the implementation of approximative exchange-correlation functionals. A second but not less important ingredient to DFT applications consists in the development and verification of numerical electronic structure codes. Here, we consider the latter subject.

The task to solve the Kohn-Sham equations of DFT is well defined, but its solution is far from being trivial. The mathematical complexity of this system of non-linear integrodifferential equations is mirrored by the complexity of its solvers, numerical codes with typically several 105 source lines. Codes of such a size and complexity bear two unavoidable problems: approximations and coding errors. Fortunately, a large number of codes has been developed in past decades by independent teams. These codes partly rely on completely different approximations for, e.g., the representation of the wave function. Thus, many of them can be considered as virtually independent implementations with little risk to suffer from the same source of inaccuracy. Hence, deviations from the numerically exact solution can be considered as code-specific and uncorrelated, at least for codes with different

**Citation:** Richter, M.; Kim, S.-J.; Koepernik, K.; Rosner, H.; Möbius, A. Accuracy and Precision in Electronic Structure Computation: Wien2k and FPLO *Computation* **2022**, *10*, 28. https://doi.org/10.3390/ computation10020028

Academic Editor: Henry Chermette

Received: 6 January 2022 Accepted: 7 February 2022 Published: 11 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

wave function representations. As noteworthy exceptions, we mention standard algebra routines that are frequently taken from common software packages and methods for *k*-space integrations. The latter are usually implemented independently in each code but may lead to method-specific convergence behavior with the number of *k*-points.

In 2016, a comparison was published regarding the volume-dependent total energy, *E*(*V*), calculated by means of 15 different DFT codes for 71 elements of the periodic table [1]. The goal of that work was to evaluate the achievable precision of Kohn-Sham solvers under well-defined common choices of exchange-correlation functional and treatment of relativistic effects. Thereby, *precision* was defined as the scatter among computed equationof-state data, in distinction to *accuracy* which would describe the deviation of the mean DFT result from experiment [1].

In the present work which is again solely focused on the comparison of computations, we will use the terms accuracy and precision in a slightly different sense. Motivated by the idea that the total energy is a variational quantity, lower total energy will be identified with higher numerical accuracy in the common meaning of a smaller systematic error [2]. This choice comes with the known caveat that numerical problems like basis set overcompleteness can reduce the resulting total energy in an unpredictable manner [3] and have to be excluded by appropriate data analysis. Furthermore, our definition only makes sense if all calculations make use of one and the same exchange-correlation approximation. With this condition fulfilled, the numerical accuracy of a modern DFT code is mainly determined by the completeness of the basis set for the representation of the Bloch states.

Regarding precision as a description of random errors [2], we will associate this term with numerical noise along the *E*(*V*) curve. This choice is motivated by the following consideration. DFT computations are usually quasi-deterministic. Thus, it is not possible to perform a statistical analysis for a single point on *E*(*V*). However, different points on *E*(*V*) are known to show quasi-random deviations from averaging approximants to *E*(*V*), such as polynomial fits. The default numerical settings of advanced codes are nowadays tight enough to keep these quasi-random contributions to *E*(*V*) well below its systematic deviations from a hypothetical, exact solution of the Kohn-Sham problem. In other words, the numerical precision of *E*(*V*) is much higher than its numerical accuracy. However, this may not be the case for derivatives of the total energy which are of interest on their own. For the sake of brevity, we will drop the qualifier "numerical" for the considered accuracy and precision in most places below.

Seven of the codes taken into account in the abovementioned comparison are socalled all-electron codes which do not rely on approximative pseudopotentials but include the complete nuclear potential and the electronic core wave functions in the calculation. Considering the scatter of results only among these seven codes used with their numerically best settings, the data produced by Wien2k [4] turned out to be close to the barycenter [1]. This finding confirms the status of Wien2k as the "gold standard" of DFT codes: a welldeserved result of several decades of dedicated work by the developer team of Wien2k, formed and headed by Karlheinz Schwarz.

The predictive computation of physical properties by means of DFT codes has been widely established in past decades. If applied to a large set of materials, such a property screening is faced with resource considerations. Disregarding so-called *O*(*N*) methods – they are not among the most commonly used DFT codes and are often geared towards certain subsets of compounds like insulators –, the required CPU time and main storage grow with the third and second power, respectively, of the Hamilton matrix rank. Hence, a code with a small matrix, i.e., a small basis set, is usually superior to a code with a large basis set in terms of computation time and main storage demand. However, small basis sets frequently limit the achievable accuracy. Thus, it is important to gauge the quality of codes with small basis sets against a confirmed quality standard such as Wien2k both in terms of accuracy and precision. In the current study, we compare Wien2k and the full-potential local-orbital (FPLO) code [5] with each other. They are based on completely different concepts to represent the Bloch wave function and the FPLO code is equipped with a much smaller basis set than Wien2k. However, both these codes belong to the

all-electron category of electronic structure codes. This allows a direct comparison of their total energies.

To establish a procedure for comparing the quality of two codes concerning their predictive power beyond total energies, we advance from the previous equation-of-state comparison [1] to the consideration of total energy derivatives which give access to more detailed thermodynamic information. For example, the first derivative of the total energy with respect to the magnetization equals the external field necessary to stabilize this magnetization, and a zero second derivative indicates a metamagnetic transition. Yet more subtle is the detection of an electronic topological transition (ETT), also called Lifshitz transition, on the sole basis of total-energy data. For hexagonal Osmium, a previous experimental study claimed visibility of an ETT as a kink in the ratio of lattice parameters *c*/*a* vs. pressure *p* [6]. Theoretical equation-of-state data should in principle allow to identify such a peculiarity. In later work, however, the existence of this particular ETT was doubted, see Ref. [7] and references therein.

Here, we demonstrate for a particular case that a conventional ETT can be detected in the third derivative of *E*(*a*) as obtained from current DFT codes, *a* denoting the lattice parameter. The reader may object that ETTs can simply be found by careful inspection of band structure and density of states. This was indeed our starting point in seeking to distinguish the traces of ETTs from numerical noise. We achieved a substantial reduction of the latter by modifying the default numerical settings for both considered codes. This allowed us to obtain even *d*3*E*/*da*3 relations with high resolution. The obtained precision enabled us to perform convincing Birch-Murnaghan fits using the minimum possible number of data points.

Our exercise focuses on the comparatively simple case of face centered cubic (fcc) Aluminium. This metal was already the target of very early Fermi surface investigations by electronic structure theory [8]. Later, Fermi surface changes under different types of pressure were predicted [9] and an ETT was observed by nuclear magnetic resonance experiments on Al under hydrostatic pressure [10].

Details of the computational methods are provided in the following Section 2. Results and the related discussion are presented in Section 3, and Section 4 summarizes this work.

#### **2. Methods**

The Wien2k calculations were performed in the version 19.1 of this code [4]. The parameter *Rk*max was varied in the range from 4 to 12, and the muffin-tin sphere radius *R*mt was chosen as 2.39 Bohr radii for Aluminium. As convergence conditions, 10−6 for density and 10−8 Rydberg for the total energy were used in almost all calculations. The only exception was one data point for *a* = 0.3646 nm, *Rk*max = 8, for which only a density convergence of 10−5 could be achieved.

FPLO was used in the version FPLO-18.00-52 [11]. Self-consistency of the iterations was considered to be achieved upon stability of both density and total energy. The default values, 10−6 for density and 10−8 Hartree for energy, were taken as related thresholds if not indicated otherwise. A specific modification of the numerical mesh for orbital-potential matrix elements is explained and discussed in Section 3.3.

All calculations, performed with either of the codes, were conducted in a scalar relativistic mode according to Koelling and Harmon [12] and using the generalized gradient approximation (GGA) in the parameterization by Perdew, Burke, and Ernzerhof (PBE) [13]. Integrations in reciprocal space were done with a linear tetrahedron method including Blöchl corrections [14]. It should be noted that broadening techniques are less suited for the detection of ETTs.

The *k*-point numbers of the specific calculations are defined in Section 3 using the following notation: A mesh with *n* × *n* × *n* intervals in the full Brillouin zone, the irreducible part of which was used in a self-consistent calculation, is denoted as scK*n*; a mesh with *m* × *m* × *m* intervals in the full Brillouin zone, the irreducible part of which was used in a single-step calculation with an input density from a calculation scK*n*, is denoted as K*m*-scK*n*.

Calculations were performed for fcc Aluminium (space group 225) with lattice parameters *a* between 0.36 nm and 0.41 nm in steps of 0.0005 nm. In certain cases, the step size was chosen as 0.0001 nm.

In our evaluation of the *E*(*a*) data, numerical differentiation plays an important role. This approach, however, requires a trade-off between the conflicting demands for the best possible resolution and the lowest possible noise. Here, we applied two approaches. The first one, called "direct evaluation" henceforth, uses generalizations of the midpoint formula of differentiation. For a given set of data points with equidistant arguments,

$$d^2 \mathbb{E}/da^2|\_{a\_i} = \left(\mathbb{E}(a\_{i+1}) - 2\mathbb{E}(a\_i) + \mathbb{E}(a\_{i-1})\right) / (a\_i - a\_{i-1})^2$$

and

$$d^3 \mathbb{E}/da^3|\_{(a\_l+a\_{l+1})/2} = \left(\mathbb{E}(a\_{l+2}) - \mathfrak{Z}(a\_{l+1}) + \mathfrak{Z}(a\_l) - \mathbb{E}(a\_{l-1})\right) / (a\_l - a\_{l-1})^3$$

were evaluated for (*ai* − *ai*−1) = 0.0005 nm. This value is a carefully chosen compromise between a high noise level for small step sizes and a suppression of physical singularities for large step sizes. A considerable further noise reduction can be achieved by starting from third-order polynomial fits within a moving window. This second approach, called "moving window differentiation" below, was carried out on equidistant data sets with (*ai* − *ai*−1) = 0.0001 nm and 33 data points in each fit.

#### **3. Results and Discussion**

*3.1. Effect of the k-Mesh on the Density of States*

The density of states (DOS), *g*, is a primary quantity obtained in any electronic structure calculation on extended systems. It has widely been used for the interpretation of numerous experimental results like photon- or electron-spectroscopy data, thermodynamic properties, or electron-nucleus interaction data. Here, we focus on the DOS at the Fermi level *g*(*e*F). It is related to the low-temperature specific heat and also to the Knight shift. Last but not least, singularities in the dependence of *g*(*e*F) on the volume (or lattice parameter, or pressure) signal Lifshitz transitions. The aim of this subsection is to sensitize the reader to the importance of choosing a dense *k*-mesh in order to achieve a precise DOS. Henceforth, we use *g*(*e*F) as abbreviation for *g*(*e*F(*a*), *a*).

Figure 1 shows *g*(*e*F) of fcc Aluminium vs. lattice parameter *a* for different *k*-meshes. Consider the upper panel of the figure first. The default *k*-mesh scK12 (red line) provides a smooth behavior of *g*(*e*F) in the whole range of the lattice parameter. It does, however, not meet the expectation of a root-like dependence for a nearly-free electron situation. For the next denser mesh scK24, a local maximum of *g*(*e*F) appears close to *a* = 0.394 nm. This maximum could be related to an ETT, but the yet denser mesh scK50 produces a much more structured curve with several local maxima or kinks. The majority of these singularities can be assigned to numerical minigaps arising from wrong band connections. Such minigaps are more numerous but smaller in the scK100 results, and they are almost invisible in the scK200 curve.

The latter two data sets each exhibit two kinks that indeed originate from ETTs. They are marked by the vertical lines and named ETT1 and ETT2. Furthermore, data obtained with a very dense *k*-mesh K400-scK200 are presented in the lower panel of the figure and in its inset. The related curve is piece-wise quasi-analytic. Here, another but less prominent Lifshitz transition, ETT0, is recognizable close to *a* = 0.392 nm, see the inset of Figure 1.

**Figure 1.** Total density of states at the Fermi level vs. lattice parameter, obtained with the FPLO code using different *k*-meshes as specified in the legends. The positions of two prominent ETTs are indicated by vertical lines in the upper panel and by arrows in the lower panel. Another, less prominent ETT is indicated by an arrow in the inset, which has an *a*-axis identical to that of the main graph. The upper panel shows data points which were obtained at a distance of 0.0005 nm. They are connected by straight lines. The lower panel presents data obtained at the same distance completed by data at a distance of 0.0001 nm in the vicinity of the ETTs. The data points are connected by lines in the lower main panel and depicted as circles in the inset. The K400-scK200 data were computed with tightened convergence criteria of 10−7 (density) and 10−9 Hartree (total energy) in the scK200-step.

All three ETTs detected above can be identified in the schematic Fermi surface plots shown in Figure 2. Therein, the well-known third-zone monster of fcc Aluminium is shown in the left-hand part within the complete Brillouin zone and as the lowest detail in the right-hand part, both for the equilibrium situation.

Consider the right-hand part of Figure 2. Upon reduction of *a*, the *β*-orbit merges with its neighbor and the monster become multiply connected (arrow in the second lowest detail). This ETT0 is barely visible in the DOS, compare the inset of Figure 1. A slight further reduction of *a* disconnects the Fermi surface and lets the *α*-orbits vanish, as shown in the third detail. This transition was named ETT1 above. Finally, the ellipsoid around symmetry point U vanishes in another Lifshitz transition, named ETT2. The latter two transitions are well-resolved in the DOS, provided it is computed with appropriate precision, see Figure 1.

To avoid the occurrence of the abovementioned artificial singularities of the DOS, resulting from wrong band connections, broadening techniques are frequently used for *k*-integration and even for the calculation of the DOS. Such approaches, however, do not allow for the identification of van Hove singularities as indicators of ETTs. Hence, the use of a linear tetrahedron method with a very dense *k*-mesh is the only safe way to identify an ETT in the DOS. An alternative way is to calculate Fermi surfaces using an equally dense mesh. We stress that a band structure at symmetry lines is usually not sufficient for this

aim since changes of the Fermi surface topology frequently happen at general *k*-points, compare Figure 2.

**Figure 2.** Schematic Fermi surface of Aluminium as obtained with the FPLO code. **Left part**: Thirdzone equilibrium-volume Fermi surface in the fcc Brillouin zone with symmetry points. Blue (red) color indicates the occupied (unoccupied) side of the Fermi surface. **Right part**: variation of a detail of the Fermi surface close to one of the square edges as indicated by a brown rectangle in the left part. The lattice parameter shrinks from bottom to top. The lowest detail shows the equilibriumvolume situation with three extremal orbits (*α*, *β*, *γ*) indicated in the common notation [9] . Our denominations ETT0, ETT1, and ETT2 refer to the Lifshitz transitions detected in Figure 1.

The limited precision of DOS data evaluated by using a coarse *k*-mesh also results in large errors in the predicted values of Sommerfeld parameters, Pauli susceptibilities, and Knight shifts. To illustrate this remark, Table 1 compiles the relative differences between the well-converged value of *g*(*e*F) obtained from K400-scK200 and *g*(*e*F) data obtained from less dense meshes for three specific lattice parameters. It shows that the default mesh, which is in many cases sufficient for fairly accurate self-consistent calculations, provides an error of up to 30% in *g*(*e*F) at the considered lattice parameters. Even the rather dense mesh scK100 yields an error of up to 2%.

**Table 1.** Relative differences 1 − *g*(*e*F)other/*g*(*e*F)400 between *g*(*e*F)400, computed from K400-scK200, and *g*(*e*F)other, computed from other *k*-meshes, for lattice parameters close to the equilibrium structure and close to ETT1 and ETT2. All data were obtained with the FPLO code.

The data presented in Table 2 corroborate the above note that *k*-meshes of usual density are sufficient for the computation of integral properties like the equilibrium lattice parameter *a*0. For each of the considered *k*-meshes, *a*0 was obtained as the minimum of a parabola through the three total energy values computed for *a* = 0.40415, 0.40420, 0.40425 nm. These lattice parameter values are close to the GGA equilibrium which we identify with *a*0,scK400 = 0.4042144 nm, denoting *a*0 computed with scK400. Already for scK36, the relative error 1 − *a*0,scKxx/*a*0,scK400 is below 2 × 10−5. Note, this precision is somewhat better than the stability of the Aluminium sample lattice parameter in an experiment without explicit and accurate temperature control. For this situation, temperature stability of 1 K and thermal expansion of 2.3 × 10−5 K−1 are assumed. After increasing the density of *k*-points by a factor of 8, that is for scK72, the deviation is already reduced to 2.5 × 10−6, far better than the achievable precision of any single X-ray experiment.

**Table 2.** Relative differences 1 − *a*0,scKxx/*a*0,scK400 between equilibrium lattice parameter *a*0,scK400, computed from scK400, and *a*0,scKxx, computed from other *k*-meshes. All data were obtained with the FPLO code.

#### *3.2. Preliminary Considerations—Wien2k*

Before a serious numerical investigation by means of any DFT code for a periodic system can be started, two important decisions regarding the numerical settings have to be taken. The first one, choosing the method and mesh density for *k*-space integration, was discussed in the previous subsection. The second one is the choice of the basis set for the Bloch wave function representation. In the Wien2k code, the extend of the basis set is encoded in the parameter *Rk*max, the product of the smallest muffin-tin sphere radius and the largest *k*-vector of the plane wave basis set used in the interstitial space. In addition, the number of local orbitals depends on *Rk*max.

In order to check the impact of *Rk*max on the accuracy, we performed calculations for a number of integer values around the default value *Rk*max = 7, using a dense *k*-mesh K100-scK50. Figure 3 shows the total energy vs. *Rk*max of fcc Aluminium for *a* = 0.404 nm, which is close to the GGA equilibrium. The default value of *Rk*max is indicated by a red circle. We note that the total energy diminishes monotonically with growing *Rk*max between 4 ≤ *Rk*max ≤ 12. Hence, the accuracy improves with growing *Rk*max at least till *Rk*max = 12. The energy seems to be converged within about 0.1 meV for this largest tested basis set. The default setting *Rk*max = 7 provides a total energy about 7 meV above the lowest value.

**Figure 3.** Total energy vs. *Rk*max of fcc Aluminium for *a* = 0.404 nm obtained with the Wien2k code. All data were computed with K100-scK50. The lowest obtained total energy value, *E*0 = −6,607,524.417 meV/atom, was chosen as reference. The inset shows part of the data using a zoomed energy axis. Its *Rk*max-axis agrees with that of the main plot. Red circles indicate the result for the default value of *Rk*max and black circles all other calculated points.

Next, we turn to the influence of *Rk*max on the precision by evaluating total energy derivatives with respect to the lattice parameter. Figure 4 shows *d*2*E*/*da*2 (left panel) and *d*3*E*/*da*3 (right panel) vs. *a*, as obtained using Wien2k for a lattice parameter step size of 0.0005 nm. Note that, in order to avoid strong overlap, the individual data sets are offset against each other as detailed in the figure caption. For both the second and third

derivatives, the numerical noise is smallest for *Rk*max = 8 or 9, slightly above the default value *Rk*max = 7. For both *Rk*max < 8 and *Rk*max > 9, the precision decreases. We stress that, in the case of third derivatives, the noise level is larger than the absolute mean value of the derivative for all values of *Rk*max. This fact will necessitate data smoothing in the search for ETTs, see Section 3.5.

**Figure 4.** Derivatives of the total energy of fcc Aluminium with respect to the lattice parameter vs. lattice parameter. The total energies were obtained with the Wien2k code for different values of *Rk*max as detailed in the legend. **Left panel**: second derivative, *d*2*E*/*da*2; the data are offset against the original data by (11 − *Rk*max) · 106 meV/nm2. **Right panel**: third derivative, *d*3*E*/*da*3; the data are offset against the original data by (11 − *Rk*max) · 109 meV/nm3. The derivatives are obtained by direct evaluation, see Section 2, for a lattice parameter step size of 0.0005 nm. All data were computed with K100-scK50.

One possible explanation for, above *Rk*max = 9, the noise unexpectedly increasing with increasing number of plane waves could be numerical problems arising from near completeness. In such a situation, the newly added basis vectors contain more and more redundant directions in Hilbert space. Shifting (a tiny) wave function weight from one to another of these directions will not alter the charge density in any appreciable manner. However, which directions are finally chosen to contribute to the wave function is increasingly determined by unavoidable numerical noise from various places of the code in a quasi-random fashion, the less important the added directions are. Such noise will reflect in derivatives of the total energy. We have no hypothesis that could explain the observed increase of noise if a smaller basis than the default one is chosen.

#### *3.3. Preliminary Considerations—FPLO*

The precision of the total energy computation, viz. the noise of *E*(*V*), can be evaluated in two ways. In the previous subsection, we used numerical differentiation for this aim. Now, we consider the comparison to a global fit of the *E*(*V*) data set: the subtraction of a Birch-Murnaghan fit, *E*BM−fit(*V*), from the original data removes a large and smooth contribution.

Figure 5 compares the results of this procedure when applied to two FPLO data sets obtained with different treatment of orbital-potential matrix elements. In both cases, the rather dense *k*-mesh K200-scK50 and default settings concerning all other parameters were used. The black and red curves present the lattice parameter dependences of the difference *E*(*a*) − *E*BM−fit(*a*) obtained with default and modified integration meshes for the orbital-potential matrix elements, respectively. Along the whole range of considered lattice parameters, these energy differences amount to less than 10−7 parts of the total energy.

**Figure 5.** Difference between the total energy calculated with the FPLO code and the total energy obtained with a related third-order Birch-Murnaghan fit vs. lattice parameter. All data were obtained with K200-scK50. The black line shows data obtained with default numerical settings and the red line shows data obtained with a modified mesh for the computation of orbital-potential matrix elements as explained in the text. The arrows indicate the positions of two ETTs, identified in Section 3.1. The vertical line marks the equilibrium lattice parameter derived from the Birch-Murnaghan fit to data obtained with the modified integration mesh.

The general behavior of both *E*(*a*) − *E*BM−fit(*a*) curves is typical of a fourth-order polynomial, as expected for the residual of a third-order fit. However, the black curve exhibits two well-visible discontinuities at *a* = 0.381 ... 0.3815 nm and *a* = 0.403 ... 0.4035 nm, as well as another, barely visible one at *a* = 0.3705 ... 0.371 nm. The energy residual *E*(*a*) − *E*BM−fit(*a*) jumps at these lattice parameter values by less than 0.05 meV, which is about two orders of magnitude smaller than the energy variations usually encountered. Hence, these discontinuities are not visible in most presentations. They are, moreover, not important in most applications. In the present study, however, they are essential.

A possible idea would be to relate the discontinuities in *E*(*a*) − *E*BM−fit(*a*) to Lifshitz transitions. However, the two ETTs identified in Section 3.1 occur at other lattice parameters, were no singularity in *E*(*a*) − *E*BM−fit(*a*) is visible in Figure 5. The true origin of these discontinuities lies in a time-saving approximation in the evaluation of the orbitalpotential matrix elements. This numerical integration is performed on a three-dimensional

grid, which combines a radial mesh with angular meshes at each of the radial points. In FPLO, computing time is saved by using angular meshes with lower density for smaller radii. Since, the transition from smaller to larger meshes with varying radius happens at reasonably chosen fixed distances from the nuclei, the number of grid points changes with varying lattice parameter, resulting in the observed discontinuities.

The remedy is to choose the most dense angular mesh at all radii, which increases the CPU time only by about 10%. This modification has been applied to compute the data shown in the red curve of Figure 5. There, no visible deviation from a smooth behavior is observed and the noise level is well below 0.01 meV. For this reason, we used the modified mesh for the evaluation of the orbital-potential matrix elements for all calculations presented in this work, except for the default data in Figure 5.

Finally, we remark that such small precision-related effects in the total energy as visualized in Figure 5 can only be detected by considering the difference between total energy and a Birch-Murnaghan or another appropriate polynomial fit. The parameter values of the Birch-Murnaghan fits used here are given in Table 3. They are equivalent to each other, i.e., the precision of these parameters is not reduced by applying the default integration mesh instead of the modified one as long as the fit is performed with a large enough number of data points.

**Table 3.** Parameter values obtained in the third-order Birch-Murnaghan fits used for Figure 5. *V*0 denotes the equilibrium atomic volume, *B*0 and *B* 0 are the bulk modulus and its pressure derivative.

#### *3.4. Comparison of Total Energies*

We now compare the accuracies reached by the Wien2k and FPLO codes. This will be done by considering the equation of state for two different basis sets of each code. Figure 6 shows results obtained with default basis sets (Wien2k: *Rk*max = 7; FPLO: local states 2*s*, 2*p*, 3*s*, 3*p*, 3*d*, 4*s*, 4*p*) and with basis sets for enhanced accuracy (Wien2k: *Rk*max = 9; FPLO: local states 2*s*, 2*p*, 3*s*, 3*p*, 3*d*, 4*s*, 4*p*, 4*d*, 4 *f* , 5*s*, 5*p*). For the Wien2k code, the enhanced basis set was defined according to the precision considerations in Section 3.2; for the FPLO code, the enhanced basis set as defined in the appendix to Ref. [1] was used. We remark that the *Rk*max settings for Wien2k in the mentioned reference were 6.5 for default and 10.0 for enhanced accuracy calculations [1].

We learn from Figure 6 that the characterization of Wien2k as a highly accurate code is once more confirmed. For the default settings, the ground-state energy obtained by FPLO is about 10 meV per atom above the corresponding value obtained by the Wien2k code; considering the enhanced settings, this difference is reduced to about 4 meV per atom. At lattice spacings much smaller than the equilibrium spacings, the order is, however, reversed for the enhanced settings: the FPLO energy is found to be about 4 meV below the Wien2k value (see upper inset). A larger value of *Rk*max, which could be chosen at the price of reduced precision, would not change this picture, since it would reduce the total energy by merely 0.5 meV compared with the value for *Rk*max = 9. We remark that the different sequences of the total energy values obtained by Wien2k and FPLO for equilibrium spacing, on the one hand, and for much smaller values of *a*, on the other hand, are not very surprising: DFT codes are usually optimized to achieve high accuracy close to the equilibrium state. For the case of high pressure, numerical settings may have to be adapted specifically.

**Figure 6.** Comparison of total energies of fcc Aluminium vs. lattice parameter calculated with two different codes and two different basis sets for each code as detailed in the legend. All data were obtained with K100-scK50. The lowest total energy, *E*0 = −6,607,523.903 meV/atom, was chosen as reference. The insets show zoomed regions of the main plot. The *a*-axis of the lower inset agrees with that of the main plot.

Table 4 presents a comparison of the parameter values obtained by means of Birch-Murnaghan fits to all four data sets. These fits yield virtually the same results: *a* = 0.4040 . . . 0.4042 nm, bulk modulus *B*0 = 79 GPa, and its pressure derivative *B* 0 = 4.3.

**Table 4.** Comparison of third-order Birch-Murnaghan fits applied to the data shown in Figure 6. Rows denoted with *Rk*max values refer to Wien2k, rows denoted with a basis type refer to FPLO. *V*0 denotes the equilibrium atomic volume, *B*0 and *B* 0 are the bulk modulus and its pressure derivative, respectively.

It is very remarkable that two completely independent codes yield total energies with an absolute difference of less than 5 meV and a relative difference of less than 10−6 along the whole range of lattice parameters. These findings confirms the reliability of *both* codes because the probability of accidental error compensation is extremely low. Note finally that Aluminium is usually considered to be a nearly-free-electron metal. Thus, it is expected to be best described by plane-wave related methods. The above results prove that, even in this case, a local orbital based code achieves almost the celebrated accuracy of an augmented plane wave code.

We would like to conclude this evaluation of accuracies with a comparison of the related numerical efforts. One anticipates that a better accuracy requires a larger effort. This expectation is indeed confirmed by the typical computation times presented in Table 5. We observed that, as a rough estimate, the Wien2k calculations consume about one order of magnitude more CPU resources than the FPLO calculations. This difference originates mainly from the following: Wien2k uses about twice as many basis functions as contained in the enhanced basis of FPLO.

**Table 5.** Typical computation times on Intel Xeon CPUs with 2.6 GHz clock frequency to calculate the data shown in Figure 6. Rows denoted with *Rk*max values refer to Wien2k, rows denoted with a basis type refer to FPLO. The runtime of each of the cores for one iteration cycle is called *t*core. It is given for serial mode (1 core) and, in the case of Wien2k, for parallel mode with 16 cores.

#### *3.5. Search for ETTs by Means of Total Energy Derivatives*

In Section 3.1, carefully inspecting the lattice parameter dependence of the density of states at the Fermi level, we detected three Lifshitz transitions. Two of them, termed ETT1 and ETT2, are clearly visible in *g*(*e*F), see Figure 1, provided the *k*-mesh is K100 or better; another one, ETT0, is hardly visible. Close to a three-dimensional ETT related to the bottom of a band at *e*c, the total energy contains a contribution proportional to (*e*F − *e*c)2.5 [15]. In other situations, the sign can be changed. Assume linear dependence of *e*F − *e*c on *a* close to the transition. Consequently, the ETT should be signaled by a divergence of the third derivative *d*3*E*/*da*3.

We now ask whether the respective optimum total energy precision of Wien2k and FPLO is sufficient for the identification of such divergences. The upper panel of Figure 7 shows *d*3*E*/*da*3 vs. *a* as obtained by means of both codes with K100-scK50 and subsequent direct numerical differentiation (d). In both cases, the step size of *a* was set to 0.0005 nm. Despite the use of rather dense *k*-meshes and although the settings of both codes were optimized with respect to low numerical noise, there is no sign of a divergence visible in the upper panel data close to any of the expected critical lattice parameters, which are indicated by vertical lines.

The fluctuations of the Wien2k data (black curve) are considerably stronger than those of the FPLO data (red), which points to a better precision of the latter. We remark that FPLO data generated with enhanced basis set as described in Section 3.4 (not included in Figure 7) exhibit up to four times stronger noise of *d*3*E*/*da*3 than data obtained using the default basis (red curve).

The lower panel of Figure 7 shows data obtained for the setting K200-scK100 using yet denser *k*-meshes than above (blue, Wien2k and orange, FPLO), and for K400-scK200 (green, FPLO). The two former data sets were computed by using a moving window differentiation (m) with 33 equidistant data points (distance 0.0001 nm). In this way, a clear singularity is seen at the expected position of ETT2, which is yet more pronounced in the FPLO K400-scK200 data set which was obtained by direct evaluation of the derivative (distance 0.0005 nm). On the one hand, the moving window differentiation reduces the numerical noise considerably; on the other hand, it broadens the singularity which is more pronounced in the latter data. Note that the individual data are offset against each other as detailed in the figure caption.

**Figure 7.** Third derivative of the total energy with respect to lattice parameter *a* in relation to electronic topological transitions. **Upper panel**: Comparison of Wien2k data obtained with *Rk*max = 8 (same data as shown in the right panel of Figure 4, black curve) with FPLO data (red). Both calculations were performed using the K100-scK50 setting with total energy convergence of 10−8 Rydberg. The derivatives were obtained by direct evaluation, indicated in the legends by "d". **Lower panel**: Comparison of Wien2k data obtained with *Rk*max = 8 (blue, with 107 meV/nm3 offset) with FPLO data (orange, with 0.5 × 107 meV/nm3 offset). Both data sets were obtained using the K200-scK100 setting with total energy convergence of 10−8 Rydberg. The derivatives were calculated using a moving window with 33 data points, indicated by "m", see Section 2. Another FPLO data set, K400-scK200 with total energy convergence of 10−9 Hartree and direct evaluation of the derivative, is shown in green. The vertical lines indicate the positions of ETT1 and ETT2 as obtained in Section 3.1.

None of the presented data sets shows any peculiarity exceeding numerical fluctuations at the position of ETT1 or of ETT0. We remind that the root-like singularities in *g*(*e*F) caused by ETT0 and ETT1 are less strong than that one caused by ETT2, see Figure 1. This difference may result from the differing natures of the respective topological changes: ETT0 and ETT1 are the opening and the closing of connections between Fermi surface parts, whereas ETT2 is the vanishing of an isolated part of the Fermi surface, see Figure 2.

To conclude this section, we have shown that an electronic topological transition can in principle be detected in a higher derivative of the total energy. To this end an extreme precision is required, such that the total energy curve appears very smooth as compared to usual code applications. Our calculations have shown that ETT2 is indeed visible as a singularity in *d*3*E*/*da*3(*a*), whereas we could not detect the other ETT in this way. The other kinks observed are not related to topological transitions. In consequence, the total energy is not suited to detect such transitions. This implies that similar analyses of experimental data, which have a naturally lower precision than a highly tuned computer code, will be challenging. These strict statements, however, do no longer hold if the ETT is of lower than three-dimensional nature, in which case the transition can become first order and more easily detectable [16]. The differentiation procedures used here are very sensitive and might be helpful both in detecting first or second order phase transitions or in disproving

their existence, such as in the case of *β*-Zirconium. There, contrary to earlier claims, no evidence for isostructural phase transitions was recently found [17].

#### *3.6. Birch-Murnaghan Fits Using Minimum Number of Data Points*

We have shown in Section 3.1 that the calculation of the equilibrium lattice constant is well converged already for moderately large *k*-meshes. Now, we would like to find out whether it is possible to utilize precise total energy data for evaluating the equilibrium volume, the bulk modulus, and its pressure-derivative with the minimum number of data points needed. For this aim, third-order Birch-Murnaghan fits were performed with only four total energy values obtained using FPLO with K100-scK50 setting. The equidistant fit points are distributed within a fit region of 0.0015 ... 0.012 nm around the equilibrium lattice parameter.

Figure 8 shows results for the three experimentally accessible fit parameters vs. the width of the fit region. For comparison, the results of two Birch-Murnaghan fits to 25 and 13 data points, respectively, are shown as well. As expected from the results presented in Table 2, the equilibrium volume is almost independent of the width of the fit region. Furthermore, the variation of the bulk modulus, 0.1 GPa within the considered range, is negligible. Only the value of *B* 0 seems to be unstable if the fit region is chosen smaller than 0.004 nm. This observation is consistent with the fact that the stability of numerical derivatives worsens with their order increasing.

**Figure 8.** Results of third-order Birch-Murnaghan fits to few data points. The equilibrium volume *V*0 (upper panel), the equilibrium bulk modulus *B*0 (middle panel), and its pressure-derivative *B* 0 (lower panel) were obtained from FPLO total energy data for K100-scK50 setting. The circles connected with black lines show results obtained from four *E*(*a*) data points, centered at the equilibrium lattice parameter and equally distributed over the fit region. The red (blue) lines show fit results obtained from 25 (13) *E*(*a*) data points equally distributed within 0.012 (0.006) nm around the equilibrium lattice parameter.

The possibility to use a minimum number of precise data points to numerically evaluate quantities derived from the total energy can be useful e.g. for screening studies. However, this chance comes at the price of a sufficiently dense *k*-mesh. As long as specific

effects of band structure or DOS, such as ETTs, Weyl points or magnetic instabilities, are not important, it can be useful to apply broadening techniques with a less dense *k*-mesh. They offer the advantages to reduce the related effort and to smoothen the selfconsistent iteration.

#### **4. Summary**

Electronic structure calculation packages such as the famous Wien2k code, developed and maintained by Karlheinz Schwarz and his group for decades, are meanwhile indispensable tools for the interpretation of experimental data as well as for the prediction of materials properties. In the latter case, there is no standard to gauge the result but by other codes.

Doing so, we have confirmed the classification of Wien2k as the "gold standard" concerning total energy accuracy. This high numerical accuracy is achieved, as in the case of similar codes [1], by using a large and elaborated basis set for the representation of the Bloch wave function. The FPLO code, taken in this study for comparison, uses a smaller basis set. In consequence, it can only reach slightly lower total energy accuracy than Wien2k, but it has the advantage of considerably smaller numerical effort.

Either of the two compared codes allows for the detection of one of the known electronic topological transitions of Aluminium under pressure by appropriate evaluation of *d*3*E*/*da*3(*a*). The other two discussed ETT, which are of different topological nature than the former one, are only visible in the density of states and in the Fermi surface. A surprising finding is the high total energy precision of FPLO, demonstrated by comparing *d*3*E*/*da*3 with related Wien2k data.

The total energy precision can be improved in both codes by refining the *k*-mesh. While this behavior is expected, further strong improvements of the numerical precision are obtained in the case of Wien2k by optimizing the value of *Rk*max and in the case of FPLO by modifying a real-space integration mesh. Thus, it also became possible to perform meaningful third-order Birch-Murnaghan fits to the minimum number of four data points of the obtained very precise *E*(*V*).

Concerning both codes, a counterintuitive finding has to be mentioned: tuning the respective basis sets toward optimum numerical accuracy results in a reduced numerical precision. This feature presumably originates from incipient linear dependence of the basis functions.

Finally, we note that our approach to quantify the numerical precision by means of total energy derivatives is not only applicable to the two codes considered in the present study, but can easily be extended to other methods.

**Author Contributions:** Conceptualization and methodology, M.R., A.M. and H.R.; software, K.K.; investigation, all authors; writing—original draft preparation, M.R. and K.K.; writing—review and editing, all authors. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** We are are grateful to Ulrike Nitzsche for her efforts in establishing and maintaining the IFW high-performance computational environment.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Density Functional Theory Study of Metal and Metal-Oxide Nucleation and Growth on the Anatase TiO2(101) Surface**

**Leila Kalantari \*, Fabien Tran and Peter Blaha**

Institute of Materials Chemistry, Vienna University of Technology, Getreidemarkt 9/165-TC, A-1060 Vienna, Austria<EMAIL_ADDRESS>(F.T<EMAIL_ADDRESS>(P.B.)

**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** Experimental studies have shown the possible production of hydrogen through photocatalytic water splitting using metal oxide (MO*y*) nanoparticles attached to an anatase TiO2 surface. In this work, we performed density functional theory (DFT) calculations to provide a detailed description of the stability and geometry of M*x*O*y* clusters M = Cu, Ni, Co, Fe and Mn, *x* = 1–5, and *y* = 0–5 on the anatase TiO2(101) surface. It is found that unsaturated 2-fold-coordinated O-sites may serve as nucleation centers for the growth of metal clusters. The formation energy of Ni-containing clusters on the anatase surface is larger than for other M clusters. In addition, the Ni*n* adsorption energy increases with cluster size *n*, which makes the formation of bigger Ni clusters plausible as confirmed by transition electron microscopy images. Another particularity for Ni-containing clusters is that the adsorption energy per atom gets larger when the O-content is reduced, while for other M atoms it remains almost constant or, as for Mn, even decreases. This trend is in line with experimental results. Also provided is a discussion of the oxidation states of M5O*y* clusters based on their magnetic moments and Bader charges and their possible reduction with oxygen depletion.

**Keywords:** DFT; anatase TiO2(101) surface; adsorption energy; Bader charge

#### **1. Introduction**

The over-exploitation of the fossil energies leads to a significant increase of CO2 in the atmosphere, resulting in severe climate problems. It is therefore absolutely necessary to replace them by alternative energy sources. When using sun light as energy source, however, materials with specialized properties are necessary, and TiO2 is such a versatile material that has numerous applications in catalysis, photocatalysis, and solar energy [1]. TiO2 crystallizes in three major different structures: rutile (tetragonal), anatase (tetragonal) and brookite (rhombohedral). Other structures exist as well, as for example cotunnite, that has been synthesized at high pressures and is one of the hardest polycrystalline materials known. However, only rutile and anatase play an important role in the applications of TiO2. A well-suited model compound for photocatalytic water splitting is anatase. It has a slightly larger band gap than rutile (∼3.0 eV for rutile and ∼3.2 eV for anatase) [2–6] and shows better performance [7], since it also has a longer electron-hole pair life time, which makes anatase more suitable for photocatalytic applications [8]. To improve the water splitting photocatalyst activity, the use of co-catalysts is necessary and the most widely used are Pt and Pd for photocatalytic reduction, and IrO2 and RuO2 for the oxidation of water. However, these are expensive materials based on rare noble metals. To achieve large-scale industrial applications of photocatalytic water splitting, the development of new co-catalysts based on cheap and widely available elements remains an important issue. With regard to novel abundant co-catalysts for photocatalytic water splitting, research efforts focus on 3*d* transition metals (M) oxides, which are known for their excellent catalytic properties and applications in industry, research and nature.

**Citation:** Kalantari, L.; Tran, F.; Blaha, P. Density Functional Theory Study of Metal and Metal-Oxide Nucleation and Growth on the Anatase TiO2(101) Surface. *Computation* **2021**, *9*, 125. https://doi.org/10.3390/computation 9110125

Academic Editor: Henry Chermette

Received: 12 October 2021 Accepted: 5 November 2021 Published: 19 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Recently, small clusters of M oxides were produced by deposition of M-acetylacetonate precursors on TiO2 and subsequent calcination in air. In particular, the Ni (and Cu) nanoparticles showed pronounced activities for water splitting, but little is known about the atomic structure of these nanoparticles [9]. This inspired the present work, where we used density functional theory (DFT) [10,11] to study possible adsorption sites of a 3*d* transition metal atom (Cu, Ni, Co, Fe, and Mn), which can be the nucleation sites for bigger M clusters (2 to 5 atoms) or M oxides clusters on TiO2(101) anatase, the lowestenergy surface of anatase [12]. There are a couple of previous theoretical works on the M adsorption on anatase, studying single M atom adsorption and diffusion [13], or the adsorption of small Ni [14] or Cu and Cu-oxide [15,16] clusters. A nice review related to this topic is given by Zhou and Dong [17]. However, note that most of these calculations used a generalized gradient approximation (GGA) in DFT for the correlated 3*d* electrons, which may not be accurate enough.

#### **2. Methods**

Spin-polarized DFT calculations of the adsorption of M*x*O*y* clusters (M = Cu, Ni, Co, Fe or Mn, *x* = 1, 2, 3, 5, and *y* = 0, 1, 2, 3, 4, 5) on the TiO2(101) surface were performed using the full-potential augmented plane wave plus local orbitals method [18,19] as implemented in the WIEN2k code [20,21]. We employed the PBEsol exchange-correlation functional [22], which is a GGA and yields lattice parameters of bulk anatase TiO2, *a* = *b* = 3.77 Å and *c* = 9.54 Å that are in good agreement with experimental data [23]. In order to treat the correlated 3*d* electrons of the M atoms a Hubbard *U* correction was used [24]. We used *U*eff = *U* − *J* of 5 eV for Cu, Ni, and Co and 4 eV for Mn and Fe. These values are deliberately chosen smaller than what is typically used in strongly correlated oxides, since in our case also less ionic (metallic) clusters were investigated and the correlation may not be always so strong. We note that small changes of *Ueff* would not affect our basic results and conclusions.

The (101) surface was modeled by a symmetric slab with a thickness of three layers of TiO2 (see Figure 1, where only one layer is shown for better visibility) and a vacuum region of 16 Å between the slabs. In order to reduce the interactions between the adsorbed clusters in neighboring cells a 3 × 1 supercell was used for all surfaces with adsorbed M atoms. Moreover, such a supercell is also helpful to release the interfacial strain energy and find more stable configurations. A 2 × 2 × 1 Γ-centered mesh of **k**-points is used during relaxation (a 5 × 4 × 1 mesh for the final results) and a basis-set size corresponding to *R*min MT *K*max = 7, where *R*min MT is the smallest atomic sphere radius and *K*max the largest reciprocal lattice vector, was used. All surface models were relaxed until all residual forces were below 1 mRy/bohr.

Consistent *R*MT and *R*min MT *K*max values have been used for calculating the adsorption energy *E*ads of clusters on anatase surface, which is defined as

$$\begin{array}{rcl} E\_{\text{ads}} &=& \frac{1}{\mathbf{x} + \mathbf{y}} [E(\text{TiO}\_2(101)) + \mathbf{x}E(\text{M}) + \mathbf{y}\frac{1}{2}E(\text{O}\_2) \\ &- E(\text{M}\_{\text{x}}\text{O}\_{\text{y}}/\text{TiO}\_2(101))], \end{array} \tag{1}$$

where *E*(M) and *E*(O2) are the total energies of the free M atom and an O2 molecule, and *E*(TiO2(101)) is for the bare TiO2(101) surface. Bader's quantum theory of atoms in molecules [25] was used to calculate the atomic magnetic moments and charges.

**Figure 1.** Side and top views of the top-most layer of the TiO2(101) surface (in the 3 × 1 supercell) with a single adsorbed Ni atom (middle) and two Ni atoms (right structure). Red, blue and gray spheres correspond to O, Ti and Ni atoms, respectively. Selected Ni-O bond distances (in Å) are also given. O3*c* and O3*c* are threefold coordinated oxygen with two Ti6*c* and one Ti5*c*, and two Ti5*c* and one Ti6*c* nearest neighbors, respectively.

#### **3. Results and Discussion**

#### *3.1. Perfect Anatase TiO*2*(101) Surface*

As shown in Figure 1, the anatase TiO2(101) surface has a stepped structure. Threefold coordinated O atoms (O3*c* and O3*c*) and sixfold coordinated Ti atoms (Ti6*c*) are fully saturated and have bulk-like coordination, whereas the twofold O2*c* and fivefold Ti5*c* atoms are under-coordinated. The O2*c* atoms are located at the ridges of the saw-tooth-like structure, and after optimization they relax inwards by ∼0.16 Å with respect to bulk TiO2. The three-fold coordinated O3*c* and O3*c* atoms relax outwards by ∼0.07 Å, while the Ti5*c* and Ti6*c* atoms relax inwards by ∼0.15 Å and ∼0.16 Å respectively, so that the surface exhibits a slightly buckled structure.

#### *3.2. Pure Metal Clusters on the Anatase TiO*2*(101) Surface*

Upon calcination of the M-acetylacetonate precursors, single M atoms may adsorb, diffuse and form larger clusters and finally oxidize on the surface. It is thus natural to study as a first step the possible adsorption sites of a single M atom or a small cluster. We have considered several possible initial configurations of the adsorbed metal atoms and clusters M*x* with *x* = 1, 2, 3, and 5. The most stable structures, which are the same for all M (Cu, Ni, Co, Fe or Mn), are presented in Figures 1–3 in the case of Ni clusters with *x* = 1, 2 and 5. The nearest-neighbor distances are shown in Table 1, while the adsorption energy *E*ads and magnetic moment are available in Table 2.

**Figure 2.** Side and top views of the top-most layer of the TiO2(101) surface with an adsorbed Ni5 cluster in a chain-like structure. The color coding is as in Figure 1.

**Figure 3.** Side view of the starting (**left**) and relaxed structures of Ni5 (**middle**) and Fe5 clusters (**right**) adsorbed on the TiO2(101) surface in a crystalline structure. The color coding is as in Figure 1.

**Table 1.** Shortest M-O bonds (in Å) between a M atom (of the pure metal clusters adsorbed on the TiO2(101) surface) and O atoms (O2*c* and O3*c*) of TiO2. M5(1) and M5(2) refer to chain-like and crystalline structures, respectively, and Mn5(2a) corresponds to a bcc like structure with antiferromagnetic configuration. The results for different M atoms are separated by semi-colons. "(2)" means two bonds of same length.

A single adsorbed M atom prefers to adsorb at the bridge site between two unsaturated oxygens (O2*c*) atoms, similar as found in previous studies [13,14]. The corresponding M-O bond length (it is the shortest for Ni, 1.84 Å) follows in general the trend found in the bulk M oxides (except for the Cu-O and Ni-O distances, which are interchanged), but in all cases is much shorter (and therefore indicates a stronger bonding) than in the bulk MO (e.g., 2.08 Å in NiO). The particular short bond distance for Ni goes hand in hand with the O2*c*-Ni-O2*c* bond angle, which reaches 177◦ and forms a nearly linear bond, while for the other metals this varies from 167◦ (Fe) to 150◦ (Co). On the other hand the two larger M-O3*c* bond distances are also almost equal (but very large) for Ni, but for all other M form a short and long bond in an asymmetric position (see Table 1). All M atoms except Cu possess a magnetic moment (see Table 2), where as expected Mn has the largest moment, while it is the smallest for Ni. It should be noted that PBEsol without *U* would lead to a

non-magnetic Ni atom and the moments for the other M would be smaller but still sizeable. The adsorption energy on the TiO2(101) surface (Table 2) is the largest for Ni (3.46 eV) and the smallest for Cu (2.24 eV). Comparing our results with literature the adsorption energy of Ni is 0.4 eV larger than in references [13,14] and the Ni position is symmetric with respect to the O3*c* atoms. On the other hand, our calculated adsorption energy for Fe is 0.7 eV smaller than that in reference [13], while it agrees well for Cu and Co. These differences could arise due to the neglect in references [13,14] of a Hubbard *U* correction for the M atom or the use of pseudopotentials.

**Table 2.** Adsorption energy (in eV/atom) and atomic magnetic moment (in *μ*B) for pure metal clusters adsorbed on the TiO2(101) surface. M5(1) and M5(2) refer to the chain-like and crystalline structures, respectively, and Mn5(2a) corresponds to the bcc structure with antiferromagnetic configuration. The cohesive energies for bulk MO are 2.75, 4.20, 3.72, 3.41, and 3.34 eV for Cu, Ni, Co, Fe, and Mn, respectively.

If 2 or 3 M atoms are adsorbed on the surface, the most stable structures are still obtained when the M atoms sit at the O2*c* bridge sites (the structure on the right side in Figure 1 shows 2 M atoms adsorbed on the surface), very similar to a single atom adsorption. Interestingly, for all M atoms except Ni the adsorption energy per M atom decreases slightly (Figure 4 and Table 2). For Ni it increases by 0.65 eV. The reason seems to be that the Ni atom is smaller than the others and an almost linear chain fits perfectly for Ni but not in the other cases, where this additional strain reduces *E*ads slightly. Most magnetic moments remain unchanged, except for Ni, where it increases to 0.97/0.98 *μB*. In our 3 × 1 supercell, 3 M atoms form an infinite linear chain of M-O2*c* pairs.

For 5 adsorbed M atoms, we tried several different starting structures in the 3 × 1 supercell. The most promising structures are based on the chain-like structure with 2 additional atoms (labelled M5(1), see Figure 2) or a structure with a more "crystalline type" nucleated on a single M bridge atom (labelled M5(2), see Figure 3). The additional

2 atoms "destroy" the symmetric chain-like 3-atom structure above the O2*c* atoms forming a quite irregular array of M atoms with bonds to O2*c*, O3*c* and O3*c* atoms. Still this structure seems to be the most stable for Fe and Mn although *E*ads is smaller than for less M atoms. The crystalline M5(2) structures behave very differently for the bcc metals Mn and Fe and the fcc metals Co, Ni and Cu. While for the former the cluster structure stays basically close to the starting one and form a "bcc" like cluster (Figure 3), for the latter the cluster rearranges and can be considered as a nucleus for an fcc like cluster. In fact, their *E*ads is larger than for the corresponding chain-like M5(1) structures and for Ni and Cu even the highest within the series.

**Figure 4.** Adsorption energy (in eV/atom) for different metal clusters M*x*. M5(1) and M5(2) correspond to chain-like and crystalline structures, respectively.

For Mn5 in a crystalline bcc structure, we also started from different antiferromagnetic configurations (labelled Mn5(2a)), but the adsorption energies were not much larger than in the ferromagnetic case.

The partial density of states (PDOS) for the TiO2(101) surface and for the energetically most favorable configuration of one metal atom on the surface are shown in Figure 5. For the TiO2(101) surface, both the valance band maximum (VBM) and the conduction band minimum (CBM) are composed of a mixture of O-2*p* and Ti-3*d* states, with the CBM having mostly 3*d* character while the VBM has more O-2*p* character, which indicates a mixed ionic and covalent bonding. After adsorption of a M atom on the anatase TiO2(101) surface, the Fermi level moves to the conduction band and the system becomes formally metallic. In the case of adsorbed Cu and Ni atoms, almost all the *d* states of Cu/Ni are concentrated in the band gap below the Fermi energy, and they slightly hybridize with the O-*p* and Ti-*d* states of the Ti and O atoms which are the closest to the Cu/Ni atom. With an adsorbed Co, Fe, or Mn atom the *d* states are not only in the band gap but also show a peak at the VBM for spin-up and inside the conduction band for spin-down states.

**Figure 5.** PDOS for (**a**) the bare TiO2(101) surface and (**b**–**f**) for one metal atom adsorbed on the TiO2(101) surface. The Fermi level is indicated by a vertical dashed line. All Ti and O PDOS are from surface atoms. The *d* partial DOS of the metals are scaled down by a factor of 3. The plots are aligned at the TiO2 CBM and the Fermi energy is set at 0 eV.

#### *3.3. MxOy on the Anatase TiO*2*(101) Surface*

Turning now to the adsorption of metal-oxides clusters M*x*O*y* on the anatase TiO2(101) surface, several 2D and 3D (starting) configurations among the numerous possible ones were considered. Some of the most stable structures are shown in Figures 6 and 7, and the adsorption energies are shown in Table 3 for the Ni*x*O*y* clusters. It can be seen that *E*ads is almost constant with the increase of cluster size indicating that bigger Ni clusters can easily be formed on the anatase TiO2(101) surface, which is in agreement with the experimental result [9]. The most favorable structure is Ni10O9 with a NaCl like structure (shown in Figure 7).

It was found experimentally [9] that the oxidation state of the M atoms changed after a hydrogen evolution reaction (HER) experiment. In order to study the M oxidation state theoretically, we started from the chain-like M5 clusters (Figure 6) and added oxygen atoms systematically until a fully oxidized M cluster has been obtained. Figure 8 compares the adsorption energy of different M5O*y* M oxide clusters on the surface and also includes for comparison the cohesive energy of the bulk pure metals and bulk metal monoxides. It can

be seen that, among all metal-oxide clusters, Ni*x*O*y* have the highest adsorption energy. All adsorption energies are lower than the cohesive energies of the corresponding bulk monoxide and the pure bulk metal. Thus they form non-crystalline clusters on the surface. It can be seen that Ni clusters are a bit more stable with less oxygen, while for the other metals the adsorption energy remains almost constant, or, in particular for Mn, is reduced with less oxygen. This is in agreement with experimental transition electron microscopy (TEM) results [9], where cluster formation was seen in particular for Ni. In addition, XPS on calcinated NiO*x*-TiO2 samples shows mainly Ni2+, but after HER experiments a substantial increase of Ni0 compared to Ni2+ appears. On the other hand, XPS on MnO*x*-TiO2 samples shows only Mn2+, even after the HER experiments.

**Figure 6.** Top view of the Ni5 cluster (**left**) and, top and side view of Ni5O5 cluster (**middle** and **right**) adsorbed on TiO2(101) surface in chain-like structure. Red, blue, gray and green spheres correspond to O, Ti, Ni and added O, respectively. From Ni5 to Ni5O5, the O atoms are added in the order O*a*, O*b*, O*c*, O*d* and O*e*.

**Figure 7.** Ni10O9 cluster with a NaCl-like structure adsorbed on the TiO2(101) surface. The color coding is as in Figure 6.

Theoretically, the oxidation state is more difficult to estimate. One can use for instance Bader's method [26], which uses the electron density to calculate the gradient vector field and searches for surfaces of zero flux. The charge enclosed within this zero flux surface (basin) can be used to define the total electronic charge of an atom. It must be mentioned that one never gets the full nominal charges, but these Bader charges are always smaller than anticipated. In addition, in particular for open structures like surfaces, these charges may sometimes correlate more with the specific geometry and the distances and number of other atoms around the M atom. For transition-metal atoms, the magnetic moment can also be a good indication of the charge state in some cases. For instance, only Cu2+ ions possess a magnetic moment, while both Cu+ and neutral Cu are non-magnetic. Similarly, one

can quite well distinguish, e.g., Fe2+ and Fe3+ based on their smaller or larger magnetic moments (provided they remain in high-spin states). The magnetic moments and the Bader charges of each atom in M5O*y* clusters are shown in bar charts in Figures 9 and 10. The largest magnetic moments for each M atom increase as expected from 0.6 *μB* for Cu2+ to 1.6, 2.6, 3.7 and 4.9 *μB* in the series Ni, Co, Fe and Mn. This is due to the number of spin-down 3*d* electrons that gets reduced when dealing with high-spin states and fully occupied spin-up states.

**Figure 8.** Adsorption energy (in eV/atom) of M*x*O*y* clusters on the TiO2(101) surface. Also shown is the cohesive energy of bulk MO (left-most data) and bulk M (right-most data).

For the charge state of Cu clusters, both methods mentioned above give quite reliable results. Bader charges of about one and magnetic moments of 0.6 *μB* indicate Cu2+, while Bader charges of 0.6 indicate Cu+ and even smaller Bader charges neutral Cu, both with vanishing magnetic moments. Apparently, in the fully oxidized cluster all Cu are in the 2+ charge state, while removal of an O*e* atom (see Figure 6) leads to a drastic reduction of the charge state of the neighboring Cu4 and Cu5 atoms.

For Ni, however, the situation is more complicated. Even in the fully oxidized cluster some Ni moments are completely quenched because the Ni-O bond length (1.77 Å) is for these Ni atoms much smaller than for the magnetic ions (1.97 Å) and it is difficult to estimate the charge states based on the magnetic moment. On the other hand the Bader charges show a very systematic decrease with oxygen removal, and the experimentally observed increase of neutral Ni correlates reasonably with these Bader charges. Going to the extreme case of Mn, the magnetic moments stay almost constant around 5 *μB* indicative of Mn2+ in a high-spin state in agreement with experiment. The Bader charges are less convincing, but nevertheless they never indicate a neutral charge state of a Mn atom.

In order to get more information about the electronic structure, the PDOS of various Ni-O clusters is shown in Figure 11. For a Ni5 cluster (Figure 11b) we see a wide range of Ni-*d* states within the gap of TiO2, but also some spin-up states mixed into the O-2*p* valence band and some spin-down states in the empty conduction band, indicating the different bonding situation of the different Ni atoms. Upon oxidation a more pronounced formation of upper and lower Hubbard bands can be seen and the remaining Ni-*d* states in the gap sharpen until for the Ni10O9 cluster (Figure 11f) quite sharp Ni-*d* states of (distorted) e*g* and t2*g* character appear similar as in bulk NiO.

**Figure 9.** Atomic magnetic moments (in *μ*B) of M atoms in M*x*O*y* clusters adsorbed on TiO2(101) surface.

**Figure 10.** Bader charges of M atoms in M*x*O*y* clusters adsorbed on TiO2(101) surface.

**Figure 11.** PDOS for (**a**) the bare anatase TiO2(101) surface, showing the surface atoms and the bulk-layer states, (**b**) Ni5, (**c**) Ni5O, (**d**) Ni5O3, (**e**) Ni5O5 clusters in chain-like structures and (**f**) a Ni10O9 cluster with NaCl-like structure adsorbed on the surface. All Ti and O PDOS are from surface atoms only. The Ni PDOS is shaded. The plots are aligned at the TiO2 CBM and the Fermi energy is set at 0 eV.

#### **4. Summary and Conclusions**

Spin-polarized DFT calculations were carried out to investigate possible stable structures of five different transition metal clusters, namely Cu, Ni, Co, Fe and Mn adsorbed on the anatase TiO2(101) surface. The adsorption site for single atoms is in all cases the bridge site between two unsaturated oxygen (O2*c*) atoms. We found that among all tested metals Ni has the largest adsorption energy, and for all metals except Ni the adsorption energies decrease with increased number of adsorbed atoms, so that probably only Ni prefers to form bigger clusters on the TiO2 surface in agreement with experimental TEM results [9]. Clusters with five M atoms can form either a fcc-like structure (Cu, Ni, Co) or still remain in a chain-like structure (Fe, Mn) as is common for smaller clusters.

We also studied several oxidized metals clusters M*x*O*y*, and in order to get more information on the charge state of the M atom a systematic study of M5O*y* (*y* = 0–5) clusters on the TiO2(101) surface was made. We found that Ni forms more stable structures with reduced oxidation, while for the other metals the adsorption energy is almost constant or is reduced. An analysis of the Bader charges and magnetic moments allows to reveal the corresponding charges of the M atoms, and we found, in agreement with the experimental results [9], that Ni can be more easily reduced during a HER experiment, while Mn atoms will keep their 2+ oxidation state.

**Author Contributions:** L.K.: investigation, formal analysis, writing—original draft, software, visualization. F.T.: writing—review and editing. P.B.: supervision, writing—review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the TU-D doctoral college (TU Vienna).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available in insert article.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

12. Hebenstreit, W.; Ruzycki, N.; Herman, G.S.; Gao, Y.; Diebold, U. Scanning tunneling microscopy investigation of the TiO2 anatase (101) surface. *Phys. Rev. B* **2000**, *62*, R16334–R16336. [CrossRef]

## *Article* **Hexatetra-Carbon: A Novel Two-Dimensional Semiconductor Allotrope of Carbon**

**Mosayeb Naseri 1,2, Jaafar Jalilian 3, Dennis R. Salahub 2,\*, Maicon Pierre Lourenço 4 and Ghasem Rezaei 3**

**Abstract:** Employing first-principles calculations based on density functional theory (DFT), we designed a novel two-dimensional (2D) elemental monolayer allotrope of carbon called hexatetracarbon. In the hexatetra-carbon structure, each carbon atom bonds with its four neighboring atoms in a 2D double layer crystal structure, which is formed by a network of carbon hexagonal prisms. Based on our calculations, it is found that hexatetra-carbon exhibits a good structural stability as confirmed by its rather high calculated cohesive energy −6.86 eV/atom, and the absence of imaginary phonon modes in its phonon dispersion spectra. Moreover, compared with its hexagonal counterpart, i.e., graphene, which is a gapless material, our designed hexatetra-carbon is a semiconductor with an indirect band gap of 2.20 eV. Furthermore, with a deeper look at the hexatetra-carbon, one finds that this novel monolayer may be obtained from bilayer graphene under external mechanical strain conditions. As a semiconductor with a moderate band gap in the visible light range, once synthesized, hexatetra-carbon would show promising applications in new opto-electronics technologies.

**Keywords:** density functional theory; hexatetra-carbon; electrical properties

#### **1. Introduction**

The discovery of the interesting behavior of graphene [1] has motivated further theoretical and experimental investigations in order to find possible new stable two-dimensional (2D) materials, especially mono-elemental 2D monolayer materials [2–7]. As carbon exhibits a large number of different allotropes, such as graphite, diamond, C60 fullerene [8], nanotube [9], carbon nano-cone [10], and nanochain [11], 2D carbon mono-elemental monolayer materials beyond graphene have attracted significant attention from both theoretical and experimental fields of study. However, the electronic structure of graphene limits its application in designing electronic nano-devices due to its semi-metallic gapless nature. Therefore, finding new 2D mono-elemental monolayer materials with a semiconductor behavior is technologically important.

In recent years, the possibility of free standing 2D carbon allotropes beyond graphene has been explored. In addition, a number of new carbon allotrope monolayers has been proposed [12–15]. Although these 2D materials, such as graphdiyne [12], penta-graphene [13], and phagraphene [14] are metastable compared with graphene, only a few have been successfully synthesized. Moreover, these new 2D carbon monolayers exhibit very interesting properties, such as anisotropic Dirac cones, inherent ferromagnetism, high catalytic activity, and potential superconductivity related to the high density of states at the Fermi level.

**Citation:** Naseri, M.; Jalilian, J.; Salahub, D.R.; Lourenço, M.P.; Rezaei, G. Hexatetra-Carbon: A Novel Two-Dimensional Semiconductor Allotrope of Carbon. *Computation* **2022**, *10*, 19. https://doi.org/ 10.3390/computation10020019

Academic Editor: Henry Chermette

Received: 27 December 2021 Accepted: 20 January 2022 Published: 25 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

This demonstrates that structural properties and the crystal configuration of 2D carbon allotropes effectively influence the electronic, optical, and other chemo-physical properties of these materials [12–19]. For instance, graphdiyne [12] is a predicted flat one-atom-thick allotrope of carbon with a Dirac cone in its band structure, in which the Dirac points are located at the *K* point. On the other hand, penta-graphene is a new 2D carbon allotrope with semiconducting properties that stabilizes in a buckled structure composed entirely of carbon pentagons, and resembles the Cairo pentagonal tiling with an intrinsic indirect band gap of about 3.25 eV [13]. Moreover, phagraphene, which is a monolayer sheet of carbon with a structure composed of 5-6-7 carbon rings has distorted Dirac cones [14].

On the other hand, carbon nanostructures have tetragons, pentagons, and hexagons as their main basic building blocks. For instance, cubane (C8H8) is a synthetic hydrocarbon molecule formed by eight carbon atoms positioned at the corner of a cube. It is attached to its three neighboring carbon atoms and a hydrogen atom with tetragonal top and side views [20], 2D graphene, and 3D graphite, which are formed by carbon hexagons. In addition, penta-graphene consists entirely of pentagons of carbon atoms. Moreover, the C60 molecule is formed by 12 pentagons, which are separated by 20 hexagons with a soccer ball shape [8].

In this paper, using first-principles calculations based on DFT, we propose a novel 2D elemental monolayer allotrope of carbon, which is called hexatetra-carbon due to its hexagonal and tetragonal top and side views. In a crystal network of the newly proposed 2D carbon allotrope, each carbon atom binds with its four neighboring atoms in a 2D double layer crystal structure, which is formed by hexagonal prisms. The meta-stability of our newly designed monolayer is shown by the cohesive energy and phonon calculations. Furthermore, evaluating the electrical properties of hexatetra-carbon shows that it exhibits a semiconductor behavior with a moderate band gap.

The current paper is organized as follows: The computational method used in our study is presented in Section 2. In Section 3, the details of the structural properties are outlined and the stability of the proposed structure is discussed. The electrical characteristics of the designed hexatetra-carbon monolayer and its potential applications are described in Section 4. Finally, in Section 5, the paper is concluded.

#### **2. Computational Methods**

To obtain accurate structural and electrical properties of the proposed 2D material, the full-potential linearized-augmented plane wave (FP-LAPW) scheme was utilized [21]. This scheme is based on the DFT implemented in the WIEN2k computational package [22], in which the generalized gradient approximation (GGA) parameterized by Perdew–Burke–Ernzerhof (PBE) was used [23]. Moreover, since the GGA method underestimates band gaps, to obtain more reliable band gaps, the screened short range hybrid functional exchange correlation implemented in the WIEN2k code [24] was employed for band gap calculations. Moreover, the Monkhorst–Pack scheme [25] was used to sample the Brillouin zone by a 24 × 24 × 1 k-mesh, where we chose RKmax = 7, Gmax = 14 Ry1/2, and lmax = 10. To avoid interlayer interactions, a large vacuum distance of 20 Å along the non-periodic direction was utilized. With regards to dynamic stability, an evaluation of all the calculations was conducted with the Quantum Espresso (QE) package [26]. Furthermore, the Martin–Troullier norm-conserving pseudopotential [27] was used to treat the core electrons, while the valence electronic wave functions were expanded using an energy cut-off of 80 Ry. However, for the investigation of structural properties of bilayer graphene, we considered the van der Waals correction in our calculations.

#### **3. Structural Properties and the Stability of Hexatetra-Carbon Monolayer**

The design of the 2D hexatetra-carbon monolayer was initiated by examining the effect of vertical compressive strains on the structural, electronic, and optical properties of bilayer graphene in its AA-stacking configuration. Recently, several research works have reported the influence of strain on the different physical characteristics of bilayer graphene.

Moreover, it is well-known that a crystal structure of a bilayer graphene stabilizes at the interlayer distance of about *h*0 = 3.4 Å layer separation, in which the long-range van der Waals (vdW) interaction plays an important role in the structural properties of the materials. In this case, when evaluating the effect of vertical compressive strains on the different properties of AA-stacked bilayer graphene, with the aim of preserving the stability of the initial AA-stacked bilayer graphene, researchers mostly consider vertical strains that influence a maximum variation of 20% in interlayer separation, i.e., *h* = *h*0 ± 20% *h*0 (in which the long-ranged van der Waals (vdW) interaction should be taken into account). In these conditions, not only the stability of the AA-stacked bilayer graphene is preserved, but also its crystal structure configuration is retained. However, one may be curious to know the effect of higher values of compressive vertical strains, in which the layer separation distance of h < 2.0 Å for the long-range van der Waals (vdW) interaction has no important effect. On this basis, we found a new unprecedented 2D allotrope of carbon called hexatetracarbon during our structure searches. Figure 1 shows the optimized structure of our designed hexatetra-carbon monolayer from different views. A unit cell of hexatetra-carbon monolayer consists of four carbon atoms with optimized lattice constants of a = b = 2.65 Å. In a crystal network of the hexatetra-carbon, C atoms are distributed in two different exactly planar atomic planes with a vertical distance of about 1.55 Å. Specifically, each C atom of the hexatetra-carbon bonds to four neighboring C atoms to form a tetra-coordinated carbon structure, in which the C–C bond lengths are about 1.55 Å, longer than those in graphene (1.43 Å) and almost equal to the C1–C2 bond lengths in penta-graphene [13].

**Figure 1.** Crystal structure of 2D hexatetra-carbon from different views.

Furthermore, before an evaluation could be performed for the electrical properties and potential applications of the hexatetra-carbon monolayer, it is imperative to analyze the stability of the monolayer structure. This was carried out by first evaluating its energetic stability through a calculation of its cohesive energy, which is given by *Ecoh* = *Ehexatetra*−*carbon*−4*Ec* 4 (*Ec*, *Ehexatetra*−*carbon* are the total energies of a single C atom, and a unit cell of the hexatetracarbon monolayer, respectively). Based on our calculation, the 2D hexatetra-carbon monolayer has a cohesive energy of about −6.86 eV/atom. For comparison, we also calculated the cohesive energy of graphene, penta-graphene, and T-carbon [28], which is a 3D carbon allotrope obtained by replacing each atom in diamond with a carbon tetrahedron at the same theoretical level. The values were about −8.01, −7.02, −6.34 eV/atom, respectively, which confirm that the proposed 2D hexatetra-carbon monolayer shows good energetic stability. Figure 2 shows a unit cell relative energy of the predicted hexatetra-carbon monolayer compared with different carbon allotropes under strain conditions.

**Figure 2.** Cohesive energy of different carbon allotropes.

Next, we examined the dynamic stability of the 2D hexatetra-carbon monolayer by calculating its phonon dispersion spectrum. As shown in Figure 3, there are no imaginary phonon modes in the whole Brillouin zone, indicating that the hexatetra-carbon is a local minimum on the potential energy surface and can be considered as a metastable allotrope of carbon compared with the other carbon allotropes. Specifically, we obtained the highest frequency of 1267 cm−1, which is higher than those obtained for silicon [29] (580 cm−1), MoS2 monolayer [30] (473 cm−1), and TiC monolayer [31] (810 cm−1). However, it is lower than the highest phonon frequency of graphene (about 1650 cm−1) and penta-graphene (about 1600 cm−1), indicating robust C–C bonds in the predicted monolayer.

**Figure 3.** Phonon dispersion spectrum of 2D hexatetra-carbon monolayer.

From the above indications that the 2D hexatetra-carbon may be meta-stable, we next systematically study its mechanical properties. Using Young's modulus (*Y*), Equation (1) can be calculated [32], where ε is the strain in the vicinity of the optimum lattice vector (ε = (*a* − a0)/a0) and *E* is the total energy. In addition, *V*0 is the equilibrium volume of the 2D material evaluated by *V*0 = 3 √3*d*2*h* 2 , where *d* is the adjacent carbon distance in the hexagon ring in the *a* and *b* plane, and *h* is the thickness of the 2D material along the *c* vector (Figure 1).

$$Y = \frac{1}{V\_0} \left(\frac{\partial^2 E}{\partial \varepsilon^2}\right)\_{\varepsilon = 0} \tag{1}$$

In this work, the estimated Young's modulus value for hexatetra-carbon is 1859.7 GPa. It is smaller than those estimated by Raman spectroscopy [33] for the single- and bilayer graphene, 2400 ± 400 GPa and 2000 ± 500 GPa, respectively. With regards to carbon fibers, the Young's modulus value is 235–427 GPa [34]. Here, the value of hexatetra-graphene is greater than those obtained for other materials, such as imogolites (∼320–370 GPa) [35] and GaS (∼270 GPa) [36], as well as the nanotubes and MoS2 monolayer (265 ± 13 GPa), which are similar to the bilayer [37].

#### **4. Electronic Properties**

To analyze the electrical properties of the 2D hexatetra-carbon monolayer, we calculated the band structure of the designed 2D monolayer. As shown in Figure 4, an indirect bandgap of about 2.20 eV, which is calculated using the hybrid functional level of theory can be seen in the band structure. Moreover, the valence band maximum (VBM) of this 2D monolayer material is located between K and Γ points, while the conduction band minimum (CBM) is located at the Γ point. Therefore, the 2D hexatetra-carbon monolayer, which is a semiconductor with a moderate band gap, is not similar to the monolayer and bilayer graphene due to the fact that they are both semi-metals.

**Figure 4.** Band structures of (**a**) graphene, (**b**) bilayer graphene, (**c**) penta-graphene, and (**d**) hexatetra-carbon.

To gain a deeper insight into the bonding nature of hexatetra-carbon and its structural analogy with respect to AA-stacked graphene, we investigated the 2D valence charge density distribution of these two nanostructures (see Figure 5d–g). The crystal structure of AA-stacked bilayer graphene, the cubane molecule, and hexatetra-carbon are shown in Figure 5a–c. It is clear that the hexatetra-carbon structure has a hexagonal face similar to graphene (xy plane). In addition, interlayer bonds along the z direction between the two graphene monolayers show a tetragonal side view, which is similar to the cubane molecule.

As shown in Figure 5d, there are three sp2 sigma covalent bonds (1.43 Å) between the C–C atoms in each monolayer of AA-stacked bilayer graphene, as well as a weak pz-pz interaction between the monolayers. In comparison, the in-plane bond length of C–C for hexatetra-carbon is 1.56 Å, which is longer than those obtained in AA-stacked bilayer graphene, i.e., the C–C in-plane orbital overlap decreases for the hexatetra-carbon (see Figure 5f). Therefore, to retain its structural stability, the hexatetra-carbon nanostructure compensates this orbital variation by creating interlayer sigma bonds between the neighboring carbons, which are located in the different planes. Moreover, these interlayer bonds would restrict pz, resulting in the semiconducting nature of 2D hexatetra-carbon.

Furthermore, we calculated the variation of the cohesive energy of two graphene sheets as they approach each other. As seen in Figure 6, when a graphene sheet moves towards another fixed graphene sheet, they tend to stabilize at the vertical distance of about 3.45 Å and form the AA-stacked bilayer graphene structure, which is the most stable structure in this situation. However, by applying an additional external vertical mechanical strain to the AA-stacked bilayer graphene through decreasing its interlayer distance, its lattice parameter increases. In addition, another local energy minimum occurs for the two graphene sheets with a vertical interlayer distance of about 1.55 Å. In this case, the hexatetra-carbon includes four carbon atoms in its unit cell, and the lattice parameter of about 2.65 Å structure is formed. In other words, by applying the external vertical strain, a type of phase transition occurs.

**Figure 5.** Crystal structure of (**a**) AA-stacked bilayer graphene, (**b**) cubane molecule, and (**c**) hexatetracarbon (**d**,**e**). Top and side views of valence charge density distribution for AA-stacked graphene and (**f**,**g**) for hexatetra-carbon obtained by WIEN2K code [21].

**Figure 6.** Cohesive energy of AA-stacked bilayer graphene and hexatetra-carbon versus interlayer distance calculated by Quantum Espresso (due to the interruption in accessing WIEN2k while following up on a reviewer's comment, we have used Quantum Espresso for Figure 6).

#### **5. Summary**

In conclusion, utilizing first-principles calculations based on DFT, we proposed a new 2D monolayer crystalline allotrope of carbon called hexatetra-carbon, which shows a semiconducting nature with a moderate band gap. The stability of 2D hexatetra-carbon monolayer is confirmed by its rather high cohesive energy, as well as the absence of imaginary phonon modes in its phonon dispersion spectrum. Comparing the cohesive energy of the hexatetra-carbon monolayer with those obtained for other 2D allotropes confirms its high bonding properties. Moreover, our analysis on the charge densities of both bilayer graphene and hexatetra-carbon monolayer indicates that by applying a vertical strain on bilayer graphene, the hexatetra-carbon structure can be obtained. Therefore, one may be optimistic that the 2D hexatetra-carbon monolayer can be experimentally achieved in the foreseeable future.

**Author Contributions:** All authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Natural Sciences and Engineering Research Council of Canada Discovery Grant number RGPIN-2019-03976. Computational resources were provided by Compute Canada/WestGrid.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Article* **Length-Gauge Optical Matrix Elements in WIEN2k**

**Oleg Rubel 1,\* and Peter Blaha 2,\***

**Abstract:** Hybrid exchange-correlation functionals provide superior electronic structure and optical properties of semiconductors or insulators as compared to semilocal exchange-correlation potentials due to admixing a portion of the non-local exact exchange potential from a Hartree–Fock theory. Since the non-local potential does not commute with the position operator, the momentum matrix elements do not fully capture the oscillator strength, while the length-gauge velocity matrix elements do. So far, length-gauge velocity matrix elements were not accessible in the all-electron full-potential WIEN2k package. We demonstrate the feasibility of computing length-gauge matrix elements in WIEN2k for a hybrid exchange-correlation functional based on a finite difference approach. To illustrate the implementation we determined matrix elements for optical transitions between the conduction and valence bands in GaAs, GaN, (CH3NH3)PbI3 and a monolayer MoS2. The non-locality of the Hartree– Fock exact exchange potential leads to a strong enhancement of the oscillator strength as noticed recently in calculations employing pseudopotentials (Laurien and Rubel: arXiv:2111.14772 (2021)). We obtained an analytical expression for the enhancement factor for the difference in eigenvalues not captured by the kinetic energy. It is expected that these results can also be extended to other non-local potentials, e.g., a many-body *GW* approximation.

**Keywords:** semiconductors; oscillator strength; density functional theory; hybrid exchange-correlation functional; non-local potential

This paper is dedicated to the 80th birthday of Professor Karlheinz Schwarz, the founder of the WIEN2k DFT package.

#### **1. Introduction**

Calculations of linear optical properties of solids require matrix elements for electric dipole transitions. Momentum matrix elements

$$\mathbf{p}\_{mn}(\mathbf{k}) = \langle m\_\prime \mathbf{k} | -i\nabla\_\mathbf{r} | n\_\prime \mathbf{k} \rangle \tag{1}$$

are widely used in full-potential codes with periodic boundary conditions [1] when optical properties are computed with local potentials (e.g., LDA (see end of the paper for the full list of abbreviations) or GGA XC functionals) and referred to in the literature as a velocity gauges. (Atomic units will be used throughout the paper.)

Starace [2] emphasised the limitations of Equation (1) when representing matrix elements for electric dipole transitions. Instead, the more general velocity matrix elements should be used

$$\mathbf{v}\_{mn}(\mathbf{k}) = \langle m, \mathbf{k} | i[\hat{H}, \mathbf{r}] | n, \mathbf{k} \rangle \tag{2}$$

with the velocity operator [2]

$$\hat{\mathbf{v}} = i[\hat{H}, \mathbf{r}] = \hat{\mathbf{p}} + i[\hat{V}\_{\text{NL}}(\mathbf{r}, \mathbf{r}'), \mathbf{r}]\_{\prime} \tag{3}$$

**Citation:** Rubel, O.; Blaha, P. Length-Gauge Optical Matrix Elements in WIEN2k. *Computation* **2022**, *10*, 22. https://doi.org/ 10.3390/computation10020022

Academic Editor: Henry Chermette

Received: 29 December 2021 Accepted: 21 January 2022 Published: 26 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

which contains an additional commutator term [*V*ˆ NL(**r**,**r** ),**r**] to account for non-local potentials (e.g., the Hartree-Fock exchange in hybrid XC functionals). With |*n*, **k** and *En*(**k**) being eigenstates of the Hamiltonian *H*ˆ , the alternative (length gauge) matrix elements can be expressed as [2]

$$
\langle m, \mathbf{k} | i[\hat{H}, \mathbf{r}] | n, \mathbf{k} \rangle = i[E\_m(\mathbf{k}) - E\_n(\mathbf{k})] \langle m, \mathbf{k} | \mathbf{r} | n, \mathbf{k} \rangle. \tag{4}
$$

Since the position operator is not well defined for periodic systems, the following substitution is used instead: [3,4]

$$\mathbf{r} = \lim\_{q \to 0} (e^{iq\mathbf{r}} - 1) / iq,\tag{5}$$

which leads to a practical expression for the velocity matrix elements in the long wavelength limit [5,6]

$$\psi\_{mn}^{(a)}(\mathbf{k}) = \lim\_{q \to 0} \frac{1}{q} \langle m, \mathbf{k} + \mathbf{q}\_a | \mathbf{e}^{i\mathbf{q}\_a \cdot \mathbf{r}} | n, \mathbf{k} \rangle [E\_m(\mathbf{k} + \mathbf{q}\_a) - E\_n(\mathbf{k})]. \tag{6}$$

Here, *α* = *x*, *y*, *z* is a Cartesian direction, *m* and *n* are band indices, and **q***x* = *q* **x**ˆ, where **x**ˆ is a unit vector in the direction of the *x* axis.

WIEN2k [7,8] is one of the most used full-potential all-electron DFT codes for solids. It offers many XC functionals to open the band gap, including hybrids with a non-local Hartree-Fock potential [9]. So far, however, WIEN2k has implemented only momentum matrix elements to compute optical properties as a part of the optic module [1]. Laurien and Rubel [10] showed that neglecting the second term in Equation (3) when using hybrid functionals can lead to an underestimation of the squared magnitude of matrix elements for electric dipole transitions between conduction and valence band edges by ca. 30%.

Here, we present a scheme for the calculation of the length-gauge optical matrix elements in WIEN2k based on a finite difference Equation (6) with the help of overlap matrix elements

$$M\_{mn}(\mathbf{k}, \mathbf{q}) = \langle \mu\_{\mathbf{k}, n} | \mu\_{\mathbf{k} + \mathbf{q}, m} \rangle \equiv \langle \psi\_{\mathbf{k}, n} | \mathbf{e}^{-i\mathbf{q} \cdot \mathbf{r}} | \psi\_{\mathbf{k} + \mathbf{q}, m} \rangle \tag{7}$$

that come from the wien2wannier module [11]. This development opens an avenue for the calculation of optical properties (frequency-dependent dielectric tensor, absorption spectrum, optical conductivity, refractive index, reflectively, loss function) in the independent particle approximation with hybrid functionals in WIEN2k.

#### **2. Methods**

DFT [12,13] calculations were performed with the WIEN2k package (version 21.1) and the Yukawa screened hybrid (YSH) functional [9]. It was shown that by choosing an appropriate screening length *λ* in the Yukawa potential the YSH functional gives very similar results as the common HSE06 XC functional [14,15]. Important parameters are summarized in Table 1. Experimental structural parameters were used for all solids (Figure 1) with internal atomic positions optimized at the PBE level when permitted by symmetry. Spin-orbit coupling (SOC) was included in all calculations. The structure of (CH3NH3)PbI3 was represented by a pseudo-cubic cell taken from Ref. [16], scaled to experimental lattice parameters at 350 K [17,18], followed by a subsequent relaxation of atomic positions while retaining the experimental lattice parameters. The pseudo-cubic structure means that the following constrains *a* = *b* = *c*, *α* = *β* = *γ* = 90◦ are applied to lattice parameters, while the formal symmetry of the structure (spacegroup P1) is not cubic.

**Table 1.** Structural and calculation parameters.

**Figure 1.** Crystal structures: (**a**) GaAs, (**b**) GaN, (**c**) quasi-cubic (CH3NH3)PbI3, and (**d**) monolayer MoS2.

Velocity-gauge optical matrix elements **p***mn*(**k**) were calculated using the optic module [1] in WIEN2k. Length-gauge optical matrix elements **v***mn*(**k**) were obtained with the forward

$$
\psi\_{mn}^{(\boldsymbol{\alpha})}(\mathbf{k}) \approx \frac{1}{q} \langle u\_{\mathbf{k},n} | u\_{\mathbf{k}+\mathbf{q},m} \rangle [E\_m(\mathbf{k}+\mathbf{q}\_n) - E\_n(\mathbf{k})] \tag{8}
$$

and central

$$v\_{mn}^{(a)}(\mathbf{k}) \approx \frac{1}{q} \langle u\_{\mathbf{k} - \frac{1}{2}\mathbf{q}\_{\mu,n}} | u\_{\mathbf{k} + \frac{1}{2}\mathbf{q}\_{\mu,m}} \rangle [E\_m(\mathbf{k} + \frac{1}{2}\mathbf{q}\_a) - E\_n(\mathbf{k} - \frac{1}{2}\mathbf{q}\_a)] \tag{9}$$

finite difference method. The overlap matrix elements *u***k**,*n*|*u***k**,*m* between the cell-periodic parts of the Bloch functions were generated by the wien2wannier module [11] (case.mmn output file). The length-gauge optical matrix elements for GaAs computed with YSH were additionally verified using VASP [22], HSE06 and projector augmented-wave potentials [22–24]. Sample scripts that illustrate a detailed workflow can be found in the Supporting Information section.

The logarithmic percent change

$$\Delta = \ln \left( \frac{\sum |v\_{mn}^{(\alpha)}|^2}{\sum |p\_{mn}^{(\alpha)}|^2} \right) 100\% \tag{10}$$

was used to evaluate differences between matrix elements. This approach has the following advantages: (i) does not require a reference, (ii) is more suitable for large changes (greater than a few percents), (iii) it has additive properties, and (iv) in the limit of small changes it reduces to the classical ratio of the relative change to the reference.

#### **3. Results**

#### *3.1. Finite Difference Calibration and Validation for Local XC*

We selected GaAs and the local PBE XC potential [25] to prove the feasibility of computing the length-gauge optical matrix elements using the finite difference methods given by Equations (8) and (9). The local potential is selected here since both the lengthand velocity-gauge should lead to identical results under these circumstances. It is also important to get a feeling for the step size *q* at which the finite difference approximation converges to the accurate result given by the momentum matrix element. Here GaAs serves as an important benchmark since the conduction band and light holes are very sharp and non-parabolic (see Figure 2a).

Results presented in Table 2 suggest that both, the forward and the central finite differences reproduce the values of the momentum matrix element within a 3% error. The central finite difference converges faster (at the wave vector shift of *q* ≈ 0.003 rad bohr−1) and will be used to derive matrix elements for the remaining part of this paper. The numerical noise of the finite difference starts to show up at *q* < 10−6 rad bohr−1.

**Table 2.** Length-gauge velocity matrix elements |*v* (*x*) *mn*| 2 (at.u.) in GaAs calculated using the finite difference approximations (forward vs. central) Equations (8) and (9) with various step sizes *q*. These values are compared with the velocity-gauge momentum matrix elements |*p* (*x*) *mn*| 2 (at.u.) from the optic module. The local (GGA-PBE) XC functional wa used, which made the velocity and the length gauges identical. The band degeneracy is given as a superscript in parentheses and the meaning of the subscripts is made clear in Figure 2a.

#### *3.2. Validation for Non-Local XC*

After validating our approach with the local potential, we applied it to the non-local YSH XC functional. Again we evaluated the velocity-gauge (momentum) matrix elements and length-gauge (velocity) matrix elements in GaAs. Now we did not expect the two matrix elements to agree given the arguments presented in Section 1. To cross-check our results, we also computed the velocity matrix elements with VASP, which should be comparable with our *v* (*α*) *mn* values.

Our YSH calculations for GaAs gave a band gap of *E*YSH g = 1.24 eV vs. *E*exp g = 1.52 eV [26] and previously reported *E*HSE06 g = 1.33 eV [27], which was a significant improvement over PBE. The band structure is shown schematically in Figure 2a where bands are labeled according to the convention. The results presented in Table 3 confirmed the agreement between WIEN2k and VASP for length-gauge matrix elements within less than a 2% deviation. The total length-gauge oscillator strength between valence and conduction band corresponded to *m*0 ∑ |*v* (*x*) *cv* | 2 ≈ 21 eV, which agreed well with the 20 eV quoted by (Yu and Cardona [28], Section 2.6) for III-V semiconductors. The momentum (velocity-gauge) matrix elements significantly underestimated the strength of optical transitions, which had previously been reported and quantified in Ref. [10]. The values of |*pcv*| 2 were almost identical to those obtained with PBE (Table 2), even though the momentum matrix elements were derived from YSH wave functions.

**Figure 2.** Schematic band structures of materials studied with SOC: (**a**) GaAs, (**b**) GaN, (**c**) quasicubic (CH3NH3)PbI3, and (**d**) monolayer MoS2. The band indices 'c', 'hh', 'lh', and 'so' stand for conduction, valence heavy-hole, light-hole, and split-off bands, respectively. The scale of band splittings is exaggerated.

**Table 3.** Length-gauge |*v* (*x*) *mn*| 2 and velocity-gauge |*p* (*x*) *mn*| 2 matrix elements (at.u.) in GaAs calculated using WIEN2k (with YSH) and VASP (with HSE06). Due to the non-local potential the velocity and the length gauges are *not* identical. The band degeneracy is given as a superscript in parentheses and the subscripts are explained in Figure 2a. The logarithmic deviation between ∑ |*p* (*x*) *vc* | 2 and ∑ |*v* (*x*) *vc* | 2 is given in parentheses (Δ as per Equation (10)).

#### *3.3. Illustrative Applications*

In the previous subsection we showed that calculations of optical properties for GaAs with the non-local hybrid XC functional (YSH or HSE06) require length-gauge optical matrix elements. If the momentum matrix elements had been used instead, the strength of optical transitions would have been underestimated by 23%. Next, we showed that a similar enhancement of the strength of direct optical transitions was also observed in other semiconductors, such as GaN, (CH3NH3)PbI3, and monolayer MoS2. The corresponding band structures are shown schematically in Figure 2b–d.

The band gap of GaN is well reproduced with YSH: *E*YSH g = 3.19 eV vs *E*exp g = 3.30 eV [29]. Compared to GaAs, the optical matrix elements in Table 4 showed an even larger disparity between the length-gauge velocity and the momentum matrix elements.

**Table 4.** Length-gauge |*v* (*α*) *mn*| 2 and velocity-gauge |*p* (*α*) *mn*| 2 matrix elements (at.u.) in GaN calculated using the YSH XC functional.

The monolayer MoS2 had a direct band gap at the *K* = (1/3, 1/3, 1/3) point. Due to large excitonic effects [30], direct comparison of the YSH band gap *E*YSH g = 2.22 eV with experiment was not possible. Thus, we used a many-body result *EG*0*W*0 g = 2.53 eV [31,32] as a reference. Similarly to other materials, the monolayer MoS2 showed strong enhancement of the matrix elements (Table 5) with the YSH XC functional. Spin selection rules disabled half of the in-plane *v* (*x*) *cv* matrix elements, while the out-of-plane matrix elements *v* (*z*) *cv* were zero for transitions at the band edges due to symmetry arguments.

**Table 5.** Length-gauge |*v* (*x*) *mn*| 2 and velocity-gauge |*p* (*x*) *mn*| 2 matrix elements (at.u.) in monolayer MoS2 calculated using YSH XC functional.

The pseudo-cubic (CH3NH3)PbI3 had a direct gap at *R* = (1/2, 1/2, 1/2). The calculated band gap *E*YSH g = 1.03 eV is an improvement relative to the PBE band gap (0.46 eV), but it is still far from the experimental 1.5–1.6 eV for the tetragonal phase [33,34]. This underestimation was due to the lack of stochastic thermal distortions of the PbI6 octahedra, which further opened the gap by ca. 0.5 eV at room temperature [35,36]. Table 6 captures the matrix elements and their anisotropy caused by the reduced (pseudo-cubic) symmetry of the unit cell. Among all materials studied here, this material showed the lowest enhancement of the velocity matrix elements compared to the momentum matrix elements.

**Table 6.** Length-gauge |*v* (*α*) *mn*| 2 and velocity-gauge |*p* (*α*) *mn*| 2 matrix elements (at.u.) in pseudo-cubic (CH3NH3)PbI3 calculated using the YSH XC functional.

*a* The true enhancement should be about 22% due to inaccuracies in |*p* (*x*,*y*,*z*) *mn* | 2 values. See text below for more details.

It should be mentioned that the momentum matrix elements calculated with the optic module in the presence of SOC had an inaccuracy that progressively increased for heavier elements. The discrepancy between |*p* (*x*,*y*,*z*) *mn* | 2 values calculated at the PBE level (including SOC) with the optic module and using the finite difference overlap matrix reached ca. 12% in the case of (CH3NH3)PbI3. The discrepancy fully vanished when SOC was excluded. After crosschecking the matrix elements with VASP we concluded that the finite difference results were correct. Since the optic module overestimated |*p* (*x*,*y*,*z*) *mn* | 2 values at PBE with SOC, the same applied to YSH with the SOC results presented in Table 6. After including this error, the true enhancement of YSH matrix elements for (CH3NH3)PbI3 were about 22% (10% average enhancement in Table 6 and 12% optic error for this material). Additional calculations with VASP and HSE06 XC functional with SOC produced a very similar result (23% enhancement of the matrix elements).

#### **4. Discussion**

YSH length-gauge |*v*YSH *m*,*n* | 2 matrix elements were systematically greater than the momentum matrix elements |*p*YSH *m*,*n* | 2. The enhancement ranged from 22 to 36% in the following order: (CH3NH3)PbI3, GaAs, GaN, and MoS2 (from the smaller to higher enhancement). This trend prompted the hypothesis that the enhancement was related to the localization of states involved in the optical transition. (CH3NH3)PbI3 has the most extended 5p-I and 6p-Pb states, while MoS2 had the most localized 4d-Mo and 3p-S states at the band edges.

To gain further insight into the difference between |*p*YSH *m*,*n* | 2 and |*v*YSH *m*,*n* | 2 we wote the momentum matrix element in the length gauge. The corresponding operator was expressed as the commutator

$$
\hat{\mathbf{p}} = \mathbf{i}[\hat{T}, \mathbf{r}]\_\prime \tag{11}
$$

where *T*ˆ is the kinetic energy operator. Following the same logic that leads to Equation (4), we derived an equivalent expression for the momentum matrix element in the length-gauge

$$\mathbf{p}\_{m,n}(\mathbf{k}) = i[T\_m(\mathbf{k}) - T\_n(\mathbf{k})] \langle m, \mathbf{k} | \mathbf{r} | n, \mathbf{k} \rangle. \tag{12}$$

After dividing Equation (4) by (12) we obtain

$$\mathbf{v}\_{m,n}(\mathbf{k}) = \mathbf{p}\_{m,n}(\mathbf{k}) \frac{E\_m(\mathbf{k}) - E\_n(\mathbf{k})}{T\_m(\mathbf{k}) - T\_n(\mathbf{k})}.\tag{13}$$

Thus, the 10 to 36% enhancement of the absolute squared magnitude of velocity matrix elements vs momentum matrix elements in calculations with YSH was directly related to the difference in eigenvalues not captured by the kinetic energy. In contrast, we expected the difference in eigenvalues to be fully captured by the kinetic energy, i.e., [*Em*(**k**) − *En*(**k**)]/[*Tm*(**k**) − *Tn*(**k**)] = 1, when a local potential was employed. Interestingly, Equation (13) predicts an *isotropic* renormalization factor shared by all Cartesian directions (*α* = *x*, *y*, *z*). Indeed, materials with anisotropic *v* (*α*) *mn*—GaN (Table 4) and (CH3NH3)PbI3 (Table 6)—showed a material-dependent yet isotropic enhancement factor, which was an indirect proof of the validity of Equation (13).

The renormalization of momentum matrix elements should have had implications for optical properties calculated with non-local potentials (hybrid or quasi-particle *GW*). The velocity matrix elements entered the frequency-dependent dielectric tensor (an imaginary part of the inter-band contribution) that took the following form in the independent-particle approximation [5]

$$
\varepsilon\_{a\beta}^{\prime\prime}(\omega) \propto \sum\_{\upsilon,\varepsilon} \int\_{\mathbf{k} \in \mathcal{BZ}} \frac{v\_{\upsilon,\varepsilon}^{(a)}(\mathbf{k}) v\_{\varepsilon,\upsilon}^{(\beta)}(\mathbf{k})}{\omega^2} \delta[E\_c(\mathbf{k}) - E\_\upsilon(\mathbf{k}) - \omega] \,d\mathbf{k}.\tag{14}
$$

However, length-gauge **v***c*,*v* matrix elements are more difficult to compute than **p***c*,*v*, especially at the quasi-particle *GW* level of theory where the finite difference method seems the only available technique [3]. Equation (13) opens a convenient possibility to use renormalized momentum matrix elements instead

$$
\epsilon\_{a\beta}^{\prime\prime}(\omega) \approx \sum\_{\upsilon,\varepsilon} \int\_{\mathbf{k} \in \mathcal{BZ}} \frac{p\_{\upsilon,\varepsilon}^{(a)}(\mathbf{k}) p\_{\varepsilon,\upsilon}^{(\delta)}(\mathbf{k})}{[T\_{\varepsilon}(\mathbf{k}) - T\_{\upsilon}(\mathbf{k})]^2} \delta[E\_{\varepsilon}(\mathbf{k}) - E\_{\upsilon}(\mathbf{k}) - \omega] \, d\mathbf{k},\tag{15}
$$

provided that eigenstates are consistent with the potential, and their kinetic energy is known. The last expression should be valid not only for hybrid XC functionals but also for the quasi-particle *GW* level of theory.

Finally, we would like to comment on the renormalization of optical transition matrix elements proposed by Levine and Allan [37]

$$v\_{v, \mathcal{c}}^{GW} = v\_{v, \mathcal{c}}^{\text{LDA/PBE}} \frac{(E\_{\mathcal{c}} - E\_v)^{GW}}{(E\_{\mathcal{c}} - E\_v)^{\text{LDA/PBE}}} \tag{16}$$

that is further used in the literature [5,38]. If we apply Equation (16) to the Γ*so* − Γ*c* transition in GaAs, one would expect the absolute squared magnitude of the velocity matrix element to increase by the ratio of [(*E*YSH *c* − *E*YSH *so* )/(*E*PBE *c* − *E*PBE *so* )]2 which amounts to +151%. This result contradicts the +23% difference between |*p*PBE *so*,*c* | 2 and |*v*YSH *so*,*c* | 2 we observed (compare Tables 2 and 3). At the same time, the dipole matrix element |Γ*so*|**r**|Γ*c*|2 becomes 127% smaller in YSH relative to PBE and counterbalances (partly) the effect of the gap opening. We further identified contributions of the muffin-tin spheres and of the interstitial volume to the value of the dipole matrix element Γ*so*|**r**|Γ*c* at PBE and YSH levels of theory: 25% Ga, 44% As, and 31% interstitial. All contributions are in phase with each other, and the proportions remain unchanged from PBE to YSH. Equation (16), in contrast, implies the equality of dipole matrix elements *m*, **k**|**r**|*n*, **k**LDA/PBE = *m*, **k**|**r**|*n*, **k***GW* (see Equation (4). Note that *v*LDA/PBE *m*,*n* = *p*LDA/PBE *m*,*n* ≈ *pGW m*,*n* [10]) leads to a gross overestimation of *vGW m*,*n* matrix elements making them inconsistent with the band curvature [10].

#### **5. Conclusions**

Strong material-dependent enhancement of the oscillator strength (22–36% in the absolute squared magnitude) was observed in the electronic structure calculations of semiconductors with a hybrid XC functional. The origin of the enhancement was traced to the non-local Hartee–Fock exchange potential. The enhancement of the absolute squared magnitude of velocity matrix elements |*v*YSH *m*,*n* | 2 vs momentum matrix elements |*p*YSH *m*,*n* | 2 in calculations with non-local potentials was directly related to the difference in eigenvalues not captured by the kinetic energy, i.e., [*Em*(**k**) − *En*(**k**)]2/[*Tm*(**k**) − *Tn*(**k**)]2. This enhancement is isotropic and can be readily included in a calculation of the dielectric function. Our

enhancement factor was much more accurate than that previously proposed by Levine and Allan (*EGW* g /*E*LDA g )2, which leads to nonphysically large *vGW m*,*n* matrix elements.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/computation10020022/s1, a SI-Computation.tar.gz file with WIEN2k workflows (bash scripts), structure files, sample k-point files, and a python script for data processing.

**Author Contributions:** Conceptualization, O.R. and P.B.; investigation, O.R.; validation, P.B. and O.R.; writing—original draft, O.R.; writing—review and editing, O.R. and P.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** O.R. acknowledges funding provided by the Natural Sciences and Engineering Research Council of Canada under the Discovery Grant Programs RGPIN-2020-04788. Calculations were performed using the Compute Canada infrastructure supported by the Canada Foundation for Innovation under John R. Evans Leaders Fund.

**Institutional Review Board Statement:** Not applicable.

**Acknowledgments:** Authors are thankful to Robert Laskowski (A\*STAR, Singapore) and Fabien Tran (VASP Software GmbH, Vienna) for helpful discussions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **References**

## *Article* **Performance Enhancement of APW+***lo* **Calculations by Simplest Separation of Concerns**

**Long Zhang 1,2,3, Anton Kozhevnikov 4, Thomas Schulthess 4, Hai-Ping Cheng 1,2,3 and Samuel B. Trickey 1,2,3,\***

	-<EMAIL_ADDRESS>(T.S.)

**Abstract:** Full-potential linearized augmented plane wave (LAPW) and APW plus local orbital (APW+*lo*) codes differ widely in both their user interfaces and in capabilities for calculations and analysis beyond their common central task of all-electron solution of the Kohn–Sham equations. However, that common central task opens a possible route to performance enhancement, namely to offload the basic LAPW/APW+*lo* algorithms to a library optimized purely for that purpose. To explore that opportunity, we have interfaced the Exciting-Plus ("EP") LAPW/APW+*lo* DFT code with the highly optimized SIRIUS multi-functional DFT package. This simplest realization of the separation of concerns approach yields substantial performance over the base EP code via additional task parallelism *without* significant change in the EP source code or user interface. We provide benchmarks of the interfaced code against the original EP using small bulk systems, and demonstrate performance on a spin-crossover molecule and magnetic molecule that are of size and complexity at the margins of the capability of the EP code itself.

**Keywords:** LAPW method; APW+*lo* method; all-electron DFT

#### **1. Dedication**

Much credit for the widespread use of full-potential linearized augmented plane wave (LAPW) methodology to solve the Kohn–Sham (KS) [1] equation for solids goes to Karlheinz Schwarz. The history of that contribution is related in Section 3 of Ref. [2]. It suffices to say here that Heinz started using the original APW in his thesis work, then came to Gainesville and Quantum Theory Project (QTP) in 1969 to work with Prof. J.C. Slater, the inventor of the APW method. That is how the last author of this paper became a collaborator with Heinz and a friend.

Years later, when linearized methods removed the explicit energy dependence difficulty in the APW basis, Heinz undertook development of the code that became WIEN [3]. Again there was a collaboration involving QTP, University of Florida, and the last author. Apparently that was the first publicly available FLAPW code. By now, it has evolved to WIEN2k [2,4]. During that evolutionary period, methodological developments led to revival of the APW scheme via the APW plus local orbitals (APW+*lo*) combination (summarized below). That has been instantiated in several other codes as well as WIEN2k, e.g., ELK [5], FlEUR [6], *exciting* [7], and Exciting-plus [8]. Here, we are pleased to contribute to the further advance of this important methodology and particularly pleased to be able to do so in honor of Heinz' birthday.

**Citation:** Zhang, L.; Kozhevnikov, A.; Schulthess, T.; Cheng, H.-P.; Trickey, S.B. Performance Enhancement of APW+*lo* Calculations by Simplest Separation of Concerns. *Computation* **2022**, *10*, 43. https://doi.org/ 10.3390/computation10030043

Academic Editor: Henry Chermette

Received: 7 February 2022 Accepted: 4 March 2022 Published: 16 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **2. Motivating Physical Systems**

The materials physics problem class driving our effort is molecular magnetism, in particular, the contriving of condensed aggregates of molecular magnets into materials of relevance to quantum information systems, notably, quantum computing [9,10]. The computational challenge is to predict promising molecular magnets [11,12], and promising aggregates of them as well as parametrizing spin Hamiltonians and aiding interpretation of experimental data. The molecules themselves are large and complicated.

A pertinent example is the molecular magnet NiCl2-[SC(NH2)2]4, dichloro-tetrakisthourea-nickel, commonly called "DTN" and its Co analogue, DTC [13]. DTN is important in this context because of its multi-ferroicity, the coexistence of ferromagnetism and ferroelectricity [14], in contrast with the absence of multi-ferroicity in DTC. The DTN cubic molecular crystal structure has two Ni atoms as magnetic centers in the unit cell. Each Ni has four S atoms and two Cl atoms as nearest neighbors, forming an octahedral structure (like the BO6 octahedra in ABO3 perovskites). The unit cell has 70 atoms and 444 electrons,

Spin-crossover systems are a closely related, highly relevant class, as they are candidate linkers for quantum information systems [15–17]. The electronic structure challenge is to calculate the low-spin, high-spin energy difference and provide the potential surface to calculate the vibrational entropic contributions. A particular significant example is the so-called [Mn(taa)] molecule ([*Mn*3+(*pyrol*)3(*tren*)]), a meridional pseudo-octahedral chelate complex of a single Mn as the magnetic center and the hexadentate tris[(E)-1-(2 azolyl)-2-azabut-1-en-4-yl]amine ligand. It has 53 atoms and 224 electrons. Calculating its spin-crossover energy with low-computational-cost, commonly used density functional methodology without user intervention and tuning has proven to be a formidable task [18].