11institutetext: Mail.ru Group, National Research University Higher School of
Economics

# Semantic clustering of Russian web search results: possibilities and
problems

Andrey Kutuzov

###### Abstract

The present paper deals with word sense induction from lexical co-occurrence
graphs. We construct such graphs on large Russian corpora and then apply the
data to cluster the results of Mail.ru search according to meanings in the
query. We compare different methods of performing such clustering and
different source corpora. Models of applying distributional semantics to big
linguistic data are described.

##  1 Introduction

The presented paper deals with the problem of semantic clustering of search
engine results page (SERP). The problem arises from the obvious fact that many
user queries are ambiguous in some way. Thus, search engines strive to
diversify their results and to present such results that are related to as
many query interpretations as possible. For example, Google search for the
Russian word ‘максим’ returns:

  1. 1.

five results related to a popular singer,

  2. 2.

two results for a magazine,

  3. 3.

one result for http://lib.ru, Maxim Moshkow’s electronic library,

  4. 4.

one result for a proper name.

However these results are not sorted by their meaning and are returned simply
according to their relevance ranking, which for many of them seems to be
almost equal. The obvious way to cluster the results is by the words their
snippets share. Unfortunately, often snippets for results belonging to one
query sense do not have a single content word in common (except for the query
itself, which is useless). Cf. two snippets for the first query meaning from
the example above:

  1. 1.

‘МакSим начинает самостоятельно заниматься своей карьерой, пишет новые песни.
В этот период певица выступает как малобюджетный проект, …’

  2. 2.

‘МакSим презентовала видеоклип «Я буду жить», получивший широкую огласку еще
до момента появления видео в сети.’

They do not have a single common word, but still belong to one meaning
(popular singer).

Moreover, snippets for different query senses can share some words. Cf. two
snippets from the same search engine results page. They share the word ‘автор’
(‘author’), however the first snippet relates to the first meaning, while the
second snippet shows the third one:

  1. 1.

‘МакSим (Марина Абросимова) – одна из самых популярных и коммерчески успешных
певиц в России, являющаяся автором и исполнителем…’

  2. 2.

‘Работает с 1994 года. Книги и тексты, разбитые по жанрам и авторам.’

That means that there is a need for more sophisticated way to cluster search
results. We should somehow learn which senses the query has and with which
words these meanings are (probabilistically) associated. One of the possible
ways to solve this problem is by extracting co-occurrence statistics from
large corpora. The idea behind this is that word meaning is in fact the sum
(or the average) of its uses. So, meaning is a function of distribution (cf.
[1]). Thus, if we know with which words the query typically co-occurs and how
these neighbors are related to each other, then we know the ‘sense set’ of the
query. After that we can somehow measure semantic similarity of each search
snippet on the SERP with each of the senses and map them to each other. This
information can then be used to either rank the results, or mark them with
appropriate labels.

The structure of the paper is as follows. In Section 2 we briefly overview
work previously done on the subject. Section 3 describes the process of
building co-occurrence graphs from large Russian corpora. In Sections 4 and 5
we conduct an experiment on clustering SERPs with ambiguous queries from
Mail.ru search engine with the help of the methods described before. The
results are evaluated in Section 6. Section 7 draws conclusions concludes and
provides suggestions for further research.

##  2 Related Work

As stated in the previous Section, we are inspired by a fundamental hypothesis
than meaning depends on the distribution [1] and that frequency of linguistic
phenomena (in our case, word co-occurrence) is important for determining these
phenomena’s place in the system of language [2]. Our work is also based on the
idea that the senses of ambiguous lexical units should be induced from the
data itself, not from a dictionary. No dictionary is perfect or comprehensive,
because ‘senses as identified in the dictionary identify points on a continuum
of possibilities for how the word is used’ [5]. The only robust source of
words’ meanings in the text is the text itself. That’s why we shift our focus
away from selecting the most suitable senses from a pre-defined inventory
towards discovering senses automatically from the raw data, which is natural
text.

One of the first notes on practical application of this idea to word sense
disambiguation and word sense induction is found in [3], where vector
representations of word similarity derived from co-occurrence data are used.
Broad review of contemporary (by 2012) state of the field is provided in [4].

The main source of methods for our present research is [6], which describes
workflow for clustering web search results using graph analysis over co-
occurrence networks. Specifically, we use the notion of query graph,
consisting of query terms and words from search engine results page augmented
with nearest neighbors and relations from a reference corpus. For partitioning
query graph and clustering query senses we employed Curvature algorithm [6]
and Hyperlex algorithm proposed in [7].

##  3 Building Co-Occurrence Graph

The first thing we had to do was to select a text corpus to build the graph
upon. It is well known that the larger the corpus is the more co-occurrence
information it contains. However, increasing corpus size also leads to
exponentially growing computation time. Thus, for the sake of time and because
of the preliminary nature of our research, we restricted ourselves to three
Russian corpora of smaller but still decent size:

  1. 1.

Open Corpora111http://opencorpora.org (1 million tokens), further ‘OC’;

  2. 2.

Disambiguated fragment of Russian National Corpus222http://ruscorpora.ru (1
million tokens), further ‘RNC’;

  3. 3.

Corpus of random search queries from Mail.ru search engine333http://go.mail.ru
(2 million tokens), further ‘QC’.

The first two items are academic corpora of Russian texts, supposedly
representing (written) language in general. They differ in that the first one
consists of full texts published under various free and open licenses, while
the second one is a random sample of sentences from the larger Russian
National Corpus. Both of them come with morphological annotation.

The third corpus was taken for comparison. It is important in view of the aim
of our research (to test semantic SERP clustering). Our intuition was that
perhaps query corpus provides more ‘real-life’ sense inventory. It is two
times as big as its counterparts, because ‘connectivity’ between its members
is lower (see Table 2) and we had to compensate for this.

At the same time, it turned out that the first two corpora mixed into one give
better results, thus below we will often refer to such ‘meta-corpus’ as ‘Mix
corpus’.

Before constructing the graph itself, we preprocessed the corpora, namely:

  1. 1.

Removed from QC all queries which did not contain Cyrillic characters (as
apparently they are not Russian),

  2. 2.

Processed QC with Freeling analyzer [8] to extract lemmas and morphological
information for all tokens,

  3. 3.

Removed stop words,

  4. 4.

Removed all tokens except nouns, as we restrict ourselves to inducing only
nominal senses (the same strategy was applied in [6]).

Sizes of preprocessed corpora are given in Table 1.

Table 1: Sizes of corpora participating in the experiment Corpus | Size (tokens)
---|---
OC | 490671
RNC | 294849
Mix | 785520
QC | 1035483

Average query length in QC is 2.47 noun tokens per query.

After the corpus has been built, the process of constructing co-occurrence
graph is rather straightforward: we create an empty graph and then populate it
with vertexes denoting word types in the text (lemmas). After that for each
lemma we find all its immediate neighbors in the corpus, that is, words to the
left and to the right (sentence boundaries not crossed, queries considered to
be ‘sentences’ as well). If two lemmas were neighbors at least one time, we
draw an edge between corresponding neighbors.

Finally, we have an undirected graph in which noun lemmas are vertexes and co-
occurrence relations are edges. For each edge we also calculate Dice
coefficient [9]. It measures the ‘strength’ of the collocation, based on
absolute frequency (c) of both words (w and w’) and collocation (w,w’):

$Dice(w,w^{\prime})=\frac{2c(w,w^{\prime})}{c(w)+c(w^{\prime})}$ (1)

One can also think about the graph as a matrix of Dice coefficient values for
all possible pairs of lemmas in the corpus.

Table 2 gives an overview of the basic features of the graphs.

Table 2: Parameters of the graphs Corpus |  Vertexes |  Edges |  Average degree |  Average path length |  Clustering coefficient
---|---|---|---|---|---
OC |  21881 |  257846 |  23.57 |  3.26 |  0.166
RNC |  22467 |  163914 |  14.6 |  3.53 |  0.136
Mix |  31984 |  395225 |  24.7 |  3.29 |  0.186
QC |  85548 |  291033 |  6.8 |  4.07 |  0.16

One can see that the average degree of QC is lower in comparison with the
other corpora (because queries are typically shorter that sentences in natural
texts). That is one of the reasons for our decision to use a larger query
corpus.

It should also be noted that all corpora comply to ‘small world’ definition
[10], because their average path length is approximately the same as in a
random graph with the same number of vertexes ($N_{V}$) and average degree
($A_{D}$), while clustering coefficient is significantly higher than it should
be in the random graph. For example if Mix corpus were a random one, its
average path length would be equal to 3.24 ($=\frac{log(N_{V})}{log(A_{D})}$),
very close to the actual value. However, in this case, its clustering
coefficient should be 0.0015 ($=\frac{2\times A_{D}}{N_{V}}$), which is
significantly lower than the actual value. The same is true for all other
corpora.

‘Small world’ nature of our graphs means that vertexes in them tend to bundle
into clusters, which is typical of many real-world networks. This finding
supports the idea of extracting senses from such clusters. It also
additionally proves the applicability of graph sense induction methods to our
corpora, as English-language graphs in the related publications also showed
such properties.

##  4 Building Query Graph

We experimented with clustering search engine results page on a set of sixty
ambiguous one-word Russian queries, taken from Analyzethis homonymous queries
analyzer444http://analyzethis.ru/?analyzer=homonymous. Analyzethis is a search
engines evaluation initiative, offering various search performance analyzers,
including one for ambiguous or homonymous queries. We crawled Mail.ru search
for these queries, getting titles and snippets (10 for each result).

The procedure of semantic clustering starts with building the so called query
graph. Here we closely follow [6].

First, we lemmatize all snippets and titles and remove stop words and the
query word itself. Then we construct a graph $G_{q}$ with all nouns from
snippets and titles as vertexes. Then we use one of the large corpora graphs
(those that we built in Section 3) to find words strongly connected to the
query word and add these words to the query graph. We consider a connection
‘strong’ if it falls under the following constraints:

$\left\\{\begin{aligned} \frac{c(q,w)}{c(q)}\geq 0.01\\\ Dice(q,w)\geq
0.005\\\ \end{aligned}\right.$ (2)

where c is absolute frequency in the corpus, q is the query and w is the word
under analysis. Thresholds 0.01 and 0.005 were determined empirically while
experimenting on the above mentioned ambiguous queries set. These thresholds
produced most convincing sense clustering. However, the issue of choosing the
thresholds is a subject for thorough evaluation in future.

Thus, we now have $G_{q}$ with no edges and vertex set consisting of words
from the search result and strong neighbors of the query word. After that, for
each pair of words (w,w’) in $G_{q}$ we check if they co-occur in the large
corpus. If they do and $Dice(w,w^{\prime})\geq 0.005$, we connect these words
in $G_{q}$ with an edge with weight = Dice(w,w’). Finally, we delete
disconnected vertexes (those with the degree equal to 0).

##  5 Processing Query Senses and Results

With query graph at hand, we are ready to find which senses the query has.
What we need is an optimal partition of the query graph, in which words
related to different senses are in different parts of the graph. We apply two
techniques for that, namely, Curvature from [6] and Hyperlex from [7].

Figure 1: Query graph for ‘амур’ (Curvature)

###  5.1 Curvature

Curvature algorithm aims at finding vertexes from $G_{q}$ with low local
clustering coefficient. Our hypothesis is that these are words which serve as
‘links’ between different senses or ‘uses’ of the query. Then we remove
vertexes with clustering coefficient below a certain threshold. It leads to
the graph disjointing into several components related to different senses.
Vertexes in these components represent lexical inventory of each sense.
Disconnected vertexes are removed from the final graph.

Let us illustrate the process with the example of ‘амур’ (‘Amur’) query.
Figure 1 shows its query graph. It is already disconnected into two components
and the meaning of love god (associated with words ‘лук’, ‘стрела’ and
‘юноша’) is separated. However, other ‘senses’ of the query remain hidden in
the giant component. Vertexes shown as triangles have low clustering
coefficient and are thus marked for deletion.

So, we delete ‘triangular’ vertexes. Note that we chose threshold 0.3 – all
vertexes with clustering coefficient below this are removed. It is also
important that we do not delete vertexes with clustering coefficient = 0. This
is because neighbors of such vertexes are not connected to anything except
this vertex. If we remove it, a lot of disconnected vertexes will appear. Such
clusters (consisting of only one word) do not make much sense. For example,
the word ‘лук’ on Figure 1 is characterized by clustering coefficient = 0. If
we remove it, then the whole component representing ‘love god’ meaning
disappears.

Figure 2 shows the query graph after removing vertexes with low clustering
coefficient. We now have 6 components (note that the labels for these clusters
are introduced by us, not by the algorithm):

  1. 1.

River (all vertexes except enumerated below)

  2. 2.

Love god (‘юноша, лук, стрела’)

  3. 3.

Hockey club (‘клуб, болельщик’)

  4. 4.

Movie (‘любовь, фильм’)

  5. 5.

Dictionary-1 (‘календарь, словарь, википедия, энциклопедия, академик, сюжет’)

  6. 6.

Dictionary-2 (‘значение, описание, характеристика’)

Figure 2: Disjointed query graph (Curvature)

First 4 components clearly represent different meanings of the word ‘амур’.
The last two are rather ‘uses’, typical contexts. However they can still be
useful in clustering as they allow to keep encyclopedic results together.

###  5.2 Hyperlex

Hyperlex algorithm described in [7] introduces the notion of ‘hubs’ within the
graph, meaning most inter-connected vertexes and employs the graph’s maximum
spanning tree. Just like the previous algorithm, it takes as an input the
query graph $G_{q}$ we prepared in Section 4 and the query itself.

First we create a list L with all vertexes from $G_{q}$ sorted in decreasing
order by their absolute frequency in the large corpus. Then for each item of
this list we check if the corresponding vertex complies to the following
constraints:

  1. 1.

Vertex normalized degree is greater than or equal to 0.05,

  2. 2.

Average Dice coefficient of vertex edges is greater than or equal to 0.007.

If the constraints are met, we add this word to the hub list, considering it
to be a kind of a connector. Simultaneously, we remove this vertex and its
neighbors from the list L and continue iterating. In case we meet a word which
does not satisfy the requirements above, we check whether the list of hubs has
at least two elements. If it does, we stop iterating, if not, we continue to
the next item. Note that it differs from the original Hyperlex algorithm,
where one should stop no matter how long the hub list is. In our Russian
material it sometimes caused the hub list to remain empty or contain only one
item, which is useless.

After we have the list of hubs, we augment $G_{q}$ with query vertex and
connect this vertex to all hubs putting infinite (or very high) Dice
coefficient on the corresponding edges. Then, we produce a maximum spanning
tree from this graph. Maximum spanning tree is an attempt to keep all the
vertexes connected while eliminating cycles and using as few edges as possible
with as high weights (in our case it is Dice coefficient) on them as possible.
In the spanning tree, there is only one path between any two vertexes and this
path lies through edges with maximum Dice. Because the query vertex and the
hubs are connected by edges with infinite Dice, they are sure to be the center
of the spanning tree and directly linked.

Figure 3: Query graph for ‘амур’ after augmenting it with the query vertex
(Hyperlex)

At last we remove the query vertex from the spanning tree, producing
disjointed subtrees with hubs as roots. These subtrees represent query
meanings. Note that we also delete all disconnected vertexes (those with
degree = 0).

Let us present an example of Hyperlex at work with the same query ‘амур’. Our
corpus is Mix. Initial state of the query graph $G_{q}$ is the same as in
Figure 1.

We add the word ‘амур’ to the graph $G_{q}$ and connect it to vertexes
selected as hubs: ‘область, фильм, команда’. The result is presented on Figure
3 with query vertex drawn as a diamond. For reference, vertexes which were
introduced from the corpus and not from search results (‘сюжет’, ‘океан’, etc)
are drawn as triangles.

Now we produce maximum spanning tree from $G_{q}$ with Dice coefficient as
weight measure. The tree is visualized on Figure 4. Note that it has much
fewer edges than the initial $G_{q}$.

Finally, we remove the query vertex and all vertexes that become disconnected
after this removal. As a result, we have a disjointed graph shown in Figure 5.
The number of the components has grown from 2 to 4 (once again, labels are
assigned by us):

  1. 1.

Love god (‘юноша, лук, стрела’)

  2. 2.

Movie (‘любовь, фильм’)

  3. 3.

Hockey club (‘клуб, игра, болельщик, команда, цвет’)

  4. 4.

River (all the remaining vertexes)

Figure 4: Maximum spanning tree for ‘амур’ query graph (Hyperlex)

One can see that Hyperlex successfully extracted the same four important
meanings as the previous algorithm. At the same time, unlike Curvature, it
managed to avoid two ‘encyclopedic’ clusters (obviously in common for too many
queries) and leave their vertexes in the ‘river’ cluster. Also, Hyperlex is
better because it describes ‘hockey club’ cluster in a richer way, using 5
relevant words instead of 2.

One can again note that in fact what we call ‘senses’ are not senses like
meanings in the dictionaries. We agree with Jean Veronis who argues that co-
occurrence networks reflect ‘uses’ rather than senses. So, what we have are
typical environments where the word is used, and these environments are only
loosely connected to what a lexicographer would call ‘senses’ or ‘meanings’.
However, we are fine with that, as we assume that clustering SERP according to
‘typical uses’ is at least equally important as clustering according to
‘proper senses’. Perhaps, these senses are in fact less related to real-life,
as even linguists sometimes have trouble matching the ‘senses’ found in a
dictionary and the occurrences found in a corpus [7]. Additionally, as has
already been stated, dictionary senses are always limited and by design cannot
cover new semantic trends and subtle meanings quickly appearing and
disappearing in the modern world. Thus, theoretically typical uses are more
relevant for clustering than academic dictionary senses. To strictly prove it
for the Russian material, one needs manually clustered data set (see Section
6), and we leave it for further research.

###  5.3 Mapping Results to Senses

Once we possess the sense inventory for the query, we can combine it with
bags-of-words for each search result to finally perform SERP clustering. We do
that in a rather straightforward way.

Given a set of senses represented by a lemma set each and a set of results
(snippet and title) represented by lemma sets as well, for each pair of result
($r$) and sense ($s$) we calculate similarity measure $sim$. It is a simple
number of lemmas in common for both sets divided by the number of lemmas in
the result:

$sim(r,s)=\frac{r\cap s}{length(r)}$ (3)

Then we choose the sense with maximum similarity and link this sense to the
result. Thus, each result receives some sense, and is ‘understood’.

In the future we plan to explore other means of calculating similarity measure
as well, for example, counting tokens not types or considering weights on
edges in the intersection.

Figure 5: Maximum spanning tree after removing query vertex (Hyperlex)

##  6 Evaluation of SERP clustering

Generally, evaluation of clustering is a rather harsh task. Perhaps, the best
way to do this is to employ human assessment, but for the time being we
limited ourselves to simple evaluation of the correctness of cluster number
(that is, number of meanings).

Analyzethis service provides data about how many senses of an ambiguous query
are there in the SERP. Thus we consider it to be an expert opinion and check
how strong is our deviation from this ‘gold standard’. For example, if
Analyzethis believes that there are three senses present on the SERP, and our
clustering algorithm puts all the results into one cluster, this signals that
the algorithm is not optimal. The same is true if the number of clusters is,
for example, eight. The less our deviation from Analyzethis assessment is, the
better. So, in fact we check that the employed algorithms do not produce
senseless results (too many or too few meanings). We once again note that in
order to evaluate the contents of the clusters themselves, one needs manually
clustered SERPs for ambiguous queries. To our knowledge, there is no such a
data set for Russian. We are working on creating it.

For the time being, we compared the number of clusters for each of ambiguous
queries in four different settings (two corpora and two word sense induction
methods). Then we calculated average deviation of our clustering number from
that of Analyzethis. Table 3 provides the results of this comparison. Note
that the average number of senses per query in Analyzethis data set was 2.65.

Table 3: Evaluation of SERP clustering (average deviation from Analyzethis assessment in number of senses and in percent from the average number of senses in the set) Corpus | Curvature | Hyperlex
---|---|---
Mix | 1.636 (61%) |  1.288 (49%)
Query | 1.742 (66%) | 1.379 (52%)

It is clear that Hyperlex consistently outperforms Curvature, and that Mix
corpus does the same with the query corpus. Hyperlex victory comes as no
surprise, as it uses maximum spanning tree notion, which seems to allow deeper
grasping of graph structure. The victory of Mix corpus (which is smaller than
the query corpus) is much less expected. We believe that there are two reasons
for this:

  1. 1.

As we have already mentioned, the query corpus is less ‘dense’ because of low
length of queries. Thus, there are fewer edges and less data for algorithms.

  2. 2.

Query corpus was lemmatized with Freeling while Mix corpus consists of
manually annotated corpora. Glitches and outright errors of Freeling could
impact graph quality. This can be fixed in the future either by improving
Freeling or by using another lemmatizer.

Thus, at the moment, using Mix corpus and Hyperlex algorithm of word sense
induction seems to be the best option. However, things surely can be different
if we employ larger corpora (which we plan to do in the future).

##  7 Conclusion and future work

We showed that state-of-the-art methods of word sense induction and search
results clustering based on semantic graphs do work for Russian data.

Application of such methods can lead to search engine results presentation
getting closer to actual semantics of the results, not simply term frequency
ranking. For a user, it would mean the possibility to immediately grasp which
results in the SERP are actually related to the query sense, and which other
senses exist. The power of this approach can be increased by wider employment
of Semantic Web paradigm: semantically marked up web pages are represented by
generally better and clearer snippets. Such snippets, in turn, should provide
better data for graph-based word sense induction algorithms.

We plan to experiment on more types of query graph processing and launch a
full-scale human evaluation of results. Also, it seems profitable to use not
only separate words, but also compound phrases, as well as to construct graphs
with not only immediate neighbors, but also with second-order co-occurrences
(neighbors of neighbors). Additionally, experiments with larger query corpora
may lead to new and inspiring insights in this field.

## References

  * [1] Harris, Z.S.:  Distributional structure.  Springer (1970)
  * [2] Bybee, J.:  Frequency of use and the organization of language.  Oxford University Press, USA (2006)
  * [3] Schütze, H., Pedersen, J.O.:  Information retrieval based on word senses.  In: Proceedings 4th Annual Symposium on Document Analysis and Information Retrieval (SDAIR 1995). (1995) 161–175
  * [4] Navigli, R.:  A quick tour of word sense disambiguation, induction and related approaches.  In: SOFSEM 2012: Theory and practice of computer science.  Springer (2012) 115–129
  * [5] Kilgarriff, A.:  Dictionary word sense distinctions: An enquiry into their nature.  Computers and the Humanities 26(5-6) (1992) 365–387
  * [6] Marco, A.D., Navigli, R.:  Clustering and diversifying web search results with graph-based word sense induction.  Computational Linguistics 39(3) (2013) 709–754
  * [7] Véronis, J.:  Hyperlex: lexical cartography for information retrieval.  Computer Speech & Language 18(3) (2004) 223–252
  * [8] Padró, L., Stanilovsky, E.:  Freeling 3.0: Towards wider multilinguality.  In Calzolari, N., Choukri, K., Declerck, T., Doğan, M.U., Maegaard, B., Mariani, J., Odijk, J., Piperidis, S., eds.: Proceedings of the Eight International Conference on Language Resources and Evaluation (LREC’12), Istanbul, Turkey, European Language Resources Association (ELRA) (may 2012)
  * [9] Smadja, F., McKeown, K.R., Hatzivassiloglou, V.:  Translating collocations for bilingual lexicons: A statistical approach.  Computational Linguistics 22(1) (1996) 1–38
  * [10] Watts, D.J., Strogatz, S.H.:  Collective dynamics of ‘small-world’networks.  Nature 393(6684) (1998) 440–442


#  Spinney: post-processing of first-principles calculations of point defects
in semiconductors with Python

Marco Arrigoni  Georg K. H. Madsen  Institute of Materials Chemistry, TU Wien,
A-1060 Vienna, Austria

###### Abstract

Understanding and predicting the thermodynamic properties of point defects in
semiconductors and insulators would greatly aid in the design of novel
materials and allow tuning the properties of existing ones. As a matter of
fact, first-principles calculations based on density functional theory (DFT)
and the supercell approach have become a standard tool for the study of point
defects in solids. However, in the dilute limit, of most interest for the
design of novel semiconductor materials, the “raw“ DFT calculations require an
extensive post-processing. Spinney is an open-source Python package developed
with the aim of processing first-principles calculations to obtain several
quantities of interest, such as the chemical potential limits that assure the
thermodynamic stability of the defect-laden system, defect charge transition
levels, defect formation energies, including electrostatic corrections for
finite-size effects, and defect and carrier concentrations. In this paper we
demonstrate the capabilities of the Spinney code using c-BN, GaN:Mg, TiO2 and
ZnO as examples.

###### keywords:

first-principles; point defects; charged defects; thermodynamic stability;
density functional theory; Python

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program Title: Spinney
Licensing provisions: MIT
Programming language: Python 3
Extetrnal libraries: NumPy [1], SciPy [2], Pandas [3], Matplotlib [4], ASE [5]
Nature of the problem:
Post-processing of first-principles calculations in order to obtain important
properties of defect laden systems in the dilute-limit: chemical potential
values ensuring thermodynamic stability, thermodynamic charge transition
levels, defect formation energies and corrections thereof using state-of-the-
art corrections schemes for electrostatic finite-size effects, equilibrium
defect and carriers concentrations.
Solution method:
Flexible low-level interface for allowing the post-processing of the raw fist-
principles data provided by any computer code. High-level interface for
parsing and post-processing the first-principles data produced by the popular
computer codes VASP and WIEN2k.
Additional comments:
An extensive documentation is available at: https://spinney.readthedocs.io
The code is hosted on GitLab: https://gitlab.com/Marrigoni/spinney

##  1 Introduction

Point defects in semiconductor and insulator materials (from now on we will
use the general term “semiconductor” to denote any material with a non-zero
fundamental gap) play a fundamental role in a wide range of technological
applications such as electronics, optoelectronics, electrochemistry and
catalysis to name a few [6, 7, 8, 9, 10, 11]. The introduction of controlled
amounts of extrinsic atomic species, usually in very low concentrations, can
have greatly beneficial effect on the properties of the host material and has
fueled the enormous growth of the semiconductor industry. At the same time
however, the unintentional introduction of impurities can have a detrimental
effect on the performance and lifetime of a device. In addition to extrinsic
species, intrinsic point defects, whose existence is guaranteed at
thermodynamic equilibrium, play a fundamental role in defining the properties
of any material and their response to doping [12, 13, 14].

The desire to control materials properties through point defects has lead to
the need of a detailed understanding of the atomic-scale structure and
electronic properties of such defects. While it is sometime possible to
experimentally obtain details of the atomic structure of surface defects [15,
16], acquiring such information remains generally challenging. First-
principles density functional theory (DFT)[17, 18] calculations give the
possibility of both atomic-scale structural characterization as well as the
calculation of thermodynamic quantities[19, 20, 21], and have become a
standard tool for complementing experimental investigations.

The vast majority of such calculations employ the supercell approach, in which
a point defect is embedded in a simulation cell containing multiple
repetitions of the primitive cell of the host material (_i.e._ the supercell).
However, for many practical applications the dilute limit, in which defect-
defect interactions can be considered negligible, is of interest. The modeling
of such a limit poses a challenge for the supercell approach, as commonly used
supercell sizes entail much larger defect concentrations. Moreover, point
defects in semiconductors can be ionized and therefore realize different
charge states [22]. While most of the spurious defect-defect interactions
present in the supercell method can be minimized by employing a moderate
supercell size, the long-ranged nature of the Coulomb potential means that
correcting spurious electrostatic interactions is fundamental in order to
obtain reliable defect formation energies [23].

The computer code Spinney was developed with the aim of obtaining the most
important properties of the defect-laden systems from the results of DFT
supercell calculations by implementing non-trivial post-processing operations.
While other codes aimed at similar properties support only calculations
performed by specific DFT packages,[24, 25, 26] Spinney offers a flexible low-
level interface operating on built-in Python or NumPy data structures. This
allows users of general DFT softwares to extract the required information,
feed it to the relevant Spinney subroutines and obtain the defect properties
of interest without being bound to use code-specific formats and data
structures. At the same time, Spinney provides user-friendly interfaces for
selected popular DFT codes, able to automatically parse the DFT calculations
outputs, extract the needed information and produce properties of interest. At
the current stage, such interfaces exist for the VASP [27, 28] and WIEN2k [29,
30] first-principles codes.

In section 2 we review the theoretical background regarding the point defects
in the dilute limit and describe the formalism which is employed by Spinney.
Section 3 gives a demonstration on how Spinney can be used to calculate common
properties of defect-laden system using important semiconductor and oxide
materials as explicit examples.

##  2 Background

###  2.1  Defect formation energy and thermodynamic charge transition levels
in the dilute limit

A key quantity that characterizes a point defect is its formation energy. The
most appropriate thermodynamic potential for the formation of point defects is
the grand potential. Therefore it is natural to define the formation energy,
$\Delta E_{d}(D^{(q)})$, of a given point defect $D$ in the charge state $q$
as the change in grand potential after the introduction of the defect in the
pristine host material [31]:

$\Delta E_{d}(D^{(q)})=\Delta
E_{f}(D^{(q)})-\sum_{i}n_{i}\Delta\mu_{i}+q\mu_{e}+E_{corr}$ (1)

where $\Delta E_{f}(D^{(q)})=E(D^{(q)})-E_{\mathrm{bulk}}-\sum_{i}n_{i}E_{i}$
is the defect formation energy with respect to the reference states of the
parent elements. $E(D^{(q)})$ is the energy of the supercell containing the
point defect, $E_{\mathrm{bulk}}$ is the energy of the supercell describing
the pristine material, $n_{i}$ is the number of atoms of type $i$ which need
to be removed ($n_{i}<0$) or added ($n_{i}>0$) to the system in order to
create the point defect, $\Delta\mu_{i}$ is the chemical potential of the
$i$-th element with respect to the standard state of the element and $\mu_{e}$
is the chemical potential of the electron.

As a common approximation, the Gibbs free energy of the solid has been
replaced by the ground-state DFT energy, $E$. For defects inducing large
atomic distortions on the host material, harmonic [32, 33, 34], and even
anharmonic contributions [35], can have a non-negligible effect at high
temperatures. However, the calculation of such terms is generally time-
consuming and is customarily neglected. On the other hand, large errors arise
if thermal contributions are neglected in the chemical potential of gas
species. Therefore, one needs to explicitly account for the pressure and
temperature dependence of gas species if these appear in equation (1).
Furthermore, $E(D^{(q)})$ must contain a correction to the "raw" DFT energy
due to the finite-size-errors which arise in the supercell method [21], which
is denoted as $E_{corr}$ in equation (1). More detailed discussions of the
chemical potentials and finite size corrections are given in Sections 2.2 and
2.3.

A thermodynamic charge transition level is defined as the value that the
electron chemical potential must have in order for two different charge states
of a defect, $q$ and $q^{\prime}$, to have the same defect energy. It is
customarily to express the electron chemical potential in terms of the valence
band maximum of the host material, $\epsilon_{V}$, and the Fermi level,
$E_{F}$, which varies between zero and the band gap of the material:
$\mu_{e}=\epsilon_{V}+E_{F}$. With this convention and equation (1), we can
then write the charge transition levels as

$\epsilon_{0}(q/q^{\prime})=\frac{E(D^{(q)})-E(D^{(q^{\prime})})}{q^{\prime}-q}-\epsilon_{V}+\Delta
E_{corr}(q/q^{\prime}),$ (2)

The expression emphasizes the fact that charge transition levels do not depend
on the chemical potentials of the elements but they do depend on the valence
band maximum and the corrections for finite size effects.

###  2.2  Determination of equilibrium chemical potential values

Chemical potentials quantify the energy cost necessary for exchanging atomic
species between the system and a particle reservoir. The values of $\mu_{i}$
and $\mu_{e}$ are thus important, not only because they determine the defect
formation energy, equation (1), but also because they connect the first-
principles calculations with the experimental growth conditions.

Thermodynamic equilibrium constrains the possible values the chemical
potentials can assume. Considering, without loss of generality, a binary
compound, $\mathrm{M_{x}N_{y}}$, the values of $\Delta\mu_{M}$ and
$\Delta\mu_{N}$ are restricted by the following conditions:

$x\Delta\mu_{M}+y\Delta\mu_{N}=\Delta h_{f}(\mathrm{M_{x}N_{y}}),$ (3a)
$\nu\Delta\mu_{M}+\mu\Delta\mu_{N}\leq\Delta
h_{f}(\mathrm{M}_{\nu}\mathrm{N}_{\mu}),$ (3b) $\Delta\mu_{M}\leq
0\quad,\quad\Delta\mu_{N}\leq 0.$ (3c)

$\Delta\mu_{\zeta}$ represents the change in chemical potential from the
standard state of species $\zeta$:
$\Delta\mu_{\zeta}=\mu_{\zeta}-\mu^{\circ}_{\zeta}$, and $\Delta h_{f}$ is the
enthalpy of formation (per formula unit) of the compound. Equation (3a)
represents the thermodynamic stability condition for $\mathrm{M_{x}N_{y}}$ and
shows that only one of the elemental chemical potentials is an independent
variable. Equation (3b) constrains the value of the chemical potentials so
that $\mathrm{M_{x}N_{y}}$ is stable and do not segregate into other phases.
Likewise, equation (3c) demands that segregation of the parent compounds are
avoided.

As mentioned above, thermal contributions to the free energies of gas phases
must be taken into account. Suppose that in the standard state species N is in
the gas phase Nγ, we can model the dependence of $\Delta\mu_{N_{\gamma}}$ from
$T$ using a set of parameters fitted to experimental data (Shomate equation)
[36] and model the dependence on $p$ by using an ideal gas model:

$\Delta\mu_{N}(T,p)=\frac{1}{\gamma}\left(h_{N_{\gamma}}(0,p^{\circ})+\Delta
g_{N_{\gamma}}(T)+k_{B}T\ln\left(\frac{p}{p^{\circ}}\right)\right)$ (4)

where $p^{\circ}=$ 1 bar is the standard pressure and $\Delta
g_{N_{\gamma}}(T)=g_{N_{\gamma}}(T,p^{\circ})-g_{N_{\gamma}}(0,p^{\circ})$.
$h_{N_{\gamma}}(0,p^{\circ})$ represents the enthalpy per molecule at standard
pressure and 0 K and it is usually approximated by the DFT-calculated
electronic energy. This model is implemented in Spinney for the most common
binary gas molecules: O2, N2, H2, Cl2 and F2. The model can also be used for
other gases taking the needed parameters for calculating $\Delta
g_{N_{\gamma}}$ from the NIST-JANAF themochemical tables [36].

###  2.3  Corrections for electrostatic finite-size effects

Due to the periodic boundary conditions, the introduction of a point defect in
a supercell entails the presence of periodic images of the defect itself. Such
an artificial array of defects furthermore generally corresponds to high
defect concentrations, far from the dilute limit. While most of the resulting
defect-defect interactions can be minimized by using a large enough supercell,
long range electrostatic interactions due to charged defects cannot be
neglected for any realistic supercell size. In such a case, spurious defect-
defect interactions can considerably affect the predicted energy of the point
defect and must be corrected for.

As a point defect is introduced in the host material, there will be a
redistribution of the charge density of the latter. In the ideal case of an
isolated point defect, this defect-induced charge density can be described by
$\rho_{isol}$. On the other hand, when a supercell is employed, the
application of periodic boundary conditions will yield a different charge
density: $\rho_{per}$. Moreover, for charged point defects, periodic boundary
conditions require the introduction of a neutralizing background, usually
taken as an homogeneous jellium of density $-\frac{q}{V}$, where $q$ is the
defect charge state and $V$ is the supercell volume, in order to ensure
convergence of the electrostatic energy [22]. The defect-induced electrostatic
potential can generally be obtained from the induced charge density by solving
Poisson equation, $\nabla^{2}\phi(\mathbf{r})=-4\pi\rho(\mathbf{r})$, with the
proper boundary conditions. This will yield the potentials $\phi_{isol}$ and
$\phi_{per}$ for $\rho_{isol}$ and $\rho_{per}$, respectively.

Assuming that the defect-induced charge density is completely contained in the
supercell and is not affected by the presence of periodic boundary conditions,
the corrective term can be calculated as:

$E_{corr}=\frac{1}{2}\int_{V}\left(\phi_{isol}(\mathbf{r})-\phi_{per}(\mathbf{r})\right)\rho_{d}(\mathbf{r})\,d\mathbf{r}$
(5)

Multiple correction schemes for electrostatic finite-size effects have been
developed [23, 37, 38, 39, 40, 41, 42, 43] in order to estimate $E_{corr}$.
Among the proposed approaches, the methods applying corrections _a posteriori_
using simple models for the defect-induced charge density have become perhaps
the most popular due to their flexibility, speed and reliability.

Spinney implements the _a posteriori_ schemes proposed by Freysoldt,
Neugebauer and Van der Walle (FNV) [41] and the development thereof proposed
by Kumagai and Oba (KO) [42] where the correction energy is be expressed as:

$E_{corr}=-E_{lat}+q\Delta\phi$ (6)

where $E_{lat}$ is the Madelung energy of $\rho_{d}$ embedded in the host
material and the jellium background when periodic boundary conditions are
present and $\Delta\phi$ is an alignment term for the electrostatic potential.
The scheme assumes that the defect-induced charge density is spherical and
introduces a model charge density from which the terms are calculated. In the
FNV scheme this is generally taken as a linear combination of a Gaussian and
an exponential functions. Using such a model for $\rho_{d}$ gives quite some
flexibility for model the defect-induced charge density and allows to
calculate analytically the Madelung energy for isotropic systems but analytic
calculations for anisotropic systems are not possible. On the other hand, KO
showed the loss of flexibility in using a point-charge (PC) model is generally
small for most system and a PC model allows the analytic calculation of the
Madelung energy in anisotropic systems [44, 45]. As a matter of fact, the KO
approach has been successfully applied to a wide range of materials, showing
that, once the defect-induced charge density is well localized within the
supercell, a convergent defect formation energy can be obtained using
supercells of moderate size [42].

The potential-alignment term, $\Delta\phi$, is obtained by comparing the
electrostatic potential far from the defect position of the pristine and the
defect containing supercells with the one of the model charge density. It can
be decomposed into a sum of different terms [41, 46, 42]:

$\Delta\phi=V_{d,q}-V_{q/b}\equiv-\Delta V_{d,q/b},$ (7)

where $V_{d,q}$ is the electrostatic potential produced by the model charge
distribution and $V_{q/b}=V_{\mathrm{defect},q}-V_{\mathrm{bulk}}$ is the
defect-induced potential as calculated in the first-principles simulation. FNV
perform the comparison using planar averages[41], whereas KO compare the
atomic-site potentials in a region far from the defect.[42]. Provided large
enough supercells, which allow for a sufficient sampling of the atomic-site
potentials, the latter method has been shown to provide better convergence,
especially in the case of ionic materials where defect-induced atomic
relaxations are large [42].

###  2.4  Equilibrium defect and carrier concentrations in the dilute limit

Equilibrium defect and carrier concentrations can be calculated in the dilute
limit once the relevant defect formation energies have been obtained. In the
dilute limit the defects do not interact and thus the energy required for
forming $n$ defects of type $D$ in charge state $q$ is given by: $\Delta
E_{d}(nD^{(q)})=n\Delta E_{d}(D^{(q)})$. A necessary condition for the
thermodynamic equilibrium is that the system grand potential, $\Phi$, is at an
extremum:

$\Phi(nD^{(q)})=\Phi(\mathrm{bulk})+\Delta E_{d}(nD^{(q)})+TS_{conf}(n)$ (8)

where $S_{conf}$ is the contribution of the configurational entropy to the
grand potential of the defect-laden system. Assume that there are
$g_{D^{(q)}}$ possible configurations in which $D^{(q)}$ has the same $\Delta
E_{d}(D^{(q)})$. Let $N$ be the number of unit cells forming the system and
$\gamma_{D^{(q)}}$ be the number of equivalent sites in the unit cell that the
defect can occupy, then (for $n\ll N\gamma_{D^{(q)}}$) the number of possible
ways to place $n$ non-interacting defects on $N\gamma_{D^{(q)}}$ sites is:

$\Omega_{D^{(q)}}=g_{D^{(q)}}\binom{N\gamma_{D^{(q)}}}{n}$ (9)

yielding the configurational entropy $S_{conf}=k_{B}\ln\Omega_{D^{(q)}}$.
Finding the extremum of equation (8) with respect to $n$, using Stirling’s
approximation for $S_{conf}$, gives the equilibrium defect concentration:

$c_{D^{(q)}}=\frac{n}{N}=\frac{\gamma_{D^{(q)}}g_{D^{(q)}}}{\exp{\left(\frac{\Delta
E_{d}(D^{(q)})}{k_{B}T}\right)}+g_{D^{(q)}}}$ (10)

which is usually approximated by the limit value for $\Delta
E_{d}({D^{(q)}})\gg k_{B}T$ as:

$c_{D^{(q)}}=g_{D^{(q)}}\gamma_{D^{(q)}}\exp{\left(-\frac{\Delta
E_{d}({D^{(q)}})}{k_{B}T}\right)}$ (11)

In case of more than one type of defect in the crystal, the equilibrium
concentration of each defect is still given by formula (10), assuming the
dilute-limit holds.

Equilibrium formation energies depend also on the electron chemical potential
$\mu_{e}$ which needs to be evaluated before computing equilibrium
concentrations. The equilibrium value of $\mu_{e}$ is fixed by the condition
that any actual solid will be characterized by a null net charge at
equilibrium. Spinney computes $\mu_{e}$ by finding the roots of the equation
describing the charge-neutrality condition:

$\sum qc_{D^{(q)}}(\mu_{e})+p_{0}(\mu_{e})-n_{0}(\mu_{e})=0,$ (12)

where $n_{0}$ is the concentration of free electrons:

$n_{0}=\int_{\epsilon_{C}}^{\infty}\frac{n(\epsilon)}{e^{(\epsilon-\mu_{e})/k_{B}T}+1}d\epsilon,$
(13)

and $p_{0}$ the concentration of free holes:

$p_{0}=\int_{-\infty}^{\epsilon_{V}}\frac{n(\epsilon)}{e^{(\mu_{e}-\epsilon)/k_{B}T}+1}d\epsilon.$
(14)

$n(\epsilon)$ is the density of states and $\epsilon_{V},\epsilon_{C}$ are the
eigenvalues of the valence band maximum and conduction band minimum,
respectively. Once the roots of equation (12) have been found, equilibrium
defect concentrations are obtained using equation (10) and carrier
concentrations using equations (13) and (14).

Often dopants are introduced in the material in conditions which are far from
the thermodynamic equilibrium assumed in the previous discussion, however, the
thermodynamic formalism is still generally used to assess the properties of
doped materials. If there is any indication that actual doping concentrations
will noticeably differ from those predicted by thermodynamic equilibrium,
Spinney allows to specify an effective doping concentrations
$N_{\mathrm{eff}}$ quantifying the amount of ionized dopant species. The
equilibrium electron chemical potential is then obtained by a modified version
of equation (12):

$\sum qc_{D^{(q)}}(\mu_{e})+p_{0}(\mu_{e})-n_{0}(\mu_{e})=N_{\mathrm{eff}}$
(15)

##  3 Implementation and Examples

###  3.1  General implementation features

Figure 1:  Typical workflow for first-principles calculations of point
defects. The properties which the Spinney code can calculate are included
within the dashed box.

Figure 1 illustrates the typical workflow for first-principles point defect
calculations in solids. According to the problem of interest, one selects the
most important native defects, eventual doping and the environmental
conditions most pertinent for the applications of the material. Different
environmental conditions entail different thermodynamic limits for the
chemical potentials of the atomic species forming the system. In order to
evaluate these limits, competing phases to the system under investigation must
be considered. Multiple online repositories offer large databases of chemical
compounds and the corresponding _ab-initio_ calculated electronic energy and
can be used for the identification of the most relevant phases. The most time-
consuming part of the whole process consists in performing the first-
principles calculations on the defect-containing supercells. Once the relevant
calculations have been completed some post-processing of the first-principles
data is required to obtain many of the energetic properties of the defect-
laden system. At this stage Spinney comes into play: first-principles data are
collected and fed to the appropriate routines which will output the quantities
necessary for calculating defect formation energies, thermodynamic charge
transition levels and equilibrium defects and carriers concentrations. The
modular design of the code allows to calculate each of the defect properties
independently. These results can then be used as input data for additional
steps in the pipeline.

Spinney is written entirely in Python 3 and a rapid execution speed is
achieved by extensive exploitation of NumPy arrays and the supported
vectorized operations [1]. As mentioned in the introduction, the basic
routines that allow for calculating the defect properties illustrated in
Figure 1 accept built-in Python’s data structures or NumPy arrays as input
data. This allows the users of general first-principles codes to parse the
data obtained in the calculations, format them and feed the result to the
appropriate routine in order to calculate the desired property. At the same
time Spinney implements higher-level routines which can automatically parse
the native output files of the popular DFT codes WIEN2k [29, 30] and VASP [27,
28]. For VASP, generally the OUTCAR and vasprun.xml files of the pristine and
defect containing supercells are sufficient. The FNV alignment scheme
furthermore requires the appropriate LOCPOT files with the calculated
electrostatic potentials. The size of these files can be considerable, adding
noticeable overhead which makes the correction scheme often I/O bound. For
WIEN2k, the files which are generally required are the case.struct and
case.scf. The KO alignment scheme also requires the case.vcoul files.
Additional information, which is brief, easily accessible and does not require
_ad hoc_ parsers, such as defect positions, values of the dielectric tensor
and valence band maximum, must be provided by the user. The main output of
Spinney’s basic routines are either scalars or tabular data which can be
accessed also as panda dataframes [3], allowing for seamless construction of
databases for point defect in solids. Defect properties are also summarized by
multiple plots enabled by the Matplotlib library [47].

###  3.2  Finding chemical potentials thermodynamic limits

Taking Nb-doped TiO2 anatase as an example, we will illustrate how to
investigate the limit values for atomic chemical potentials using Spinney. In
particular, the Spinney’s module used for this analysis is
thermodynamics.chempots.

The Ti-O system has a complex chemistry and many phases are know to exist.
Niobium doping of anatase is a promising method for obtaining alternative
transparent-conducting oxides [48] and it is also employed for improving the
photocatalytic properties of the material [49]. The inequalities and equality
constraints in equations (3) define the feasible region the three atomic
chemical potentials. In order to predict this region, we used the Materials
Project’s online repository [50], from which we obtained the first-principles
energies of more than 100 compounds in the Ti-O-Nb system. Formation energies
for all these compounds where then calculated using as reference state the HCP
structure for Ti, the BCC one for Nb and the triplet molecular state for O2.
As discussed in Ref. [51] the oxide materials formation energies were
calculated including the term correcting for the binding energy of the O2
molecule proposed in Ref. [52].

Case |  |  $\Delta\mu_{\text{Ti}}$ (eV) |  $\Delta\mu_{\text{O}}$ (eV) |  $\Delta\mu_{\text{Nb}}$ (eV)
---|---|---|---|---
1. | min | -10.33 | -4.46 | $-\infty$
max | -1.65 | -0.12 | 0
2. | min | -10.33 | -0.12 | $-\infty$
max | -10.33 | -0.12 | -10.35
3. | min | -1.78 | -4.46 | 0
max | -1.65 | -4.39 | 0
Table 1: Minimum and maximum values of $\Delta\mu$ calculated for the Ti-O-Nb
system. The first row (1.) takes into account the whole feasible region. The
second row (2.) considers the intersection of the feasible region with the
$\Delta\mu_{\text{O}}=\Delta\mu_{\text{O}}^{\mathrm{max}}$ plane. The third
row (3.) considers the intersection of the feasible region with the
$\Delta\mu_{\text{Nb}}=0$ plane.

Often defect properties, such as formation energies, are reported considering
only the extreme values of the atomic chemical potentials, for which the
system is in equilibrium with other phases. For example, assume one is
interested in the O-rich limit and desire to find out the possible values that
the chemical potentials $\Delta\mu_{\text{Ti}}$ and $\Delta\mu_{\text{Nb}}$
can achieve. Figure 2a shows the intersection of the feasible region with the
$\Delta\mu_{\text{O}}=\Delta\mu_{\text{O}}^{\mathrm{max}}$ plane, where
$\Delta\mu_{\text{O}}^{\mathrm{max}}$, which is reported in Table 1,
represents the maximum value that $\Delta\mu_{\text{O}}$ can achieve within
the feasible region described by equations (3). As Figure 2a and Table 1 show,
thermodynamic equilibrium requires that the maximum value
$\Delta\mu_{\text{Nb}}$ can achieve is $-10.35$ eV in the O-rich limit. For
higher values, Nb-doping would cause the segregation of $\mathrm{Nb_{2}O_{5}}$
in the $C2/m$ space group. Once $\Delta\mu_{\text{O}}$ has been fixed,
$\Delta\mu_{\text{Ti}}$ is fixed as well due to the equality constraint.
Figure 2b shows another extreme condition where the feasible region intersects
the $\Delta\mu_{\text{Nb}}=0$ plane. In this Nb-rich limit,
$\Delta\mu_{\text{Ti}}$ and $\Delta\mu_{\text{O}}$ can vary over a narrow
range, quantified in Table 1, along the line describing the equality
constraint. Outside this range, Ti2O and Ti3O5 would start to precipitate.

(a)  (b)

Figure 2:  Intersection of the feasible region of anatase TiO2:Nb with: (a)
the $\Delta\mu_{\text{O}}=\Delta\mu_{\text{O}}^{\mathrm{max}}=-0.12$ eV plane;
(b) the $\Delta\mu_{\text{Nb}}=0$ eV plane. The gray area represents the
inequalities constraints and the bold black line the equality constraint (see
equations (3c)).

###  3.3  Corrections for electrostatic finite-size effects

In this section we take as an illustrative example two intrinsic point defects
in two wide-band-gap materials: the B vacancy in the charge state -3,
$\square_{\mathrm{B}}^{-3}$, in cubic BN and the oxygen vacancy in the charge
state +2, $\square_{\mathrm{O}}^{+2}$, in ZnO. Spinney’s modules for
calculating electrostatic finite-size effects corrections employing the KO and
FNV schemes are located in spinney.defects.

Figure 3 compares the calculated defect formation energy as a function of the
supercell size for $\square_{\mathrm{B}}^{-3}$ in cubic BN (left-hand side)
and for $\square_{\mathrm{O}}^{+2}$ in ZnO (right-hand side) using the two
correction schemes for electrostatic effects implemented in Spinney. It can be
observed that both correction schemes predict the same value of the defect
formation energy for large enough supercells. In particular, Figure 3a shows
that for an isotropic system the predicted value of the defect formation
energy, after applying the correction term, converges to the limit value for
an infinite large supercell $E^{\prime}(\infty)$, described by the heuristic
equation: $E^{\prime}(N_{at})=E^{\prime}(\infty)+aN_{at}^{-1/3}+bN_{at}^{-1}$,
where $E^{\prime}(N_{at})$ represents the _uncorrected_ energy of a supercell
containing $N_{at}$ atoms [53].

(a)  (b)

Figure 3:  Relative defect formation energy as a function of the supercell size. (a) $\square_{\mathrm{B}}^{-3}$ in c-BN; (b) $\square_{\mathrm{O}}^{+2}$ in ZnO. The zero of the energy scale represents the convergent value of the defect formation energy. For $\square_{\mathrm{B}}^{-3}$ in c-BN the uncorrected defect formation energies are fitted as: $E^{\prime}(N_{at})=E^{\prime}(\infty)+aN_{at}^{-1/3}+bN_{at}^{-1}$. In the scheme of Freysoldt _et al._ we used a simple Gaussian charge density $N\exp(-r^{2})$ where $N$ is the normalization and $r$ is given in Å. System | Supercell |  $-E_{lat}$ (eV) |  $q\Delta\phi_{KO}$ (eV) |  $q\Delta\phi_{FNV}$ (eV)
---|---|---|---|---
| $2\times 2\times 2$ | 3.67 | 0.03 | 0.53
$\square_{\mathrm{B}}^{-3}$ | $3\times 3\times 3$ | 2.44 | 0.03 | 0.17
c-BN | $4\times 4\times 4$ | 1.83 | 0.09 | 0.08
| $5\times 5\times 5$ | 1.47 | 0.01 | 0.04
| $3\times 3\times 2$ | 0.80 | -0.23 | -0.25
$\square_{\mathrm{O}}^{+2}$ | $4\times 4\times 2$ | 0.66 | -0.11 | -0.14
ZnO | $5\times 5\times 3$ | 0.50 | -0.05 | -0.05
| $6\times 6\times 4$ | 0.40 | -0.02 | -0.02
Table 2: Contribution to the correction energy for electrostatic finite-size
effects calculated with the method of Kumagai and Oba and with the one of
Freysoldt _et al._ for two type of point defects. For c-BN, a supercell
expansion of the conventional cubic cell is used so that the $2\times 2\times
2$ cell contains 64 atoms. The $3\times 3\times 2$ ZnO cell contains 108
atoms. For $\square_{\mathrm{B}}^{-3}$ the B-rich, $\Delta\mu_{\mathrm{B}}=0$,
and for $\square_{\mathrm{O}}^{+2}$ the O-rich, $\Delta\mu_{\mathrm{O}}=0$,
limits were considered. All first principles calculations were performed using
the Perdew, Burke and Ernzerhof (PBE) [54] functional and the projector
augmented-wave (PAW) method [55]. For these two systems, $E_{lat}$ calculated
by the two methods differs in all case by at most 0.01 eV and is therefore
reported only once.

Table 2 compares the two terms entering equation (6) calculated in the two
correction schemes. The method of KO generally allows for a more robust
potential alignment procedure since it was found in Ref. [42] that the atomic-
site potentials are able to converge much faster far from the defect in ionic
materials; while this is not generally the case for the plane-averaged
electrostatic potential. Potential-vs-distance plots like those in Figure 4 or
Figure 1 of Ref. [42] represent therefore a valuable tool for assessing the
accuracy of the correction scheme and can be readily obtained with Spinney.

Figure 4:  Local electrostatic potentials used to compute the potential
alignment term, $\Delta\phi$, in the scheme of Kumagai and Oba, plotted with
respect to the distance from $\square_{\mathrm{B}}^{-3}$ in c-BN using a
$4\times 4\times 4$ supercell expansion. The atomic site potentials are used:
$V_{q/b}$ is the difference between the pristine and defect containing
supercell potentials obtained from the first-principles calculations,
$V_{PC,\,q}$ is the potential generated by the point charge model and $\Delta
V_{PC,\,q/b}$ is the difference between these two potentials. The latter term
is used to calculate $\Delta\phi$ and its value is well converged in the
sampling region far from the defect. The dashed horizontal line indicate the
zero of the potential.

We finally compared the defect formation energies of
$\square_{\mathrm{B}}^{-3}$ in c-BN calculated obtained with VASP and WIEN2k.
The computed values are in very good agreement between the two codes. The
difference in $\Delta E_{d}$ is about 0.2 to 0.3 eV for the $3\times 3\times
3$ and $4\times 4\times 4$ cells which amounts to about $2\%$.

###  3.4  Thermodynamic charge-transition levels

In these last two sections, Mg-doped GaN is taken as an illustrative example.
GaN is a very important semiconductor material which finds applications in
photodetectors, light emitting diodes (LEDs) in the blue and ultraviolet
region, laser diodes and bipolar transistors [56, 57, 58]. The material shows
an intrinsic n-type conductivity. On the other hand, for improving the
properties of GaN-based electronic devices such as LEDs and lasers the
synthesis of p-type GaN is desirable. Obtaining p-type GaN is generally
challenging but Mg is one of the most successful dopants used for this purpose
[59, 60].

This section reports thermodynamic charge transition levels for GaN:Mg. The
next section, section 3.5, will use these results for calculating defects and
carriers equilibrium concentrations. The employed Spinney’s modules are
located in spinney.defects. All calculation were performed using the VASP code
and PAW pseudopotentials with the PBE exchange-correlation functional. A
supercell containing 96 atoms was used and corrections for electrostatic
finite-size effects for charged defects were included using the KO scheme. All
intrinsic point defects (Ga and N vacancies, Ga and N intersitials, GaN and
NGa antisites) and the MgGa substitutional impurity were considered.

(a)  (b)

Figure 5:  Defect formation energies as a function of the Fermi level. The
white are represents the fundamental band gap predicted by the PBE functional.
This area is extended by the gray regions to the HSE fundamental gap, aligning
the valence band maxima of the two functionals. (a) Ga-rich conditions, (b)
N-rich conditions.

Figure 5 shows the calculated $\Delta E_{d}$ as a function of the Fermi level,
$E_{F}$, in Ga-rich ($\Delta\mu_{\mathrm{Ga}}=0$ eV, left-hand side) and
N-rich ($\Delta\mu_{\mathrm{N}}=0$ eV, right-hand side) conditions. Using the
same approach illustrated in section 3.2, we found that in both limits, the
atomic chemical potential values are constrained by the formation of
$\mathrm{Mg_{3}N_{2}}$. In the former limit, at equilibrium with
$\mathrm{Mg_{3}N_{2}}$, $\Delta\mu_{\text{N}}=-0.93$ eV and
$\Delta\mu_{\text{Mg}}=-0.67$ eV. In the latter limit, since
$\mathrm{Mg_{3}N_{2}}$ is still the competing phase,
$\Delta\mu_{\text{Ga}}=-0.93$ eV, while $\Delta\mu_{\text{Mg}}=-1.29$ eV.
These chemical potential values were used to calculate defect formation
energies in Figure 5. The calculated charge transition levels for selected
point defects in GaN:Mg are reported in Table 3. As charge transition levels
are defined for values of $E_{F}$ within the band gap, they are highly
affected by the well-known DFT band-gap error. It has been shown that
quantitative predictions of the charge transition levels in better agreement
with more accurate functionals, such as hybrid ones, can be obtained by
aligning the valence bad maximum predicted by local/semilocal functionals with
the one of the more accurate functional [61, 62]. Figure 5 displays such
alignment: the PBE band gap (white area) is extended (through the gray
regions) so that the conduction band maxima of the PBE functional and the
hybrid functional of Heyd, Scuseria and Ernzerhof (HSE) [63] are aligned. The
offset between the PBE and HSE valence band maxima has been taken from Ref.
[64], which considers the same system and the same functionals. While, after
such an alignment has been performed, PBE charge transition levels predicted
for $\square_{\mathrm{N}}$ agree discretely well with HSE calculations (_cf._
Ref. [64]), this is not the case for $\mathrm{Mg}_{\mathrm{Ga}}$ which is
predicted by PBE to be a quite deep acceptor, while HSE predicts it to be much
shallower [65]. This feature is not surprising as the
$\mathrm{Mg}_{\mathrm{Ga}}$ defect is an uncommon shallow acceptor, and HSE
calculations have found that it is characterized by a strongly localized hole
on a neighboring N atom [65], which cannot correctly be described by PBE due
to the well known self-interaction error.

Defect | $q/q^{\prime}$ |  $\epsilon(q/q^{\prime})$ (eV)
---|---|---
$\square_{\mathrm{N}}$ | 2/3 | 0.51
1/2 | 0.59
0/1 | 2.77
-1/0 | 3.41
$\mathrm{Mg}_{\mathrm{Ga}}$ | -1/0 | 1.03
Table 3: Thermodynamic charge transition levels of the most relevant defects
in GaN:Mg calculated using the PBE functional. The zero of the Fermi level is
set to the top of the HSE valence band.

###  3.5  Equilibrium defects and carriers concentrations

Once defect formation energies have been computed for the defects of interest,
equilibrium defect concentrations in the dilute limit can be calculated with
the formalism presented in section 2.4. Bulk GaN is usually growth at high
temperatures and the thermodynamic conditions can be described as Ga-rich due
to the very high nitrogen equilibrium pressure at these temperaratures [66].
Intrinsic GaN shows high electron concentrations in the range of 1017-1020
cm-3 [66, 67] and there is general consensus among experimental and
theoretical studies that they arise from the ionization of nitrogen vacancies
which are present at high concentration at the growth conditions (see Ref.
[66] and references therein).

Figure 6 shows the equilibrium concentrations of the N vacancy in Ga-rich
conditions calculated for intrinsic GaN considering a high-temperature range
representing the experimental growth conditions of bulk samples. From the
calculations, only $\square_{\mathrm{N}}$ in the illustrated charge states
assumes concentrations larger than 1010 cm-3 in the considered temperature
range. This indicates that the high concentrations of free electrons do indeed
stem from the ionization of N vacancies and, in particular, from singly
ionized donors, whose predicted concentration ranges from 1018 to
2$\times$1019 cm-3 in the temperature range between 1000 and 1500 K.

Figure 6:  Equilibrium concentrations of $\square_{\mathrm{N}}$ in intrinsic
GaN for high-temperature (1000-1500 K) and Ga-rich conditions.

Figure 7 shows the calculated carrier concentrations as a function of the
effective concentration of Mg doping, which represents the amount of activated
single-acceptor impurities at a given temperature (_cf._ with equation (15)).
Ga-rich conditions and a growth temperature of 1000 K were considered. The
concentration of the most relevant donor species, $\square_{\mathrm{N}}^{+}$,
is also shown. From the picture it is clear that the acceptor doping is
compensated by an increase of the intrinsic donor $\square_{\mathrm{N}}^{+}$
but hole concentrations larger than 1016 cm-3 can be obtained for an effective
Mg doping larger than 1019 cm-3. This model predicts a monotone increase in
the hole concentration as the amount of doping increases. In practice, the
dilute limit theory breaks down for large dopant concentrations and
segregation of Mg, with a decline of the p-type conductivity, has indeed been
observed for large Mg concentrations (see [68] and references therein).

Figure 7:  Equilibrium concentrations of carriers and
$\square_{\mathrm{N}}^{+}$ as a function of the effective concentration of Mg
calculated for Ga-rich conditions and a growth temperature of 1000 K. For the
carrier concentrations, negative values represent electrons and positive ones
holes.

##  4 Conclusion

We have presented Spinney a Python 3 package for post-processing first
principles calculations of point defects in semiconductors. Based on the
theory of solid state solutions in the dilute limit, the package is able to
calculate the most relevant energetic properties of point defects, including
formation energies and thermodynamic charge transition levels, and applying
state-of-the-art correction schemes for electrostatic finite-size effects. The
package can be used as a Python module, making it easy to integrate with other
computational frameworks. In this contribution we have shown in detail how the
package can be used to analyse and predict the properties of materials of
technological relevance.

##  5 Acknowledgements

M.A would like to thank Peter Blaha for the fruitful discussions. The authors
acknowledge support from the Austrian Science Funds (FWF) under project CODIS
(FWF-I-3576-N36). We also thank the Vienna Scientific Cluster for providing
the computational facilities (1523306: CODIS).

## References

  * [1] S. van der Walt, S. C. Colbert, G. Varoquaux, The numpy array: A structure for efficient numerical computation, Computing in Science & Engineering 13 (2011) 22–30.
  * [2] P. Virtanen, _et al._ , Scipy 1.0–fundamental algorithms for scientific computing in python (2019).  arXiv:1907.10121.
  * [3] W. McKinney, Data structures for statistical computing in python, Proceedings of the 9th Python in Science Conference (2010) 51–56.
  * [4] J. D. Hunter, Matplotlib: A 2d graphics environment, Computing in Science & Engineering 9 (2007) 90–95.
  * [5] A. H. Larsen, _et al._ , The atomic simulation environment—a python library for working with atoms, J. Phys.: Condens. Matter 29 (2017) 273002.
  * [6] S. Pizzini, Physical Chemistry of Semiconductor Materials and Processes, John Wiley & Sons, Ltd, 2015.  doi:10.1002/9781118514610.
  * [7] E. E. H. Matthew D. McCluskey, Dopants and Defects in Semiconductors, CRC Press, 2012.
  * [8] H. J. Queisser, E. E. Haller, Defects in semiconductors: Some fatal, some vital, Science 281 (5379) (1998) 945–950.  arXiv:https://science.sciencemag.org/content/281/5379/945.full.pdf, doi:10.1126/science.281.5379.945.
URL https://science.sciencemag.org/content/281/5379/945

  * [9] J. Maier, Nanoionics: ion transport and electrochemical storage in confined systems, Nature Materials 4 (11) (2005) 805–815.  doi:10.1038/nmat1513.
URL https://doi.org/10.1038/nmat1513

  * [10] J. A. del Alamo, Nanometre-scale electronics with iii-v compound semiconductors, Nature 479 (7373) (2011) 317–323.  doi:10.1038/nature10677.
URL https://doi.org/10.1038/nature10677

  * [11] X. Yu, T. J. Marks, A. Facchetti, Metal oxides for optoelectronic applications, Nature Materials 15 (4) (2016) 383–396.  doi:10.1038/nmat4599.
URL https://doi.org/10.1038/nmat4599

  * [12] D. B. Laks, C. G. Van de Walle, G. F. Neumark, S. T. Pantelides, Role of native defects in wide-band-gap semiconductors, Phys. Rev. Lett. 66 (1991) 648–651.  doi:10.1103/PhysRevLett.66.648.
URL https://link.aps.org/doi/10.1103/PhysRevLett.66.648

  * [13] A. Zunger, Practical doping principles, Applied Physics Letters 83 (1) (2003) 57–59.  arXiv:https://doi.org/10.1063/1.1584074, doi:10.1063/1.1584074.
URL https://doi.org/10.1063/1.1584074

  * [14] A. Walsh, A. Zunger, Instilling defect tolerance in new compounds, Nature Materials 16 (10) (2017) 964–967.  doi:10.1038/nmat4973.
URL https://doi.org/10.1038/nmat4973

  * [15] M. Setvin, C. Franchini, X. Hao, M. Schmid, A. Janotti, M. Kaltak, C. G. Van de Walle, G. Kresse, U. Diebold, Direct view at excess electrons in ${\mathrm{tio}}_{2}$ rutile and anatase, Phys. Rev. Lett. 113 (2014) 086402.  doi:10.1103/PhysRevLett.113.086402.
URL https://link.aps.org/doi/10.1103/PhysRevLett.113.086402

  * [16] M. Reticcioli, M. Setvin, X. Hao, P. Flauger, G. Kresse, M. Schmid, U. Diebold, C. Franchini, Polaron-driven surface reconstructions, Phys. Rev. X 7 (2017) 031053.  doi:10.1103/PhysRevX.7.031053.
URL https://link.aps.org/doi/10.1103/PhysRevX.7.031053

  * [17] P. Hohenberg, W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136 (1964) B864–B871.  doi:10.1103/PhysRev.136.B864.
URL https://link.aps.org/doi/10.1103/PhysRev.136.B864

  * [18] W. Kohn, L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140 (1965) A1133–A1138.  doi:10.1103/PhysRev.140.A1133.
URL https://link.aps.org/doi/10.1103/PhysRev.140.A1133

  * [19] C. G. Van de Walle, J. Neugebauer, First-principles calculations for defects and impurities: Applications to III-nitrides, Journal of Applied Physics 95 (8) (2004) 3851–3879.  doi:10.1063/1.1682673.
  * [20] D. A. Drabold, S. Estreicher, Theory of Defects in Semiconductors, Springer-Verlag Berlin Heidelberg, 2007.
  * [21] C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer, G. Kresse, A. Janotti, C. G. Van de Walle, First-principles calculations for point defects in solids, Rev. Mod. Phys. 86 (2014) 253–305.  doi:10.1103/RevModPhys.86.253.
URL https://link.aps.org/doi/10.1103/RevModPhys.86.253

  * [22] M. Leslie, N. J. Gillan, The energy and elastic dipole tensor of defects in ionic crystals calculated by the supercell method, Journal of Physics C: Solid State Physics 18 (5) (1985) 973–982.  doi:10.1088/0022-3719/18/5/005.
  * [23] G. Makov, M. C. Payne, Periodic boundary conditions in ab initio calculations, Phys. Rev. B 51 (1995) 4014–4022.  doi:10.1103/PhysRevB.51.4014.
URL https://link.aps.org/doi/10.1103/PhysRevB.51.4014

  * [24] E. Péan, J. Vidal, S. Jobic, C. Latouche, Presentation of the pydef post-treatment python software to compute publishable charts for defect energy formation, Chemical Physics Letters 671 (2017) 124 – 130.  doi:https://doi.org/10.1016/j.cplett.2017.01.001.
URL http://www.sciencedirect.com/science/article/pii/S0009261417300015

  * [25] A. Goyal, P. Gorai, H. Peng, S. Lany, V. Stevanović, A computational framework for automation of point defect calculations, Computational Materials Science 130 (2017) 1 – 9.  doi:https://doi.org/10.1016/j.commatsci.2016.12.040.
URL http://www.sciencedirect.com/science/article/pii/S0927025617300010

  * [26] D. Broberg, B. Medasani, N. E. Zimmermann, G. Yu, A. Canning, M. Haranczyk, M. Asta, G. Hautier, Pycdt: A python toolkit for modeling point defects in semiconductors and insulators, Computer Physics Communications 226 (2018) 165 – 179.  doi:https://doi.org/10.1016/j.cpc.2018.01.004.
URL http://www.sciencedirect.com/science/article/pii/S0010465518300079

  * [27] G. Kresse, J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B 47 (1993) 558–561.  doi:10.1103/PhysRevB.47.558.
URL https://link.aps.org/doi/10.1103/PhysRevB.47.558

  * [28] G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54 (1996) 11169–11186.  doi:10.1103/PhysRevB.54.11169.
URL https://link.aps.org/doi/10.1103/PhysRevB.54.11169

  * [29] P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka, J. Luitz, F. T. R. Laskowski, L. D. Marks, WIEN2k, An Augmented Plane Wave + Local Orbitals Program for Calculating Crystal Properties, Karlheinz Schwarz, Techn. Universität Wien, Austria, 2018.
  * [30] P. Blaha, K. Schwarz, F. Tran, R. Laskowski, G. K. H. Madsen, L. D. Marks, Wien2k: An apw+lo program for calculating the properties of solids, The Journal of Chemical Physics 152 (7) (2020) 074101.  doi:10.1063/1.5143061.
  * [31] S. B. Zhang, J. E. Northrup, Chemical potential dependence of defect formation energies in gaas: Application to ga self-diffusion, Phys. Rev. Lett. 67 (1991) 2339–2342.  doi:10.1103/PhysRevLett.67.2339.
URL https://link.aps.org/doi/10.1103/PhysRevLett.67.2339

  * [32] T. S. Bjørheim, M. Arrigoni, D. Gryaznov, E. Kotomin, J. Maier, Thermodynamic properties of neutral and charged oxygen vacancies in BaZrO3 based on first principles phonon calculations, Phys. Chem. Chem. Phys. 17 (2015) 20765–20774.  doi:10.1039/C5CP02529J.
URL http://dx.doi.org/10.1039/C5CP02529J

  * [33] T. S. Bjørheim, M. Arrigoni, S. W. Saeed, E. Kotomin, J. Maier, Surface segregation entropy of protons and oxygen vacancies in BaZrO3, Chemistry of Materials 28 (5) (2016) 1363–1368.  doi:10.1021/acs.chemmater.5b04327.
URL https://doi.org/10.1021/acs.chemmater.5b04327

  * [34] M. Arrigoni, T. S. Bjørheim, E. Kotomin, J. Maier, First principles study of confinement effects for oxygen vacancies in BaZrO3 (001) ultra-thin films, Phys. Chem. Chem. Phys. 18 (2016) 9902–9908.  doi:10.1039/C6CP00830E.
URL http://dx.doi.org/10.1039/C6CP00830E

  * [35] A. Glensk, B. Grabowski, T. Hickel, J. Neugebauer, Breakdown of the arrhenius law in describing vacancy formation energies: The importance of local anharmonicity revealed by ab initio thermodynamics, Phys. Rev. X 4 (2014) 011018.  doi:10.1103/PhysRevX.4.011018.
URL https://link.aps.org/doi/10.1103/PhysRevX.4.011018

  * [36] J. Malcolm W. Chase, NIST-JANAF thermochemical tables, Fourth edition. Washington, DC : American Chemical Society ; New York : American Institute of Physics for the National Institute of Standards and Technology, 1998., 1998, issued as: Journal of physical and chemical reference data; monograph no. 9, 1998.;Includes bibliographies.
URL https://search.library.wisc.edu/catalog/999842910902121

  * [37] P. A. Schultz, Local electrostatic moments and periodic boundary conditions, Phys. Rev. B 60 (1999) 1551–1554.  doi:10.1103/PhysRevB.60.1551.
URL https://link.aps.org/doi/10.1103/PhysRevB.60.1551

  * [38] P. A. Schultz, Charged local defects in extended systems, Phys. Rev. Lett. 84 (2000) 1942–1945.  doi:10.1103/PhysRevLett.84.1942.
URL https://link.aps.org/doi/10.1103/PhysRevLett.84.1942

  * [39] S. Lany, A. Zunger, Assessment of correction methods for the band-gap problem and for finite-size effects in supercell defect calculations: Case studies for zno and gaas, Phys. Rev. B 78 (2008) 235104.  doi:10.1103/PhysRevB.78.235104.
URL https://link.aps.org/doi/10.1103/PhysRevB.78.235104

  * [40] I. Dabo, B. Kozinsky, N. E. Singh-Miller, N. Marzari, Electrostatics in periodic boundary conditions and real-space corrections, Phys. Rev. B 77 (2008) 115139.  doi:10.1103/PhysRevB.77.115139.
URL https://link.aps.org/doi/10.1103/PhysRevB.77.115139

  * [41] C. Freysoldt, J. Neugebauer, C. G. Van de Walle, Fully ab initio finite-size corrections for charged-defect supercell calculations, Phys. Rev. Lett. 102 (2009) 016402.  doi:10.1103/PhysRevLett.102.016402.
URL https://link.aps.org/doi/10.1103/PhysRevLett.102.016402

  * [42] Y. Kumagai, F. Oba, Electrostatics-based finite-size corrections for first-principles point defect calculations, Phys. Rev. B 89 (2014) 195205.  doi:10.1103/PhysRevB.89.195205.
URL https://link.aps.org/doi/10.1103/PhysRevB.89.195205

  * [43] T. R. Durrant, S. T. Murphy, M. B. Watkins, A. L. Shluger, Relation between image charge and potential alignment corrections for charged defects in periodic boundary conditions, The Journal of Chemical Physics 149 (2) (2018) 024103.  arXiv:https://doi.org/10.1063/1.5029818, doi:10.1063/1.5029818.
URL https://doi.org/10.1063/1.5029818

  * [44] G. Fischerauer, Comments on real-space green’s function of an isolated point-charge in an unbounded anisotropic medium, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 44 (6) (1997) 1179–1180.  doi:10.1109/58.656617.
  * [45] R. Rurali, X. Cartoixà, Theory of defects in one-dimensional systems: Application to al-catalyzed si nanowires, Nano Letters 9 (3) (2009) 975–979.  doi:10.1021/nl802847p.
  * [46] H.-P. Komsa, T. T. Rantala, A. Pasquarello, Finite-size supercell correction schemes for charged defect calculations, Phys. Rev. B 86 (2012) 045112.  doi:10.1103/PhysRevB.86.045112.
  * [47] J. D. Hunter, Matplotlib: A 2d graphics environment, Computing in Science & Engineering 9 (3) (2007) 90–95.  doi:10.1109/MCSE.2007.55.
  * [48] Y. Furubayashi, T. Hitosugi, Y. Yamamoto, K. Inaba, G. Kinoda, Y. Hirose, T. Shimada, T. Hasegawa, A transparent metal: Nb-doped anatase tio2, Applied Physics Letters 86 (25) (2005) 252101.  arXiv:https://doi.org/10.1063/1.1949728, doi:10.1063/1.1949728.
URL https://doi.org/10.1063/1.1949728

  * [49] H. Su, Y.-T. Huang, Y.-H. Chang, P. Zhai, N. Y. Hau, P. C. H. Cheung, W.-T. Yeh, T.-C. Wei, S.-P. Feng, The synthesis of nb-doped tio2 nanoparticles for improved-performance dye sensitized solar cells, Electrochimica Acta 182 (2015) 230 – 237.  doi:https://doi.org/10.1016/j.electacta.2015.09.072.
URL http://www.sciencedirect.com/science/article/pii/S0013468615304849

  * [50] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, K. a. Persson, The Materials Project: A materials genome approach to accelerating materials innovation, APL Materials 1 (1) (2013) 011002.  doi:10.1063/1.4812323.
URL http://link.aip.org/link/AMPADS/v1/i1/p011002/s1&Agg=doi

  * [51] M. Arrigoni, G. K. H. Madsen, A comparative first-principles investigation on the defect chemistry of TiO2 anatase, The Journal of Chemical Physics 152 (4) (2020) 044110.  doi:10.1063/1.5138902.
  * [52] L. Wang, T. Maxisch, G. Ceder, Oxidation energies of transition metal oxides within the $\mathrm{GGA}+\mathrm{U}$ framework, Phys. Rev. B 73 (2006) 195107.  doi:10.1103/PhysRevB.73.195107.
URL https://link.aps.org/doi/10.1103/PhysRevB.73.195107

  * [53] C. W. M. Castleton, A. Höglund, S. Mirbt, Managing the supercell approximation for charged defects in semiconductors: Finite-size scaling, charge correction factors, the band-gap problem, and the ab initio dielectric constant, Phys. Rev. B 73 (2006) 035215.  doi:10.1103/PhysRevB.73.035215.
URL https://link.aps.org/doi/10.1103/PhysRevB.73.035215

  * [54] J. P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77 (1996) 3865–3868.  doi:10.1103/PhysRevLett.77.3865.
URL https://link.aps.org/doi/10.1103/PhysRevLett.77.3865

  * [55] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50 (1994) 17953–17979.  doi:10.1103/PhysRevB.50.17953.
URL https://link.aps.org/doi/10.1103/PhysRevB.50.17953

  * [56] O. Ambacher, Growth and applications of group III-nitrides, Journal of Physics D: Applied Physics 31 (20) (1998) 2653–2710.  doi:10.1088/0022-3727/31/20/001.
URL https://doi.org/10.1088%2F0022-3727%2F31%2F20%2F001

  * [57] H. Morkoç, S. Strite, G. B. Gao, M. E. Lin, B. Sverdlov, M. Burns, Large-band-gap sic, iii-v nitride, and ii-vi znse-based semiconductor device technologies, Journal of Applied Physics 76 (3) (1994) 1363–1398.  arXiv:https://doi.org/10.1063/1.358463, doi:10.1063/1.358463.
URL https://doi.org/10.1063/1.358463

  * [58] J. C. Johnson, H.-J. Choi, K. P. Knutsen, R. D. Schaller, P. Yang, R. J. Saykally, Single gallium nitride nanowire lasers, Nature Materials 1 (2) (2002) 106–110.  doi:10.1038/nmat728.
URL https://doi.org/10.1038/nmat728

  * [59] H. Amano, M. Kito, K. Hiramatsu, I. Akasaki, P-type conduction in mg-doped GaN treated with low-energy electron beam irradiation (LEEBI), Japanese Journal of Applied Physics 28 (Part 2, No. 12) (1989) L2112–L2114.  doi:10.1143/jjap.28.l2112.
URL https://doi.org/10.1143%2Fjjap.28.l2112

  * [60] I. Akasaki, H. Amano, M. Kito, K. Hiramatsu, Photoluminescence of mg-doped p-type gan and electroluminescence of gan p-n junction led, Journal of Luminescence 48-49 (1991) 666 – 670.  doi:https://doi.org/10.1016/0022-2313(91)90215-H.
URL http://www.sciencedirect.com/science/article/pii/002223139190215H

  * [61] A. Alkauskas, P. Broqvist, A. Pasquarello, Defect energy levels in density functional calculations: Alignment and band gap problem, Phys. Rev. Lett. 101 (2008) 046405.  doi:10.1103/PhysRevLett.101.046405.
URL https://link.aps.org/doi/10.1103/PhysRevLett.101.046405

  * [62] A. Alkauskas, A. Pasquarello, Band-edge problem in the theoretical determination of defect energy levels: The o vacancy in zno as a benchmark case, Phys. Rev. B 84 (2011) 125206.  doi:10.1103/PhysRevB.84.125206.
URL https://link.aps.org/doi/10.1103/PhysRevB.84.125206

  * [63] J. Heyd, G. E. Scuseria, M. Ernzerhof, Hybrid functionals based on a screened coulomb potential, The Journal of Chemical Physics 118 (18) (2003) 8207–8215.  doi:10.1063/1.1564060.
URL https://doi.org/10.1063/1.1564060

  * [64] J. L. Lyons, C. G. Van de Walle, Computationally predicted energies and properties of defects in gan, npj Computational Materials 3 (1) (2017) 12.  doi:10.1038/s41524-017-0014-2.
URL https://doi.org/10.1038/s41524-017-0014-2

  * [65] J. L. Lyons, A. Janotti, C. G. Van de Walle, Shallow versus deep nature of mg acceptors in nitride semiconductors, Phys. Rev. Lett. 108 (2012) 156403.  doi:10.1103/PhysRevLett.108.156403.
URL https://link.aps.org/doi/10.1103/PhysRevLett.108.156403

  * [66] P. Perlin, T. Suski, H. Teisseyre, M. Leszczynski, I. Grzegory, J. Jun, S. Porowski, P. Bogusławski, J. Bernholc, J. C. Chervin, A. Polian, T. D. Moustakas, Towards the identification of the dominant donor in gan, Phys. Rev. Lett. 75 (1995) 296–299.  doi:10.1103/PhysRevLett.75.296.
URL https://link.aps.org/doi/10.1103/PhysRevLett.75.296

  * [67] G. Kamler, J. Zachara, S. Podsiadło, L. Adamowicz, W. Gębicki, Bulk gan single-crystals growth, Journal of Crystal Growth 212 (1) (2000) 39 – 48\.  doi:https://doi.org/10.1016/S0022-0248(99)00890-8.
URL http://www.sciencedirect.com/science/article/pii/S0022024899008908

  * [68] P. P. Paskov, B. Monemar, 2 \- point defects in group-iii nitrides (2018) 27 – 61doi:https://doi.org/10.1016/B978-0-08-102053-1.00002-8.
URL http://www.sciencedirect.com/science/article/pii/B9780081020531000028