# Improving the Optimization in Model Predictive Controllers: Scheduling Large
Groups of Electric Vehicles

Leoni Winschermann, Marco E.T. Gerards, Johann Hurink
Department of Electrical Engineering, Mathematics and Computer Science
University of Twente
Enschede, The Netherlands
<EMAIL_ADDRESS>
###### Abstract

In parking lots with large groups of electric vehicles (EVs), charging has to
happen in a coordinated manner, among others, due to the high load per vehicle
and the limited capacity of the electricity grid. To achieve such
coordination, model predictive control can be applied, thereby repeatedly
solving an optimization problem. Due to its repetitive nature and its
dependency on the time granularity, optimization has to be (computationally)
efficient.

The work presented here focuses on that optimization subroutine, its
computational efficiency and how to speed up the optimization for large groups
of EVs. In particular, we adapt Focs, an algorithm that can solve the
underlying optimization problem, to better suit the repetitive set-up of model
predictive control by adding a pre-mature stop feature. Based on real-world
data, we empirically show that the added feature speeds up the median
computation time for 1-minute granularity by up to 44%. Furthermore, since
Focs is an algorithm that uses maximum flow methods as a subroutine, the
impact of choosing various maximum flow methods on the runtime is
investigated. Finally, we compare Focs to a commercially available solver,
concluding that Focs outperforms the state-of-the-art when making a full-day
schedule for large groups of EVs.

###### Index Terms:

algorithm, optimization, computational performance, real data, electric
vehicle

##  I Introduction

With the on-going transition to electric mobility, the number of electric
vehicles (EVs) in the Netherlands is increasing rapidly [1]. Especially near-
work locations have been found to be good locations to charge large groups of
EVs, partially due to the long stays [2] and the potential availability of
solar energy [3]. However, due to limited grid capacity [4], high
synchronicity [5] and possible temporal mismatches with said solar energy
production [6], charging has to occur in a coordinated manner. For successful
coordination, a controller needs information on the charging sessions, for
example the (expected) energy demand and departure time of an EV in the
parking lot. In practice, the EVs currently do not communicate such
information [7]. Therefore, next to the sheer amount of charge required per
vehicle, the biggest challenge in managing the charging processes of the EVs
is the absence of information.

To bridge this information gap, model predictive control (Mpc) can be applied
(e.g., [8]). Based on the information available at the current point in time
and an underlying system model, Mpc generates and solves a suitable offline
optimization problem and an (in time) initial portion of this solution is
realized. This process is repeated frequently over the given time horizon.
Specifically to solve this offline problem, the Flow-Based Offline Charging
Scheduler (Focs) has been developed [9]. Focs is a polynomial time algorithm
which in theory has a complexity of $\mathcal{O}(n^{2}\mu)$, where $n$ is the
number of EVs to be charged and $\mu$ is the complexity of the maximum flow
algorithm embedded in Focs.

Notably, research has shown that for some algorithms theoretical and empirical
runtime do not align. For example, the simplex algorithm, while theoretically
exponential in runtime, performs well in practice. The other way around,
Ellipsoid methods for linear programs are known to be polynomial on paper,
while performing poorly in practice. This motivates us to investigate the
empirical runtime of the novel algorithm Focs based on real world EV data.

The main contributions of the work presented here are:

  * •

Introduction and validation of a pre-mature stop feature that increases
computational efficiency in model predictive controllers using Focs.

  * •

Extensive comparative runtime analysis of Focs, and its validation using real
world EV charging data.

The remainder of the paper is organized as follows. Section II discusses the
mathematical model of the considered EV scheduling problem (Section II-A),
briefly introduces Focs (Section II-B) and proposes a pre-mature stop feature
for Focs (Section II-C) expected to increase its efficiency for Mpc
applications. Then, Section III specifies the experimental setup and gives
some information on the implementation of the used algorithms. Furthermore,
the section describes the real world data set used for the experiments
(Section III-C). Section IV presents the experimental results, followed by a
discussion in Section V and a conclusion in Section VI.

##  II EV scheduling problem and algorithms

In this section, we describe the considered EV scheduling problem in terms of
notation, constraints and considered objective function. Furthermore, we
briefly introduce Focs, an algorithm that can be used to determine an optimal
solution for said problem. For more details on Focs, we refer to [9].

###  II-A Problem definition

We consider jobs $1,...,n$, where each job $j$ corresponds to a pending
charging session (or EV) with associated energy requirement $e_{j}$, arrival
time $a_{j}$, departure time $d_{j}$ and job-specific charging power limit
$p_{j}^{\max}$. We denote the set of all jobs by $[n]$. We consider a problem
instance feasible if

$e_{j}\leq p_{j}^{\max}(d_{j}-a_{j})\hskip 10.0pt\forall j\in[n].$

Moreover, we discretize the time horizon $[\min_{j}a_{j},\max_{j}d_{j}]$ into
atomic intervals by using all arrival and departure times of EVs as
breakpoints. The resulting (ordered) sequence of breakpoints we denote as
$t_{1},\dots,t_{m+1}$, resulting in $m\leq 2n-1$ atomic intervals
$I_{i}=[t_{i},t_{i+1}]$ for $i=1,...,m$.

To relate jobs to atomic intervals, we denote by $J(i)$ the jobs that are
available in interval $I_{i}$, i.e.,

$J(i)=\\{j|(a_{j}\leq t_{i})\land(t_{i+1}\leq d_{j})\\}.$

Similarly, $J^{-1}(j)$ is defined as the set of indices $i$ of the intervals
$I_{i}$ where $j$ is available. Finally, as atomic intervals are in general
not unit-sized, we introduce maximum energy limits
$e_{i,j}^{\max}=p_{j}^{\max}(t_{i+1}-t_{i})$ per job $j$ and interval $I_{i}$.

As decision variables, let $e_{i,j}$ be the energy that EV $j$ charges in
$I_{i}$. Note that preemption is allowed, meaning the charging can be
suspended for a few intervals to be continued at a later interval. In this
work, we do not consider V2G, implying that $e_{i,j}\geq 0$.

Summarizing, the EV model is constrained by

$\displaystyle\sum_{i\in J^{-1}(j)}e_{i,j}$ $\displaystyle\geq e_{j}$
$\displaystyle\forall j\in[n]$ (1a) $\displaystyle e_{i,j}$ $\displaystyle\geq
0$ $\displaystyle\forall j\in[n],i\in J^{-1}(j)\hskip 3.0pt$ (1b)
$\displaystyle e_{i,j}$ $\displaystyle\leq e_{i,j}^{\max}$
$\displaystyle\forall j\in[n],i\in J^{-1}(j).$ (1c)

From a grid perspective, the aggregated power level resulting from an EV
schedule is of interest. For a given schedule, the average aggregated power
level $p_{i}$ in atomic interval $I_{i}$ is given by

$p_{i}=\frac{\sum_{j\in J(i)}e_{i,j}}{t_{i+1}-t_{i}}.$

The two most frequently considered objectives for the overall grid usage are
the $\ell_{\infty}$ and $\ell_{2}$-norms of the aggregated power. Focs was
developed for minimization of objective functions

$F\left(p_{1},\dots,p_{m}\right)$ (2)

that are convex, differentiable and for which increasing any value $p_{i}$
increases the value of the objective function. The complexity of Focs is
independent on the exact form of the objective function. In this work, we
specifically consider the square of the $\ell_{2}$-norm given by

$F\left(p_{1},\dots,p_{m}\right)=\sum_{i=1}^{m}L_{i}\cdot\left(p_{i}^{2}\right)$
(3)

where $L_{i}$ is a normalization term to account for the length of the
interval.

###  II-B The general working of Focs

Focs is a recursive algorithm that finds an optimal solution to the problem
described above. To this end, it makes use of various properties of optimal
solutions to this problem. First, any optimal aggregated power profile that
minimizes an objective function (2) is a step function in the sense that
within each atomic interval, the aggregated power is constant. Based on this,
Focs iteratively identifies a critical interval where the highest speed is
needed, determines the minimum speed for this interval, and fixes the schedule
for this interval, ultimately minimizing the maximum peak. The selected
critical interval and the charge provided in that interval are then removed
from the problem instance.

To determine the critical intervals, Focs computes multiple iterations of
maximum flows. Figure 1 depicts the basic structure of the considered flow
network. Here, at each iteration, Focs adapts the capacities of edges adjacent
to sink node $t$. The iterations until a critical interval is identified, form
a round $r$. The edge capacities
$g_{r,k}:\\{I_{i}|i\in[m]\\}\rightarrow\mathbb{R}_{\geq 0}$ in round $r$ and
iteration $k$ represent a fill-level that is increased over the iterations
until critical interval(s) are found.

Figure 1: Schematic of flow network structure of EV charging schedule.

###  II-C Adapting Focs for model predictive control

Focs was developed to determine the optimal solution of the offline EV
scheduling problem described in Section II-A. To this end, it solves the
instance over the entire planning horizon. However, for Mpc, often only the
next (few) time intervals and their control actions are of interest. Based on
this, we propose an adaptation to Focs. For convenience, we assume in the
following that only the first time interval is of interest. The proposed
concept extends naturally to include the first few intervals.

$j$ | $a_{j}$ | $d_{j}$ | $e_{j}$ | $p_{j}^{\max}$
---|---|---|---|---
1 | 0 | 2 | 4 | 2
2 | 1 | 3 | 2 | 1

Figure 2: Intermediate states of Focs with pre-mature stop for an example
instance, tracked over rounds and iterations.

Focs recursively identifies intervals that require the highest aggregated
power in an optimal schedule (so-called critical intervals), determines the
partial schedule corresponding to that interval, and repeats this with the
remaining problem instance. In particular, each of those recursion steps fixes
a part of the optimal schedule. Therefore, if we are only interested in the
schedule of the first interval, we can pre-maturely terminate Focs once the
schedule for that interval was found, without losing any information on the
first interval.

We illustrate this concept in Figure 2 for an instance with two EVs and three
intervals. The annotations at the edges of the flow network are of the form
$x/y$, where $x$ is the flow going through the edge, and $y$ is the edge
capacity. Moreover, the green shaded blocks are the identified critical
intervals. Those are fixed, and in the last step combined into the algorithm
output. After finding the schedule for the first interval (i.e., when we find
the green block for the first interval in the fifth step in Figure 2), we
terminate the algorithm. Note that even though the schedule for the third
interval has not yet been found and thereby the returned schedule (last power
profile in Figure 2) is not feasible considering the entire problem instance,
there exists an optimal schedule that complements the found partial solution.
In particular, the aggregated power in the first interval is optimal. This
follows by optimality of Focs, and the following lemmas paraphrased from [9].

###### Lemma 1

(Lemma 2 in [9]) If a critical interval $I_{i}$ and its aggregated power are
found, and if there is at least one interval $I_{i^{\prime}}$ left that at
that point has not yet been scheduled, then no scheduled charging work can
feasibly be moved from $I_{i}$ to $I_{i^{\prime}}$ for any such
$I_{i^{\prime}}$.

###### Lemma 2

(Lemma 3 in [9]) The aggregated power of the subsequently identified (sets of)
critical intervals by Focs is strictly decreasing. In particular, if Focs
fixes the schedule for $I_{i}$ before that of $I_{i^{\prime}}$, then

$p_{i}>p_{i^{\prime}}$

in the optimal solution.

For the precise formulation and proof of the lemmas, we refer the reader to
[9].

Notably, for EV parking lots near office locations, the optimal schedule over
a day typically resembles a concave curve. The pre-mature stop feature we
propose here terminates the algorithm once the first interval has been
scheduled, where intervals are scheduled by order of their aggregated power in
the optimal solution. Therefore, the time of day at which we apply the pre-
mature stop feature impacts its effectiveness. In particular for an office
parking lot, we expect the greatest reduction in computation time with pre-
mature stop compared to the unrestricted Focs to be around noon. In the
following, we refer to the implementation of Focs with pre-mature stop as
Focs-pm.

##  III Methods

In this section, we provide some details on the implementation of the three
tested optimization models (Quadratic program in Gurobi (QP), Focs and Focs
with pre-mature stop (Focs-pm)), after which we describe the experimental
setup and real-world data set used.111The source code with all implementations
and the experimental setup can be found under
https://github.com/lwinschermann/FlowbasedOfflineChargingScheduler.

###  III-A QP implementation

Based on the problem constraints (1) and the quadratic objective function (3),
we use the commercial solver Gurobi [10] to solve this problem. Gurobi is used
for comparison with state-of-the art solvers. Note that while Focs leads to
optimal solutions also for more general objective functions, for practical
reasons we compare its runtime to a Gurobi model and for that we restrict us
to a quadratic objective function.

###  III-B Focs implementation

The Focs implementation used for the empirical runtime analysis of Focs has
been developed specifically as a proof of concept for this paper. The code is
written in python, heavily relying on the networkx package [11]. This package
is chosen for its user friendly flow models, and because various maximum flow
algorithms are readily available within the package. An overview of the
considered algorithms is provided in Table I.

TABLE I: Complexity table for four maximum flow algorithms. Here, $N$, $M$ and $U$ are respectively the number of nodes, edges and the maximum capacity in a given flow network. Max flow algorithm | Complexity [12] | Networkx complexity | Complexity on Focs network
---|---|---|---
shortest_augmenting_path() | $\mathcal{O}(NMU)$ | $\mathcal{O}(N^{2}M)$ | $\mathcal{O}(n^{4})$
edmonds_karp() |  $\mathcal{O}(M^{2}N)$ or $\mathcal{O}(m^{2}\log U)$ | $\mathcal{O}(M^{2}N)$ | $\mathcal{O}(n^{5})$
preflow_push() | $\mathcal{O}(MN\log(\frac{N^{2}}{M}))$ | $\mathcal{O}(N^{2}\sqrt{M})$ | $\mathcal{O}(n^{3})$
dinitz() |  $\mathcal{O}(N^{2}M)$ or $\mathcal{O}(MN\log U)$ | $\mathcal{O}(N^{2}M)$ | $\mathcal{O}(n^{4})$

Next to that, we also implemented the adapted Focs with pre-mature stop as
introduced in Section II-C by adding a feature that stops the algorithm once
the relevant parts of the schedule are computed. In the following, we refer to
Focs-pm for results achieved with this pre-mature stop.

###  III-C Dataset

The experiments described in this paper use data collected at a real life EV
charging parking lot in Utrecht, the Netherlands [13]. The data was collected
between September 1st 2022 and August 31st 2023, resulting in a total of 13694
charging sessions. Each session specifies (among others) an arrival time,
departure time and total energy charged during the session. At the parking
lot, most charging stations have two sockets each. If only one socket is
occupied, the maximum charging power is 22 kW. Else, it is 11 kW per socket.
Therefore, for the experiments, we set the maximum charging power per session
to either 11 kW or 22 kW, depending on whether the average observed charging
power is below or above 11 kW. From the 13694 sessions in the data, a maximum
of 113 sessions was recorded in a single day.

###  III-D Numerical experiments

The main focus of this paper is on optimization in Mpc for EV charging parking
lots and in particular its efficiency. To this end, we present numerical
experiments based on real world data, evaluating the empirical runtime of Focs
and Focs-pm compared to the benchmark results by Gurobi.

The main parameters to classify the experiments are the

  * •

instance size $n$, i.e., the number of EV charging sessions;

  * •

time granularity, i.e., either 1 minute, 15 minutes, 30 minutes or 1 hour;

  * •

used maximum flow method in Focs (and Focs-pm).

Given a set of parameters for the experiment, we randomly sample $n$ sessions
from the dataset described in Section III-C and apply all three optimization
models. We record the CPU runtime using the function time.process_time() in-
built in the python package time. In particular, we record the CPU runtime for
building the models (QP, Focs, Focs-pm) and solving them. For a given model
and experiment, the total runtime is the sum of CPU time taken to build and
solve the model.

As pointed out in Section II-C, Focs-pm is most effective when the first
interval is expected to have one of the highest power peaks in the optimal
profile. Therefore, next to solving the entire instance over one day, we also
repeat the experiment assuming optimization starting at noon. In that case,
due to the characteristics of an office parking lot, the first interval is
expected to have quite a high power. To get meaningful results, we repeat this
process 500 times, starting at the instance sampling. The values reported in
this work are the median runtimes over those 500 runs. All experiments are run
on an Intel Xeon E5-2630 v3 processor.

##  IV Results

In this section, we present and analyze the results of the runtime
experiments. In particular, we are interested in the efficiency and therefore
usability of the various algorithms for Mpc in (large) EV parking lots. To
this end, we focus on the dependency of runtime on the instance size $n$. In
Section IV-A we first discuss the pre-mature stop feature for Focs introduced
in Section II-C, and how it reflects in the empirical results. Then, in
Section IV-B, we focus on the impact that the different maximum flow methods
have on the performance of Focs. Lastly, we compare the impact of various time
granularities under invariant maximum flow methods in Section IV-C.

###  IV-A Focs-pm

As discussed in Section II-C, an Mpc approach is mostly interested in the next
(few) control action(s). Therefore, we propose Focs-pm as an adapted Focs that
speeds up computation by stopping the optimization pre-maturely when the
control action for the first upcoming intervals is known. While by design we
expect the feature to speed up computation, this section investigates the
extend of the effect, and validates it empirically.

Figure 3 shows the total runtimes of Focs and Focs-pm relative to instance
size $n$ for the four time granularities, and maximum flow method
shortest_augmenting_path(). All instances are solved for the full day. First,
we note that the runtimes of Focs and Focs-pm are very similar for full day
instances. Comparing the values over all instance sizes, the improvement of
Focs-pm relative to Focs is at most 8%, 11%, 12% and 15% for time
granularities of 1 minute, 15 minutes, 30 minutes and 1 hour respectively. The
average improvements are 0%, 2%, 6% and 9% respectively. We see that
especially for finer time granularities, the computational gain of the feature
is rather small for this use case.

(a) 60s time granularity

(b) 15m time granularity

(c) 30m time granularity

(d) 1h time granularity

Figure 3: Runtimes relative to instance size, for various time granularities,
solved for the whole day and using shortest_augmenting_path().

In the analysis above, we considered instances for a full day. In a next step,
we present results of experiments that consider instances starting at noon.

(a) 60s time granularity

(b) 15m time granularity

(c) 30m time granularity (d) 1h time granularity

Figure 4: Runtimes for instances solved at noon using
shortest_augmenting_path().

Figure 4 shows the total runtimes of Focs and Focs-pm relative to instance
size $n$ for the four time granularities, and maximum flow method
shortest_augmenting_path(). However, as opposed to Figure 3, instances are
solved starting from noon. As expected, for situations with the same instance
size and time granularity, solving the (smaller) problem at noon is faster
than solving the entire day. Furthermore, the added value of the pre-mature
stop feature is even more pronounced in the noon experiments. For time
granularities of 1 minute, 15 minutes, 30 minutes and 1 hour, the experiments
show runtime improvements of at most 44%, 39%, 37% and 38% respectively for up
to 400 EVs. That maximum improvement of the median runtime corresponds to
instance sizes 400, 390, 290 and 340, all rather large instances. The average
improvements amount to 33%, 28%, 26% and 23% respectively.

###  IV-B Maximum flow method comparison

The efficiency of an algorithm in the context of an Mpc among others depends
on the implementation. As mentioned in Section II-B, Focs relies on solving
multiple maximum flow problems. Hereby, in theory it does not matter which
concrete algorithm is used. However, since according to Table I the available
maximum flow algorithms are non-linear, this section investigates their impact
on the scalability of Focs for large groups of EVs.

As mentioned in Section III-B, networkx provides a number of functions that
compute maximum flows. In particular, our experiments use the
shortest_augmenting_path(), edmonds_karp(), preflow_push() and dinitz()
functions from the networkx package. Table I provides an overview of those
functions, and the following three properties. First, per function, we cite
the theoretical complexity of its namesake based on [12]. Note that the
implementations of shortest_augmenting_path(), edmonds_karp(), preflow_push()
and dinitz() likely do not strictly follow the original algorithms ([14],
[15], [16]/[17] and [18] respectively). Since their publication, various
improvements and variants have been introduced that are usually referred to by
the same name. Secondly, Table I presents the theoretical complexity as
reported by networkx222last accessed 19 December 2023. Finally, the last
column takes the complexity reported by networkx and bounds the theoretical
complexity specifically for flow networks of Focs instances. We note that the
number of nodes (compare flow network in Figure 1) is bound by 2 plus the
number of jobs and intervals, being $2+n+(2n-1)=3n+1$. The number of edges is
roughly bound by $n+n(2n-1)+(2n-1)=2n^{2}+2n-1$. We derive the theoretical
complexity in terms of the input size of Focs instances by substituting this
in the bound reported by networkx. Since these maximum flow solvers are only
used in the solving of Focs (and therefore do not influence the model
building), we restrict the discussion of the results to the solving part of
the runtime.

Figure 5 presents the solving time of Focs for a full day and quarterly time
granularity. Among the methods using augmenting paths, we see a clear
difference in performance between methods. While the time taken to solve Focs
using shortest_augmenting_path() seems to grow almost linearly for instance
sizes up to 400, solving times with edmonds_karp() and dinitz() clearly
increase non-linearly. Finally, the preflow-push method behaves similarly to
shortest_augmenting_path().

If we compare this to the theoretical complexity of the various maximum flow
methods embedded in networkx (see Table I), the most striking observation is
that while preflow_push() has the smallest theoretical runtime on Focs
networks ($\mathcal{O}(n^{3})$), the theoretically quartic runtime of
shortest_augmenting_path() ($\mathcal{O}(n^{4})$) outperforms the other tested
methods. Edmonds_karp() shows the highest theoretical complexity, and
empirically is the second slowest method, overtaken only by dinitz().

Overall, shortest_augmenting_path() appears to be the dominating method in our
experiments with respect to runtime. The experiments further indicate that the
proper choice of maximum flow method greatly increases the usability of Focs
for EV scheduling.

Figure 5: Runtimes for solving only for quarterly granularity, comparing
maximum flow algorithms applied in Focs.

###  IV-C Runtime comparison

In the following, we compare Focs to the QP implementation for comparison with
(commercial) optimization tools already available for Mpcs. The total runtime
(i.e., the sum of the model building and solving), relative to instance size
$n$ for the four time granularities is depicted in Figure 6. Note that for
Focs, we assume that the flow method shortest_augmenting_path() is applied.
All instances are solved for the full day.

(a) 60s time granularity

(b) 15m time granularity

(c) 30m time granularity

(d) 1h time granularity

Figure 6: Runtimes relative to instance size, for various time granularities,
solved for the whole day and using shortest_augmenting_path().

First, we note that across subfigures, the quadratic program solved in Gurobi
initially either outperforms Focs or runs at similar speed. However, we
observe a gradually steeper increase in runtime relative to instance size for
QP compared to the Focs model. As a result, the median runtime of Focs
outperforms QP starting from instance sizes of approximately 90, 120, 90 and
80 for time granularities of 1 minute, 15 minutes, 30 minutes and 1 hour
respectively. Furthermore, we observe a sudden steep increase in the runtime
of the QP implementation, at instance sizes from around 90, 130, 170 and 230
in Figures 6a, 6b, 6c and 6d respectively. This behaviour possibly reflects
that a certain memory threshold is reached after which solving becomes more
costly for Gurobi due to swapping memory.

Overall, as expected, a finer time granularity for a given instance size
results in a larger runtime for both the QP and Focs. Notably, the slope of
the runtime for Focs seems to be almost linear for the instance sizes
considered, except in Figure 6a. There, a non-linear increase can be observed.
To give an indication of this growth, the median runtimes for Focs and QP with
instance size 400 under time granularity 1 minute are respectively 36 and 62
seconds.

In a next step, we disaggregate the runtimes into model building and solving
components (see Figure 7).

(a) Focs

(b) QP

Figure 7: Runtimes disaggregated into model building and solving for 15 m
granularity.

Figure 7a presents the results for Focs, whereas Figure 7b presents the
results for QP. Both figures are stackplots, depicting the total runtime as
sum of the model building (bottom part, orange shading) and solving (top part,
blue shading). The results correspond to experiments with time granularity of
15 minutes and shortest_augmenting_path() applied within Focs.

Notably, Focs spends most of its computation time on solving. The almost
negligible contribution of model building can be attributed to the mostly
dictionary-based implementation of flows in the networkx package. For QP, on
the other hand, model building makes up a considerable portion of the runtime.
Moreover, Figure 7b clearly attributes the earlier observed sudden steep
increase in the runtime of the QP (discussed above) to its solving. This can
be seen in the continued smooth increase in the bottom (orange) area, whereas
approximately at instance size 130 the contribution of the solving time to the
total runtime increases drastically.

##  V Discussion

In the previous section, we have presented the empirical results of our
experiments. Next, we want to shortly discuss some of the limitations of this
research, and put the results in perspective to physical parking lots for EV
charging with Mpc.

First, we consider the implementation of Focs as used for the experiments.
While the overall Focs method builds on augmentation of flows, within the used
implementation the number of considered intervals is reduced between
consecutive iterations and a new maximum flow with the empty flow as
initialization is executed. A more efficient implementation would use the flow
from the previous iteration, and augment on top of that. Especially if Focs
uses a maximum flow method that relies on augmenting paths, this would speed
up the optimization.

Moreover, in Mpc settings, both efficiency and quality of the optimization are
crucial. Therefore, it is important to note that for the conclusions of
Section IV-B we have not considered the exact form of the found schedule but
focused solely on the runtime of the various flow methods. However, in the
setting of EV scheduling, the fairness, robustness and user acceptability of
the found schedule itself are of great importance. To illustrate, consider the
example in Figure 8.

$j$ | $a_{j}$ | $d_{j}$ | $e_{j}$ | $p_{j}^{\max}$
---|---|---|---|---
1 | 0 | 2 | 2 | 2
2 | 0 | 1 | $\frac{1}{2}$ | 2
3 | 1 | 2 | $\frac{1}{2}$ | 2
4 | 0 | 2 | 2 | 2

Figure 8: Example of an instance where shortest_augmenting_path() results in a
schedule with bias for one job over the other.

Here, jobs 1 and 4 share the same properties. The example provides two
schedules. In the schedule on the left those two jobs are scheduled
consecutively in disjoint intervals, whereas on the right they are scheduled
in parallel. Both schedules are optimal solutions to the problem with
constraints (1) minimizing objective function (3). One may argue, however,
that the schedule on the right is fairer to the users and more robust to early
departure of either job. Such considerations were left out of scope of this
work, but should be considered in future work considering optimization in Mpc
settings.

Finally, in the comparison of Focs to QP, we do not take the adaptivity of the
models in an Mpc setting into account. While the model building for QP takes
considerable time, Gurobi is known for efficiently altering and re-solving its
models, using previous results in its initialization. This was left out of
scope for the comparison. While we do recommend future work in that direction,
we first suggest to conduct theoretical research on the efficient
initialization of Focs and to implement the subsequent results to allow for a
reasonable comparison.

##  VI Conclusion

This work investigated potential algorithms for an optimization subroutine in
model predictive controllers for EV scheduling. To this end, we introduced a
pre-mature stop feature for Focs, an algorithm suitable for this subroutine.
The proposed feature can decrease runtime in Mpc settings. It is particularly
effective close to peak hours, for example around noon at an office parking
lot when parking lot occupancy is expected to be at its peak. Using a 15
minute granularity, the proposed variant of Focs shows a median runtime
improvement of on average 28% compared to the original unadapted version to
calculate a charging schedule at noon. Moreover, this work has illustrated the
considerable impact of the chosen maximum flow method for Focs. Finally, for
large EV parking lots, the presented runtimes show the competitiveness of Focs
with commercial state-of-the-art optimization software.

Future work may dive deeper into the integration of Focs into an Mpc, and into
the properties of the found schedules depending on the used maximum flow
method.

## Acknowledgements

This research is conducted within the SmoothEMS met GridShield project
subsidized by the Dutch ministries of EZK and BZK (MOOI32005).

## References

  * [1] R. Wolbertus and R. van den Hoed, _Plug-in (Hybrid) Electric Vehicle Adoption in the Netherlands: Lessons Learned_. Cham: Springer International Publishing, 2020, pp. 121–143.
  * [2] N. Sadeghianpourhamami, N. Refa, M. Strobbe, and C. Develder, “Quantitive analysis of electric vehicle flexibility: A data-driven approach,” _International Journal of Electrical Power & Energy Systems_, vol. 95, pp. 451–462, 2018.
  * [3] G. J. Osório, M. Gough, M. Lotfi, S. F. Santos, H. M. Espassandim, M. Shafie-khah, and J. P. Catalão, “Rooftop photovoltaic parking lots to support electric vehicles charging: A comprehensive survey,” _International Journal of Electrical Power & Energy Systems_, vol. 133, p. 107274, 2021.
  * [4] J. W. Eising, T. van Onna, and F. Alkemade, “Towards smart grids: Identifying the risks that arise from the integration of energy and transport supply chains,” _Applied Energy_ , vol. 123, pp. 448–455, 2014.
  * [5] K. Turitsyn, N. Sinitsyn, S. Backhaus, and M. Chertkov, “Robust broadcast-communication control of electric vehicle charging,” in _2010 First IEEE International Conference on Smart Grid Communications_ , 2010, pp. 203–207.
  * [6] M. Cañigueral and J. Meléndez, “Flexibility management of electric vehicles based on user profiles: The arnhem case study,” _International Journal of Electrical Power & Energy Systems_, vol. 133, p. 107195, 2021.
  * [7] G. Vecchio and L. Tricarico, ““may the force move you”: Roles and actors of information sharing devices in urban mobility,” _Cities_ , vol. 88, pp. 261–268, 2019.
  * [8] Q. Zhou and C. Du, “A quantitative analysis of model predictive control as energy management strategy for hybrid electric vehicles: A review,” _Energy Reports_ , vol. 7, pp. 6733–6755, 2021.
  * [9] L. Winschermann, M. E. T. Gerards, A. Antoniadis, G. Hoogsteen, and J. Hurink, “Relating electric vehicle charging to speed scaling with job-specific speed limits,” 2023.
  * [10] Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,” 2023. [Online]. Available: https://www.gurobi.com
  * [11] A. A. Hagberg, D. A. Schult, and P. J. Swart, “Exploring network structure, dynamics, and function using networkx,” in _Proceedings of the 7th Python in Science Conference_ , G. Varoquaux, T. Vaught, and J. Millman, Eds., Pasadena, CA USA, 2008, pp. 11 – 15.
  * [12] O. Cruz-Mejía and A. N. Letchford, “A survey on exact algorithms for the maximum flow and minimum-cost flow problems,” _Networks_ , vol. 82, no. 2, pp. 167–176, 2023.
  * [13] B. Nijenhuis, L. Winschermann, N. Bañol Arias, G. Hoogsteen, and J. L. Hurink, “Protecting the distribution grid while maximizing ev energy flexibility with transparency and active user engagement,” in _CIRED Porto Workshop 2022: E-mobility and power distribution systems_ , vol. 2022, 2022, pp. 209–213.
  * [14] L. R. Ford and D. R. Fulkerson, “Maximal flow through a network,” _Canadian Journal of Mathematics_ , vol. 8, p. 399–404, 1956.
  * [15] J. Edmonds and R. M. Karp, “Theoretical improvements in algorithmic efficiency for network flow problems,” _J. ACM_ , vol. 19, no. 2, p. 248–264, apr 1972\.
  * [16] A. Karzanov, “Determining the maximal flow in a network by the method of preflows,” _Doklady Mathematics_ , vol. 15, p. 434–437, 02 1974.
  * [17] A. V. Goldberg and R. E. Tarjan, “A new approach to the maximum-flow problem,” _J. ACM_ , vol. 35, no. 4, p. 921–940, oct 1988.
  * [18] Y. Dinitz, “Algorithm for solution of a problem of maximum flow in networks with power estimation,” _Soviet Math. Dokl._ , vol. 11, pp. 1277–1280, 01 1970.


# radioactivedecay: A Python package for radioactive decay calculations

Alex Malins  Thom Lemoine

††margin:  DOI:
https://doi.org/10.21105/joss.03318https://doi.org10.21105/joss.03318 Software
• https://github.com/openjournals/joss-
reviews/issues/3318https://doi.orgReview •
https://github.com/radioactivedecay/radioactivedecayhttps://doi.orgRepository
• https://doi.org/10.5281/zenodo.6334651https://doi.orgArchive Editor:
https://github.com/kellyrowlandhttps://doi.orgKelly Rowland
Reviewers: • https://github.com/munkmhttps://doi.org@munkm •
https://github.com/shyamdhttps://doi.org@shyamd •
https://github.com/kellyrowlandhttps://doi.org@kellyrowland Submitted: 30
April 2021
Published: 16 March 2022 License
Authors of papers retain copyright and release the work under a Creative
Commons Attribution 4.0 International License
(http://creativecommons.org/licenses/by/4.0/https://doi.orgCC BY 4.0).

## Summary

radioactivedecay is a Python package for radioactive decay modelling. It
contains functions to fetch decay data, define inventories of nuclides and
perform decay calculations. The default nuclear decay dataset supplied with
radioactivedecay is based on ICRP Publication 107, which covers 1252
radioisotopes of 97 elements. The code calculates an analytical solution to a
matrix form of the decay chain differential equations using double or higher
precision numerical operations. There are visualization functions for drawing
decay chain diagrams and plotting activity decay curves.

## Statement of Need

Calculations for the decay of radioactivity and the ingrowth of progeny
underpin the use of radioisotopes in a wide range of research and industrial
fields, spanning from nuclear engineering, medical physics, radiation
protection, environmental science and archaeology to non-destructive testing,
mineral prospecting, food preservation, homeland security and defence.
radioactivedecay is an open source, cross-platform package for decay
calculations and visualization. It supports decay chains with branching decays
and metastable nuclear isomers. It includes a high numerical precision decay
calculation mode, which resolves numerical problems with using double-
precision floating-point numbers to calculate decay chains involving
radionuclides with disparate half-lives (Bakin et al., 2018).

This set of features distinguishes radioactivedecay from other commonly-used
decay packages, such as Radiological Toolbox (Hertel et al., 2015) and PyNE
(Scopatz et al., 2012). Radiological Toolbox is a closed-source Windows
application, so it is not easily scriptable and its use of double-precision
arithmetic makes it susceptible to numerical round-off errors. PyNE uses
approximations to help mitigate numerical issues, however these may
potentially affect accuracy. Moreover as of v0.7.5, PyNE does not correctly
model metastable nuclear isomers within decay chains, which means, for
example, it cannot simulate the production of ${}^{99m}\textrm{Tc}$ from
${}^{99}\textrm{Mo}$ for medical imaging applications.

## Theory and Implementation

radioactivedecay implements the solution to the decay differential equations
outlined by Amaku et al. (2010). If vector $\mathbf{N}$ contains the number of
atoms of each radionuclide in a system, its elements $N_{i}$ can be ordered
such that no progeny (either first or subsequent generation) of radionuclide
$i$ has itself an index lower than $i$. Ordering in this manner is possible
because natural radioactive decay processes do not increase the mass number of
the decaying radionuclide, and there are no cyclic decay chains where
radionuclide $i$ can decay to other radionuclides then reform itself (Ladshaw
et al., 2020). Note metastable nuclear isomers have distinct indices from
their ground states in $\mathbf{N}$.

The radioactive decay chain differential equations expressed in matrix form
are:

$\frac{\mathrm{d}\mathbf{N}}{\mathrm{d}t}=\varLambda\mathbf{N}.$ (1)

$\varLambda$ is a lower triangular matrix with elements:

$\varLambda_{ij}=\begin{cases}0&\textrm{for }i<j,\\\ -\lambda_{j}&\textrm{for
}i=j,\\\ b_{ji}\lambda_{j}&\textrm{for }i>j.\end{cases}$ (2)

$\lambda_{j}$ is the decay constant of radionuclide $j$, and $b_{ji}$ is the
branching fraction from radionuclide $j$ to $i$. $\varLambda$ is
diagonalizable so its eigendecomposition can be used to rewrite Equation 1 as:

$\frac{\mathrm{d}\mathbf{N}}{\mathrm{d}t}=C\varLambda_{d}C^{-1}\mathbf{N}.$
(3)

$\varLambda_{d}$ is a diagonal matrix whose elements are the negative decay
constants, i.e. $\varLambda_{dii}=-\lambda_{i}$. Matrix $C$ and its inverse
$C^{-1}$ are both lower triangular matrices that are calculated as:

$C_{ij}=\begin{cases}0&\text{for }i<j,\\\ 1&\text{for }i=j,\\\
\frac{\sum_{k=j}^{i-1}\varLambda_{ik}C_{kj}}{\varLambda_{jj}-\varLambda_{ii}}&\text{for
}i>j,\end{cases}\quad\text{and}\quad C^{-1}_{ij}=\begin{cases}0&\text{for
}i<j,\\\ 1&\text{for }i=j,\\\ -\sum_{k=j}^{i-1}C_{ik}C^{-1}_{kj}&\text{for
}i>j.\end{cases}$ (4)

The analytical solution to Equation 3 given an initial condition of
$\mathbf{N}(0)$ at $t=0$ is:

$\mathbf{N}(t)=Ce^{\varLambda_{d}t}C^{-1}\mathbf{N}(0).$ (5)

$e^{\varLambda_{d}t}$ is a diagonal matrix with elements
$e^{\varLambda_{d}t}_{ii}=e^{-\lambda_{i}t}$. radioactivedecay evaluates
Equation 5 upon each call for a decay calculation.

Matrices $C$ and $C^{-1}$ are independent of time so they are pre-calculated
and imported from files into radioactivedecay. $C$ and $C^{-1}$ are stored in
sparse matrix data structures to minimize memory use and maximize efficiency
when computing the matrix multiplications in Equation 5. For decay
calculations with double-precision floating-point operations, $C$ and $C^{-1}$
are stored in SciPy (Virtanen et al., 2020) Compressed Sparse Row (CSR) matrix
data structures. Conversely, they are stored in SymPy (Meurer et al., 2017)
SparseMatrix data structures for high numerical precision calculations.

The high numerical precision decay calculation mode resolves numerical issues
arising from using double-precision floating-point numbers for decay
calculations for chains containing nuclides with disparate half-lives. One
example is the decay chain for ${}^{254}\textrm{Es}$, which contains
${}^{238}\textrm{U}$ (4.468 billion year half-life) and ${}^{214}\textrm{Po}$
($t_{1/2}$ is 164.3 $\mu$s half-life). This a 20 orders of magnitude
difference in half-life. Loss of numerical precision inevitably occurs when
evaluating the off-diagonal elements of $C$ and $C^{-1}$ in Equation 4 with
double-precision floating-point numbers (which hold approximately 15 decimal
places of numerical precision). Note loss of precision also occurs in the
converse scenario, i.e. when a decay chain contains radionuclides with similar
half-lives. However this scenario does not occur in the ICRP Publication 107
decay dataset, as the relative difference between half-lives of any two
radionuclides in the same decay chain is always greater than 0.1%.

The default operation of the high precision decay mode is to evaluate Equation
5 using floating-point numbers with 320 significant figures of precision. This
is sufficient precision to ensure accurate results for any physically relevant
decay calculation users may wish to perform. Moreover, computations in the
high precision mode are still fast, taking less than one second on a notebook
equipped with an Intel Core i5-8250U processor.

## Decay & Atomic Mass Datasets

The default dataset supplied with radioactivedecay uses decay data from ICRP
Publication 107 (Eckerman & Endo, 2008) and atomic masses from the Atomic Mass
Data Center (AMDC) (Huang et al., 2021; Kondev et al., 2021; Wang et al.,
2021). Endo et al. (2005) and Endo & Eckerman (2007) describe the development
of the ICRP Publication 107 decay dataset. Raw data from ICRP 107 and AMDC
were converted into dataset files suitable for radioactivedecay in a Jupyter
https://github.com/radioactivedecay/datasetshttps://doi.orgnotebook. Along
with SciPy and SymPy versions of the sparse matrices $C$ and $C^{-1}$, the
dataset files contain radionuclide half-lives, decay constants, progeny,
branching fractions, decay modes and atomic masses. Although there is a
default dataset, radioactivedecay allows the import and use other decay data.

## Main Functionality

Figure 1: Examples of the plotting capabilities of radioactivedecay: (a) Decay
chain diagram for molybdenum-99. (b) Graph showing the decay of 1 kBq of
${}^{99}\textrm{Mo}$ along with the ingrowth of ${}^{99m}\textrm{Tc}$ and a
trace quantity of ${}^{99}\textrm{Tc}$.

The main functionality of radioactivedecay is based around Nuclide, Inventory
and InventoryHP classes. The Nuclide class is used for fetching atomic and
decay data about a single nuclide, such as its atomic mass, half-life, decay
modes, progeny and branching fractions. It creates diagrams of the nuclide’s
decay chain (ex. Figure 1(a)) using the NetworkX library (Hagberg et al.,
2008).

An Inventory can contain multiple nuclides, each with an associated quantity
(the number of atoms of the nuclide). Nuclides can be stable or radioactive.
The decay() method calculates the decay of the radioactive nuclides in an
Inventory, adding any ingrown progeny automatically. The numbers(),
activities(), masses(), and moles() methods output the inventory of nuclides
as different quantities using the atomic data stored in the decay dataset.
Additional activity_fractions(), mass_fractions(), and mole_fractions()
methods provide the relative amounts of each nuclide in the inventory with
respect to different quantities. Plots can be made of the variation of nuclide
activities, masses and moles over time (ex. Figure 1(b)) using Matplotlib
(Hunter, 2007).

The InventoryHP class is the high numerical precision complement of the
Inventory class. It has the same API as the Inventory class, but uses SymPy
high numerical precision routines for all calculations.

## Validation

Decay calculations with radioactivedecay v0.4.2 were cross-checked against
Radiological Toolbox v3.0.0 (Hertel et al., 2015) and PyNE v0.7.5 (Scopatz et
al., 2012) (see Jupyter notebooks in the
https://github.com/radioactivedecay/comparisonshttps://doi.orgcomparisons
repository). Radiological Toolbox employs the ICRP Publication 107 decay data.
Fifty radionuclides were randomly selected and a decay calculation was
performed for 1 Bq of each for a random decay time within a factor of
$10^{-3}$ to $10^{3}$ of the half-life. Differences between decayed activities
reported by each code were within 1% of each other in 64% of cases.
Discrepancies greater than 1% were attributed to rounding differences,
erroneous results from Radiological Toolbox, or numerical issues relating to
decay chains containing radionuclides with disparate half-lives.

A dataset was prepared for radioactivedecay with the same Evaluated Nuclear
Structure Data File (ENSDF, 2019) decay data as used by PyNE v0.7.5. Bugs in
PyNE v0.7.5 cause incorrect decay calculation results for chains containing
metastable nuclear isomers, ${}^{183}\textrm{Pt}$, ${}^{172}\textrm{Ir}$ or
${}^{152}\textrm{Lu}$. Thus the affected chains were not used for the
comparisons. The decay of 1 Bq of every radionuclide was calculated for
multiple decay times varying from $0$ to $10^{6}$ times the radionuclide’s
half-life. The absolute difference between the decayed activities reported by
each code was less than $10^{-13}$ Bq. Relative differences depended on the
magnitude of the activity. Relative errors of greater than 0.1% only occurred
when the calculated activity was less than $2.5\times 10^{-11}$ Bq, i.e. 10
orders of magnitude smaller than the initial activity of the parent
radionuclide. The discrepancies between the two codes were attributed to
methodological differences for computing decay chains with radionuclides with
large disparities between half-lives, and numerical issues arising from
double-precision floating-point operations.

## Limitations

radioactivedecay does not model neutronics, so cannot evaluate radioactivity
produced from activations or induced fission. It does not support external
sources of radioactivity input or removal from an inventory over time. Caution
is required if decaying backwards in time, as this can cause floating-point
overflows when computing the exponential terms in Equation 5.

There are also some limitations associated with the ICRP Publication 107 decay
dataset. It does not contain data for the radioactivity produced from
spontaneous fission decay pathways and the minor decay pathways of some
radionuclides. More details on limitations are available in the
https://radioactivedecay.github.io/overview.html#limitationshttps://doi.orgdocumentation,
Endo et al. (2005), and Endo & Eckerman (2007).

## Acknowledgements

We thank Mitsuhiro Itakura, Kazuyuki Sakuma, colleagues in JAEA’s Center for
Computational Science & Systems, Wolfgang Kerzendorf & Bernardo Gameiro for
their support for this project. We thank Kenny McKee, Daniel Jewell, Ezequiel
Passaro, Hunter Ratliff & Jayson Vavrek for helpful suggestions, and Bjorn
Dahlgren, Anthony Scopatz & Jonathan Morrell for their work on radioactive
decay calculation software. We also thank the editors and reviewers at the
Journal of Open Source software for constructive comments.

## References

Amaku, M., Pascholati, P. R., & Vanin, V. R. (2010). Decay chain differential
equations: Solution through matrix algebra. _Computer Physics Communications_
, _181_(1), 21–23. https://doi.org/10.1088/0952-4746/26/3/N02

Bakin, R. I., Kiselev, A. A., Shvedov, A. M., & Shikin, A. V. (2018).
Computational Errors in the Calculation of Long Radioactive Decay Chains.
_Atomic Energy_ , _123_(6), 406–411. https://doi.org/10.1007/s10512-018-0360-2

Eckerman, K. F., & Endo, A. (2008). ICRP 107: Nuclear Decay Data for
Dosimetric Calculations. _Annals of the ICRP_ , _38_(3), 119.
https://doi.org/10.1016/j.icrp.2008.10.003

Endo, A., & Eckerman, K. F. (2007). _JAEA-Data/Code 2007-021: Nuclear Decay
Data for Dosimetry Calculation - Data for Radionuclides with Half-lives Less
than 10 Minutes_. https://doi.org/10.11484/jaea-data-code-2007-021

Endo, A., Yamaguchi, Y., & Eckerman, K. F. (2005). _JAERI 1347: Nuclear Decay
Data for Dosimetry Calculation; Revised data of ICRP Publication 38_.
https://doi.org/10.11484/jaeri-1347

ENSDF. (2019). _From ENSDF database as of October 4, 2019. Version available
at:_ https://www.nndc.bnl.gov/ensarchivals/

Hagberg, A. A., Schult, D. A., & Swart, P. J. (2008). Exploring network
structure, dynamics, and function using NetworkX. _Proceedings of the 7th
Python in Science Conference (SciPy2008)_ , 11–15.
http://conference.scipy.org/proceedings/SciPy2008/paper_2/

Hertel, N. E., Eckerman, K. F., & Sun, C. (2015). Radiological Tookbox 3.0.0.
_Transactions of the American Nuclear Society_ , _113_ , 977–980.
https://www.ans.org/pubs/transactions/article-38022/

Huang, W. J., Wang, M., Kondev, F. G., Audi, G., & Naimi, S. (2021). The AME
2020 atomic mass evaluation (i). Evaluation of input data, and adjustment
procedures. _Chinese Physics C_ , _45_(3), 030002.
https://doi.org/10.1088/1674-1137/abddb0

Hunter, J. D. (2007). Matplotlib: A 2D graphics environment. _Computing in
Science & Engineering_, _9_(3), 90–95. https://doi.org/10.1109/MCSE.2007.55

Kondev, F. G., Wang, M., Huang, W. J., Naimi, S., & Audi, G. (2021). The
NUBASE2020 evaluation of nuclear physics properties. _Chinese Physics C_ ,
_45_(3), 030001. https://doi.org/10.1088/1674-1137/abddae

Ladshaw, A., Wiechert, A. I., Kim, Y., Tsouris, C., & Yiacoumi, S. (2020).
Algorithms and algebraic solutions of decay chain differential equations for
stable and unstable nuclide fractionation. _Computer Physics Communications_ ,
_246_ , 106907. https://doi.org/10.1016/j.cpc.2019.106907

Meurer, A., Smith, C. P., Paprocki, M., ?ertik, O., Kirpichev, S. B., Rocklin,
M., Kumar, A., Ivanov, S., Moore, J. K., Singh, S., Rathnayake, T., Vig, S.,
Granger, B. E., Muller, R. P., Bonazzi, F., Gupta, H., Vats, S., Johansson,
F., Pedregosa, F., … Scopatz, A. (2017). SymPy: symbolic computing in Python.
_PeerJ Computer Science_ , _3_ , e103. https://doi.org/10.7717/peerj-cs.103

Scopatz, A. M., Romano, P. K., Wilson, P. P. H., & Huff, K. D. (2012). PyNE:
Python for nuclear engineering. _Transactions of the American Nuclear Society_
, _107_ , 985–987. https://www.ans.org/pubs/transactions/article-14978

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T.,
Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., Walt,
S. J. van der, Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A.
R. J., Jones, E., Kern, R., Larson, E., … Vazquez-Baeza, Y. (2020). SciPy 1.0:
fundamental algorithms for scientific computing in Python. _Nature Methods_ ,
_17_(3), 261–272. https://doi.org/10.1038/s41592-019-0686-2

Wang, M., Huang, W. J., Kondev, F. G., Audi, G., & Naimi, S. (2021). The AME
2020 atomic mass evaluation (II). Tables, graphs and references. _Chinese
Physics C_ , _45_(3), 030003. https://doi.org/10.1088/1674-1137/abddaf


Boron nanotube structure explored by evolutionary computations

Tomasz Tarkowski 1 and Nevill Gonzalez Szwacki 2*

1 Institute of Theoretical Physics, Faculty of Physics, University of Warsaw,
Pasteura 5, PL-02093 Warsaw, Poland<EMAIL_ADDRESS>
2 Institute of Experimental Physics, Faculty of Physics, University of Warsaw,
Pasteura 5, PL-02093 Warsaw, Poland<EMAIL_ADDRESS>
* Correspondence<EMAIL_ADDRESS>Tel.: +48-22-5532797

###### Abstract

In this work, we explore the structure of single-wall boron nanotubes with
large diameters (about 21 Å) and a broad range of surface densities of atoms.
The computations are done using an evolutionary approach combined with a
nearest neighbors model Hamiltonian. For the most stable nanotubes, the number
of 5-coordinated boron atoms is about $63\%$ of the total number of atoms
forming the nanotubes, whereas about $11\%$ are boron vacancies. For hole
densities smaller than about 0.22, the boron nanotubes exhibit randomly
distributed hexagonal holes and are more stable than a flat stripe structure
and a quasi-flat B36 cluster. For larger hole densities ($>0.22$) the boron
nanotubes resemble porous tubular structures with hole sizes that depend on
the surface densities of boron atoms.

Keywords: genetic algorithms; nanotubes; nanowires; model Hamiltonian

Copyright (c) 2022 Tomasz Tarkowski and Nevill Gonzalez Szwacki. License: CC
BY 4.0 (https://creativecommons.org/licenses/by/4.0/).

##  1 Introduction

Boron nanotubes (BNTs) have been synthesized for the first time in 2004 [1]. A
magnesium-substituted mesoporous silica template (Mg-MCM-41) was applied by
Ciuparu et al. [1] to prepare pure boron single-wall nanotubes at
$870{}^{\circ}\mathrm{C}$ with uniform diameters ($36\pm 1$ Å) using the
mixture of $\mathrm{BCl}_{3}$ and $\mathrm{H}_{2}$ as gas sources. The authors
attributed the Raman peaks at $210\mathrm{~{}cm}^{-1}$ and between $300$ and
$500\mathrm{~{}cm}^{-1}$ to typical tubular structures, where the first
($210\mathrm{~{}cm}^{-1}$) corresponds to the characteristic radial breathing
mode. In 2010 Liu et al. [2] reported the first large-scale fabrication of
single crystalline multilayered BNTs using boron $(99.99\%)$ and boron oxide
powders $(99.99\%)$ as source materials. The as-synthesized BNTs had lengths
of several micrometers and diameters in a range from 10 to 40
$\mathrm{~{}nm}$. The nanotubes were cataloged as multilayered single
crystalline BNTs with an interlayer spacing of about $3.2$ Å. Moreover, these
nanotubes were experimentally proven to have metallic properties regardless of
their chirality. The metallic behavior of BNTs makes them attractive in the
design of novel electronic nanodevices, such as field-effect transistors,
light-emitting diodes, and field-emission displays or for photosensitive
device applications [3, 4].

The experimental studies were preceded by several theoretical investigations,
most often based on density functional theory (DFT). Quasiplanar [5, 6],
tubular [7, 8], convex, and spherical [9] boron clusters have been
computationally explored. Moreover, the existence of quasiplanar boron
clusters [10, 11] implies that boron fullerenes [12, 13] and caped BNTs [14,
15] may exist because larger in size quasiplanar or planar clusters will tend
to remove dangling edge bonds by forming closed tubular or polyhedral
structures.

Although some reports affirm that BNTs with diameters smaller than 1.7 nm [16,
17] or 2 nm [18] show semiconducting behavior due to the band opening at the
Fermi level through curvature-induced out-of-plane buckling of certain atoms,
later calculations based on second-order Møller–Plesset perturbation theory
[19] and dispersion-corrected DFT calculations [20] showed that the surface
buckling is more likely an artifact of standard DFT approaches.

Thus, BNTs are found to be metallic and independent of diameter and chirality
(armchair or zigzag) [21]. Moreover, there have been a lot of morphologies of
BNTs in contrast to only one morphology of carbon nanotubes (CNTs) [22]. A
more detailed description of the structure and properties of 1D-boron
structures can be found in Refs. [23, 4].

Crystal structure prediction (CSP) is a long-standing challenge in physical
and materials sciences. Several methods have been applied for materials design
[24, 25]. Evolutionary algorithms, such as genetic algorithms (GAs) [25] that
use human evolution mechanisms (such as crossover and mutation), have been
used, for instance, to obtain the most stable structures of prototype
nanotubes composed of particles interacting through a Lennard-Jones potential
[26].

Boron two-dimensional structures [27] conceptually form foundations for BNTs
by “cutting” adequate stripes and “gluing” them along the edges. Our
previously developed methodology based on the application of the floating-
point representation is not the only possible strategy for the CSP of boron
nanostructures [28]. Alternatively, one can use binary (for monoatomic
materials) or integer (for alloys) representations with a fixed crystal
lattice. To obtain the total energy of the system, several approaches can be
used. In our previous work, we used a genetic algorithm [29] combined with DFT
to predict the structure of boron nanowires [28]. In this work, we will use a
GA combined with a DFT-based nearest neighbors model Hamiltonian to predict
the structure of BNTs. We assume that the nanotubes have a perfect cylindrical
shape and their structure is closely related to one-atom-thick sheets of boron
atoms arranged on a hexagonal lattice [27].

##  2 Computational approach and results of simulations

###  2.1 Unit cell definition

The information about the position of atoms and vacancies (“holes”) on a
hexagonal lattice wrapped cylindrically around a nanotube axis is encoded
using a binary representation. The type of considered nanotubes is limited to
an armchair $(n,n)$ chirality as described in Ref. [8]. A cylindrical
coordinate system can be conveniently used for the nanotube description. It is
assumed that the nanotube has periodic boundary conditions (PBCs) in the $z$
direction. The unit cell size is described with two numbers
$n_{\phi},n_{z}\in\mathbb{N}_{+}$, where $n_{\phi}>1$. The lattice constant
(interatomic distance) of the unwrapped hexagonal layer is equal to $a$ (to
plot the BNTs, we have assumed $a=1.675$ Å). The positions of atoms and
vacancies are not encoded in the genotype but are specified by the following
position vectors:

$\vec{r}_{i}=\left(\frac{a\sqrt{6}}{4\sqrt{1-\cos(\pi/n_{\phi})}},\frac{\pi}{n_{\phi}}\cdot\left\lfloor\frac{i}{n_{z}}\right\rfloor,(i\bmod
n_{z})\cdot a+\left\\{\begin{array}[]{ll}0&{\rm if}\ \lfloor
i/n_{z}\rfloor\,{\rm is\ even}\\\ a/2&{\rm
otherwise}\end{array}\right.\right)_{(\rho,\phi,z)}.$ (1)

The use of binary representation, $g=(x_{i})_{i=0}^{c-1}\in\mathbb{B}^{c}$,
where $c=2n_{\phi}n_{z}$, justifies the way the vectors $\vec{r}_{i}$ are
indexed (more details about indices are given in Appendix A). It also implies
the number of atoms in the unit cell equal to $N_{\rm
a}=\\#\\{i\in\iota_{c}\mid x_{i}\\}$ and the number of vacancies equal to
$N_{\rm v}=c-N_{\rm a}$.

The unit cell of the BNT is defined in terms of predicate theory. Predicate
$Q\equiv Q_{0}\wedge Q_{1}\wedge Q_{2}$, where predicates $Q_{0}$, $Q_{1}$,
and $Q_{2}$ state that atoms within the unit cell are connected, that there is
an atomic neighborhood at the boundary of two unit cells along the nanotube
axis, and that there is an atomic neighborhood at the unit cell boundary
perpendicular to the circumference of the nanotube, respectively. The
predicate Q is sufficient to ensure the nanotube connectivity and its closure
at the circumference (see also Appendix A for a more precise formulation of
the predicate Q). It does not restrict any possible valid result of structures
that are equivalent by translational or rotational symmetry.

###  2.2 Fitness function

The fitness function values are calculated using a nearest neighbors model
Hamiltonian. The optimization task consists in finding such configuration of
atoms and vacancies which maximizes the fitness function—the binding energy
per atom, $E_{\rm b}$, calculated within the nearest neighbors model
Hamiltonian. In this model, $E_{\rm b}$ can be defined as:

$E_{\rm
b}(n_{1},n_{2},n_{3},n_{4},n_{5},n_{6})=\frac{1}{N_{\textnormal{a}}}\sum_{i=1}^{6}n_{i}e_{i},$
(2)

where $n_{i}$ is the number of atoms in the unit cell with $i$ nearest
neighbors, while $e_{i}$ is the $E_{\rm b}$ for a structure consisting of
atoms having only $i$ nearest neighbors (e.g., two-atom molecule, single
atomic chain, or a honeycomb structure, _hc_ , for $i$ equal to $1$, $2$ or
$3$, respectively). Obviously, the total number of atoms in the unit cell
$N_{\rm a}=\sum_{i=1}^{6}n_{i}$. The values of $e_{i}$ are taken from our
previous work [30] and they are equal to $e_{1}=1{.}7803$, $e_{2}=5{.}1787$,
$e_{3}=5{.}6504$, $e_{4}=6{.}2522$, $e_{5}=6{.}5718$, and $e_{6}=6{.}5116$ in
units of ${\rm eV}/{\rm atom}$. Since our model Hamiltonian includes only
nearest neighbors interactions and was developed for the case of a one-atom-
thick layer of boron atoms [31, 30], it does not take into account the
curvature of the BNT. Therefore, the $e_{i}$ parameters can be applied to
nanotubes with large radii, where the effect of curvature on the properties of
the nanotube is expected to be negligible.

In the case of carbon, the elastic properties of the nanotubes do not depend
on the diameter already at about 18 Å [32]. The same is true for electronic
properties such as the energy gap [33]. Making the assumption that the pure
boron case situation is the same, $n_{\phi}=23$ is taken which corresponds to
$\rho\approx 10{.}628\,\text{\AA}$. This gives the size of the potential
solution space (described by a general formula of $2^{2n_{\phi}n_{z}}$) equal
to $2^{46n_{z}}$.

The Bernoulli ${\rm B}(1,0{.}5)$ probability distribution [$B(n,p)$ is a
discrete probability distribution, where the Bernoulli random variable $n$ can
have only 0 or 1 as the outcome values, $p$ is the probability of success
($n=1$), and $1-p$ is the probability of failure ($n=0$)] was chosen for the
creation of the first generation of nanotubes, rejecting genotypes not
satisfying the predicate $Q$. A one-point recombination with bit-flipping
mutation ($p=1/c$) applied stochastically with probabilities equal to $p_{\rm
r}=p_{\rm m}=0{.}5$ is employed during the evolution.

The generation size $\mu$ is set to $4n_{\phi}n_{z}$ while parent multiset
size $2k=2n_{\phi}n_{z}$. Stochastic universal sampling (SUS) with linear
ranking selection ($s=2$) is used for the parent selection and selection to
the next generation mechanisms. The genetic search is terminated after
reaching a _plateau_ for the fitness function or, more precisely, when after
$100$ generations the fitness function maximum has not improved by $\Delta
E=1\,{\rm meV}$.

Some of the obtained nanotubes can be equivalent through rotations around the
$z$-axis and translations along this axis. There is a relationship between
natural numbers and the above-defined nanotubes. Encoding
$g\in\mathbb{B}^{2n_{\phi}n_{z}}$ corresponds to certain number
$n=\sum_{i=0}^{2n_{\phi}n_{z}-1}2^{i}x_{i}$, where $\mathbb{B}\ni{\rm
true}\equiv 1\in\mathbb{N}$ and $\mathbb{B}\ni{\rm false}\equiv
0\in\mathbb{N}$. There are classes of abstractions of natural numbers
regarding rotations around the $z$-axis (nanotube axis) and translations along
this axis, which form equivariant nanotubes. Rotations and translations in the
set $\iota_{2n_{\phi}n_{z}}$ can be defined as follows:

$\displaystyle R_{n_{\phi},n_{z}}^{\Delta n_{\phi}}(n)$
$\displaystyle=n^{\prime},\,x_{(i+2\Delta n_{\phi}n_{z})\bmod
2n_{\phi}n_{z}}^{\prime}=x_{i},$ (3) $\displaystyle T_{n_{\phi},n_{z}}^{\Delta
n_{z}}(n)$ $\displaystyle=n^{\prime},\,x_{\left((i+\Delta n_{z})\bmod
n_{z}\right)+\lfloor i/n_{z}\rfloor n_{z}}^{\prime}=x_{i},$ (4)

where $\Delta n_{\phi},\Delta n_{z}\in\mathbb{Z}$. The following identities
hold:

$\displaystyle R_{n_{\phi},n_{z}}^{n_{\phi}}$ $\displaystyle={\rm id},$ (5)
$\displaystyle T_{n_{\phi},n_{z}}^{n_{z}}$ $\displaystyle={\rm id},$ (6)
$\displaystyle R_{n_{\phi},n_{z}}^{\Delta n_{\phi}}$
$\displaystyle=\left(R_{n_{\phi},n_{z}}^{1}\right)^{\Delta n_{\phi}},$ (7)
$\displaystyle T_{n_{\phi},n_{z}}^{\Delta n_{z}}$
$\displaystyle=\left(T_{n_{\phi},n_{z}}^{1}\right)^{\Delta n_{z}},$ (8)
$\displaystyle R_{n_{\phi},n_{z}}^{-1}$
$\displaystyle=R_{n_{\phi},n_{z}}^{n_{\phi}-1},$ (9) $\displaystyle
T_{n_{\phi},n_{z}}^{-1}$ $\displaystyle=T_{n_{\phi},n_{z}}^{n_{z}-1}.$ (10)

It may also turn out that the nanotube obtained as a result of evolution has a
shorter elementary unit cell along $\hat{e}_{z}$ than it would result from the
$n_{z}$ value. This occurs if, and only if, the maximum divisor $d_{n_{z}}$ of
the number $n_{z}$, for which the condition

$T_{n_{\phi},n_{z}}^{n_{z}/d_{n_{z}}}(n)=n,$ (11)

is satisfied for $n$ representing a given nanotube, is greater than $1$. If
maximum divisor $d_{n_{\phi}}$ of $n_{\phi}$, for which analogous condition

$R_{n_{\phi},n_{z}}^{n_{\phi}/d_{n_{\phi}}}(n)=n$ (12)

is met, is greater than $1$ then the nanotube has a nontrivial rotational
symmetry about the axis of the tube. The values of $d_{n_{z}}$ and
$d_{n_{\phi}}$ are assumed to denote the maximum values satisfying conditions
given by Eqs. 11 and 12, respectively.

###  2.3 Results of simulations

The most stable evolutionary obtained nanotubes are labeled as
$\tau_{n_{\phi},n_{z}}$. The computational complexity for a given value
$n_{\phi}$ grows exponentially with unit cell length $n_{z}$. For
$n_{\phi}=23$, under the assumption that a single structure is proceeded for
$10\,\mu{\rm s}$ even for $n_{z}=2$ with sequential computations, the time
five orders of magnitude longer than the Universe age [34] is needed to check
all possible potential solutions. The application of a genetic algorithm is
therefore advisable.

$n_{z}$ | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
---|---|---|---|---|---|---|---|---
|  |  |  |  |  |  |  |
| $\tau_{23,2}$ | $\tau_{23,3}$ | $\tau_{23,4}$ | $\tau_{23,5}$ | $\tau_{23,6}$ | $\tau_{23,7}$ | $\tau_{23,8}$ | $\tau_{23,9}$
$E_{\rm b}$ | $6{.}512$ | $6{.}547$ | $6{.}551$ | $6{.}550$ | $6{.}548$ | $6{.}550$ | $6{.}551$ | $6{.}548$
$c$ | 92 | 138 | 184 | 230 | 276 | 322 | 368 | 414
$N_{\rm a}$ | 92 | 122 | 164 | 203 | 247 | 285 | 327 | 370
$N_{\rm v}$ | 0 | 16 | 20 | 27 | 29 | 37 | 41 | 44
$n_{4}$ | 0 | 4 | 2 | 5 | 4 | 6 | 5 | 6
$n_{5}$ | 0 | 88 | 116 | 152 | 166 | 210 | 236 | 252
$n_{6}$ | 92 | 30 | 46 | 46 | 77 | 69 | 86 | 112
$d_{n_{\phi}}$ | 23 | 1 | 1 | 1 | 1 | 1 | 1 | 1
$d_{n_{z}}$ | 2 | 1 | 1 | 1 | 1 | 1 | 1 | 1
$\Delta t$ | 2:20:42:33 | 0:04:43:57 | 0:03:32:40 | 0:05:33:33 | 0:11:15:37 | 1:00:12:09 | 2:13:54:34 | 6:12:22:56
Figure 1: Fragments of boron sheets that if rolled into nanotubes give the
most stable BNTs that were obtained in our evolutionary computations. The unit
cell of each BNT is defined by two parameters $(n_{\phi},n_{z})$, where
$n_{\phi}=23$ and $n_{z}$ is specified above each stripe. Bellow each stripe
are listed several parameters that characterize each nanotube (see text). The
binding energy per atom, $E_{\rm b}$, is given in units of ${\rm eV}/{\rm
atom}$. The holes in the nanotubes are vacancies in the hexagonal structure
wrapped on a cylinder to form a nanotube. The density of holes can be
calculated using the formula $\eta=N_{\rm v}/c$, where $c=N_{\rm a}+N_{\rm
v}$. The values of the calculation time, $\Delta t$, are presented using the
days:hours:minutes:seconds format. One evolutionary process was performed on
one logical thread of a multithreaded CPU [35]. The unit cell of the
$\tau_{23,8}$ BNT is shown in the inset of Fig. 2. Figure 2: Binding energy
per atom, $E_{\rm b}$, as a function of hole density, $\eta$. The data points
in black correspond to BNTs with maximum values of $E_{\rm b}$ obtained for a
given $n_{z}$ parameter ($n_{\phi}=23$) and shown with a different label. The
data points in gray show the full evolution and correspond to nonoptimal BNTs.
The $n_{z}\in\\{3,4,\dots,9\\}$ to label correspondence is provided above the
data points. For comparison purposes, shown are also the values of $E_{\rm b}$
for the $hc$ and $\alpha$ sheets (with $E_{\rm b}$ equal to $e_{3}$ and
$6{.}603\,{\rm eV/atom}$, respectively [30]) as well as $E_{\rm b}$ values
taken from Ref. [30] corresponding to the BTC and ${\rm B}_{36}$ structures.
The inset on the left (with gray background) shows an enlarged view of a
fragment of the main picture. The histogram on the right presents (from
bottom) the values of $n_{5}/c$ (gray), $n_{6}/c$ (white), $n_{4}/c$ (black),
and $N_{\rm v}/c$ (white) for $\tau_{23,n_{z}}$ nanotubes. In the center, the
unit cell of the $\tau_{23,8}$ nanotube is shown.

The most energetically favorable nanotubes are presented in Fig. 1. The
evolutionary obtained structure $\tau_{23,2}$ is the case of a fully filed
with boron atoms lattice. The removal of any atom from a fully filled lattice
is energetically favorable since raises the number of 5-coordinated atoms
(raise of the beneficial value of $n_{5}$). The density of vacancies can not
be too high since a situation with holes that are too close to each other may
be detrimental for $E_{\rm b}$ since may increase the number of 4-coordinated
atoms (raise of the detrimental value of $n_{4}$). Looking at Fig. 1, we may
conclude that, in general, the most energetically favorable are those
nanotubes for which the boron holes are separated by fragments of the boron
double chain (BDC) or boron triple chain (BTC). Less favorable are those
situations in which the holes are next to each other. The highest $E_{\rm b}$
in our simulation belongs to the $\tau_{23,8}$ nanotube ($6{.}55108\,{\rm
eV}/{\rm atom}$) after which the $\tau_{23,4}$ nanotube ($6{.}55102\,{\rm
eV}/{\rm atom}$) is located. Recalling that DFT calculations precision is of
the order of $1\,{\rm meV}/{\rm atom}$ [36], one can assume that both
nanotubes have equal $E_{\rm b}$ (within our model). In Fig. 1, we also show
the time taken to obtain within our approach, the most stable examples of BNTs
for a given value of $n_{z}$. Although our methodology is relatively simple,
the time taken to obtain for instance $\tau_{23,9}$ exceeds 6 days. This tells
us that computations that would involve more accurate fitness functions may
require prohibitive computational times.

The relation between $E_{\rm b}$ and hole density, $\eta=N_{\rm v}/c$, is
presented in Fig. 2. If one were to draw a curve that is the average of the
$E_{\rm b}$ values shown in Fig. 2 as a function of hole density, the maximum
of this curve would correspond to the hole density of $1/9=0.111$ which is the
hole density of the $\alpha$-sheet [27] (the optimal flat neutral boron 2D
structure). The histogram shown in the inset of Fig. 2 presents the number of
4, 5, and 6 coordinated boron atoms as well as boron vacancies relative to the
total number of possible boron sites ($c=N_{\rm a}+N_{\rm v}$) for the most
stable BNTs. The average over the 7 relevant cases studied
($n_{z}\in\\{3,4,\dots,9\\}$) gives us that the number of 4-, 5-, and
6-coordinated boron atoms is about $2\%$, $63\%$, and $24\%$, respectively, of
the total number of atoms forming the nanotubes, whereas about $11\%$ are
boron vacancies. This gives us a very interesting result that more than $85\%$
of the boron atoms in BNTs with the largest $E_{\rm b}$ values are highly
coordinated atoms.

In Fig. 3, we show three examples of BNTs for three substantially different
hole densities $\eta$ equal to 0, 0.122, and 0.239. For hole densities smaller
than about 0.22, the boron nanotubes exhibit randomly distributed hexagonal
holes and are more stable than a flat stripe structure and a quasi-flat B36
cluster (see Fig. 2). For larger hole densities ($>0.22$) the boron nanotubes
resemble porous tubular structures with hole sizes that depend on the surface
densities of boron atoms. Finally should be noted that almost all the studied
BNTs in this work have $E_{\rm b}$ values between $5.65~{}{\rm eV}/{\rm atom}$
and $6.60~{}{\rm eV}/{\rm atom}$ which correspond to $E_{\rm b}$ of the $hc$
and $\alpha$ sheets, respectively [30].

Figure 3: Examples of BNTs with three different hole densities. For each
nanotube, the $n_{\phi}$ and $n_{z}$ parameters are provided. To obtain the
same height of the nanotubes the 2, 5, and 3 values of $n_{z}$ are multiplied
by 15, 6, and 10, respectively (for a total height of 50.25 Å). The binding
energy, $E_{\rm b}$, and hole density, $\eta$, is also given.

##  3 Summary

In conclusion, we have presented an adaptation of the genetic algorithm for
the structure prediction of tubular BNTs in a broad range of possible hole
densities. The evolution of the structures is done using a binary
representation and a fitness function calculated using a model Hamiltonian.
For the most stable BNTs, the number of 5-coordinated boron atoms is about
$63\%$ of the total number of atoms forming the nanotubes, whereas about
$11\%$ are boron vacancies. Moreover, more than $85\%$ of the boron atoms in
BNTs with the largest $E_{\rm b}$ values are 5- and 6-coordinated atoms. For
hole densities smaller than about 0.22, the BNTs exhibit randomly distributed
hexagonal holes and are more stable than the BTC structure and a quasi-flat
B36 cluster. For larger hole densities ($>0.22$) the BNTs resemble porous
tubular structures with hole sizes that depend on the surface densities of
boron atoms. These prototype genetic algorithm calculations may serve for the
exploration of the structure of nanotubes based on a mixture of boron and
other atoms (e.g. carbon and/or nitrogen).

## Appendix A

Let us define some helper functions describing the indices of six lattice
sites (atoms and vacancies) that are neighbors of a given site of the wrapped
hexagonal lattice. If a given lattice site has an index $i$, then the sites up
left, up right, left, right, down left, and down right have indices equal to
$u_{l}^{*}(i)$, $u_{r}^{*}(i)$, $l^{*}(i)$, $r^{*}(i)$, $d_{l}^{*}(i)$, and
$d_{r}^{*}(i)$, respectively. Those functions are formulated in the following
way:

$\displaystyle u_{l}^{*}$ $\displaystyle=l^{*}\circ u_{r}^{*},$ (13)
$\displaystyle u_{r}^{*}(i)$ $\displaystyle=(i+n_{z}+o^{*}(i))\bmod
2n_{\phi}n_{z},$ (14) $\displaystyle l^{*}(i)$
$\displaystyle=i-1+\left\\{\begin{array}[]{ll}0&{\rm if}\ i\bmod n_{z}\neq
0\\\ n_{z}&{\rm otherwise}\end{array}\right.,$ (17) $\displaystyle r^{*}(i)$
$\displaystyle=i+1-\left\\{\begin{array}[]{ll}0&{\rm if}\ i\bmod n_{z}\neq
n_{z}-1\\\ n_{z}&{\rm otherwise}\end{array}\right.,$ (20) $\displaystyle
d_{l}^{*}$ $\displaystyle=l^{*}\circ d_{r}^{*},$ (21) $\displaystyle
d_{r}^{*}(i)$ $\displaystyle=(i-n_{z}+o^{*}(i))\bmod 2n_{\phi}n_{z},$ (22)

where the offset $o^{*}$ is defined as:

$o^{*}(i)=\left\\{\begin{array}[]{ll}0&{\rm if}\ \lfloor i/n_{z}\rfloor\ {\rm
is\ even}\\\ 1&{\rm if}\ \lfloor i/n_{z}\rfloor\ {\rm is\ odd}\wedge i\bmod
n_{z}\neq n_{z}-1\\\ 1-n_{z}&{\rm otherwise}\end{array}\right..$ (23)

Let ${\rm id}(i)=i$. The following relations hold:

$\displaystyle u_{l}^{*}$ $\displaystyle={d_{r}^{*}}^{-1},$ (24)
$\displaystyle l^{*}$ $\displaystyle={r^{*}}^{-1},$ (25) $\displaystyle
d_{l}^{*}$ $\displaystyle={u_{r}^{*}}^{-1},$ (26)
$\displaystyle{u_{l}^{*}}^{2n_{\phi}}$
$\displaystyle={d_{l}^{*}}^{2n_{\phi}}={l^{*}}^{2},$ (27)
$\displaystyle{l^{*}}^{n_{z}}$ $\displaystyle={r^{*}}^{n_{z}}={\rm id},$ (28)
$\displaystyle{u_{r}^{*}}^{2n_{\phi}}$
$\displaystyle={d_{r}^{*}}^{2n_{\phi}}={r^{*}}^{2}.$ (29)

Functions $u_{l}^{*}$, $u_{r}^{*}$, $r^{*}$, and $d_{r}^{*}$ are used to
ensure neighborhood along the nanotube axis direction and at its
circumference. To ensure connectivity, an additional set of functions is
needed. These functions do not reflect PBCs. If a given atom or vacancy does
not possess a suitable neighbor (i.e., it lies at the unit cell boundary),
then the appropriate function returns an unchanged index:

$\displaystyle u_{l}(i)$ $\displaystyle=i+\left\\{\begin{array}[]{ll}0&{\rm
if}\ \lfloor i/n_{z}\rfloor=2n_{\phi}-1\\\ 0&{\rm if}\ \lfloor i/n_{z}\rfloor\
{\rm is\ even}\wedge i\bmod n_{z}=0\\\ n_{z}-1\hphantom{()0}&{\rm if}\ \lfloor
i/n_{z}\rfloor\ {\rm is\ even}\wedge i\bmod n_{z}\neq 0\\\ n_{z}&{\rm
otherwise}\end{array}\right.,$ (34) $\displaystyle u_{r}(i)$
$\displaystyle=i+\left\\{\begin{array}[]{ll}0&{\rm if}\ \lfloor
i/n_{z}\rfloor=2n_{\phi}-1\\\ 0&{\rm if}\ \lfloor i/n_{z}\rfloor\ {\rm
odd}\wedge i\bmod n_{z}=n_{z}-1\\\ n_{z}+1\hphantom{()0}&{\rm if}\ \lfloor
i/n_{z}\rfloor\ {\rm odd}\wedge i\bmod n_{z}\neq n_{z}-1\\\ n_{z}&{\rm
otherwise}\end{array}\right.,$ (39) $\displaystyle l(i)$
$\displaystyle=i-\left\\{\begin{array}[]{ll}1&{\rm if}\ i\bmod n_{z}\neq 0\\\
0\hphantom{(n_{z}+1)}&{\rm otherwise}\end{array}\right.,$ (42) $\displaystyle
r(i)$ $\displaystyle=i+\left\\{\begin{array}[]{ll}1&{\rm if}\ i\bmod n_{z}\neq
n_{z}-1\\\ 0\hphantom{(n_{z}+1)}&{\rm otherwise}\end{array}\right.,$ (45)
$\displaystyle d_{l}(i)$ $\displaystyle=i-\left\\{\begin{array}[]{ll}0&{\rm
if}\ \lfloor i/n_{z}\rfloor=0\\\ 0&{\rm if}\ \lfloor i/n_{z}\rfloor\ {\rm is\
even}\wedge i\bmod n_{z}=0\\\ (n_{z}+1)\hphantom{0}&{\rm if}\ \lfloor
i/n_{z}\rfloor\ {\rm is\ even}\wedge i\bmod n_{z}\neq 0\\\ n_{z}&{\rm
otherwise}\end{array}\right.,$ (50) $\displaystyle d_{r}(i)$
$\displaystyle=i-\left\\{\begin{array}[]{ll}0&{\rm if}\ \lfloor
i/n_{z}\rfloor=0\\\ 0&{\rm if}\ \lfloor i/n_{z}\rfloor\ {\rm is\ odd}\wedge
i\bmod n_{z}=n_{z}-1\\\ (n_{z}-1)\hphantom{0}&{\rm if}\ \lfloor
i/n_{z}\rfloor\ {\rm is\ odd}\wedge i\bmod n_{z}\neq n_{z}-1\\\ n_{z}&{\rm
otherwise}\end{array}\right..$ (55)

To formulate the $Q$ predicate, elements of graph theory is used. The set of
vertices, $V$, consists of indices of all atoms taken from a single unit cell:

$V=\left\\{i\in\iota_{2n_{\phi}n_{z}}\mid x_{i}\right\\},$ (56)

while the set of (undirected) edges, $E$, which vertices fulfill the
neighborhood condition in the nanotube’s fragment (unit cell) without PBCs, is
now defined as:

$E=\left\\{\left\\{i,j\right\\}\subseteq V\mid
j\in\left\\{u_{l}(i),u_{r}(i),l(i),r(i),d_{l}(i),d_{r}(i)\right\\}\setminus\left\\{i\right\\}\right\\}.$
(57)

By taking a pair of both objects one obtains graph $H=\left(V,E\right)$. The
$Q$ predicate is defined as $Q\equiv Q_{0}\wedge Q_{1}\wedge Q_{2}$, where:

$\displaystyle Q_{0}\equiv$ $\displaystyle\,\omega(H)=1,$ (58) $\displaystyle
Q_{1}\equiv$ $\displaystyle\left(\exists
i\in\iota_{2n_{\phi}n_{z}}\colon\lfloor i/n_{z}\rfloor\ {\rm is\ odd}\wedge
i\bmod n_{z}=n_{z}-1\wedge x_{i}\wedge(x_{u_{r}^{*}(i)}\vee
x_{d_{r}^{*}(i)})\right)\vee$ $\displaystyle\left(\exists
i\in\iota_{2n_{\phi}n_{z}}\colon\lfloor i/n_{z}\rfloor\ {\rm is\ even}\wedge
i\bmod n_{z}=n_{z}-1\wedge x_{i}\wedge x_{r^{*}(i)}\right),$ (59)
$\displaystyle Q_{2}\equiv$ $\displaystyle\left(\exists
i\in\iota_{2n_{\phi}n_{z}}\colon\lfloor i/n_{z}\rfloor=2n_{\phi}-1\wedge
x_{i}\wedge(x_{u_{l}^{*}(i)}\vee x_{u_{r}^{*}(i)})\right),$ (60)

and $\omega(H)$ is the number of connected components of graph $H$.

## Acknowledgments

This work is a result of the projects funded by the National Science Centre of
Poland (Twardowskiego 16, PL-30312 Kraków, Poland, http://www.ncn.gov.pl/)
under grants number UMO-2013/11/B/ST3/04273 and UMO-2016/23/B/ST3/03575.

## References

  * [1] D. Ciuparu, R. F. Klie, Y. Zhu, L. Pfefferle, Synthesis of pure boron single-wall nanotubes, The Journal of Physical Chemistry B 108 (13) (2004) 3967–3969.  doi:10.1021/jp049301b.
URL https://doi.org/10.1021/jp049301b

  * [2] F. Liu, C. Shen, Z. Su, X. Ding, S. Deng, J. Chen, N. Xu, H. Gao, Metal-like single crystalline boron nanotubes: synthesis and in situ study on electric transport and field emission properties, Journal of Materials Chemistry 20 (11) (2010) 2197.  doi:10.1039/b919260c.
URL https://doi.org/10.1039/b919260c

  * [3] J. Tian, Z. Xu, C. Shen, F. Liu, N. Xu, H.-J. Gao, One-dimensional boron nanostructures: Prediction, synthesis, characterizations, and applications, Nanoscale 2 (8) (2010) 1375.  doi:10.1039/c0nr00051e.
URL https://doi.org/10.1039%2Fc0nr00051e

  * [4] Y. Tian, Z. Guo, T. Zhang, H. Lin, Z. Li, J. Chen, S. Deng, F. Liu, Inorganic boron-based nanostructures: Synthesis, optoelectronic properties, and prospective applications, Nanomaterials 9 (4) (2019) 538.  doi:10.3390/nano9040538.
URL https://doi.org/10.3390%2Fnano9040538

  * [5] I. Boustani, New convex and spherical structures of bare boron clusters, Journal of Solid State Chemistry 133 (1) (1997) 182–189.  doi:10.1006/jssc.1997.7424.
URL https://doi.org/10.1006%2Fjssc.1997.7424

  * [6] I. Boustani, Systematic ab initio investigation of bare boron clusters: Determination of the geometry and electronic structures of $B_{n}~{}(n=2\mathrm{-}14)$, Physical Review B 55 (24) (1997) 16426–16438.  doi:10.1103/physrevb.55.16426.
URL https://doi.org/10.1103/physrevb.55.16426

  * [7] I. Boustani, A. Quandt, Nanotubules of bare boron clusters: Ab initio and density functional study, Europhysics Letters (EPL) 39 (5) (1997) 527–532.  doi:10.1209/epl/i1997-00388-9.
URL https://doi.org/10.1209%2Fepl%2Fi1997-00388-9

  * [8] A. Gindulytė, W. N. Lipscomb, L. Massa, Proposed Boron Nanotubes, Inorganic Chemistry 37 (1998) 6544–6545.  doi:10.1021/ic980559o.
  * [9] I. Boustani, A. Rubio, J. A. Alonso, Ab initio study of $b_{32}$ clusters: competition between spherical, quasiplanar and tubular isomers, Chemical Physics Letters 311 (1-2) (1999) 21–28.  doi:10.1016/s0009-2614(99)00767-8.
URL https://doi.org/10.1016%2Fs0009-2614%2899%2900767-8

  * [10] B. Kiran, S. Bulusu, H.-J. Zhai, S. Yoo, X. C. Zeng, L.-S. Wang, Planar-to-tubular structural transition in boron clusters: $b_{20}$ as the embryo of single-walled boron nanotubes, Proceedings of the National Academy of Sciences 102 (4) (2005) 961–964.  doi:10.1073/pnas.0408132102.
URL https://doi.org/10.1073/pnas.0408132102

  * [11] Z. A. Piazza, H.-S. Hu, W.-L. Li, Y.-F. Zhao, J. Li, L.-S. Wang, Planar hexagonal $\textnormal{B}_{36}$ as a potential basis for extended single-atom layer boron sheets, Nature Communications 5 (2014) 3113.  doi:10.1038/ncomms4113.
  * [12] N. Gonzalez Szwacki, A. Sadrzadeh, B. I. Yakobson, ${\mathrm{B}}_{80}$ Fullerene: An Ab Initio Prediction of Geometry, Stability, and Electronic Structure, Physical Review Letters 98 (2007) 166804.  doi:10.1103/PhysRevLett.98.166804.
  * [13] H.-J. Zhai, Y.-F. Zhao, W.-L. Li, Q. Chen, H. Bai, H.-S. Hu, Z. A. Piazza, W.-J. Tian, H.-G. Lu, Y.-B. Wu, Y.-W. Mu, G.-F. Wei, Z.-P. Liu, J. Li, S.-D. Li, L.-S. Wang, Observation of an all-boron fullerene, Nature Chemistry 6 (2014) 727–731.  doi:10.1038/nchem.1999.
  * [14] N. Gonzalez Szwacki, Boron Fullerenes: A First-Principles Study, Nanoscale Research Letters 3 (2007) 49.  doi:10.1007/s11671-007-9113-1.
  * [15] L. A. Chernozatonskii, P. B. Sorokin, B. I. Yakobson, New boron barrelenes and tubulenes, JETP Letters 87 (9) (2008) 489–493.  doi:10.1134/s0021364008090087.
URL https://doi.org/10.1134/s0021364008090087

  * [16] A. K. Singh, A. Sadrzadeh, B. I. Yakobson, Probing Properties of Boron $\alpha$-Tubes by Ab Initio Calculations, Nano Letters 8 (5) (2008) 1314–1317.  doi:10.1021/nl073295o.
URL https://doi.org/10.1021/nl073295o

  * [17] X. Yang, Y. Ding, J. Ni, Ab initio prediction of stable boron sheets and boron nanotubes: Structure, stability, and electronic properties, Physical Review B 77 (4) (jan 2008).  doi:10.1103/physrevb.77.041402.
URL https://doi.org/10.1103/physrevb.77.041402

  * [18] H. Tang, S. Ismail-Beigi, First-principles study of boron sheets and nanotubes, Physical Review B 82 (11) (sep 2010).  doi:10.1103/physrevb.82.115412.
URL https://doi.org/10.1103/physrevb.82.115412

  * [19] N. Gonzalez Szwacki, C. Tymczak, The symmetry of the boron buckyball and a related boron nanotube, Chemical Physics Letters 494 (1-3) (2010) 80–83.  doi:10.1016/j.cplett.2010.05.086.
URL https://doi.org/10.1016/j.cplett.2010.05.086

  * [20] R. N. Gunasinghe, C. B. Kah, K. D. Quarles, X.-Q. Wang, Dispersion corrections in the boron buckyball and nanotubes, Applied Physics Letters 98 (26) (2011) 261906.  doi:10.1063/1.3604018.
URL https://doi.org/10.1063/1.3604018

  * [21] J. Kunstmann, V. Bezugly, H. Rabbel, M. H. Rümmeli, G. Cuniberti, Unveiling the Atomic Structure of Single-Wall Boron Nanotubes, Advanced Functional Materials 24 (26) (2014) 4127–4134.  doi:10.1002/adfm.201304146.
  * [22] E. Caliskan, M. Kirca, Tensile characteristics of boron nanotubes by using reactive molecular dynamics simulations, Computational Materials Science 209 (2022) 111368.  doi:10.1016/j.commatsci.2022.111368.
URL https://doi.org/10.1016/j.commatsci.2022.111368

  * [23] T. Kondo, Recent progress in boron nanomaterials, Science and Technology of Advanced Materials 18 (1) (2017) 780–804.  doi:10.1080/14686996.2017.1379856.
URL https://doi.org/10.1080/14686996.2017.1379856

  * [24] T. M. Dieb, S. Ju, J. Shiomi, K. Tsuda, Monte carlo tree search for materials design and discovery, MRS Communications 9 (2) (2019) 532–536.  doi:10.1557/mrc.2019.40.
URL https://doi.org/10.1557/mrc.2019.40

  * [25] W. Paszkowicz, Genetic algorithms, a nature-inspired tool: Survey of applications in materials science and related fields, Materials and Manufacturing Processes 24 (2) (2009) 174–197.  doi:10.1080/10426910802612270.
URL https://doi.org/10.1080/10426910802612270

  * [26] T. E. B. Davies, M. M. Popa, C. V. Ciobanu, G. Maroulis, T. E. Simos, Genetic algorithm for structural optimization of tubular nanostructures, in: AIP Conference Proceedings, AIP, 2007.  doi:10.1063/1.2827001.
URL https://doi.org/10.1063/1.2827001

  * [27] H. Tang, S. Ismail-Beigi, Novel Precursors for Boron Nanotubes: The Competition of Two-Center and Three-Center Bonding in Boron Sheets, Physical Review Letters 99 (2007) 115501.  doi:10.1103/PhysRevLett.99.115501.
  * [28] T. Tarkowski, N. Gonzalez Szwacki, The structure of thin boron nanowires predicted using evolutionary computations (2022).  doi:10.48550/ARXIV.2211.11901.
URL https://arxiv.org/abs/2211.11901

  * [29] T. Tarkowski, Quilë (apr 2022).  doi:10.5281/zenodo.6484743.
URL https://zenodo.org/record/6484447

  * [30] T. Tarkowski, N. Gonzalez Szwacki, M. Marchwiany, Structure of porous two-dimensional boron crystals, Physical Review B 104 (2021) 195423.  doi:10.1103/PhysRevB.104.195423.
  * [31] T. Tarkowski, J. A. Majewski, N. Gonzalez Szwacki, Energy decomposition analysis of neutral and negatively charged borophenes, FlatChem 7 (2018) 42–47.  doi:10.1016/j.flatc.2017.08.004.
  * [32] E. Hernández, C. Goze, P. Bernier, A. Rubio, Elastic properties of single-wall nanotubes, Applied Physics A 68 (1999) 287–292.  doi:10.1007/s003390050890.
  * [33] P. Umari, O. Petrenko, S. Taioli, M. M. D. Souza, Communication: Electronic band gaps of semiconducting zig-zag carbon nanotubes from many-body perturbation theory calculations, The Journal of Chemical Physics 136 (18) (2012) 181101.  doi:10.1063/1.4716178.
URL https://doi.org/10.1063/1.4716178

  * [34] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, S. Basak, R. Battye, K. Benabed, J.-P. Bernard, M. Bersanelli, P. Bielewicz, J. J. Bock, J. R. Bond, J. Borrill, F. R. Bouchet, F. Boulanger, M. Bucher, C. Burigana, R. C. Butler, E. Calabrese, J.-F. Cardoso, J. Carron, A. Challinor, H. C. Chiang, J. Chluba, L. P. L. Colombo, C. Combet, D. Contreras, B. P. Crill, F. Cuttaia, P. de Bernardis, G. de Zotti, J. Delabrouille, J.-M. Delouis, E. Di Valentino, J. M. Diego, O. Doré, M. Douspis, A. Ducout, X. Dupac, S. Dusini, G. Efstathiou, F. Elsner, T. A. Enßlin, H. K. Eriksen, Y. Fantaye, M. Farhang, J. Fergusson, R. Fernandez-Cobos, F. Finelli, F. Forastieri, M. Frailis, A. A. Fraisse, E. Franceschi, A. Frolov, S. Galeotta, S. Galli, K. Ganga, R. T. Génova-Santos, M. Gerbino, T. Ghosh, J. González-Nuevo, K. M. Górski, S. Gratton, A. Gruppuso, J. E. Gudmundsson, J. Hamann, W. Handley, F. K. Hansen, D. Herranz, S. R. Hildebrandt, E. Hivon, Z. Huang, A. H. Jaffe, W. C. Jones, A. Karakci, E. Keihänen, R. Keskitalo, K. Kiiveri, J. Kim, T. S. Kisner, L. Knox, N. Krachmalnicoff, M. Kunz, H. Kurki-Suonio, G. Lagache, J.-M. Lamarre, A. Lasenby, M. Lattanzi, C. R. Lawrence, M. Le Jeune, P. Lemos, J. Lesgourgues, F. Levrier, A. Lewis, M. Liguori, P. B. Lilje, M. Lilley, V. Lindholm, M. López-Caniego, P. M. Lubin, Y.-Z. Ma, J. F. Macías-Pérez, G. Maggio, D. Maino, N. Mandolesi, A. Mangilli, A. Marcos-Caballero, M. Maris, P. G. Martin, M. Martinelli, E. Martínez-González, S. Matarrese, N. Mauri, J. D. McEwen, P. R. Meinhold, A. Melchiorri, A. Mennella, M. Migliaccio, M. Millea, S. Mitra, M.-A. Miville-Deschênes, D. Molinari, L. Montier, G. Morgante, A. Moss, P. Natoli, H. U. Nørgaard-Nielsen, L. Pagano, D. Paoletti, B. Partridge, G. Patanchon, H. V. Peiris, F. Perrotta, V. Pettorino, F. Piacentini, L. Polastri, G. Polenta, J.-L. Puget, J. P. Rachen, M. Reinecke, M. Remazeilles, A. Renzi, G. Rocha, C. Rosset, G. Roudier, J. A. Rubiño-Martín, B. Ruiz-Granados, L. Salvati, M. Sandri, M. Savelainen, D. Scott, E. P. S. Shellard, C. Sirignano, G. Sirri, L. D. Spencer, R. Sunyaev, A.-S. Suur-Uski, J. A. Tauber, D. Tavagnacco, M. Tenti, L. Toffolatti, M. Tomasi, T. Trombetti, L. Valenziano, J. Valiviita, B. Van Tent, L. Vibert, P. Vielva, F. Villa, N. Vittorio, B. D. Wandelt, I. K. Wehus, M. White, S. D. M. White, A. Zacchei, A. Zonca, Planck 2018 results - VI. Cosmological parameters, Astronomy & Astrophysics 641 (2020) A6.  doi:10.1051/0004-6361/201833910.
  * [35] [link].
URL https://www.amd.com/en/product/8436

  * [36] R. Hoffmann, P. von Ragué Schleyer, H. F. Schaefer III, Predicting Molecules–More Realism, Please!, Angewandte Chemie International Edition 47 (38) (2008) 7164–7167.  doi:10.1002/anie.200801206.