# Simulating Charged Defects at Database Scale

Jimmy-Xuan Shen  Materials Science Division, Lawrence Livermore National
Laboratory, Livermore, California 94550, U.S.A  Laboratory for Energy
Applications for the Future, Lawrence Livermore National Laboratory,
Livermore, California 94550, U.S.A  Lars Voss  Materials Engineering Division,
Lawrence Livermore National Laboratory, Livermore, California 94550, U.S.A
Joel B. Varley  Materials Science Division, Lawrence Livermore National
Laboratory, Livermore, California 94550, U.S.A  Laboratory for Energy
Applications for the Future, Lawrence Livermore National Laboratory,
Livermore, California 94550, U.S.A

###### Abstract

Point defects have a strong influence on the physical properties of materials,
often dominating the electronic and optical behavior in semiconductors and
insulators. The simulation and analysis of point defects is therefore crucial
for understanding the growth and operation of materials especially for
optoelectronics applications. In this work, we present a general-purpose
Python framework for the analysis of point defects in crystalline materials,
as well as a generalized workflow for their treatment with high-throughput
simulations. The distinguishing feature of our approach is an emphasis on a
unique, unitcell, structure-only, definition of point defects which decouples
the defect definition and the specific supercell representation used to
simulate the defect. This allows the results of first-principles calculations
to be aggregated into a database without extensive provenance information and
is a crucial step in building a persistent database of point defects that can
grow over time, a key component towards realizing the idea of a “defect
genome” that can yield more complex relationships governing the behavior of
defects in materials. We demonstrate several examples of the approach for
three technologically relevant materials and highlight current pitfalls that
must be considered when employing these methodologies, as well as their
potential solutions.

##  I Introduction

Understanding the physical characteristics of point defects (including
dopants) in semiconductors and insulators is a crucial analysis tool in
materials and semiconductor research since a small number of defects and
dopants can significantly change the electrical and optical properties of a
given material [1]. Since defects are usually formed during the initial growth
process, defect properties are difficult to isolate and experiment on. Hence,
first-principles calculations are often used to analyze these properties.
Specifically, calculations using density-functional theory (DFT) [2, 3] have
been widely adopted as the standard method. The analysis of point defects from
first principles usually begins with an analysis of the defect thermodynamics
by calculating the formation energy of the defect at a given set of chemical
conditions. These formation energies can elucidate which defects are likely to
form and different growth conditions and provide information on the doping
behavior of these defects [3]. Recently, more advanced analysis capabilities
of quantum recombination rates, for radiative [4], and non-radiative
recombination [5, 6], have been developed to help understand the interactions
of electrons and photons around point defects.

The computational cost for performing first-principles defect calculations is
unusually high due to a combination of factors. First, the periodic simulation
cell containing the defect must be large enough to avoid spurious interactions
between the atomic relaxation of the defect and its periodic images. Second,
since the defect represents a localized system embedded in a periodic crystal,
the system must be simulated using a level of theory that can accurately
describe the electronic structure of both the defect and the host crystal.
This usually means performing hybrid functional calculations in a plane-wave
basis set, which is notoriously 10 to 100 times more expensive than
traditional local or semi-local functional calculations.

Recently, automated defect analysis tools have been developed to facilitate
the simulation and analysis of defects in a high-throughput manner. Notable
among these frameworks are pyCDT [7], pyDefect [8], doped[9] and DASP [10].
pyCDT, pyDefect and doped are Python packages aimed at automating the
thermodynamic analysis of defects, while DASP also includes the simulation of
various quantum recombination rates.

While these tools accomplish their state goal of automating defect analysis,
their adoption by the broader defects research community has been limited and
building a large-scale persistent defects database is still out of reach. This
is a key limitation towards the realization of a curated “Defect Genome” that
can be refined over time and complement the more general Materials Genome
Initiative for accelerating materials engineering and discovery. [11] We
believe that significant improvements can be made in three areas that will
lead to broader adoption: (I) improved integration with existing workflow
automation and data aggregation tools, (II) conceptual separation between the
code for standard defect analysis and high-throughput calculations, and (III)
simulation-independent definition of defects which allows for calculations to
be grouped in a database without explicit provenance information. In this
work, we present the advances we have made in these areas and demonstrate how
a persistent database of defect properties can be built.

The core software driving our defect analysis is pymatgen-analysis-defects
[12], a part of the pymatgen [13] suite materials informatics. To address (I)
and (II), we developed pymatgen-analysis-defects to be a general-purpose tool
for the analysis of defects in crystalline materials and delegated all of the
code for high-throughput workflow automation to the atomate2 package built
upon the jobflow framework for automating general high-performance computing
tasks [14]. This allows us to codify common defect analysis tasks within
pymatgen-analysis-defects and provide a consistent interface for users to
perform these tasks without the added complication of high-throughput
calculations. For high-throughput calculations, we have developed a
generalized defect workflow as part of the atomate2 package which is capable
of orchestrating thousands of DFT defect calculations and organizing their
results in a standardized format.

Finally, for (III), we formalized a structure-only definition for point
defects which allows us to group calculation results without explicit
provenance information. A common critique of high-throughput defect
calculations is that structure minimization algorithms can miss the global
energy minimum, either due to a lack of local symmetry breaking [9] or not
finding the correct spin multiplicity [15]. Having a persistent definition
that is independent of the simulation details allows us to revisit previous
calculations and find lower energy configurations by using different
calculation settings or initial conditions. This allows us to update any
database incrementally as more calculations are performed, such as those that
increase the accuracy or sample additional structures. Over time, users can
improve both the breadth of their defects databases — by including
calculations for additional charge state, and the quality of data — by finding
relaxations to local or global energy minima that were previously missed by
the automated workflow.

##  II Background and Methods

The thermodynamic properties of a point defect ($X$) are controlled by the
defect formation energy ($E^{f}[X^{q}]$), which measures the energy cost of
creating a defect with a charge state $q$ at a given set of chemical
conditions. The formation energy is determined by the chemical potentials of
everything that goes into forming the defect [7, 3]:

$E^{f}[X^{q}]=E_{\rm tot}[X^{q}]-E_{\rm tot}[{\rm
bulk}]+\sum_{i}n_{i}\mu_{i}+qE_{\rm F}+\Delta^{q}\,,$ (1)

where $E_{\rm tot}[X^{q}]$ is the total energy of the defect simulation
supercell, $E_{\rm tot}[{\rm bulk}]$ is the total energy of the bulk cell of
the same size; $n_{i}$ is the change in the number of atoms of type $i$
required to form the defect, while $\mu_{i}$ is the chemical potential of atom
type $i$; $E_{\rm F}$ is the chemical potential of the electron (often called
the Fermi level), and $\Delta^{q}$ is the finite-size correction required to
account for the fact that we are simulating a charged defect in a finite
periodic cell.

While all of the different chemical potentials are interconnected, it is
conceptually easier to separate them into two categories: the chemical
potential of the atoms and the chemical potential of electrons. The chemical
potential of electrons (often called the Fermi level) accounts for how all of
the external conditions, including the presence of other defects, affect the
energy costs of adding or removing an electron from the defect system. Since
the electrons in the system can be manipulated after the defect is formed, the
Fermi level is often considered a free variable and is shown as the x-axis
whenever $E^{f}[X^{q}]$ is plotted for a given set of the other chemical
potentials. The chemical potentials of the different atoms added to or removed
from the defect are less dynamic and assumed to be fixed after the defect is
formed. So the formation energy of each charge state of a defect is typically
shown as a linear function of $E_{\rm F}$ with a slope of $q$ with a
$y$-intercept determined by the chemical environment during crystal growth.

The core goal of any defect analysis code is to calculate the formation energy
of a given defect in the different charge states. However, since $E_{\rm
tot}[X^{q}]$ is calculated by creating a point defect in a periodic supercell
and $\Delta^{q}$ presumably eliminates the finite size effects, we can have
many different supercell calculations that all correspond to the same defect
charge state. Traditionally, users will have to track the provenance of the
defect and charge state for all the calculations. This can quickly become
untenable as the number of defects and charge states increases. Additionally,
since the most stable spin and atomic configuration of a defect might not be
found by the first structure relaxation, it is often necessary to revisit
defect calculations to improve upon the data over time. This necessitates a
way for us to compare and group defects without explicit provenance
information. Toward that end, we have developed a structure-only definition of
point defects that will allow persistent, incrementally growing databases of
defect simulation results to be constructed (cf. Sec. II.1). To compute the
different atomic chemical potentials ($\mu_{i}$), we combine explicit DFT
calculations of the pure elemental phases and experimentally corrected values
of the formation enthalpy to arrive at a more accurate measure of the
elemental chemical potentials with fewer calculations.

For the calculation of the finite-size correction term ($\Delta^{q}$), we use
the finite-size correction method of Freysoldt, Neugebauer, and Van de Walle
(FNV) [16]. Previous implementations [16, 7] of this method required the
position of the defect to be supplied by the user. Since we require our data
to be provenance-agnostic, we have developed a method to calculate the
position of a point defect without prior knowledge of how the defect was
defined (cf. Sec. II.3). This eliminates a persistent point of user
intervention in the workflow of defect analysis and allows for more software-
driven high-throughput approaches.

The combinations of these advances allow us to arrive at an analysis framework
to define, simulate, and analyze defects in a fully automated manner. All of
the code responsible for the definition and analysis of point defects is
available via the pymatgen-analysis-defects add-on package to pymatgen. [17]
The automation workflow is distributed independently as part of the atomate2
package which allows us to orchestrate and organize thousands of quantum
chemistry calculations (cf. II.4). We will detail these advancements in the
remainder of this section.

###  II.1 Structure-only Definition of Defect

For most DFT calculations, the physical system you are simulating is solely
defined by the atomic positions and the lattice vectors of the unit cell. This
allows large databases of DFT calculations to be built by simply storing the
atomic structures and aggregating all of the calculations related to a given
structure into a single entry. While the unit cell is not unique, modern
symmetry analysis tools allow different representations of the same periodic
structure to be matched to each other so it is trivial to group calculations
that represent the same physical system. However, when simulating point
defects, we are using an arbitrary supercell to simulate the defect in the
dilute limit, so multiple symmetry-distinct supercell structures can all
correspond to the same point defect. This inherently breaks the one-to-one
mapping between structures and calculations, and thus the structure matching
along is not enough facilitate the building of a persistent database of defect
calculations. A core feature of pymatgen-analysis-defects is the ability to
define defects in a structure-only manner, which allows us to group
calculations for the same defect without explicit provenance information.
Ultimately, this definition frees us from the need to explicitly track the
relationship and dependencies between sets of defect calculations and makes
database building significantly easier.

In pymatgen-analysis-defects, a defect is defined by the combination of:

  * •

The type of defect: whether it is a vacancy, interstitial, or substitutional
defect.

  * •

A unique primitive unit cell of the bulk material.

  * •

The fractional position of the defect in the unit cell.

This definition allows us to create a representative structure
($\mathcal{S}_{X}$) of the defect $X$ in the unit cell. As long as the
primitive cell is fixed, two defects $X$ and $Y$ are considered to be the same
if and only if $\mathcal{S}_{X}$ and $\mathcal{S}_{Y}$ can be structurally
matched to each other, i.e. they are equivalent under the symmetry operations
of the primitive cell. While there might be multiple equivalent sites in a
primitive cell for a point defect to form, the structures representing these
different versions will be symmetrically equivalent to each other. As
examples, all of the symmetry-distinct native defects of wurtzite GaN are
shown in Fig. 1. The vacancy and antisite defects are trivially constructed,
however, generating interstitial defects has always been challenging.
Following the method from Ref. 18, the symmetry-distinct (as defined using
Spglib [19]) local minima in the charge density were identified as possible
interstitial insertion sites, which resulted in two symmetry-distinct versions
of the Ga and N interstitial defects. Note that this method and be easily
combined with charge density API [20] from the Materials to quickly generate
defect structures without additional DFT calculations.

Figure 1: Defect structures for native point defects in GaN, note that two
symmetry distinct initial interstitial sites were identified.

These representative defect structures can be used to define additional
quantities such as predicted oxidation states of the elements and degeneracy
factors which are required to calculate various thermodynamic quantities
related to point defects [3]. Additionally, once a defect has been created, we
can generate a nearly cubic supercell structure containing that point defect
to minimize the effects of periodic images. The resulting calculations can be
associated with the original defect so the thermodynamic quantities from the
DFT calculations amount to one representation of a concrete concept such as
“energy of VGa in the $q=-1$ charge state”. For the 4-atom primitive cell of
GaN, we can create a 128-atom supercell that has orthogonal lattice vectors
that differ in length by less than 10%. For the remainder of this work,
explicit supercell defect DFT calculations will be performed using this GaN
supercell structure.

###  II.2 Chemical Potential Corrections Based on Experimental Data

The chemical potential is a measure of the thermodynamic cost of adding or
removing an atom from the defect when the formations of different competing
phases are considered. If the formation of the defect is only competing with
the pure elemental phases, then the chemical potential is simply the energy
per atom from the DFT calculation of the elemental phase ($\mu^{\rm
DFT}_{i}$). However, competing compounds with lower formation enthalpy can
lower the chemical potential ($\mu_{i}=\mu^{\rm DFT}_{i}+\Delta\mu_{i}$) since
they create a more stable environment for the defect. For a set of atomic
species (indexed by $i$), and a set of competing compounds (indexed by
$\alpha$) these correction terms must satisfy the following constraints:

$\displaystyle\Delta\mu_{i}$ $\displaystyle\leq 0$ $\displaystyle{}\forall{}$
$\displaystyle\text{ elements }i\,$
$\displaystyle\sum_{\alpha}n^{\alpha}_{i}\Delta\mu_{i}$ $\displaystyle<\Delta
H_{f}(\alpha)\qquad$ $\displaystyle{}\forall{}$ $\displaystyle\text{ compounds
}\alpha$ (2) $\displaystyle\sum_{\alpha}n^{\rm bulk}_{i}\Delta\mu_{i}$
$\displaystyle=\Delta H_{f}({\rm bulk})$

Where $n^{\alpha}_{i}$ is the number of atoms of species $i$ in a formula unit
of compound $\alpha$, and the $\Delta H_{f}(\alpha)$’s are the formation
enthalpies of one formula unit of the compound $\alpha$. Eq. II.2 defines a
manifold in chemical-potential space where you are constrained to the facet
representing the host material but bound on different sides by the formation
of competing phases. Since the total energy from DFT can change under
different calculation settings, we are usually forced to calculate all of the
formation enthalpies using the same level of theory. However, since the
formation energy data from the Materials Project is specifically corrected to
fit experimental formation enthalpies [21], we can use the formation enthalpy
data from the Materials Project to calculate the chemical potential
corrections ($\Delta\mu_{i}$) only and expect similar if not better results.
Hence, the chemical potential in Eq. 1 can be expressed as:

$\mu_{i}=\mu^{\rm DFT}_{i}+\Delta\mu^{\rm MP}_{i}$ (3)

where we combine the explicitly calculated energy of the elemental phase in a
given level of theory with the formation enthalpy correction due to competing
phases from the Materials Project. This approach is distinct from other
approaches based on correcting inaccuracies in chemical potential references,
such as fitted elemental-phase reference energies (FERE) [22], which instead
corrects the $\mu^{\rm DFT}_{i}$ terms based on a consistent calculation
framework and thermodynamic databases from which to derive the corrections.
From a database-building perspective, this simultaneously provides a more
accurate prediction of relative stability between the competing phases (since
we are using experimentally-fitted formation enthalpies) and allows users to
build a database using only the DFT energy from the elemental phases alone,
which significantly reduces the computational cost.

###  II.3 Identify Defect Position With No Prior Knowledge

The total energy of the defective supercell ($E^{f}[X^{q}]$ in Eq. 1) refers
to the energy of the defect simulation supercell after atomic relaxation.
Since atomic relaxation is often unpredictable, it is difficult to know the
exact position of the defect in the supercell afterward. This is especially
problematic when we attempt to calculate the finite-size charge correction for
native defects in a supercell since the position of the defect is an input
that helps determine the potential alignment contribution [16]. Usually, one
can arrive at a decent estimate of the defect position by keeping track of the
position during the defect-creating processes. But for cases where the
relaxation of the atomic positions is severe, such as the case of split
interstitials and split vacancies, the relaxed defect position will differ
significantly from the defect creation site. Additionally, as we aim to build
databases without explicit provenance information, we need to find a way to
identify the position of the defect in the supercell after structure
relaxation without any prior knowledge of the defect creation site. In
pymatgen-analysis-defects, we included a method to identify the position of
the defect in the supercell after structure relaxation by calculating a site-
specific distortion field based on SOAP vectors [23, 24]. For each site in the
relaxed supercell, we calculate the SOAP vectors ($\bm{v}^{\rm
defect}_{\sigma}$) for each site $\sigma$ and compare them with the SOAP
vectors of the sites in the pristine bulk structure ($\bm{v}^{\rm
bulk}_{\sigma^{\prime}}$). The distortion on site $\sigma$ is defined as:

$\delta_{\sigma}=1-\max_{\sigma^{\prime}}\left\\{\frac{\bm{v}^{\rm
defect}_{\sigma}\cdot\bm{v}^{\rm bulk}_{\sigma^{\prime}}}{\left|\bm{v}^{\rm
defect}_{\sigma}\right|\left|\bm{v}^{\rm
bulk}_{\sigma^{\prime}}\right|}\right\\}$ (4)

To identify the set of sites with the largest distortions, we first sort the
distortion values ($\delta_{\sigma}$) in descending order and cut off the
descending series at the largest drop in distortion value. The position of the
defect is then given by taking the average position of the most distorted
sites, weighted by their distortion value.

Using this site finder algorithm, we can calculate the FNV finite-size
correction with only the electrostatic potentials from the DFT calculations
and the dielectric constant of the material. The alignment-like term in the
original Freysoldt correction paper requires the short-range potential
($V^{\rm sr}_{q/0}$) to be shifted by a constant value ($C$) so that the
potential is effectively zero far away from the defect. Fig. 2 shows the
planar-averaged values of the different potentials involved in calculating the
FNV charge-state correction. The difference in electrostatic potentials
between the defect and pristine cells (black curve) and the short-range
potential (green curve) are both directly computed from the files storing the
electrostatic potential information (e.g. LOCPOT files for VASP calculations).
The position of the defect was automatically determined using the finder
algorithm and can be validated by the fact that the potential difference and
short-range potential profile both peak at the origin.

Figure 2:  Planar-averaged potentials used during the calculations of the FNV
charge-state correction for the $V_{\rm Ga}$ defect in GaN. The long-range
(black) term is computed from a Gaussian model, the potential difference (red)
and short-range terms are both shifted so the origin is at the defect
positions identified by the finder algorithm. The sampling region (shaded
red), over which the $C$ offset (dashed) is computed in a small region far
from the origin in each crystal direction.

This method removed the usual need for user input from previous methods [16,
7] and allow for the entire process of defect formation energy calculation to
be codified within our Python framework. This makes the investigation of
singular defects much more approachable and significantly reduces the book-
keeping complexity incurred in other high-throughput defect analysis
approaches. Additionally, the ability to easily identify the positions of
native vacancies and interstitials in structures using only the atomic
positions is likely to have uses beyond the present context.

###  II.4 Atomate2 Workflow for Defect Formation Energies

The release of atomate [25] and aiida [26] marked a major milestone in the
open-source materials workflow development, as they codified many standard
materials science workflow under a well-defined automation framework. However,
complex workflows such as defect formation energy calculation were difficult
to develop due to the fundamental limitations such as processing of the
volumetric data and the ability to dynamically create additional DFT
calculations based on the results of previous calculations. atomate2 directly
addresses these challenges and allows us to define dynamic and composable
workflows that can also natively handle the storage of large volumetric data.
With these advances, we were able to develop a flexible and composable high-
throughput defect simulation framework that is capable of handling a variety
of different simulation needs.

The flowchart for the workflow to calculate the defect formation energy of a
single defect is shown in Fig. 3. The two primary inputs for the workflow are
the uniquely defined Defect object (as discussed in Sec. II.1) and the Relaxer
construct that contains all of the specific calculation settings for the
relaxations procedure. Because we have abstracted away the computational
details of the charge state relaxation using the Relaxer, any functional
support by VASP [27] or even entirely different DFT simulation codes (such as
CP2K [28]) can be used with minimal modification to the code.

In the workflow for a single defect, we first create a supercell of the bulk
material and relax the lattice parameters and atomic positions. Then, we
create a defect supercell with the same lattice parameters as the bulk
supercell and relax only the atomic positions of each charge state. The charge
states in this case, are derived from a combination of the allowed formal
oxidation states of the species in the bulk formula, and the abundances of
specific charge states for these species in the ICSD [29]. Once these
individual charge state relaxations are completed, their raw results are
stored in a database along with a serialized version of the original defect
object used to generate these calculations. This allows the final data
aggregation step to only examine the Defect objects and the individual
calculation settings to quickly find the set of calculations corresponding to
charge states of a particular defect simulated using a given level of theory.
The final step is to combine the computed energies and electrostatic potential
data from the bulk calculation and different charge state calculations to
arrive at a prediction of the formation energy of the defect. Since we have
designed our framework specifically to eliminate the need to track data
provenance, this final step can be performed independently from the simulation
process. This becomes more important as a database scales and bookkeeping for
the provenance of specific calculations becomes more difficult and mistake-
prone.

Figure 3:  Flowchart for the workflow to calculate the defect formation
energy. The Relaxer and Defect objects are the required inputs. A single bulk
supercell calculation is required to obtain the electrostatic potential (e.g.
a LOCPOT file for VASP calculations) data and to determine the supercell
lattice parameters. The outputs from the different charge state supercell
calculations can be combined with the bulk calculation to compute the defect
formation energy.

The composable nature of atomate2 allows us to connect many defect workflows
to avoid repeating the Bulk supercell relaxation step. We can perform this
step once and initialize the Spawn charge state calculations step without any
additional DFT calculations.

##  III Results & Discussions

Using the framework described in Sec. II, we have calculated the formation
energies of all of the native defects in GaN at different levels of theory and
under different atomic chemical conditions. The formation energy diagrams for
native defects of GaN, computed using HSE06 [30] and PBEsol functionals are
shown in Fig. 4. We elected not to fit the exact-exchange fraction to the band
gap and used a default of 25% to make the data directly comparable between
material systems. Nonetheless, the results of our HSE06 calculations match
well with previously reported results with a different exact exchange that
better represents the experimental band gap of GaN [31] with one notable
exception of the NGa defect. For NGa, the formation energy of the $q=0$ charge
state in HSE06 was too high to be stable relative to the $q=-1$ and $q=+1$
charge states and the estimated charge state range did not agree with the
charge states found by Lyons et al. [31]. We use this example to illustrate
the first shortcoming in typical high-throughput defect approaches that
generate structures that may adopt insufficiently perturbed initial structures
or an incomplete set of defect charge states (based on the constituents’
oxidation states) that are explored. The first problem can be fixed by
performing the HSE06 calculation with larger distortions to the environment
around the defect, which lowered the formation energy of the $q=0$ charge
state by 0.58 eV and stabilized the neutral charge state in $n$-type
conditions. The second problem is more subtle and caused by the unusual
oxidation behavior of this particular defect. Since Ga is found in the $+3$
oxidation state and N is found in the $-3$ oxidation state the relevant charge
states for the NGa defect are expected to range from $-6$ to $0$, however,
Lyons et al. [31] found multiple positive charge states for this defect to be
stable. This represents one example of how “standard" types of point defects
may be interpreted as defect complexes for particular oxidation states, e.g.
$V_{\rm Ga}$-Ni complexes that can adopt higher positive charge states for
Fermi levels closer to the VBM. We note that this type of consideration can be
factored in a priori for the prototypical point defect types considered here.
While formal treatment of defect complexes, as well as the integration with
informed symmetry-breaking algorithms [9] are fully supported by the pymatgen-
analysis-defects framework, they should only be considered on a case-by-case
basis due to the increased computational cost. While these issues only
affected a single high-formation-energy defect in those studied in GaN, they
highlight the general difficulty in automating defect calculations and why our
approach to incrementally updating the data is required.

We compared our PBEsol and HSE06 calculations and discovered two general
methods to improve the overall speed and reliability of defect calculations.
First, we found that the use of a PBEsol calculation to precondition the HSE06
relaxation can significantly reduce the computational cost of the defect
calculations. However, we acknowledge that this must be considered with care
for certain charge states, as a bias toward charge delocalization from
semilocal functionals can lead to relaxed structures in local minima which may
require additional perturbations or relaxation steps to access global minima
for a given defect configuration and charge state [9]. Second, we found that
constraining the atoms far from the defect to their bulk positions can both
reduce the computational cost and improve the reliability of the defect
calculations. Both of these approaches have been implemented as optional
settings that can be easily toggled on or off by the user.

Figure 4:  Formation energy diagrams for native defects in GaN computed using
PBEsol (a) and HSE06 (b). Both sets for formation energy diagrams are
calculated using Ga-rich chemical conditions.

We recommend performing HSE06 defect calculations using a Relaxer that
consists of a fast PBEsol relaxation to precondition the structure and
wavefunctions, followed by a slower HSE06 relaxation which typically produces
more accurate descriptions of the localized defect states and energies. We
again stress that the extent to which the defect structures are relaxed with
(semi)local functionals before switching to hybrid functionals is sensitive to
the particular oxidation state and symmetry of the defect, with electronic
closed-shell defects being much more amenable to extended relaxations with
semilocal functionals, such as the $V_{\rm Ga}^{-3}$. For 8 native defects in
GaN, we identified a total of 70 distinct charge state calculations that
needed to be performed based on our analysis of the formal oxidation state of
the defects. Each electronic optimization step of an HSE06 calculation is
usually around 10 to 200 times more expensive than the corresponding PBEsol
calculations so we can simply count to number of HSE06 electronic steps taken
to determine the relative computational cost of the different approaches.
Using only HSE06 to perform the atomic relaxations, we needed an average of
132.3 electronic steps to complete each charge state calculation. Using PBEsol
atomic to precondition the HSE06 relaxation, we needed an average of 84.2
HSE06 electronic steps to complete each charge state calculation. This
represents a 36% reduction in the number of electronic steps required with a
less than 0.05 eV difference in the final formation energy predictions.

A well-known problem when performing defect calculations is the fact that the
perturbation of the crystal caused by the defect can cause another phase
transition in the host materials so no matter how large the supercell is, the
defect distortion caused by the defect will still cause periodic image
interactions. This can be especially problematic for interstitial defects
which can cause large distortions in the surrounding lattice. These
distortions are often not physically meaningful since the bulk phase is
stabilized by finite temperature effects not captured by the ground-state DFT
calculations. To avoid this problem, we first calculate the largest sphere
that can be inscribed in the supercell and then constrain all atoms outside of
this sphere to their bulk positions. This ensures that the lattice distortions
caused by the defect cannot interact with itself across periodic boundaries.

Using the same methods for generating and analyzing native defects in GaN, we
have also computed the formation energies of the native defects in
$\beta$-Ga2O3, a lower-symmetry, more complicated monoclinic structure that
has two symmetry-distinct Ga sites and three symmetry-distinct O sites. This
leads to 10 distinct “simple” native defects in Ga2O3, excluding
interstitials: two each for the $V_{\rm Ga}$ and OGa defects and three each
for the $V_{\rm O}$ and GaO defects. with a number of interstitial sites also
possible in the $\beta$-gallia structure. Using the same charge-density method
as above [18], we identified 5 symmetry distinct interstitial sites by
analyzing the electronic charge density and they coincide with recently
published work on interstitials in Ga2O3. [32]. Since the number of candidate
interstitial sites can be large, we recommend limiting the insertions to the
two sites with the lowest average charge density in the surrounding region
since that is known to be a good predictor of stability. [33] The formation
energy diagrams for native defects in Ga2O3 are again computed with two
different levels of theoretical treatment, this time with the PBEsol and HSE06
functionals and shown in Fig. 5. As for GaN, we used a standard fraction of
exact exchange of 25% for HSE, which is also known to underestimate the band
gap of $\beta$-Ga2O3. Once again, our fully automated framework largely
replicated recently published computational results [34] without any human
intervention. We note that the qualitative orderings of the defects are
consistent with other results, and summarize the computed transition levels in
the Supporting Information. However, it has been shown that $\beta$-Ga2O3
exhibits several highly distorted vacancy configurations (so called split-
vacancies) that are lower in energy and that the automated relaxation
framework was unable to initially capture. [35, 34, 36] This echoes a key
current deficiency in point defect workflows similar to the case with GaN, and
highlights the importance of efficient databasing and the ability to apply
additional structural perturbations that may be outside the typical tolerance
of automated symmetry-breaking algorithms.

Figure 5:  Formation energy diagrams for native defects in Ga2O3, computed
using PBEsol (a) and HSE06 (b). Both sets for formation energy diagrams are
calculated using Ga-rich chemical conditions.

As a third case study, we also computed the formation energies of the
different native defects of cubic SrTiO3 (STO). While the cubic perovskite
phase is structurally simple, STO thermodynamically prefers a lower-symmetry
structure at 0K and also exhibits a more complex chemical stability region
given the larger number of elemental constituents. The perturbation caused by
the defect can lead to spurious transformations to the cubic phase when the
defect is simulated in a finite-sized periodic cell. To alleviate this, we
define a sphere centered at the defect with the largest radius allowed by the
periodic boundary conditions and freeze all atoms outside this sphere. This
leads to faster convergence of the structural relaxation and more reliable
results.

The atomic structure of STO only has a single symmetry-distinct site for each
atomic species, resulting in three distinct vacancies and six distinct
antisite defects. Two symmetry-distinct interstitial sites were identified
using the electronic charge density, which resulted in another six
interstitial native point defects. Since the Sr-Ti-O chemical space is more
complex, the benefit of using thermodynamic databases like the Materials
Project to determine the relevant boundaries for the chemical potentials is
more apparent. A total of nine distinct compounds are involved in defining the
stability region of STO in chemical potential pace, which would usually
require nine additional hybrid DFT calculations on top of the calculations for
the total energies of the elemental phases. Using the method outlined in Sec.
II.2, we can skip these nine calculations and use the formation enthalpy data
from the Materials Project to compute the chemical potential limits. The
stability region for STO in chemical potential space is shown in Fig. 6 (a),
and we have selected the Ti-rich and O-rich growth conditions to represent two
limiting conditions for the growth of STO. This contrasts previous studies
[37, 38] which only considered TiO2 and SrO as the competing phases and only
considered a single interstitial site.

Figure 6:  (a) The stability region of STO in the (Sr, Ti) chemical potential
space, the Ti-rich (teal), and O-rich (magenta) growth conditions are
indicated. (b-d) The formation energy diagrams of native point defects of STO
under Ti-rich conditions. (e-g) The formation energy diagrams of native point
defects of STO under O-rich conditions.

The formation energy diagrams of native defects of STO under O-rich conditions
are shown in Fig. 6 (b-d) and the ones for Ti-rich conditions are shown in
Fig. 6 (e-g). Similar to Ref. 37 and Ref. 38, we find that the three vacancies
$V_{\rm O}$, $V_{\rm Ti}$ and $V_{\rm Sr}$ are the most stable native defects
under most conditions. Under O-rich conditions, the +2 charge state of $V_{\rm
O}$ is the most stable defect across most Fermi levels in the band gap. Under
Ti-rich conditions, the $-2$ charge state of $V_{\rm Sr}$ becomes more stable
for $n$-type conditions, and the $-4$ charge states of $V_{\rm Ti}$ become
stable just below the conduction band. The relevant defects, charge states,
and formation energies of these defects change under different chemical
conditions are all consistent with previous studies [37, 38] with errors
typical for defect calculations. [39, 40, 41]

As we have seen from the examples above using GaN, Ga2O3, and STO, our fully
automated framework can reproduce previous computational results without any
human intervention. Instead of a static database of the results at the end,
our strategy results in a dynamic database that can be incrementally improved
over time as more accurate calculations are performed. In each example above,
we were able to identify a set of defects that are likely to be relevant under
a given set of chemical conditions and have a reasonable estimate of the
charge state behavior of these defects.

A limitation of any high-throughput approach is the trade-off between
generality and accuracy. Here, we have elected to use a single mixing
parameter of 25% for all calculations, which is a choice that simplifies the
comparison of energetics within a given level of theory and computational
parameters. However, a fixed mixing parameter prevents us from using the
common strategy for tuning the mixing parameter to match the experimental band
gap, a typical choice for more direct comparisons of calculated defect
transition levels with experiments. Choices of “incorrect" exact exchange will
result in some inaccuracies in the formation energies and defect transition
levels that are largely predictable and quantifiable given that localization
is qualitatively consistent between the calculations, and unlikely to affect
the overall trends in the defect behavior. [42] This re-emphasizes the need to
separate the responsibilities of the defect analysis and high-throughput
defect simulation since high-throughput defect analysis is often more
concerned with the relative stability of the defects and not the exact
thermodynamic properties of a specific defect. The data shown in Figs. (4-6)
are available as tables in the Supplementary Material.

##  IV Conclusion

We have demonstrated a general-purpose codebase (pymatgen-analysis-defects)
for the analysis of point defects in crystalline materials. Our code enforces
a strong connection between a symmetry-distinct defect structure and core
concepts in defect analysis, allowing for further extension of the code base
in an object-oriented fashion. Although not discussed in the manuscript, we
have also implemented methods for configuration-coordinate diagram analysis as
well as tools for simulating radiative and non-radiative recombination at
defect centers. This framework codifies some of the more complex computational
workflows in materials informatics and the strong integration with existing
analysis and automation software allows for more consistent development and
maintenance by the community as the core pymatgen codebase has demonstrated
over the past decade.

Since the code is designed for minimal user intervention, we were able to
integrate this with existing workflow automation tools to build a fully
automated high-throughput defect analysis framework. We have demonstrated that
a fully automated approach can largely reproduce previously published
computational results without any human intervention and identified current
drawbacks with the current implementation. Specifically, we find that
undersampling of global and local minima structures of particular oxidation
states, particularly in the case of defect complexes or other defects that
exhibit large structural distortions, may be missed in initial searches but
can be readily revisited with our emphasis on a robust database-centered
representation of point defects. These issues are currently being resolved in
future versions that will handle symmetry-breaking and defect complexes in a
more rigorous manner. With more computational resources, it is now possible to
build a persistent database of defect properties that can be improved over
time as more accurate calculations are performed, although that is beyond the
scope of the present work.

##  V Supplementary Material

The transition levels for GaN, Ga2O3 and SrTiO3 are presented as tables in the
Supplementary Material.

##  VI Acknowledgements

This work was performed under the auspices of the US DOE by Lawrence Livermore
National Laboratory under contract DE-AC52-07NA27344 and supported by LLNL
LDRD funding under project number 22-SI-003.

##  VII Author Declarations

###  VII.1 Conflict of Interest

The authors have no conflicts to disclose.

###  VII.2 Author Contributions

JXS contributed to the software development and the manuscript. JBV
contributed to the manuscript. LV supervised the project and reviewed the
manuscript.

###  VII.3 Data Availability

The data that supports the findings of this study are available within the
article and its supplementary material. The code used to generate the data is
released under the modified BSD license and is available at
https://github.com/materialsproject/pymatgen-analysis-defects.

## References

  * McCluskey [2012] H. McCluskey,  _Dopants and Defects in Semiconductors_ (Taylor & Francis, Andover, England, UK, 2012).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, “Self-consistent equations including exchange and correlation effects,” Phys. Rev. 140, A1133–A1138 (1965).
  * Freysoldt  _et al._ [2014] C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer, G. Kresse, A. Janotti, and C. G. Van de Walle, “First-principles calculations for point defects in solids,” Rev. Mod. Phys. 86, 253–305 (2014).
  * Dreyer  _et al._ [2020] C. E. Dreyer, A. Alkauskas, J. L. Lyons, and C. G. Van de Walle, “Radiative capture rates at deep defects from electronic structure calculations,” Phys. Rev. B 102, 085305 (2020).
  * Wickramaratne  _et al._ [2016] D. Wickramaratne, J.-X. Shen, C. E. Dreyer, M. Engel, M. Marsman, G. Kresse, S. Marcinkevičius, A. Alkauskas, and C. G. Van de Walle, “Iron as a source of efficient Shockley-Read-Hall recombination in GaN,” Appl. Phys. Lett. 109 (2016), 10.1063/1.4964831.
  * Alkauskas, Yan, and Van de Walle [2014] A. Alkauskas, Q. Yan, and C. G. Van de Walle, “First-principles theory of nonradiative carrier capture via multiphonon emission,” Phys. Rev. B 90, 075202 (2014).
  * Broberg  _et al._ [2018] D. Broberg, B. Medasani, N. E. R. Zimmermann, G. Yu, A. Canning, M. Haranczyk, M. Asta, and G. Hautier, “Pycdt: A python toolkit for modeling point defects in semiconductors and insulators,” Comput. Phys. Commun. 226, 165–179 (2018).
  * Kumagai  _et al._ [2021] Y. Kumagai, N. Tsunoda, A. Takahashi, and F. Oba, “Insights into oxygen vacancies from high-throughput first-principles calculations,” Phys. Rev. Mater. 5, 123803 (2021).
  * Mosquera-Lois  _et al._ [2023] I. Mosquera-Lois, S. R. Kavanagh, A. Walsh, and D. O. Scanlon, “Identifying the ground state structures of point defects in solids,” npj Comput. Mater. 9, 1–11 (2023).
  * Huang  _et al._ [2022] M. Huang, Z. Zheng, Z. Dai, X. Guo, S. Wang, L. Jiang, J. Wei, and S. Chen, “Dasp: Defect and dopant ab-initio simulation package,” J. Semicond. 43, 42101 (2022).
  * Yan  _et al._ [2024] Q. Yan, S. Kar, S. Chowdhury, and A. Bansil, “The Case for a Defect Genome Initiative,” Adv. Mater. n/a, 2303098 (2024).
  * Shen and Varley [2024] J.-X. Shen and J. Varley, “pymatgen-analysis-defects: A Python package for analyzing point defects in crystalline materials,” Journal of Open Source Software 9, 5941 (2024).
  * Ong  _et al._ [2013a] S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson, and G. Ceder, “Python materials genomics (pymatgen): A robust, open-source python library for materials analysis,” Computational Materials Science 68, 314–319 (2013a).
  * Rosen  _et al._ [2024] A. S. Rosen, M. Gallant, J. George, J. Riebesell, H. Sahasrabuddhe, J.-X. Shen, M. Wen, M. L. Evans, G. Petretto, D. Waroquiers, G.-M. Rignanese, K. A. Persson, A. Jain, and A. M. Ganose, “Jobflow: Computational Workflows Made Simple,” Journal of Open Source Software 9, 5995 (2024).
  * Xiong  _et al._ [2023] Y. Xiong, C. Bourgois, N. Sheremetyeva, W. Chen, D. Dahliah, H. Song, S. M. Griffin, A. Sipahigil, and G. Hautier, “High-throughput identification of spin-photon interfaces in silicon,” arXiv (2023), 10.48550/arXiv.2303.01594, 2303.01594 .
  * Freysoldt, Neugebauer, and de Walle [2009] C. Freysoldt, J. Neugebauer, and C. V. de Walle, “Fully ab initio finite-size corrections for charged-defect supercell calculations,” Phys. Rev. Lett. 102, 016402 (2009).
  * Ong  _et al._ [2013b] S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson, and G. Ceder, “Python Materials Genomics (pymatgen): A robust, open-source python library for materials analysis,” Comput. Mater. Sci. 68, 314–319 (2013b).
  * Shen, Horton, and Persson [2020] J.-X. Shen, M. Horton, and K. A. Persson, “A charge-density-based general cation insertion algorithm for generating new Li-ion cathode materials,” npj Comput. Mater. 6 (2020), 10.1038/s41524-020-00422-3.
  * Togo and Tanaka [2018] A. Togo and I. Tanaka, “Spglib: a software library for crystal symmetry search,” arXiv (2018), 1808.01590 .
  * Shen  _et al._ [2022] J.-X. Shen, J. M. Munro, M. K. Horton, P. Huck, S. Dwaraknath, and K. A. Persson, “A representation-independent electronic charge density database for crystalline materials,” Sci. Data 9, 1–7 (2022).
  * Wang  _et al._ [2021] A. Wang, R. Kingsbury, M. McDermott, M. Horton, A. Jain, S. P. Ong, S. Dwaraknath, and K. A. Persson, “A framework for quantifying uncertainty in DFT energy corrections,” Sci. Rep. 11, 1–10 (2021).
  * Stevanović  _et al._ [2012] V. Stevanović, S. Lany, X. Zhang, and A. Zunger, “Correcting density functional theory for accurate predictions of compound enthalpies of formation: Fitted elemental-phase reference energies,” Phys. Rev. B 85, 115104 (2012).
  * Bartók, Kondor, and Csányi [2013] A. P. Bartók, R. Kondor, and G. Csányi, “On representing chemical environments,” Phys. Rev. B 87, 184115 (2013).
  * De  _et al._ [2016] S. De, A. P. Bartók, G. Csányi, and M. Ceriotti, “Comparing molecules and solids across structural and alchemical space,” Phys. Chem. Chem. Phys. 18, 13754–13769 (2016).
  * Mathew  _et al._ [2017] K. Mathew, J. H. Montoya, A. Faghaninia, S. Dwarakanath, M. Aykol, H. Tang, I. heng Chu, T. Smidt, B. Bocklund, M. Horton, J. Dagdelen, B. Wood, Z.-K. Liu, J. Neaton, S. P. Ong, K. Persson, and A. Jain, “Atomate: A high-level interface to generate, execute, and analyze computational materials science workflows,” Comput. Mater. Sci. 139, 140–152 (2017).
  * Huber  _et al._ [2020] S. P. Huber, S. Zoupanos, M. Uhrin, L. Talirz, L. Kahle, R. Häuselmann, D. Gresch, T. Müller, A. V. Yakutovich, C. W. Andersen, F. F. Ramirez, C. S. Adorf, F. Gargiulo, S. Kumbhar, E. Passaro, C. Johnston, A. Merkys, A. Cepellotti, N. Mounet, N. Marzari, B. Kozinsky, and G. Pizzi, “AiiDA 1.0, a scalable computational infrastructure for automated reproducible workflows and data provenance,” Sci. Data 7, 1–18 (2020).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, “Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set,” Phys. Rev. B 54, 11169–11186 (1996).
  * Kühne  _et al._ [2020] T. D. Kühne, M. Iannuzzi, M. D. Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöß, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack, and J. Hutter, “Cp2k: An electronic structure and molecular dynamics software package - quickstep: Efficient and accurate electronic structure calculations,” The Journal of Chemical Physics 152, 194103 (2020).
  * Bergerhoff  _et al._ [1983] G. Bergerhoff, R. Hundt, R. Sievers, and I. D. Brown, “The inorganic crystal structure data base,” J. Chem. Inf. Comput. Sci. 23, 66–69 (1983).
  * Krukau  _et al._ [2006] A. V. Krukau, O. A. Vydrov, A. F. Izmaylov, and G. E. Scuseria, “Influence of the exchange screening parameter on the performance of screened hybrid functionals,” J. Chem. Phys. 125 (2006), 10.1063/1.2404663.
  * Lyons and Van de Walle [2017] J. L. Lyons and C. G. Van de Walle, “Computationally predicted energies and properties of defects in GaN,” npj Comput. Mater. 3, 1–10 (2017).
  * Frodason  _et al._ [2023] Y. K. Frodason, J. B. Varley, K. M. H. Johansen, L. Vines, and C. G. Van de Walle, “Migration of Ga vacancies and interstitials in $\beta\text{$-$}{\mathrm{Ga}}_{2}{\mathrm{O}}_{3}$,” Phys. Rev. B 107, 024109 (2023).
  * Shen  _et al._ [2023] J.-X. Shen, H. H. Li, A. Rutt, M. K. Horton, and K. A. Persson, “Topological graph-based analysis of solid-state ion migration,” npj Comput. Mater. 9, 1–5 (2023).
  * Ingebrigtsen  _et al._ [2019] M. E. Ingebrigtsen, A. Yu. Kuznetsov, B. G. Svensson, G. Alfieri, A. Mihaila, U. Badstübner, A. Perron, L. Vines, and J. B. Varley, “Impact of proton irradiation on conductivity and deep level defects in $\beta$-Ga2O3,” APL Mater. 7 (2019), 10.1063/1.5054826.
  * Varley  _et al._ [2011] J. B. Varley, H. Peelaers, A. Janotti, and C. G. Van de Walle, “Hydrogenated cation vacancies in semiconducting oxides,” J. Phys.: Condens. Matter 23, 334212 (2011).
  * Kyrtsos, Matsubara, and Bellotti [2017] A. Kyrtsos, M. Matsubara, and E. Bellotti, “Migration mechanisms and diffusion barriers of vacancies in ${\mathrm{Ga}}_{2}{\mathrm{O}}_{3}$,” Phys. Rev. B 95, 245202 (2017).
  * Liu  _et al._ [2014] B. Liu, V. R. Cooper, H. Xu, H. Xiao, Y. Zhang, and W. J. Weber, “Composition dependent intrinsic defect structures in SrTiO3,” Phys. Chem. Chem. Phys. 16, 15590–15596 (2014).
  * Baker  _et al._ [2017] J. N. Baker, P. C. Bowes, D. M. Long, A. Moballegh, J. S. Harris, E. C. Dickey, and D. L. Irving, “Defect mechanisms of coloration in Fe-doped SrTiO3 from first principles,” Appl. Phys. Lett. 110 (2017), 10.1063/1.4978861.
  * Janotti  _et al._ [2014] A. Janotti, J. B. Varley, M. Choi, and C. G. Van de Walle, “Vacancies and small polarons in ${\mathrm{SrTiO}}_{3}$,” Phys. Rev. B 90, 085202 (2014).
  * Varley, Janotti, and Van de Walle [2014] J. B. Varley, A. Janotti, and C. G. Van de Walle, “Hydrogenated vacancies and hidden hydrogen in SrTiO3,” Phys. Rev. B 89, 075202 (2014).
  * Choi, Oba, and Tanaka [2009] M. Choi, F. Oba, and I. Tanaka, “Role of Ti Antisitelike Defects in ${\mathrm{SrTiO}}_{3}$,” Phys. Rev. Lett. 103, 185502 (2009).
  * Alkauskas, Broqvist, and Pasquarello [2011] A. Alkauskas, P. Broqvist, and A. Pasquarello, “Defect levels through hybrid density functionals: Insights and applications,” Phys. Status Solidi B 248, 775–789 (2011).


# Visual Discourse Parsing

Arjun R Akula
University of California, Los Angeles
<EMAIL_ADDRESS>
&Song-Chun Zhu
University of California, Los Angeles
<EMAIL_ADDRESS>

###### Abstract

Text-level discourse parsing aims to unmask how two segments (or sentences) in
the text are related to each other. We propose the task of Visual Discourse
Parsing, which requires understanding discourse relations among scenes in a
video. Here we use the term scene to refer to a subset of video frames that
can better summarize the video. In order to collect a dataset for learning
discourse cues from videos, one needs to manually identify the scenes from a
large pool of video frames and then annotate the discourse relations between
them. This is clearly a time consuming, expensive and tedious task. In this
work, we propose an approach to identify discourse cues from the videos
without the need to explicitly identify and annotate the scenes. We also
present a novel dataset containing 310 videos and the corresponding discourse
cues to evaluate our approach. We believe that many of the multi-discipline
Artificial Intelligence problems such as Visual Dialog and Visual Storytelling
would greatly benefit from the use of visual discourse cues.

##  1 Introduction

Discourse structure aids in understanding a piece of text by linking it with
other text units (such as surrounding clauses, sentences, etc.) from its
context Carlson et al. (2003); Soricut and Marcu (2003); LeThanh et al.
(2004). A text span may be linked to another span through semantic
relationships such as contrast relation, causal relation, etc. Marcu and
Echihabi (2002); Duverle and Prendinger (2009). Text-level discourse parsing
algorithms aim to unmask such relationships in text, which is central to many
downstream natural language processing (NLP) applications such as information
retrieval, text summarization, sentiment analysis Wang and Lan (2015) and
question answering Chai and Jin (2004); Akula et al. (2013); Akula (2015).

Figure 1: (a) Video with 9 frames; (b) Textual description of video; (c) RST
discourse structure (represented as sequence of words) of description in (b).
The notations LEFT and RIGHT represent the direction (nuclearity) of the
rhetorical relations; (d) RST discourse structure of the video using 3 frames,
i.e. scenes. The symbols N and S indicate the nucleus and satellite of each
rhetorical relation.

Recently, there has been a lot of focus on multi-discipline Artificial
Intelligence (AI) research problems such as visual storytelling Huang et al.
(2016) and visual dialog Das et al. (2016). Solving these problems requires
multi-modal knowledge that combines computer vision (CV), NLP, and knowledge
representation & reasoning (KR), making the need for commonsense knowledge and
complex reasoning more essential. In an effort to fill this need, we introduce
the task of Visual Discourse Parsing.

Task Definition. The concrete task in Visual Discourse Parsing is the
following - given a video, understand discourse relationships among its
scenes. Specifically, given a video, the task is to identify a scene’s
relation with the context. Here we use the term scene to refer to a subset of
video frames that can better summarize the video. We use Rhetorical Structure
Theory (RST) Mann and Thompson (1988) to capture discourse relations among the
scenes.

Consider for example, nine frames of a video shown in Figure 1. We can
represent the discourse structure of this video using only 3 out of 9 frames,
i.e. there are only 3 scenes for this video. The discourse structure in the
Figure 1 interprets the video as follows: the event “person going to the
bathroom and cleaning his stains” is caused by the event “the person spilling
coffee over his shirt”; the event “the person used his handkerchief to dry the
water on his shirt” is simply an elaboration of the event “person going to the
bathroom and cleaning his stains”.

It is a time consuming, expensive and tedious task to collect a dataset for
learning discourse cues from videos. This is because one needs to manually
identify the scenes from a large pool of video frames and annotate the
discourse relations between them. To this end, we propose an approach to
identify discourse cues from the videos without the need to explicitly
identify and annotate the scenes. We also present a novel dataset containing
310 videos and the corresponding discourse cues to evaluate our approach.

##  2 Approach

Algorithm 1 presents our approach for learning a model to identify discourse
structure from the videos. Firstly, we generate natural language text
descriptions from the videos automatically using video captioning methods such
as Yu et al. (2016), Venugopalan et al. (2016) and Pasunuru and Bansal (2017).
Secondly, we obtain discourse structures of the above text descriptions using
text-level discourse parsers such as Duverle and Prendinger (2009) and Ji and
Eisenstein (2014). We represent the discourse structure as a sequence of words
(see Figure 1).

Next we use an end-to-end trainable architecture for learning to predict the
text-level discourse structures (i.e. sequence of words) from the videos (i.e.
sequence of video frames). This gives us a model to map videos to their
corresponding text-level discourse structures. Finally we use saliency methods
Ramanishka et al. (2017) to replace the textual descriptions in the discourse
structure, i.e. elementary discourse units (EDUs), with the scenes. For
example, in Figure 1, (b) is the text description of video shown in (a). Text-
level discourse structure of (b) is shown in (c) as a sequence of words. We
then map (c) to (d) using saliency methods.

Algorithm 1 Our approach for Visual Discourse Parsing.

1:Generate natural language descriptions of videos.

2:Find discourse structure of text descriptions.

3:Train a sequence to sequence model to learn the mapping from video frames to
the sequence of words in the discourse structure.

4:Use saliency methods to map EDUs in text-level discourse structure to
scenes.

The quality of natural language descriptions and the their text-level
discourse structures is crucial for learning a robust model. The state-of-the-
art video captioning and text-level discourse parsing approaches, as we found
in our experiments, may generate a lot of noise in their outputs. While
developing our corpus, we manually performed these two steps. These manual
annotations are still much easier to perform compared to the tedious task of
directly annotating discourse relations between video frames.

In the step 3 of our algorithm, we use the standard machine translation
encoder-decoder RNN model Sutskever et al. (2014). As RNN suffers from
decaying of gradient and blowing-up of gradient problem, we use LSTM units,
which are good at memorizing long-range dependencies due to forget-style gates
Hochreiter and Schmidhuber (1997). The sequence of video frames are passed to
the encoder. The last hidden state of the encoder is then passed to the
decoder. The decoder generates the discourse structure as a sequence of words.
Let the input sequence of video frames be
$\mathbf{x}=\left(x\\_1,...,x\\_p\right)$ and the output sequence of words as
$\mathbf{y}=\left(y\\_1,...,y\\_n\right)$. The distribution of the output
sequence w.r.t. the input sequence is:

$p(y_{1},...,y_{n}\mid x_{1},...,x_{p})=\prod_{t=1}^{n}p(y_{t}\mid h_{t}^{d})$
(1)

where $h_{t}^{d}$ is the hidden state at the $t^{th}$ time step of the
decoding LSTM.

Soft Attention: We further improve our encoder-decoder model using an
attention based sequence-to-sequence model Bahdanau et al. (2014). The
attention weights act as an alignment mechanism by re-weighting the encoder
hidden states that are more relevant for decoder time step.

RNN Type | #Hidden Units | Bidirectional | #Layers | Relations | Edges | Relations+Edges | Bleu4
---|---|---|---|---|---|---|---
LSTM | 256 | NO | 1 | 0.3 | 0.51 | 0.21 | 3.8
LSTM | 512 | NO | 1 | 0.52 | 0.62 | 0.42 | 6.2
LSTM | 1024 | YES | 1 | 0.49 | 0.51 | 0.42 | 4.3
LSTM | 1024 | NO | 1 | 0.35 | 0.51 | 0.21 | 4
LSTM | 512 | NO | 2 | 0.35 | 0.51 | 0.21 | 5.4
LSTM | 512 | NO | 3 | 0.56 | 0.62 | 0.42 | 6.8
LSTM | 512 | NO | 4 | 0.56 | 0.62 | 0.42 | 6.8
GRU | 512 | NO | 1 | 0.3 | 0.51 | 0.21 | 4.2
Table 1: Evaluation using sequence-to-sequence model without Attention. RNN Type | #Hidden Units | Bidirectional | #Layers | #Attention Type | Relations | Edges | Relations+Edges | Bleu4
---|---|---|---|---|---|---|---|---
LSTM | 512 | NO | 1 | general | 0.63 | 0.69 | 0.53 | 10.6
LSTM | 512 | NO | 1 | dot | 0.52 | 0.65 | 0.45 | 8.9
LSTM | 512 | NO | 1 | concat | 0.52 | 0.65 | 0.45 | 9.9
LSTM | 512 | NO | 2 | general | 0.52 | 0.65 | 0.45 | 8.8
LSTM | 512 | NO | 3 | general | 0.5 | 0.65 | 0.39 | 7.8
Table 2: Evaluation using sequence-to-sequence model using Attention.

##  3 Experiments

We developed a new dataset containing 310 videos. These videos are shot at
various settings such as playing sports (Table Tennis, Frisbee, Tennis,
Rugby), bus stop, dining hall, elevator, classroom, library, garden, study
room, etc. On average, the length of each video is about 19 seconds. We first
manually generated descriptions of each video and then annotated the discourse
structure of these descriptions - with the help of 5 graduate students. Each
video is annotated by at least 2 students. We solved the disagreements found
in the annotations together. As the training data is not large enough, we
chose short videos and described each video using only three sentences, i.e.
discourse structure (RST tree) of each video contains only two relations and
two edges. This reduces the total number of parameters that need to be learned
from the end-to-end training (in step 3 of Algorithm 1).

We evaluate our approach by using the following three metrics:

  1. (a)

BLEU score We used the BLEU score Papineni et al. (2002) to evaluate the
translation quality of the discourse structure generated from the videos. We
computed our BLEU score on the tokenized predictions and ground truth.

  2. (b)

Relations Accuracy Each video, in our dataset, contains two discourse
relations. The Relations Accuracy metric is defined as the total number of
relations correctly predicted by the model.

  3. (c)

Edges Accuracy Each video, in our dataset, contains two edges. The Edges
Accuracy metric is defined as the total number of edges (i.e. RST node
nuclearity directions) correctly predicted by the model.

  4. (d)

Relations+Edges Accuracy Here, we compute the correctness of the complete
discourse structure, i.e. the predicted discourse structure will be considered
correct only if all the relations and the edges are correctly predicted by the
model.

We use 210 videos for training, 30 videos for validation and the rest 70
videos for testing. We fix our sampling rate to 5fps to bring uniformity in
the temporal representation of actions across all videos. These sampled frames
are then converted into features using VGGNet Simonyan and Zisserman (2014).
We initialize the decoder embedding with Google pre-trained word2vec word
embeddings Mikolov et al. (2013). We tune all hyperparameters using our
validation data: learning rate, weight initializations, hidden states. We use
a 1024-dimension RNN hidden state size. We use Adam optimizer Kingma and Ba
(2014) and apply a dropout of 0.5.

Table 1 presents our results, reporting several configurations of the encoder-
decoder RNN model, using the four evaluation metrics. The metrics Relations,
Edges and Relations+Edges accuracies are scaled between 0 and 1. It may be
noted that RNN configuration using LSTM unit with 512 hidden units and 3,4
encoder layers outperformed other RNN configurations. Table 2 reports the
performance of our attention based sequence to sequence model. Attention based
models gave better accuracies than the simple sequence to sequence models. In
particular, the LSTM unit with 512 hidden units and single encoding layer
outperformed other configurations. These results are encouraging considering
the small size of our training dataset.

##  4 Related Work

Text-level Discourse Parsing
Several works Carlson et al. (2003); LeThanh et al. (2004); Marcu and Echihabi
(2002); Duverle and Prendinger (2009) have been proposed in the past to parse
text documents into RST style discourse representation. Carlson et al. (2003)
proposed a data driven approach using RST Discourse Treebank. Reitter (2003)
proposed a chart-parsing-style techniques to derive discourse trees from
documents. LeThanh et al. (2004) proposed a three-stage approach to find
discourse parse. It segments the elementary discourse units (EDUs) into trees
for each successive level of document: first at the sentence level, then
paragraph and finally at the document level. However, using a three-stage
approach exponentially increased their search space, making it computationally
intractable to find the optimal discourse tree. Duverle and Prendinger (2009)
proposed end-to-end parsing algorithms. However, in none of these works, the
problem of extracting discourse information from videos has been addressed.
Video Captioning
There are few works on video captioning which come closer to our line of
research. Video captioning is the task of describing the content of a video.
Guadarrama et al. (2013); Thomason et al. (2014); Yu et al. (2016) proposed a
multi-stage approach to identify dependency roles such as subject and object
in order to generate a description of the video. More recently, Venugopalan et
al. (2016) proposed a sequence-to-sequence model with encoder and decoder
RNNs. Pasunuru and Bansal (2017) proposed a multi-task learning approach to
generate captions from the video. However, all of these works mainly focus on
generating text from the video but not on understanding relationships between
the scenes (or frames).

##  5 Conclusions

We introduced a new AI task - Visual Discourse Parsing, where the AI agent
needs to understand discourse relations among scenes in a video. We presented
an end-to-end learning approach to identify the discourse structure of videos.
Central to our approach is the use of text descriptions of videos to identify
discourse relations. In the future, we plan to extend this dataset to include
longer videos that need more than three sentences to describe. We also intend
to experiment with multi-task learning approaches. Our results indicate that
there is significant scope for improvement.

## References

  * Akula et al. (2013) Arjun Akula, Rajeev Sangal, and Radhika Mamidi. 2013.  A novel approach towards incorporating context processing capabilities in nlidb system.  In Proceedings of the sixth international joint conference on natural language processing. pages 1216–1222.
  * Akula (2015) Arjun R Akula. 2015.  A novel approach towards building a generic, portable and contextual nlidb system.  International Institute of Information Technology Hyderabad .
  * Bahdanau et al. (2014) Dzmitry Bahdanau, Kyunghyun Cho, and Yoshua Bengio. 2014.  Neural machine translation by jointly learning to align and translate.  arXiv preprint arXiv:1409.0473 .
  * Carlson et al. (2003) Lynn Carlson, Daniel Marcu, and Mary Ellen Okurowski. 2003.  Building a discourse-tagged corpus in the framework of rhetorical structure theory.  In Current and new directions in discourse and dialogue, Springer, pages 85–112.
  * Chai and Jin (2004) Joyce Y Chai and Rong Jin. 2004.  Discourse structure for context question answering.  In Proceedings of the Workshop on Pragmatics of Question Answering at HLT-NAACL 2004.
  * Das et al. (2016) Abhishek Das, Satwik Kottur, Khushi Gupta, Avi Singh, Deshraj Yadav, José MF Moura, Devi Parikh, and Dhruv Batra. 2016.  Visual dialog.  arXiv preprint arXiv:1611.08669 .
  * Duverle and Prendinger (2009) David A Duverle and Helmut Prendinger. 2009.  A novel discourse parser based on support vector machine classification.  In Proceedings of the Joint Conference of the 47th Annual Meeting of the ACL and the 4th International Joint Conference on Natural Language Processing of the AFNLP: Volume 2-Volume 2. Association for Computational Linguistics, pages 665–673.
  * Guadarrama et al. (2013) Sergio Guadarrama, Niveda Krishnamoorthy, Girish Malkarnenkar, Subhashini Venugopalan, Raymond Mooney, Trevor Darrell, and Kate Saenko. 2013.  Youtube2text: Recognizing and describing arbitrary activities using semantic hierarchies and zero-shot recognition.  In Proceedings of the IEEE international conference on computer vision. pages 2712–2719.
  * Hochreiter and Schmidhuber (1997) Sepp Hochreiter and Jürgen Schmidhuber. 1997.  Long short-term memory.  Neural computation 9(8):1735–1780.
  * Huang et al. (2016) Ting-Hao Kenneth Huang, Francis Ferraro, Nasrin Mostafazadeh, Ishan Misra, Aishwarya Agrawal, Jacob Devlin, Ross Girshick, Xiaodong He, Pushmeet Kohli, Dhruv Batra, et al. 2016.  Visual storytelling.  In Proceedings of the 2016 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies. pages 1233–1239.
  * Ji and Eisenstein (2014) Yangfeng Ji and Jacob Eisenstein. 2014.  Representation learning for text-level discourse parsing.  In Proceedings of the 52nd Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers). volume 1, pages 13–24.
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba. 2014.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980 .
  * LeThanh et al. (2004) Huong LeThanh, Geetha Abeysinghe, and Christian Huyck. 2004.  Generating discourse structures for written texts.  In Proceedings of the 20th international conference on Computational Linguistics. Association for Computational Linguistics, page 329\.
  * Mann and Thompson (1988) William C Mann and Sandra A Thompson. 1988.  Rhetorical structure theory: Toward a functional theory of text organization.  Text-Interdisciplinary Journal for the Study of Discourse 8(3):243–281.
  * Marcu and Echihabi (2002) Daniel Marcu and Abdessamad Echihabi. 2002.  An unsupervised approach to recognizing discourse relations.  In Proceedings of the 40th Annual Meeting on Association for Computational Linguistics. Association for Computational Linguistics, pages 368–375.
  * Mikolov et al. (2013) Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg S Corrado, and Jeff Dean. 2013.  Distributed representations of words and phrases and their compositionality.  In Advances in neural information processing systems. pages 3111–3119.
  * Papineni et al. (2002) Kishore Papineni, Salim Roukos, Todd Ward, and Wei-Jing Zhu. 2002.  Bleu: a method for automatic evaluation of machine translation.  In Proceedings of the 40th annual meeting on association for computational linguistics. Association for Computational Linguistics, pages 311–318.
  * Pasunuru and Bansal (2017) Ramakanth Pasunuru and Mohit Bansal. 2017.  Multi-task video captioning with video and entailment generation.  arXiv preprint arXiv:1704.07489 .
  * Ramanishka et al. (2017) Vasili Ramanishka, Abir Das, Jianming Zhang, and Kate Saenko. 2017.  Top-down visual saliency guided by captions.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR). volume 1, page 7.
  * Reitter (2003) David Reitter. 2003.  Simple signals for complex rhetorics: On rhetorical analysis with rich-feature support vector models.  In LDV Forum. volume 18, pages 38–52.
  * Simonyan and Zisserman (2014) Karen Simonyan and Andrew Zisserman. 2014.  Very deep convolutional networks for large-scale image recognition.  arXiv preprint arXiv:1409.1556 .
  * Soricut and Marcu (2003) Radu Soricut and Daniel Marcu. 2003.  Sentence level discourse parsing using syntactic and lexical information.  In Proceedings of the 2003 Conference of the North American Chapter of the Association for Computational Linguistics on Human Language Technology-Volume 1. Association for Computational Linguistics, pages 149–156.
  * Sutskever et al. (2014) Ilya Sutskever, Oriol Vinyals, and Quoc V Le. 2014.  Sequence to sequence learning with neural networks.  In Advances in neural information processing systems. pages 3104–3112.
  * Thomason et al. (2014) Jesse Thomason, Subhashini Venugopalan, Sergio Guadarrama, Kate Saenko, and Raymond J Mooney. 2014.  Integrating language and vision to generate natural language descriptions of videos in the wild.  In Coling. volume 2, page 9.
  * Venugopalan et al. (2016) Subhashini Venugopalan, Lisa Anne Hendricks, Raymond Mooney, and Kate Saenko. 2016\.  Improving lstm-based video description with linguistic knowledge mined from text.  arXiv preprint arXiv:1604.01729 .
  * Wang and Lan (2015) Jianxiang Wang and Man Lan. 2015.  A refined end-to-end discourse parser.  In Proceedings of the Nineteenth Conference on Computational Natural Language Learning-Shared Task. pages 17–24.
  * Yu et al. (2016) Haonan Yu, Jiang Wang, Zhiheng Huang, Yi Yang, and Wei Xu. 2016.  Video paragraph captioning using hierarchical recurrent neural networks.  In Proceedings of the IEEE conference on computer vision and pattern recognition. pages 4584–4593.