# Partitioning and Placement of Deep Neural Networks on Distributed Edge
Devices to Maximize Inference Throughput

Arjun Parthasarathy  Crystal Springs Uplands School
Email<EMAIL_ADDRESS>Bhaskar Krishnamachari  University of
Southern California
Email<EMAIL_ADDRESS>
(November 2022)

###### Abstract

Edge inference has become more widespread, as its diverse applications range
from retail to wearable technology. Clusters of networked resource-constrained
edge devices are becoming common, yet no system exists to split a DNN across
these clusters while maximizing the inference throughput of the system. We
present an algorithm which partitions DNNs and distributes them across a set
of edge devices with the goal of minimizing the bottleneck latency and
therefore maximizing inference throughput. The system scales well to systems
of different node memory capacities and numbers of nodes. We find that we can
reduce the bottleneck latency by 10x over a random algorithm and 35% over a
greedy joint partitioning-placement algorithm. Furthermore we find empirically
that for the set of representative models we tested, the algorithm produces
results within 9.2% of the optimal bottleneck latency.

##  I Introduction

Deep Neural Networks (DNNs) have greatly accelerated machine learning across
different disciplines, such as Computer Vision [5] and Natural Language
Processing [21]. Edge Inference is becoming an increasingly popular field with
multiple facets [33], as sensor-driven computation in IoT systems necessitates
DNN inference in the field. IoT applications for edge inference range from
retail to wearable technology [6, 4].

The edge can come in multiple configurations [19, 28], and there are multiple
approaches to facilitate edge inference. For cloud-edge hybrid inference, one
such approach is model compression [11], which deals exclusively with DNN
optimization but does not address the system’s runtime configuration. In this
paper, we focus on clusters of resource-constrained edge devices. These edge
clusters are becoming increasingly common due to their low-cost and
scalability at the edge [24]. Unlike a cloud data center, the edge brings
system resource limitations and communication bottlenecks between devices.

With this in mind, we address the following problem: How can we take advantage
of multi-device edge clusters to enable high-performance DNN inference while
respecting computational resource constraints and taking into account the
heterogeneity of communication links?

To partition a deep learning model, we first split the model into components
that are executed sequentially. Each partition is assigned to a different edge
device, and once each node performs inference with its piece of the model,
that intermediate inference result is sent to the next node with the
corresponding partition in the sequence. This inference pipeline is shown in
Figure 1.

Figure 1: Partitioning and Distributing a Model Across Edge Devices to Create
an Inference Pipeline

In an edge cluster, although we have a lower computational power in each node,
we can take advantage of this inference pipelining to increase system
throughput. Since each node can perform inference with its partition
individually, prior nodes in the pipeline can send their finished inference
results to the subsequent nodes in the pipeline and accept new batches.

We define the throughput metric of a system as the number of inference cycles
it can perform per unit time. As we showed in our previous work DEFER [23], we
can achieve higher throughput with distributed edge inference as opposed to
inference on a single device because of pipelining. The throughput is defined
as the reciprocal of the bottleneck latency. For nodes
$[k]=\\{1,2,\dots,k\\}$, the bottleneck latency $\beta$ is defined as

$\begin{array}[]{c}S=\\{k\in[k]\mid c_{k},\gamma_{k}\\}\\\ \beta=\max_{s\in
S}{s}\end{array}$ (1)

where $c_{k}$ is the compute time of the operations on node $k$, and
$\gamma_{k}$ is the communication time between node $k-1$ and $k$. We use
ResNet50 [12], which is a representative model for our use case. On a
Raspberry Pi 4, the inference speed was found to be 225 ms [25]. Next, we
found the amount of data transferred between each layer of the model. On
average, 10.2 Mbits of data was transferred between layers. Given an average
WiFi bandwidth of 6 Mbps for a low-end edge network, this gives us a
communication time of $1.7s$. This is 7.5x slower than the compute time. In
reality, many models are larger than ResNet50 and will therefore be split
across devices, so each device will have less operations to execute. This
means that communication time will outweigh compute time as the bottleneck.
Therefore, we can simplify the expression for bottleneck latency to the
following:

$\beta=\max_{k\in[k]}{\gamma_{k}}$ (2)

Since throughput is defined as $\frac{1}{\beta}$, by minimizing the bottleneck
latency we maximize inference throughput. Additionally, we assume that all
nodes are homogeneous in RAM. If the devices are not the same capacity, then
the algorithm will take the smallest memory capacity across all nodes in the
cluster, and take that as the capacity of each node. In this paper, we
primarily analyze image and text models due to their prevalence on the edge
for visual analytics applications [34, 22].

Our main contribution in this paper is a novel partitioning and placement
algorithm for DNNs across a cluster of edge devices distributed spatially
within the same WiFi network. The algorithm finds the candidate partition
points, finds the optimal partition sizes to transfer the least amount of
data, and finds the arrangement of nodes with the highest bandwidth. Together,
these aim to minimize the resulting bottleneck latency according to the
throughput metric. We found that our algorithm results in a 10x improvement
over a random partitioning/placement algorithm, and a 35% reduction in
bottleneck latency for systems with 50 compute nodes. We empirically observe
an average approximation ratio of 1.092 for the bottleneck latency (i.e. it is
9.2% more than the optimal bottleneck latency, on average).

##  II Related Work

Early works on the topic of partitioning DNN models divided them into head and
tail models with the former distilled to enable running on a resource-
constrained device and reduce data transfer [20]. Some prior works on DNN edge
inference mathematically perform DNN model slicing by layer [35, 36], after
calculating layer impact during the training stage; these do not account for
communication demands on the edge. Others abstract model layers into certain
“execution units,” [17, 7] which they then choose to slice based on certain
resource requirements. Li _et al._ [16] regressively predict a layer’s latency
demand and optimize communication bandwidth accordingly. DeeperThings [29]
performs layer fusion on CNNs to optimize data transfer. These works are
optimized for a hybrid edge-cloud pipeline and do not address the demands of a
cluster of edge devices. Couper [13] uses a similar partitioning scheme to
minimize inter-partition data transfer, but does not address the communication
bottleneck associated with an edge cluster. Hu _et. al_ [14] optimize the
partitioning of a CNN onto a set of devices by taking compute time as a
bottleneck, while employing compression to deal with communication
constraints, and do not consider placement. Our paper builds on and
differentiates itself from these works by addressing the bandwidth limitation
of an edge cluster, and aims to maximize inter-node bandwidth during the
placement stage to minimize bottleneck latency.

##  III Partitioning and Placement Algorithm

We are given two graphs:

  1. 1.

An unweighted DAG $G_{m}$ representing the computation graph of a DNN, where
each vertex represents a layer in the model. This DAG can be found using
common ML libraries such as Tensorflow [1] and Keras [8].

  2. 2.

A weighted complete graph $G_{c}$ representing the communication graph of a
cluster of homogeneous physical compute nodes, where each vertex represents a
physical compute node and each edge represents the bandwidth between those
nodes. The graph is complete because we assume that these edge devices will
communicate over the same WiFi network.

Our goal is to optimally partition the model and place these partitions on a
set of edge devices. We do so as follows.

###  III-A Converting a Complex DAG to a Linear DAG

First, we need to distill $G_{m}$ into a linear DAG. The vertices where it is
possible to partition the model are called “candidate partition points.” We
illustrate this in Figure 2.

For $v\in V$, edges $e\in E$ and source vertex $s$ of $G_{m}$, find the
longest path from $s$ to $v$. This can be done by topologically sorting the
DAG and for each vertex in the resulting list, relaxing each neighbor of that
vertex. We call the length of this longest path the topological depth of that
vertex in the graph. Let $LP(v)$ denote the length of longest path from $s$ to
$v$.

To verify that all paths from vertex $v_{prev}$ go through vertex $v$, use a
modified DFS by recursing on the incident edges of each vertex. If we
encounter a vertex with a greater topological depth than $v$, return false. If
we reach vertex $v$, return true. Let $AP(v_{prev},v)$ denote the result of
this algorithm.

Given the previously found candidate partition point $p_{k-1}$ and the current
vertex $u$, the next candidate partition point $p_{k}=u$ iff:

  1. 1.

$LP(u)\neq LP(v)\forall v\in\\{V-u\\}$

  2. 2.

$AP(p_{k-1},u)=\text{true}$

with $p_{0}=s$.

The time complexity of LP is $O(V+E)$. AP runs in polynomial time by returning
upon reaching a vertex with a greater topological depth. Therefore, this
algorithm runs in polynomial time.

Figure 2: Partition points for ResNet50 and InceptionResNetV2 models

Figure 2 shows the candidate partition points at certain sections of the DAG
of ResNet50 [12] and InceptionResNetV2 [30]. Each rectangle represents a model
layer in the DAG.

Figure 3: Histogram of Number of Candidate Partition Points

We then calculated the number of partition points for the set of Keras
pretrained models. As shown in Figure 3, almost all the models have at least
25 candidate partition points. There are some model architectures, like NASNet
[37], which do not allow partitioning under our scheme.

Figure 4: Portion of NASNet’s layer DAG

As shown in Figure 4, NASNet cannot be partitioned because there is no single
point that splits the model into a distinct execution unit that does not have
any dependencies to a previous or subsequent layer. If we run our LP
algorithm, we find that there is no single layer that has distinct topological
depth from other layers. We found that 64 of the 66 (97 %) pretrained Keras
models [15] could be partitioned under our scheme, and only the NASNet
variants could not.

###  III-B Optimal model partitioning and placement

Our goal is to maximize throughput of the system. As previously discussed,
this means we need to minimize the bottleneck latency. Latency is defined as
$\frac{\text{data}}{\text{bandwidth}}$. Given a tuple of partition points
$P_{opt}$, their transfer sizes $T$, and a set of bandwidths $B$ between
compute nodes, the latency between each set of compute node is defined as

$\gamma_{k}=\frac{T_{opt,k}}{B_{k}}\forall 0\leq k<\lvert P_{opt}\rvert$ (3)

The bottleneck latency for the system is then given by Equation 2. For the
purposes of explanation, we separate the problems of optimizing the partitions
(thereby optimizing transfer size) and optimizing placement (thereby
optimizing bandwidth between nodes). We show empirically that this results in
the the smallest bottleneck latency. In Section V, we compare this formulation
to an algorithm that tries to jointly optimize transfer size and bandwidth.

####  III-B1 Finding optimal partitions

Our heuristic for finding optimal partitions is the “transfer size” of the
partition; i.e how much data will be transferred from that partition to the
next. Given the tuple of candidate partition points
$P=(p_{0},p_{1},\dots,p_{k})$, we now need to find a set of model partitions
which minimizes the sum of transfer sizes. Assuming a batch size of 1, the
transfer size $t_{k}$ of candidate partition point $p_{k}$ is defined as

$t_{k}=\frac{\eta}{\lambda}$ (4)

Given a floating point array representing the output of the model layer
$l_{k}$ (which is the same layer given by the candidate partition point
$p_{k}$), $\eta$ represents the size of that array.

$\lambda\approx 1.44*2.1$ represents the total compression ratio given by
multiplying the average ZFP compression ratio [18] by the average LZ4
compression ratio [9].

To better illustrate our algorithm, we classify the transfer size $t_{k}$ into
3 transfer size classes (“low”, “medium”, or “high”) based on the distribution
of the transfer sizes.

$\begin{array}[]{lr}C=\\{L,M,H\\}&t_{k}\subseteq C\end{array}$ (5)

The optimal set of partitions is the scheme which minimizes the sum of the
transfer sizes of said partitions.
Let $G_{p}$ represent a DAG, where each vertex is represented by a possible
partition. The vertices are defined as follows:

$\begin{array}[]{l}p_{i},p_{i+1},\dots,p_{j}\\})<\kappa\mid\\{p_{i},p_{i+1},\dots,p_{j}\\}\\}\\\
\forall{0\leq i<\lvert P_{opt}\rvert,0\leq j<\lvert
P_{opt}\rvert-i}\end{array}$ (6)

The set of vertices represents every possible contiguous subarray of candidate
partition points, where $\omega(P)$ finds whether the memory use of partition
$P$ is within the memory capacity $\kappa$ of the compute node. We quantize
the models using TFLite [32] quantization to reduce their memory footprint.
However, when calculating the memory footprint of a partition, we do not
consider this quantization. This means that we are conservative on partition
size and in turn provide extra space on each device for the memory overhead
from containerization. Each partition is a set of layers that fall between the
partition points $p_{i}$ and $p_{j}$.

The set of edges is defined as follows:

$E=\\{(u,v)\in V,\rho(u_{\lvert u\rvert-1})=\rho(v_{0}-1)\mid(u,v)\\}$ (7)

The function $\rho(\upsilon)$ finds the index of element $\upsilon$ in
$P_{opt}$. There is an edge between vertices if the last partition point of
$u$’s partition is adjacent in $P$ to the first partition point of $v$’s
partition. For example, if $u=[1,2]$, $v=[3,4]$, and $P=(1,2,3,4)$, then
$(u,v)$ is an edge. Each edge has a weight $w(u,v)$ which corresponds to its
transfer size class.

Figure 5: Example partition graph, where the partition points are
$P=\\{1,2,3,4,5\\}$

Figure 5 shows an example partition graph, where edges that are the same color
will have the same weight. In the figure, “root” vertices have in-degree 0,
“leaf” vertices have out-degree 0, and “intermediate” vertices have neither.

Algorithm 1 Optimal Partitioning

// Map to store memoized paths

$pathFrom\leftarrow\textsc{NEW-MAP}($)

procedure MIN-COST-PATH($G$, $v$)

if $v.children=\emptyset$ then

return $v,0$

end if

$partitionLastLayer\leftarrow v[v.length-1]$

if $partitionLastLayer\notin pathFrom$ then

$paths\leftarrow[]$

for $c\in v.children$ do

$path,cost\leftarrow\textsc{MIN-COST-PATH}(G,c)$

$paths\leftarrow\textsc{APPEND}(paths,(path,cost))$

end for

$pathFrom[partitionLastLayer]=\textsc{MIN}(paths)$

end if

$best\leftarrow pathFrom[partitionLastLayer]$

$minPath,minCost\leftarrow best$

$chosenNode\leftarrow minPath[0]$

// Path starting at v and going to a leaf

$newPath\leftarrow\textsc{APPEND}([v],...minPath)$

$newCost\leftarrow minCost+w(v,chosenNode)$

return $newPath,newCost$

end procedure

procedure PARTITION(G)

$roots\leftarrow\textsc{GET-ROOT-VERTICES}(G)$

for $r\in roots$ do

$path,cost\leftarrow\textsc{MIN-COST-PATH}(G,r)$

$paths\leftarrow\textsc{APPEND}(paths,(path,cost))$

end for

$minPath,minCost\leftarrow\textsc{MIN}(paths)$

return minPath

end procedure

$\Theta\leftarrow\textsc{PARTITION}(G_{p})$

Algorithm 1 finds the shortest path in the graph from a root to a leaf. Since
edges which bridge the same candidate partition points (and have the same
color as shown in Figure 5) will have the same subsequent paths, we can
memoize the shortest path. On line 2, we store a map on which tells us for
each candidate partition point what the shortest path is from that point.
Using memoization, Algorithm 1 takes $O(N)$ to find the shortest path, but
$O(N^{2})$ to construct the partition graph. Therefore the runtime of
Algorithm 1 is $O(N^{2})$, where N is the number of nodes.

Let $\Theta$ represent the set of chosen partitions. For each subarray in
$\Theta$, we take the last element of the subarray, add that to the list of
partition points $Q$, and add its corresponding transfer size to the list $S$.
The resulting list is then sorted based on the topological depth of each
partition point, so that the partitions are executed in the order they appear
in the model.

####  III-B2 Finding optimal model placement

With the set of optimal partitions $Q$ and their corresponding transfer sizes
$S$, we now need to “match” them to the vertices of $G_{c}$. We know from
Equation 5 that $S\subseteq C$. Let $c(e)$ return the bandwidth class of a
given edge of $G_{c}$. We use the following threshold function to classify
each edge:

$\tau(X,t)=\left\\{\begin{array}[]{lr}c(e)=C_{\arg X-1},&\text{if }e<t\\\
c(e)=X,&\text{if }e\geq t\end{array}\right\\}{\forall e\in E_{c}}$ (8)

If the edge is greater than or equal to the threshold, it will be classified
as class $X$, otherwise it will be classified as the class in $C$ right below
$X$. In order for our algorithm to work, we set the number of transfer size
classes equal to the number of bandwidth classes.

Figure 6: Example communication graph with different bandwidth classes

Figure 6 shows an example communication graph, where the nodes are colored
black and different bandwidth classes of edges are shown.

Given the array of transfer sizes $S$ and array of communication graph edges
$E_{c}$, the lower bound on bottleneck latency we can achieve is given by
Theorem 1.

###### Theorem 1

The lowest bottleneck latency we can achieve is:

$\min(\beta)=\frac{\max{S}}{\max{E_{c}}}$ (9)

Therefore, if we achieve $\min(\beta)$, then we have found the optimal minimum
bottleneck latency.

We prove Theorem 1 as follows:

Given the highest transfer size ($\max S$), then it must be matched with the
highest bandwidth ($\max E_{c}$) to have the lowest bottleneck latency. There
are two cases in which the system would have another bottleneck latency:

  1. 1.

$\beta=\frac{\max{S}}{e}\forall e\in E_{c}-max(E_{c})$ (10)

  2. 2.

$\begin{array}[]{l}\beta=\frac{s}{e}\\\ \forall\\{s\in S-max(S),\\\ e\in
E_{c}-max(E_{c})\mid\beta\geq\frac{\max{S}}{\max{E_{c}}}\\}\end{array}$ (11)

In Equation 10, the latency of the system would be higher than 9, since the
transfer size is being matched with a lower bandwidth edge. In Equation 11,
some other transfer size $s$ and bandwidth $e$ may result in a higher
bottleneck latency, in which case Equation 9 still holds. Therefore, Theorem 1
holds.

We run tests in Section V to see how often we get this optimal solution.
Algorithm 2 performs the matching between $S$ and $G_{c}$ to try to reach the
optimal latency as outlined above. Let $N$ represent the array of nodes that
we choose from $G_{c}$, with length $\lvert S\rvert$.

Algorithm 2 Finding K-Paths

procedure SUBGRAPH-K-PATH($X$, $k$, $s$, $u$)

// Sort by weight in descending order

$edgeList\leftarrow\textsc{SORT}(G_{c},\\{e\in E_{c}\mid w(e)\\},reverse)$

$low\leftarrow 0$

$high\leftarrow edgeList.length$

$bestPath\leftarrow[]$

while $low<high$ do

$median\leftarrow\frac{(low+high)}{2}$

$\tau(X,edgeList[median])$

$\tau(X,threshold)$

// Induced subgraph of $G_{c}$ w/ class $X$ edges

$G^{X}_{c}\leftarrow\\{E^{X}_{c}=\\{e\in E_{c}\land c(e)=X\mid e\\}\mid
V^{X}_{c},E^{X}_{c}\\}$

$result\leftarrow\textsc{K-PATH}(G^{X}_{c},k,s,u)$

if $result=\text{FALSE}$ then

$low\leftarrow median+1$

else

$high\leftarrow median$

$bestPath\leftarrow result$

end if

end while

for $N\in bestPath$ do

$\textsc{DEL}(G_{c},N)$

end for

end procedure

In algorithm 2, we use the color-coding k-path algorithm [2], which finds a
path of length $k$ (where $k$ is the number of vertices) in $G^{X}_{c}$ if a
$k$-path exists and does so in polynomial time if $k<\log(|V^{X}|)$. As we
show in Section V, this is possible using node memory capacities which reflect
real-world devices. We use a binary search to find the maximum threshold for
which a $k$-path exists. On line 3, we sort in descending order so that we can
find the maximum viable edge-weight threshold with a binary search. As $N$
starts to be filled in, the $k$-paths have to be found between certain nodes
in order for $N$ to be a contiguous path of nodes. We modify the $k$-path
algorithm to start at $s$ and stop once it reaches $u$. We make the algorithm
more efficient by stopping a particular iteration if we reach $u$ before we
have a path of length $k$. If $s$ is $null$, find any find any $k$-path that
ends at $u$. Similarly, if $u$ is $null$, find any $k$-path that starts at
$s$.

Algorithm 3 K-Path Matching

1:procedure K-PATH-MATCHING(S, C)

2: for $X\in C$ do

3: $x\\_paths\leftarrow\textsc{FIND-SUBARRAYS}(S,X)$

4: $x\\_paths\leftarrow\textsc{SORT}(x\\_paths,\\{p\in x\\_paths\mid
p.length\\})$

5: for $i\leftarrow 0\textbf{ to }x\\_paths.length$ do

6: $startIdx\leftarrow\textsc{INDEX-OF}(x\\_paths[i][0],S)$

7: $x\\_len\leftarrow x\\_paths.length$

8: $startV\leftarrow N[startIdx]$

9: $endV\leftarrow N[startIdx+x\\_paths[i].length+1]$

10: $path\leftarrow\textsc{SUBGRAPH-K-PATH}(X,x\\_len+1,startV,endV)$

11: $N[startIdx:startIdx+path.length]\leftarrow path$

12: end for

13: end for

14:end procedure

Algorithm 3 performs the $k$-path matching of partitions onto vertices of
$G_{c}$.

FIND-SUBARRAYS() on line3 returns a list of subarrays of a certain class, by
iterating over the list of transfer sizes $S$.

Lines 2-13 match paths for all bandwidth classes, starting with class $H$. On
line 8, $startV$ represents the vertex before the current $k$ path. If this is
equal to $null$, meaning that the algorithm hasn’t reached the iteration of
finding the previous $k$-path, then SUBGRAPH-K-PATH will find a path that
starts at any vertex. Similarly, on line 9 $endV$ represents the vertex of the
current $k$-path. If this is equal to $null$, meaning that the algorithm
hasn’t reached the iteration of finding that subsequent $k$-path, then
SUBGRAPH-K-PATH will find a path that ends at any vertex.

By starting with the longest $H$-subarrays and working to the shortest $L$
subarrays, we are greedily finding the best bandwidth paths to match with the
highest transfer size terms of $S$. We continue this process until we have
found $k$-path matchings for all subarrays of $S$. In some cases, a high
number of bandwidth classes will prevent the algorithm from returning a
result, because it has very few edges to choose from during each iteration of
the matching. In this case, we can re-run the algorithm with fewer bandwidth
classes.

##  IV Evaluation Methodology

We simulated a set of randomly placed edge devices using a random complete
graph. For each evaluation, we created a random complete graph by drawing the
positions of the nodes from a uniform distribution with the range
$x,y\in(-150,-1)\cup(1,150)$ to simulate a WiFi router with a range of 150m.
Between each set of nodes, we calculated the edge weight using Shannon’s
capacity equation, assuming that the SNR decreases proportionally to the
inverse square of the device’s distance from the router.

$\begin{array}[]{l}r(x,y)=\log_{2}\left(1+\frac{a}{(\sqrt{x^{2}+y^{2}})^{2}}\right)=\log_{2}\left(1+\frac{a}{x^{2}+y^{2}}\right)\\\
x,y\in(-B,-1)\cup(1,B)\end{array}$ (12)

In Equation 12, we found $a=283230$ by assuming that the bandwidth at 80 m
from the router was 5.5 Mbps, which matches the characteristics of a low-power
edge network. We exclude $x,y\in(-1,1)$ to satisfy the domain of Equation 12,
and to simplify the creation of our geometric graphs for our simulations. From
a practical standpoint, this means that we assume that no devices will be
within 1m of the router.

We used the following configuration to test the algorithm:

  1. 1.

Model: MobileNetV2 [27], EfficientNetB1 [31], ResNet50, or InceptionResNetV2,
chosen for their diversity in model topologies and memory sizes. Due to space
constraints, we only show ResNet50 and InceptionResNetV2 in some results
figures.

  2. 2.

Number of Nodes: 5, 10, 15, 20, or 50 randomly placed edge devices.

  3. 3.

Number of Bandwidth Classes: 2, 5, 8, 11, 14, 17, or 20 bandwidth classes,
which provide granularity in how to classify the transfer sizes and edge
bandwidths.

  4. 4.

Node Memory Capacity: 64, 128, 256, or 512 MB of RAM for a compute node.

For each test, we used a different random communication graph generated using
the procedure above. With each algorithm result, we then calculated the
bottleneck latency according to Equation 3. The resulting bottleneck latency
from each configuration of model, node capacity, number of nodes, and number
of bandwidth classes was run 50 times and averaged.

We compare the resulting bottleneck latency of our algorithm to that of the
following two algorithms:

  1. 1.

Random Algorithm: Select a random node and a random partition that can be
accommodated on that node.

  2. 2.

Joint-Optimization Algorithm: Let $Q$ and $N$ represent the optimal set of
partitions and optimal arrangement of nodes, respectively, chosen under this
algorithm For each node $n$, do the following:

    1. (a)

At each step choose the partition with the smallest transfer size that will
fit within the node. Add this partition to the set of chosen partitions $p$.

    2. (b)

Starting at $n$, find the neighbor in the communication graph whose edge $e$
has the highest bandwidth, and add that to the path of chosen nodes $c$. Then,
find the highest bandwidth edge from $e$, and so on.

    3. (c)

Compare the bottleneck latency found with $p$ and $c$ to the smallest
bottleneck found with all nodes $n$ thus far, and update $Q$ and $N$ with $p$
and $c$ if the current bottleneck is smaller.

For each of these algorithms, we used the same configuration and methodology
as above to find the bottleneck latency. These algorithms don’t use bandwidth
classes, so we didn’t need to include that as part of the configuration.

##  V Results

Figure 7: Color Map of Bottleneck Latency ($s$) based on Model, Node Capacity,
Number of Nodes, and Number of Bandwidth Classes - Optimal
Partitioning/Placement

Figure 7 shows a color map of the resulting bottleneck latency based on the
factors described in Section IV. In Figure 7, the color map was only generated
for the node capacities which were too small for the models to fit on a single
device of that capacity. All models were able to fit on a single 512 MB
device. The lack of bottleneck latency values for InceptionResNetV2 with 5
nodes and 64MB node capacity indicates that the model could not be partitioned
with these physical constraints. For each model, the lowest bottleneck latency
for a given node capacity comes from the combination of the most number of
bandwidth classes and number of nodes. The lowest bottleneck latency comes
with the highest node capacity. These results follow from the fact that a
larger node and number of nodes allows the partitioning algorithm to have
greater choice in selecting the smallest transfer sizes. Similarly, a high
number of bandwidth classes allows the placement algorithm to better perform
the $k$-path matching.

Figure 8: Comparison of Algorithm 3 with Random Algorithm - based on Model,
Node Capacity, Number of Nodes

In Figure 8, the optimal algorithm reduces bottleneck latency by $\approx 10$x
on average for this selection of models. The difference is the smallest for
ResNet50, with the optimal algorithm producing a $\approx 2$x lower bottleneck
latency. The models with the greatest variance in transfer size will result in
the largest difference in bottleneck latency between the optimal random
algorithms. Overall, we see that the optimal algorithm produces a significant
reduction in bottleneck latency compared to the random algorithm.

Figure 9: Comparison of Algorithm 3 with Joint Optimization - based on Model,
Node Capacity, Number of Nodes

In Figure 9, the joint optimization algorithm tends to perform better for a
smaller number of nodes. Since each of these algorithms use the same optimal
partitioning logic, we can only compare the models based on their differing
placement logic. As the number of nodes increases, our $k$-path algorithm
performs better. This makes sense, because the difference in the greedy
strategy of the joint optimization algorithm and the matching strategy of our
algorithm only becomes more apparent as the communication graph grows bigger
and there are more options for node paths. In particular, for 50 nodes, our
algorithm outperforms the joint optimization algorithm by 35%. We hypothesize
that this trend would continue for more complex models which have a greater
number of candidate partition points and a greater variance in transfer size,
necessitating the $k$-path matching strategy to minimize bottleneck latency.

Figure 10: Histogram of Average Approximation Ratio for Keras Pretrained
Models

We then ran our algorithm 1000 times for the set of Keras pretrained models.
We used a configuration of 50 nodes and 64 MB node memory capacity. For each
trial, we then divided the resulting bottleneck latency by the optimal
bottleneck latency as given by Theorem 1. We took the average of this ratio
across all trials for each model. The results are presented in Figure 10. We
see that the bottleneck latencies for 75% of these models are within 9% of the
optimal solution. The average approximation ratio across all these models is
$\approx 1.092$ or within 9.2% of the optimal bottleneck latency.

##  VI Conclusion

We have presented a framework to partition and place a model across a set of
resource-constrained edge devices, with the goal of maximizing inference
throughput. Additionally, we show that for different models and node
configurations, we can outperform a greedy joint-optimization algorithm. We
further show empirically that our algorithm is within 9.2% of the optimal
bottleneck latency for the models we tested.

###  VI-A Future Work

With minor edits, we could extend our framework to work with geographically-
distributed edge devices for a truly scalable edge inference solution.

With software changes, we could potentially run the average image model on a
cluster of micro-controllers. We could use RiotOS [26] without any
containerization and perform optimizations to run with limited device memory.
Some devices we could potentially take advantage of are the Raspberry Pi Pico
[10] and Arduino Uno [3].

##  VII Acknowledgements

We would like to acknowledge the helpful input and pointers provided by Prof.
Anil Vullikanti from the University of Virginia, particularly in directing us
to the color-coding $k$-path algorithm.

## References

  * [1] Abadi, M., et al. TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  Software available from tensorflow.org.
  * [2] Alon, N., Yuster, R., and Zwick, U. Color-coding.  Journal of the ACM (JACM) 42, 4 (1995), 844–856.
  * [3] Arduino.  Arduino uno.  https://store-usa.arduino.cc/products/arduino-uno-rev3, 2003.
  * [4] Biswas, A., Jain, A., et al. Survey on edge computing–key technology in retail industry.  In Computer Networks and Inventive Communication Technologies. Springer, 2021, pp. 97–106.
  * [5] Chai, J., Zeng, H., Li, A., and Ngai, E. W. Deep learning in computer vision: A critical review of emerging techniques and application scenarios.  Machine Learning with Applications 6 (2021), 100134.
  * [6] Chen, J., and Ran, X. Deep learning with edge computing: A review.  Proceedings of the IEEE 107, 8 (2019), 1655–1674.
  * [7] Cho, E., Yoon, J., Baek, D., Lee, D., and Bae, D.-H. Dnn model deployment on distributed edges.  In International Conference on Web Engineering (2021), Springer, pp. 15–26.
  * [8] Chollet, F., et al. Keras.  https://keras.io, 2015.
  * [9] Collet, Y. Lz4.  https://github.com/lz4/lz4, 2011.
  * [10] Foundation, R. P. Raspberry pi pico.  https://www.raspberrypi.com/products/raspberry-pi-pico/, 2021.
  * [11] Gholami, A., Kim, S., Dong, Z., Yao, Z., Mahoney, M. W., and Keutzer, K. A survey of quantization methods for efficient neural network inference.  arXiv preprint arXiv:2103.13630 (2021).
  * [12] He, K., et al. Deep residual learning for image recognition.
  * [13] Hsu, K.-J., Bhardwaj, K., and Gavrilovska, A. Couper: Dnn model slicing for visual analytics containers at the edge.  In Proceedings of the 4th ACM/IEEE Symposium on Edge Computing (2019), pp. 179–194.
  * [14] Hu, D., and Krishnamachari, B. Fast and accurate streaming cnn inference via communication compression on the edge.  In 2020 IEEE/ACM Fifth International Conference on Internet-of-Things Design and Implementation (IoTDI) (2020), IEEE, pp. 157–163.
  * [15] Keras.  Keras applications.  https://keras.io/api/applications/, 2021.
  * [16] Li, E., Zhou, Z., and Chen, X. Edge intelligence: On-demand deep learning model co-inference with device-edge synergy.  In Proceedings of the 2018 Workshop on Mobile Edge Communications (2018), pp. 31–36.
  * [17] Li, M., Gao, J., Zhou, C., Zhuang, W., et al. Slicing-based ai service provisioning on network edge.  arXiv preprint arXiv:2105.07052 (2021).
  * [18] Lindstrom, P. Fixed-rate compressed floating-point arrays.  IEEE Transactions on Visualization and Computer Graphics 20 (08 2014).
  * [19] Luo, Q., Hu, S., Li, C., Li, G., and Shi, W. Resource scheduling in edge computing: A survey.  IEEE Communications Surveys & Tutorials 23, 4 (2021), 2131–2165.
  * [20] Matsubara, Y., Baidya, S., Callegaro, D., Levorato, M., and Singh, S. Distilled split deep neural networks for edge-assisted real-time systems.  In Proceedings of the 2019 Workshop on Hot Topics in Video Analytics and Intelligent Edges (2019), pp. 21–26.
  * [21] Min, B., Ross, H., Sulem, E., Veyseh, A. P. B., Nguyen, T. H., Sainz, O., Agirre, E., Heinz, I., and Roth, D. Recent advances in natural language processing via large pre-trained language models: A survey.  arXiv preprint arXiv:2111.01243 (2021).
  * [22] Nayak, S., Patgiri, R., Waikhom, L., and Ahmed, A. A review on edge analytics: Issues, challenges, opportunities, promises, future directions, and applications.  arXiv preprint arXiv:2107.06835 (2021).
  * [23] Parthasarathy, A., and Krishnamachari, B. Defer: Distributed edge inference for deep neural networks.  In 2022 14th International Conference on COMmunication Systems & NETworkS (COMSNETS) (2022), IEEE, pp. 749–753.
  * [24] Premkumar, S., and Sigappi, A. A survey of architecture, framework and algorithms for resource management in edge computing.  EAI Endorsed Transactions on Energy Web 8, 33 (2021), e15–e15.
  * [25] PyTorch.  Real time inference on raspberry pi 4.  https://pytorch.org/tutorials/intermediate/realtime_rpi.html, 2021\.
  * [26] Riot.  Riotos.  https://www.riot-os.org/, 2013.
  * [27] Sandler, M., Howard, A., Zhu, M., Zhmoginov, A., and Chen, L.-C. Mobilenetv2: Inverted residuals and linear bottlenecks.  In Proceedings of the IEEE conference on computer vision and pattern recognition (2018), pp. 4510–4520.
  * [28] Sonkoly, B., Czentye, J., Szalay, M., Németh, B., and Toka, L. Survey on placement methods in the edge and beyond.  IEEE Communications Surveys & Tutorials 23, 4 (2021), 2590–2629.
  * [29] Stahl, R., Hoffman, A., Mueller-Gritschneder, D., Gerstlauer, A., and Schlichtmann, U. Deeperthings: Fully distributed cnn inference on resource-constrained edge devices.  International Journal of Parallel Programming 49, 4 (2021), 600–624.
  * [30] Szegedy, C., Ioffe, S., Vanhoucke, V., and Alemi, A. A. Inception-v4, inception-resnet and the impact of residual connections on learning.  In Thirty-first AAAI conference on artificial intelligence (2017).
  * [31] Tan, M., and Le, Q. Efficientnet: Rethinking model scaling for convolutional neural networks.  In International conference on machine learning (2019), PMLR, pp. 6105–6114.
  * [32] Tensorflow.  Tensorflow lite.  https://www.tensorflow.org/lite, 2019.
  * [33] Wu, R., Guo, X., Du, J., and Li, J. Accelerating neural network inference on fpga-based platforms—a survey.  Electronics 10, 9 (2021), 1025.
  * [34] Zhang, Q., Sun, H., Wu, X., and Zhong, H. Edge video analytics for public safety: A review.  Proceedings of the IEEE 107, 8 (2019), 1675–1696.
  * [35] Zhang, Z., Li, Y., Guo, Y., Chen, X., and Liu, Y. Dynamic slicing for deep neural networks.  In Proceedings of the 28th ACM Joint Meeting on European Software Engineering Conference and Symposium on the Foundations of Software Engineering (2020), pp. 838–850.
  * [36] Zhang, Z., Ng, L. K., Liu, B., Cai, Y., Li, D., Guo, Y., and Chen, X. Teeslice: slicing dnn models for secure and efficient deployment.  In Proceedings of the 2nd ACM International Workshop on AI and Software Testing/Analysis (2022), pp. 1–8.
  * [37] Zoph, B., Vasudevan, V., Shlens, J., and Le, Q. V. Learning transferable architectures for scalable image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition (2018), pp. 8697–8710.


Spin-flip pair-density functional theory: A practical approach to treat
static and dynamical correlations in large molecules

Oinam Romesh Meitei
Department of Chemistry, Virginia Tech, Blacksburg, Virginia 24061, USA

Nicholas J. Mayhall
Department of Chemistry, Virginia Tech, Blacksburg, Virginia 24061, USA

We present a practical approach to treat static and dynamical correlation
accurately in large multi-configurational systems. The static correlation is
accounted for using the spin-flip approach which is well known for capturing
static correlation accurately at low-computational expense.
Unlike previous approaches to add dynamical correlation to spin-flip models
which use perturbation theory or coupled-cluster theory, we explore the
ability to use the on-top pair-density functional theory approaches recently
developed by Gagliardi and co-workers (JCTC, 10, 3669, 2014).
External relaxations are carried out in the spin-flip
calculations though a restricted active space framework for which a
truncation scheme for the orbitals used in the external excitation is
The performance of the approach is demonstrated by computing energy gaps
between ground and excited states for diradicals, triradicals and linear
polyacene chains ranging from naphthalene to dodecacene. Accurate results
are obtained using the new approach for these challenging open-shell
molecular systems.

§ INTRODUCTION

An accurate computational description of large molecular systems with multi-configurational
characters or strongly correlated systems still remains a challenge
due to the lack of a general approach that delivers high accuracy while remaining
computationally affordable.

The most widely used approach for treating multi-configurational systems is the
complete active space self-consistent field (CASSCF) method in which the full
configurational interaction (CI) is solved within a chosen active
space with the orbitals and the CI coefficients optimized
<cit.>. The active space is typically selected with a
prior knowledge of the important
orbitals contributing to the target chemical problem or the electronic
state (although recent advances are making it possible to algorithmically select active space orbitals <cit.>). When supplied a well-defined active space, the method has
been shown to reliably recover static correlation. However in
doing so, the bulk of the dynamical correlation is absent, precluding any hope for quantitative predictions.
A general strategy adopted to recover the dynamical correlation is to use the
complete active space second-order perturbation theory (CASPT2) approach.<cit.>
As dynamical correlation is well-described by perturbation theory,
CASPT2 provides significant accuracy improvements compared to CASSCF, albeit with increased computational cost and memory requirements. The CASSCF method
itself is limited to 20 electrons in 20 orbital space with the current
state-of-the-art computational resources.<cit.>
With a motivation to enable computation using larger active spaces, several
approaches have been developed such as
various flavors of selected configuration interaction (SCI)
approaches<cit.>, the density matrix renormalization group (DMRG)<cit.>, generalized
active-space self-consistent field (GASSCF) method<cit.>, the 2-RDM driven CASSCF
(v2RDM-CASSCF) method<cit.> or the full CI quantum Monte Carlo
(FCIQMC) method<cit.>.
With these reduced cost methods, larger active spaces can be used which begin to account for dynamical correlation as well, but generally a large amount of dynamical correlation is still missing.

An alternative approach for treating multi-configurational systems is the
spin-flip (SF) method<cit.>.
The SF method is a relatively simple approach which
uses a well-defined single-determinant reference on which spin-flipping
excitations are carried out to access electronic states starting from a
high-spin configuration. The SF methods have been further extended to remain
spin-pure, including external relaxations though a restricted active space
(RAS-$n$SF)<cit.>. Recently, we have also proposed a redox SF approach that
simultaneously accounts for both spin and spatial degeneracies which the SF
approach alone cannot handle<cit.>.
The SF approaches have favorable computational scaling for practical
applications to large multi-configurational systems.
In addition to allowing traditionally single-reference methods to be applied to multi-reference problems,
(e.g., coupled-cluster theory), spin-flip also reduces the dependence on user input, since the active-space is defined automatically once the reference spin state has
been chosen.
As a result, prior knowledge of the orbitals that are important for the target
electronic state is not necessarily required.
Here, the active space is simply chosen as the singly occupied
orbitals. The external excitation space in RAS-$n$SF are simply the doubly
occupied and the virtual orbitals although knowledge of the important
contributing orbitals can be utilized as well.
Another advantage of the SF approach is the treatment of the
electronic states including the ground state on equal footing.
However, because the RAS-$n$SF approach excludes all double excitations outside of the active space,
it can only deliver qualitative accuracy due to the lack of dynamical correlation.<cit.>

An appealing approach to recover dynamical correlations for practical
applications to large molecular systems is the combination of multi-configurational
wavefunction methods with a density functional theory (DFT) based method.
Here, the idea is to use a qualitatively correct multi-configurational wavefunction to capture the
static correlation and include a DFT-based description of the dynamical correlation.
The idea was first introduced by Lie and Clementi where they demonstrated that
the DFT results can be substantially improved by adding the correlation energy
obtained from a multi-configurational wavefunction<cit.>.
However, the approach does not separate the static and dynamical correlation,
thereby double counting the correlation energy.
Several approaches have emerged in this line which attempts to overcome the
deficiencies and improve over the previous works<cit.>.
Notable example include the
approach by Colle and Solvetti<cit.>, CAS-DFT by Savin and co-workers<cit.>, modified
CAS-DFT approaches by Grafenstein and
Cremer<cit.>, Yamaguchi and
co-workers<cit.> and CAS-DFT using on-top
pair density by Gusarov and co-workers<cit.>. A nice
review of the available approaches are provided in Ref .
The main challenge in combining multi-configurational wavefunction with a
DFT-based description of the dynamical correlation is due to the double counting of
electron correlation. The multi-configurational wavefunction generally includes
some part of the dynamical correlation within the active space. The second
complication is the “symmetry dilemma” in KS-DFT<cit.>. The spin densities of the
multi-configurational wavefunctions are not compatible with standard density
functionals for low spin states.

Recently, a new
approach has been proposed by Gagliardi and co-workers to include DFT-based
correlation energy to a multi-configurational wavefunction
which addresses both the double counting of correlation energy as well as the
“symmetry dilemma”.<cit.> The method referred to as multi-configurational-pair
density functional theory (MC-PDFT) only uses the multi-configurational
wavefunction to compute the classical Coulomb and the kinetic energy while
the rest of the exchange-correlation energy and correction to the kinetic
energy is obtained from DFT using the on-top pair density functionals. The
success of the method has been demonstrated in conjunction with CASSCF<cit.>,
GASSCF<cit.>, DMRG<cit.> and v2RDM-CASSCF<cit.>.

In this work, we propose an alternate approach to combine MC-PDFT with the
RAS-$n$SF approaches and our recently developed SF-IP/EA approach. The goal is
to recover dynamical correlation energy in the SF approaches using MC-PDFT for
practical applications to large multi-configurational systems at a considerably
low computational cost. A brief description of the RAS-$n$SF approaches and
the MC-PDFT method is provided in Section <ref> and
<ref>. A practical scheme is also presented to further reduce the
computational cost in the MC-PDFT calculation. Computational details are
provided in Section <ref>. The performance of the approach
is demonstrated by computing energy gaps between spin states of challenging
biradicals, triradicals and polyacenes ranging from naphthalene to
dodecacene. The results are presented in Section <ref>,
<ref> and <ref>. Finally, a summary is provided in Section

§ METHODS

§.§ RAS-$n$-SF and combined IP/EA approach

We give a short overview of the RAS-$n$-SF and
the combined IP/EA approaches. For further details of the methods, we refer the reader to
Refs. , , and .

The spin-flip (SF) approach, proposed by Krylov, provides an efficient way to model
a large number of multi-configurational problems using only a single reference
The key idea in SF approach is that while the various S$_z$ multiplets of high
spin states have identical electronic energies, states with maximum S$_z$ have a
single determinant representation, while lower S$_z$ multiplets are highly
multi-configurational. SF approach leverages this degeneracy, by optimizing the orbitals
for the single configurational high spin (S$_z$ = S) state, then uses
spin-flipping excitations to access the target S$_z$ manifold of states. Because
both high-spin and low-spin states appear in the target S$_z$ space,
the high-spin states and low-spin states are treated on an equal footing.

As an
example, the spin-degeneracy in the valence bonding orbitals
$\ket{\sigma^{2}\pi^{4}}$ with the anti-bonding $\sigma^{*}$ and $\pi^{*}$
orbitals upon bond dissociation of N$_2$ molecule can be resolved  using the
SF approach. Using a high-spin heptet state
($\ket{\sigma^1\pi^2\pi^{*2}\sigma^1}$) with $m_s=3$ as the reference state,
the ground singlet state can be accessed using a 3SF operator.

\begin{equation}
\ket{\Psi}=\sum_{i<j<k,\bar{a}<\bar{b}<\bar{c}}c_{ijk}^{\bar{a}\bar{b}\bar{c}}\hat{a}_{ijk}^{\bar{a}\bar{b}\bar{c}}\ket{\Psi^{Reference}}
\end{equation}

Systems with spin as well as spatial degeneracies can be treated using a
combination of the spin-flip method and the electron addition (EA) or
elimination (IP) method. Unlike in the SF approach, the combined SF-IP/EA
approach works with the closest well-defined high spin state obtained either
by oxidizing or reducing the the system. While only spatial degeneracies arise in the previous example,
if the cationic system is considered instead (N$_{2}^{+}$) both spin
and spatial degeneracy occur. In this case, the system is multi-configurational even in the high
spin state as one of the degenerate bonding ($\sigma^{2}\pi^{3}$) and anti-bonding
($\sigma^{*}\pi^{*}$) orbitals must be left unoccupied upon dissociation. The
ambiguity is resolved by simultaneously using the 2SF-IP operator on
the neutral N$_2$ molecule.

\begin{equation}
\ket{\Psi}=\sum_{ijk\bar{a}\bar{b}}c_{ijk}^{\bar{a}\bar{b}}\hat{a}_{ijk}^{\bar{a}\bar{b}}\ket{\Psi^{Reference}}
\end{equation}

The SF and IP/EA excitations are carried out only in the singly occupied
orbital space to ensure spin-pure state. External relaxation effects are taken
into account using RAS(S) by allowing the full set of singles excitations
defined by (h, p, hp)<cit.>. The approach without the external effects is
denoted by CAS-$n$SF. Overall, the
RAS-$n$-SF and the combined IP/EA approaches provides qualitatively accurate
descriptions of static correlation as demonstrated in earlier works with the
only disadvantage being the neglect of dynamical correlations.

Truncation schemes used for the RAS1 and RAS3 active space in the
RAS(S) excitation in the natural orbital basis. The truncations are
performed separately for the RAS1 and RAS3 space.

Scheme     Abbreviation     Description [2mm]
I          RAS-SS           A state specific approach, truncate the active
space for each state separately. Each state has different number of
orbitals in the active space[2mm]

II         RAS-SA           A state average (SA) approach, truncate based on
an average density of the participating states. Each state has the same
number of orbitals in the active space.[2mm]

III        RAS-eff          An effective SA approach, the active space is
defined by the maximum of the number of orbitals in each state for a
given threshold (separately for RAS1 and RAS3).[2mm]

§.§ Multi-configurational pair-density functional theory

We shortly review the MC-PDFT approach which combines multi-configurational
methods with DFT-based methods without incurring double counting of
correlation energy and the symmetry dilemma in the context of KS-DFT<cit.>. By using
1- and 2-RDMs from a multi-configurational calculation, the electronic energy
in the MC-PDFT framework is given by

\begin{equation}
E = \sum_{pq}h_{pq}D_{pq} + \frac{1}{2}\sum_{pqrs}g_{pqrs}D_{pq}D_{rs} + E_{ot}[\rho,\Pi]\label{eq:mc_pdft_eq}
\end{equation}

where $h_{pq}$ and $g_{pqrs}$ are the one- and two-electron integrals
respectively, $D_{pq}$ is the 1-RDM and $E_{ot}$ is the on-top energy with
$\rho$ and $\Pi$ being that total density and the on-top pair density
respectively. The last two terms are the classical Coulomb term and an on-top
density functional term which replaces the two-electron contribution in the
usual electronic energy expression for many-electron systems. Here, the
one-electron contribution containing the kinetic and electron-nuclear
potential energy as well as the classical electrostatic contribution are
directly taken from the multi-configurational wavefunction. The remainder
exchange and the correlation contributions are folded into the on-top pair
density functional.

The total density $\rho$ and the on-top pair density, $\Pi$ are defined in
terms of 1-RDM and 2-RDM obtained from the multi-configurational calculations
respectively as

\begin{equation}
\rho(r) = \sum_{pq}\phi_p(r)\phi_q(r)D_{pq} \label{eq:rdm1}
\end{equation}

\begin{equation}
\Pi(r) = \sum_{pqrs}\phi_p(r)\phi_q(r)\phi_r(r)\phi_s(r)D_{pqrs} \label{eq:rdm2}
\end{equation}

The on-top pair density functional used within the MC-PDFT formalism is simply
obtained with “translated” existing exchange-correlation functionals
employed in standard KS-DFT. The derivation of the “translated” and the
corresponding “fully translated” functionals are provided in Ref. .

In the present RAS-$n$SF(-IP/EA)-PDFT strategies, the appropriate RDMs entering in
Equation <ref> are obtained from the RAS-$n$SF(-IP/EA) wavefunction
presented in the previous section. Because the computation of the 2-RDM in the combined spaces of
RAS1, RAS2, and RAS3 would prevent application to large systems,
we propose the following truncation scheme which improves efficiency while  producing only small affects on the final energies:

* Perform RAS-$n$SF(-IP/EA) using the full space for the external RAS(S)

* Diagonalize 1-RDMs of the RAS1 and RAS3 spaces separately to obtain natural orbitals for
the RAS subspaces. The RAS-$n$SF methods are invariant with respect to orbital rotations within a
RAS subspace.

* Define a threshold for the orbital occupation number to separately
truncate the RAS1 and RAS3 active space for the RAS(S) excitation in the
natural orbital basis. The truncation schemes are tabulated in Table <ref>

* Using the truncated RAS1/RAS3 spaces repeat the RAS-$n$SF(-IP/EA) calculation
to obtain the required one- and two-RDMs for the RAS-$n$SF(-IP/EA)-PDFT calculation.

Schematic representations of the poly-radicals and polyacenes considered in this work.

§ COMPUTATIONAL DETAILS

In order to determine the accuracy of the RAS-$n$SF(-IP/EA)-PDFT approach
presented in this work, we compare to existing accurate methods for the
computation of the doublet-quartet gap of triradical and singlet-triplet gaps
in aromatic biradicals and linear polyacenes and to experimental results
wherever available. The structures of 1,3,5-trimethylenebenzene (TMB) and its
negative ion (TMB$^-$) were optimized with UB3LYP/cc-pVDZ. The structure of
meta-xylylene was taken from Ref. , the $ortho$-, $para$- and $meta$-benzyne
radicals from Ref. ,  and the
polyacenes from naphthalene to dodecacene from
Ref. . A schematic representation of the poly-radicals
and polyacenes considered in this work is provides in Figure <ref>.

The doublet state of neutral 1,3,5-trimethylenebenzene triradical was obtained
by performing a single SF operation on the quartet reference state with
$m_{s}=\tfrac{3}{2}$. The singlet and triplet state of its negative ion was
obtained by performing 1SF-EA operation on the neutral high-spin quartet
state obtained by oxidizing the target anionic state. For the benzyne radicals
as well as the polyacenes, the singlet states are obtained by performing only
a single SF operation on the triplet reference state as done in Refs.  and .

Augmented Dunning's correlation consistent aug-cc-pVDZ basis set was employed
for all the calculations<cit.>.
The RAS-$n$SF-PDFT and the RAS-$n$SF-IP/EA-PDFT calculations were performed
with tPBE, ftPBE, tBLYP and ftBLYP on-top density functional.
Cholesky decomposition was used in all the two-electron integral calculations
with a decomposition threshold of 10$^{-6}$ a.u.<cit.>
The RAS-$n$SF(-IP/EA) and the RAS-$n$SF(IP/EA)-PDFT
calculations were performed using OpenMolcas<cit.> and the truncation scheme described in Section <ref> was computed using an in-house python plug-in code to OpenMolcas.

Illustration of the operation of 1SF and 1SF-EA on the high spin reference
state for (a) TMB, (b) TMB$^-$ and (c) MX. The orbitals shown corresponds to the
non-bonding orbitals and the spin-flips are indicated by
red in color.

§ RESULTS AND DISCUSSION

§.§ All $\pi$-polyradicals

The 1,3,5-trimethylenebenezene (TMB) presents as a prototypical high-spin
all-$\pi$ triradical<cit.>.
The triradical system has an open-shell quartet ground
state with the three unpaired electrons each occupying the nondisjoint
$\pi$ non-bonding orbitals (NBO).
Here, the three NBOs are nearly degenerate and most of the low-lying states have
heavily multiconfigurational wavefunctions.
theoretical studies are available for the energy gap between the ground state
and the energetically lowest doublet state<cit.>. The doublet state is accesible
with just a single spin-flip operation on the high spin quartet state
($m_s=\tfrac{3}{2}$), see the sketch in Figure <ref>(a).

The negative ion of TMB (TMB$^-$) on the other hand cannot be described solely by the
spin-flip approach. Four electrons occupy the nearly degenerate three NBOs and
it is ambigous to which of the three orbitals should be doubly occupied.
The ground state is a triplet state and thus is basically a
$\pi\pi$-diradical<cit.>. The ground state as well as the low-lying singlet states
can be accessed using 1SF-EA operation on the quartet reference state obtained
by oxidizing the anion. This is illustrated in Figure <ref>(b).

Another prototypical all $\pi$-diradical is the $meta$-xylylene (MX) with two
electrons distributed in two nearly degenerate NBOs. The singlet-triplet gap
of MX has been intensively studied both theoretically as well as
experimentally<cit.>. The diradical has an open-shell triplet ground state. The
low-lying singlet state can be accessed using a single SF operation on the
triplet reference state ($m_s=1$), see the illustration in Figure

Convergence of energy gap w.r.t. truncation of natural orbitals,
doublet-quartet gap of TMB (top), singlet-triplet gap of TMB$^-$ (center)
and MX (bottom). Natural orbitals with occupation number $<$ threshold
$or$ $> 2.0-$threshold are truncated in the RAS-$n$SF(-IP/EA)-PDFT
computation following the truncation schemes (in parenthesis), SS, SA, eff as outlined
in Table <ref>.

Figure <ref> plots the convergence of the different truncation schemes
described in the last part of Section <ref>.
The natural orbitals computed from the 1-RDM of the RAS3 active space with
occupation number smaller than a given threshold are truncated in the RAS(S)
excitation scheme. The doubly occupied orbitals in the natural orbital basis
constructed from the 1-RDM of the RAS1  space is truncated with
occupation number larger than the difference of two and the given threshold.
In the three truncation schemes namely, SS, SA, and eff, the distinction
between the SS and the eff schemes vanishes at the convergence.
The RAS(S)-$n$SF-PDFT converges with a threshold of 0.0001 for all the
truncation schemes with the exception of SS truncation scheme in
TMB$^-$. Here, the difference in the energy gap between spin states of TMB$^-$
obtained using the SS and eff truncation scheme is 0.13 kcal/mol at the
tightest threshold considered
which is well within the target accuracy of the presented approach.
Furthermore, the SA truncation scheme is slightly slower to converge for TMB$^-$ and MX, with the largest difference at the minimum threshold considered being 0.24 kcal/mol. Hereafter, we will consider the eff truncation
with a threshold of 0.0001 for the RAS3 orbitals and 1.9999 for the RAS1

The energy gap between spin-states of TMB, TMB$^-$ and MX computed using
CAS-SF, CASSCF, RAS(S)-SF and the corresponding MC-PDFT method using the tPBE
functional is provided in Table <ref>. Literature values for
the energy gaps from previous electronic structure calculations as well the
experimental energy gap for MX are also reported in therein.

For the TMB and TMB$^-$ radicals, experimental values
were not available and so comparisons are made to the available literature
values. First, we consider the vertical doublet-quartet (DQ) energy gap
of TMB radical. The application of MC-PDFT on CAS-1SF which does not account
for any external relaxation effect decreases the DQ energy gap by almost 15.0
kcal/mol. CAS-1SF-tPBE strongly underestimates the DQ energy gap as compared
to the previously reported gaps from CASPT2 and DDCI methods. Contrary to
this, CASSCF-tPBE increase the DQ energy gap computed with CASSCF (with
CAS(2,2)) by 7.3 kcal/mol. Interestingly, the CASSCF-tPBE energy gap is very close to a
previously reported CASSCF value with a larger active space (CAS(9,9)). On the
other hand RAS(S)-1SF-tPBE reduces the DQ energy gap from the bare RAS(S)-1SF
by almost 10 kcal/mol. Comparing to the available DDCI
and CASPT2 literature, the inclusion of dynamical correlation in the MC-PDFT
framework greatly improves the bare RAS(S)-1SF doublet-quartet energy

Next, we consider the vertical singlet-triplet (ST) gap for the anion
of TMB. Here, CAS-1SF as well as CAS-1SF-tPBE predicts a singlet ground state
as opposed to the experimental evidences for a high spin triplet ground state
for the anion<cit.>. The ST gap of bare CASSCF with
CAS(2,2) active space is increased by about 2.0 kcal/mol with the
corresponding MC-PDFT method. In contrast to this, RAS(S)-1SF-tPBE decreases
the ST gap computed using the bare RAS(S)-1SF by 1.7 kcal/mol. Here, the ST
gap agrees well with an existing CASPT2 result available in literature based
on CASSCF orbital with CAS(10,9).

For the MX diradical, the experimental ST energy gap is available and
so comparisons can be made to the adiabatic ST energy gap from the various
methods tabulated in Table <ref>. The CAS-1SF-tPBE strongly
underestimates the ST energy gap which collectively with the results obtained
for the TMB and TMB$^-$ radicals suggest that the one- and two-RDMs from the bare
CAS-1SF without any external relaxation effect is not accurate to be used in
the MC-PDFT equation. The CAS-1SF-tPBE even produces qualitatively wrong
ordering of spin-states in the TMB$^-$ radical.
On the other hand, the CASSCF-tPBE based on CAS(2,2) active space
underestimates the experimental ST gap of MX by 2.2 kcal/mol whereas the
RAS(S)-1SF-tPBE underestimates the energy gap by almost 4 kcal/mol. In
literature, the different wavefunction
based methods tend to overestimate the ST gap of MX, see Table
<ref>. LCCQMC based on a stochastic approach from Alavi
et. al. provides the ST gap of MX very close to experimental result<cit.>. The difficulty
in achieving an ST gap using wavefunction-based approaches comparable to the
experimental result have been demonstrated in earlier works<cit.>. An extensive
study using different wavefunction methods and DFT methods for the ST gap of
MX radical can be found in Ref. . The authors therein pointed out
the importance of including dynamical correlation with high accuracy<cit.>. In our
case, while the bare RAS(S)-1SF which also includes some dynamical correlation<cit.>
overestimates the experimental ST gap by 4.3 kcal/mol, RAS(S)-1SF-tPBE
underestimates the ST by 3.9 kcal/mol.

Energy gaps between spin states of TMB[Optimized geometry with UB3LYP/cc-pVDZ] (vertical doublet-quartet
gap), TMB$^-$<ref> (vertical singlet-triplet
gap) and MX[Geometry from Ref. , UB3LYP/6-311G(d,p)] (adiabatic singlet-triplet gap) obtained
using RAS(S)-1SF and RAS(S)-1SF-tPBE. The energy gap for TMB$^-$ was
obtained using 1SF-EA operation. All results with aug-cc-pVDZ basis set
and natural orbitals truncated with a threshold of 0.0001. Units are in kcal/mol.

TMB          TMB$^-$          MX [2mm]
CAS-1SF              17.6     -1.9     4.3[1mm]
CAS-1SF-tPBE         2.8      -6.1     1.7[2mm]
CASSCF[CAS(3,3) for TMB, CAS(4,3) for TMB$^-$ and CAS(2,2) for MX]      8.3     2.9      2.1[1mm]
CASSCF-tPBE<ref>          15.6     5.0      7.7[2mm]
RAS(S)-1SF           19.6     5.5     14.2[1mm]
RAS(S)-1SF-tPBE       9.8     3.8      6.0[2mm]
expt.                       9.6[1mm]
expt.$-$ZPE                 9.9[2mm]
15.7[CASSCF(9,9)/6-31G(d,p) from
Ref. ] ,
11.2[DDCI/6-31G(d) using localized orbitals and
complete virtual orbitals (fully variational) from
Ref. ],
13.6[CASPT2/6-31G(d,p) with CASSCF(9,9)
orbitals from Ref. ]
3.5[CASSCF(10,9)/ANO-L and CASPT2/ANO-L using
the same CASSCF orbitals from Ref. ],
1.9[state averaged MS-CASPT2/ANO-L from Ref. ]

12.9[CASSCF/6-31G(d) with CAS(8,8) from Ref. ],
11.7[CASPT2/6-31G(d) with CASSCF orbitals with CAS(8,8) from
Ref. ],
11.3[EOM-SF-CCSD/6-31G(d) from Ref. ],
11.8[Multireference second-order Møller-Plesset/aug-cc-pVTZ from
Ref. ],
9.5[LCCQMC/6-311++g(d,p) from Ref. ]

Adiabatic singlet-triplet energy gap for $ortho$-, $meta$- and
$para$-benzyne using the RAS-1SF and RAS-1SF-PDFT method. All results
with aug-cc-pVDZ basis set and natural orbitals
truncated with a threshold of 0.0001. Units are in

$o$-benzyne     $m$-benzyne     $p$-benzyne [2mm]
CAS-1SF            16.7     2.3     0.6 [1mm]
CAS-1SF-tPBE       30.2    14.1     3.4 [2mm]
CASSCF[using CAS(2,2) active space]     27.6    10.1     1.5 [1mm]
CASSCF-tPBE<ref>        39.0    23.4     5.4 [2mm]
RAS(S)-1SF            38.0     22.1     3.8 [1mm]
RAS(S)-1SF-tPBE       34.9     20.0     5.8 [1mm]
RAS(S)-1SF-ftPBE      33.0     17.7     5.4 [1mm]
RAS(S)-1SF-tBLYP      35.5     19.0     6.0 [1mm]
RAS(S)-1SF-ftBLYP     33.7     17.0     5.6 [2mm]
expt.[expt. value from
Ref. ]                 37.5     21.0     3.8 [1mm]
expt.$-$ZPE[ZPE correction with SF-DFT/6-311G* from
Ref. ]
38.1     20.0     3.3 [2mm]
35.1[CASSCF using CAS(8,8) active space from
Ref. ],
32.6[CASPT2 using CASSCF(8,8) orbitals from
Ref. ],
36.8[ic-MRCCSD(T) based on CAS(2,2) from Ref. ]
37.3[SF-CCSD(T) from
Ref. ]

16.1<ref>, 19.0<ref>,

3.8<ref>, 5.8<ref>,

§.§ Benzyne radicals

The $ortho$-, $meta$- and $para$-benzyne isomers have been used as
benchmarks for new theoretical approaches and at the same time challenging
because of the strongly correlated
biradical electrons<cit.>. The benzyne
isomers have a closed-shell singlet state
with the diradical character correlating with the distance between the
unpaired electrons. As a result, the energy gap between the ground singlet and
the lowest excited triplet state decreases following the $ortho$-, $meta$- and
$para$- sequence.

The adiabatic singlet-triplet gaps of the benzyne isomers computed using bare
CAS-1SF, CASSCF, RAS(S)-1SF and the corresponding MC-PDFT methods are presented
in Table <ref> along with the available experimental values and
selected literature values for comparison.
The RAS(S)-1SF and RAS(S)-1SF-PDFT ST energy gaps were computed using
truncated natural orbitals outlined in Section <ref>.
A threshold of 0.0001 and 1.9999 was used for the virtual and doubly
occupied space respectively following the results from Section
<ref> in the external RAS(S) excitation scheme.
The different truncation schemes, namely, SS, SA and eff resulted in similar
ST gaps and so only results from the eff truncation scheme is presented. The
largest difference in the ST gaps between the different truncation schemes was
only 0.22 kcal/mol. Also presented in Table <ref> are the results
from employing different on-top functionals, viz., tPBE, tBLYP and the fully
translated variants. Overall the results obtained from using tPBE performs
better than the other functionals and so we focus our discussion only to using
the tPBE functional.

The CAS-1SF and the CASSCF methods are based on minimal active spaces, i.e.,
CAS-1SF does not account for any external relaxations while the CASSCF method
is based on CAS(2,2) active space. This is reflected in the ST gaps obtained
from CAS-1SF-tPBE as compared to the experimental results (ZPE included) for
$ortho$- and $meta$-benzyne radicals. In spite of the small active space,
CASSCF-tPBE on the other hand performs well for the $ortho$-benzyne radical,
within 1.0 kcal/mol of the experimental value whereas the ST gap is
overestimated by 3.4 kcal/mol for the $meta$-benzyne radical.
Contrary to this, RAS(S)-1SF-tPBE underestimates the experimental ST gap by
3.2 kcal/mol for $ortho$-benzyne while it exactly matches the result for the
$meta$-benzyne. Here, the ST gap of $ortho$-benzyne obtained from RAS(S)-1SF
is within 0.1 kcal/mol of the experimental value. We note that the
corresponding bare RAS(S)-1SF method accounts for some dynamical correlations
but we emphasise again that the MC-PDFT method avoids double counting of
electron correlation. Comparing to literature, we find that for the $ortho$-
and $meta$-benzyne radicals, RAS(S)-1SF-tPBE is
more accurate than CASPT2 (based on CASSCF orbitals with CAS(8,8) active
space)<cit.> which is often a method of choice for including dynamical correlations
in modelling strongly correlated systems.

The $para$-benzyne radical presents as a more challenging case and has been
the focus in comparing the ST gap obtained from various highly accurate
theoretical methods to the available experimental result<cit.>. The discrepency
between experimental and theoretical results have been discussed in earlier
works on the diradical<cit.>. Köhn and co-workers have pointed out the possibility
of wrong assignments of the singlet and triplet states in the
experiment<cit.>. The
ST gap for the diradical presented in Table <ref> further supports
this possiblity. The ST gap from CAS-1SF-tPBE in which the bare CAS-1SF does
not account for any external relaxation agrees well within 0.1 kcal/mol of the
experimental result in contrast to RAS(S)-1SF-tPBE.
This would mean that CAS-1SF provides a more accurate one- and two-RDMS than
the RAS(S)-1SF in the MC-PDFT equation in contrast to the results obtained for
the other radicals discussed before. Furthermore, the ST gap form RAS(S)-1SF agrees well with other
literature values obtained from CASPT2 (based on CASSCF orbitals with CAS(8,8)
active space)<cit.> and ic-MRCCSD(T)<cit.>. The ST gap also agrees well with the

Adiabatic singlet-triplet gap with increasing length of

Adiabatic singlet-triplet energy gap ($E_{triplet}-E_{singlet}$)
for polyacenes. aug-cc-pVDZ basis set are used upto Hexacene and cc-pVDZ
basis set for the remainder. RAS(S)-1SF-PDFT results are obtained using tPBE
functional. Units are in kcal/mol.

Acene     2cCAS-1SF
2cCASSCF[using CAS(2,2)]
2cRAS(S)-1SF     expt.     expt.[1mm]

$bare$     -PDFT       $bare$     -PDFT       $bare$     -PDFT           $-$ZPE[ZPE correction
with B3LYP/6-31G(d,p) from Ref. ][2mm]
Naphthalene     52.7     57.4     64.7     76.8     56.2      62.3
61.1<cit.>     64.1, 64.3
Anthracene      31.4     36.7     48.7     51.9     33.8     40.4
43.1<cit.>     45.0, 45.4
Tetracene       26.4     26.6     41.8     29.1     27.7     28.5
29.5<cit.>     31.4
Pentacene       14.5     18.1     29.4     24.1     15.4     19.2      19.8<cit.>
20.9 [1mm]
Hexacene        10.8     12.2     23.2     13.7     11.5     13.3      12.5<cit.>
13.9 [1mm]
Heptacene       5.9     8.6     14.1     13.5     10.6      7.8        [1mm]
Octacene      3.4     5.6     10.1     6.4     10.0     4.8          [1mm]
Nonacene      1.3     4.1     5.0     6.7    7.1      3.3          [1mm]
Decacene      -0.2     1.9     3.2     2.3     5.6      2.0          [1mm]
Undecacene      -1.8     2.6     0.1     3.4     4.8      2.3          [1mm]
Dodecacene      -3.0     1.2     -1.1     0.9    3.9      2.1          [1mm]

§.§ Polyacenes

In the previous sections, we have shown that the RAS-$n$SF(-IP/EA)-PDFT
approach can describe challenging medium-sized biradical and triradical
molecular systems with good accuracy. In this section, we demonstrate the
performance of our approach for larger systems: the polyacenes, ranging from
naphthalene to dodecacene. The polyacenes have a singlet ground state, with
the open shell character increasing as the number of the benzene ring
increases. The singlet-triplet (ST) gap in this case correlates with the acene
length which exponentialy decreases.

Table <ref> reports the adiabatic singlet-triplet gaps of the
polyacenes obtained from RAS-1SF and RAS-1SF-PDFT. Here, the on-top tPBE
functional was employed. Using the triplet state
(m$_s$=1) as the reference, the singlet state is accessed by performing only a
single SF operation on the reference state. The external RAS(S) excitations
were carried out from all the valence $\pi$ orbitals.
The adiabatic ST gaps are plotted in Figure <ref> along with
selected literature values, namely, CCSD(T) extrapolated according to focal
point analysis<cit.>, Monte Carlo (JSD-LRDMC)<cit.>, v2RDM-CASSCF-PDFT<cit.> and GAS-PDFT (WFP-3
partitioning)<cit.>. Several other literature values are available, a detailed
analysis of the available ST gaps in literature are provided in
Ref.  and

The RAS(S)-1SF-PDFT adiabatic ST gaps have a good agreement with the corresponding
experimental values as compared to the RAS(S)-1SF values.
Although the bare RAS(S)-1SF includes some dynamical correlation<cit.>,
RAS(S)-1SF-PDFT improves over the bare RAS(S)-1SF method.

Here, a delicate balance between the dynamic and the static
correlation exist in the singlet and triplet state<cit.>. From Table
<ref>, it can be seen that the larger acenes agree very well with
the experimental value whereever available while for anthracene the difference is
$\sim$5 kcal/mol. However, we note that the experimental values were obtained
from measurements in solutions or as solids unlike the theoretical results
from gas phase calculations<cit.>.
Table <ref> also presents ST gaps obtained from bare
CAS-1SF, CASSCF(2,2) and the corresponding MC-PDFT approaches. In both the
cases, including the PDFT correction improves upon the bare method for the ST
gaps. However, the RAS(S)-1SF-PDFT performs better in comparison to CAS-1SF-PDFT
and CASSCF-PDFT for the polyacene chains. This shows that RAS(S)-1SF provides a more accurate 1- and
2-RDMs entering the MC-PDFT equation (Equation <ref>, <ref>
and <ref>) than CAS-1SF and CASSCF based on minal CAS(2,2) active space.

For the acenes where experimental values are not
available, the RAS-1SF-PDFT values are lower than the available literature values
shown in Figure <ref> although in some cases good agreements can be
seen. The CCSD(T) as well as the Monte Carlo values tend to overestimate the
available experimental values while it can be seen that the largest reported
acenes with these methods agrees well with the RAS-1SF-PDFT values.
The values reported for the
v2RDM-CASSCF-PDFT differ between 1.8 to 5.7 kcal/mol while the
difference is between 0.3 to 3.1 kcal/mol with the GAS-PDFT

Another feature that can be observed in Figure <ref> is the
smoothly decaying exponential curve for the RAS-1SF-PDFT ST gaps which can be
fitted to the form $aexp(-bx)+c$. The fit then can be used to estimate an
extrapolated ST gap for polyacene with infinite length. Using the fitting

\begin{equation}
E_{ST}(x)=138.14 e^{(-0.40 x)}+0.10
\end{equation}

The adiabatic ST gap estimated for an infinitely long linear polyacene chain
is 0.10 kcal/mol.
The value is in very good agreement with the ST gap of 0.18
kcal/mol from CCSD(T)/cc-pV$\infty$Z reported in Ref.  and with pp-RPA/cc-pVDZ with a
value between 0.0 and 2.3 kcal/mol reported in Ref. . The later however is
estimated using vertical ST gaps.
The ST gap is
directly related to the HOMO-LUMO gap and the obtained result suggest a
closure for the HOMO-LUMO gap for inifinitly long polyacene chains.
Note that
the result is in contrast to Ref.  where the ST gap for
the infinite chain is reported to be 4.87 kcal/mol with the v2RDM-CASSCF-PDFT/cc-pVTZ
method and Ref.  with a value of 5.06 and 5.37
kcal/mol obtained as their best estimate and from SF-CCSD/6-31+G(d,p) method respectively.
With GAS-PDFT/6-31+G(d,p) (WFP-3 partitioning), the ST gap is reported to be 1.9 kcal/mol
which is roughly half way between our value and the ones obtained from
v2RDM-CASSCF-PDFT and SF-CCSD methods<cit.>.

Although the theoretical results presented above show qualitative consistency,
the basis sets used to compute the ST gaps differ for the various approaches and so comparison with the available literature values does not account for basis set effects.

§ SUMMARY

In this work, we have presented a practical approach for treating large
multiconfigurational molecular system with a low computational cost.
The new method, spin-flip pair-density functional theory (SF-PDFT) uses
a spin-flip or a redox spin-flip operator to account for the static
correlation while the dynamical correlation is described with DFT using the
MC-PDFT approach.
The SF-PDFT method improves upon the result of bare spin-flip approach,
thereby capturing the missing dynamical correlations in the spin-flip
approach. In cases where the spin-flip results are already close to the
available experimental results, the SF-PDFT only changes the spin-flip results
slightly ascertaining the relaibility of SF-PDFT. The reason is because in
such cases, the static correlation dominates the electron correlation.
The method yields good accuracy for energy gaps between ground and the lowlying
excited states for challenging open-shell molecular systems. The applicability
range of the method was demonstrated by computing the singlet-triplet gap in
linear polyacene chains ranging from naphthalene to dodecacene. The SF-PDFT
predicts a closing singlet-triplet gap for an infinitly long linear polyacene

§ ACKNOWLEDGEMENTS

The authors are grateful for financial support provided by the U.S. Department
of Energy (Award No. DE-SC0018326) and computational infrastructure from the
Advanced Research Computing at Virginia Tech.
N. M. appreciates helpful discussions with Prof. Eugene DePrince III.