Grand canonically optimized grain boundary phases in hexagonal close-packed
titanium

Enze Chen1,2,3,4,∗, Tae Wook Heo2, Brandon C. Wood2, Mark Asta1,3, Timofey
Frolov2,∗

1 Department of Materials Science and Engineering, University of California,
Berkeley, CA 94720, USA
2 Materials Science Division, Lawrence Livermore National Laboratory,
Livermore, CA 94550, USA
3 Materials Sciences Division, Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, USA
4 Present address: Department of Materials Science and Engineering, Stanford
University, Stanford, CA 94305, USA
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Grain boundaries (GBs) profoundly influence the properties and performance of
materials, emphasizing the importance of understanding the GB structure and
phase behavior. As recent computational studies have demonstrated the
existence of multiple GB phases associated with varying the atomic density at
the interface, we introduce a validated, open-source GRand canonical Interface
Predictor (GRIP) tool that automates high-throughput, grand canonical
optimization of GB structures. While previous studies of GB phases have almost
exclusively focused on cubic systems, we demonstrate the utility of GRIP in an
application to hexagonal close-packed titanium. We perform a systematic high-
throughput exploration of tilt GBs in titanium and discover previously
unreported structures and phase transitions. In low-angle boundaries, we
demonstrate a coupling between point defect absorption and the change in the
GB dislocation network topology due to GB phase transformations, which has
important implications for the accommodation of radiation-induced defects.

## Introduction

Grain boundaries (GBs) are interfacial defects in crystalline materials that
have long been studied for their influence on materials properties and
performance [1]. Given their ability to exist in multiple stable and
metastable states, which have been termed GB phases [2] or complexions [3], it
is desirable to obtain an atomic-level understanding of the GB structures and
possible phase transition pathways between them [4, 5]. The structure–property
relationships of these interfacial phases are believed to have a profound
influence on an array of phenomena, such as diffusion [6] and GB migration [7]
in materials.

Recent experiments have provided direct [8] and indirect [9] evidence for GB
phase stability, coexistence, and transitions in metallic systems; however,
given the vast five-dimensional space characterizing the macroscopic degrees
of freedom (DOF) for GBs, it is not yet clear where these phases may appear.
Atomistic simulations provide a powerful tool to guide such searches and
unveil the microscopic mechanisms underlying the formation of GB phases [10].
Previous atomistic modeling studies have discovered a diverse array of GB
phases present in face-centered cubic (FCC) [4, 11, 12], body-centered cubic
(BCC) [13, 14], diamond cubic [15, 16], and other cubic systems [17, 18]. One
notable feature shared by the aforementioned works is the ability to add or
remove atoms from the GB region in the simulation cell, i.e., grand canonical
optimization (GCO), which was required to access new ground states and
metastable states. While the exchange of atoms at an interface could naturally
occur in real polycrystalline materials due to diffusion, irradiation, and
mechanical deformation at finite temperature, this variation is omitted in the
majority of computational simulations employing the $\gamma$-surface method
[19]. The $\gamma$-surface method is the traditional technique for simulating
GBs where only relative translations are allowed between two bulk slabs before
the atoms are relaxed using conjugate gradient energy minimization to their
equilibrium positions at $0\text{\,}\mathrm{K}$. It is often adopted for its
simplicity, but the deficiencies exposed by the previous studies suggest that
more DOF must be considered during optimization in order to find the true
ground-state structure in certain GBs. A few alternative approaches from the
literature for atomistic modeling of GBs include high-temperature molecular
dynamics (MD) simulations [4, 15, 20], Monte Carlo sampling [21, 22], and
evolutionary algorithms [17, 23, 22, 11, 24] that can access a greater
diversity of structures and atomic densities at the interface, with a
concomitant trade-off in computational complexity. These algorithms have
employed GCO for GB structures in a variety of systems, although seldom in a
high-throughput manner [21, 20], and it remains unclear if the ubiquity of GB
phases that they have yielded extends to lower-symmetry crystalline systems
that are ubiquitous in nature and engineering applications.

A particular system of immense technological relevance is the hexagonal close-
packed (HCP) crystal structure, which is considerably more complex than cubic
systems, as it displays anisotropy in its crystalline lattice vectors and a
basis containing more than one atom. This structure is adopted by elemental
metals such as Mg, Zr, and Ti, the last of which ($\alpha$-Ti) will be the
focus of this work. Ti alloys are important structural alloys for aerospace,
biomedical, and energy applications, particularly where high specific strength
and strong corrosion resistance are desired [25]. The importance of GBs in the
$\alpha$-Ti system is highlighted in recent studies that used grain refinement
to mitigate low-temperature oxygen embrittlement in $\alpha$-Ti [26] and 3D
electron backscatter diffraction to map out the complete GB character
distribution in this system [27]. In a previous study [28], we used an
evolutionary algorithm to discover a ground-state structure for the
$\\{11\bar{2}4\\}$$[1\bar{1}00]$ twin boundary (TB) in $\alpha$-Ti that was in
closer agreement to density functional theory (DFT) calculations and high-
resolution transmission electron microscopy results than previously reported
structures. In addition to these experimental works, there are also several
atomistic simulation studies in the literature that systematically model
symmetric tilt grain boundaries (STGBs) in $\alpha$-Ti along the $[0001]$ [29,
30], $[1\bar{1}00]$ [31, 32, 33, 34], and $[1\bar{2}10]$ [35, 34] tilt axes.
Despite the simplicity of STGBs—only a tilt axis and tilt angle ($2\theta$)
are required to describe the crystallographic misorientation between two bulk
crystals—they include coherent TBs as an important subclass, several of which
are experimentally observed in deformation microstructures and thus important
for mechanical behavior in $\alpha$-Ti [25, 26, 28]. STGBs are also model
systems to study the geometric relationships of defects at the interface [35];
however, as the previous studies utilized the $\gamma$-surface method, it is
important to clarify the effects of GCO on STGB structure in $\alpha$-Ti and
more broadly whether interfacial phases exist in HCP metals.

Herein, we perform GCO of low-index STGBs in $\alpha$-Ti using an open-source
GRand canonical Interface Predictor (GRIP) tool that we developed to
rigorously sample microscopic DOF at the GB. We use this tool along with
empirical potentials to perform GB structure search, discovering new ground-
state structures and GB phases. We further employ high-temperature MD
simulations to explore the $\\{21\bar{3}0\\}$$[0001]$ STGB and demonstrate GB
phase (meta)stability and phase transitions through a novel dislocation-
pairing mechanism. We conclude by discussing the broader implications of these
results on GB phase behavior in HCP systems and how the GRIP tool can benefit
future studies for diverse crystal structures.

## Results

### Grand canonical optimization—the GRIP tool

GB structure prediction is a long-standing challenge in materials modeling
that requires rigorous and often advanced sampling of possible interfacial
structures. Previous studies of GBs in HCP metals generated the interfaces
using the common $\gamma$-surface method [19], which is not guaranteed to
yield the true ground state configuration in general [15, 11]. In the
traditional approach, conjugate gradient minimization from different starting
points representing distinct relative transitions of the grains across the
boundary simply allows the atoms to fall into a nearby local minimum, which
may be far away from the ground state. For example, complex GB core
configurations may exist that require significant rearrangement of the
constituent atoms [15, 11].

The other significant limitation of the $\gamma$-surface method is that it is
not grand canonical: All GBs created using this method are composed of the
same number of atoms derived from the constituent perfect half-crystals. This
poses a substantial constraint because many other structures, including true
ground states, can be realized out of a different number of atoms at the
interface [15, 4, 13]. For STGBs and a fixed reconstruction area, the number
of distinct atomic densities that can give rise to different GB structures is
given by the total number of atoms in one atomic plane parallel to the
boundary, which we denote $N^{\mathrm{bulk}}_{\mathrm{plane}}$. This quantity
is the limit because removing a full plane of atoms from a crystal will return
the exact same configuration up to a relative grain translation.

Figure 1: Success rate as a function of search parameters in GRIP. The success
rate of finding the ground state is sensitive to the search parameters,
including temperature and duration of the MD sampling. The optimal parameters
are not known a priori and vary for each particular boundary and
reconstruction, demonstrating why rigorous sampling of the parameter space is
critical. (a) The U shape of the $E_{\mathrm{gb}}$ vs. $T$ plot illustrates
the inefficient frozen dynamics of the GB structure at low temperatures and
the generation of disordered liquid-like GBs at very high temperatures for a
representative GB. (b) The fraction of ground-state structures out of all
structures sampled at each temperature are plotted and fitted with a Gaussian
distribution. (c) A sufficiently large number of MD steps is required to
obtain the ground state, even at an optimal $T$. (d) $\Delta E_{\mathrm{gb}}$
is plotted against the number of MD steps. The least-squares regression line
shows the convergence in energy at longer duration.

Here we address these shortcomings through the development of an open-source
tool GRIP to perform grand canonical GB structure search. During the
optimization, we systematically explore all possible microscopic DOF by
sampling different relative grain translations and atomic densities (see
Methods for details). The latter is accomplished by randomly removing a
fraction of atoms between $0$ and $N^{\mathrm{bulk}}_{\mathrm{plane}}$ from
the boundary plane. For a fixed translation and number of GB atoms, we
optimize the GB structure using dynamic sampling (performed here using MD
simulations) at different temperatures within a wide window between room
temperature and $1200\text{\,}\mathrm{K}$ (approximately
$T_{\alpha\rightarrow\beta}$ for Ti), and for different durations up to
$0.6\text{\,}\mathrm{ns}$. At the end of each MD run, we perform conjugate
gradient energy minimization at $0\text{\,}\mathrm{K}$ until convergence
before calculating the GB energy, $E_{\mathrm{gb}}$ (see Equation 1 in
Methods). The random and diverse GB structure initialization coupled with the
extensive dynamic sampling for each microscopic DOF done by hundreds of
parallel calculations ensure a rigorous GB structure exploration.

To further underscore the need for rigorous sampling, Figure 1 shows the
structural diversity and success rates from randomly sampling MD simulation
parameters, namely temperature ($T$) and duration (MD steps). Figure 1a shows
the range of $E_{\mathrm{gb}}$ as a function of $T$ when the duration is
fixed, and only at intermediate temperatures does the algorithm find the
ground-state structure. For this representative boundary, the U shape of the
$E_{\mathrm{gb}}$ vs. $T$ plot illustrates the inefficient frozen dynamics of
the GB structure at low temperatures and the generation of disordered liquid-
like GBs at very high temperatures. From this data, we can compute the
probability of finding the ground state—calculated as the fraction of ground-
state structures out of all sampled structures at each $T$—which peaks at
approximately $1150\text{\,}\mathrm{K}$ and is zero for very low and very high
temperatures (Figure 1b). These panels illustrate the existence of an optimal
$T$ that is sensitive to the structural DOF of each system and outside of
which the GB may fail to be optimized. Analogously, simply choosing an optimal
$T$ (e.g., $1000\text{\,}\mathrm{K}$) is insufficient, as too few MD steps
will never achieve the ground state, as shown in Figure 1c. The energy range,
$\Delta E_{\mathrm{gb}}$, is plotted in Figure 1d to show how the spread
generally decreases as the MD duration increases; however, we emphasize that
the optimal parameters are not known a priori. These optimal parameters can
vary significantly not only with the boundary character described by the five
macroscopic DOF, but also for larger area reconstructions of the same
boundary. The uniform sampling of possible optimization parameters and GB DOFs
implemented in GRIP allows for robust high-throughput optimization of large GB
datasets.

As motivated in the Introduction, we showcase the performance of GRIP in the
following sections through a detailed analysis of STGBs in HCP $\alpha$-Ti.
Importantly, however, we note that we also comprehensively benchmark our tool
by reproducing well-studied literature results for tilt and twist GBs in
elemental cubic metals [4, 14] and more challenging covalently-bonded, lower-
symmetry systems [15, 21] (Supplementary Figure S1). Even in the thoroughly
studied BCC W system [14], we discover a new ground-state structure with a
different GB atomic density and markedly different dislocation network in the
GB (Supplementary Figure S2). Such results, while not discussed further in
this work, underscore the opportunities of having a robust method for
exploring GB phase space across disparate chemical systems. The new ground
states also position GRIP as a tool capable of advancing the state of the art
in GB structure prediction through its extensive dynamic sampling of the
relevant DOF.

### Survey of new GB phases in HCP $\alpha$-Ti

Figure 2: Calculation of the number of atoms per plane,
$N^{\mathrm{bulk}}_{\mathrm{plane}}$, in HCP. The accurate calculation of
$N^{\mathrm{bulk}}_{\mathrm{plane}}$ ensures that GB structures with all
possible atomic densities are explored. Because HCP crystals have two basis
atoms, two different cases are possible when (a) all atoms inside the plane
have the same $z$-position, or (b) they have two distinct $z$-positions
resulting in two structurally different surface terminations. In both cases,
$N^{\mathrm{bulk}}_{\mathrm{plane}}$ is calculated as the total number of
atoms found inside the region spanned by the hexagonal interplanar spacing,
$d_{hkl}$ (boxed). The distinct $z$-positions of atoms belonging to the same
plane are indicated by the dashed magenta lines.

The two-atom basis of HCP Ti presents additional considerations during
optimization, and Figure 2 illustrates one nuance in having two possible cases
of calculating $N^{\mathrm{bulk}}_{\mathrm{plane}}$. For the orientation shown
in Figure 2a, all atoms found inside the planar region have the same
$z$-positions indicated by the dashed magenta line. Such orientations are
analogous to cubic systems and have only one distinct surface termination. For
the second case shown in Figure 2b, the atoms belonging to the same plane can
have two distinct $z$-positions, giving rise to two structurally different
surface terminations. These two distinct terminations are possible because HCP
has two basis atoms. In all orientations, the thickness of the planar region,
$d_{hkl}$, corresponds to the smallest normal component of a lattice vector
connecting two atoms on the same sublattice with different $z$-positions. We
further note that this definition of a plane of atoms works for both cases,
allowing us to uniformly apply it in calculating GB atomic density, $n$ (see
Equation 2 in Methods).

Figure 3: Representative GB structure searches using GRIP. (a) GB energy
$E_{\mathrm{gb}}$ vs. GB atomic density $n$ reveals two GB phases of $\Sigma
49$$\\{11\,2\,\overline{13}\,0\\}$$[0001]$ which are shown in panels (b) and
(c). (b) The $n=0$ GB phase, composed of edge dislocations, is metastable at 0
K. (c) The $n=0.75$ GB phase is the ground state, highlighting the importance
of GCO. (d) $E_{\mathrm{gb}}$ vs. $n$ for $\Sigma
13$$\\{31\bar{4}0\\}$$[0001]$. On the right, two orthogonal projections are
shown for each minimum-energy structure obtained using (e) GRIP at $n=0.5$,
(f) the $\gamma$-surface method, and (g) GRIP at $n=0$. The atoms are colored
according to the common neighbor analysis (CNA) in OVITO [36, 37].

Figure 3 shows the results of the GRIP searches for two representative
boundaries evaluated using a modified embedded-atom method (MEAM) potential
[38], illustrating the need for grand canonical structure optimization for GBs
in HCP Ti. Panels (a) and (e) show plots of $E_{\mathrm{gb}}$ vs. $n$, which
was introduced for cubic crystals in our previous work [4, 11]. Each point on
the plot corresponds to a particular GB structure obtained after energy
minimization. The thorough exploration enabled by GRIP generates hundreds of
distinct structures covering different densities and energies.

For the $\\{11\,2\,\overline{13}\,0\\}$$[0001]$ GB, the structure search
identifies two GB phases with different atomic densities $n=0$ and $n=0.75$.
The structures are illustrated in panels (b) and (c), respectively. The $n=0$
phase does not require insertion or removal of atoms and is metastable at
$0\text{\,}\mathrm{K}$. It is composed of well-separated cores of edge
dislocations with Burgers vectors $\mathbf{b}_{\text{I}}=\frac{1}{3}\langle
1\bar{2}10\rangle$, as identified in green by the dislocation extraction
algorithm (DXA) in OVITO [39, 37]. The newly predicted ground state of this
boundary has $n=0.75$ and thus cannot be generated by using the simplistic
$\gamma$-surface approach or sampling different terminations. Its structure is
significantly different from the $n=0$ state, where the dislocation cores
overlap and the boundary structure appears completely flat. The energy of the
ground state (corresponding to $n=0.75$) is $1\text{\,}\mathrm{\char
37\relax}$ lower than that of the metastable phase ($n=0$). We finally note
that this GB is a single-surface termination type of boundary where all atoms
belonging to a bulk plane have the same $z$-coordinate.

Figure 4: GB energy map of $[0001]$ STGBs. GRIP finds new ground states and
multiple GB phases within the entire misorientation range, demonstrating the
need for GCO. Each subplot is analogous to the boundary in Figure 3a, marking
the minimum-energy structures at different GB atomic densities. The squares
mark the ground state for each tilt angle, the larger circles mark a
metastable state, and the blue triangles mark the $\gamma$-surface structure.
The green and purple marker colors correspond to different GB phases
(commensurate with the dislocation core colors in Figure 6).

The GRIP search for the $\\{31\bar{4}0\\}$$[0001]$ GB shown in Figure 3d
illustrates that the microscopic descriptor $n$ properly captures all possible
distinct GB configurations, even for orientations with two distinct surface
terminations. Similar to the previous boundary, the prediction of the ground-
state structure at $n=0.5$ also requires an insertion (or removal) of half of
the atoms in one $\\{31\bar{4}0\\}$ plane; however, different from the first
example, this particular ground-state structure can also be generated using
the $\gamma$-surface approach that considers two possible surface
terminations. The different surface terminations are obtained in a
straightforward manner by removing half of a plane that contains two layers of
atoms, as visualized in Supplementary Figure S3. We emphasize that while
sampling terminations may suffice in some cases, it is clearly restricted to
atomic densities of $0$ and $0.5$, thereby performing very limited
optimization of the atomic structure. In our search, for example, the GRIP
tool finds a GB structure at $n=0$ with
$E_{\mathrm{gb}}=$0.509\text{\,}\mathrm{J}\text{\,}{\mathrm{m}}^{-2}$$ (green
circle), approximately $18\text{\,}\mathrm{\char 37\relax}$ lower in energy
than the best $\gamma$-surface structure (blue triangle). For comparison, the
ground-state structure and the two metastable states at $n=0$ are shown in
panels (e), (f), and (g).

Our structure searches performed for 150 GBs with three different tilt axes
show that the need for GCO and presence of multiple GB phases is a general
phenomenon in HCP Ti. Figure 4 summarizes the results from GRIP for the family
of $[0001]$ STGBs studied. Each subplot is equivalent to the gray boundary in
Figure 3a, denoting the minimum-energy structures at different $n$ and the
square marks the ground state. Evidently as many of the minima are located at
$n=0.5$, GCO is necessary to find the ground state in multiple $[0001]$ STGBs.
Similar to $\\{31\bar{4}0\\}$, the $\gamma$-surface method often performs
poorly for these GBs, getting higher energies and different structures than
the GRIP tool. The color map of the ground states reveals three distinct
intervals that correspond to different GB structural units. The low-angle GBs
in the intervals $\theta\leq$6.58\text{\,}\mathrm{\SIUnitSymbolDegree}$$ and
$\theta\geq$23.41\text{\,}\mathrm{\SIUnitSymbolDegree}$$ are composed of
isolated $\mathbf{b}_{\text{I}}$ edge dislocations (green markers) that for
the lowest angles do not require GCO. The near-energy-degenerate minima at
$n=0.5$ are composed of the same type of dislocations, with the extra atoms
accommodated by dislocation climb, resulting in GB structures with unevenly
spaced GB dislocations. Different GB dislocations stabilize at
$\theta\approx$10.89\text{\,}\mathrm{\SIUnitSymbolDegree}$$ with twice the
Burgers vector of $\mathbf{b}_{\text{II}}=\frac{1}{3}\langle
2\bar{4}20\rangle$ (purple markers). We investigate the transition between
these two states in this GB in detail in the next section. For high-angle GBs
in the interval
$$13.90\text{\,}\mathrm{\SIUnitSymbolDegree}$\leq\theta\leq$21.79\text{\,}\mathrm{\SIUnitSymbolDegree}$$,
the ground states at $n=0.5$ are composed of structural units that match the
dislocation core structures of $\mathbf{b}_{\text{II}}$, as outlined in
Supplementary Figure S4. Additional energy maps for select $[1\bar{1}00]$ and
$[1\bar{2}10]$ STGBs with low and high tilt angles out of 134 total studied
are shown in Figure 5, and results for the embedded-atom method (EAM)
potential [40] are presented in Supplementary Figures S5 and S6, and
Supplementary Note 1. We also perform select DFT calculations using the
optimized GRIP structures as inputs to confirm the stability of the GB
dislocation core structures and the relative energies between phases (see
Supplementary Figure S7). Taken together, these results demonstrate the
ubiquitous need for grand canonical sampling in locating the ground-state
structures for multiple tilt axes in HCP Ti.

Figure 5: GB energy maps of (a) $[1\bar{1}00]$ and (b) $[1\bar{2}10]$ STGBs.
Five different misorientations are selected for each tilt axis as
representative boundaries. Minima at $n=0.5$ indicate that GCO is broadly
required to find the ground states in all families of STGBs studied in
$\alpha$-Ti.

### Phase transitions and coexistence

The multiple GB phases predicted by GRIP opens up an opportunity to explore GB
phase transformations in HCP Ti; specifically, we focus on low-angle STGBs and
investigate transformations that change the topology of the dislocation
network arrangement. By elucidating the transformation pathways, we are able
to predict the structure of a nucleus with a distinct dislocation network
topology embedded inside a different parent dislocation network. We use point
defects to drive the transformation and we study the coupling between defect
absorption and changes in the dislocation network topology. While low-angle GB
phase transformations due to solutes and temperature have been previously
reported by experimental observations and simulations in a few metals [8, 41,
14, 42], the questions of transition states and the role of intrinsic point
defects have not been investigated.

Figure 6: GB phases of $\Sigma 7$$\\{21\bar{3}0\\}$$[0001]$ STGB. (a) The plot
of $E_{\mathrm{gb}}$ vs. $n$ reveals two GB phases, with the ground state at
$n=0$ and a metastable phase at $n=0.5$, shown in (b) and (c), respectively.
Both states are composed of edge dislocations indicated by green and purple
lines (identified using DXA in OVITO [39, 37]). The Burgers vectors of the
metastable boundary is twice that of the ground state. The non-HCP atoms of
the dislocation cores are colored according to the CNA.

We select the $\\{21\bar{3}0\\}$$[0001]$ STGB, which marks the transition
between the two different GB dislocation types at a misorientation angle of
$\theta\approx$10.9\text{\,}\mathrm{\SIUnitSymbolDegree}$$. The structure
search performed on this GB is illustrated in Figure 6, where we identify two
distinct GB phases corresponding to atomic densities $n=0$ (green circle) and
$n=0.5$ (orange square). Both structures correspond to GB energy cusps with
respect to $n$ and the structures of the two phases at $0\text{\,}\mathrm{K}$
are shown in Figure 6b and 6c, respectively. The ground state ($n=0$) is
composed of an array of $\mathbf{b}_{\text{I}}=\frac{1}{3}\langle
1\bar{2}10\rangle$ edge dislocations, while the second phase is composed of
$\mathbf{b}_{\text{II}}=\frac{1}{3}\langle 2\bar{4}20\rangle$ edge
dislocations with Burgers vector twice that of the ground state and
consequently half the line density within the GB plane. The optimized
dislocation core structures are consistent with the ‘T’ and ‘A’ structural
units, respectively, reported by Wang and Ye using constrained molecular
statics [29]. We perform MD simulations of each structure at temperatures as
high as $1150\text{\,}\mathrm{K}$ for up to $20\text{\,}\mathrm{ns}$ to
confirm that they are dynamically stable and indeed represent two GB phases.
The other two energy cusps at $n=0.33$ and $n=0.67$ are the mixed states
expected from the lever rule, where the GB region is patterned by weighted
fractions of $\mathbf{b}_{\text{I}}$ and $\mathbf{b}_{\text{II}}$ dislocations
corresponding to the proportions between $n=0$ and $n=0.5$, as shown in
Supplementary Figure S8.

Figure 7: Topological GB dislocation network transformation in
$\\{21\bar{3}0\\}$$[0001]$. (a) Adding Ti atoms (light blue) to the left half
of the ground-state structure and performing MD simulations at
$T=$600\text{\,}\mathrm{K}$$ triggers a dislocation-pairing transformation
$2\mathbf{b}_{\text{I}}\rightarrow\mathbf{b}_{\text{II}}$ in a quasi-2D
geometry. (b) Analogously, adding Ti atoms (dark red) to the left half of the
metastable structure triggers a dislocation-unpairing transition
$\mathbf{b}_{\text{II}}\rightarrow 2\mathbf{b}_{\text{I}}$. (c) Topological
transition of the GB dislocation network upon defect absorption. The view of
the GB plane shows a paired-dislocation GB island (nucleus) inside the parent
ground state.

Because the two GB phases are composed of different numbers of atoms, first-
order transitions between the two structures can be triggered by changing the
concentration of point defects [4]. The requisite high, local non-equilibrium
concentrations of vacancies or interstitials may occur, for example, as a
result of radiation damage, rapid cooling from high temperatures, or
deformation by creep. To mimic these conditions, we insert extra atoms into
interstitial sites in the ground-state structure, triggering a local
transformation of the GB structure illustrated in Figure 7. During the
transformation, the $\mathbf{b}_{\text{I}}$ dislocations of the ground-state
structure pair up into $\mathbf{b}_{\text{II}}$ dislocations and absorb the
extra atoms. Analogously, adding Ti atoms to the left half of the metastable
structure and performing high-temperature MD triggers a dislocation-unpairing
transition ($\mathbf{b}_{\text{II}}\rightarrow 2\mathbf{b}_{\text{I}}$) as
shown in Figure 7b. The transformed states remain stable at finite temperature
and the transformation can be reversed by introducing vacancies near the GB,
which we show in Supplementary Figure S9. Effectively, this sequence of states
and partial transformations illustrate the possibility of GB transformation-
mediated creep. Indeed, such a bicrystal can grow (shrink) by periodically
alternating its GB structure and absorbing only half a plane of atoms
(vacancies) at a time. If only one GB phase were present, the whole plane of
atoms would have to be absorbed in concert through disconnection motion before
returning to the original GB structure.

The simulated heterogeneous states containing two different GB phases show
stable coexistence in the closed system at high temperatures. While not
visible in Figure 7a and 7b, the two phases are separated by a line defect
called a GB phase junction, which is a dislocation as well as a force monopole
[43]. The Burgers vector of this junction is non-zero because the GB phases
have different dimensions [44]. The structure of this defect becomes more
apparent when considering nucleation in fully 3D. To illustrate the shape of
the nucleus during such a transformation, we increase the cross-section of the
GB and place interstitial atoms of Ti (light blue) in a relatively small
section. During the subsequent high-temperature simulation at
$T=$1000\text{\,}\mathrm{K}$$, the extra atoms diffuse to the boundary core
and locally trigger the pairing transition. The equilibrium structure of the
obtained nucleus is illustrated in Figure 7c. The transformation changes the
dislocation network topology as the dislocations of the parent structure shown
in green ($\mathbf{b}_{\text{I}}$) pair up on the nucleus boundary to form
three individual purple segments ($\mathbf{b}_{\text{II}}$). GB phase
nucleation by absorption of point defects has been previously investigated in
high-angle boundaries [4, 14]. The important distinction of the transformation
studied here is that it occurs in a low-angle GB; therefore, the core
structure of the GB phase junction is represented by a collection of
dislocation nodes where two dislocations pair up into one. To the best of our
knowledge, such transition states facilitating the change in the GB
dislocation network topology in a pure metal by point defect absorption have
not been previously reported.

## Discussion

In this work, we perform grand canonical GB structure search to discover new
GB phases in an HCP metal, $\alpha$-Ti. While GBs in $\alpha$-Ti have been
investigated extensively [29, 30, 31, 35, 32, 33, 34], prior simulations were
restricted to a fixed number of atoms derived from perfect surface
terminations with no point defects (see Supplementary Note 1). By rigorously
exploring atomic densities at GBs, we show that the minimum-energy structures
can be found for atomic densities inaccessible to the $\gamma$-surface method
for both high-angle and low-angle GBs across the misorientation range. The
ubiquitous need for GCO and presence of multiple GB phases with different
atomic densities is consistent with phenomena previously illustrated in
elemental cubic metals with FCC [4, 11, 12] and BCC [13, 14] crystal
structures.

Subsequent high-temperature MD simulations guided by this detailed sampling of
phase space yield the discovery of a novel GB phase transformation mechanism.
The two phases shown in Figure 6 are composed of periodic arrays of edge
dislocations with distinct localized cores that contain different atomic
densities in the GB. In the transition between these phases, dislocations of a
less dense GB ($\mathbf{b}_{\text{I}}$, $n=0$) pair up to form a new
dislocation core ($\mathbf{b}_{\text{II}}$, $n=0.5$), leading to a doubling of
the Burgers vector and absorption of interstitial atoms. Previous studies of
GB phase transitions in low-angle GBs revealed defect absorption by individual
dislocation cores without the change of the Burgers vector [11, 14] and other
studies demonstrated the change in the dislocation network topology due to
temperature [45] and solute segregation [41]. It is also well established and
expected that individual dislocations absorb point defects by climb [46]; yet
here, we demonstrate a different mechanism where the dislocation network
topology and the number of constituent atoms are coupled. This coupling
suggests an important mechanism for point-defect absorption in polycrystalline
materials with non-equilibrium concentrations of point defects produced by
rapid quenching, irradiation, or additive manufacturing approaches that can
yield dense dislocation cellular walls [47]. This work thus provides important
insights into the ways in which low-angle GBs and dislocation arrays interact
with point defects [48]. The work here, focused on an HCP metal, may be
particularly relevant for engineering materials with such structures that
experience radiation damage, such as Zr-based nuclear fuel cladding [49].

Herein, we further extend the notion of the number of atoms in a GB plane
($N^{\mathrm{bulk}}_{\mathrm{plane}}$) to non-cubic, multi-basis crystals like
HCP. Previous studies on elemental cubic metals calculated this quantity as
the total number of atoms located in one planar cut parallel to the GB, i.e.,
all these atoms are equidistant in the $z$-direction. This is not always the
case for HCP metals or any multi-basis crystal, as visualized in Figure 2.
Generally, $N^{\mathrm{bulk}}_{\mathrm{plane}}$ includes all atoms located
inside a region with height equal to the minimum normal component of a lattice
vector. In this work, we show that if the GB structure search considers only
those terminations of the surface by complete atomic layers, it is restricted
to sampling states with $n=0$ or $n=0.5$ solely. Importantly, we demonstrate
that this restriction misses lower-energy GB structures with intermediate
values of $n$. A thorough search must consider all different atomic densities,
as generalized by the framework presented here.

We implement this framework for handling the structural DOF and predicting new
GB phases in the open-source GRIP tool, written in Python with minimal
dependencies (see Code availability). The algorithm rapidly samples the
configurational space described by relative translations and different atomic
densities and moves the system toward equilibrium. The relevant DOF—e.g.,
atomic density, reconstructions, temperature—are specified by the user in a
single input file and the code exhaustively explores the GB phase space by
sampling as many structures as possible in parallel. The energy calculations
presented here use empirical interatomic potentials (IAPs) to perform the
dynamic sampling, but other techniques such as DFT can be used as well, as
those calculations are decoupled from the structure optimization steps;
however, the use of IAPs enables us to access low-angle GBs and larger
reconstructions with thousands of atoms, as demonstrated here in simulations
up to $3\,\times\,13$ reconstructions to validate the dislocation character.
This methodology can thus take advantage of the increasing availability of
computational resources and the advent of high-fidelity machine-learned IAPs
to enable quantum-accurate atomistic simulations of large systems with
extended defects [50]. Advances in sampling and structure generation
algorithms will further expand the diversity of results and the modular
structure of the code enables different techniques to be easily plugged in. Of
particular interest would be extensions to multicomponent systems, which could
be handled using a Monte Carlo approach [21, 22] for compositional DOF and
would enable grand canonical sampling of GB structures in technologically
relevant alloy chemistries.

## Methods

### GB structure search

Figure 8: Simulation cell setup for GRIP. The cell is oriented such that the
$y$-axis is the tilt axis direction, the $x$-axis is the orthogonal in-plane
direction, and the $z$-axis is the out-of-plane normal direction. Periodic
boundary conditions are maintained in the GB plane ($xy$-plane). GRIP begins
by translating the upper crystal and removing atoms from the GB (magenta
line). During MD, atoms in the GB regions are free to move while atoms in the
buffer and semi-rigid regions in the upper slab are constrained to move
together, and atoms in the lower two regions are fixed. During relaxation,
atoms in both buffer regions are free to move while the semi-rigid region is
still constrained. For the $\gamma$-surface method, the upper slab is only
allowed to translate as a whole before relaxation is applied.

We perform atomic-level optimization of GB structures using the open-source,
Python-based GRand canonical Interface Predictor (GRIP) code (see Code
availability), which rigorously explores structural DOF through dynamic
sampling. Bicrystal slabs can be automatically generated using the Atomic
Simulation Environment (ASE) [51] library or supplied as external files, which
we created for $\alpha$-Ti using a combination of ASE and Pymatgen [52].
Figure 8 shows the orientation of the simulation cell, such that the $y$-axis
is the tilt axis direction, the $x$-axis is the orthogonal in-plane direction,
and the $z$-axis is the out-of-plane normal direction. We ensure periodicity
in the GB plane ($xy$-plane) and an integer multiple of the interplanar
spacing that totals at least $3.5\text{\,}\mathrm{nm}$ in the $z$-direction
for each slab to minimize cell size effects. While any size cell can be used
in principle, for computational tractability in this high-throughput study, we
choose to simulate only $[0001]$, $[1\bar{1}00]$, and $[1\bar{2}10]$ STGBs
where all Miller-Bravais indices for the plane and the in-plane $x$-direction
are less than or equal to 15, resulting in 16, 40, and 94 STGBs for each of
the tilt axes, respectively (150 total).

For an individual GB, each iteration of the algorithm has three stages. During
the first stage, the initial configuration is created by uniformly sampling a
specific set of GB DOF. Specifically, the algorithm randomly samples an
$m\,\times\,n$ replication of the unit GB cell (here, up to $3\,\times\,3$),
randomly translates the upper slab in the $xy$-plane, and removes a randomly
chosen (from the user-specified interval) fraction of atoms from the GB. To
further increase the structural diversity of the initial GB configurations, we
have implemented random swaps of atoms on crystal lattice sites and
interstitial sites in the GB region. The algorithm identifies interstitial
sites near the GB as the vertices of the Voronoi diagram of the GB region.

During the second stage, it performs dynamic sampling to optimize the GB
structure consistent with the imposed DOF. In this study, we used standard
finite-temperature MD simulations using the Large-scale Atomic/Molecular
Massively Parallel Simulator (LAMMPS) [53] in the canonical ($NVT$) ensemble
with a Langevin thermostat and a time step of $2\text{\,}\mathrm{fs}$ as our
dynamic sampling technique. The temperature and duration of the MD are also
randomly sampled based on the user-specified ranges and only the atoms in the
GB region are allowed to move freely during the dynamical sampling phase. It
is straightforward to substitute this MD optimization with other, more
sophisticated sampling techniques implemented in LAMMPS or other codes. For
this study, we choose a GB region of $1\text{\,}\mathrm{nm}$ thickness on each
side, a temperature between $300\text{\,}\mathrm{K}$ and
$1200\text{\,}\mathrm{K}$, and a duration up to $0.6\text{\,}\mathrm{ns}$.
Finally, the temperature is quickly ramped down to $100\text{\,}\mathrm{K}$
for $2\text{\,}\mathrm{ps}$. With some probability (here,
$5\text{\,}\mathrm{\char 37\relax}$) the algorithm skips the dynamic sampling
for one iteration and jumps to the third stage.

In the third stage, each GB structure following MD sampling is fully relaxed
at $0\text{\,}\mathrm{K}$ using a conjugate gradient minimization scheme,
where atoms in the GB and buffer regions can move freely while the semi-rigid
region is constrained to move together. Here, we specified the buffer region
to be $0.6\text{\,}\mathrm{nm}$ beyond each side of the GB region, and larger
values lowered $E_{\mathrm{gb}}$ by no more than $1\text{\,}\mathrm{\char
37\relax}$. The convergence criteria are ${10}^{-15}\text{\,}\mathrm{f}$or
relative energy ($\mathrm{d}E/E$ in successive iterations) and
${10}^{-15}\text{\,}\mathrm{eV}\text{\,}{\mathrm{\text{Å}}}^{-1}$ for forces,
with a maximum of ${10}^{5}\text{\,}\mathrm{e}$valuations for each criterion.
The algorithm repeats these stages on each processor independently until
termination, saving each relaxed structure to disk and periodically deleting
duplicates. Duplicates are defined as structures with the same value of
$E_{\mathrm{gb}}$ and $n$ to three decimal places, and the algorithm will keep
the structure with a smaller reconstruction and relative translations.

For each relaxed structure, the GB energy, $E_{\mathrm{gb}}$, is computed
according to:

$E_{\mathrm{gb}}=\frac{E^{\mathrm{gb}}_{\mathrm{total}}-N^{\mathrm{gb}}_{\mathrm{total}}E^{\mathrm{bulk}}_{\mathrm{coh}}}{A^{\mathrm{gb}}_{\mathrm{plane}}}$
(1)

where $E^{\mathrm{gb}}_{\mathrm{total}}$ and
$N^{\mathrm{gb}}_{\mathrm{total}}$ are the total energy and number,
respectively, of atoms in the GB and buffer regions,
$E^{\mathrm{bulk}}_{\mathrm{coh}}$ is the cohesive energy per atom in bulk
$\alpha$-Ti, and $A^{\mathrm{gb}}_{\mathrm{plane}}$ is the area of the GB
plane. We also track the fraction of atoms in one plane or GB atomic density,
$n$, according to:

$n=\frac{N_{\mathrm{total}}\mod
N^{\mathrm{bulk}}_{\mathrm{plane}}}{N^{\mathrm{bulk}}_{\mathrm{plane}}}\in[0,1)$
(2)

where $N_{\mathrm{total}}$ is the total number of atoms in the simulation cell
and $N^{\mathrm{bulk}}_{\mathrm{plane}}$ is the number of atoms in one plane
of the bulk structure. Previous calculations of
$N^{\mathrm{bulk}}_{\mathrm{plane}}$ simply counted the number of atoms at a
single $z$ value in the bulk [4, 11]; however, due to the 2-atom basis of the
HCP crystal structure, atoms associated with one plane may be offset in the
$z$-direction, as we show in Figure 2. Therefore, we calculate
$N^{\mathrm{bulk}}_{\mathrm{plane}}$ as the number of atoms within a region
equal to the minimum normal component of a lattice vector; in HCP $\alpha$-Ti,
this is equivalent to the interplanar spacing of the hexagonal lattice
($d_{hkl}$) given by [54]:

$\frac{1}{d_{hkl}^{2}}=\frac{4}{3}\left(\frac{h^{2}+hk+k^{2}}{a^{2}}\right)+\left(\frac{l}{c}\right)^{2}$
(3)

where $h$, $k$, and $l$ are the Miller indices, and $a$ and $c$ are the HCP
lattice constants. We note this extended definition of
$N^{\mathrm{bulk}}_{\mathrm{plane}}$ reduces to taking a planar slice for
unary, single-basis systems like elemental BCC and FCC metals, consistent with
previous studies [4, 11].

We compare the results of structure optimization using two different
interatomic potentials, an embedded-atom method (EAM) potential for Ti–Al from
Zope and Mishin [40] and a modified embedded-atom method (MEAM) potential for
Ti from Hennig, et al. [38] For each STGB and potential, we also optimize the
structure using the $\gamma$-surface method [19] for comparison, using a
$2\,\times\,4$ replication of the same bicrystals and translating the top slab
in increments of $0.025\text{\,}\mathrm{nm}$ in the $x$\- and $y$-directions
prior to a conjugate gradient energy minimization.

### High-temperature MD simulations

To study GB phase stability and transitions, we perform high-temperature MD
simulations using methods adapted from previous work [4]. Briefly, we
replicate the STGB in the $x$\- and $y$-directions until the simulation cell
is around $10\text{\,}\mathrm{nm}$ in the $x$-direction and
$3\text{\,}\mathrm{nm}$ in the $y$-direction along the tilt axis. We freeze
the bottom $1\text{\,}\mathrm{nm}$ layer of atoms and constrain the top
$1\text{\,}\mathrm{nm}$ layer to be semi-rigid throughout the
$20\text{\,}\mathrm{ns}$ simulation. We use periodic boundary conditions
(PBCs) in the $y$-direction and both PBCs and open surfaces with
$1\text{\,}\mathrm{nm}$ of vacuum in the $x$-direction. We scan a range of
temperatures between $600\text{\,}\mathrm{K}$ and $1200\text{\,}\mathrm{K}$.

To induce a phase transition, we either insert additional Ti atoms at
interstitial sites in the GB region or delete Ti atoms from a region near the
top of the GB region. These MD simulations are performed in the canonical
($NVT$) ensemble between $600\text{\,}\mathrm{K}$ and
$1200\text{\,}\mathrm{K}$ for up to $20\text{\,}\mathrm{ns}$, using the MEAM
potential and associated structures. For clarity of visualization, we relax
all structures at $0\text{\,}\mathrm{K}$ using a conjugate gradient
minimization scheme.

### DFT calculations

To validate select GB structures, we perform additional density functional
theory (DFT) calculations using the Vienna Ab initio Simulation Package (VASP)
[55, 56, 57, 58] with projector augmented-wave potentials [59] and the
generalized gradient approximation exchange correlation functional of Perdew,
Burke, and Ernzerhof [60]. The semi-core $3p$ states are treated as valence
states (Ti_pv potential). We use Monkhorst-Pack [61] $\mathbf{k}$-point grids
with a density of 5000 $\mathbf{k}$ points per reciprocal atom and apply
Methfessel-Paxton smearing [62] with a width of $0.1\text{\,}\mathrm{eV}$. The
plane wave cutoff energy is $500\text{\,}\mathrm{eV}$ and the convergence
criteria are set at ${10}^{-5}\text{\,}\mathrm{eV}$ for energy and
$0.02\text{\,}\mathrm{eV}\text{\,}{\mathrm{\text{Å}}}^{-1}$ for forces. We
create the input structure by extracting a section near the GB region of the
optimized structure from GRIP of approximately $4.5\text{\,}\mathrm{nm}$ in
thickness (200–300 atoms) and adding $1\text{\,}\mathrm{nm}$ of vacuum on top.
The axes are rescaled to match equilibrium DFT values and atomic positions are
fully relaxed while the cell shape and volume are fixed to maintain stresses
in the GB plane. The energy of the GB is computed as the difference in total
energy of a structure with the GB and a second bulk structure in the same
orientation with the same number of atoms and vacuum but without a GB, divided
by the planar area.

## Data availability

The data that support the findings of this study, including input and relaxed
STGB structures, are available at Zenodo at publication time. Other data are
available from the corresponding author upon reasonable request.

## Code availability

The GRand canonical Interface Predictor (GRIP) tool that implements the GB
structure optimization algorithm described here can be found at
https://github.com/enze-chen/grip at publication time.

## Acknowledgments

This work was performed under the auspices of the U.S. Department of Energy
(DOE) by the Lawrence Livermore National Laboratory (LLNL) under Contract No.
DE-AC52-07NA27344. Computing support for this work comes from the LLNL
Institutional Computing Grand Challenge program and Bridges-2 at Pittsburgh
Supercomputing Center through allocation DMR110087 from the Advanced
Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS)
program. T.F. acknowledges support from the U.S. DOE, Office of Science under
an Office of Fusion Energy Sciences Early Career Award. E.C. and M.A.
acknowledge funding through the DOE, Office of Science, Office of Basic Energy
Sciences, Materials Sciences and Engineering Division, under Contract No. DE-
AC02-05-CH11231 within the Materials Project program (KC23MP). E.C. also
acknowledges a fellowship through the National Science Foundation Graduate
Research Fellowship Program under Grant No. DGE-2146752 and support from the
Computational Chemistry & Materials Science Institute at LLNL. Part of this
work was funded by the Laboratory Directed Research and Development (LDRD)
program at LLNL under projects with a tracking code 22-ERD-002. Part of this
work was also supported by the U.S. DOE, Office of Energy Efficiency and
Renewable Energy, Hydrogen and Fuel Cell Technologies Office, through the
Hydrogen storage Materials Advanced Research Consortium (HyMARC). The authors
thank Tomas Oppelstrup and Daryl Chrzan for helpful discussions. All figures
are produced using matplotlib[63] and OVITO [37].

## Author contributions

T.F. and E.C. designed the study which was initially supervised by T.W.H. and
B.C.W. T.F. and E.C. designed the GCO algorithm. E.C. performed the
simulations and analyzed the results. E.C. and T.F. drafted the manuscript
with feedback from M.A. All authors discussed the results and contributed to
the writing of the manuscript.

## Declaration of competing interests

The authors declare no competing financial or non-financial interests.

## Materials and correspondence

Correspondence and requests for materials should be addressed to E.C. or T.F.

## References

  * [1] A. P. Sutton and R. W. Balluffi.  Interfaces in Crystalline Materials.  Oxford University Press, New York, 1st ed., (1995).
  * [2] T. Frolov and Y. Mishin.  Phases, phase equilibria, and phase rules in low-dimensional systems.  J. Chem. Phys. 143, 044706 (2015).
  * [3] P. R. Cantwell, M. Tang, S. J. Dillon, et al.  Grain boundary complexions.  Acta Mater. 62, 1–48 (2014).
  * [4] T. Frolov, D. L. Olmsted, M. Asta, and Y. Mishin.  Structural phase transformations in metallic grain boundaries.  Nat. Commun. 4, 1–7 (2013).
  * [5] P. R. Cantwell, T. Frolov, T. J. Rupert, et al.  Grain boundary complexion transitions.  Annu. Rev. Mater. Res. 50, 465–492 (2020).
  * [6] T. Frolov, S. V. Divinski, M. Asta, and Y. Mishin.  Effect of interface phase transformations on diffusion and segregation in high-angle grain boundaries.  Phys. Rev. Lett. 110, 255502 (2013).
  * [7] J. Wei, B. Feng, R. Ishikawa, et al.  Direct imaging of atomistic grain boundary migration.  Nat. Mater. 20, 951–955 (2021).
  * [8] T. Meiners, T. Frolov, R. E. Rudd, G. Dehm, and C. H. Liebscher.  Observations of grain-boundary phase transformations in an elemental metal.  Nature 579, 375–378 (2020).
  * [9] S. Rajeshwari K., S. Sankaran, K. C. Hari Kumar, et al.  Grain boundary diffusion and grain boundary structures of a Ni–Cr–Fe-alloy: Evidences for grain boundary phase transitions.  Acta Mater. 195, 501–518 (2020).
  * [10] T. Frolov, M. Asta, and Y. Mishin.  Phase transformations at interfaces: Observations from atomistic modeling.  Curr. Opin. Solid St. M. 20, 308–315 (2016).
  * [11] Q. Zhu, A. Samanta, B. Li, R. E. Rudd, and T. Frolov.  Predicting phase behavior of grain boundaries with evolutionary search and machine learning.  Nat. Commun. 9, 467 (2018).
  * [12] T. Brink, L. Langenohl, H. Bishara, and G. Dehm.  Universality of grain boundary phases in fcc metals: Case study on high-angle [111] symmetric tilt grain boundaries.  Phys. Rev. B. 107, 054103 (2023).
  * [13] T. Frolov, W. Setyawan, R. J. Kurtz, et al.  Grain boundary phases in bcc metals.  Nanoscale 10, 8253–8268 (2018).
  * [14] T. Frolov, Q. Zhu, T. Oppelstrup, J. Marian, and R. E. Rudd.  Structures and transitions in bcc tungsten grain boundaries and their role in the absorption of point defects.  Acta Mater. 159, 123–134 (2018).
  * [15] S. von Alfthan, P. D. Haynes, K. Kaski, and A. P. Sutton.  Are the structures of twist grain boundaries in silicon ordered at 0 K?  Phys. Rev. Lett. 96, 055505 (2006).
  * [16] J. Zhang, C.-Z. Wang, and K.-M. Ho.  Finding the low-energy structures of Si[001] symmetric tilted grain boundaries with a genetic algorithm.  Phys. Rev. B. 80, 174102 (2009).
  * [17] A. L.-S. Chua, N. A. Benedek, L. Chen, M. W. Finnis, and A. P. Sutton.  A genetic algorithm for predicting the structures of interfaces in multicomponent systems.  Nat. Mater. 9, 418–422 (2010).
  * [18] A. B. Mazitov and A. R. Oganov.  Grain boundaries in minerals: Atomic structure, phase transitions, and effect on strength of polycrystals.  Zapiski Rmo (Proc. Russ. Mineral. Soc.) 150, 92–102 (2021).
  * [19] Y. Mishin and D. Farkas.  Atomistic simulation of [001] symmetrical tilt grain boundaries in NiAl.  Philos. Mag. A 78, 29–56 (1998).
  * [20] J. Han, V. Vitek, and D. J. Srolovitz.  Grain-boundary metastability and its statistical properties.  Acta Mater. 104, 259–273 (2016).
  * [21] A. D. Banadaki, M. A. Tschopp, and S. Patala.  An efficient Monte Carlo algorithm for determining the minimum energy structures of metallic grain boundaries.  Comput. Mater. Sci. 155, 466–475 (2018).
  * [22] S. Yang, N. Zhou, H. Zheng, S. P. Ong, and J. Luo.  First-order interfacial transformations with a critical point: Breaking the symmetry at a symmetric tilt grain boundary.  Phys. Rev. Lett. 120, 085702 (2018).
  * [23] S. E. Restrepo, S. T. Giraldo, and B. J. Thijsse.  A genetic algorithm for generating grain boundaries.  Model. Simul. Mater. Sci. Eng. 21, 055017 (2013).
  * [24] C. Yang, M. Zhang, and L. Qi.  Grain boundary structure search by using an evolutionary algorithm with effective mutation methods.  Comput. Mater. Sci. 184, 109812 (2020).
  * [25] G. Lütjering and J. C. Williams.  Titanium.  Springer, Berlin, 2nd ed., (2007).
  * [26] Y. Chong, R. Gholizadeh, T. Tsuru, et al.  Grain refinement in titanium prevents low temperature oxygen embrittlement.  Nat. Commun. 14, 404 (2023).
  * [27] M. N. Kelly, K. Glowinski, N. T. Nuhfer, and G. S. Rohrer.  The five parameter grain boundary character distribution of $\alpha$-Ti determined from three-dimensional orientation data.  Acta Mater. 111, 22–30 (2016).
  * [28] M. S. Hooshmand, R. Zhang, Y. Chong, et al.  Twin-boundary structural phase transitions in elemental titanium.  arXiv:2103.06194 (2021).
  * [29] Y. C. Wang and H. Q. Ye.  On the tilt grain boundaries in hcp Ti with [0001] orientation.  Philos. Mag. A 75, 261–272 (1997).
  * [30] Z. Zheng, D. S. Balint, and F. P. E. Dunne.  Investigation of slip transfer across HCP grain boundaries with application to cold dwell facet fatigue.  Acta Mater. 127, 43–53 (2017).
  * [31] D. Farkas.  Grain-boundary structures in hexagonal materials: Coincident and near coincident grain boundaries.  Metall. Mater. Trans. A 25, 1337–1346 (1994).
  * [32] J. Wang and I. J. Beyerlein.  Atomic structures of $[0\bar{1}10]$ symmetric tilt grain boundaries in hexagonal close-packed (hcp) crystals.  Metall. Mater. Trans. A 43, 3556–3569 (2012).
  * [33] C. Ni, H. Ding, M. Asta, and X. Jin.  Computational study of $\langle 1\bar{1}00\rangle$ symmetric tilt grain boundaries in Mg and Ti.  Scripta Mater. 109, 94–99 (2015).
  * [34] M. A. Bhatia and K. N. Solanki.  Energetics of vacancy segregation to symmetric tilt grain boundaries in hexagonal closed pack materials.  J. Appl. Phys. 114, 244309 (2013).
  * [35] J. Wang and I. J. Beyerlein.  Atomic structures of symmetric tilt grain boundaries in hexagonal close packed (hcp) crystals.  Model. Simul. Mater. Sci. Eng. 20, 024002 (2012).
  * [36] A. Stukowski.  Structure identification methods for atomistic simulations of crystalline materials.  Model. Simul. Mater. Sci. Eng. 20, 045021 (2012).
  * [37] A. Stukowski.  Visualization and analysis of atomistic simulation data with OVITO—the Open Visualization Tool.  Model. Simul. Mater. Sci. Eng. 18, 015012 (2009).
  * [38] R. G. Hennig, T. J. Lenosky, D. R. Trinkle, S. P. Rudin, and J. W. Wilkins.  Classical potential describes martensitic phase transformations between the $\alpha$, $\beta$, and $\omega$ titanium phases.  Phys. Rev. B. 78, 054121 (2008).
  * [39] A. Stukowski, V. V. Bulatov, and A. Arsenlis.  Automated identification and indexing of dislocations in crystal interfaces.  Model. Simul. Mater. Sci. Eng. 20, 085007 (2012).
  * [40] R. R. Zope and Y. Mishin.  Interatomic potentials for atomistic simulations of the Ti–Al system.  Phys. Rev. B. 68, 024102 (2003).
  * [41] K. E. Sickafus and S. L. Sass.  Grain boundary structural transformations induced by solute segregation.  Acta Metall. 35, 69–79 (1987).
  * [42] N. J. Peter, T. Frolov, M. J. Duarte, et al.  Segregation-induced nanofaceting transition at an asymmetric tilt grain boundary in copper.  Phys. Rev. Lett. 121, 255502 (2018).
  * [43] I. S. Winter, R. E. Rudd, T. Oppelstrup, and T. Frolov.  Nucleation of grain boundary phases.  Phys. Rev. Lett. 128, 035701 (2022).
  * [44] T. Frolov, D. L. Medlin, and M. Asta.  Dislocation content of grain boundary phase junctions and its relation to grain boundary excess properties.  Phys. Rev. B. 103, 184108 (2021).
  * [45] D. L. Olmsted, D. Buta, A. Adland, et al.  Dislocation-pairing transitions in hot grain boundaries.  Phys. Rev. Lett. 106, 046101 (2011).
  * [46] Y. Komem, P. Pétroff, and R. W. Balluffi.  Direct observation of grain boundary dislocation climb in ion-irradiated gold bicrystals.  Philos. Mag. 26, 239–252 (1972).
  * [47] Z. Shang, C. Fan, S. Xue, et al.  Response of solidification cellular structures in additively manufactured 316 stainless steel to heavy ion irradiation: an _in situ_ study.  Mater. Res. Lett. 7, 290–297 (2019).
  * [48] I. Beyerlein, M. Demkowicz, A. Misra, and B. Uberuaga.  Defect-interface interactions.  Prog. Mater. Sci. 74, 125–210 (2015).
  * [49] C. Yan, R. Wang, Y. Wang, X. Wang, and G. Bai.  Effects of ion irradiation on microstructure and properties of zirconium alloys—A review.  Nucl. Eng. Technol. 47, 323–331 (2015).
  * [50] R. Freitas and Y. Cao.  Machine-learning potentials for crystal defects.  MRS Commun. 12, 510–520 (2022).
  * [51] A. H. Larsen, J. J. Mortensen, J. Blomqvist, et al.  The atomic simulation environment—a Python library for working with atoms.  J. Phys. Condens. Matter 29, 273002 (2017).
  * [52] S. P. Ong, W. D. Richards, A. Jain, et al.  Python Materials Genomics (pymatgen): A robust, open-source python library for materials analysis.  Comput. Mater. Sci. 68, 314–319 (2013).
  * [53] A. P. Thompson, H. M. Aktulga, R. Berger, et al.  LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales.  Comput. Phys. Commun. 271, 108171 (2022).
  * [54] B. D. Cullity and S. R. Stock.  Elements of X-Ray Diffraction.  Pearson, 3rd ed., (2001).
  * [55] G. Kresse and J. Hafner.  _Ab initio_ molecular dynamics for liquid metals.  Phys. Rev. B. 47, 558–561 (1993).
  * [56] G. Kresse and J. Hafner.  _Ab initio_ molecular-dynamics simulation of the liquid-metal–amorphous-semiconductor transition in germanium.  Phys. Rev. B. 49, 14251–14269 (1994).
  * [57] 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).
  * [58] G. Kresse and J. Furthmüller.  Efficiency of _ab initio_ total energy calculations for metals and semiconductors using a plane-wave basis set.  Comput. Mater. Sci. 6, 15–50 (1996).
  * [59] G. Kresse and D. Joubert.  From ultrasoft pseudopotentials to the projector augmented-wave method.  Phys. Rev. B. 59, 1758–1775 (1999).
  * [60] J. P. Perdew, K. Burke, and M. Ernzerhof.  Generalized gradient approximation made simple.  Phys. Rev. Lett. 77, 3865–3868 (1996).
  * [61] H. J. Monkhorst and J. D. Pack.  Special points for Brillouin-zone integrations.  Phys. Rev. B. 13, 5188–5192 (1976).
  * [62] M. Methfessel and A. T. Paxton.  High-precision sampling for Brillouin-zone integration in metals.  Phys. Rev. B. 40, 3616–3621 (1989).
  * [63] J. D. Hunter.  Matplotlib: A 2D graphics environment.  Comput. Sci. Eng. 9, 90–95 (2007).
  * [64] Y. Mishin, M. J. Mehl, D. A. Papaconstantopoulos, A. F. Voter, and J. D. Kress.  Structural stability and lattice defects in copper: _Ab initio_ , tight-binding, and embedded-atom calculations.  Phys. Rev. B. 63, 224106 (2001).
  * [65] X. Zhou, R. A. Johnson, and H. N. G. Wadley.  Misfit-energy-increasing dislocations in vapor-deposited CoFe/NiFe multilayers.  Phys. Rev. B. 69, 144113 (2004).
  * [66] F. H. Stillinger and T. A. Weber.  Computer simulation of local order in condensed phases of silicon.  Phys. Rev. B. 31, 5262–5271 (1985).
  * [67] A. R. Oganov and C. W. Glass.  Crystal structure prediction using _ab initio_ evolutionary techniques: Principles and applications.  J. Chem. Phys. 124, 244704 (2006).
  * [68] A. O. Lyakhov, A. R. Oganov, H. T. Stokes, and Q. Zhu.  New developments in evolutionary structure prediction algorithm USPEX.  Comput. Phys. Commun. 184, 1172–1182 (2013).
  * [69] G. J. Ackland.  Theoretical study of titanium surfaces and defects with a new many-body potential.  Philos. Mag. A 66, 917–932 (1992).
  * [70] N. A. Zarkevich and D. D. Johnson.  Titanium $\alpha$–$\omega$ phase transformation pathway and a predicted metastable structure.  Phys. Rev. B. 93, 020104 (2016).
  * [71] B. Waters, D. S. Karls, I. Nikiforov, et al.  Automated determination of grain boundary energy and potential-dependence using the OpenKIM framework.  Comput. Mater. Sci. 220, 112057 (2023).
  * [72] E. Torres.  Atomistic study of the structure and deformation behavior of symmetric tilt grain boundaries in $\alpha$-zirconium.  Comput. Mater. Sci. 197, 110600 (2021).

Supplementary Information for
Grand canonically optimized grain boundary phases in hexagonal close-packed
titanium

Enze Chen1,2,3,4,∗, Tae Wook Heo2, Brandon C. Wood2, Mark Asta1,3, Timofey
Frolov2,∗

1 Department of Materials Science and Engineering, University of California,
Berkeley, CA 94720, USA
2 Materials Science Division, Lawrence Livermore National Laboratory,
Livermore, CA 94550, USA
3 Materials Sciences Division, Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, USA
4 Present address: Department of Materials Science and Engineering, Stanford
University, Stanford, CA 94305, USA
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
Figure S1: Validation of the GRIP tool. We perform grand canonical
optimization of the (a) $\Sigma 5(210)[001]$ tilt GB in Cu (EAM potential
[64]), (b) $\Sigma 5(310)[001]$ twist GB in W (EAM [65]), (c) $\Sigma
5(310)[001]$ tilt and (d) twist GB in Si (Stillinger-Weber [66]). The
$E_{\mathrm{gb}}$ vs. $n$ plots and low-energy structures match those in the
literature for Cu [4], W [14], and Si [21, 15]. The atoms are colored
according to CNA [36], where gray are bulk-coordinated atoms and brown are
non-bulk-coordinated atoms (in the GB).

Figure S2: Discovering a new ground state in W$(118)[1\bar{1}0]$. A previous
study [14] using the evolutionary algorithm USPEX [67, 68] and an EAM
potential [65] found a ground state for the W$(118)[1\bar{1}0]$ STGB at
$n=0.33$ with
$E_{\mathrm{gb}}=$0.225\text{\,}\mathrm{J}\text{\,}{\mathrm{m}}^{-2}$$, as
pictured on the top. Using GRIP and the same EAM potential, we find a
structure with lower energy
$E_{\mathrm{gb}}=$0.220\text{\,}\mathrm{J}\text{\,}{\mathrm{m}}^{-2}$$ at
$n=0.67$, as shown on the bottom. This structure has a different dislocation
network where the $\langle 001\rangle$-type edge dislocations (magenta, as
identified using DXA) overlap in the GB plane instead of residing in parallel.

Figure S3: Different planar terminations and corresponding GB structures.
$\\{31\bar{4}0\\}$$[0001]$ is a GB where the two basis atoms have different
$z$ coordinates normal to the GB plane, and two initial configurations are
shown where the lower slabs terminate at different basis atoms (gold for the
first, blue for the second). We use the $\gamma$-surface method to optimize
the GB structure in both cases. The first configuration ($n=0$) produces a
higher-energy structure shown on the left
($E_{\mathrm{gb}}=$0.625\text{\,}\mathrm{J}\text{\,}{\mathrm{m}}^{-2}$$),
while the second configuration ($n=0.5$) produces the lower-energy structure
shown on the right
($E_{\mathrm{gb}}=$0.351\text{\,}\mathrm{J}\text{\,}{\mathrm{m}}^{-2}$$),
which matches the ground-state structure from GRIP. The final structures and
colors correspond to those in Figure 3 in the main manuscript.

Figure S4: GB structural units in $[0001]$ STGBs. GRIP results for select GBs
at $n=0.0$ and $n=0.5$ viewed down $[0001]$. Low-angle GBs adopt a dislocation
core configuration at $n=0.0$ that is accommodated at $n=0.5$ through
dislocation climb. At
$\theta\approx$10.89\text{\,}\mathrm{\SIUnitSymbolDegree}$$, there is a
transition at $n=0.5$ to a different structural unit. At even higher angles
$\theta\geq$22.41\text{\,}\mathrm{\SIUnitSymbolDegree}$$, the GB structural
units transform back into the motifs at low angles. The lower-energy structure
at each tilt angle, as evaluated using the MEAM potential [38], is outlined in
black.

Figure S5: Energy map of $[0001]$ STGBs simulated using an EAM potential [40].
The profiles are qualitatively similar to those in Figure 4 in the main text,
which was generated using a MEAM potential [38]. See Figure 4 for additional
descriptions of features.

Figure S6: Comparison of optimized GB structures. A few STGBs are provided for
each tilt axis and two projections are shown for each ground state obtained
using the EAM [40] and MEAM [38] potentials. The GB atomic density ($n$) is
shown in the lower-right corner and may not be equal for both potentials. The
atoms are colored according to the common neighbor analysis (CNA) in OVITO
[36, 37], where gray are HCP-coordinated atoms, gold are FCC-coordinated
atoms, and blue have a different coordination (in the GB).

Figure S7: DFT validation for select STGBs. The structures from GRIP are used
as inputs to VASP for further relaxation (see Methods in the main manuscript
for details). The GB dislocation core structures remain stable and the
relative energies between different phases are consistent with those from
GRIP.

Figure S8: Mixed state at $n=0.33$ and $n=0.67$ for $\\{21\bar{3}0\\}$. At
intermediate values of $n\in(0,0.5)$, the localized dislocation cores in the
GB alternate between $\mathbf{b}_{\text{I}}=\frac{1}{3}\langle
1\bar{2}10\rangle$ (green) and $\mathbf{b}_{\text{II}}=\frac{1}{3}\langle
2\bar{4}20\rangle$ (purple) dislocations to obtain the minimum-energy
structure. The amount of each phase follows the conventional lever rule for
phase fractions. Dislocations are identified using the dislocation extraction
algorithm (DXA) in OVITO [39, 37].  Figure S9: (Reverse) phase transformation
through vacancy absorption. In the main text, we demonstrate interstitial-
induced phase transformation and coexistence in $\\{21\bar{3}0\\}$$[0001]$,
where every two $\mathbf{b}_{\text{I}}$ dislocations pair up to form one
$\mathbf{b}_{\text{II}}$ dislocation. Here, by injecting vacancies (outlined
in dark red) and performing MD simulations at $1200\text{\,}\mathrm{K}$, we
reverse the transformation, whereby the $\mathbf{b}_{\text{II}}$ dislocations
absorb the vacancies and unpair.

## Supplementary Note 1

Figure S10: $E_{\mathrm{gb}}$ vs. $\theta$ for STGBs in $\alpha$-Ti. All
results are shown for the (a) $[0001]$, (b) $[1\bar{1}00]$, and (c)
$[1\bar{2}10]$ tilt axes. Each point corresponds to the minimum-energy
structure for that tilt angle. The orange and blue lines correspond to our
calculations using the EAM [40] and MEAM [38] potentials, respectively. Solid
lines with stars are for GRIP and translucent lines are for the
$\gamma$-surface method. The other data (gray in (a) and green/purple in (b)
and (c)) are referenced from the literature [30, 34, 35, 32]. Select TBs
corresponding to energy cusps are labeled in panels (b) and (c).

To better characterize where our algorithm improves upon existing studies, we
compare plots of $E_{\mathrm{gb}}$ vs. the tilt angle ($\theta$) for all
ground-state structures in Figure S10. Using the EAM potential from Zope and
Mishin [40] (solid orange), we match or improve upon the results from Bhatia
and Solanki [34] (green) and Wang and Beyerlein [35, 32] (purple) for the
$[1\bar{1}00]$ and $[1\bar{2}10]$ STGBs, as seen in panels (b) and (c). Our
results for $[0001]$ STGBs are in good agreement with those from Zheng, et al.
[30] (gray), although more precise comparisons are not possible as they used a
different EAM potential [69]. As noted in the main manuscript, previous
studies may have inconsistently sampled different terminations [33] or deleted
overlapping atoms [34], so we also perform standard $\gamma$-surface
calculations with perfect bulk slabs for each STGB and plot the results in
corresponding translucent colors in Figure S10. For both the EAM (orange) and
MEAM (blue) potentials, the GRIP data (solid lines) are lower bounds for the
$\gamma$-surface values (translucent), which is consistent with expectations.
We attribute the large differences in $\gamma$-surface sampling, i.e., the
sharp peaks in panels (b) and (c), to our definition of planar terminations
that may not have been similarly enforced in previous works [35, 32, 34, 33].
The majority of remaining discrepancies between the GRIP and $\gamma$-surface
data occur for the family of $[0001]$ STGBs, with the largest difference at
$\theta\approx$16.1\text{\,}\mathrm{\SIUnitSymbolDegree}$$, which is the
$\\{31\bar{4}0\\}$ STGB shown in Figure 3 in the main manuscript. We find for
the other two tilt axes closer agreement for $E_{\mathrm{gb}}$, even when many
STGBs require GCO, as seen in Figure 4 in the main text, suggesting that
structural differences in these STGBs are small (see Figure S6). Consistent
with the only existing study that used an evolutionary algorithm to study a
few GBs in an HCP metal (Mg) [24], we also observe a zigzag distribution of GB
dislocations instead of a straight line in several $[1\bar{1}00]$ and
$[1\bar{2}10]$ STGBs.

Figure S11: Faceting of the boundary in $[1\bar{1}00]$ STGBs. The
$\\{11\bar{2}4\\}$$[1\bar{1}00]$
($\theta\approx$38.5\text{\,}\mathrm{\SIUnitSymbolDegree}$$) TB when simulated
with the MEAM potential adopts a thick interfacial structure that is a
strained version of a metastable bulk polymorph [28, 70]. At nearby
misorientation angles, the interfacial phase is partly preserved and the
boundaries with the surrounding $\alpha$-Ti slabs are faceted.

Different empirical potential formalisms are expected to result in different
GB properties, but previous benchmark studies on STGBs in cubic metals [71]
and HCP $\alpha$-Zr [72] found largely similar $E_{\mathrm{gb}}$ vs. $\theta$
profiles using the $\gamma$-surface method. In contrast, we find notable
differences for multiple tilt axes in $\alpha$-Ti when comparing the EAM [40]
(orange) and MEAM [38] (blue) parameterizations. In the family of $[0001]$
STGBs, the values for $E_{\mathrm{gb}}$ from the MEAM potential are
consistently lower than those produced by the EAM potential, with a noticeable
energy cusp (local minimum) at
$\theta\approx$16.1\text{\,}\mathrm{\SIUnitSymbolDegree}$$. Moreover, this
cusp is only present for the GCO data and not the $\gamma$-surface
results—which in fact peaks—further demonstrating the advantages of the GRIP
algorithm. Likewise, in the family of $[1\bar{1}00]$ STGBs (Figure S10b), only
using the MEAM potential do we recover a low-energy $\\{11\bar{2}4\\}$ TB at
$\theta\approx$38.5\text{\,}\mathrm{\SIUnitSymbolDegree}$$ that we have
extensively characterized using transmission electron microscopy and DFT [28].
There, we found the TB to adopt a thick body-centered orthorhombic (BCO)
structure that is a strained version of a metastable bulk polymorph of Ti
[70]; however, what this high-throughput study also reveals is faceting around
the interfacial phase at nearby tilt angles to accommodate the strained BCO
phase, as shown in Figure S11. Stabilization of the BCO phase may be
responsible for the significantly lower energy of $[1\bar{1}00]$ STGBs
simulated with the MEAM potential vs. the EAM potential shown in Figure S10b,
and we leave this analysis to future work.