The band structure problem
To properly describe a defect in an otherwise perfect crystal, using periodic boundary conditions, the defect is usually embedded within a sufficiently large supercell of the host material. The construction of a supercell leads to the back-folding of the bands of the primitive unit cell into the correspondingly reduced 1st Brillouin zone of the supercell (Fig. 1). The introduction of a defect (here a carbon dimer substitution) then leads to a perturbation of the valence and conduction bands and, depending on the defect, to the appearance of one or more localized (i.e., weakly dispersing) states in the band gap. The description of such a system using a tight-binding Hamiltonian obtained by jointly fitting the defect bands and the perturbed valence/conduction bands becomes a formidable task. Directly fitting the TB bands to the respective DFT bands (i.e., minimizing \(| {\epsilon }_{nk}^{DFT}-{\epsilon }_{nk}^{TB}{| }^{2}\) for all n and k) requires their disentanglement: the respective eigenvalues have to be from the same band n at the same k. Tweaking tight-binding parameters to shift the “wrong” eigenstate to match a DFT energy will actually worsen the overall fit. The disentanglement of the bands becomes cumbersome quite quickly with increasing supercell size, and Wannierization reaches its limits already for defective 2D materials35.
Using instead the PDOS offers several advantages. Most importantly, the key problem of the disentanglement of the bands is overcome by projecting the density of states onto the orbitals and atomic sites. Furthermore, the PDOS at the defect site and its close neighbors converges quickly as a function of supercell size. Even in the limit of large and complex defect supercells (where DFT calculations are restricted to a single k-point), the PDOS is still accessible and reliable.
The PDOS has significant advantages for the fitting process, but the perturbation of the host material by the defect also influences the tight-binding Hamiltonian which presents another challenge. The reduction of the symmetry from the primitive cell to the defect supercell leads to many new parameters. In principle, all site energies and the parameters for all possible hoppings between the first, second, and third nearest neighbors within the supercell would need to be refitted. In the case of a defective hBN supercell of 162 atoms, we would be faced with more than 1000 parameters. Although this number can be reduced by symmetry35, in the present work we aim for a minimal parameter set. We therefore employ two distance dependencies for the pristine parameters, resulting in a simpler fitting problem.
Specifically, we used a neural network that predicts these tight-binding parameters for defects from the PDOS. The work flow is shown in Fig. 2. The starting points are optimized geometries of the pristine cell and the defect supercell, calculated using DFT. The pristine tight-binding model is obtained by fitting to the ab-initio band structure in the primitive cell. A distance dependence of the hopping parameters can be obtained by refitting to ab-initio calculations of strained unit cells. Once the pristine tight-binding parameters have been established, we used the defect geometry obtained from a DFT relaxation to build the defect Hamiltonian as a perturbation of the strain-dependent pristine Hamiltonian.
Oval nodes indicate the starting points, cornered nodes represent steps including calculations and rounded nodes represent data. DFT related calculations are displayed in orange, TB in green and ML in purple. Although some similarities are present to Schattauer et al.35, instead of refitting the pristine parameters, we treat the defect as a perturbation to reduce the number of fitting parameters. Furthermore, we introduce an additional step for the generation of the pristine model, where we calculate a distance (strain) dependence of the TB parameters for the long-range description. Last but not least, instead of using the band structure, we use the PDOS to avoid the cumbersome disentanglement of the bands.
Random sampling of the tight-binding parameters that define this Hamiltonian generates the data sets on which we train, and the neural network learns the mapping between the tight-binding PDOS (input) and the tight-binding parameters (output). A key benefit is that the training data is solely acquired within the tight-binding framework, making it computationally inexpensive. After the training, a DFT PDOS is used to predict tight-binding parameters that accurately describe the Hamiltonian defect. Before discussing the details of the tight-binding model, we will introduce the machine learning approach to obtain the parameters.
Machine learning
The key task in our approach is to find tight-binding parameters that give rise to localized defect states within the band gap and simultaneously describe the bands properly. Using standard fitting tools to directly fit the tight-binding PDOS to DFT provides a first hint on the complexity of this task (Supplementary Information). Without prior knowledge of the defect parameters, standard fitting tools struggle to find the localized peak and simultaneously describe the bands properly. Another problem is the scalability of this procedure. While with a small number of parameters, one can still manually fit the parameters, this becomes increasingly difficult the more parameters need to be fitted.
In recent years, state-of-the-art fitting tools have been complemented by machine learning algorithms, enabling the exploration of a complex parameter space with many parameters. It is the goal of this article to make this connection between machine learning and fitting TB parameters via the PDOS: The neural network learns the direct correspondence between the PDOS and the parameters used for the generation, minimizing the chance of being stuck at a local minimum and making sure that the defect state and the bands are properly described.
The PDOS D(r; E) as a fitting property has multiple advantages, but it also inherits some difficulties. The considered defects are deep defects within hBN, a material with a large band gap. For such a material, D(r; E) is mostly zero within the band gap, with narrow defect peaks with large contributions at the defect site, resulting in few activations in the first layer.
To prepare the data for the training, we first simplify the problem by considering D(r; E) only at a small set of points ri and only within a restricted energy range around the band gap to train only on the 1D PDOS data. In addition, we include the sum ∑iD(ri; E) for training, since we get this additional information free. We then normalize each Di(E) = D(ri; E) (between 0 and 1) separately as normalization of the data is essential for training and allows access of all features equally at every site (Fig. 3).
Upper part: scheme of the data preparation. PDOS D(ri; E) at different sites ri is normalized separately. Lower part: Scheme of the neural network, displaying the order of the different parts of the network.
The separately normalized Di(E) are then stacked before entering a 1D convolution (conv1D) layer. The aim is to capture the correlated nature of the PDOS, but also to compress the data and simultaneously reduce the fraction of zeros (Supplementary Fig. 1). Here, the stride defines how much the data is squished while the kernel size determines the broadening of the PDOS features. The output of the conv1D layer is pooled along the channel direction and averages the outputs for different kernels to maintain the features of the input.
The pooling of the conv1D is followed by a scaled dot-product attention as proposed by Vaswani et al.38. The aim is to include long-range dependencies between the defect states and bands. This is achieved by separating the input in query (Q), key (K), and value (V). A matrix-multiplication of the weighted Q and K results in a dependency between each of the entries of the respective inputs. Finally, the matrix is multiplied by the weighted value matrix to incorporate the dependencies into the output. The short formulation of the process reads as38
$$\,\text{Attention}(Q,K,V)=\text{softmax}\,\left(\frac{Q{K}^{T}}{\sqrt{{d}_{k}}}\right)V\,,$$
(2)
here the softmax function associates the matrix entries with a probability and \(1/\sqrt{{d}_{k}}\) stabilizes the training.
The final component of the network architecture is a multi-layer perceptron, which consists of fully connected linear layers. For each layer the ReLu function is used as an activation function. The architecture of the network is depicted in Fig. 3. The final output is evaluated by the standard mean-square error (MSE)
$${\text{loss}}=\frac{1}{{\text{N}}\,}\sum _{n}{({y}_{n}-{\hat{y}}_{n})}^{2}\,,$$
(3)
where \({\hat{y}}_{n}\) is the generated and yn the expected output of the network. We tested including the Gaussian relation for ϵdefect and σ for a given d, but we found no improvement for the final results.
This network is trained using the data produced with the tight-binding approximation. The data set consists of PDOSs, serving as input, and the corresponding labels [\({\epsilon }^{d},\sigma ,{t}_{i}^{d}\)] as output. After training, the network can predict the tight-binding parameters from DFT PDOS obtained from a single calculation.
The network has been developed in parallel with the tight-binding model. Since such architectures solved the parametrization problem at hand, we did not increase its complexity further. For more complex materials and more parameters to be fitted, this can be adjusted. Possible adjustments include further separation of the training at different levels of the neural network. For example, a separation of the conv1D layer with padding in reflection mode enables more data from the bands. On can also use multi-head attention or include more scaled dot-product attention mechanism similar to a transformer neural network.
Defect tight-binding model
The tight-binding Hamiltonian is an approximation of the full many-body Hamiltonian projected onto the localized atomic orbital basis \(\left\vert n\right\rangle\) located at site n. The Hamiltonian can be formulated as follows
$${\hat{{\mathcal{H}}}}_{0}=\sum _{n}{\epsilon }_{n}^{\,\text{prist}\,}\left\vert n\right\rangle \left\langle n\right\vert +\sum _{n,m}\left\vert n\right\rangle {t}_{nm}\left\langle m\right\vert \,,$$
(4)
where diagonal matrix entries \(\left\langle n\right\vert {\hat{{\mathcal{H}}}}_{0}\left\vert n\right\rangle\) are described with the onsite energies \({\epsilon }_{n}^{\,\text{prist}\,}\), whereas the off diagonal elements, \(\left\langle n\right\vert {\hat{{\mathcal{H}}}}_{0}\left\vert m\right\rangle\), are referred to as hopping matrix elements, tnm, from site m to n. Due to the localized nature of the orbitals (fast decay), the hopping parameters can be restricted to the n-th nearest neighbor (usually n ≤ 3), limiting at the same time the number of hopping parameters and the accuracy of the model. The Hamiltonian can be obtained by a fit of the parameters to experiments or ab-initio calculations.
In the case of a monolayer of hBN, we employ a pristine model including the pz orbitals of both N and B with nonzero hopping parameters up to the 3rd nearest neighbors. Furthermore, we use different 2nd nearest neighbor parameters for N-N and B-B hoppings, \({{\rm{t}}}_{{2}^{nd}}^{\,\text{NN}\,}\) and \({{\rm{t}}}_{{2}^{nd}}^{\,\text{BB}\,}\), in order to reproduce the asymmetry in the highest valence and lowest conduction band. Together with the two onsite energies, a total of six tight-binding parameters have been obtained and are shown in table 1. To best describe the valence band maximum and the conduction band minimum, these have been fitted to an ab-initio band structure along \(\overline{KM}\) (Fig. 4a).
a DFT (blue) and the TB (orange) band structure fitted to reproduce the DFT band structure along \(\overline{MK}\). Red dotted line indicates range for later PDOS calculations. b Defect and its nearest neighbors and the adaptations used for the pristine hopping parameters. The short-range (SR) impact on the hopping parameters is described with an additional parameter between the nearest neighbors (1NN) of the defect. The long-range impact effects all other hopping parameters via a distance dependence between the hopping neighbors.
In the more general case of a defective monolayer of hBN, the perturbation introduced by the defect is naturally expected to significantly influence the physics of the system, including changes in the positions of neighboring atoms. The defect has different chemical properties which we account for by new parameters. A new onsite parameter ϵdefect is complemented by three hopping parameters to its nearest neighbors similar to the pristine model, all having the same distance to the defect, as there is no Jahn-Teller distortion for the investigated defects.
Although this introduction of four additional TB parameters for the defect site is straightforward, it is not sufficient because the defect also has a subtle and delicate influence on the chemical properties of its surroundings. Instead of refitting all the hoppings and site energies in the defect supercell, we limit the number of additional fitting parameters to two by shifting the onsite energy of each atom as a function of its distance from the defect site. Additionally, the hopping between pristine neighbors may vary due to geometric deformation. This is taken into account by making the corresponding hopping parameters dependent on the hopping distance.
For the onsite energies, we employ a dependence on the distance from the defect. We use a Gaussian function to describe the perturbation potential introduced by the defect. The function is defined as follows
$$\epsilon ({\epsilon }^{{\rm{defect}}},\sigma ,d)={\epsilon }^{{\rm{prist}}}+\left[{\epsilon }^{{\rm{defect}}}-{\epsilon }^{{\rm{prist}}}\right]\cdot {e}^{-\frac{{d}^{2}}{2{\sigma }^{2}}}\,,$$
(5)
where d is the distance to the defect site and σ is the variance, which is directly proportional to the width of the potential. The dependency is modeled such that the height of the potential is represented by the difference of the onsite energy of the defect (ϵdefect) and the host atom (ϵprist). This approach is similar to the approach of Lambin et al.32 who used a Gaussian potential to investigate the NC substitution defect in graphene in the TB approximation. Accounting for the defect and its perturbation on the onsite parameters (height and width of the Gaussian function), this results in a total of five parameters for the defect itself.
For the defect’s influence on the hopping parameters, we distinguish between short-range and long-range effects. The first model contribution accounts for the local change of the chemical environment introduced by the defect. This is done by means of an additional hopping parameter between its nearest neighbors and is referred to as the short-range (SR) model (Fig. 4b). Because no Jahn-Teller distortion is present for the studied monomers, only one extra parameter needs to be fitted for the defective supercell.
The second contribution to the model accounts for the long-range (LR) component of the perturbation, and it is an expansion of the pristine tight-binding model by a distance dependence between the hopping neighbors. With it we aim at capturing the perturbation on the geometry of the host crystal (Fig. 4b). The four hopping parameters were fitted to strain-dependent band structures, while the onsite energies were kept fixed for later adjustments. For the studied defects, at the DFT level, the maximum change in the distances of the host atoms is 3.7% for the CN and 3.5% for the CB. A biaxial strain within a range of 4% is therefore sufficient to capture the distance dependence for which the band structure was still well described. The strain-dependent fit of the hopping parameters can be found in Supplementary Fig. 4.
The best description of the long-range distance dependence is found to be a quadratic dependence
$$t(d)=\alpha \cdot {(d-{d}_{{prist}})}^{2}+\beta \cdot (d-{d}_{{prist}})+{t}_{{prist}}\,,$$
(6)
which is constructed to reproduce the pristine hopping parameters for the pristine distance between atoms. It adds two more parameters to the respective pristine hopping terms. An advantage of our method is that this distance dependency is naturally well suited for applications where strain plays a role.
To benchmark our methods, we investigated the carbon monomers and the carbon dimer. Both show deep defect states respectively and are perfect candidates to test the workflow. The defect states within the band gap are also pz orbitals and therefore match the pz model for the pristine crystal. In the following, we will benchmark the tight-binding method by fitting to the PDOS and comparing the band structures. We will investigate, in particular the importance of including SR and/or LR contributions to the TB model.
Benchmark results
For both carbon monomers, we want to answer the question of how important adjustments of the pristine hopping parameters are. To do so, we analyze the individual contributions of all components of our model. Therefore, we use the complete model (SR + LR), along with models obtained by inclusion of either short or long range (SR, LR) or the exclusion of both (w/o).
We first check the performance of the neural network with respect to the number of PDOS used as an input for the neural network. To identify the trends, we take into account the density of states projected to the defect site, its first, second and, finally, its third nearest neighbors. Our analysis shows that the usage of the PDOS at three sites is sufficient to describe the monomers (Supplementary Information). We can use these results and generate new data sets with a denser sampling of the PDOS to obtain the best parameters (Table 2).
These parameters are then used for the evaluation of the hopping contributions. We compare the sum of the PDOS values and their cosine similarities up to the 6th nearest neighbor (defect site and 3 atom sites of each element). Since the defect peak outweighs the differences in the bands, we separate the evaluation into different parts of the PDOS, namely the valence band (−0.4–0.1 eV) and the conduction band (4.6–4.85 eV). The cosine similarities of the respective contributions are shown in Fig. 5 and reveal no major differences with high similarities of above 0.96. The valence band is described slightly worse for the models including an additional parameter within the SR contribution. The opposite is observed for the conduction band, where the additional parameter improves the cosine similarity. For all contributions, the cosine similarity of the full PDOS does not change considerably. The CB, is described equally well for all different hopping models.
Cosine similarity of the different hopping description, including the complete model (SR+LR), either just the short-range (SR) or the long-range (LR) model, and without any adaptation. The cosine similarity is calculated for the valence (vb), conduction band (cb), and the full PDOS of CN in (a) and CB in (b). The corresponding PDOS of the respective bands can be found in Fig. 6.
The good agreement of similarities is also reflected in the direct comparison of the summed PDOS of the bands (Fig. 6). Again, we observe a slightly worse description for the valence band when including the SR contribution for the CN, which is reflected in a shift of the valence band maximum. In general, we find that the differences for all contributions are negligible.
Comparison of the summed PDOS at the defect and its 6 first nearest neighbors of the CN in (a) and CB in (b) for DFT and the different hopping approximations.
This indicates that the perturbation of both monomers to the host material is sufficiently accounted for without any hopping contribution. Therefore, considering new defect tight-binding parameters and their impact on the pristine onsite energies is enough to describe the carbon monomers. This is a key advantage for the monomers, as it enables a tight-binding model with few fitting parameters, but also for more complex defects like the carbon dimer.
The carbon dimer consists of two carbon atoms that substitute for one pair of neighboring nitrogen and boron atoms. This results in a shift of the two respective monomer states towards the band edges, an effect that is similar to that of two atoms forming a molecule with a bonding and an anti-bonding state. This requires the introduction of an additional hopping parameter \({\,\text{t}}_{{1}^{st}}^{\text{CC}\,}\). However, our fitting attempts show that this is not sufficient to obtain a good fit for the positions of both defect peaks at the same time. In order to achieve this, one needs to take into account the additional symmetry breaking for the carbon dimer that leads to new distances in the vicinity of the defect. To properly account for all changes, more hopping parameters than those already used for the respective monomers would be necessary. However, we have obtained a distance dependence for the carbon hopping matrix elements that describe the evolution of the respective peaks under compressive and tensile strain (Supplementary Information). These can be used to properly account for the different distances of the carbon atoms from its neighbors. Our calculations for varying distances between carbon monomers indicate that for distances where the Gaussian perturbations of the respective defects overlap, the parameters are not transferable (Supplementary Information). Thus, we have to fit the respective Gaussian dependencies (\({\epsilon }_{\,\text{N}}^{{\rm{defect}}},{\sigma }_{{\rm{N}}},{\epsilon }_{{\rm{B}}}^{{\rm{defect}}},{\sigma }_{{\rm{B}}}\)) and the new hopping parameter (\({\,\text{t}}_{{1}^{st}}^{\text{CC}\,}\)).
In conclusion, the analysis of the long- and short-range contributions showed no major improvements for the monomers. However, to accurately describe strained conditions, the long-range part should be included, since the distance dependence also impacts the peak position. Therefore, further benchmark calculations for the two monomers are performed with the LR contribution. For the fitting of the carbon dimer, we use the model described above, in which we only need to fit five parameters.
The PDOS on the full energy range that we trained provides additional insight into the descriptive power of the predicted tight-binding parameters (Fig. 7). Both carbon monomers show a very good agreement between their respective tight-binding and DFT PDOS. In particular, for CN, the second nearest neighbor has a larger contribution to the defect peak than the first, whereas the opposite is observed for CB. Both are properly described with the respective tight-binding Hamiltonians.
Projected densities of states obtained from DFT (upper subplots) and from the fitted TB parameters (lower subplots) of the (a) CN and (b) of the CB for the defect and its first 6 nearest neighbors including the sum of all. c Projected densities of states of the CBN for the defects and the neighbors depicted in the upper-right sub-panel. The defect peaks are scaled by 0.3 for better visibility in comparison to the band edge contributions to the PDOS. The inset shows the defect peak without rescaling to assess the different contributions of the neighbors.
For the carbon dimer, we use the simplified model described in the previous section, and we use PDOS at both defect sites to obtain the parameters. We observe some differences for the contributions of different atomic sites for the peaks, but the results are satisfactory as we were able to properly describe the main features by only fitting five parameters.
In summary, our results show that the machine learning algorithm is capable of predicting tight-binding models for carbon defects in hBN via the PDOS. Although the tight-binding descriptions are relatively simple, they are able to reproduce the DFT PDOS fairly well. The remaining question is whether the PDOS as a fitting observable is enough to also reproduce the band structure of the defective supercell.
We compare the tight-binding and DFT band structures in Fig. 8, where we use the same tight-binding parameter as for the previous PDOS. As expected, the narrow defect peaks observed in the PDOS result in a correct description of the undispersive defect states in the band gap. The introduction of the defect results in the splitting of the bands from which the defect states emerge. This feature is captured by all tight-binding Hamiltonians as well. However, for CB, we observe an underestimation of this effect which is enhanced by the difference of the conduction band minimum already observed in the pristine fit. This difference is related to the simple pristine model with six parameters. It is not able to capture simultaneously the shape of the conduction band and the position of its minimum (Supplementary Fig. 5). More hopping parameters might fix this, but the aim of this work is to use a pristine model with its limitations to describe a defect. An improvement of the pristine model should also result in a better description of the defective model.
DFT (blue) and tight-binding (orange) band structures of a CN, b CB and c CBN. The tight-binding band structures are calculated with the parameters obtained from the fit to the PDOS (LR). The energy axis has been cut.
The conduction band of the carbon dimer shows a split off band which is captured with the tight-binding model. Although the splitting itself is similar, we observe a difference between the conduction band minima, similar to the CB. Nevertheless, our very simple model reproduces the main features of the band structure, including the valence band, the defect states and the split off conduction band.
Overall, we have demonstrated that the use of the PDOS to obtain tight-binding parameters also shows good agreement for the band structures, thus enabling a novel method to obtain a tight-binding model for defective systems. The observed differences are related to the simplicity of tight-binding models (pristine and defective), rather than to the fitting method.
