\section{Introduction}
Understanding and controlling microstructural evolution in metals and metallic alloys 
is one of the central tasks of materials science and engineering. The dynamics of grain
growth and the evolution of grain shape in metallic microstructures strongly depends on the individual mobility of different grain boundaries (GBs), i.e., on the change of the interface energy with its structural parameters \cite{salama2020role, nino2023influence, conry2022enginering, bhattacharya2021grain}.

Nowadays, numerical modeling methods for microstructure and microstructure evolution are available, which explicitly include the variation of interface energy with the geometric degrees of freedom of the grain boundary \cite{vakili2020multiphasefield,
steinbach2011phasefield,
moelans2008introduction,
lee2000montecarlo,
pauza2021computer}. To some extent, this variation can be captured by analytical models based on the Read-Shockley-Wolf (RSW) model \cite{wolf1989read}, see e.g. \cite{bulatov2014grain,dette2017efficient, chirayutthanasak2022anisotropic,sarochawikasit2021grain}. The improvement of such models, and even more a purely numerical treatment of energy as a function of geometry,  
rely on comprehensive data bases of grain boundary energies, which can be generated in a systematic fashion via atomistic simulations.
However, being comprehensive is a quite challenging task. On the one hand, the parameter space of GBs is five-dimensional, defined by the misorientation axis and angle, as well as the grain boundary plane inclination. On the other hand, the grain boundary energy does not vary in a monotonous manner, but exhibits deep cusps at certain misorientations or boundary plane inclinations.
Thus, a standard high-throughput sampling of the parameter space on a regular grid has a twofold drawback:
It is both very time consuming and likely to miss the most important features in the energy landscape.
To provide an efficient data base, however, these should be included in the sampling. This either requires a sampling strategy based on prior knowledge or at least reasonable assumptions on the topology of the energy landscape  \cite{olmsted2009survey,bulatov2014grain,baird2021five, homer2022examination}, and sometimes even the manual addition of the relevant data  \cite{kim2011identification}.
Based on a symmetry analysis and prior knowledge, Olmsted et al.~\cite{olmsted2009survey} designed a strategy to create an energy data base starting with several 1D subspaces and extending to higher dimensions on from there. Bulatov \cite{bulatov2014grain} used 
this data  to interpolate in between the sampled regions by an extended RSW model.
Homer et al.~\cite{homer2015grain} focused on 2D inclination subspaces of coincidence-site lattice based grain boundaries and in a first step reduced the size of the subspace of interest as far as possible by exploiting their point symmetries \cite{patala2013representation,patala2013symmetries}.
Randomly chosen structures from the reduced subspaces were then simulated to further explore it. 
Although impressive progress has been made in these publications, even tackling the complete 5D parameter space \cite{homer2022examination}, none of the mentioned approaches solves the problem of how to find the cusps in the energy landscape automatically.

It is tempting to replace the necessary a priori knowledge by the use of modern machine learning methods, which have become more and more popular and effective in material science \cite{butler2014machine}. Zhang et al.
\cite{zhang2018strategy}, for example, studied how machine learning can be applied accurately to sparse datasets.
As an example they studied the prediction of the band gap of binary semiconductors.
First machine learning approaches for grain boundary energies have been proposed e.g.~by Restrepo et al.~\cite{restrepo2014artificial}, who successfully trained an artifical neural network to predict GB energies by training it with the data collected in \cite{kim2011identification}.
However, also here information concerning position and energy of the cusps was already part of the training data.
\emph{Active learning} \cite{settles2012active} provides a promising remedy for this drawback.
In contrast to traditional design of experiment approaches, where the sampling design is fixed beforehand, active learning starts with a comparatively small dataset and then successively proposes where to further explore the parameter space, based on the analysis of the existing data, until a learning goal is reached. 
Such a sequential procedure hopefully results in a better detection of important features than mere high-throughput sampling with a regular sampling design.

Recently,  Kroll et al.~\cite{kroll2022efficient}  have proposed a method along these lines.
It combines a statistical sampling of the parameter space via a sequential design of experiment approach with a Kriging interpolator to estimate the energy function. Using the jackknife variance, the choice of the next point in the sequential design is a compromise between sampling the region of largest  fluctuations and avoiding a clustering of data points. In this way, the cusps of the energy can be found within a small number of iterations and refined as desired.

To turn this approach into an active learning method, one needs to answer the question: 
when is the sequential sampling good enough to stop the atomistic simulations? The answer sounds simple -- when all relevant cusps of the grain boundary energy  have been found and the predicted energies in unsampled regions are accurate enough. However, to express this answer in measurable quantities and implement it in terms of an automated stopping criterion is not equally obvious, for the following reasons:
\begin{itemize}
\item There is no way to calculate the absolute error of the predicted values, since the true energy distribution is unknown.
\item There is no way to determine a priori the number of cusps, even in a low-dimensional subspace of the 5D parameter space.
\end{itemize}
Thus, what is needed is a measure for the convergence of the energy prediction, which is based only on the already calculated data, and a definition of a \emph{sufficient number of cusps} to describe all relevant features of the grain boundary energy variation. In this work, we develop  such a measure, which addresses both aspects and use it to define a stopping criterion for a sequential sampling algorithm. 

In addition to the new  stopping criterion for sequential sampling of grain boundary energies, this work extends the methodology of \cite{kroll2022efficient}  for the 1D subspace of symmetrical tilt grain boundaries (energy as a function of misorientation) to the two-dimensional subspace of energies as function of GB plane inclination. This creates the additional challenges of 
how to choose an initial design for the fundamental zone (FZ) of the 2D subspace, and how to properly interpolate the energy on a suitable grid. The fundamental zone is the minimum area, from which the whole inclination subspace can be constructed by symmetry operations such as rotation and reflections. The concept of fundamental zones itself is described in \citep{homer2015grain}.
The different instances of interpolation in 2D will be explained below.

In the following Section~\ref{sec:meth} the active learning algorithm is explained, starting with a general description of the overall procedure in Section \ref{sec21}. 
The different types of  grids  used in various steps of the  algorithm in 2D
are introduced in Section \ref{sec:grid_generation} and further illustrated in Appendix~\ref{app:grids}.
The stopping criterion itself and the difference between its application to 1D and 2D samplings is elucidated in Section \ref{sec:stopping-criterion}. Roughly speaking it monitors two quantities: the development of the energy profile and the number of cusps.
The results part of this paper in Section \ref{sec:results} starts with a  validation of the stopping criterion by post-processing the data of the 1D STGB subspaces from \cite{kroll2022efficient}. Here it is demonstrated  that the criterion makes the sampling more efficient,  which means that  we can achieve the same accuracy with fewer atomistic simulations
(Section~\ref{sec:1D_Results}). In Section \ref{sec:2D_subspaces} 
we demonstrate the advantages of the new algorithm for sampling  2D subspaces of grain boundary plane inclinations. In particular, we investigate in Section~\ref{sec:Delta_E_stat}  the impact of the choice of the energy threshold (i.e., the desired accuracy)  on the quality  of the prediction and the number of necessary steps to reach it.
Section \ref{sec:contribution} analyses how the two  quantities 
monitored by the stopping criterion evolve throughout the sampling
and it is demonstrated that both criteria are indispensable.
Finally, the active learning  procedure requires an initial design and the 
influence of  its size towards the quality of the sampling is discussed in Section~\ref{sec:effect}.

\section{Methodology}\label{sec:meth}
\subsection{Basic steps of the procedure}  \label{sec21}
\begin{figure}
     \centering
     \begin{subfigure}{0.9\textwidth}
         \centering
         \includegraphics[width=\textwidth]{Plots/General_overview_and_grids_V3.pdf}
     \end{subfigure}
    \caption{Flowchart of the overall
    procedure for the sampling of grain boundary energy subspaces. The method can be divided into three parts: initial design, sequential design and final Kriging interpolation. The stopping criterion (grey box) represents the new part of the algorithm compared to the method proposed in \cite{kroll2022efficient}. The second addition is the extension to 2D energy subspaces of grain boundary inclinations. The example shown for the final interpolation is the fundamental zone of such a subspace for [110]7.5$^{\circ}$ GBs in fcc nickel.} 
    \label{fig:Flowchart}
\end{figure}
The overall sampling approach is schematically shown in Figure \ref{fig:Flowchart}. It consists of three parts, the \emph{initial design}, the \emph{sequential design} step including 
the stopping criterion, and the \emph{final prediction}.
In a nutshell, the building blocks of the overall method can be summarized as follows:
\begin{itemize}
    \item[(1)]\label{it:summary:1} The \emph{initial design} 
    defines  ${N_{\mathrm{init}}}$ points 
    in the fundamental zone, for which the corresponding GB energies are calculated from molecular statics, as explained in Appendix~\ref{sec:atomistics}. These data points can be obtained in a regular high-throughput scheme. 
    \item[(2a)]\label{it:summary:2} The \emph{sequential design} consists of an ongoing sequence of data generation  and data evaluation. Starting with the initial design, the energy distribution in the fundamental zone is predicted by a  Kriging estimator (see Appendix~\ref{app:kriging} for details). Then, in each following iteration the next sampling point is chosen from ${N_{\mathrm{cand}}}$ candidate points via a jacknife method (basically a cross-validation method).
 Here we compute 
    for any  candidate point the  Kriging prediction  from  the  currently available observations and  from  the reduced sets obtained by deleting exactly one of the current observations at a time.
    The Kriging estimators obtained this way are combined in order to define a jackknife estimate of the uncertainty of the prediction at any of the ${N_{\mathrm{cand}}}$ candidate locations, and the next sampling point is then chosen among the maximizers of this estimate.
    Finally, an atomistic simulation is conducted at this new point,
    and an updated Kriging model is fitted to the augmented  data set (previous design + new point).
    We refer to \cite{kroll2022efficient} for more details.
    
     \item[(2b)]
    Next, the validity of the new {\it stopping criterion},  which will be explained in detail in Section \ref{sec:stopping-criterion} below, is checked. If the criterion is satisfied, sequential sampling is terminated, otherwise continued. 
   This means, that, in contrast to \cite{kroll2022efficient}, the number of atomistic simulations, say ${N_{\mathrm{seq}}}$, is not fixed  a priori, but determined by the stopping criterion based on the results of the  conducted experiments.
    \item[(3)]\label{it:summary:3} After sequential sampling has stopped, 
    \emph{final prediction} on a dense grid of points is performed via Kriging on the basis of the complete dataset (initial design + sequential design). 
 \end{itemize}

In Appendix~\ref{app:kriging}, we give a brief summary of  Kriging prediction, which is used in steps (2a) and (2b) of the overall procedure
and will be used as a black box method 
in the remaining part of the paper.  We also refer to \cite{rasmussen2006gaussian}, Section \ref{sec:grid_generation}, or \cite{stein1999interpolation} for further details.

The overall procedure extends the one from \cite{kroll2022efficient} mainly in two directions.
First, the methodology from \cite{kroll2022efficient} is extended also to the 2D case which requires the definition of suitable sampling grids beyond the 1D case where the definition of regular grids is obvious.
Second, and maybe more important, in step (2b) we add the stopping criterion to the overall procedure which is ought to release us from the task of specifying the number of sequential steps in advance.
These two ingredients of the overall procedure will be described explicitly in the following sections.
 

\subsection{Grid generation}
\label{sec:grid_generation}
All three steps outlined in Section~\ref{sec21} involve Kriging interpolation of the available data which is carried out on a pre-defined grid.
For step~(1) of the algorithm, the initial design, a space-filling $s^2$-grid is used.
This grid is especially suited to achieve a homogeneous distribution of the sampling points and can be motivated from the use of Kriging as the interpolation method of our choice,
see Appendix~\ref{sec:s2-grid} for a more detailed discussion.
Such an $s^2$-grid with a larger sample size was also used to define the reference designs for the evaluation of the performance of different predictions. Of course, in a real application these reference designs are not available.
Since the construction of a large $s^2$-grid is rather time consuming, a regular equally distant angular grid (see Appendix~\ref{sec:angular-grid}) is used for the Kriging interpolation in step~(2a) and (2b) of the algorithm, where it is required to evaluate the stopping criterion.
In step~(3), the final interpolation is carried out on a very fine version of this regular equally distant angular grid.
Note that, as the number of points increases, the advantages of the $s^2$-grid become less pronounced, and for a grid with a very high density the regular equally distant angular grids perform nearly in the same way.

Also the candidate points for the jackknife method, which is applied to find the next point during the sequential design step, are supposed to lie on a pre-defined space-filling grid. Here, the space-filling property is desirable to make measurements in all areas of the FZ at least potentially possible. For this purpose, we have taken, mainly for computational reasons, the reduced angular grid (see Appendix~\ref{sec:reduced-angular-grid}).
In \citep{kroll2022efficient} some heuristics concerning the choice of ${N_{\mathrm{cand}}}$ were briefly discussed, but only a fixed value of ${N_{\mathrm{cand}}} = 75$ was considered in the simulations.
In this work, after a more detailed investigation of its impact (see Appendix~\ref{app:Ncand}), ${N_{\mathrm{cand}}}$ is increased to 200.

\subsection{Stopping criterion}
\label{sec:stopping-criterion}
\begin{figure}
\centering
    \includegraphics[width=0.57\textwidth]{Plots/stopping_criterion_flowchart.pdf}
     \caption{Flowchart, which represents the code behind the stopping criterion. After performing Kriging, the interpolated dataset is analysed for minima. Sharp minima are processed directly, while for extended minima, so-called valleys (see text), first the central point is determined. After deleting minima, which are too close, the topological stopping criterion (expression A) can be checked. For the statistical one, the energy difference between the datapoints in the current and previous step is analysed as indicated schematically in the lower panel and enters expression B.}
\label{fig:matrix}
\end{figure}
In this section we develop a stopping criterion which combines a topological aspect with a statistical one. 
Before it is applied, the available data is interpolated by performing Kriging on a very fine grid with equally distant points. In the 1D subspace of symmetrical tilt grain boundaries, this is a 1D grid of equally distant misorientation angles. In the 2D space that defines the boundary plane normal it is a 2D grid of equally distant polar and azimuthal angles (the regular equally distant angular grid, see Appendix~\ref{sec:angular-grid}.)
To check the topological aspect of the stopping criterion, the cusps in the energy landscape are identified by comparing the energies of neighbouring points. This requires a list of neighbours (defined by the misorientation angle) for the 1D case and a matrix of neighbours (defined by the azimuthal angle in the column and the polar angle in the row) as indicated in Figure \ref{fig:matrix}.

In 1D subspaces, each element of the list is analysed by checking the previous and the following element (2-connectivity).
If both neighbouring elements have a higher energy then the analysed element is classified as a minimum (cusp). If several neighbouring elements have the same energy and are surrounded by elements with a higher energy, they form what is called a valley. To define the actual position of the energy minimum, the centre of the valley is calculated, i.e.~the position of the cusp is chosen as the mean value of the misorientation angles of the elements in the valley. 

To identify energy minima in the 2D energy subspaces, the difference between the energy of each element  to  each of its eight neighbours (8-connectivity) is analysed. 
Again, extended regions of low energy, the above mentioned valleys, can occur and their centres have to be identified.  For this purpose a recursive method is applied, which is also illustrated  in Figure \ref{fig:matrix} (calculate center
of valleys). The middle points between each element of a valley and its nearest neighbours (if still part of the valley) are added to a new list of reduced number of points. This will be repeated with every element in the new list to create a further reduced list. The process repeats until the list only contains one element, which will be the centre of the valley. 

After identification of all minima, the minima which are closer to each other than 
2$^{\circ}$ for 1D and 2D are  reduced to the minimum with the lowest energy.
Here, for 2D subspaces the angular distance is calculated as
\begin{equation}
    \Delta \alpha = \sqrt{\Delta \varphi^2 + \Delta \vartheta^2}.
\end{equation}

Finally, the stopping criterion is checked, i.e., whether the number of cusps is constant with respect to the previous ${N_{\mathrm{iter, cusps}}}$ steps. 

The statistical aspect of the stopping criterion is defined similarly: The change in energy at each point compared to the previous sequential step is calculated. If the maximum difference of the energy towards the previous iterations in the whole subspace is lower than a threshold $\Delta E_{\mathrm{stat}}$ for the next $N_{\mathrm{iter}, \Delta E}$ iterations, also the statistical aspect of the stopping criterion is fulfilled.
If both subcriteria are met, the criterion itself is fulfilled and the sampling stops. A Voronoi tesselation can be applied to the cusps to divide the subspaces into cells. Strictly speaking only the maximum difference in the overall subspace is needed to evaluate energetic aspect, but the calculation of the maximum difference in each Voronoi cell provides a closer look on what is happening in the cells while sampling the subspace.

To summarize, both the overall energy as well as the number of cusps are monitored. Convergence is reached once both, the energy and the number of cusps are stable for a certain number of iterations.  The relevant parameters are listed in Table~\ref{table}.
\begin{table}\footnotesize
	\begin{tabular}{ll}
    \multicolumn{2}{l}{\bfseries{Pre-defined parameters}}\\
    \\
	${N_{\mathrm{init}}}$ & number of initial design points\\
	${N_{\mathrm{cand}}}$ & number of candidate points from which the next point is chosen in each\\
    & sequential step\\	
  $\Delta E_{\mathrm{stat}}$ & 
  threshold to which $\Delta E_{\mathrm{prev}}$ is compared as part of the stopping criterion\\
    $N_{\mathrm{iter}, \Delta E}$ & number of iterations during which $\Delta E_{\mathrm{prev}}$ must stay below the threshold 
  \\ & $\Delta E_{\mathrm{stat}}$ (statistical stopping criterion)\\
	${N_{\mathrm{iter, cusps}}}$ & number of iterations for which the number of division must not change\\
  &   (topological stopping criterion)\\
    \\
\multicolumn{2}{l}{\bfseries{Runtime variables}}\\
    \\
	${N_{\mathrm{seq}}}$ & current number of sequential step\\
    ${N_{\mathrm{cusps}}}$ & number of divisions/cusps in the FZ\\
	$\Delta E_{\mathrm{prev}}$ & maximum absolute deviation between Kriging interpolators from two\\
    & consecutive steps\\
  		$\Delta E_{\mathrm{ref}}$ & in current step, maximum absolute deviation of the Kriging interpolator\\
    & from a reference data set\\
  ${N_{\mathrm{stop}}}$ &\makecell[l] {number of iterations required to fullfil the stopping criterion}
\end{tabular}
 \caption{Overview of code-related abbreviations used in the text.}
\label{table}
\end{table}

In this work, $N_{\mathrm{iter}, \Delta E}$ and ${N_{\mathrm{iter, cusps}}}$ are chosen as 3.
Moreover, as mentioned already above, a fixed value of ${N_{\mathrm{cand}}}=200$ was considered.

\section{Results}\label{sec:results}

\subsection{$1D$ STGB subspaces}
\label{sec:1D_Results}

Recently it was demonstrated in \cite{kroll2022efficient}  that a sequential design is able to identify the cusps in the 1D energy landscape of symmetrical tilt grain boundaries in body-centered cubic iron. 
However, in the cited paper the optimal number of sequential steps was determined \emph{a posteriori} by analysing the maximum error with respect to a reference data base.
In the following, the same data is used to validate the new stopping criterion.

The topological aspect of the stopping criterion is illustrated in Figure~\ref{fig:division_boundaries_1D}
for the 
$[100]$ and $[110]$ STGB subspaces,
which also shows the location of the design points of the sequential sampling.
On the bottom line, the  locations of the initial design points are displayed.
Going up along the $y$-axis, the evolution of the positions of the cusps 
(vertical dashed lines) and the location of the sequentially chosen design points ($\blacktriangle$) and initial chosen design points ($\blacktriangle$ located at the $x$-axis) can be tracked. For example, for 
the [100]-subspace the algorithm chooses 
the  $88.5^\circ$ misorientation angle 
in the first iteration, $0.75^\circ$ in the second,  $85.15^\circ$, in the third and so on.
It can be seen that the sequential algorithm allocates a large number of design points in neighbourhoods of the (unknown) cusps. Moreover, it does not select new sampling locations from the same region over several iterations but rather visits neighbourhoods of other cusps. In addition, the plots illustrate that the number of cusps increases with the number of sequential steps performed and finally converges. 
For example, for the [110] subspace the algorithm starts with $5$ cusps, and after $8$ and $10$ iterations it detects $6$ and $7$ cusps, respectively.

\begin{figure}[htb]
     \centering
    \begin{subfigure}{0.49\textwidth}
         \centering
    \includegraphics[width=\textwidth]{Plots/100_cusp_positions_plot.pdf}     
     \end{subfigure}
     \begin{subfigure}{0.49\textwidth}
         \centering
          \includegraphics[width=\textwidth]{Plots/110_cusp_positions_plot.pdf} 
     \end{subfigure}
     \caption{Evolution of the sequential design for the $[100]$ 
     (left panel) and $[110]$ STGB subspace (right panel). 
     $x$-axis: misorientation angle; left $y$-axis:  current sequential step; right $y$-axis:   energy. The solid blue line is the grain boundary energy as a function of misorientation after the final sampling step. The \textcolor{black}{$\blacktriangle$} on the $x$-axis display the initial design, the other \textcolor{black}{$\blacktriangle$}
     indicate the positions of the new design points calculated by 
     the sequential algorithm. The vertical dashed lines mark the positions of the cusps (they start at the sequential step where a new cusp was discovered).}
\label{fig:division_boundaries_1D}
\end{figure}

\begin{figure}[htb]
     \centering
    \begin{subfigure}{0.49\textwidth}
         \centering
         \label{fig:1D:max:error:100}
    \includegraphics[width=\textwidth]{Plots/100_max_error_ref_prev_with_div_numb_comparision.pdf}
     \end{subfigure}
     \begin{subfigure}{0.49\textwidth}
         \centering
          \includegraphics[width=\textwidth]{Plots/110_max_error_ref_prev_with_div_numb_comparision.pdf}   
     \end{subfigure}
     \caption{
     Evolution of the two contributions to the stopping criterion: maximum absolute error with respect to  the previous sequential step
     (left $y$-axis, \textcolor{red}{\textbullet})
and the number of cusps, ${N_{\mathrm{cusps}}}$ (right $y$-axis, \textcolor{blue}{$\blacksquare$}).
For the sake of comparison the maximum absolute error  with respect to  a reference data base (left $y$-axis, $\blacktriangle$) 
     is also displayed.
     Left panel: $[100]$ subspace with ${N_{\mathrm{init}}} = 16$; Right panel: $[110]$ subspace with ${N_{\mathrm{init}}} = 31$.
}
\label{fig:STGB_max_error_comparison}
\end{figure}
In \cite{kroll2022efficient}, the quality of sampling was evaluated by the maximum error of the Kriging interpolator with respect to a reference data set. This was done for development purposes - in a practical application, i.e.,
when sampling a completely unexplored subspace, such a reference data set does of course not exist.
Thus, in this work, the maximum error with respect to the previous iteration, $\Delta E_{\mathrm{prev}}$, is introduced as an alternative error measure, which can be computed from the observed data only and is monitored by the stopping criterion.

Figure~\ref{fig:STGB_max_error_comparison} shows the development of 
this error measure  during the  sequential sampling  (red circles; the size of the error can be determined from the left $y$-axis). Similarly, we display the development of the number of cusps/minima (blue squares; the number of cusps can be determined from right $y$-axis). Moreover, the error with respect to the reference data base is represented by the black triangles for the sake of a qualitative comparison.
For example, for the $[100]$-subspace we observe from the left part of Figure~\ref{fig:STGB_max_error_comparison}  
 that after $5$ iterations
$\Delta E_{\mathrm{prev}} \approx 0.3 \si{J}/\si{m}^2$ and $5$ cusps have been detected.
 
The comparison of the evolution of $\Delta E_{\mathrm{prev}}$ and ${N_{\mathrm{cusps}}}$ with $\Delta E_{\mathrm{ref}}$ shows the strength of the new criterion.  While the latter keeps decreasing, $\Delta E_{\mathrm{prev}}$ increases again when a new cusp has been found, and even if ${N_{\mathrm{cusps}}}$ remains constant, the sequential design continues until the desired threshold for  $\Delta E_{\mathrm{prev}}$ is reached. For instance, 
in the [110] and the [111] subspaces $\Delta E_{\mathrm{ref}}$ remains constant for several iterations (after $6$ iterations for [110] and $2$ iterations for [111]), while $\Delta E_{\mathrm{prev}}$ still varies.
At the same time, Figure~\ref{fig:STGB_max_error_comparison} also illustrates the importance of the topological part of the stopping criterion. After $8$ respectively $10$ iterations in the [110] subspace the number of cusps increases, accompanied by an increase in $\Delta E_{\mathrm{prev}}$, but $\Delta E_{\mathrm{ref}}$ does not change. Nevertheless, monitoring only $\Delta E_{\mathrm{prev}}$ is not sufficient. To see this, consider the [100] subspace
and note that after $9$ iterations $\Delta E_{\mathrm{prev}}$ is lower than the threshold $\Delta E_{\mathrm{stat}} = 150 \si{mJ}/\si{m}^2$ for several iterations, but the number of cusps still increases, which means that the stopping criterion is not yet fulfilled here.
\begin{figure}[htb]
     \centering
    \includegraphics[width=\textwidth]{Plots/STGBs_N_stop_3_3.pdf}     
     \caption{
   Number  of sequential steps (${N_{\mathrm{stop}}}$) required until the algorithm terminates for different values of $\Delta E_{\mathrm{stat}}$ that specify the desired accuracy (measured by $\Delta E_{\mathrm{prev}}$ in each iteration).
    $\blacktriangle$ mark the $[100]$ subspace with ${N_{\mathrm{init}}} = 16$, \textcolor{red}{\textbullet} the $[110]$ subspace with ${N_{\mathrm{init}}} = 31$ and \textcolor{blue}{$\blacksquare$} the $[111]$ subspace with ${N_{\mathrm{init}}} = 21$.}
\label{fig:STGB_N_stop_3_3}
\end{figure}
In Figure~\ref{fig:STGB_N_stop_3_3} the impact of the required accuracy $\Delta E_{\mathrm{stat}}$ on the stopping criterion is studied for three 1D subspaces  (STGB subspaces with a fixed rotation axis of $[100]$, $[110]$ and $[111]$). 
More precisely, for various values of the input parameter $\Delta E_{\mathrm{stat}}$ the figure displays the number of sequential steps, denoted by ${N_{\mathrm{stop}}}$, which is required until termination of the algorithm.
Note that the algorithm eventually stops for any choice of $\Delta E_{\mathrm{stat}}$ and for any subspace. For example, for a statistical accuracy $\Delta E_{\mathrm{stat}}= 100 \si{mJ}/\si{m}^2$ the algorithm stops sequential sampling after $14$, $8$, and $6$ iterations for the $[100]$, $[111]$, and $[110]$ subspace, respectively.
Clearly, the required number of iterations is a decreasing function of the  desired accuracy $\Delta E_{\mathrm{stat}}$.

The benefit of the stopping criterion becomes clear when the point of termination and resulting accuracy is compared to the empirically chosen number of sequential iterations in \cite{kroll2022efficient}. In that work, the number of iterations was set to $20$, which 
corresponds to a sampling where  $\Delta E_{\mathrm{stat}}$ is not more than $25 \si{mJ}/\si{m}^2$ for the [110] and [111] subspaces  and not more than  $10 \si{mJ}/\si{m}^2$ for the [100] subspace.
On the basis of the new stopping criterion the algorithm terminates 
sampling much earlier and still achieves the same precision.
In this regard, the stopping criterion is not only a tool to automatise, but also to optimise the sampling procedure.

\subsection{$2D$ inclination subspaces}
\label{sec:2D_subspaces}
In this section the algorithm with the new stopping criterion is applied to  $2$D inclination subspaces.
Different samplings of the inclination space of the $\Sigma 3[111]60^\circ$ grain boundaries in bcc Fe, as well as the of the $\Sigma 5[100]36.87^\circ$ and $\Sigma 7[111]38.21^\circ$ boundaries in fcc Al and $[110]7.5^\circ$ boundaries in fcc Ni are  considered.

As for the 1D STGB subspaces, a reference data base was generated for each subspace to be able to judge the quality of the the sequential sampling and the stopping criterion. 
The results of the Kriging interpolation for this data base are displayed  in the left part of  Figure \ref{fig:All_Kriging_plots}, which shows the energy distribution in the fundamental zone.
The complete space of grain boundary inclination represents the surface of a sphere, with the normal vector of the GB plane pointing to a point on this sphere. 
The individual figures show how the complexity and the amount of substructure of the energy plot increases with the $\Sigma$ value, i.e.~with decreasing symmetry (the $\Sigma$ value for the non-periodic $[110] 7.5^\circ$ small angle grain boundaries is infinite).  The size of the fundamental zone, but also the density of cusps increases from the $\Sigma 3$ to the $[110] 7.5^\circ$ boundaries. This makes a sampling very challenging, calling once more for a reliable stopping criterion.

In the right column of the figure we show the prediction based on sequential sampling. We also display  the Voronoi cells around the cusps 
to illustrate the increasing 
complexity of fundamental  zones.
The results demonstrate the advantages of our approach. We obtain qualitatively very similar energy plots, but, compared to the reference design, the number of atomistic simulations to obtain these plots is substantially reduced.  For example, for the left panel in part (d)  of Figure \ref{fig:All_Kriging_plots} a reference design with $100$ points is used, while the right panel is obtained from a sequential design with $54$ points, which yields a substantial reduction of computational costs.
In the following, the two contributions to the new criterion are evaluated for the 2D cases. Furthermore the quality of the energy prediction is also governed by two aspects which will be discussed in more detail: the threshold value $\Delta E_{\mathrm{stat}}$ which defines the convergence of the energy
and the number of initial design points ${N_{\mathrm{init}}}$.\\

\begin{figure}
     \centering
    \includegraphics[width=0.7\textwidth]{Plots/All_Kriging_plots.pdf}     
     \caption{ Left panels: 
     predicted energies  on the basis of a 
     reference data bases.
     Right  panels: predicted energies from the sequential sampling. 
     $\Delta E_{\mathrm{stat}}$ equals $50 \si{mJ}/\si{m}^2$ in each sampling and the 
     black lines indicate the boundaries of 
     Voronoi cells.\\
     (a): $\Sigma 3$ subspace; ${N_{\mathrm{ref}}} = 150$, ${N_{\mathrm{init}}} = 25$ and ${N_{\mathrm{seq}}} = 13$. (b): $\Sigma 5$  subspace;  ${N_{\mathrm{ref}}} = 100$, ${N_{\mathrm{init}}} = 20$ and ${N_{\mathrm{seq}}} = 6$. (c): $\Sigma 7$ subspace;  ${N_{\mathrm{ref}}} = 100$, ${N_{\mathrm{init}}} = 28$ and ${N_{\mathrm{seq}}} = 3$. (d): $[110] \, 7.5^\circ$ subspace; ${N_{\mathrm{ref}}} = 100$, ${N_{\mathrm{init}}} = 40$ and ${N_{\mathrm{seq}}} = 14$.}
\label{fig:All_Kriging_plots}
\end{figure}

\begin{figure}
     \centering
         \centering
    \includegraphics[width=\textwidth]{Plots/all_subspaces_N_stat_iter_3_N_geo_iter_3.pdf}
     \caption{
     Number ${N_{\mathrm{stop}}}$ of sequential steps when the algorithm is terminating.
     $\Sigma 3$ subspace (${N_{\mathrm{init}}} = 50$): $\blacktriangle$;   $\Sigma 5$ subspace (${N_{\mathrm{init}}} = 20$): \textcolor{red}{\textbullet};
     $\Sigma 7$ subspace with ${N_{\mathrm{init}}} = 28$: \textcolor{blue}{$\blacksquare$}; $[110]7.5^\circ$ subspace (${N_{\mathrm{init}}} = 40$): \textcolor{SpringGreen4}{$\star$}.
     }
\label{fig:2D_iend}
\end{figure}
\subsection{$\Delta E_{\mathrm{stat}}$ and the speed of convergence}
\label{sec:Delta_E_stat}
The algorithm with the new  stopping criterion is applied to the data of the 2D inclination subspaces with different values for $\Delta E_{\mathrm{stat}}$.
In Figure \ref{fig:2D_iend}, the number of sequential steps, ${N_{\mathrm{stop}}}$, when the algorithm is terminating (such that the stopping criterion is satisfied) is displayed
as a function  of $\Delta E_{\mathrm{stat}}$.  Similar to the 1D subspace, ${N_{\mathrm{stop}}}$ decreases 
with an increasing $\Delta E_{\mathrm{stat}}$, and, if two neighbouring $\Delta E_{\mathrm{stat}}$ yield the same value  ${N_{\mathrm{stop}}}$, the lower value of $\Delta E_{\mathrm{stat}}$ marks the actual accuracy. In other words, a lower $\Delta E_{\mathrm{stat}}$ does not lead to a gain in speed, because the fluctuations in the energy from one step to the next are small, anyhow. 
The values of  ${N_{\mathrm{stop}}}$ for the $\Sigma 3$, $\Sigma 5$ and $\Sigma 7$ subspaces differ only  by $3$, but  for  the $[110]7.5^\circ$ ${N_{\mathrm{stop}}}$ is  significantly larger. This is an effect of the complex energy landscape of  the $[110]7.5^\circ$  subspace rather than of the size of the fundamental zone, which becomes apparent when the two contributions to the stopping criterion are analysed.

\subsection{The impact of $\Delta E_{\mathrm{prev}}$ and ${N_{\mathrm{cusps}}}$ on  the stopping criterion}
\label{sec:contribution}
\begin{figure}
     \centering
    \begin{subfigure}{0.49\textwidth}
         \centering
         \caption{}
    \includegraphics[width=\textwidth]{Plots/Sig3_111_60_50_25_nu_MLE_Ncand_200_max_error_ref_prev_comparision_plus_div_number_equally_distant_grid.pdf}     
     \end{subfigure}
     \begin{subfigure}{0.49\textwidth}
         \centering
         \caption{}
          \includegraphics[width=\textwidth]{Plots/Sig5_100_36_87_20_10_nu_0_5_Ncand_200_max_error_ref_prev_comparision_plus_div_number_equally_distant_grid.pdf}   
     \end{subfigure}
         \begin{subfigure}{0.49\textwidth}
         \centering
         \caption{}
    \includegraphics[width=\textwidth]{Plots/Sig7_111_38_21_28_14_nu_0_5_Ncand_200_max_error_ref_prev_comparision_plus_div_number_equally_distant_grid.pdf}     
     \end{subfigure}
     \begin{subfigure}{0.49\textwidth}
         \centering
         \caption{}
          \includegraphics[width=\textwidth]{Plots/LAGB_110_7_5_40_20_max_error_ref_prev_comparision_plus_div_number_equally_distant_grid.pdf}   
     \end{subfigure}
     \caption{Maximum absolute error 
     $\Delta E_{\mathrm{ref}}$ with respect to a reference data base,  (left $y$-axis, $\blacktriangle$); maximum absolute error $\Delta E_{\mathrm{prev}}$ with respect to the previous sequential step (left $y$-axis, \textcolor{red}{\textbullet}); number of cusps, ${N_{\mathrm{cusps}}}$ (right $y$-axis, \textcolor{blue}{$\blacksquare$}).\\
     (a) $\Sigma 3$ subspace with ${N_{\mathrm{init}}} = 50$, (b) $\Sigma 5$ subspace with ${N_{\mathrm{init}}} = 20$, (c) $\Sigma 7$ subspace with ${N_{\mathrm{init}}} = 28$ and (d) $[110]7.5^\circ$ subspace with ${N_{\mathrm{init}}} = 40$.}    
\label{fig:2D_error}
\end{figure}

In Section \ref{sec:1D_Results} it has been demonstrated for 1D subspaces that $\Delta E_{\mathrm{prev}}$ is a reasonable criterion to control the sequential sampling procedure. By monitoring both $\Delta E_{\mathrm{prev}}$ and ${N_{\mathrm{cusps}}}$ simultaneously,
a further improvement has been shown. A corresponding comparison in Figure~\ref{fig:2D_error} confirms these findings for 2D subspaces.
It shows that $\Delta E_{\mathrm{prev}}$ contains more information about the state of the sampling than $\Delta E_{\mathrm{ref}}$ and the role of ${N_{\mathrm{cusps}}}$ becomes even more important now.
This is particularly visible in the $[110]7.5^\circ$ subspace with the highest amount of substructure in the energy landscape. 
Between the first and the 9th iteration ${N_{\mathrm{cusps}}}$ changes several times, while $\Delta E_{\mathrm{prev}}$ is already smaller than $0.11 \si{J}/\si{m}^2$ between the second and the 9th iteration. Therefore, without controlling ${N_{\mathrm{cusps}}}$  the algorithm  would terminate even though several cusps are still to be discovered in the next iterations. On the other hand ${N_{\mathrm{cusps}}}$ does not change after $9$ iterations, while $\Delta E_{\mathrm{prev}}$ goes up to $0.15 \si{J}/\si{m}^2$ after 11 iterations. This shows that both aspects of the stopping criterion are necessary to optimize the automatized active learning procedure.

\subsection{The effect of the size of the initial design}
\label{sec:effect}
An aspect which has not been discussed so far, is the influence of the size of the initial design on the quality of the sampling.
\begin{figure}
    \centering
    \begin{subfigure}{0.46\textwidth}
       \caption{}
        \label{fig:Ncusps_allDeltaEstat}
    \includegraphics[width=0.95\textwidth]{Plots/N_cusp_in_dependence_of_Delta_E_stat.pdf}
     \end{subfigure}
     \begin{subfigure}{0.46\textwidth}
    \caption{}
     \label{fig:maxDeltaEref_allDeltaEstat}
          \includegraphics[width=0.95\textwidth]{Plots/max_Delta_E_ref_in_dependence_of_Delta_E_stat.pdf}   
     \end{subfigure}
    \begin{subfigure}{0.46\textwidth}
        \caption{}\label{fig:Ncusps_DeltaEstat0_05}
   \includegraphics[width=0.95\textwidth]{Plots/N_cusps_bar_diagram_Delta_E_stat_0_05.pdf}
    \end{subfigure}
    \begin{subfigure}{0.46\textwidth}
         \caption{}\label{fig:maxDeltaEref_DeltaEstat0_05}
         \includegraphics[width=0.95\textwidth]{Plots/max_Delta_E_ref_diagram_Delta_E_stat_0_05.pdf}   
    \end{subfigure}
     \caption{
    {Results of sequential and regular sampling (with the same total number of sampling points) for 2D spaces: (a) number 
     of cusps 
     detected by the sequential algorithm
     as a function of $\Delta E_{\mathrm{stat}}$; (b) 
     maximum absolute error with respect to the reference data base
     of the sequential algorithm with ${N_{\mathrm{init}}}$ initial design points 
     as a function of $\Delta E_{\mathrm{stat}}$
     compared to a regular sampling with the same total number of sampling points;
     (c) number
     of cusps for different initial designs
     for different subspaces
  ($\Delta E_{\mathrm{stat}} = 50\si{mJ}/\si{m}^2$);
 (d) 
  maximum absolute error with respect to the reference data base 
  for different initial designs for different subspaces
($\Delta E_{\mathrm{stat}} = 50\si{mJ}/\si{m}^2$).}}  
\end{figure}
Figure~\ref{fig:Ncusps_allDeltaEstat} shows the number of cusps which have been found at the end of the sampling as a function of the threshold value $\Delta E_{\mathrm{stat}}$ for different sizes of the initial design.
We observe  that consistently more cusps are found for a larger initial design, with the exception of the $\Sigma 3$ subspace. This exception is due to the rather smooth and flat energy distribution in the fundamental zone of the $\Sigma 3$ grain boundaries. Generally, the number of identified cusps increases with the size and the complexity of the subspace. 
In  Figure \ref{fig:maxDeltaEref_allDeltaEstat} we compare  the 
maximum error $\Delta E_{\mathrm{ref}}$ with respect to the reference data base for a sequential design and a regular design with same number of atomistic simulations.
A remarkable observation   
   is that sequential sampling   yields a smaller  maximum error 
   than the regular high-throughput sampling counterpart. This is particularly visible for the most complex subspace of the $[110]7.5^\circ$ boundaries.  
   These findings show that the algorithm is most efficient for sampling complex subspaces with a large  number of cusps.
While the statements above are largely independent of $\Delta E_{\mathrm{stat}}$, we additionally observe from 
Figures \ref{fig:Ncusps_allDeltaEstat} and \ref{fig:maxDeltaEref_allDeltaEstat} that the sequential design shows increasing improvement (i.e., an increase in ${N_{\mathrm{cusps}}}$ and a decrease in $\Delta E_{\mathrm{ref}}$) for $\Delta E_{\mathrm{stat}}$ smaller than $50\si{mJ}/\si{m}^2$. 

For the value $\Delta E_{\mathrm{stat}} = 50\si{mJ}/\si{m}^2$, Figures \ref{fig:Ncusps_DeltaEstat0_05} and \ref{fig:maxDeltaEref_DeltaEstat0_05} show a detailed comparison of the different samplings with the regular one. 
The numbers above the bars refer to the sequential sampling (blue) and represent the number of initial plus the number of sequential design points, after the stopping criterion has been fulfilled. The total number of sample points is the sum of both, i.e., ${N_{\mathrm{total}}} = 
{N_{\mathrm{init}}} + {N_{\mathrm{seq}}}$. The results for  equivalent regular high-throughput sampling with the same ${N_{\mathrm{total}}}$ as the sequential sampling are shown as red bars.
We observe from the diagrams that energy convergence can also be reached with a small initial design, leading to a small total number of calculations. 
However, the number of determined cusps  is higher and the error with respect to the reference data base is lower for the larger initial designs (with the above mentioned exception for $\Sigma 3$, here the sampling with ${N_{\mathrm{init}}} = 25$ identifies the highest number of cusps.) To choose the optimal ${N_{\mathrm{init}}}$, the size of the different subspaces has to be considered, which goes along with an increase in complexity of the energy function. The maximum polar angle always equals $\vartheta_{\mathrm{FZ}}=90^\circ$. As shown in Figure~\ref{fig:All_Kriging_plots} the maximum azimuthal angle equals $\varphi_{\mathrm{FZ}, \Sigma 3} = 30^\circ$, $\varphi_{\mathrm{FZ}, \Sigma 5} = 45^\circ$, $\varphi_{\mathrm{FZ}, \Sigma 7} = 60^\circ$ and $\varphi_{\mathrm{FZ}, [110]7.5^\circ} = 90^\circ$. The point densities $\rho$ of the inital design can be calculated with the following equation:
\begin{equation*}
    \rho = \frac{{N_{\mathrm{init}}}}{\varphi_{\mathrm{FZ}}}\text{.}
\end{equation*}
In this study we obtained good results for point densities between 
$0.44$ and $0.66$ per degree.

The analysis of Figures \ref{fig:Ncusps_DeltaEstat0_05} and \ref{fig:maxDeltaEref_DeltaEstat0_05} also shows that the sequential design outperforms a regular high-throughput sampling in particular in low-symmetry cases, provided that the threshold for the stopping criterion, $\Delta E_{\mathrm{stat}}$, is chosen small enough. In this work, the variation of grain boundary energies in the fundamental zones angles from $300$ ($\Sigma 7$) to $600 \si{mJ}/\si{m}^2$ ([110] low angle GBs) for the low-symmetry cases, and a good result is obtained for $\Delta E_{\mathrm{stat}} \leq 50\si{mJ}/\si{m}^2$, i.e., for less than 8--16\% of this variation.

\section{Discussion and Conclusion}

This work introduces an algorithm for an automated sampling of grain boundary energies in the spirit of an active learning technique. It is based on the sequential sampling strategy introduced in \cite{kroll2022efficient}, which has been developed for 1D STGB subspaces and was now successfully extended in two directions. On the one hand  the new algorithm can be used for  2D applications. On the other hand, and more importantly,
the new algorithm is able to decide on the basis of the collected data, when it is
reasonable to stop sequential sampling.
The major difference to other methods published so far, e.g.~\cite{kim2011identification} and \cite{restrepo2014artificial}, is that no prior knowledge concerning the location of the cusps is needed. The proposed algorithm rather learns the locations of the cusps automatically and terminates when no new cusps are found or major changes in the energy landscape arise over several iterations. This is a breakthrough in terms of the discovery of cusps in multidimensional subspaces.

Our results  demonstrate that the sequential sampling technique is not only able to outperform a regular high-throughput sampling in 1D subspaces, it also has the possibility to outperform high-throughput samplings of high and low symmetry, coincidence-site lattice based and even and particularly non-periodic low angle grain boundaries 2D inclination subspaces.

To arrive at a practical scheme, two quantities to evaluate the quality of the sampling were defined: the maximum deviation of the energy between two following sequential steps and the number of identified minima in the energy landscape. Both quantities can be calculated on the fly, and in combination they are well-suited to evaluate the sampling based on the available data. This allows to define a stopping criterion for automated sequential runs. In our  applications  it  is demonstrated that the total number of sequential steps depends on the variability of the energy landscape. A larger number of cusps requires more sequential steps. The benefit of monitoring this number  becomes particularly clear when looking at low-angle grain boundaries with a rather volatile energy distribution in the fundamental zone and thus a high frequency in the variation of energy versus angles.

It was also shown that the sampling can be further refined  by a careful choice of the
size of the initial design. 
We observed that a larger number of initial design points patched with a smaller number of sequential design points   outperforms a regular high-throughput sampling. This superiority  is independent of the pre-defined threshold for the energy convergence.

With the help of a reasonable metric, the algorithm presented in this paper can  be extended to even higher dimensional subspaces, or the whole parameter space. To some extent, the required data for such an analysis is already available from \cite{baird2021five} and \cite{homer2022examination}, who examined the topology of the 5D parameter space.

To sum up, the  maximum deviation of the energy between two following sequential steps and the number of identified minima in the energy landscape  are useful measures  to describe the quality of the sampling on the fly while giving advanced information about the subspace itself and its complexity. We have developed a  sequential sampling algorithm with a stopping criterion, which is based on a simultaneous monitoring  of these two quantities.  As a result  we obtain a very efficient active learning procedure for the exploration of grain boundary subspaces. This approach  outperforms high-throughput sampling especially for subspaces of non-periodic grain boundaries after only a small number of iterations, which opens the door to explore the energy landscape of low angle grain boundaries precisely.

\section*{Acknowledgements}
This research has been supported by the German Research Foundation (DFG), project number 414750139.

\printbibliography

\newpage

\pagenumbering{arabic}
\renewcommand*{\thepage}{A\arabic{page}}

