# MultiBinding Sites United in Covalent-Organic Frameworks (MSUCOF) for H2
Storage and Delivery at Room Temperature

Marcus Djokic  Department of Chemical Engineering and Material Science,
Michigan State University, East Lansing, MI, 48824, USA.  Jose L. Mendoza-
Cortes<EMAIL_ADDRESS>Department of Chemical Engineering and Material
Science, Michigan State University, East Lansing, MI, 48824, USA.

###### Abstract

The storage of hydrogen gas (H2) has presented a significant challenge that
has hindered its use as a fuel source for transportation. To meet the
Department of Energy’s ambitious goals of achieving $50$ g L-1 volumetric and
$6.5$ wt % gravimetric uptake targets, materials-based approaches are
essential. Designing materials that can efficiently store hydrogen gas
requires careful tuning of the interactions between the gaseous H2 and the
surface of the material. Metal-Organic Frameworks (MOFs) and Covalent-Organic
Frameworks (COFs) have emerged as promising materials due to their
exceptionally high surface areas and tunable structures that can improve gas-
framework interactions. However, weak binding enthalpies have limited the
success of many current candidates, which fail to achieve even $10$ g L-1
volumetric uptake at ambient temperatures. To overcome this challenge, We
utilized quantum mechanical (QM) based force fields (FF) to investigate the
uptake and binding enthalpies of 3 linkers chelated with 7 different
transition metals (TM), including both precious metals (Pd and Pt) and first
row TM (Co, Cu, Fe, Ni, Mn), to design 24 different COFs in-silico. By
applying QM-based FF with grand canonical Monte Carlo (GCMC) from 0-700 bar
and 298 K, We demonstrated that Co-, Ni-, Mn-, Fe-, Pd-, and Pt-based MSUCOFs
can already achieve the Department of Energy’s hydrogen storage targets for
2025. Surprisingly, the COFs that incorporated the more affordable and
abundant first-row TM often outperformed the precious metals. This promising
development brings us one step closer to realizing a hydrogen-based energy
economy.

###### keywords:

H2 storage | metalation | chelation | transition metals | covalent organic framework | imine | triazine | structure-property relations

##  1 INTRODUCTION

Storing hydrogen at operational temperatures poses a significant challenge to
its commercial viability, particularly in applications such as transportation.
The U.S. Department of Energy (DOE) has set ambitious targets for hydrogen
storage by 2025, including achieving 5.5 weight percent (wt %) and 40 g L-1 at
233-358 K, with an ultimate target of 6.5 wt % and 50 g L-1.1 Additionally,
during the DOE’s Hydrogen Shot summit in 2021, a further goal was announced:
to reduce the cost of clean hydrogen to $1 per 1 kg within the next decade.2
To achieve these targets, it is evident that materials-based approaches are
preferred. Rather than attempting to compress or cryogenically cool hydrogen
storage tanks, adsorbents such as Metal-Organic Frameworks (MOFs) and
Covalent-Organic Frameworks (COFs) can be introduced to the vessel, improving
gas uptakes without incurring high utility costs. Of these adsorbents, COFs
show particular promise due to their high surface area and composition of
lighter, more abundant elements compared to MOFs. 3, 4, 5

While many MOFs can offer high uptakes, operating at cryogenic temperature (77
K) swings to ambient conditions (298 K) is often prohibitively expensive. To
make the process more economical, it is preferable to perform the adsorption
process at ambient conditions. However, this is challenging for many MOFs and
COFs, as adsorption energies decay at higher temperatures due to entropic
changes. 6, 7 As a result, MOFs and COFs can lose their loading capacity at
higher temperatures, and those with promising uptakes at 77 K may not perform
as well at room temperature.

Efforts to improve gas uptake in porous materials such as MOFs and COFs have
led to the development of several strategies. One approach is reticular
chemistry, which allows for the selection of different nodes and linkers to
create customized structures. Building on this concept, Han et al.
investigated how various factors, including pore size, impregnation,
catenation, open metal sites, and functionalized linkers, affect hydrogen
storage in MOFs and COFs.3 In our study, we demonstrate that COFs with smaller
pores and tri-topic metalated linkers exhibit state-of-the-art performance. To
achieve optimal results, it is essential that metal sites are easily
accessible, enabling them to interact more efficiently with adsorbates. We
propose that by incorporating multiple metal sites into a linker, a greater
number of H2 molecules can interact directly with the sites via long-range
interactions without being shielded. Furthermore, cooperative interactions may
occur if these metalated sites are in close proximity. By exploring these
ideas, we aim to advance the field’s understanding of how customizable
structures and metalated linkers can enhance gas uptake in porous materials.

The conventional approach for predicting gas uptake in COFs and MOFs at
cryogenic temperatures (77 K) relies on pore volume, void fraction, and
surface area.8, 9 However, these features overlook crucial gas-framework
interactions, such as non-covalent and orbital interactions. To achieve
optimal hydrogen storage at room temperature and pressures up to 700 bar, it
is crucial to consider binding energies ($\Delta H_{bind}^{0}$) in the range
of 7-15 kJ/mol, which fall between chemi- and physisorption.10, 11, 6, 3 Since
COFs typically consist of many repeating binding sites, $\Delta H_{bind}^{0}$
is a good approximation to the experimentally obtained initial isosteric heat
of adsorption (Q${}_{st}^{0}$). Therefore, this heuristic of 7-15 kJ/mol
should also apply to the Q${}_{st}^{0}$. The results in this paper further
support this correlation that the 7-15 kJ/mol Q${}_{st}^{0}$ range results in
optimal performance.

##  2 MATERIALS AND METHODS

###  2.1 Design of Triazine COFs

Figure 1: New linkers with tri-topic connectivity for N-containing-COFs where
a high density of metalation sites are present: Dipyrazino [2,3-f:2’,3’-h]
quinoxaline (DPQ) and 2,4,6-tri(pyrimidin-2-yl)-1,3,5-triazine (TPT). The a)
pristine linkers are shown compared to b) the metalated case. Note that DPQ
sites can accommodate TMCl2, while TPT sites accommodate TMCl.

To create new frameworks with a high density of metalation sites and
cooperative effect, we introduced two linkers (Figure 1): Dipyrazino
[2,3-f:2’,3’-h] quinoxaline (DPQ) and 2,4,6-tri(pyrimidin-2-yl)-1,3,5-triazine
(TPT). These linkers were inspired by the well-known 4,4’-bipyridine (BPY) and
4,4’-diethylamido-2,2’-bipyridine (BPYAM) linkers and were chosen to increase
the density of metalation sites. To generate triazine-COFs, we constrained the
connectivity to a semi-regular net with a 3-4 connection. We selected only the
ctn ($I-43d$ space group) topologies for these frameworks, as they have been
shown to be the most stable.4, 12, 13, 14 We minimized the frameworks using
molecular mechanics without any symmetry constraints, and the coordinates are
reported in the SI. We note that related triazine frameworks have been
reported experimentally.15, 16, 17, 18

Figure 2: Approach suggested for the incorporation of the hydrazine
functionality with known boroxine and chelation chemistry to form a) MSUCOF-1
and MSUCOF-1-TMCl2, b) MSUCOF-2 and MSUCOF-2-TMCl2, or by in-situ triazine
formation to form c) MSUCOF-3 and MSUCOF-1-TMCl

To generate new N-containing COFs, we utilized four experimental building
blocks (1-4) as shown in Figure 2. These include methanetetrayltetrakis
(benzene-4,1-diyl) tetraboronic acid (1), dipyrazino [2,3-f:2’,3’-h]
quinoxaline-2,3,6,7,10,11-hexaol (2), diquinoxalino [2,3-a:2’,3’-c]
phenazine-2,3,8,9,14,15-hexaol (3), and 5,5’,5”,5”’-methanetetrayl-tetrakis
(pyrimidine-2-carbonitrile) (4). We used the known boronic chemistry to create
MSUCOF-1 (1+2) and MSUCOF-2 (1+3). 3 serves as an extended version of 2 with
the addition of benzene rings to each tri-topic end of 2, in effect creating a
MSUCOF-2 with a larger pore size, Psize, than MSUCOF-1 (as seen in Figure 3).
Linker 1 has already been used for imine COFs.19, 4 By using the new triazine
formation method, we design MSUCOF-3 (self-condensation of 4). Ultimately,
MSUCOF-1 and MSUCOF-2 have DPQ sites, while MSUCOF-3 has TPT sites. Upon metal
chelation, TM-Cl${}_{\text{x}}$ (where TM = Co, Cu, Fe, Mn, Ni, Pd, and Pt)
will occupy the DPQ and TPT sites to form metalated COFs using known methods.
The intention is to form interaction between the s-orbitals from H2 and the
d-orbitals from the TM.10

Figure 3: FF-Optimized structures of the pristine and metalated cases for a)
MSUCOF-1, b) MSUCOF-2, and c) MSUCOF-3.

###  2.2 Volumetric and Gravimetric uptake calculation

When analyzing porous materials, there are multiple methods to compare their
properties. Gas absorption is typically assessed by calculating volumetric
uptakes (g of H2 L-1) and gravimetric uptakes (g of H2 g-1). Gravimetric
uptakes are particularly important in vehicles since they measure the added
weight contribution of the porous material. Volumetric uptakes, on the other
hand, evaluate the energy density of hydrogen storage relative to the
container size. To accurately evaluate gas adsorption, it is essential to
distinguish between absolute, excess, and net adsorption. Absolute uptakes
measure the total number of H2 molecules loaded into the structure, whereas
excess uptakes quantify the loading with respect to H2 that would occupy the
pore volume if there were no COF/MOF present, as shown in Equation 1. 20

$\centering n^{exc.}=n^{abs.}-V_{P}*\rho^{exp.}_{H_{2}}\@add@centering$ (1)

Here, $n^{exc.}$ is the quantity of excess gas molecules, $n^{abs.}$ is the
absolute total quantity of gas molecules in the adsorbent system, $V_{p}$ is
the pore volume, and $\rho^{exp.}_{H_{2}}$ is the experimental hydrogen
density at the same temperature and pressure as the adsorbent system. In our
calculations, we used the NIST database to obtain hydrogen gas densities.21 A
major pitfall to excess uptakes is it fails to take into account the space
contribution of the adsorbent itself, which makes the adsorbent appear to be
more promising. A positive excess uptake may still result in a worse
performance than an tank full of hydrogen alone at the same pressure.
Additionally, excess uptake is directly sensitive to how the pore volume is
determined (e.g. geometric pore volume, probe-occupiable pore volume, van der
Waals volume, Connolly volume, and more)22, 23, 24, which pores are
accessible25, and whether the pore volume can change as is the case with
flexible MOFs/COFs 26, 27. Net uptakes, shown in Equation 2, try to address
these issues by measuring uptake with respect to the empty unit cell at the
same pressure, thereby directly assessing whether the adsorbent improves
overall storage capacity.

$\centering n^{net}=n^{abs.}-V_{U.C.}*\rho^{exp.}_{H_{2}}\@add@centering$ (2)

Replacing $V_{p}$ in Equation 1 with the unit cell volume of the MOF/COF,
denoted as $V_{U.C.}$, yields the net uptake, $n^{net}$. A positive net uptake
indicates an improvement in gaseous loading when the adsorbent is added to a
storage vessel. Nevertheless, net, absolute, and excess uptakes do not
consider the kinetics of the adsorption and desorption processes. An adsorbent
with high uptakes may not be practical if it cannot release hydrogen quickly
and inexpensively. To address this, usable working capacity (WC) quantifies
the amount of gas that can be withdrawn from the adsorbent in a reasonable
time frame. Specifically, WC can be calculated by subtracting the absolute
uptake at 5 bar, $n^{abs.}_{5}$, from the maximum uptake at the highest
pressure tested, $n^{abs.}_{700}$ (in this study, 700 bar).

$\centering WC=n^{abs.}_{700}-n^{abs.}_{5}\@add@centering$ (3)

###  2.3 QM-Based Force Fields

To accurately capture the interaction between molecular hydrogen and
transition metal chelation complexes (framework-TMCl${}_{\text{x}}$), a Morse
potential (shown in Equation 4) was fitted and used to parameterize the
interaction. The fitting process involved considering various first-row
elements (Co, Cu, Fe, Mn, and Ni) and late precious metals (Pd and Pt) using
density functional theory (DFT). The first row TM were previously fit11 using
B3LYP-D3 in Amsterdam Density Functional and GULP. In this paper, this same
methodology was used to expand to Pd and Pt, where the Morse fitted parameters
are summarized in Table 1. These parameters were incorporated as non-bonding
terms in the QM-modified force field, allowing for more accurate modeling of
the porous framework and gas interaction. It is important to note that the
Morse potential is specifically tailored to the adsorption of molecular
hydrogen and the transition metal chelation complexes in the framework,
ensuring a more precise representation of the system.

$\centering
U_{ij}^{\text{Morse}}(r_{ij})=D_{0}[(1-\text{e}^{-\alpha(r_{ij}-r_{0})})^{2}-1]\@add@centering$
(4)

Element | D0 (kcal/mol) |  $\boldsymbol{\alpha}$ (Å-1) | r0 (Å)
---|---|---|---
H2-Co | 0.879 | 0.850 | 2.985
H2-Cu | 0.818 | 1.462 | 2.931
H2-Fe | 1.092 | 1.180 | 3.155
H2-Mn | 0.994 | 0.990 | 3.015
H2-Ni | 1.154 | 1.210 | 3.207
H2-Pd | 0.641 | 0.974 | 3.268
H2-Pt | 1.652 | 0.920 | 3.282
H2-Cl | 0.146 | 1.446 | 3.725
Table 1: Table of fitted Morse potential force field parameters between
chelated atoms and H2

###  2.4 Grand Canonical Monte Carlo (GCMC)

In this study, we implemented the Ab-initio parameterized Force Fields in
Grand Canonical Monte Carlo (GCMC) simulations. To perform these simulations,
we employed the Metropolis algorithm code implemented in the Sorption Module
of Materials Studio28. Our investigations were conducted at 298 K and
pressures ranging from 1 bar to 700 bar. During each simulation step, hydrogen
gas underwent various processes in a predefined ratio (2:1:1:1:1:0.1) for
translation, rotation, creation, annihilation, regrowth, respectively. The
translation and rotation steps allowed the hydrogen atoms to freely move
within the rigid COF pores, while the creation and annihilation steps
facilitated the addition and removal of hydrogen molecules. The regrowth
process accommodated the conformational changes in asymmetric molecules, which
for molecular hydrogen, resulted in translation or rotation. Equilibrium was
achieved after 1,000,000 simulation steps, and the following 3,000,000 steps
were utilized to calculate the ensemble average loading and isosteric heat of
adsorption (Qst) for each temperature and pressure point.

###  2.5 Equation of State

$\displaystyle P=\frac{RT}{\nu-b}-\frac{a}{\nu^{2}}$
$ln\frac{f}{P}=(b-\frac{a}{RT})\frac{P}{RT}$ (5)

To make these results applicable to real-world data, an equation of state
(EOS) is required to convert the computed fugacity to pressure. The Van der
Waals29(Equation 5) and Peng-Robinson30 (Equation 6) equations of state (EOS)
were both considered. In each equation, empirical terms $a$ and $b$ are used,
with values of 0.2476 L2 bar/mol2 and 0.02661 L/mol, respectively, for
hydrogen. Other variables include $T$ for temperature, $\nu$ for molar volume,
$R$ for the gas constant, $P$ for pressure, $f$ for fugacity, and $Z$ for
compressibility factor. Note that for simplicity, terms A, B, and C are used
which represent various groupings of the previously mentioned variables into
single terms. Figure 4 demonstrates the differences between an ideal gas
(y=x), Van der Waals EOS (blue), and the Peng-Robinson EOS (orange). In the
low pressure regime <100 bar, Hydrogen gas will act like an ideal gas. At the
700 pressure range however, it strays drastically from the ideal pressures.
Between the two EOS, the calculated pressures showed a 0.735% difference
between the two EOS. Though both seem to be fairly interchangeable under these
conditions, we chose to use the Van der Waals EOS for fugacity and pressure
conversions.

$\displaystyle P=\frac{RT}{\nu-b}-\frac{a(T)}{\nu(\nu+b)+b(\nu-b)}$
$\displaystyle Z^{3}-(1-B)Z^{2}+(A-3B^{2}-2B)Z-(AB-B^{2}-B^{3})=0$
$\displaystyle A=\frac{aP}{R^{2}T^{2}},B=\frac{bP}{RT},C=\frac{P\nu}{RT}$
$ln\frac{f}{P}=Z-1-ln(Z-B)-\frac{A}{2\sqrt{2}B}ln(\frac{Z+2.414B}{Z-0.414B})$
(6)

Figure 4: Calculated Van der Waals (purple) and Peng-Robinson (green) fugacity
to pressure conversion showing good agreement up to the 700 bar pressure
evaluated. It is clear that these pressure and fugacity conversion stray from
the ideal gas approximation (light blue y=x line) above values of 100 bar for
H2 gas at 298 K.

##  3 RESULTS AND DISCUSSION

The H2 uptakes for the new N-triazine-COFs were calculated from 1 to 700 bar
at 298 K (Figure 5). The un-metalated COFs of this family give total
volumetric uptakes of 28 to 33 g L-1 while the metalated cases give uptakes of
at least 35 g L-1 at 700 bar. Of the three COF families tested, the MSUCOF-3
family gives the best performance in volumetric units of up to 73 g L-1, where
MSUCOF-3-PtCl2 obtains approximately triple the uptake of the pristine
MSUCOF-3. This is due to the high interaction of the MSUCOF-3 framework with
H2 while also having the smallest Psize of 7.7-9.4 Å, which avoids wasting
useful space. Congruently, this is close to the postulated range of the ideal
Psize for H2 storage of 6 to 12 Å. Thus these new COFs can overcome the 2025
DOE volumetric target of 40 g/L and gravimetric target of 5.5 weight %. The
excess gravimetric uptakes of the pure N-triazine-COFs are very similar to
other top-performing COFs. The metalated triazine-COFs perform as well as the
previously discovered promising COF-301-PdCl2 with total 60 g L-1 and 4.2
excess wt % uptakes at 100 bar and 298 K.31 In comparison, MSUCOF-3-PtCl2 at
100 bar and 298 K achieve a total uptake of 60.4 g L-1 and 4.06 excess wt %.

Figure 5: We show the high pressure H2 volumetric isotherms at 298 K for (a)
MSUCOF-1-TMCl2, (b) MSUCOF-2-TMCl2, and (c) MSUCOF-3-TMCl.

For the unmetalated triazine-COFs, we calculated the isosteric heat of
adsorption (Qst) to be around 3-4 kJ mol-1, as shown in Figure 6. However,
upon metalation, the Qst values increased significantly, ranging from 4-26 kJ
mol-1. Notably, the MSUCOF-3 family exhibited the highest Qst values, likely
due to its small pore size and multiple binding sites per pore. This trend is
reminiscent of the COF-301-PdCl2 case, where the pore size (Psize) was 6.3 Å
and the space between neighboring framework layers (Dlayers) was 9.9 Å. We
propose that the compact nature and multiple binding sites of MSUCOF-3 may
induce a cooperative effect, resulting from the overlap of the energy surfaces
of adjacent metalated sites.

A common trend observed in Qst values is that the maximum values are achieved
at lower pressures. As gas molecules are introduced into the structure at
higher pressures, the interactions between the framework and H2 become
obscured by the presence of other H2 molecules in the pore. The initial few
layers of adsorption can only accommodate a limited number of molecules,

Figure 6: We show the high pressure H2 isosteric heat of adsorption, Qst, at
298 K for (a) MSUCOF-1-TMCl2, (b) MSUCOF-2-TMCl2, and (c) MSUCOF-3-TMCl.

while the bulk pore can efficiently incorporate more hydrogen molecules. As a
result, the initial Qst value, Q${}_{st}^{0}$, at 1 bar gradually transitions
to a predominantly bulk-filled interaction as the pressure increases up to 700
bar (refer to Table 2).

To investigate whether other transition metals could perform similarly to Pd
or Pt, we calculated the H2 uptake for MSUCOF-1-TMCl2, MSUCOF-2-TMCl2, and
MSUCOF-3-TMCl2, where TM includes Co, Ni, Fe, Cu, Mn, Pd, and Pt.
Interestingly, several of the tested early transition metals exhibited high
uptakes that were comparable to those of Pd.

Figure 7: We show the high pressure H2 net volumetric isotherms, which
compares uptakes with respect to H2 loading in an empty unit cell of the same
volume, at 298 K for (a) MSUCOF-1-TMCl2, (b) MSUCOF-2-TMCl2, and (c)
MSUCOF-3-TMCl.

Among the early transition metals tested, Co, Ni, Mn, and Fe showed the best
performance across the new MSUCOF families. Notably, MSUCOF-3-CoCl exhibited
the highest volumetric uptake of 55.7 g L-1 at 700 bar, followed by
MSUCOF-3-PdCl with 47.8 g L-1, MSUCOF-1-CoCl2 with 39.7 g L-1, and
MSUCOF-1-PdCl2 with 37.8 g L-1. These results demonstrate that early
transition metals can achieve comparable performance to heavy transition
metals, which is particularly significant because

Figure 8: We show the high pressure H2 gravimetric isotherms at 298 K for (a)
MSUCOF-1-TMCl2, (b) MSUCOF-2-TMCl2, and (c) MSUCOF-3-TMCl.

early transition metals are relatively abundant and less expensive. In
particular, CoCl2, NiCl2, and FeCl2 chelation reactants can be obtained at
prices ($ mol-1) that are one to two orders of magnitude cheaper than PtCl2 or
PdCl2 (see SI Section Economic Analysis of MSUCOFs). Additionally, MnCl2 and
CuCl2 can be purchased at prices three orders of magnitude cheaper than the
precious transition metal chlorides, which is even more impressive.

To quantify the superiority of the COFs over an empty H2 tank, it is essential
to calculate the net volumetric uptakes. A positive value of net uptake
indicates that the COF can accommodate more hydrogen in the storage vessel
compared to bulk H2 at the same pressure. The peaks on the net volumetric
chart reveal the optimal pressure (Popt) at which the COF performs best.
Across the three families of COFs, the Popt occurs at approximately 100 bar,
suggesting that the performance diminishes at higher pressures. This is likely
due to a higher proportion of H2 bulk pore filling as adsorption sites begin
to saturate. Yet as illustrated in figure 7, MSUCOF-3-TMCl (TM = Pt, Co, Mn,
Ni, Fe, Pd) demonstrates extra H2 uptakes of 8.55-34.0 g L-1 at 700 bar, which
implies that these highlighted COFs remain exceptional even in applications up
to 700 bar.

Improving volumetric uptakes often comes at the expense of gravimetric
uptakes, as seen in many strategies. For example, COFs with smaller pores,
such as MSUCOF-3, tend to be denser than those with larger pores, such as
MSUCOF-1 and MSUCOF-2. Although the MSUCOF-3 family exhibits the most
impressive volumetric uptakes among these three, it also has the lowest
gravimetric uptake of up to 6.7 wt % (Figure 8). In contrast, the MSUCOF-2
family boasts an exceptional gravimetric uptake of up to 19.2 wt %. The
incorporation of TM through chelation significantly increases volumetric at
the expense of slightly diminished gravimetric uptakes. Heavier precious
metals add significant weight to the COFs with similar volumetric performance
as first row TM. Therefore, pristine COFs, as well as chelated COFs with the
most promising first-row TM, achieve the maximum gravimetric uptakes. This
further demonstrates how Co, Ni, Mn, and Fe may be more favorable than Pd and
Pt, as they often reach similar uptakes while contributing much less weight
and cost.

Table 2 reveals some interesting trends. Firstly, metalation of COFs often
resulted in up to 50% improvement in volumetric working capacities compared to
the pristine COFs. However, the addition of chelation complexes increased the
weight of COFs, resulting in a reduced gravimetric working capacity. Secondly,
COFs with Q${}_{st}^{0}$ values in the 7-15 kJ/mol range tend to exhibit the
best WC. For instance, although MSUCOF-3-PtCl demonstrated impressive absolute
and net uptakes, its WC was inferior to that of the three COF families tested,
owing to an overly strong interaction with hydrogen gas, which leads to a
large number of loaded molecules strongly binding to the COF, even at low
pressures. Third, the Qst tends to decrease to a bulk-filling interaction as
gas is loaded into the COF. COFs that maintain Qst within this ideal range
even at 700 bar exhibit the most impressive WC. Although MSUCOF-1-PtCl2
achieved DOE’s gravimetric and volumetric targets, several COFs in the
MSUCOF-3 family achieved higher volumetric WC. Importantly, common features
like S${}_{\text{A}}$ and V${}_{\text{P}}$ had little to no bearing on
predicting which COFs would perform the best.

|  |  |  | 1 bar | 700 bar |  | WC (700-5)
---|---|---|---|---|---|---|---
COF | SA | VP | Psize | Q${}_{\mathbf{st}}^{\mathbf{0}}$ | Qst | Popt |  H2 Upt. |  H2 Upt.
| (m2/g) | (cm3/g) | (Å) | (kJ/mol) | (kJ/mol) | (bar) | (g/L) | (wt %)
MSUCOF-1 | 5060 | 5.06 | 18.5 | 3.08 | 3.19 | - | 32.4 | 15.2
MSUCOF-1-CoCl2 | 4330 | 2.97 | 15.2 | 8.97 | 5.70 | 205 | 38.6 | 11.4
MSUCOF-1-CuCl2 | 4060 | 2.92 | 16.8 | 6.10 | 4.52 | 113 | 34.4 | 10.1
MSUCOF-1-FeCl2 | 4400 | 2.99 | 15.2 | 8.70 | 5.49 | 139 | 37.0 | 11.0
MSUCOF-1-MnCl2 | 4380 | 3.00 | 15.2 | 8.46 | 5.52 | 169 | 37.0 | 11.0
MSUCOF-1-NiCl2 | 4200 | 2.96 | 15.2 | 9.65 | 5.66 | 169 | 37.2 | 11.0
MSUCOF-1-PdCl2 | 3630 | 2.58 | 16.8 | 7.09 | 5.14 | 139 | 36.9 | 9.69
MSUCOF-1-PtCl2 | 2890 | 2.09 | 16.8 | 19.3 | 8.43 | 169 | 41.4 | 8.71
MSUCOF-2 | 4980 | 6.59 | 22.7 | 3.04 | 3.11 | 1.5 | 33.2 | 18.9
MSUCOF-2-CoCl2 | 4380 | 4.25 | 18.8 | 8.60 | 4.72 | 169 | 36.7 | 14.4
MSUCOF-2-CuCl2 | 3910 | 4.20 | 21.1 | 5.64 | 3.99 | 92 | 34.5 | 13.6
MSUCOF-2-FeCl2 | 4440 | 4.58 | 18.7 | 8.23 | 4.54 | 139 | 35.6 | 14.2
MSUCOF-2-MnCl2 | 4420 | 3.29 | 18.7 | 8.28 | 4.59 | 169 | 36.1 | 14.9
MSUCOF-2-NiCl2 | 4310 | 4.25 | 18.9 | 8.94 | 4.64 | 92 | 35.8 | 14.2
MSUCOF-2-PdCl2 | 3480 | 3.78 | 21.2 | 7.03 | 4.44 | 169 | 36.2 | 12.9
MSUCOF-2-PtCl2 | 2890 | 3.12 | 21.1 | 19.8 | 6.77 | 169 | 38.7 | 11.4
MSUCOF-3 | 5020 | 1.90 | 9.64 | 3.85 | 4.05 | - | 28.4 | 6.35
MSUCOF-3-CoCl | 3770 | 0.94 | 7.80 | 16.2 | 13.7 | 139 | 44.3 | 5.22
MSUCOF-3-CuCl | 3650 | 0.91 | 7.84 | 10.5 | 9.17 | 139 | 36.7 | 4.38
MSUCOF-3-FeCl | 3810 | 0.95 | 7.87 | 16.2 | 13.1 | 113 | 42.2 | 5.10
MSUCOF-3-MnCl | 3840 | 0.95 | 7.73 | 15.9 | 13.2 | 139 | 43.1 | 5.20
MSUCOF-3-NiCl | 3690 | 0.93 | 7.85 | 17.2 | 13.6 | 139 | 41.7 | 4.95
MSUCOF-3-PdCl | 2970 | 0.75 | 7.82 | 13.1 | 11.4 | 139 | 42.7 | 4.20
MSUCOF-3-PtCl | 2190 | 0.55 | 7.78 | 26.8 | 21.4 | 75 | 32.6 | 2.29
Table 2: Tabulated properties of the three MSUCOF families showing the
calculated: surface areas (SA), pore volume (VP), pore size (Psize), initial
isosteric heat of adsorption at 1 bar (Q${}_{\mathbf{st}}^{\mathbf{0}}$),
isosteric heat of adsorption at 700 bar (Qst), optimal pressure (Popt),
volumetric working capacity, and gravimetric working capacity.

##  4 SUMMARY.

The need for lightweight adsorbents that can operate under ambient conditions
is critical for hydrogen storage in fuel cell vehicles. COFs have emerged as a
promising candidate for this application, owing to their similar properties to
MOFs, while being composed of fewer heavy elements. By incorporating
transition metals into COFs and placing them in accessible areas within the
pore, their interaction with gaseous components can be optimized. By carefully
selecting linkers, researchers can generate COFs that are tailored for
hydrogen storage applications. Increasing the number of binding sites per
linker not only enhances the number of sites for H2 physisorption or
chemisorption, but also enables overlapping interactions between these sites,
resulting in a compounding gas-framework interaction. In this study, the
authors have developed 24 new promising COFs using tri-topic linkers chelated
with first-row transition metals and precious metals. MSUCOF-3-CoCl,
MSUCOF-3-FeCl, MSUCOF-3-MnCl, MSUCOF-3-NiCl, MSUCOF-3-PdCl, and MSUCOF-1-PtCl2
have all surpassed DOE’s 2025 40 g L-1 working capacity target at ambient
temperatures. Among these, MSUCOF-3-CoCl is particularly exceptional,
achieving an impressive 55.7 absolute g L-1 and 6.66 absolute wt % uptakes,
without the need for expensive Pd or Pt.

##  5 ACKNOWLEDGMENTS

This work was supported in part by computational resources and services
provided by the Institute for Cyber-Enabled Research at Michigan State
University. We also thank Prof. Srimanta Pakhira for initial helpful
discussions.

##  6 ASSOCIATED CONTENT

Supporting Information Available

##  7 AUTHOR INFORMATION

Corresponding Authors
<EMAIL_ADDRESS>
Notes
The authors declare no competing financial interests.

## References

  * DOE 2017 DOE DOE Technical Targets for Onboard Hydrogen Storage for Light-Duty Vehicles. 2017
  * Dinh 2022 Dinh, H. N. Hydrogen Production and Hydrogen Shot: Options for Producing Low-Carbon Hydrogen at Scale. 2022,
  * Han et al. 2009 Han, S. S.; Mendoza-Cortés, J. L.; Iii, W. A. G. Recent advances on simulation and theory of hydrogen storage in metal–organic frameworks and covalent organic frameworks. 2009, _38_ , 1460–1476
  * Geng et al. 2020 Geng, K.; He, T.; Liu, R.; Dalapati, S.; Tan, K. T.; Li, Z.; Tao, S.; Gong, Y.; Jiang, Q.; Jiang, D. Covalent Organic Frameworks: Design, Synthesis, and Functions. 2020, _120_ , 8814–8933
  * Furukawa and Yaghi 2009 Furukawa, H.; Yaghi, O. M. Storage of Hydrogen, Methane, and Carbon Dioxide in Highly Porous Covalent Organic Frameworks for Clean Energy Applications. 2009, _131_ , 8875–8883
  * Bhatia and Myers 2006 Bhatia, S. K.; Myers, A. L. Optimum Conditions for Adsorptive Storage. 2006, _22_ , 1688–1700
  * Han and Goddard 2008 Han, S. S.; Goddard, W. A. I. High H2 Storage of Hexagonal Metal-Organic Frameworks from First-Principles-Based Grand Canonical Monte Carlo Simulations. 2008, _112_ , 13431–13436
  * Ahmed et al. 2019 Ahmed, A.; Seth, S.; Purewal, J.; Wong-Foy, A. G.; Veenstra, M.; Matzger, A. J.; Siegel, D. J. Exceptional hydrogen storage achieved by screening nearly half a million metal-organic frameworks. 2019, _10_ , 1568
  * Ahmed and Siegel 2021 Ahmed, A.; Siegel, D. J. Predicting hydrogen storage in MOFs via machine learning. 2021, _2_ , 100291
  * Kubas 2007 Kubas, G. J. Fundamentals of H2 Binding and Reactivity on Transition Metals Underlying Hydrogenase Function and H2 Production and Storage. 2007, _107_ , 4152–4205
  * Pramudya and Mendoza-Cortes 2016 Pramudya, Y.; Mendoza-Cortes, J. L. Design Principles for High H2 Storage Using Chelation of Abundant Transition Metals in Covalent Organic Frameworks for 0–700 bar at 298 K. 2016, _138_ , 15204–15213
  * El-Kaderi et al. 2007 El-Kaderi, H. M.; Hunt, J. R.; Mendoza-Cortés, J. L.; Côté, A. P.; Taylor, R. E.; O’Keeffe, M.; Yaghi, O. M. Designed Synthesis of 3D Covalent Organic Frameworks. 2007, _316_ , 268–272
  * Delgado-Friedrichs et al. 2006 Delgado-Friedrichs, O.; O’Keeffe, M.; Yaghi, O. M. Three-periodic nets and tilings: edge-transitive binodal structures. 2006, _62_ , 350–355
  * Schmid and Tafipolsky 2008 Schmid, R.; Tafipolsky, M. An accurate force field model for the strain energy analysis of the covalent organic framework COF-102. 2008, _130_ , 12600–12601
  * Kuhn et al. 2008 Kuhn, P.; Antonietti, M.; Thomas, A. Porous, Covalent Triazine-Based Frameworks Prepared by Ionothermal Synthesis. 2008, _47_ , 3450–3453
  * Singh and Nagaraja 2022 Singh, G.; Nagaraja, C. Rational design of Cu(I)-anchored porous covalent triazine framework (CTF) for simultaneous capture and conversion of CO2 at ambient conditions. 2022, _63_ , 102132
  * Sun et al. 2022 Sun, R.; Wang, X.; Wang, X.; Tan, B. Three-Dimensional Crystalline Covalent Triazine Frameworks via a Polycondensation Approach. 2022, _61_ , e202117668
  * Lin and EL-Mahdy 2022 Lin, K. Y.; EL-Mahdy, A. F. Covalent triazine frameworks based on triphenylpyridine building block for high-performance supercapacitor and selective CO2 capture. 2022, _281_ , 125850
  * Waller et al. 2015 Waller, P. J.; Gándara, F.; Yaghi, O. M. Chemistry of Covalent Organic Frameworks. 2015, _48_ , 3053–3063
  * Brandani et al. 2016 Brandani, S.; Mangano, E.; Sarkisov, L. Net, excess and absolute adsorption and adsorption of helium. 2016, _22_ , 261–276
  * Linstrom and Mallard 2001 Linstrom, P. J.; Mallard, W. G. The NIST Chemistry WebBook: A Chemical Data Resource on the Internet. 2001, _46_ , 1059–1063
  * Ongari et al. 2017 Ongari, D.; Boyd, P. G.; Barthel, S.; Witman, M.; Haranczyk, M.; Smit, B. Accurate Characterization of the Pore Volume in Microporous Crystalline Materials. 2017, _33_ , 14529–14538
  * Herrera et al. 2012 Herrera, L. F.; Fan, C.; Nguyen, V.; Do, D. D.; Horikawa, T.; Nicholson, D. A self-consistent method to determine accessible volume, area and pore size distribution (APSD) of BPL, Norit and AX-21 activated carbon. 2012, _50_ , 500–509
  * Connolly 1985 Connolly, M. L. Computation of molecular volume. 1985, _107_ , 1118–1124
  * Van Heest et al. 2012 Van Heest, T.; Teich-McGoldrick, S. L.; Greathouse, J. A.; Allendorf, M. D.; Sholl, D. S. Identification of Metal–Organic Framework Materials for Adsorption Separation of Rare Gases: Applicability of Ideal Adsorbed Solution Theory (IAST) and Effects of Inaccessible Framework Regions. 2012, _116_ , 13183–13195
  * Schneemann et al. 2014 Schneemann, A.; Bon, V.; Schwedler, I. I.; Senkovska, I.; Kaskel, S.; Fischer, R. A. Flexible metal–organic frameworks. 2014, _43_ , 6062–6096
  * Sakamoto et al. 2022 Sakamoto, H.; Ito, A.; Ohtani, M. Unusual ligand substitution of a metal–organic framework with distorted metal–ligand coordination. 2022, _24_ , 1690–1694
  * BIOVIA 2022 BIOVIA, D. S. Materials Studio 2022. 2022
  * Winn 1988 Winn, J. S. The fugacity of a van der Waals gas. _Journal of Chemical Education_ 1988, _65_ , 772
  * Peng and Robinson 1976 Peng, D.-Y.; Robinson, D. B. A New Two-Constant Equation of State. 1976, _15_ , 59–64
  * Mendoza-Cortes et al. 2012 Mendoza-Cortes, J. L.; Goddard, W. A. I.; Furukawa, H.; Yaghi, O. M. A Covalent Organic Framework that Exceeds the DOE 2015 Volumetric Target for H2 Uptake at 298 K. _The Journal of Physical Chemistry Letters_ 2012, _3_

See pages - of
SI_Multi_Metallic_Multibinding_Sites_in_Covalent_Organic_Frameworks_McCOF.pdf


$\displaystyle=\left\langle\widetilde{V},\widetilde{D}-\widetilde{V}\right\rangle=2\left\langle\widetilde{V}_{\alpha\gamma},\widetilde{D}_{\alpha\gamma}-\widetilde{V}_{\alpha\gamma}\right\rangle+\left\langle
Z(\widetilde{D}_{\beta\beta}),\widetilde{D}_{\beta\beta}-Z(\widetilde{D}_{\beta\beta})\right\rangle$
$\displaystyle\geq
2\left\langle\widetilde{V}_{\alpha\gamma},\widetilde{D}_{\alpha\gamma}-\widetilde{V}_{\alpha\gamma}\right\rangle=-2\sum\limits_{i\in\alpha,j\in\gamma}\frac{\lambda_{j}}{\lambda_{i}}\left(\widetilde{V}_{ij}\right)^{2}.$
(82)

On one hand, by picking $\alpha^{\prime}=\beta$ in Lemma 11 of [30], we choose
$\widehat{Z}\in\partial\operatorname{\Pi}_{S^{|\beta|}_{+}}(0)$ such that
$\widetilde{V}_{\beta\beta}=\widehat{Z}(\widetilde{D}_{\beta\beta})=\widetilde{D}_{\beta\beta}$
and the equality can be attained in the inequality of (A.2). On the other
hand, observing Proposition 3.3 of [8], we have $\Xi_{1}=E$ therein and for
any $Y\in S^{|\beta|}$, $0\in D^{*}\mathcal{N}_{S^{|\beta|}_{+}}(0,0)(Y)$,
which implies that $0\in
D^{*}\mathcal{N}_{S^{|\beta|}_{+}}(0,0)(-\widetilde{V}_{\beta\beta})$ and

$\displaystyle
P\left(\begin{array}[]{ccc}0&0&\widetilde{V}_{\alpha\gamma}-\widehat{\operatorname{\Sigma}}_{\alpha\gamma}\circ\widetilde{V}_{\alpha\gamma}\\\
0&0&\widetilde{V}_{\beta\gamma}\\\
\widetilde{V}_{\gamma\alpha}-\widehat{\operatorname{\Sigma}}_{\gamma\alpha}\circ\widetilde{V}_{\gamma\alpha}&\widetilde{V}_{\gamma\beta}&\widetilde{V}_{\gamma\gamma}\end{array}\right)P^{*}\in
D^{*}\mathcal{N}_{S^{m}_{+}}(X,U)(-V),$

where $\widehat{\operatorname{\Sigma}}\in S^{m}$ and

$\displaystyle\widehat{\operatorname{\Sigma}}_{ij}:=\frac{|\lambda_{i}(A)|+|\lambda_{j}(A)|}{\max\\{\lambda_{i}(A),0\\}+\max\\{\lambda_{j}(A),0\\}},\qquad
i\in\alpha,\,j\in\gamma.$

Therefore, since

$\displaystyle-Z\in
D^{*}\mathcal{N}_{S^{m}_{+}}(X,U)(-V)\quad\Longleftrightarrow\quad V\in
D^{*}\operatorname{\Pi}_{S^{m}_{+}}(X+U)(Z+V),$

we can conclude that

$\displaystyle
P\left(\begin{array}[]{ccc}\widetilde{V}_{\alpha\alpha}&\widetilde{V}_{\alpha\beta}&\widehat{\operatorname{\Sigma}}_{\alpha\gamma}\circ\widetilde{V}_{\alpha\gamma}\\\
\widetilde{V}_{\beta\alpha}&\widetilde{V}_{\beta\beta}&0\\\
\widehat{\operatorname{\Sigma}}_{\gamma\alpha}\circ\widetilde{V}_{\gamma\alpha}&0&0\end{array}\right)P^{*}\in\mathop{\operatorname{argmin}}_{{D\in
S^{m},V=WD}\atop{W\in\partial\operatorname{\Pi}_{S^{m}_{+}}(X+U)}}\langle
V,D-V\rangle\cap\Big{\\{}D\,\Big{|}\,V\in
D^{*}\operatorname{\Pi}_{S^{m}_{+}}(X+U)(D)\Big{\\}}.$

The desired result follows.

## References

  * [1] F. J. Aragón Artacho and M. H. Geoffroy, Characterization of metric regularity of subdifferentials, Journal of Convex Analysis, 15 (2008), pp. 365–380.
  * [2] J.-P. Aubin, Contingent derivatives of set-valued maps and existence of solutions to nonlinear inclusions and differential inclusions, in Nachbin, L. (ed.) Mathematical Analysis and Applications, Part A, Academic Press, New York, 1981, pp. 160–232.
  * [3] J. F. Bonnans and H. Ramírez C., Perturbation analysis of second-order cone programming problems, Mathematical Programming, 104 (2005), pp. 205–227.
  * [4] J. F. Bonnans and A. Shapiro, Perturbation Analysis of Optimization Problems, Springer, New York, 2000.
  * [5] F. H. Clarke, On the inverse function theorem, Pacific Journal of Mathematics, 64 (1976), pp. 97–102.
  * [6] , Optimization and Nonsmooth Analysis, Second Edition, Classics in Applied Mathematics, 5, Society for Industrial and Applied Mathematics, Philadelphia, PA, 1990.
  * [7] Y. Cui, D. Sun, and K.-C. Toh, On the R-superlinear convergence of the KKT residuals generated by the augmented Lagrangian method for convex composite conic programming, Mathematical Programming, (2019), pp. 381–415.
  * [8] C. Ding, D. Sun, and J. J. Ye, First order optimality conditions for mathematical programs with semidefinite cone complementarity constraints, Mathematical Programming, 147 (2014), pp. 539–379.
  * [9] A. L. Dontchev and R. T. Rockafellar, Implicit Functions and Solution Mappings-A View from Variational Analysis, Springer, New York, 2009.
  * [10] H. Gfrerer and J. V. Outrata, On Lipschitzian properties of implicit multifunctions, SIAM Journal on Optimization, 26 (2016), pp. 2160–2189.
  * [11] H. Gfrerer and J. V. Outrata, On (local) analysis of multifunctions via subspaces contained in graphs of generalized derivatives, Journal of Mathematical Analysis and Applications, 508 (2022), p. 125895.
  * [12] M. R. Hestenes, Multiplier and gradient methods, Journal of Optimization Theory and Applications, 4 (1969), pp. 303–320.
  * [13] T. Hoheisel, C. Kanzow, B. S. Mordukhovich, and H. Phan, Generalized Newton’s method based on graphical derivatives, Nonlinear Analysis: Theory, Methods $\&$ Applications, 75 (2012), pp. 1324–1340.
  * [14] A. D. Ioffe, Approximate subdifferentials and applications. I: The finite-dimensional theory, Transactions of the American Mathematical Society, 281 (1984), pp. 389–416.
  * [15] A. D. Ioffe and J. V. Outrata, On metric and calmness qualification conditions in subdifferential calculus, Set-Valued Analysis, (2008), pp. 199–227.
  * [16] D. Klatte and B. Kummer, Aubin property and uniqueness of solutions in cone constrained optimization, Mathematical Methods of Operations Research, 77 (2013), pp. 291–304.
  * [17] B. Kummer, Lipschitzian inverse functions, directional derivatives, and applications in $C^{1,1}$ optimization, Journal of Optimization Theory and Applications, 70 (1991), pp. 561–582.
  * [18] Y. Liu, Y. Sun, and S. Pan, Computation of graphical derivatives of normal cone maps to a class of conic constraint sets, Set-Valued and Variational Analysis, 27 (2019), pp. 783–806.
  * [19] F. J. Luque, Asymptotic convergence analysis of the proximal point algorithm, SIAM Journal on Control and Optimization, 22 (1984), pp. 277–293.
  * [20] A. Mohammadi, B. S. Mordukhovich, and M. E. Sarabi, Parabolic regularity in geometric variational analysis, Transactions of the American Mathematical Society, 374 (2021), pp. 1711–1763.
  * [21] , Variational analysis of composite models with applications to continuous optimization, Mathematics of Operations Research, 47 (2022), pp. 397–426.
  * [22] B. S. Mordukhovich, Maximum principle in problems of time optimal control with nonsmooth constraints, Journal of Applied Mathematics and Mechanics, 40 (1976), pp. 960–969.
  * [23] , Metric approximations and necessary optimality conditions for general classes of extremal problems, Soviet Mathematics Doklady, 22 (1980), pp. 526–530.
  * [24] B. S. Mordukhovich, Stability theory for parametric generalized equations and variational inequalities via nonsmooth analysis, Transactions of the American Mathematical Society, 343 (1994), pp. 609–657.
  * [25] B. S. Mordukhovich, Variational Analysis and Generalized Differentiation, I: Basic Theory, II: Applications, Springer, Berlin, 2006.
  * [26] , Variational Analysis and Applications, Springer, Cham, Switzerland, 2018.
  * [27] B. S. Mordukhovich, T. T. A. Nghia, and R. T. Rockafellar, Full stability in finite-dimensional optimization, Mathematics of Operations Research, 40 (2015), pp. 226–252.
  * [28] B. S. Mordukhovich and M. E. Sarabi, Stability analysis for composite optimization problems and parametric variational systems, Journal of Optimization Theory and Applications, 172 (2017), pp. 554–577.
  * [29] J. V. Outrata and H. Ramiriz C., On the Aubin property of perturbed second-order cone programs, SIAM Journal on Optimization, 21 (2011), pp. 798–823.
  * [30] J.-S. Pang, D. Sun, and J. Sun, Semismooth homeomorphisms and strong stability of semidefinite and lorentz complementarity problems, Mathematics of Operations Research, 28 (2003), pp. 39–63.
  * [31] R. A. Poliquin and R. T. Rockafellar, Tilt stability of a local minimum, SIAM Journal on Optimization, 8 (1998), pp. 287–299.
  * [32] M. J. D. Powell, A method for nonlinear constraints in minimization problems, in Optimization, Fletcher, R. (ed.), Academic Press, New York, 1969, pp. 283–298.
  * [33] L. Q. Qi and J. Sun, A nonsmooth version of Newton’s method, Mathematical Programming, 58 (1993), pp. 353–367.
  * [34] S. M. Robinson, Strongly regular generalized equations, Mathematics of Operations Research, 5 (1980), pp. 43–62.
  * [35] R. T. Rockafellar, Convex Analysis, Princeton University Press, 1970\.
  * [36] , Monotone operators associated with saddle functions and minimax problems, Nonlinear Functional Analysis, Part 1, F. E. Browder, ed., Symposia in Pure Math., 18 (1970), pp. 397–407.
  * [37] , Augmented Lagrangians and applications of the proximal point algorithm in convex programming, Mathematics of Operations Research, 1 (1976), pp. 97–116.
  * [38] , Monotone operators and the proximal point algorithm, SIAM Journal on Control and Optimization, 14 (1976), pp. 877–898.
  * [39] , Generalized second derivatives of convex functions and saddle functions, Transactions of the American Mathematical Society, 322 (1990), pp. 51–77.
  * [40] R. T. Rockafellar and R.-B. Wets, Variational Analysis, Springer, New York, 2009.
  * [41] D. Sun, The strong second order sufficient condition and constraint nondegeneracy in nonlinear semidefinite programming and their applications, Mathematics of Operations Research, 31 (2006), pp. 761–776.
  * [42] X.-Y. Zhao, D. Sun, and K.-C. Toh, A Newton-CG augmented Lagrangian method for semidefinite programming, SIAM Journal on Optimization, 20 (2010), pp. 1737–1765.


Efficient compilation of expressive problem space specifications to neural network solvers Supported by organization x.

Efficient compilation to neural network solvers

Matthew L. Daggitt 1 2

Wen Kokke

Robert Atkey

M.L. Daggitt et al.

Heriot-Watt University, Edinburgh, UK
University of Western Australia, Perth, Australia
University of Strathclyde, Glasgow, UK

Recent work has described the presence of the embedding gap in neural network verification.
On one side of the gap is a high-level specification about the network's behaviour, written by a domain expert in terms of the interpretable problem space.
On the other side are a logically-equivalent set of satisfiability queries, expressed in the uninterpretable embedding space in a form suitable for neural network solvers.
In this paper we describe an algorithm for compiling the former to the latter.
We explore and overcome complications that arise from targeting neural network solvers as opposed to standard SMT solvers.

§ INTRODUCTION

Neural network verification is a fast growing field and in recent years there has been a rapid increase in the number of neural network solvers.
These solvers are capable of checking whether there exists a satisfying assignment to a given set of linear constraints over the inputs and outputs of neural networks [6, 7, 8, 11, 14].
However, as with the SMTLib format for traditional SMT solvers, the input formats for these solvers are low-level and are designed to be read by machines rather than humans [1].
To address this, various tools have proposed domain specific languages (DSLs), which allow users to writer their specifications at a higher level of abstraction.
These specifications are then compiled down to low-level queries for solvers.
For example, the DNNV framework [13] uses a Python DSL, and the Caisar framework [5] uses the WhyML language [4].

There are still several ways the expressive power of the current generation of DSLs can be improved. They currently cannot express properties containing multiple applications of the same network (e.g. monotonicity of a network [9]) or properties involving two different networks (e.g. approximate equality of teacher-student networks [12]). Furthermore, they don't support conditional statements or arbitrary numbers of quantifiers.
However, we believe the most crucial missing capability is that currently specifications can only be written in terms of the direct inputs and outputs to the network. We will now elaborate further on this last point.

Recall that in a typical neural network setup, an unknown function ${\mathcal{H}: P \rightarrow Q}$ is approximated by $u \circ f \circ e$, where $e : P \rightarrow \real^m$ is the embedding function, $f : \real^m \rightarrow \real^n$ is the neural network and $u : \real^n \rightarrow Q$ is the unembedding function.
By design, neural network solver queries can only refer to the $m$ input and $n$ output variables in the embedding space, $\real^m \times \real^n$.
However, these vectors of numbers are often uninterpretable to users as the functions $e$ and $u$ usually remove all structure present in the problem space, $P \times Q$.
Consequently, specifications written in terms of these input and output variables are not easily interpretable in terms of the original problem space.
We argue that DSLs should allow domain experts to write specifications in the problem space, i.e. $P \times Q$.

Let us look at an example. Consider a neural network trained to control the brakes of a car and a specification “if the distance to the object is more than 100m the force applied to the brakes should be less than 250 newtons”.
As stated above, values are expressed in SI units in the problem space.
However if the embedding functions map inputs and outputs values to and from the range [0,1] and the domain expert is forced by the DSL to write the specification in the embedding space, then the specification might become “if the distance to the object is more than 0.7514 the force applied to the brakes should be less than 0.1123”. This new version is no longer as interpretable as the original version, as the reader cannot easily map the values in the embedding space back to SI units in the problem space.

§.§ Contributions

In this paper we describe a DSL and an efficient compilation algorithm which addresses the issues described: specifications are allowed to refer to an arbitrary number of networks, network applications, and quantifiers and supports conditional statements and higher-order functions.
Furthermore, quantified variables can be first arbitrarily linearly transformed before being used as inputs to neural networks, allowing them to represent values in the problem space rather than the embedding space.
We show the algorithm scales optimally for a selection of common queries, even when $m$ and $n$ > 1000 as is common in image processing.

The paper proceeds as follows. In Section <ref> we describe a query language for neural network solvers and our high-level DSL for writing specifications.
In Section <ref> we describe the compilation algorithm itself.

§ BACKGROUND

§.§ Query language

A generic query language for neural network solvers is shown in Figure <ref>. Given a list of $l$ neural networks $f_1 : \real^{m_1} \rightarrow \real^{n_1}, f_2 : \real^{m_2} \rightarrow \real^{n_2}, \ldots, f_l : \real^{m_l} \rightarrow \real^{n_l}$, a single query is a tree of conjunctions and disjunctions with assertions at the leaves.
An assertion is pair of a linear expression and a constant related either by an equality or an inequality.
The variables in the linear expression must be of the form $\lInputVar{i}$ for $i < \sum_{k = 1}^l m_k$ or $\lOutputVar{j}$ for $j < \sum_{k = 1}^l n_k$.
Concretely if $i = (\sum_{k = 0}^{a} m_k) + b$ for some $a \in [1, l]$ and $b \in [1, m_a]$ then $X_i$ represents the $a^{th}$ input of the $b^{th}$ network.
Likewise if $i = (\sum_{k = 0}^{a} n_k) + b$ for some $a \in [1, l]$ and $b \in [1, n_a]$ then $Y_i$ represents the $a^{th}$ output of the $b^{th}$ network.
For example, given two networks $f_1 : \real^{2} \rightarrow \real, f_2 : \real \rightarrow \real^{2}$ then $X_1$, $X_2$ and $Y_1$ are inputs and outputs of $f_1$, and $X_3$, $Y_2$ and $Y_3$ are inputs and outputs of the $f_2$.

Syntax for a generic neural network solver query language.

Note that when representing queries about multiple applications of the same network with different inputs (e.g. in monotonicity specifications), $f_c$ and $f_d$ are not necessarily distinct even when $c \neq d$.
This language is a generalisation of the standard VNNLib query format [1] supported by a wide range of solvers.
There are two main differences: firstly readability has been improved by making operators infix, and secondly, support for queries over multiple networks has been added.

In practice very few solvers support queries with arbitrary interleavings of disjunctions and conjunctions.
Instead most require that queries must be converted to disjunctive normal form (DNF) and submitted as multiple separate queries.
Nonetheless in this paper we target the non-DNF form above, as a) it is a generalisation of DNF and therefore DNF is trivial to obtain from it, b) maintaining the disjuncts is actually a harder compilation problem than doing everything in DNF, c) it avoids an exponential increase in the size of the queries in the worst case, and d) we hypothesis that some solvers may be able to reuse intermediate results between different disjunctive branches.

§.§ Specification language

Figure <ref> shows a typed lambda calculus, with builtin operations for logic, arithmetic and size-safe vector operations, for writing high-level specifications.

Syntax for the computational fragment of the Vehicle DSL

In this language a specification is a list of two types of declarations:  declarations which bring a new network into scope and  declarations which define a boolean expression who's veracity should be checked by queries to a neural network solver.
At the type-level it has a sized index and vector types, the latter of which is polymorphic in its element type, as well as  and .

Available vector-level operations include the lookup operator , where  is the value at position $i$ in vector $x$ and higher-order functions such as ,  and  which can be used to define further operations over vectors.
The language also includes vector addition $\oAdd$ and vector equality $\oEq$ as builtins, despite the fact that they are definable using a combination of polymorphism and type-classes and the operations  and .
The reason for this will be made clear in Section <ref>.

We emphasise this is a core language. A usable specification language should contain quality of life features such as auxiliary definitions,  and  declarations, polymorphism (allowing the same function to work over different types, e.g.  and ), type-classes (allowing the same symbol to represent different functions, e.g. `` can represent addition over both numbers and vectors).
However, we assume these features are dealt with in compiler passes prior to compilation to queries and therefore they are not considered in this paper.
We also assume the presence of a sound and complete typing relation.

§.§ What is not compilable

It is crucial to note that the addition of a high-level DSL on top of a neural network solver does not increase the classes of problems that can be solved.
Therefore the following classes of specifications cannot be compiled. In all the examples we assume a single network declaration:
@network f : Vector Real 1 -> Vector Real 1

Non-linear relationships between problem space and embedding space variables - in general query language's only support linear constraints, and therefore example1 is not compilable as it requires solving the non-linear relationship $X_1 = a * a * a$ between the problem space and embedding space variables.
@property example1 = forall a . f [a * a * a] ! 0 >= 0

Alternating quantifiers - example2 is not compilable as it has alternating quantifiers in it, and therefore cannot be reduced to an equi-satisfiable set of queries involving only implicit existentials.
@property example2 = forall (a : Real) . exists (b : Real) . f [b] ! 0 == a

Quantifying over functions - example3 is not compilable as it involves quantifying over functions, something not supported .
@property example3 = forall (h : Real -> Real) . f [h x] ! 0 == 0

Recent work has showed how specifications with non-linear constraints or alternating quantifiers can be detected and detailed explanations of why they are uncompilable can be generated via an evidence producing type-system [2].

§.§ What should be compilable

Conceptually the approach to compilation is simple: normalise the specification while replacing network applications with embedding space variables and solving for user variables.
We now present a set of examples, many of which are minimised versions of real-world specifications, that illustrate why it is challenging to do this in general for such an expressive DSL.
All examples assume the presence of two network declarations:
@network f : Vector Rat 1 -> Vector Rat 1
@network g : Vector Rat 2 -> Vector Rat 1

Multiple nested network applications - multiple network applications may be nested arbitrarily deeply and the same network may be applied multiple times.
@property example4 =
forall (a : Real). g [0.5 + f [a] ! 0 , g [0, a] ! 0] ! 0 >= 0

Hidden boolean structure in assertions - if statements allow the embedding of arbitrary boolean structure, including quantifiers inside what appear to be simple assertions.
@property example5 =
0.5 <= (if (forall (a : Real) . f [a] >= 1) then 0 else 1)

Multiply determined and under-determined user variables - disjunctions may result in there not being a single  consistent solution for the user quantified variables in terms of the network variables.
Furthermore, in the right-hand side of example3's disjunct, the user variables are undetermined as there is only one network input variable, but two user variables that need to be eliminated.
@property example6 =
exists (a b : Real). g [a , b] >= 5 or f [ a + b ] >= 0

Quantified vector variables that block evaluation of higher-order functions - vector variables may be quantified over and these variables can block the evaluation of higher order functions such as  and .
@property example7 =
forall (as : Vector Real 10). fold (r -> f a >= 0 and r) false as

§ COMPILATION ALGORITHM

In this section, we describe an algorithm for compiling a  in the specification language in <ref> to a set of queries in the query language in Figure <ref>.
The algorithm is implemented in Haskell [10], and while we restrict ourselves to using basic syntax as much as possible, the use of monads to cleanly represent global state is unavoidable.
For readers not well-versed in the dark arts of monadic programming, MonadWriter can be though of a write-only global variables, and MonadState can be thought of as read-and-write global variables.

§.§ Variables

For efficiency of normalisation, variables are represented using de Bruijn levels, and referred to using the Haskell type Level.
For clarity, we introduce the following type synonyms for the various categories of variables:
type UserRealVar = Level    type NetworkRealVar = Level
data UserVectorVar = Level  type NetworkVectorVar = Level
type RealVar = Level        type VectorVar = Level    type AnyVar = Level

The the global state of the algorithm is represented by a Ctx object and contains i) the stack of variables current in scope, where Levels are used to index from the bottom of the stack and ii) a mapping between vector variables and a list of variables that represent the elements of that vector:
data Ctx = Context
boundVars :: Stack AnyVar
, vectorVars :: Map VectorVar [RealVar]

As we will discuss in Section <ref>, maintaining variables for both a vector and its elements will be crucial for getting efficient performance.
We therefore define the following helper function which introduces a new set of variables representing a vector of size $n$ and its individual elements to the context and returns the vector variable:
addVectorVariableToCtx :: MonadState Ctx m => Int -> m Value
addVectorVariableToCtx n = do
GlobalCtx .. <- get
let l = length boundVars
let v = Var l []
let v1, ..., vn = Var (l + 1) [], ..., Var (l + n) []
set (GlobalCtx
boundVars = [vn...,v1,v] <> boundVars
, elementVars = v -> [vn,...,v1] <> elementVars
return a

§.§ Normalising Vehicle expressions

We assume the existence of a datatype  which faithfully represents the Vehicle language's <> type in Figure <ref> post-type checking.
Compilation operates on normalised expressions, where functions are fully applied and all user definitions have been inlined.
We therefore define the  datatype which represents the normal form of .
Normalisation takes place in an environment represented by the following datatype:
type Env = Stack Value

where the variable at Level $i$ is assigned the $i^{th}$  counting from the bottom of the stack.

Crucial to the efficiency of the algorithm, will be the ability to adjust the environment when we traverse over a universal or existential quantifier introducing a problem-space user variable.
Therefore our definition of  encodes expressions in weak-head normal form <cit.> as opposed to  full normal form.
Concretely this means that the body of lambdas and quantifiers are represented by an unnormalised (, ) pair, rather than a fully normalised .

data Value
= Lam Binder (Env, Expr)
| Var Level [Value]
| RealLiteral Real
| Add Value Value
| Mul Value Value
| BoolLiteral Bool
| Not Value
| And Value Value
| Or Value Value
| Equal Value Value
| NotEqual Value Value

| Order OrderOp Value Value
| If Value Value Value
| Forall Binder (Env, Expr)
| Exists Binder (Env, Expr)
| IndexLiteral Int
| VecLiteral [Value]
| At Value Value
| VectorEqual Value Value
| VectorNotEqual Value Value
| VectorAdd Value Value
| Map Value Value
| Fold Value Value Value
| NetworkApplication Name Value

Finally, for each of builtin functions in the  we define an evaluation function which tries to evaluate the expression if possible, otherwise returning the unevaluated expression in which case we say it's evaluation is “blocked”. For example, evalAdd would be defined as:
evalAdd : Value -> Value -> Value
evalAdd (RealLiteral x) (RealLiteral y) = RealLiteral (x + y)
evalAdd e1              e2              = Add e1 e2

We assume the existence of the eval function that uses these to take an environment  and an unnormalised  and returns a normalised :
eval : Env -> Expr -> Value.

§.§ Blocking

Compilation will proceed by normalising the expression and recursively processing “compilable” boolean values.

A  $\:e$ of type Bool is compilable if it is one of BoolLiteral, Not, And, Or, If, Equal, NotEqual, VectorEqual, VectorNotEqual, Order, Exists, Forall.

However, as illustrated by example7 in Section <ref>, we may encounter non-compilable higher-order functions whose evaluation is blocked.
Being able to unblock such expressions will be crucial for the algorithm to make progress.

A  $\:e$ is blocked by subexpressions $\{e_1, \ldots, e_n\}$ if replacing all $e_i$ with literals of the same type allows $e$ to evaluate further.

For example, `x + y' is blocked on the variables $\{$x$,$ y$\}$, and `fold (y -> x + y) 0 (map f xs)' is blocked on either $\{$xs$\}$ or $\{$map f xs$\}$.

If a blocked  $\:e$ is of type Vector or Index, there exists a set of blocking subexpressions, each one of which is either a i) a NetworkApplication, ii) an If, iii) a Var of type Vector or Index.

If $e$ is well-typed as a Vector then it must be one of the following cases:

* NetworkApplication f x - case i) holds.

* If c x y - case ii) holds.

* Var l - case iii) holds.

* Map f xs / Fold f e xs - must be blocked on xs of type Vector and so holds by induction.

* At xs i - must be blocked on at least one of xs of type Vector or i of type Index and so holds by induction.

* VectorAdd xs ys / VectorEqual xs ys / VectorNotEqual xs ys - must be blocked on at least one of xs or ys both of type Vector and so holds by induction.

Otherwise if $e$ is well-typed as a Index then it must one of the following cases:

* IndexLiteral - contradicts assumption that expression is blocked.

* If c x y - case ii) holds.

* Var l - case iii) holds.

* Fold f e xs - same as Fold case above.

* At xs i - same as At case above.

If a  $\:e$ is of type Bool and is not compilable, there exists a set of blocking subexpressions, each one of which is either i) a NetworkApplication, ii) an If or iii) a Var of type Vector, Index, or Bool.

As $e$ is well-typed and not one of the above constructors, then it must be one of:

* Var l - case iii) holds.

* Fold f e xs - must be blocked on xs of type Vector and so holds by Lemma <ref>.

* At xs i - must be blocked on at least one of xs of type Vector or i of type Index and so holds by Lemma <ref>.

Thanks to Lemma <ref>, it is possible to determine a way of unblocking non-compilable expressions of type Bool.
First, as a  is always of type Bool we know that the only variables in scope in a normalised expression are those introduced by Forall and Exists.
We can expand out Forall and Exists over the finite types Index or Bool into a finite conjunction/disjunction.
Therefore we can maintain the invariant that there are never any variables of Index or Bool in scope.
Consequently by Lemma <ref>, a non-compilable boolean expression must be blocked by a set of subexpressions that are either an i) If or a ii) NetworkApplication or iii) Var of type Vector.
We now explore how to handle each case.

§.§ Handling blocking Ifs

If we are blocked on an If, we can lift it out to the top-level. For example the value `haskell0 == (if c then x else y)' becomes `haskellif c then 0 == x else 0 == y'. This is achieved by the following function:  haskell
liftIf :: Monad m => Value -> (Value -> m Value) -> m Value
liftIf (If c x y) k = If c <$> liftIf e1 k <*> liftIf e2 k
liftIf e k = k e
\end{minted}
which takes a continuation and a value. If the value is an \haskell{If} it recursively lifts the true and false branches and reconstructs the \haskell{If} around the results, otherwise it applies the continuation.

\subsection{Handling blocking \haskell{NetworkApplication}s}
\label{sec:handling-blocked-networkapps}

If we are blocked on a \haskell{NetworkApplication}, we can replace it with the corresponding input/output embedding space variables as follows:
\begin{minted}{haskell}
elimNetworkApplication ::
(MonadState Ctx m, MonadWriter [Value] m) =>
Name -> Value -> m Value
elimNetworkApplication networkName arg = do
x <- addVectorVariableToCtx (numberOfInputs networkName)
y <- addVectorVariableToCtx (numberOfOutputs networkName)
tell [EqualVector arg x]
return y
\end{minted}
Note that we cannot conjunct the input equality to the output variable, as it will not type check and so instead the equality is added to a global list of such generated equalities using the \haskell{tell} function of \haskell{MonadWriter}.
For example, in the expression `\vehicle{0 <= f [a + 2] ! 0}' we would end up replacing `\vehicle{f [a + 2]}' with `\vehicle{y0}' and storing the equality `\vehicle{x0 == a + 2}'.

\subsection{Handling blocking \haskell{Var} of type \haskell{Vector}}
\label{sec:handling-vector-variable}

If we are blocked on a \haskell{Var} $x$ of type \vehicle{Vector A n} for some type $A$ and size $n$, we can replace it with the expression `\haskell{VecLiteral [x_1, ..., x_n]}' where $x_i$ are variables of type \vehicle{A} and represent the \emph{elements} of the vector.
The \haskell{vectorVars} field in the definition of \haskell{Ctx} in Section~\ref{sec:variables} is used to track the relationship between a vector variable and the variables representing its elements.

Note that similar to \vehicle{Index} and \vehicle{Bool} types, it would be possible to maintain the invariant that there are \emph{no} variables of type \vehicle{Vector} by expanding out a quantifier over vectors into a series of $n$ quantifiers over numeric variables.
Given this, an important question is why do we allow vector variables to appear in the expression at all?
The answer is efficiency. Consider the property `\vehicle{forall (a : Vector Real n) . f [a + c] >= 0}' where $c$ is some large constant vector.
If we immediately replace $a$ with $[a_1, ..., a_n]$ then later when solving for the user variables in terms of the embeddings space variables we will have to solve a system of $n$ identical equations:
\begin{align*}
a_1 + c_1 == X_1, \quad ..., \quad a_{n} + c_{n} == X_{n}
\end{align*}
and perform $n$ substitutions through the AST.
This would result in extremely slow compilation when $n$, the number of network inputs, is large.
However, if we choose to keep vector variables, then we can instead solve the following equation only once at the vector level:
\begin{align*}
a + c == X
\end{align*}
and perform only a single substitution through the tree.
However, this results in a subtly. In order to preserve the vector level equalities we need to not normalise the vector variables if they are arguments to operators in such an equality, i.e. \haskell{VectorAdd}, \haskell{VectorEqual}, \haskell{VectorNotEqual} and \haskell{NetworkApplication}.
In summary we therefore have the following code:
\begin{minted}{haskell}
elimVectorVariable :: MonadState Ctx m => Bool -> AnyVar -> m Value
elimVectorVariable isVecOpArg var = do
Context{..} <- get
case lookup var vectorVars of
Just elements | not isVecOpArg -> return (VecLiteral elements)
Nothing -> return (Var var [])
\end{minted}

\subsection{Unblocking}
\label{sec:unblocking}

Using these three functions we can define the function \vehicle{unblock} that locates the set of blocking subexpressions in a \valueType{} and unblocks them while keeping track of whether we are recursing over a vector operation or not.
\begin{minted}{haskell}
-- Invariant: `expr` is not a function type.
unblock ::  (MonadState Ctx m, MonadWriter [Value] m) =>
Bool -> Value -> m Value
unblock isVecOpArg expr = case expr of
-- Base cases --
Var v [] -> elimVectorVariable isVecOpArg v spine
BoolLiteral {} -> return expr
And{} -> return expr
Or{} -> return expr
Not{} -> return expr
If{} -> return expr
Forall{} -> return expr
Exists{} -> return expr
-- Inductive vector cases (VectorEqual, VectorNotEqual are similar)
NetworkApplication name x ->
liftIf (unblock True arg) (\x' -> elimNetworkApplication name x')
AddVector xs ys ->
liftIf (unblock True xs) (\xs' ->
liftIf (unblock True ys) (\ys' ->
AddVector xs' ys'))
-- Other inductive cases (all other missing cases are similar)
Add x y ->
liftIf (unblock False x) (\x' ->
liftIf (unblock False y) (\y' ->
evalAdd x' y'))
Fold f e xs ->
liftIf (unblock False xs) (\xs' -> evalFold f e xs')
\end{minted}
If we reach any non-trivial boolean structure or a variable we immediately stop recursing.
Otherwise if we reach a builtin function that does not have boolean arguments we first recurse into only the arguments that it could be blocked by, and then lift any \haskell{If}s out of those arguments.
For example, for \haskell{Add} we recurse into both the left and right arguments as it could be blocked on either.
For the higher order functions \haskell{Fold} we only recurse into the last vector argument.
Crucially this means that we never recurse into expressions which have a function type and therefore the \haskell{Lam} case is impossible.
The final case is the network application. Here we first recurse into the argument, and then replace the network application.
Note it is important to do it in this order, so that nested network applications as in \haskell{example1} in Section~\ref{sec:evaluation} are handled correctly.

Finally, we can define the function \haskell{tryUnblockBool} as follows:
\begin{minted}{haskell}
-- Invariant: `expr` is of type Bool.
tryUnblockBool :: MonadState Ctx m =>
Value -> (Value -> m a) -> m a -> m a
tryUnblockBool expr success failure = do
(unblockedExpr, networkEqualities) <- runWriterT (unblockExpr expr)
let finalExpr = foldr And unblockedExpr networkEqualities
return (if expr /= finalExpr then success finalExpr else failure)
\end{minted}
This function takes a \valueType{} of type \vehicle{Bool} and tries to unblock it, obtaining the list of noted network input equalities. Because we know the resulting expression is of type \vehicle{Bool} we can conjunct the resulting expression with any input variable equalities that were generated by replacing network applications while preserving well-typedness.
Crucially this also maintains equisatisfiability with the original expression as \haskell{unblock} never recurses inside any \haskell{Or}, \haskell{If} or \haskell{Exists} constructors.
Finally if we succeeded in unblocking the expression then we apply the continuation that was passed as the second argument to the unblocked result, otherwise we return the fallback value provided as the third argument.

\subsection{Queries}
\label{sec:queries}

We use a sparse representation of linear expressions that are generic over the variable and constant representations:
\begin{minted}{haskell}
data LinearExpr var const = LinearExpr (Map var Real) const
\end{minted}
and use it to define the following types of assertions where the linear expression is implicitly related to zero or the zero vector in the appropriate way:
\begin{minted}{haskell}
data RealEquality = RealEquality (LinearExpr RealVar Real)
data RealInequality = RealInequality Order (LinearExpr RealVar Real)
data VectorEquality = VectorEquality (LinearExpr VectorVar Vector)
data Assertion = RE RealEquality | RI RealInequality | VE VectorEquality
\end{minted}
The boolean structure will be represented by the following types:
\begin{minted}{haskell}
type Disjuncts a = NonEmptyList a
type Conjuncts a = NonEmptyList a
data MaybeTrivial a = Trivial Bool | NonTrivial a
data Tree a = Atom a | Conj (Tree a) (Tree a) | Disj (Tree a) (Tree a)
\end{minted}
Therefore the structure that is the final target for compilation:
\begin{minted}{haskell}
type Query = Tree Assertion
type QuerySet = Disjuncts Query
type Property = MaybeTrivial (Tree (Bool, QuerySet))
\end{minted}
A \haskell{Query} is a perhaps trivial tree of assertions. A \haskell{QuerySet} is a list of disjuncted queries each with a different solution for the problem-space variables in terms of embedding space variables (see \vehicle{example6} in Section~\ref{sec:challenging-examples}). A \haskell{Property} is a tree at the top level because the users may have arbitrary structure:
\begin{minted}{vehicle}
@property example8 = (forall a . f a > 0) or (exists b . f b < -5)
\end{minted}
At the leaves of the tree, in addition to a \haskell{QuerySet}, a \haskell{Bool} is also stored to mark whether the queries have been negated because there was a top-level universal quantifier.

\subsection{Compiling properties}
\label{sec:properties}

We now have sufficient machinery to start describing compilation.
The top-level boolean structure over a property compiled using \haskell{compileProperty}.
\begin{minted}{haskell}
-- Invariant: `expr` is of type Bool in the empty context.
compileProperty :: MonadState Ctx m => Expr -> m Property
compileProperty expr = go (eval emptyEnv expr)
go :: Value -> m Property
go expr = case expr of
-- Inductive cases --
And e1 e2 -> andTrivial Conj <$> go e1 <*> go e2
Or e1 e2 -> orTrivial Disj <$> go e1 <*> go e2
If c x y -> go (eliminateIf c x y)
Not e -> go =<< eliminateNot e
Equal e1 e2 -> tryUnblockBool expr go IMPOSSIBLE
NotEqual e1 e2 -> tryUnblockBool expr go IMPOSSIBLE
Order _ e1 e2 -> tryUnblockBool expr go IMPOSSIBLE
EqualsVector e1 e2 -> tryUnblockBool expr go IMPOSSIBLE
NotEqualsVector e1 e2 -> tryUnblockBool expr go IMPOSSIBLE
-- Base cases --
BoolLiteral b -> return (Trivial b)
Exists binder body -> do
querySet <- compileQueries (Exists binder body)
return (NonTrivial (Atom (False, querySet)))
Forall binder body -> do
querySet <- compileQueries (Exists binder (Not body))
return (NonTrivial (Atom (True, querySet)))
_ -> IMPOSSIBLE
\end{minted}

The first four inductive cases are relatively straight-forward. \haskell{And} and \haskell{Or} are preserved, and thanks to the assumption that \haskell{expr} is typable as \vehicle{Bool} in the empty context, we can eliminate \haskell{If} and \haskell{Not} easily (see Section~\ref{sec:elimination} for code).
The remaining inductive cases are when we hit an assertion. As we are in the empty context, all assertions must be trivial, and their evaluation cannot be blocked by a variable or a network application, and hence by Lemma~\ref{lem:blocked-bool} they must be blocked by \haskell{If} within the assertion.
Therefore we can call \haskell{tryUnblockBool} to lift the \haskell{If} out and then recurse, and failing to unblock is impossible.

There are three base cases.
\haskell{BoolLiteral} is the trivial case for when the property is vacuously true or false.
If we reach an \haskell{Exists} or \haskell{Forall} then we need to start compiling a set of queries. If we find a \haskell{Forall} we must first change it to an \haskell{Exists} and mark any queries as negated.
The function \haskell{compileQueries} is shown below.

\begin{minted}{haskell}
-- Invariant: `expr` is of type `Bool` in the current global context
compileQueries :: MonadState Ctx m =>
Value -> m (MaybeTrivial QuerySet)
compileQueries = go
go :: Value -> m (MaybeTrivial QuerySet)
go expr = case expr of
-- Base cases --
BoolLiteral b -> return (Trivial b)
Forall {} -> error "alternating quantifiers"
-- Maybe base cases --
Order op x y -> tryUnblockBool expr go (compileIneq op x y)
Equal x y -> tryUnblockBool expr go (compileEq x y)
VectorEqual xs ys -> tryUnblockBool expr go (compileVecEq xs ys)
-- Inductive cases --
If c x y -> go (eliminateIf c x y)
Not x -> go =<< eliminateNot x
NotEqual x y -> go =<< eliminateNotEqual x y
VectorNotEqual xs ys -> go =<< eliminateVectorNotEqual xs ys
Exists binder body -> eliminateExists binder body
And x y -> andTrivial (cartesianProduct Conj) <$> go x <*> go y
Or x y -> orTrivial (++) <$> go x <*> go y
_ -> tryUnblockBool expr go IMPOSSIBLE
\end{minted}

The bases cases are when we reach either a boolean literal or a positive assertion.
However, the latter three may only appear to be a base case if the relation has an \vehicle{If} hidden inside of it as demonstrated by \vehicle{example5} in Section~\ref{sec:challenging-examples}.
We therefore need to first recurse through the expression and lift all the \vehicle{If} out of the relation using \vehicle{tryUnblockBool}.
Only if that fails do we know that we can then proceed to compile the assertion.
Compiling assertions is relatively unexciting. The full code for \haskell{compileEq} and \haskell{compileIneq} can be found in the appendices. The only interesting function is \haskell{compileVecEq}, which may fail to compile.

The recursive cases can be split into three classes. The first class consists of \haskell{If}, \haskell{Not}, \haskell{NotEqual}, \haskell{VectorNotEqual}, \haskell{Exists}, all of which can be eliminated. The second class consists of \haskell{Or} and \haskell{And}, and in this case we recurse directly and then combine the results by appending them or taking the cartesian product respectively.
In the final case, we know we have a blocked, non-compilable boolean value and therefore by Lemma~\ref{lem:blocked-bool} we can unblock and then recurse.

\subsection{Elimination operations}
\label{sec:elimination}

The first three elimination algorithms for \haskell{If}, \haskell{NotEqual} and \haskell{NotVectorEqual} are direct translations as they can be expressed natively in the DSL:
\begin{minted}{haskell}
eliminateNotEqual :: Value -> Value -> Value
eliminateNotEqual x y = Or (Order Le x y) (Order Ge x y)

eliminateVectorNotEqual :: Value -> Value -> Value
eliminateVectorNotEqual x y = -- Haskell encoding of
-- `fold (\a b -> a or b) False (zipWith (\a b -> a != b) x y)`

eliminateIf :: Value -> Value -> Value -> Value
eliminateIf c x y = Or (And c x) (And (Not c) y)
\end{minted}

The elimination algorithm for \haskell{Not} is a little more complicated and may once again require the argument be unblocked, invoking Lemma~\ref{lem:blocked-bool} for soundness.
\begin{minted}{haskell}
eliminateNot :: MonadState Ctx m => Value -> m Value
eliminateNot expr = case expr of
Not x -> return x
BoolLiteral b -> return (BoolLiteral (not b))
Order op x y -> return (Order (neg op) x y)
Equal x y -> return (NotEqual x y)
NotEqual x y -> return (Equal x y)
Forall binder (env, body) -> return (Exists binder (env, Not body))
Exists binder (env, body) -> return (Forall binder (env, Not body))
EqualsVector xs ys -> NotEqualsVector xs ys
NotEqualsVector xs ys -> EqualsVector xs ys
Or x y -> And <$> eliminateNot x <*> eliminateNot y
And x y -> Or <$> eliminateNot x <*> eliminateNot y
If c x y -> If c <$> eliminateNot x <*> eliminateNot y
_ -> tryUnblockBool expr eliminateNot IMPOSSIBLE
\end{minted}
The most complicated operation is the elimination of an existential introducing a new problem-space variable.
\begin{minted}{haskell}
eliminateExists :: MonadState Ctx m =>
Binder -> (Env, Expr) -> m (MaybeTrivial QuerySet)
eliminateExists binder env body = do
v <- addVectorVariableToCtx (vectorSize binder)
normExpr <- eval (push v env) body
queries <- compileQueries normExpr
return (solveVariable queries v)
\end{minted}
First we create the problem space variable from the binder assuming, without loss of generality, that the variable is a vector.
Next we add that variable to the environment and normalise the body of the quantifier, which we then recursively compile to obtain a set of queries that are using that variable.
The final step is to then solve for the problem-space variable in terms of the embedding space variables in each of the individual queries.
\begin{minted}{haskell}
solveVariable :: MaybeTrivial Queries -> UserVar -> MaybeTrivial Queries
solveVariable queries v = do
let solve query = fmap (solveConstraints v) (findConstraints v query)
let solvedQueries = fmap (fmap solve) queries
flattenTrivial (fmap (filterTrivial . fmap filterTrivial) solvedQueries
\end{minted}
Solving for the user variable within a single query comes in two stages. First we use \haskell{findConstraints} to extract the relevant assertions that involve the variable, and then we use \haskell{solveConstraints} to eliminate it.
The last line flattens the result of type \haskell{MaybeTrivial (Disjuncts (Disjuncts (MaybeTrivial QuerySet))} back to \haskell{MaybeTrivial QuerySet}.

Before we define \haskell{findConstraints} a key thing to note is that, as demonstrated by \vehicle{example6} in Section~\ref{sec:challenging-examples}, a \haskell{Query} may contain disjunctions and therefore each side of the disjunctions may contain a different solution for the problem-space variable.
Consequently when this happens it is necessary to split the query into multiple separate queries, with each query having consistent constraints for the variable.
Furthermore, normally we would expect that every problem space variable is used as a unique network input, and therefore will generate an equality that we can solve.
However, it is possible, as demonstrated by the LHS of \vehicle{example6}, that the system of linear equations relating the problem space and embedding space variables may be underconstrained.
Therefore we may need to fall back on Fourier-Motzkin elimination~\cite{} which uses inequalities instead of equalities to eliminate a problem-space variable.
Given this we define the notion of a \haskell{ConstrainedQuery} as follows:
\begin{minted}{haskell}
data ConstrainedQuery
= Equality (Either RealEquality VectorEquality, MaybeTrivial Query)
| Inequalities (ConjunctAll RealInequality, MaybeTrivial Query)
| NoConstraints Query
\end{minted}
In particular, it represents the implicit conjunction of the strongest set of constraints over the problem space variable found so far and the remainder of the query.
We can then define \vehicle{findConstraints} as follows:
\begin{minted}{haskell}
findConstraints :: UserVar -> Query -> Disjuncts ConstrainedQuery
findConstraints var query = case query of
Query assertion -> if var `notReferencedBy` assertion
then (NoConstraints query)
else case assertion of
(RE realEq) -> Equality (Left realEq, Trivial True)
(RI realIneq) -> Inequalities ([realIneq], Trivial True)
(VE vectorEq) -> Equality (Right vectorEq, Trivial True)
Disj x y -> findConstraints x <> findConstraints y
Conj x y -> do
r1 <- findConstraints x
r2 <- findConstraints y
let (eqsLHS, remainingLHS) = partition (map (shortCircuitEq y) r1)
let (eqsRHS, remainingRHS) = partition (map (shortCircuitEq x) r2)
let ineqs = cartesianProduct merge remainingLHS remainingRHS

if null remainingLHS then eqsLHS
else if null remainingRHS then eqsRHS
else eqsLHS <> eqsRHS <> ineqs
\end{minted}
The \haskell{Conj} case is so complicated as we try to avoid unnecessarily converting the constrained to DNF, as that can lead to exponential explosion in the number of queries.
This achieved by the function \haskell{shortCircuit}, which exploits the fact that if we find an equality constraint in one side of the conjunction, we can ignore any disjuncts found on the other side.
\begin{minted}{haskell}
shortCircuit ::
Query ->
ConstrainedQuery ->
Either ConstrainedQuery ConstrainedQuery
shortCircuit disjunct constraint = case constraint of
Equality (eq, remaining) -> do
let newRemaining = andTrivial Conj remaining (NonTrivial disjunct)
Left (Equality (eq, newRemaining))
_ -> Right constraint
\end{minted}
If there are still disjuncts left on both sides of the tree, then we do have to do a full DNF conversion using \haskell{cartesianProduct} and \haskell{merge}, the latter of which is:
\begin{minted}{haskell}
merge :: ConstrainedQuery -> ConstrainedQuery -> ConstrainedQuery
merge c1 c2 = case (c1, c2) of
(NoConstraints t1, NoConstraints t2) ->
NoConstraints (Conj t1 t2)
(NoConstraints t1, Inequalities (ineqs2, t2)) ->
Inequalities (ineqs2, andTrivial Conj (NonTrivial t1) t2)
(Inequalities (ineqs1, t1), NoConstraints t2) ->
Inequalities (ineqs1, andTrivial Conj t1 (NonTrivial t2))
(Inequalities (ineqs1, t1), Inequalities (ineqs2, t2)) ->
Inequalities (ineqs1 <> ineqs2, andTrivial Conj t1 t2)
\end{minted}
Once we have a \haskell{ConstrainedQuery} we can solve for the constraints as follows:
\begin{minted}{haskell}
solveConstraints :: UserVar -> ConstrainedQuery -> MaybeTrivial QuerySet
solveConstraints var query = case query of
Equality (eq, remainingQuery) ->
fmap singleton (substituteFor var eq remainingQuery)
Inequalities (ineqs, remainingQuery) -> do
let newIneqs = fourierMotzkinElimination var ineqs
let newTree = andTrivial Conj remainingQuery newIneqs
fmap singleton newTree
NoConstraints query -> do
Context{..} <- get
case lookup var vectorVars of
Nothing -> query
Just elemVars -> foldl solveVariable (NonTrivial [tree]) elemVars
\end{minted}
If we found an equality involving the problem-space variable, we can use the equality to eliminate it from the remainder of the query.
Note that in the case where \haskell{var} is a vector variable, we also need to subsitute for its element variables as well.
If we found a set of inequalities involving the variable, we can use Fourier-Motzkin elimination to generate a new equi-satisfiable set of inequalities without that variable.
In this case, we don't have to substitute through, because the remainder of the tree is guaranteed to not refer to the user variable.
If we found no constraints, then we check whether the variable being solver for is a vector variable.
If not then the variable was never used so we can just return the query.
If it is a vector variable, then we are forced instead to solve for each of the variables representing its elements individually.

\subsection{Post-processing}

We have now completely described the function \haskell{compileProperty} that returns a set of queries referring only to the network input and output variables in the embedding space.
These before being fed into a neural network solver they require a few uninteresting post-processing steps.
Firstly, any \haskell{VectorEquality}s must be trivially converted to a conjunction of \haskell{RealEquality}s.

Secondly, note that in \haskell{elimNetworkApplication} in Section~\ref{sec:handling-blocked-networkapps}, we assigned a fresh set of embedding space variables to \emph{every} network application despite it being the case that the user might have referred to the same neural network application more than once in a specification.
In this case the solver will likely unnecessarily duplicate the entire network internally, resulting in exponentially greater amount of work.
Eliminating these duplicates can't be handled via a cache inside \haskell{elimNetworkApplication}, as the notion of two applications being the same is a semantic rather than a syntactic notion, e.g. in
\begin{minted}{vehicle}
@property example9 = forall (x y : Vector Real 1) . f x == f y and x == y
\end{minted}
the two applications \vehicle{f x} and \vehicle{f y} are syntactically different, but semantically identical.
Therefore this procedure must be performed after we know solutions for the user variables.

%To evaluate the performance of the algorithm, we evaluate it on three common example specifications.
%\item Epsilon-ball robustness - ?
%\item Monotonicity - ?
%\item Hierarchical classification - ?
%	\input{figures/verifier-compilation-scaling.tex}
%	\caption{Performance when compiling the specifications as the number of inputs to the network grows.}
%	\label{fig:scaling-experiments}
%\matthew{Add other two experiments}
%The results are shown in Figure~\ref{fig:scaling-experiments} and show that compilation time grows linearly for the first examples and quadratically for the second examples.
%These results are asymptotically optimal, as for the first two experiments the size of the queries themselves grow in proportion to $O(n)$ and $O(n^2)$ for the latter experiment.

\section{Conclusions and open questions}
\label{sec:conclusion}

In this paper we have presented an algorithm for compiling down high-level problem-space specifications to embedding space queries for neural network solvers.
Such specifications are considerably more expressive than other state of the art of the DSLs for similar tools.
% Furthermore, we have shown how when implemented in the Vehicle framework it scales optimally for a range of common specifications.

There remain important open questions in the domain of neural network DSL specification design, of which we now discuss two.
Firstly, there is the issue that the soundness of compilation depends on the use of real numbers.
Most DSLs, including this one, assume real number arithmetic, but are unsound under floating point arithmetic.
In this algorithm in particular, the soundness of \haskell{compileRealAssert}, \haskell{eliminateNotEqual} and \haskell{solveConstraints} all depend on the use of real arithmetic.

A second open problem is the ergonomic compilation to abstract-interpretation based solvers. In particular, when the embedding space queries don't conform to the supported abstract domains, how to provide meaningful explanations to the user in terms of the problem space specification.

\bibliographystyle{splncs04}
\begin{thebibliography}{10}
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\urlprefix}{URL }
\providecommand{\doi}[1]{https://doi.org/#1}

[1]
{VNNLib format}, \url{https://vnnlib.org/}, accessed on 01.12.2021

[2]
Daggitt, M.L., Atkey, R., Kokke, W., Komendantskaya, E., Arnaboldi, L.:
Compiling higher-order specifications to {SMT} solvers: How to deal with
rejection constructively. In: Krebbers, R., Traytel, D., Pientka, B.,
Zdancewic, S. (eds.) Proceedings of the 12th {ACM} {SIGPLAN} International
Conference on Certified Programs and Proofs, {CPP} 2023, Boston, MA, USA,
January 16-17, 2023. pp. 102--120. {ACM} (2023).
\doi{10.1145/3573105.3575674}, \url{https://doi.org/10.1145/3573105.3575674}

[3]
Daggitt, M.L., Kokke, W., Atkey, B., Natalia, S., Arnaboldi, L.,
Komendantskaya, E.: {Vehicle: Bridging the Embedding Gap in the Verification
of Neuro-Symbolic Programs}  \textbf{abs/2007.11206} (2024), \url{???}

[4]
Filli{\^a}tre, J.C., Paskevich, A.: {Why3—where programs meet provers}. In:
Programming Languages and Systems: 22nd European Symposium on Programming,
ESOP 2013, Held as Part of the European Joint Conferences on Theory and
Practice of Software, ETAPS 2013, Rome, Italy, March 16-24, 2013. Proceedings
22. pp. 125--128. Springer (2013)

[5]
Girard-Satabin, J., Alberti, M., Bobot, F., Chihani, Z., Lemesle, A.: Caisar: A
platform for characterizing artificial intelligence safety and robustness.
In: AISafety. CEUR-Workshop Proceedings, Vienne, Austria (Jul 2022),
\url{https://hal.archives-ouvertes.fr/hal-03687211}

[6]
Ivanov, R., Weimer, J., Alur, R., Pappas, G.J., Lee, I.: Verisig: verifying
safety properties of hybrid systems with neural network controllers. In:
Proceedings of the 22nd ACM International Conference on Hybrid Systems:
Computation and Control. pp. 169--178 (2019)

[7]
Katz, G., Barrett, C., Dill, D.L., Julian, K., Kochenderfer, M.J.: Reluplex: An
efficient smt solver for verifying deep neural networks. In: International
Conference on Computer Aided Verification. pp. 97--117. Springer (2017)

[8]
Katz, G., Huang, D.A., Ibeling, D., Julian, K., Lazarus, C., Lim, R., Shah, P.,
Thakoor, S., Wu, H., Zelji{\'c}, A., et~al.: {The Marabou framework for
verification and analysis of deep neural networks}. In: International
Conference on Computer Aided Verification. pp. 443--452. Springer (2019)

[9]
Liu, X., Han, X., Zhang, N., Liu, Q.: Certified monotonic neural networks.
Advances in Neural Information Processing Systems  \textbf{33},  15427--15438

[10]
Marlow, S., et~al.: Haskell 2010 language report  (2010)

[11]
M{\"u}ller, M.N., Makarchuk, G., Singh, G., P{\"u}schel, M., Vechev, M.: Prima:
general and precise neural network certification via scalable convex hull
approximations. Proceedings of the ACM on Programming Languages
\textbf{6}(POPL),  1--33 (2022)

[12]
Park, D.Y., Cha, M.H., Kim, D., Han, B., et~al.: Learning student-friendly
teacher networks for knowledge distillation. Advances in neural information
processing systems  \textbf{34},  13292--13303 (2021)

[13]
Shriver, D., Elbaum, S., Dwyer, M.B.: {DNNV: A Framework for Deep Neural
Network Verification}. In: Silva, A., Leino, K.R.M. (eds.) Computer Aided
Verification. pp. 137--150. Springer International Publishing, Cham (2021)

[14]
Wang, S., Zhang, H., Xu, K., Lin, X., Jana, S., Hsieh, C.J., Kolter, J.Z.:
{Beta-CROWN}: Efficient bound propagation with per-neuron split constraints
for complete and incomplete neural network verification. Advances in Neural
Information Processing Systems  \textbf{34} (2021)

\end{thebibliography}

%\section{Missing code}
%compileVecEq :: _
%compileVecEq = _
%compileRealAssert :: (MonadState Ctx m) => Value -> m RealLinearExpr
%compileRealAssert = go
%  where
%  go :: Value -> m RationalLinearExpr
%  go e = case e of
%    RatLiteral l -> return (constantExpr l)
%    VBoundVar lv [] -> singletonVarExpr 0 <$> getRealVariable lv
%    Neg v -> scaleExpr (-1) <$> go v
%    Add e1 e2 -> addExprs 1 1 <$> go e1 <*> go e2
%    Mul e1 e2 -> do
%      e1' <- go e1
%      e2' <- go e2
%      case (isConstant e1', isConstant e2') of
%        (Just c1, _) -> return (scaleExpr c1 e2')
%        (_, Just c2) -> return (scaleExpr c2 e1')
%        _ -> error "non-linear relationship"

\end{document}