probability greater or equal to $1-4n^{-2}-\exp(-n\mathbb{E}[I_{\Omega}])$,

$\displaystyle[4\beta^{-1}L\psi]_{i}$
$\displaystyle=[4\beta^{-1}Lu]_{i}-[4\beta^{-1}Lv]_{i}-\epsilon(C_{\max}+1)[4\beta^{-1}L\phi]_{i}$
$\displaystyle\leq f_{i}+\epsilon
C_{u,i}-f_{i}-\epsilon(C_{\max}+1)(1-C_{\phi,i}\epsilon)$
$\displaystyle\leq-\epsilon\left(1-C_{\phi,i}\epsilon(1+C_{\max})\right)\leq
0.$

We did not multiply the exponent by 2 in the probability estimate
$1-4n^{-2}-\exp(-n\mathbb{E}[I_{\Omega}])$ because the same set of points is
used in the estimates for $u$ and $\phi$.

Hence, by the discrete maximum principle in Lemma 4.6, the minimum of $\psi$
is achieved on $\mathcal{D}$. Since $u_{i}=v_{i}=g_{i}$ for all
$i\in\mathcal{D}$ and $\phi_{i}=0$ for $i\in\mathcal{D}$, it follows that
$\psi_{i}=0$ for all $i\in\mathcal{D}$. Therefore, $\psi_{i}\geq 0$ for all
$1\leq i\leq n$, i.e.,

$u_{i}-v_{i}\geq\epsilon(C_{\max}+1)\phi_{i}\equiv-\epsilon(C_{\max}+1)|\phi_{i}|\quad\forall
1\leq i\leq n.$ (224)

Next, we define the second comparison function

$\psi=u-v+\epsilon\left(\max_{1\leq i\leq n}C_{u,i}+1\right)\phi\equiv
u-v+\epsilon(C_{\max}+1)\phi$ (225)

and argue that $[L\psi]_{i}\geq 0$ for all $i\in\mathcal{I}$ with high
probability. Indeed, (220) and (222) imply that with probability greater or
equal to $1-4n^{-2}-\exp(-n\mathbb{E}[I_{\Omega}])$,

$\displaystyle[4\beta^{-1}L\psi]_{i}$
$\displaystyle=[4\beta^{-1}Lu]_{i}-[4\beta^{-1}Lv]_{i}+\epsilon(C_{\max}+1)[4\beta^{-1}L\phi]_{i}$
$\displaystyle\geq f_{i}-\epsilon
C_{u,i}-f_{i}+\epsilon(C_{\max}+1)(1-C_{\phi,i}\epsilon)$
$\displaystyle\geq\epsilon\left(1-C_{\phi,i}\epsilon(1+C_{\max})\right)\geq 0$

for $\epsilon$ small enough. By Lemma 4.6, the maximum of $\psi_{i}$ is
achieved on $\mathcal{D}$. Since $u_{i}=v_{i}=g_{i}$ and $\phi_{i}=0$ for
$i\in\mathcal{D}$, $\psi_{i}\leq 0$ for all $i\in\mathcal{D}$. Hence
$\phi_{i}\leq 0$ for all $1\leq i\leq n$, i.e.,

$u_{i}-v_{i}\leq-\epsilon(C_{\max}+1)\phi_{i}\equiv\epsilon(C_{\max}+1)|\phi_{i}|\quad\forall
1\leq i\leq n.$ (226)

Inequalities (224) and (226) imply that, for $\epsilon$ small enough, with
probability greater or equal to $1-4n^{-2}-\exp(-n\mathbb{E}[I_{\Omega}])$,

$|u(x_{i})-v(x_{i})|\leq\epsilon(C_{\max}+1)|\phi(x_{i})|\quad\forall 1\leq
i\leq n$ (227)

where $\phi$ is the solution to (221). Finally, we note that

$C_{\max}=\max_{1\leq i\leq
n}C_{u,i}\leq\max_{x\in\mathcal{M}}C_{u,\mu,\rho}(x)$

where $C_{u,\mu,\rho}(x)$ is defined be (104) in the statement of Theorem 3.8.
∎

##  5 Test problems

###  5.1 Illustrating error model

In Corollary 3.2.1 it was shown that bias error prefactors $\mathcal{B}_{2}$
and $\mathcal{B}_{3}$ can be eliminated when $\rho$ is uniform on
$\mathcal{M}$ and when the test function $f$ is the committor. One can expect
that these cancellations reduce the magnitude of the overall bias error
prefactor $\mathcal{B}_{1}+\mathcal{B}_{2}+\mathcal{B}_{3}$, though this might
not be necessary. The example presented in this section demonstrates that
these cancellations indeed reduce the magnitude of the prefactor, thus showing
that using uniform densities and approximating the committor function is a
viable strategy for improving the accuracy of TMDmap.

####  5.1.1 Theoretical considerations

The leading order terms in the error model in Theorem 1.1 can be briefly
stated as follows:

$\displaystyle|4\beta^{-1}L^{(n)}_{\epsilon,\mu}f(x)-\mathcal{L}f(x)|\leq\sqrt{\frac{\log
n}{n\epsilon^{2+d/2}}}|\mathcal{V}[f,\mu,\rho](x)|+\epsilon|\mathcal{B}[f,\mu,\rho](x)|.$
(228)

Here $\mathcal{V}$ and $\mathcal{B}$ are bias and variance prefactors
respectively given by

$\displaystyle\mathcal{V}$
$\displaystyle:=\frac{1}{\rho(x)^{1/2}}\left(2\|\nabla
f(x)\|\epsilon^{1/2}+11|f(x)|\right),$ (229) $\displaystyle\mathcal{B}$
$\displaystyle:=\mathcal{B}_{1}+\mathcal{B}_{2}+\mathcal{B}_{3}.$ (230)

$\mathcal{V}$ and $\mathcal{B}$ are both $O(1)$ in $n$ and $\epsilon$. While
(228) is an upper bound on the absolute error
$|4\beta^{-1}L^{(n)}_{\epsilon,\mu}f(x)-\mathcal{L}f(x)|$ due to the presence
of the _variance error_ , the _bias error_ estimate from Theorem 3.1 is
asymptotically sharp and does not involve absolute values. Thus if $\epsilon$
and $n$ are scaled such that the variance error is fixed, then with high
probability the error behaves approximately as a linear function in
$\epsilon$, i.e.

$\displaystyle
4\beta^{-1}L^{(n(\epsilon))}_{\epsilon,\mu}f(x)-\mathcal{L}f(x)\sim a+\epsilon
b.$ (231)

To fix the variance error, the number of data points $n$ must depend on
$\epsilon$ such that

$\displaystyle\frac{\log n(\epsilon)}{n(\epsilon)\epsilon^{2+d/2}}\sim O(1).$
(232)

After fixing the variance error, the slope of the linear function (231) yields
an estimate for the prefactor of the bias error. Using this method, we will
obtain a prefactor estimate for a 1D test system embedded in 2D.

####  5.1.2 1D test system

Consider a system on the unit circle $\mathbb{S}^{1}$ parameterized by arc
length (Figure 3). Let it evolve according to the overdamped Langevin dynamics
(1) with

$\displaystyle
V(\theta)=\left(4\cos^{2}\left(\frac{\theta}{2}\right)-\frac{3}{2}\right)^{2}.$
(233)

Figure 3: The potential $V(\theta)$, $0\leq\theta<2\pi$, given by (233) and
the invariant measure $\mu$ for a 1D system on $\mathbb{S}^{1}$ following the
overdamped Langevin dynamics (1) with the drift $\nabla V$.

In this case, the committor can be derived analytically. Let
$0<\theta_{1}<\theta_{2}<2\pi$ be the minima of $V$ where

$\displaystyle\theta_{1}=2\cos^{-1}\left(-\frac{\sqrt{3}}{2\sqrt{2}}\right)\approx
1.8235,\quad\theta_{2}=2\cos^{-1}\left(\frac{\sqrt{3}}{2\sqrt{2}}\right)\approx
4.4597.$ (234)

Set $A=[\theta_{1}-r,\theta_{1}+r]$ and $B=[\theta_{2}-r,\theta_{2}+r]$ where
$r=0.1$. Moreover, fix $\beta=1$. Then the committor is

$\displaystyle q(\theta)$
$\displaystyle=\frac{\int_{\theta_{1}+r}^{\theta}\exp{\left(\beta
V(\theta^{\prime})\right)}\,d\theta^{\prime}}{\int_{\theta_{1}+r}^{\theta_{2}-r}\exp{\left(\beta
V(\theta^{\prime})\right)}\,d\theta^{\prime}},\quad$
$\displaystyle\theta\in[\theta_{1}+r,\theta_{2}-r],$ $\displaystyle q(\theta)$
$\displaystyle=\frac{\int_{\theta}^{2\pi}\exp{\left(\beta
V(\theta^{\prime})\right)}\,d\theta^{\prime}+\int_{0}^{\theta_{1}-r}\exp{\left(\beta
V(\theta^{\prime})\right)}\,d\theta^{\prime}}{\int_{\theta_{2}+r}^{2\pi}\exp{\left(\beta
V(\theta^{\prime})\right)}\,d\theta^{\prime}+\int_{0}^{\theta_{1}-r}\exp{\left(\beta
V(\theta^{\prime})\right)}\,d\theta^{\prime}},\quad$
$\displaystyle\theta\in[\theta_{2}+r,2\pi),$ (235) $\displaystyle q(\theta)$
$\displaystyle=\frac{\int_{\theta}^{\theta_{1}-r}\exp{\left(\beta
V(\theta^{\prime})\right)}\,d\theta^{\prime}}{\int_{\theta_{2}+r}^{2\pi}\exp{\left(\beta
V(\theta^{\prime})\right)}\,d\theta^{\prime}+\int_{0}^{\theta_{1}-r}\exp{\left(\beta
V(\theta^{\prime})\right)}\,d\theta^{\prime}},\quad$
$\displaystyle\theta\in[0,\theta_{1}-r].$

####  5.1.3 The experimental setup

To demonstrate the error model in Theorem 1.1 and an improvement in the bias
error prefactor (10), the following experiment is conducted.

  1. 1.

For each $\epsilon$, the number $n(\epsilon)$ is chosen such that the error
behaves according to (231):

$\displaystyle\frac{n(\epsilon)}{\log(n(\epsilon))}=0.25\epsilon^{-5/2}.$
(236)

The constant prefactor of $0.25$ is chosen as the $O(1)$ term in (232) to make
the resulting number of points $n(\epsilon)$ fall into the interval
$[10^{4},3\cdot 10^{4}]$ for $\epsilon\in[0.023,0.033]$.

  2. 2.

To measure the effect of changing the sampling density, for each choice of
$\epsilon$, $n(\epsilon)$ points are drawn i.i.d. on $\mathbb{S}^{1}$ using
two sampling densities: the uniform density $\rho^{u}$ and a non-uniform
density $\rho^{n.u}$ specified in (238). The resulting datasets are defined as
follows:

$\displaystyle\mathcal{X}^{u}$
$\displaystyle=\psi\left(\\{\theta^{u}_{i}\\}_{i=1}^{n}\right),\quad$
$\displaystyle\theta^{u}_{i}$
$\displaystyle\sim\rho^{u}:=\mathcal{U}[0,2\pi],$ (237)
$\displaystyle\mathcal{X}^{n.u}$
$\displaystyle=\psi\left(\\{\theta^{n.u}_{i}\\}_{i=1}^{n}\right),\quad$
$\displaystyle\theta^{n.u}_{i}$
$\displaystyle\sim\rho^{n.u}:=\pi+0.1+\sigma\pi\\{\mathcal{N}(0,1)\\}.$ (238)

The function $\psi(\theta)=[\cos\theta,\sin\theta]^{\top}$ is the embedding
map $\psi:\mathbb{S}^{1}\rightarrow\mathbb{R}^{2}$ and
$\\{\mathcal{N}(0,1)\\}$ denotes the fractional part of the standard normal
distribution. The fractional normal distribution is shifted such that the
resulting distribution $\rho^{n.u}$ oversamples around the point
$\theta=\pi+0.1$ lying in the transition region and undersamples around the
sets $A$ and $B$. The parameter $\sigma$ modulates the variance of the
fractional normal distribution. It was set to $\sigma=0.2$. The empirical
densities of the data are visualized in Figure 4.

Figure 4: The empirical densities of $\mathcal{X}^{u}$ and
$\mathcal{X}^{n.u}$. Note that the nonuniform points $\mathcal{X}^{n.u}$ are
clustered in the region around $\theta=\pi+0.1$.

  3. 3.

To measure the effect of changing the test function,
$L^{(n(\epsilon))}_{\epsilon,\mu}f(x)$ was computed for two choices of test
function: the committor $q$ given by formula (235) and
$t(\theta)=\sin(\theta)$ (Figure 5). The evaluation point was fixed to be
$\theta=\pi$.

Figure 5:  The two choices of $f$ for evaluating the consistency error as in
(231): the committor $q$ (left) and the trigonometric polynomial
$t(\theta)=\sin\theta$ (right).

  4. 4.

For a given combination of sampling density and test function, the error in
the left-hand side of (231) was computed. This error is a function of
$\epsilon$ only.

The above experiment was conducted for 10 equispaced values of $\epsilon$ in
the range $[0.023,0.033]$. Moreover, since
$4\beta^{-1}L^{(n(\epsilon))}_{\epsilon,\mu}f(x)-\mathcal{L}f(x)$ is random
due to the i.i.d. draws from $\rho$, the experiment was repeated 50 times for
every value of $\epsilon$ to mitigate the statistical variation. Then a linear
fit to (231) was obtained to this data as a function of $\epsilon$ and the
slope of this fit was taken as the estimate of the prefactor of the bias error
for the given combination of $f$ and $\rho$. A total of four linear fits were
produced. The results are presented in Figure 6.

####  5.1.4 Summary of numerical results for reduction in bias error

We computed the error
$4\beta^{-1}L^{n(\epsilon)}_{\epsilon,\mu}f(x)-\mathcal{L}f(x)$ for $f=q,t$
and $\rho=\rho^{u},\rho^{n.u}$ as a function of $\epsilon$ ranging in
$[0.023,0.033]$. In Figure 6, a linear fit is obtained for
$4\beta^{-1}L^{n(\epsilon)}_{\epsilon,\mu}f(x)-\mathcal{L}f(x)$ for each
choice of $f$ and $\rho$, where the slope of the linear regression is then an
approximation of the bias error prefactor involved. Note that the bias error
prefactor in general may be negative since the bias error formula in Theorem
3.1 does not require taking absolute values. Here, it is of interest to
quantify the prefactor of lowest _magnitude_ to measure which combination of
sampling density and test function gives faster convergence of the bias error.
The estimates of the prefactor magnitude are presented in Table 2. We find
that the bias error prefactor of the lowest magnitude occurs when using a
uniform density as input to TMDmap and applying the resulting generator to the
committor function, presumably due to the uniform sampling and the committor
function canceling the bias prefactors $\mathcal{B}_{2}$ and $\mathcal{B}_{3}$
at the given point $x$.

| Sampling density, Test function
---|---
|  $\rho^{n.u}$,$t(\theta)$. |  $\rho^{u}$, $t(\theta)$ |  $\rho^{n.u}$, $q$ |  $\rho^{u}$, $q$
Abs. bias error prefactor | 1.024 | 0.778 | 1.148 | 0.398
Table 2: The bias error prefactor $b$ in (231) is estimated for two choices of
$\rho$ and two choices of $f$ for the 1D system described in Section 5.1.2.
The overall prefactor of lowest magnitude is obtained when using a uniform
sampling density and applying the TMDmap generator to the committor function.
This fact can be explained by the cancellations of the prefactors
$\mathcal{B}_{2}$ and $\mathcal{B}_{3}$ in Theorem 3.1. Figure 6: The
consistency error for the committor, $f=q(\theta)$, and for $f=\sin(\theta)$
at $\theta=\pi$ as a function of $\epsilon$. The number of points
$n(\epsilon)$ for each $\epsilon$ is chosen according to (236) to keep the
variance error constant. The datasets we sampled according to the uniform and
the nonuniform densities in (237) and (238) respectively.

###  5.2 Calculating the committor function

The goal of this section is to compute the committor as the solution to (7)
using TMDmap. The main focus will be on studying the discrepancy between the
numerical solution $q^{{\sf TMD}}$ and the true solution $q$ as a function of
the sampling density $\rho$ and kernel bandwidth $\epsilon$.

####  5.2.1 Theoretical considerations

Corollary 3.8.1 suggests that solving the committor problem numerically with a
uniform sampling density $\rho$ will likely produce improvements in the
accuracy of the numerical solution. Furthermore, bandwidth $\epsilon$ needs to
be commensurate with the fixed number of points $n$ to ensure that the error
formula in Theorem 1.1 holds.

####  5.2.2 Optimal bandwidth

Given fixed $n,\rho,\mu,$ and $x$, the expression in (228) may be minimized
for $\epsilon$ to obtain an optimal bandwidth $\epsilon^{*}$ given by

$\displaystyle\epsilon^{*}\sim(1+d/4)^{\frac{1}{2+d/4}}\left(\frac{\log
n}{n}\right)^{\frac{1}{4+d/2}}\left|\frac{V}{\mathcal{B}_{1}+\mathcal{B}_{2}+\mathcal{B}_{3}}\right|^{\frac{1}{2+d/4}}.$
(239)

A similar version of the optimal bandwidth formula was derived in [19] for
Laplacian eigenmaps. That work, however, stopped short of expanding further on
this formula due to the complicated dependence of the optimal bandwidth on
manifold-related parameters. In contrast, this paper offers a full enumeration
of these parameters given by $\mathcal{V},$ $\mathcal{B}_{1},$
$\mathcal{B}_{2}$ and $\mathcal{B}_{3}$. Moreover, the pointwise error in the
left-hand side of (228) splits into two regimes locally around $\epsilon^{*}$.

  1. 1.

If $\epsilon\leq\epsilon^{*}$, the error is dominated by the divergence of the
variance error term $\epsilon^{-(1+d/4)}$.

  2. 2.

If $\epsilon\geq\epsilon^{*}$, the error is dominated by the linear bias error
$\epsilon\mathcal{B}$. A reduction in the magnitude of $\mathcal{B}$ will make
the error diverge more slowly from the minimum at $\epsilon^{*}$.
Consequently, a smaller bias error prefactor will make the error more stable
to small perturbations around $\epsilon^{*}$.

It is difficult to explicitly compute prefactors in Theorem 1.1. Hence it is
not advisable to select the bandwidth using (239). The common practice is to
use a heuristic estimate for $\epsilon^{\ast}$ instead. We use the _kernel
sum_ or the _Ksum_ test [31, 43, 16]:

$\displaystyle\epsilon^{*}_{\sf
Ksum}=\underset{\epsilon}{\text{argmax}}\frac{\partial\log
S(\epsilon)}{\partial\log\epsilon},\quad
S(\epsilon):=\sum_{ij}[K^{(n)}_{\epsilon}]_{ij}.$ (240)

To implement the Ksum test numerically, the sum $S(\epsilon)$ of all entries
of the kernel matrix
$[K^{(n)}_{\epsilon}]_{ij}=\exp(-\|x_{i}-x_{j}\|^{2}/\epsilon)$ is computed at
a range of $\epsilon$ values and its logarithm is plotted against
$\log\epsilon$. The value of $\epsilon$ for which the slope of this plot is
the largest is taken to be $\epsilon^{*}_{\sf Ksum}$. Therefore, at
$\epsilon=\epsilon^{*}_{\sf Ksum}$, the kernel matrix is the most sensitive to
bandwidth for the entire input point cloud. Although this test is not exact,
it gives a good initial guess for a more exhaustive search for the optimal
$\epsilon$. We highlight $\epsilon^{*}_{\sf\sf Ksum}$ in Figures 9, 8, 12, and
11. The RMS error is flatter around $\epsilon^{*}$ and is nearly optimal at
$\epsilon^{*}_{\sf Ksum}$ in all our test problems when we use quasi-uniform
sampling densities.

####  5.2.3 Optimal sampling

The preceding discussion on choosing optimal $\epsilon$ reveals that the
sensitivity of the error to small perturbations around $\epsilon^{*}$ can be
ameliorated through reductions in the bias error prefactors. A quasi-uniform
sampling density can be used to reduce $\mathcal{B}_{2}$ and
$\mathcal{B}_{1}$, thus necessitating increased sampling of rare high-energy
configurations lying in the transition region away from the metastable states.
There are numerous ways of such rare event sampling: some examples relevant to
molecular simulation are temperature acceleration [59], importance sampling
[60], umbrella sampling [61], splitting methods [62], as well as recent work
utilizing generative models such as normalizing flows [63], Generative
Adversarial Networks (GANs) or diffusion models [64, 65]. See [66] for a
recent comprehensive survey on enhanced sampling in molecular dynamics. In the
present case, the situation is complicated because the low-dimensional
manifold $\mathcal{M}$ where $\rho$ is supported is unknown and can only be
accessed through sampling. Consequently, in this setting, metadynamics [67]
emerges as a cheap and robust option to generate samples on $\mathcal{M}$.
These samples can then be postprocessed to $\delta$-nets to enhance spatial
uniformity.

Metadynamics can be described as a modification of the Euler-Maruyama scheme
of discretizing overdamped Langevin dynamics (1). The goal of metadynamics is
to bias the system towards exploring regions of the energy landscape that are
not sufficiently sampled under the original potential $V(x)$. To do so, the
metadynamics algorithm introduces a bias potential $W(x,t)$ that is added to
the original potential $V(x)$ at each step of the Euler-Maruyama algorithm.
The bias potential is constructed as a sum of Gaussian potentials deposited
during the simulation time $t$:

${W(x,t)=\sum_{t_{j}<t}w_{0}\,\exp\left(-\frac{({\theta(x)-\theta(x(t_{j})}))^{2}}{2\sigma^{2}}\right),}$
(241)

where $\theta(\cdot)$ is a vector of collective variables used for biasing the
trajectory, $x(t_{j})$ is the set of atomic positions at time $t_{j}$, $w_{0}$
is the height of each deposited Gaussian potential and $\sigma$ controls the
width of the Gaussians. The bias potential $W(x,t)$ modifies the original
potential energy, leading to an effective potential energy $V_{\sf
eff}(x,t)=V(x)+W(x,t)$. Thus, the metadynamics fills up the free energy wells
with Gaussian bumps enabling the system to escape from metastable states and
explore the configurational space. It is summarized in Algorithm 1.

Input : Gaussian height $w_{0}$, Gaussian width $\sigma_{0}$, deposition
stride $s$, maximum number of steps $n$, timestep $\Delta t$.

Output : Visited states trajectory $\mathcal{X}(n)=\\{x_{1},\ldots x_{n}\\}$

Initialize the visited states trajectory: $\mathcal{X}(n)\leftarrow[\,]$;

Initialize the history of visited states: $H\leftarrow\emptyset$;

Initialize the potential: $V_{\sf eff}\leftarrow V$;

for  _$i\leftarrow 1$ to $N_{\text{steps}}$_ do

Sample a new state: $x_{i}\leftarrow$ Euler-Maruyama step of size $\Delta t$
with potential $V_{\sf eff}$;

Append the current state to the trajectory:
$\mathcal{X}(n)\leftarrow\mathcal{X}(n)\cup\\{x_{i}\\}$;

if  _$i\bmod s=0$_ then

Deposit Gaussian: $V_{\sf eff}\leftarrow V_{\sf
eff}+w_{0}\exp\left(-\frac{(\theta(\cdot)-\theta(x_{i}))^{2}}{2\sigma^{2}}\right)$

end if

end for

Algorithm 1 Metadynamics

The $\delta$-net algorithm outlined in Algorithm 2 is a greedy procedure used
to select a sample of points from the dataset $\mathcal{X}(n)$ such that any
pair of points is at least a distance of $\delta$ apart. The output point
cloud tends to be spatially uniform and approximates samples from the uniform
density on $\mathcal{M}$.

Input : Set of points $P$, threshold distance $\delta$

Output : $\delta$-net $D\subseteq P$

Initialize an empty set $D\leftarrow\emptyset$;

foreach  _$p\in P$_ do

if  _$\text{no point in }D\text{ is within distance }\delta\text{ from }p$_
then

Add $p$ to the $\delta$-net: $D\leftarrow D\cup\\{p\\}$;

end if

end foreach

return _$D$_ ;

Algorithm 2 $\delta$-net Algorithm

####  5.2.4 Experimental setup.

To understand the effect of sampling density $\rho$ and the bandwidth
parameter $\epsilon$ on the quality of the numerical committor $q^{{\sf TMD}}$
we will examine the root mean squared error (RMSE):

$\displaystyle\text{RMSE}\>|q-q^{{\sf
TMD}}|=\left(\sum_{i=1}^{n}|q(x_{i})-q^{{\sf TMD}}(x_{i})|^{2}\right)^{1/2}.$
(242)

To investigate the effect of $\rho$, the underlying point cloud
$\mathcal{X}(n)$ for computing $q^{{\sf TMD}}$ via the TMDmap algorithm will
be drawn from three different sampling densities via the following sampling
algorithms.

  1. 1.

The Euler-Maruyama sampling (see e.g. [68]) leads to the data being sampled
essentially through the invariant Gibbs density $\mu\propto\exp{(-\beta
V(x))}$.

  2. 2.

Metadynamics outlined in Algorithm 1 results in a more uniform distribution of
data on $\mathcal{M}$.

  3. 3.

Metadynamics with $\delta$-net first generates data through Algorithm 1 and
then post-processes them through Algorithm 2. We will illustrate that such
configurations tend to be spatially quasi-uniform on $\mathcal{M}$.

For each choice of sampling density, we will draw a fixed number of points and
feed the resulting point cloud as input to the TMDmap algorithm. To study the
effect of $\epsilon$, we will then vary $\epsilon$ and compute the RMSE as a
function of $\epsilon$ for two test systems governed by the overdamped
Langevin dynamics (1) with Müller’s potential and a two-well potential in
$\mathbb{R}^{2}$. Since both of these examples are 2D, the function
$\theta(\cdot)$ in (241) is chosen to be the identity map.

####  5.2.5 Müller’s potential

We consider a system governed by (1) with $V$ being Müller’s potential [69]:

$V(x_{1},x_{2})=\sum\limits_{i=1}^{4}D_{i}\exp\left\\{a_{i}(x_{1}-X_{i})^{2}+b_{i}(x_{1}-X_{i})(x_{2}-Y_{i})+c_{i}(x_{2}-Y_{i})^{2}\right\\}$
(243)

where

$\displaystyle[a_{1},a_{2},a_{3},a_{4}]$ $\displaystyle=[-1,-1,-6.5,0.7],$
$\displaystyle[b_{1},b_{2},b_{3},b_{4}]$ $\displaystyle=[0,0,11,0.6],$
$\displaystyle[c_{1},c_{2},c_{3},c_{4}]$ $\displaystyle=[-10,-10,-6.5,0.7],$
$\displaystyle[D_{1},D_{2},D_{3},D_{4}]$ $\displaystyle=[-200,-100,-170,15],$
$\displaystyle[X_{1},X_{2},X_{3},X_{4}]$ $\displaystyle=[1,0,-0.5,-1],$
$\displaystyle[Y_{1},Y_{2},Y_{3},Y_{4}]$ $\displaystyle=[0,0.5,1.5,1].$

We computed $q^{{\sf TMD}}$ with the three types of input density specified in
Section 5.2.4. The timestep for the Euler-Maruyama sampling was $\Delta
t=10^{-4}$, and the trajectory length was $10^{6}$ timesteps. The trajectory
was subsampled to keep every $100$th point resulting in a point cloud of size
$10^{4}$. In the metadynamics approach, a Gaussian bump with $\sigma=0.1$ and
$w_{0}=0.5$ was deposited at every $100$th timestep. The dataset was then
post-processed to a $\delta$-net. The reactant set $A$ and product set $B$
were chosen to be the balls of radius $0.1$ centered at $[-0.558,1.441]$ and
$[0.623,0.028]$ respectively. To validate the results, the numerical solution
to the committor problem was computed using the finite element method (FEM).
The computational domain for the FEM was
$\\{x\in\mathbb{R}^{2}\leavevmode\nobreak\ |\leavevmode\nobreak\ V(x)\leq
10\\}$. The homogeneous Neumann boundary condition was imposed on the outer
boundary where $V=10$ (Figure 7). This FEM solution was used as a surrogate
for the true solution in computing the RMSE in (242).

Figure 7: Landscape of Müller’s potential with different sampling densities

Since it is of interest to study the effect of uniformization on the RMSE, we
modulate $\delta$ in the $\delta$-net algorithm used for post-processing the
metadynamics dataset containing $10^{4}$ points. For each $\delta$, the RMSE
over a varying $\epsilon$ is computed and presented in Figure 8. The
$\delta$-net postprocessing tends to result in a flatter error curve around
the minimum hence making the numerical solution more robust with respect to
the choice of the bandwidth.

The TMDmap solutions obtained using datasets generated with the three sampling
types and as well as the FEM mesh consisting of nearly equilateral triangles
of nearly equal sizes are compared in Figure 9. The error curve for the
metadynamics plus $\delta$-net dataset has a flat region around its minimum
and achieves almost as small a minimum as the one for the FEM mesh point
cloud.

Figure 8: We track the RMSE for $q^{\sf TMD}$ computed for Müller’s potential
as a function of $\epsilon$ for different degrees of spatial uniformity. The
dotted line marks the $\epsilon^{*}_{Ksum}$ computed using (240). Figure 9:
Tracking the RMSE for $q^{\sf TMD}$ in Müller’s potential for different
sampling densities as a function of $\epsilon$. Note that using the mesh from
the FEM algorithm gives the least error.

####  5.2.6 Two-well potential in 2D

A 2D test system governed by (1) with a two-well potential

$V(x_{1},x_{2})=(x^{2}_{1}-1)^{2}$ (244)

at $\beta=1.0$ magnifies the advantage of using point clouds generated by
metadynamics plus $\delta$-net even stronger. The sets $A$ and $B$ are balls
of radius $0.1$ centered at $[-1,0]$ and $[1,0]$ respectively. The datasets
are sampled from the Gibbs density, metadynamics, and metadynamics with
$\delta$-net with $\delta=0.1$, $0.02$, $0.03$, and $0.04$. Four of these
datasets are displayed in Figure 10.

Figure 10: The two-well potential landscape with different sampling densities

In Figure 11, the effect of $\delta$-net is quantified by tracking the RMSE as
a function of $\epsilon$ for each value of $\delta$. We find a similar pattern
as for the test problem with Müller’s potential as in Figure 8 that moderate
levels of uniformization (at $\delta=0.02$ in the case of $V$ given by (244))
improves and stabilizes the RMSE at the optimal $\epsilon$.

Figure 11: The RMSE for $q^{{\sf TMD}}$ the test system with the two-well
potential (244) obtained using datasets generated via metadynamics and
metadynamics plus $\delta$-net with $\delta=0.01$, $0.02$, and $0.03$.

In Figure 12, the three choices of sampling density are compared. Adding
uniformization using $\delta$-nets produces the same pronounced effect as for
Müller’s potential in Figure 9. Additionally, since this is a 2D example, it
is also possible to provide to TMDmap the set of points from the nearly
regular FEM mesh used for producing the reference solution. The TMDmap
committor $q^{TMD}$ obtained using the FEM mesh dataset has a notably smaller
error for a broad range of $\epsilon$ values than it is for the other point
clouds used here. This reinforces the hypothesis that uniform sets are error-
optimal for TMDmaps. The caveat, however, is that, in higher dimensions, such
uniform meshing is infeasible to generate and thus practitioners must resort
to randomly sampled data on $\mathcal{M}$.

Figure 12: Tracking the RMSE for different sampling densities as a function of
$\epsilon$. Note that using the mesh from the FEM algorithm gives the least
error.

####  5.2.7 Summary of numerical results

Our numerical experiments in Figures 8, 9, 11, and 12 confirm that:

  1. 1.

Spatial uniformization improves not only the accuracy of the RMSE $|q-q^{{\sf
TMD}}|$ at $\epsilon^{*}$ but also the _stability_. Our error estimate
suggests that this results from the faster bias error convergence due to
quasi-uniform $\rho$.

  2. 2.

Metadynamics with $\delta$-net is a viable strategy to realize the faster bias
error convergence rate.

  3. 3.

If the original data $\mathcal{X}$ has good coverage of the manifold
$\mathcal{M}$ then _deleting_ some points can _improve_ the RMSE as long as
those deletions improve the spatial uniformity of the data. In other words,
reducing the number of data points can enhance the accuracy of the solution!

##  6 Conclusion

In this paper, we have derived sharp error estimates revealing the precise
error rates for the target measure diffusion map. Our results extend the
consistency theory for graph Laplacians to a manifold learning setting
involving density reweighting. The main advantage of incorporating the TMDmap
density reweighting is the free choice of the sampling density, enabling the
use of arbitrary enhanced sampling algorithms for generating the input data.

We have provided a principled approach for tuning the sampling density $\rho$
and kernel bandwidth parameter $\epsilon$ to improve the accuracy of the
TMDmap algorithm because our results are _asymptotically sharp_. The obtained
error formula contains a complete quantification of all prefactors involved in
the bias and variance errors. This allows us to find strategies for reducing
the approximation error when discretizing the backward Kolmogorov operator
$\mathcal{L}$ or when solving a Dirichlet boundary-value problem (BVP) with
this discretization. These formulas illuminate that (a) solving a homogeneous
BVP such as the committor problem and/or (b) sampling $\mathcal{M}$ with a
uniform density are the settings in which bias and variance errors are
attenuated. Importantly, we have demonstrated that these improvements in
accuracy are attainable in practice.

A significant aspect of our work has also been the exploration of uniform
subsampling via $\delta$-nets as a robust and simple technique for reducing
the error in the TMDmap approximation to the committor. The consistency
formula from Theorem 1.1 shows that _a priori_ quasi-uniform densities may
yield pointwise speedups in the convergence of both bias and variance. Here we
obtain a quasi-uniform point cloud with expansive coverage of $\mathcal{M}$
via a simple greedy procedure that subsamples the input dataset into a
$\delta$-net. This substantially improves the stability and accuracy of
TMDmap. Such distance-based quasi-uniform sampling has previously been used
for model reduction [17] and in molecular dynamics applications [70].

Additionally, we have probed the interplay between $\delta$ and $\epsilon$
numerically and illustrated that, for low dimensional cases, the $\delta$–net
tends to yield more favorable scaling between the optimal bandwidth and number
of points beyond the optimal scaling for general densities in Theorem 1.1.
Although our justification for the use of uniform densities is grounded in the
i.i.d. assumption, we find that the scaling between $n$ and $\epsilon$ error
is also improved with $\delta$-nets where the i.i.d. assumption no longer
holds. We aim to address this interplay of $n$, $\epsilon$, and $\delta$ in
future work.

##  7 Acknowledgements

This work was partially supported by AFOSR MURI grant FA9550-20-1-0397.

##  8 Code availability

The code and data for generating the figures in this paper are published on
GitHub [71, 72].

## References

  * [1] R. R. Coifman, S. Lafon, Diffusion maps, Applied and computational harmonic analysis 21 (1) (2006) 5–30.
  * [2] B. Nadler, S. Lafon, R. R. Coifman, I. G. Kevrekidis, Diffusion maps, spectral clustering and reaction coordinates of dynamical systems, Applied and Computational Harmonic Analysis 21 (1) (2006) 113–127.
  * [3] S. B. Kim, C. J. Dsilva, I. G. Kevrekidis, P. G. Debenedetti, Systematic characterization of protein folding pathways using diffusion maps: Application to trp-cage miniprotein, The Journal of chemical physics 142 (8) (2015).
  * [4] R. Banisch, P. Koltai, Understanding the geometry of transport: Diffusion maps for lagrangian trajectory data unravel coherent sets, Chaos: An Interdisciplinary Journal of Nonlinear Science 27 (3) (2017).
  * [5] M. Li, M. A. Erdogdu, Riemannian langevin algorithm for solving semidefinite programs, Bernoulli 29 (4) (2023) 3093–3113.
  * [6] E. P. Hsu, Stochastic analysis on manifolds, no. 38, American Mathematical Soc., 2002.
  * [7] A. Bovier, F. Den Hollander, Metastability: a potential-theoretic approach, Vol. 351, Springer, 2016.
  * [8] M. Belkin, P. Niyogi, Laplacian eigenmaps for dimensionality reduction and data representation, Neural computation 15 (6) (2003) 1373–1396.
  * [9] M. Belkin, P. Niyogi, Convergence of laplacian eigenmaps, Advances in neural information processing systems 19 (2006).
  * [10] X. Cheng, N. Wu, Eigen-convergence of gaussian kernelized graph laplacian by manifold heat interpolation, Applied and Computational Harmonic Analysis 61 (2022) 132–190.
  * [11] C. L. Wormell, S. Reich, Spectral convergence of diffusion maps: Improved error bounds and an alternative normalization, SIAM Journal on Numerical Analysis 59 (3) (2021) 1687–1734.
  * [12] R. Lai, J. Lu, Point cloud discretization of fokker–planck operators for committor functions, Multiscale Modeling & Simulation 16 (2) (2018) 710–726.
  * [13] R. Banisch, Z. Trstanova, A. Bittracher, S. Klus, P. Koltai, Diffusion maps tailored to arbitrary non-degenerate itô processes, Applied and Computational Harmonic Analysis 48 (1) (2020) 242–265.
  * [14] E. Vanden-Eijnden, W. E, Towards a theory of transition paths, Journal of statistical physics 123 (3) (2006) 503–523.
  * [15] A. Laio, M. Parrinello, Escaping free-energy minima, Proceedings of the National Academy of Sciences 99 (20) (2002) 12562–12566.  doi:10.1073/pnas.202427399.
  * [16] L. Evans, M. K. Cameron, P. Tiwary, Computing committors via mahalanobis diffusion maps with enhanced sampling data, The Journal of Chemical Physics 157 (21) (2022) 214107.
  * [17] M. Crosskey, M. Maggioni, Atlas: A geometric approach to learning high-dimensional stochastic systems near manifolds, Multiscale Modeling & Simulation 15 (1) (2017) 110–156.
  * [18] M. Hein, J.-Y. Audibert, U. v. Luxburg, Graph laplacians and their convergence on random neighborhood graphs., Journal of Machine Learning Research 8 (6) (2007).
  * [19] A. Singer, From graph to manifold laplacian: The convergence rate, Applied and Computational Harmonic Analysis 21 (1) (2006) 128–134.
  * [20] K. W. Morton, D. F. Mayers, Linear second order elliptic equations in two dimensions, 2nd Edition, Cambridge University Press, 2005, p. 194–234.  doi:10.1017/CBO9780511812248.008.
  * [21] A. Barducci, M. Bonomi, M. Parrinello, Metadynamics, Wiley Interdisciplinary Reviews: Computational Molecular Science 1 (5) (2011) 826–843.
  * [22] U. Von Luxburg, M. Belkin, O. Bousquet, Consistency of spectral clustering, The Annals of Statistics (2008) 555–586.
  * [23] S. S. Lafon, Diffusion maps and geometric harmonics, Yale University, 2004.
  * [24] M. Hein, J.-Y. Audibert, U. v. Luxburg, From graphs to manifolds–weak and strong pointwise consistency of graph laplacians, in: International Conference on Computational Learning Theory, Springer, 2005, pp. 470–485.
  * [25] X. Cheng, H.-T. Wu, Convergence of graph laplacian with knn self-tuned kernels, A Journal of the IMA 11 (2022) 889–957.  doi:https://doi.org/10.1093/imaiai/iaab019.
  * [26] X. Wang, Spectral convergence rate of graph laplacian, arXiv preprint arXiv:1510.08110 (2015).
  * [27] N. García Trillos, M. Gerlach, M. Hein, D. Slepčev, Error estimates for spectral convergence of the graph laplacian on random geometric graphs toward the laplace–beltrami operator, Foundations of Computational Mathematics 20 (4) (2020) 827–887.
  * [28] A. Singer, H.-T. Wu, Vector diffusion maps and the connection laplacian, Communications on pure and applied mathematics 65 (8) (2012) 1067–1144.
  * [29] P. Hoyos, J. Kileel, Diffusion maps for group-invariant manifolds, arXiv preprint arXiv:2303.16169 (2023).
  * [30] R. R. Coifman, I. G. Kevrekidis, S. Lafon, M. Maggioni, B. Nadler, Diffusion maps, reduction coordinates, and low dimensional representation of stochastic systems, Multiscale Modeling & Simulation 7 (2) (2008) 842–864.
  * [31] A. Singer, R. Erban, I. G. Kevrekidis, R. R. Coifman, Detecting intrinsic slow variables in stochastic dynamical systems by anisotropic diffusion maps, Proceedings of the National Academy of Sciences 106 (38) (2009) 16090–16095.
  * [32] D. Wassermann, M. Descoteaux, R. Deriche, et al., Diffusion maps clustering for magnetic resonance q-ball imaging segmentation, International Journal of Biomedical Imaging 2008 (2008).
  * [33] O. Katz, R. Talmon, Y.-L. Lo, H.-T. Wu, Alternating diffusion maps for multimodal data fusion, Information Fusion 45 (2019) 346–360.
  * [34] A. V. Little, M. Maggioni, L. Rosasco, Multiscale geometric methods for data sets i: Multiscale svd, noise and curvature, Applied and Computational Harmonic Analysis 43 (3) (2017) 504–567.
  * [35] W. K. Allard, G. Chen, M. Maggioni, Multi-scale geometric methods for data sets ii: Geometric multi-resolution analysis, Applied and computational harmonic analysis 32 (3) (2012) 435–462.
  * [36] J. C. Bremer, R. R. Coifman, M. Maggioni, A. D. Szlam, Diffusion wavelet packets, Applied and Computational Harmonic Analysis 21 (1) (2006) 95–112.
  * [37] S. Liang, S. W. Jiang, J. Harlim, H. Yang, Solving pdes on unknown manifolds with machine learning, arXiv preprint arXiv:2106.06682 (2021).
  * [38] H. Antil, T. Berry, J. Harlim, Fractional diffusion maps, Applied and Computational Harmonic Analysis 54 (2021) 145–175.
  * [39] S. W. Jiang, J. Harlim, Ghost point diffusion maps for solving elliptic pdes on manifolds with classical boundary conditions, Communications on Pure and Applied Mathematics 76 (2) (2023) 337–405.
  * [40] M. O. Williams, I. G. Kevrekidis, C. W. Rowley, A data-driven approximation of the koopman operator: Extending dynamic mode decomposition, Journal of Nonlinear Science 25 (2015) 1307–1346.
  * [41] L. Evans, M. K. Cameron, P. Tiwary, Computing committors in collective variables via mahalanobis diffusion maps, arXiv preprint arXiv:2108.08979 (2021).
  * [42] X. Cheng, H.-T. Wu, Convergence of graph laplacian with knn self-tuned kernels, Information and Inference: A Journal of the IMA 11 (3) (2022) 889–957.
  * [43] T. Berry, J. Harlim, Variable bandwidth diffusion kernels, Applied and Computational Harmonic Analysis 40 (1) (2016) 68–96.
  * [44] T. Berry, T. Sauer, Local kernels and the geometric structure of data, Applied and Computational Harmonic Analysis 40 (3) (2016) 439–469.
  * [45] A. W. Long, A. L. Ferguson, Landmark diffusion maps (l-dmaps): Accelerated manifold learning out-of-sample extension, Applied and Computational Harmonic Analysis 47 (1) (2019) 190–211.
  * [46] D. Kohli, A. Cloninger, G. Mishne, Ldle: Low distortion local eigenmaps, The Journal of Machine Learning Research 22 (1) (2021) 12914–12977.
  * [47] H. Li, O. Lindenbaum, X. Cheng, A. Cloninger, Variational diffusion autoencoders with random walk sampling, in: Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part XXIII 16, Springer, 2020, pp. 362–378.
  * [48] A. Singer, R. R. Coifman, Non-linear independent component analysis with diffusion maps, Applied and Computational Harmonic Analysis 25 (2) (2008) 226–239.
  * [49] R. J. Allen, D. Frenkel, P. R. ten Wolde, Simulating rare events in equilibrium or nonequilibrium stochastic systems, The Journal of chemical physics 124 (2) (2006).
  * [50] B. P. Vani, J. Weare, A. Dinner, Computing transition path theory quantities with trajectory stratification, The Journal of Chemical Physics (2022).
  * [51] Y. Khoo, J. Lu, L. Ying, Solving for high-dimensional committor functions using artificial neural networks, Research in the Mathematical Sciences 6 (2019) 1–13.
  * [52] Q. Li, B. Lin, W. Ren, Computing committor functions for the study of rare events using deep learning, The Journal of Chemical Physics 151 (5) (2019).
  * [53] Z. Song, M. K. Cameron, H. Yang, A finite expression method for solving high-dimensional committor problems, arXiv preprint arXiv:2306.12268 (2023).
  * [54] Z. Trstanova, B. J. Leimkuhler, T. Lelièvre, Local and global perspectives on diffusion maps in the analysis of molecular systems, Proceedings. Mathematical, Physical, and Engineering Sciences 476 (2019).
URL https://api.semanticscholar.org/CorpusID:119356107

  * [55] M. Belkin, P. Niyogi, Towards a theoretical foundation for laplacian-based manifold methods, Journal of Computer and System Sciences 74 (8) (2008) 1289–1308.
  * [56] O. G. Smolyanov, H. v. Weizsäcker, O. Wittich, Chernoff’s theorem and discrete time approximations of brownian motion on manifolds, Potential Analysis 26 (1) (2007) 1–29.
  * [57] L. C. Evans, Partial differential equations, American Mathematical Society, Providence, R.I., 2010.
  * [58] K. W. Morton, D. F. Mayers, Numerical solution of partial differential equations: an introduction, Cambridge university press, 2005.
  * [59] L. Maragliano, E. Vanden-Eijnden, A temperature accelerated method for sampling free energy and determining reaction pathways in rare events simulations, Chemical physics letters 426 (1-3) (2006) 168–175.
  * [60] Z. I. Botev, P. L’Ecuyer, B. Tuffin, Markov chain importance sampling with applications to rare event probability estimation, Statistics and Computing 23 (2013) 271–285.
  * [61] J. Kästner, Umbrella sampling, Wiley Interdisciplinary Reviews: Computational Molecular Science 1 (6) (2011) 932–942.
  * [62] R. J. Webber, D. Aristoff, G. Simpson, A splitting method to reduce mcmc variance, arXiv preprint arXiv:2011.13899 (2020).
  * [63] S. Falkner, A. Coretti, S. Romano, P. Geissler, C. Dellago, Conditioning normalizing flows for rare event sampling, arXiv preprint arXiv:2207.14530 (2022).
  * [64] J. Sohl-Dickstein, E. Weiss, N. Maheswaranathan, S. Ganguli, Deep unsupervised learning using nonequilibrium thermodynamics, in: International conference on machine learning, PMLR, 2015, pp. 2256–2265.
  * [65] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, Y. Bengio, Generative adversarial networks, Communications of the ACM 63 (11) (2020) 139–144.
  * [66] J. Hénin, T. Lelièvre, M. R. Shirts, O. Valsson, L. Delemotte, Enhanced sampling methods for molecular dynamics simulations [Article v1.0], Living Journal of Computational Molecular Science 4 (1) (2022) 1583.  doi:10.33011/livecoms.4.1.1583.
  * [67] A. Laio, M. Parrinello, Escaping free-energy minima, Proceedings of the national academy of sciences 99 (20) (2002) 12562–12566.
  * [68] P. E. Kloeden, E. Platen, P. E. Kloeden, E. Platen, Stochastic differential equations, Springer, 1992.
  * [69] K. Müller, L. D. Brown, Location of saddle points and minimum energy paths by a constrained simplex optimization procedure, Theoretica chimica acta 53 (1979) 75–93.
  * [70] J. Rydzewski, M. Chen, O. Valsson, Manifold learning in atomistic simulations: A conceptual review, arXiv preprint arXiv:2303.08486 (2023).
  * [71] S. Sule, L. Evans, M. Cameron, https://github.com/ShashankSule/TMDmaps (2023).
  * [72] M. Cameron, https://github.com/mar1akc/transition_path_theory_FEM_distmesh/tree/main (2023).


meche] Department of Mechanical Engineering, Carnegie Mellon University,
Pittsburgh PA, USA 15213 Joint First Authorship meche] Department of
Mechanical Engineering, Carnegie Mellon University, Pittsburgh PA, USA 15213
Joint First Authorship meche] Department of Mechanical Engineering, Carnegie
Mellon University, Pittsburgh PA, USA 15213 meche] Department of Mechanical
Engineering, Carnegie Mellon University, Pittsburgh PA, USA 15213 [biomed]
Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh PA,
USA 15213 [mld] Machine Learning Department, Carnegie Mellon University,
Pittsburgh PA, USA 15213 IR,NMR,UV

# MOFormer: Self-Supervised Transformer model for Metal-Organic Framework
Property Prediction

Zhonglin Cao  [  Rishikesh Magar  [  Yuyang Wang  [  Amir Barati Farimani
<EMAIL_ADDRESS>[

###### Abstract

Metal-Organic Frameworks (MOFs) are materials with a high degree of porosity
that can be used for applications in energy storage, water desalination, gas
storage, and gas separation. However, the chemical space of MOFs is close to
an infinite size due to the large variety of possible combinations of building
blocks and topology. Discovering the optimal MOFs for specific applications
requires an efficient and accurate search over an enormous number of potential
candidates. Previous high-throughput screening methods using computational
simulations like DFT can be time-consuming. Such methods also require
optimizing 3D atomic structure of MOFs, which adds one extra step when
evaluating hypothetical MOFs. In this work, we propose a structure-agnostic
deep learning method based on the Transformer model, named as MOFormer, for
property predictions of MOFs. The MOFormer takes a text string representation
of MOF (MOFid) as input, thus circumventing the need of obtaining the 3D
structure of hypothetical MOF and accelerating the screening process.
Furthermore, we introduce a self-supervised learning framework that pretrains
the MOFormer via maximizing the cross-correlation between its structure-
agnostic representations and structure-based representations of crystal graph
convolutional neural network (CGCNN) on $>$400k publicly available MOF data.
Using self-supervised learning allows the MOFormer to intrinsically learn 3D
structural information though it is not included in the input. Experiments
show that pretraining improved the prediction accuracy of both models on
various downstream prediction tasks. Furthermore, we revealed that MOFormer
can be more data-efficient on quantum-chemical property prediction than
structure-based CGCNN when training data is limited. Overall, MOFormer
provides a novel perspective on efficient MOF design using deep learning.

###### keywords:

American Chemical Society, LaTeX

##  1 Introduction

Metal-organic frameworks (MOFs) are a type of porous crystalline materials1,
2, which have been extensively researched during the past several decades.
Research interests have been induced by the porous structure and versatile
nature of MOFs on their potential applications such as gas adsorption3, 4, 5,
water harvesting and desalination6, 7, 8, and energy storage9, 10, 11. MOFs
typically consist of several building blocks, including metal-based metal
nodes and organic linkers 12, 13, 4. The assembly of those building blocks
following certain topologies generates the 2-dimensional or 3-dimensional
porous structures of MOFs. Because of the countless possible combinations of
metal nodes, organic linkers, and topologies13, 14, there is a sheer amount of
MOFs with different mechanical properties and surface chemistry. Given the
enormous variety of possible MOF structures, rapidly and inexpensively
selecting the potential top performers for each specific task can be
challenging. High-throughput screening with computational tools such as
molecular simulation 5, 15 or density functional theory (DFT) 16, 17 has been
widely used to evaluate the properties of MOFs. Without the need to
experimentally synthesize MOF structures, those computational tools accelerate
the screening process and allow researchers to screen hundreds of thousands of
hypothetical MOF structures for their performance in different applications.

Recently, machine learning (ML) models have become increasingly popular in the
field of MOF property prediction18, 19, 20, 21, 22, 23, 24. The advantage of
the ML models over the simulation methods is their instantaneous inference of
the properties of MOFs. In contrast, the simulation methods require a
computationally expensive rerun for every new MOF. In the last decade,
multiple large scale MOF dataset are released, including the CoRE MOF 201925,
hypothetical MOFs5, and QMOF26. These datasets contain the atomic structures
of MOFs and their properties like CO_2 absorption and band gap. These datasets
are large enough to train accurate data-driven ML models for the prediction of
MOF properties. Handcrafted geometrical features such as large cavity diameter
and pore limiting diameter have been used as input to a multilayer perceptron
(MLP) to predict MOFs properties 19, 23. Although the training of MLP with a
few layers can be fast, this method suffers from underwhelming accuracy due to
the simplicity of network architecture. Moreover, selecting features requires
extensive domain knowledge from the researchers and optimized 3D structures of
MOFs, thus making this method less generic. Given the aforementioned
drawbacks, a novel method that can achieve high accuracy with a more generic
input of MOFs representations should be pursued. Wang et al. 27 utilize the
crystal graph convolutional neural network (CGCNN) 28 to predict methane
absorption of MOFs. CGCNN is a prevalent model with has an architecture
designed specifically for crystalline materials. It takes the type and the 3D
coordinates of atoms in the crystalline materials as input and constructs a
crystal graph. CGCNN can extract features that encode rich chemical
information through convolution operations on the crystal graph. However, one
drawback of using CGCNN for MOFs property prediction is that it requires the
optimized 3D atomic structures of MOFs which are computationally expensive to
obtain. In addition, some large MOF structures consist of hundred or even
thousand of atoms, thus rendering crystal graph for them can be memory-
inefficient.

Enlightened by the fact that all MOFs are combinations of metal nodes, organic
linkers, and topologies, Bucior et al.29 proposed a text string representation
of MOFs called MOFid. The two core sections of a typical MOFid include the
chemical information of building blocks in the format of SMILES30 and the
topology and the catenation of the MOF structure. The building blocks are
represented by an extensively used string representation of molecules called
SMILES 30. The topology and catenation are each represented by a code adopted
from the Reticular Chemistry Structure Resource (RCSR) database31. Therefore,
MOFid is a concise text string representation of MOFs that preserves the
chemical and the majority of the structural information through topology
encoding. The MOFid text based representation enables the application of
language ML models that take text string as input for MOF property prediction.

In this work, we proposed and developed a Transformer-based language model for
MOF property prediction. Transformer and its variants have become the top
choice for the natural language processing tasks since publication in 2017 by
Vaswani et al32. The multi-head attention mechanism allows the Transformer
model to learn contextual information in a sequence without suffering from
long-range dependency33, 34. With its success in processing long sequential
data, Transformer and its other variants are also adopted for chemistry or
bioinformatics application such as molecular35, 36, 37 and protein38 property
prediction. The Transformer model in our work, named as MOFormer, takes a
modified MOFid as input to make predictions of various MOF properties. The
advantage of this method is that it does not require the 3D atomic structure
of the MOF (structure-agnostic), thus enabling a much faster and more flexible
exploration of the hypothetical MOF space. In practice, pretraining the
Transformer model in a self-supervised manner39, 40, 38, 41, 42 leverages
large quantity of unlabeled data to help the MOFormer learn a more robust
representation of the sequence, and further improve its performance in
downstream tasks. To take advantage of pre-training benefits, we also added a
self-supervised learning framework in which the MOFormer and the CGCNN model
are jointly pretrained with $>$400k MOF structures. Dimensionality reduction
tools are used to visualize the latent representation learned by both models
to provide insight into their performance characteristics. Visualization of
attention weights in MOFormer demonstrates that MOFormer learns MOF
representations based on some key atoms and topology. Lastly, we compared the
data efficiency of models to show which one is a better choice when training
data is limited.

##  2 Methods

###  2.1 MOFid tokenization and Transformer

The MOFormer is built upon the encoder part of the Transformer model that
takes tokenized MOFid as input (Figure. 1a). The MOFid tokenizer is a
customized version of the SMILES tokenizer43. The SMILES strings of all
secondary building units (SBUs) of the MOFid is tokenized by the SMILES
tokenizer, while the topology and catenation section of the MOFid is
separately tokenized based on the topology encoding adopted from RCSR31. The
tokens from both sections are then connected by a separator token “&&”. The
tokenization process followed the BERT39 to add a [CLS] token and a [SEP]
token at the beginning and the end of the sequence to symbolize the start and
the end, respectively. Since the tokenized sequences conform to a fixed length
of 512, sequences longer than the fixed length are truncated and the sequences
shorter than that are padded with special tokens [PAD].

Tokenized sequence is embedded and combined with a positional encoding
(Figure. 1a) to include information about the relative and absolute position
of each token. The position encoding is calculated by:

$\begin{split}\text{PE}_{(pos,2i)}&=\text{sin}(\frac{pos}{10000^{2i/d_{emb}}})\\\
\text{PE}_{(pos,2i+1)}&=\text{cos}(\frac{pos}{10000^{2i/d_{emb}}})\end{split}$
(1)

where $pos$ is the position of the token in the sequence, $i$ is the index of
dimension, and $d_{emb}=512$ is the embedding dimension. The Transformer model
is a deep neural network model built upon the self-attention mechanism
(detailed in the supporting information). Each of the Transformer encoder
layer consists of a multi-head attention layer followed by a simple feed-
forward multilayer perceptron (MLP). Residue connection44 and layer
normalization45 is adopted for both the attention and the MLP. In each head of
the attention layer (Figure. 1b), the input sequence embedding $X$ is
multiplied with three learnable weight vectors $W_{q}$, $W_{k}$, and $W_{v}$
to be converted to the query, key, and value vector ($Q$, $K$, $V$). The
scaled dot-product attention $A$ is then calculated by the equation:

$A=\text{softmax}(\frac{QK^{\top}}{\sqrt{d_{k}}})V$ (2)

where $d_{k}$ is the dimension of $Q$ and $K$ (detailed in the supporting
information). The randomly initialized $W_{q}$, $W_{k}$, and $W_{v}$ vectors
in each head allow the model to learn the contextual information between
tokens in different representation subspaces32. Attention from all heads are
concatenated together and then fed into the MLP for the projected output
embedding, which is of the same size of the input embedding. Given that the
self-mechanism can incorporate the information of the whole sequence into each
one of the token embeddings, theoretically any one of the embeddings can be
used as a representation of the whole sequence. Therefore, we followed the
common practice of related works35, 46, 36, 47 to use the embedding of the
first token, [CLS], for further supervised learning tasks.

Figure 1: (a) The pipeline of the MOFormer model. A MOFid of a MOF
(qmof-2521623 is used as an example) is the input to the model. The MOFid is
converted into tokenized sequence before embedded and applied with positional
encoding. The sequence is then fed into multiple Transformer encoder layers.
The learned embedding of the first token will be used as input to a MLP
regression head for downstream prediction task. (b) A schematic showing the
details of each Transformer encoder layer. Embeddings of the sequence pass
through the multi-head scaled dot-product attention layer and then a MLP.
Residue connection and layer normalization is adopted for both the attention
and the MLP. (c) The self-supervised learning framework with CGCNN and
MOFormer. The 3D structure and the MOFid of the same MOF are fed into the
CGCNN and MOFormer, respectively for representation learning. The MLP head
following each models project the representations into embeddings ($Z_{A}$ and
$Z_{B}$). A cross correlation matrix is then constructed using the embeddings.
Barlow Twins loss is applied to optimize the cross correlation matrix to be as
close as possible to an identity matrix.

###  2.2 Self-supervised pretraining with CGCNN

We introduce a self-supervised learning (SSL) paradigm for MOF representation
learning. We designed the framework by taking into consideration two
modalities of data including the text and graph information. One of the
modalities is the text string representation (MOFid) that encapsulates
building blocks’ stoichiometry and bonds (SMILES) and the topology of the MOF.
The text string information is processed by the MOFormer. One of the
limitations of text string data is the lack of availability of information
about the geometry and the neighborhood of atoms creating an information
bottleneck for the tex string input based models. The structure agnostic
nature of the text string input can prevent the transformer to achieve higher
performance than the graph based models. To mitigate these issues of the
MOFormer framework, we introduce SSL pretraining with CGCNN. We leverage
structure as a modality with CGCNN28. The CGCNN model takes as input the
structure of the MOF. The structure data consists of the atom information
along with the neighborhood information which is critical in property
prediction tasks. To implement the SSL pipeline, we take inspiration from
Crystal Twins (CT) framework48. The CT model makes use of the Barlow Twins
loss function introduced by Zbontar et al49 and SimSiam loss50 functions. In
this work, we use the Barlow Twins loss function on the emebddings generated
from the Transformer and CGCNN encoder. As shown in Figure 1C, we initially
encode both the text string representation and graph representation with their
respective encoders. The MOFormer will encode the text string representation
and the CGCNN will encode the graph representation. We generate an embedding
of size 512 from both the encoders and use it to generate the cross-
correlation matrix (Equation 3). Ideally, we want cross-correlation as close
to the identity matrix as both the representations generated from MOForemer
and CGCNN are essentially capturing the same MOF. The Barlow Twins loss
function, which we used for SSL pre-training (Equation 4), tries to force the
cross-correlation matrix to the identity matrix.

$C_{ij}\overset{\Delta}{=}\frac{\sum_{b}\boldsymbol{Z}_{b,i}^{A}\boldsymbol{Z}_{b,j}^{B}}{\sqrt{(\boldsymbol{Z}_{b,i}^{A})^{2}}{\sqrt{(\boldsymbol{Z}_{b,j}^{B})^{2}}}}$
(3)

where $b$ is the batch index and $i,j$ index the 512 dimensional output from
the projector ($Z^{A}$ and $Z^{B}$), $A$ is the graph representation and $B$
is text representation for the same MOF.

$L_{\mathrm{{BT}}}\overset{\Delta}{=}\sum_{i}(1-C_{ii})^{2}+\lambda\sum_{i}\sum_{j\neq
i}C_{ij}^{2},$ (4)

where $\boldsymbol{C}$ is the cross-correlation matrix of embeddings from the
MOFormer and CGCNN, the cross correlation matrix is given by Equation 3. The
$\lambda$ used in this work is set to $0.0051$.

Finally after pretraining the models using SSL, the encoder weights are shared
during the finetuning stage (Figure S1). The pretraining hyperparamter details
are shown in Table S3 For finetuning, we initialize the model with pretrained
weights and train the model for 200 epochs to generate the final
prediction(Hyperparameters: Table S1 and Table S2). The MOFormer and CGCNN
models are finetuned separately. We observed that using SSL pre-training
framework improved the results of both CGCNN and MOFormer consistently for all
the datasets.

###  2.3 Datasets and other featurizations

Three public MOF datasets including the CORE MOF 201925, the hypothetical
MOFs5 (hMOF), and the Boyd&Woo4 are combined to create a large dataset for the
self-supervised pretraining. The pretraining dataset only includes MOFs with
both 3D structure and MOFid available. After the removal of duplicated MOFs,
the final pretraining dataset has 413535 unique MOFs. In the downstream
prediction task, the MOFormer and the CGCNN are trained on the quantum MOF26
(QMOF) and hMOF in a supervised manner. The QMOF dataset contains 20375 MOFs
each with a label of DFT calculated band gap in eV. Only 7466 MOFs in the QMOF
dataset have MOFid available. On the other hand, the hMOF has 137652 MOFs, in
which 102858 have available MOFid. The models are trained on hMOF with the
label of CO_2 and CH_4 adsorption in mol$\cdot$kg-1 at 0.05, 0.5, and 2.5 bar
of pressure. The benchmark datasets are split into training, test, and
validation sets with a ratio of 0.7-0.15-0.15. During the training, the model
with the best validation performance is recorded and then tested with the test
set. According to the splitting rule, MOFormer has 5226-1119-1119 QMOF data
and 72000-15428-15428 hMOF data, while other models have 14262-3056-3056 QMOF
data and 96356-20647-20647 hMOF data in the training, validation and test set,
respectively. Although the MOFs with available MOFid form a subset of both
benchmark datasets, the subset with MOFid shares the same distribution and has
approximately the same mean and standard deviation compared with the original
whole dataset (Figure S2 and S3 in the supporting information). Therefore, it
is fair to compare the performance of MOFormer and other models.

We also benchmarked the MOFormer and CGCNN against other non-DL based
featurization methods such as the Smooth Overlap of Atomic Positions51, 52, 53
(SOAP) and the Stoichiometric-12054 features. SOAP is a structure-based
featurization method and the Stoichiometric-120 is a structure-agnostic
featurization method. The parameters used for creating SOAP features are
included in the supporting information Table S4. XGBoost55 model is used to
make predictions using those handcrafted features.

##  3 Results and discussion

###  3.1 QMOF

The first dataset we benchmark models on is the QMOF dataset, where the label
for each MOF is DFT calculated band gap. Lower band gap value results in
better conductivity of the MOF. Accurate prediction of band gap can help to
identify conductive MOFs which is useful in energy storage applications56, 57.
The accuracy of models follows the rank of CGCNN $>$ MOFormer $>$ SOAP $>$
Stoichiometric-120 (Table. 1). MOFormer has 21.2$\%$ lower MAE compared with
the other structure-agnostic method Stoichiometric-120. It is worth noting
that structure-agnostic MOFormer outperforms structure-based SOAP with a
smaller size of training set, indicating that MOFormer is capable of
extracting critical features from the MOFid for energy-related property
prediction. The self-supervised pretraining helps reduce the mean absolute
error (MAE) of CGCNN by 6.79$\%$ and MOFormer by 5.34$\%$. The reduced error
proves the improvement brought by the pretraining.

CGCNNscratch CGCNNpretrain SOAP MOFormerscratch MOFormerpretrain
Stoichiometric-120 MAE (eV) 0.275 $\pm$0.015 0.256$\pm$0.006 0.424$\pm$0.007
0.387$\pm$0.001 0.367$\pm$0.005 0.466$\pm$0.011

Table 1: Benchmark performance of different models on the band gap prediction
of the QMOF dataset. Mean absolute error (MAE, in the unit of eV) and standard
deviation of 3 runs of different initial seeds of each model is reported. The
left three models are structure-based and the right three models are
structure-agnostic. The best performance of each category is marked as bold.
Figure 2: (a) The binned scatter plot shows the comparison between predicted
and DFT calculated band gap for MOFs in the QMOF dataset. MOFs included in
this figure are from the randomly splitted test set. Darker color of each
hexagonal bin represents more data points are in the bin. Dashed line
represents perfect prediction. (b) - (c) Visualization58 of the MOF structure
with lowest (qmof-1c923ff) and the second lowest (qmof-6bda2bd) band gap in
the test set. The bar plot shows the comparison between predictions made by
different models.

To better understand the superior performance of MOFormer and CGCNN in QMOF,
we trained the 4 models with the same training set and then examined their
performance on the same test set consisting 1119 randomly selected MOFs. The
binned scatter plot (Figure 2a) shows the comparison between the predicted and
the DFT calculated band gap. Darker color means more data fall in the bin.
More predictions made by CGCNN and MOFormer are closer to the ground truth,
especially for MOFs with band gap lower than $\leq$ 2eV. The SOAP and
Stoichiometric-120 are more likely to overpredict the lower band gap. This
weakness of SOAP and Stoichiometric-120 can also be confirmed by the kernel
density estimation of predicted values (Figure S4). The MOFs with the top-2
lowest band gaps in this test set are the qmof-1c923ff (0.03eV) and
qmof-6bda2bd (0.039eV). Band gap prediction by MOFormer and CGCNN is much
closer to the DFT calculated value than predictions by SOAP and
Stoichiometric-120 (Figure 2b-c), especially for qmof-6bda2bd. Accurately
predicting low band gap of MOFs can lead to the discovery of conductive MOF,
rendering pretrained MOFormer and CGCNN more valuable for prescreening MOFs.

###  3.2 hMOF

The models are also benchmarked on the hMOF dataset with the label of CO_2 and
CH_4 adsorption under 0.05, 0.5, and 2.5 bar of pressure. Table 2 shows that
pretrained MOFormer is constantly outperforming the other structure-agnostic
Stoichiometric-120 by achieving 35$\%$-48$\%$ lower MAE. Pretrained CGCNN
outperforms other models for the CO_2 adsorption prediction. The pretraining
in average improves the accuracy of MOFormer by 4.3$\%$ and the CGCNN by
16.5$\%$ over all gas adsorption predictions. SOAP has surprisingly low MAE
for the gas adsorption prediction, outperforming pretrained CGCNN for 2 of 3
CH_4 adsorption predictions, and the CGCNN trained from scratch for all gas
adsorption predictions. The outstanding performance of SOAP on hMOF can be
attributed to the low variation of elements included in the hMOF dataset. Only
11 different elements are present among all 137652 hMOFs, which is very
limited compared to 79 in the QMOF dataset. Less number of elements results in
much smaller and less sparse SOAP feature vector (SOAP feature has an size of
2772 for hMOFs and 19908 for QMOFs with same parameters), thus leading to the
high prediction accuracy of the following XGBoost regressor. However, the high
accuracy of SOAP can hardly be sustainable when it is used in exploring more
diverse hypothetical MOFs. When more elements are included in the dataset, the
SOAP feature vector size and sparsity increase drastically, causing the data
and model to be too large to be accommodated by the memory of local machines
and drop of prediction accuracy (Table S5). MOFormer and CGCNN will not suffer
from such an issue since their inputs remain invariant with increasing types
of elements in the dataset, making them better choices when exploring more
diverse chemical space for MOFs.

CO_2 0.05bar CO_2 0.5 bar CO_2 2.5 bar CH_4 0.05bar CH_4 0.5bar CH_4 2.5bar
CGCNNscratch 0.126$\pm$0.005 0.391$\pm$0.017 0.818$\pm$0.050 0.028$\pm$0.001
0.121$\pm$0.006 0.333$\pm$0.017 CGCNNpretrain 0.110$\pm$0.001 0.330$\pm$0.002
0.645$\pm$0.003 0.025$\pm$0.001 0.099$\pm$0.001 0.258$\pm$0.008 SOAP
0.115$\pm$0.002 0.339$\pm$0.004 0.666$\pm$0.003 0.022$\pm$0.001
0.106$\pm$0.001 0.239$\pm$0.002 MOFormerscratch 0.178$\pm$0.002
0.558$\pm$0.001 1.000$\pm$0.013 0.034$\pm$0.000 0.174$\pm$0.002
0.385$\pm$0.003 MOFormerpretrain 0.158$\pm$0.001 0.545$\pm$0.008
0.982$\pm$0.011 0.033$\pm$0.000 0.161$\pm$0.011 0.384$\pm$0.003
Stoichiometric-120 0.282$\pm$0.002 0.983$\pm$0.005 1.895$\pm$0.003
0.050$\pm$0.001 0.269$\pm$0.001 0.631$\pm$0.002

Table 2: Benchmark performance of different models on gas adsorption
prediction of the hMOF dataset. Mean absolute error (mol kg-1) and standard
deviation of 3 runs of different initial seeds of each model is reported. The
top three models are structure-based and the bot three models are structure-
agnostic. The best performance of each category is marked as bold.

The representations of MOFs learned by the MOFormer and CGCNN after finetuning
are visualized to provide interpretablity to the models (Figure 3). Each
representation is projected to the 2D space using the dimension reduction tool
t-SNE59. t-SNE clusters more similar data points together while placing less
similar data points further away. Only MOFs which has the top-10 most common
topologies in hMOF is included in Figure 3 because they take $>99.7\%$ of the
whole dataset. We can observe that CGCNN representations cluster MOF with high
CO_2 adsorption more closely than MOFormer representations by comparing Figure
3a and 3c. This contributes to higher prediction accuracy of CGCNN. On the
other hand, MOFormer representation clusters MOFs with the same topology
closer than CGCNN representation does. For example, the MOFs with dia (green)
and tbo (brown) topologies form two clusters in the lower left corner of
MOFormer representation visualization (Figure 3b). Those MOFs are much loosely
clustered in the CGCNN representation visualization (Figure 3d). The MOFormer
representation focusing on topology can be caused by the fact that gas
adsorption is more dependent on the 3D structure of the MOF compared with its
atom composition. The only structure-related information contained in the
MOFid is the topology encoding. Therefore, more weights are on the topology
after MOFormer is finetuned to predict the gas adsorption of MOFs. The input
of CGCNN is the 3D structure of MOFs with atomic resolution, thus causing
CGCNN to rely less on the topology for gas adsorption prediction.

Figure 3: The t-SNE59 dimension reduced visualization of MOF representations
learned by (a)-(b) the MOFormer and (c)-(d) the CGCNN. Each data point in (a)
and (c) is colored by its CO_2 adsorption at 0.5 bar pressure, and in (b) and
(d) is colored by its topology. Only MOF which has the top-10 most common
topologies in hMOF dataset is shown.

###  3.3 Visualization of attention weights

Figure 4 demonstrates the attention maps between tokens of a MOFid (qmof-
ba40858) in the last MOFormer layer after finetuning on band gap prediction.
The attention map can serve as a visual interpretation of how the MOFormer
learns MOF representations60. We observe a strong attention in head 5 from all
tokens to the metal node Ytterium, and to the topology encoding pcu in head 1
and 3. The attention from metal node to the topology encoding is especially
high in head 1. Moreover, large attention weight can be observed between
tokens in the SMILES of the SBUs in head 6. Head 1, 3, 5, and 8 also show
large attention weights on the carbon and the oxygen atom, and the double bond
in the organic building block. The attention weight visualization shows that
MOFormer learns a representation that emphasizes on the contextual information
between key components in the MOFid including important atoms (e.g. Y, O, and
C) and the topology, thus leading to more accurate prediction.

Figure 4: Heatmap of the attention between tokens (MOFid of qmof-ba40858) in
different heads of the last MOFormer layer. If we index each block in the
heatmap as $\text{block}_{n,i,j}$, where $i$, $j$, and $n$ is the row, column,
and head index, respectively. The value of $\text{block}_{n,i,j}$ represents
the attention on the $j$-th token from the $i$-th token in $n$-th heatmap.

###  3.4 Data efficiency comparison

Obtaining high-quality MOF data using experimental or DFT method can be time-
consuming and expensive. A model with high data efficiency is ideal when the
training data size is limited. We compared the data efficiency of different
models on the QMOF and hMOF dataset (CO_2 adsorption at 0.5 bar pressure). For
band gap prediction (Figure 5a), the pretrained MOFormer outperforms CGCNN
when the training set size $\leq 1000$. This makes MOFormer more valuable in
predicting quantum-chemical properties when the training dataset is difficult
to build (i.e. experimentally synthesized MOFs). Both MOFormer and CGCNN
achieve higher accuracy than SOAP regardless of the training set size on the
QMOF dataset. For CO_2 adsorption prediction (Figure 5b), CGCNN constantly
achieves higher accuracy than MOFormer regardless of the training set size,
indicating its higher data efficiency. CGCNN outperforms MOFormer on hMOF
because the CO_2 adsorption correlates more with the MOF structure and the
input to CGCNN provides more structural information than MOFid. SOAP achieves
higher data efficiency than CGCNN and MOFormer on hMOF, but is eventually
caught up by CGCNN after the training set size exceeds 50000. Both Figure 5a
and 5b show that pretraining consistently improves the data efficiency of
MOFormer and CGCNN. Moreover, SOAP is shown to have diminishing improvement
with increasing training set size, but CGCNN and MOFormer do not suffer from
such an issue.

Figure 5: Data efficiency comparison between different models on the (a) QMOF
and the (b) hMOF dataset. The models are trained on a subset of the training
set, while the validation and test set are kept the same. The training set
sizes are 100, 500, 1000, 5000, 10000, and 50000 (hMOF only). Since only less
than 7500 MOFs in QMOF have available MOFid, the maximum training subset size
for MOFormer on QMOF is 5000.

##  4 Conclusion

In summary, we propose a Transformer-based model, named as MOFormer, for
structure-agnostic MOF property prediction. Taking only MOFid as input, the
MOFormer model is expected to expedite the exploration of hypothetical MOFs.
We also introduce a self-supervised learning framework to jointly pretrain the
MOFormer and CGCNN model on large unlabeled MOF dataset to enhance their
prediction accuracy in downstream tasks. Compared with another structure-
agnostic method Stoichiometric-120, MOFormer achieves $21.4\%$ higher accuracy
on band gap prediction, and $35\%$-$48\%$ higher accuracy on various gas
adsorption prediction tasks. MOFormer even outperforms structure-based SOAP
method in band gap prediction with less training data. The pretraining is
further shown to improve the accuracy of MOFormer by $5.34\%$ and $4.3\%$ in
average and CGCNN by $6.79\%$ and $16.5\%$ in average, for band gap and gas
adsorption prediction, respectively. MOFormer and CGCNN are shown to be less
likely to overpredict the band gap of MOFs compared with SOAP and
Stoichiometric-120, making them better choices for prescreening conductive
MOFs for energy applications. When used for gas adsorption prediction of MOFs,
MOFormer relies more on the topology information compared with CGCNN because
of the strong correlation between the label and the structure of MOF.
Visualization of the attention weights in the last MOFormer layer reveals that
the attention layers in MOFormer focus more on several important atoms and the
topology to learn the representation of a MOF. Lastly, MOFormer is shown to be
more data efficient than CGCNN for band gap prediction when the training set
size $\leq 1000$. As a structure-agnostic model, MOFormer can make rapid and
accurate inference on the property of MOFs (especially for quantum-chemical
properties) using arbitrarily constructed MOFid as input. Therefore, MOFormer
can serve as a tool for exploring the vast chemical space of hypothetical
MOFs.

This work is supported by the start-up fund from Mechanical Engineering
Department at CMU.

Transformer and self-attention mechanism. Details of CGCNN and MOFormer model.
Details of the self-supervised pretraining. Distribution of QMOF and hMOF
dataset. Parameters for SOAP feature vector creation and the effect of SOAP
vector length to model accuracy. Kernel density estimation of band gap
prediction from different models.

##  5 Code Availability

The python code of this work can be found on Github
https://github.com/zcao0420/MOFormer.

## References

  * James 2003 James, S. L. Metal-organic frameworks. _Chemical Society Reviews_ 2003, _32_ , 276–288
  * Zhou et al. 2012 Zhou, H.-C.; Long, J. R.; Yaghi, O. M. Introduction to metal–organic frameworks. _Chemical reviews_ 2012, _112_ , 673–674
  * 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. _Nature communications_ 2019, _10_ , 1–9
  * Boyd et al. 2019 Boyd, P. G.; Chidambaram, A.; García-Díez, E.; Ireland, C. P.; Daff, T. D.; Bounds, R.; Gładysiak, A.; Schouwink, P.; Moosavi, S. M.; Maroto-Valer, M. M., et al. Data-driven design of metal–organic frameworks for wet flue gas CO 2 capture. _Nature_ 2019, _576_ , 253–256
  * Wilmer et al. 2012 Wilmer, C. E.; Leaf, M.; Lee, C. Y.; Farha, O. K.; Hauser, B. G.; Hupp, J. T.; Snurr, R. Q. Large-scale screening of hypothetical metal–organic frameworks. _Nature chemistry_ 2012, _4_ , 83–89
  * Almassad et al. 2022 Almassad, H. A.; Abaza, R. I.; Siwwan, L.; Al-Maythalony, B.; Cordova, K. E. Environmentally adaptive MOF-based device enables continuous self-optimizing atmospheric water harvesting. _Nature communications_ 2022, _13_ , 1–10
  * Hanikel et al. 2020 Hanikel, N.; Prévot, M. S.; Yaghi, O. M. MOF water harvesters. _Nature nanotechnology_ 2020, _15_ , 348–355
  * Cao et al. 2019 Cao, Z.; Liu, V.; Barati Farimani, A. Water desalination with two-dimensional metal–organic framework membranes. _Nano letters_ 2019, _19_ , 8638–8643
  * Baumann et al. 2019 Baumann, A. E.; Burns, D. A.; Liu, B.; Thoi, V. S. Metal-organic framework functionalization and design strategies for advanced electrochemical energy storage devices. _Communications Chemistry_ 2019, _2_ , 1–14
  * Zhao et al. 2016 Zhao, Y.; Song, Z.; Li, X.; Sun, Q.; Cheng, N.; Lawes, S.; Sun, X. Metal organic frameworks for energy storage and conversion. _Energy storage materials_ 2016, _2_ , 35–62
  * Xu et al. 2017 Xu, G.; Nie, P.; Dou, H.; Ding, B.; Li, L.; Zhang, X. Exploring metal organic frameworks for energy storage in batteries and supercapacitors. _Materials today_ 2017, _20_ , 191–209
  * Sharp et al. 2021 Sharp, C. H.; Bukowski, B. C.; Li, H.; Johnson, E. M.; Ilic, S.; Morris, A. J.; Gersappe, D.; Snurr, R. Q.; Morris, J. R. Nanoconfinement and mass transport in metal–organic frameworks. _Chemical Society Reviews_ 2021,
  * Moosavi et al. 2020 Moosavi, S. M.; Nandy, A.; Jablonka, K. M.; Ongari, D.; Janet, J. P.; Boyd, P. G.; Lee, Y.; Smit, B.; Kulik, H. J. Understanding the diversity of the metal-organic framework ecosystem. _Nature communications_ 2020, _11_ , 1–10
  * Falcaro et al. 2011 Falcaro, P.; Hill, A. J.; Nairn, K. M.; Jasieniak, J.; Mardel, J. I.; Bastow, T. J.; Mayo, S. C.; Gimona, M.; Gomez, D.; Whitfield, H. J., et al. A new method to position and functionalize metal-organic framework crystals. _Nature communications_ 2011, _2_ , 1–8
  * Ozcan et al. 2020 Ozcan, A.; Semino, R.; Maurin, G.; Yazaydin, A. O. Modeling of gas transport through polymer/MOF interfaces: a microsecond-scale concentration gradient-driven molecular dynamics study. _Chemistry of Materials_ 2020, _32_ , 1288–1296
  * Canepa et al. 2013 Canepa, P.; Arter, C. A.; Conwill, E. M.; Johnson, D. H.; Shoemaker, B. A.; Soliman, K. Z.; Thonhauser, T. High-throughput screening of small-molecule adsorption in MOF. _Journal of Materials Chemistry A_ 2013, _1_ , 13597–13604
  * Nandy et al. 2021 Nandy, A.; Duan, C.; Taylor, M. G.; Liu, F.; Steeves, A. H.; Kulik, H. J. Computational discovery of transition-metal complexes: from high-throughput screening to machine learning. _Chemical Reviews_ 2021, _121_ , 9927–10000
  * Fung et al. 2021 Fung, V.; Zhang, J.; Juarez, E.; Sumpter, B. G. Benchmarking graph neural networks for materials chemistry. _npj Computational Materials_ 2021, _7_ , 1–8
  * Burner et al. 2020 Burner, J.; Schwiedrzik, L.; Krykunov, M.; Luo, J.; Boyd, P. G.; Woo, T. K. High-Performing Deep Learning Regression Models for Predicting Low-Pressure CO2 Adsorption Properties of Metal–Organic Frameworks. _The Journal of Physical Chemistry C_ 2020, _124_ , 27996–28005
  * Altintas et al. 2021 Altintas, C.; Altundal, O. F.; Keskin, S.; Yildirim, R. Machine learning meets with metal organic frameworks for gas storage and separation. _Journal of Chemical Information and Modeling_ 2021, _61_ , 2131–2146
  * Lee et al. 2021 Lee, S.; Kim, B.; Cho, H.; Lee, H.; Lee, S. Y.; Cho, E. S.; Kim, J. Computational screening of trillions of metal–organic frameworks for high-performance methane storage. _ACS Applied Materials & Interfaces_ 2021, _13_ , 23647–23654
  * Choudhary et al. 2022 Choudhary, K.; Yildirim, T.; Siderius, D. W.; Kusne, A. G.; McDannald, A.; Ortiz-Montalvo, D. L. Graph neural network predictions of metal organic framework CO2 adsorption properties. _Computational Materials Science_ 2022, _210_ , 111388
  * Moghadam et al. 2019 Moghadam, P. Z.; Rogge, S. M.; Li, A.; Chow, C.-M.; Wieme, J.; Moharrami, N.; Aragones-Anglada, M.; Conduit, G.; Gomez-Gualdron, D. A.; Van Speybroeck, V., et al. Structure-mechanical stability relations of metal-organic frameworks via machine learning. _Matter_ 2019, _1_ , 219–234
  * Nandy et al. 2021 Nandy, A.; Duan, C.; Kulik, H. J. Using Machine Learning and Data Mining to Leverage Community Knowledge for the Engineering of Stable Metal–Organic Frameworks. _Journal of the American Chemical Society_ 2021, _143_ , 17535–17547
  * Chung et al. 2019 Chung, Y. G.; Haldoupis, E.; Bucior, B. J.; Haranczyk, M.; Lee, S.; Zhang, H.; Vogiatzis, K. D.; Milisavljevic, M.; Ling, S.; Camp, J. S., et al. Advances, updates, and analytics for the computation-ready, experimental metal–organic framework database: CoRE MOF 2019. _Journal of Chemical & Engineering Data_ 2019, _64_ , 5985–5998
  * Rosen et al. 2021 Rosen, A. S.; Iyer, S. M.; Ray, D.; Yao, Z.; Aspuru-Guzik, A.; Gagliardi, L.; Notestein, J. M.; Snurr, R. Q. Machine learning the quantum-chemical properties of metal–organic frameworks for accelerated materials discovery. _Matter_ 2021, _4_ , 1578–1597
  * Wang et al. 2020 Wang, R.; Zhong, Y.; Bi, L.; Yang, M.; Xu, D. Accelerating Discovery of Metal–Organic Frameworks for Methane Adsorption with Hierarchical Screening and Deep Learning. _ACS Applied Materials & Interfaces_ 2020, _12_ , 52797–52807
  * Xie and Grossman 2018 Xie, T.; Grossman, J. C. Crystal graph convolutional neural networks for an accurate and interpretable prediction of material properties. _Physical review letters_ 2018, _120_ , 145301
  * Bucior et al. 2019 Bucior, B. J.; Rosen, A. S.; Haranczyk, M.; Yao, Z.; Ziebel, M. E.; Farha, O. K.; Hupp, J. T.; Siepmann, J. I.; Aspuru-Guzik, A.; Snurr, R. Q. Identification schemes for metal–organic frameworks to enable rapid search and cheminformatics analysis. _Crystal Growth & Design_ 2019, _19_ , 6682–6697
  * Weininger 1988 Weininger, D. SMILES, a chemical language and information system. 1. Introduction to methodology and encoding rules. _Journal of chemical information and computer sciences_ 1988, _28_ , 31–36
  * O’Keeffe et al. 2008 O’Keeffe, M.; Peskov, M. A.; Ramsden, S. J.; Yaghi, O. M. The reticular chemistry structure resource (RCSR) database of, and symbols for, crystal nets. _Accounts of chemical research_ 2008, _41_ , 1782–1789
  * Vaswani et al. 2017 Vaswani, A.; Shazeer, N.; Parmar, N.; Uszkoreit, J.; Jones, L.; Gomez, A. N.; Kaiser, Ł.; Polosukhin, I. Attention is all you need. Advances in neural information processing systems. 2017; pp 5998–6008
  * Bahdanau et al. 2014 Bahdanau, D.; Cho, K.; Bengio, Y. Neural machine translation by jointly learning to align and translate. _arXiv preprint arXiv:1409.0473_ 2014,
  * Hochreiter and Schmidhuber 1997 Hochreiter, S.; Schmidhuber, J. Long short-term memory. _Neural computation_ 1997, _9_ , 1735–1780
  * Schwaller et al. 2019 Schwaller, P.; Laino, T.; Gaudin, T.; Bolgar, P.; Hunter, C. A.; Bekas, C.; Lee, A. A. Molecular transformer: a model for uncertainty-calibrated chemical reaction prediction. _ACS central science_ 2019, _5_ , 1572–1583
  * Schwaller et al. 2021 Schwaller, P.; Probst, D.; Vaucher, A. C.; Nair, V. H.; Kreutter, D.; Laino, T.; Reymond, J.-L. Mapping the space of chemical reactions using attention-based neural networks. _Nature Machine Intelligence_ 2021, _3_ , 144–152
  * Xu et al. 2022 Xu, C.; Wang, Y.; Farimani, A. B. TransPolymer: a Transformer-based Language Model for Polymer Property Predictions. _arXiv preprint arXiv:2209.01307_ 2022,
  * Elnaggar et al. 2021 Elnaggar, A.; Heinzinger, M.; Dallago, C.; Rehawi, G.; Wang, Y.; Jones, L.; Gibbs, T.; Feher, T.; Angerer, C.; Steinegger, M., et al. ProtTrans: towards cracking the language of lifes code through self-supervised deep learning and high performance computing. _IEEE transactions on pattern analysis and machine intelligence_ 2021,
  * Devlin et al. 2018 Devlin, J.; Chang, M.-W.; Lee, K.; Toutanova, K. Bert: Pre-training of deep bidirectional transformers for language understanding. _arXiv preprint arXiv:1810.04805_ 2018,
  * Liu et al. 2019 Liu, Y.; Ott, M.; Goyal, N.; Du, J.; Joshi, M.; Chen, D.; Levy, O.; Lewis, M.; Zettlemoyer, L.; Stoyanov, V. Roberta: A robustly optimized bert pretraining approach. _arXiv preprint arXiv:1907.11692_ 2019,
  * Wang et al. 2022 Wang, Y.; Wang, J.; Cao, Z.; Barati Farimani, A. Molecular contrastive learning of representations via graph neural networks. _Nature Machine Intelligence_ 2022, _4_ , 279–287
  * Wang et al. 2022 Wang, Y.; Magar, R.; Liang, C.; Barati Farimani, A. Improving Molecular Contrastive Learning via Faulty Negative Mitigation and Decomposed Fragment Contrast. _Journal of Chemical Information and Modeling_ 2022,
  * Schwaller et al. 2018 Schwaller, P.; Gaudin, T.; Lanyi, D.; Bekas, C.; Laino, T. “Found in Translation”: predicting outcomes of complex organic chemistry reactions using neural sequence-to-sequence models. _Chemical science_ 2018, _9_ , 6091–6098
  * He et al. 2016 He, K.; Zhang, X.; Ren, S.; Sun, J. Deep residual learning for image recognition. Proceedings of the IEEE conference on computer vision and pattern recognition. 2016; pp 770–778
  * Ba et al. 2016 Ba, J. L.; Kiros, J. R.; Hinton, G. E. Layer normalization. _arXiv preprint arXiv:1607.06450_ 2016,
  * Schwaller et al. 2021 Schwaller, P.; Hoover, B.; Reymond, J.-L.; Strobelt, H.; Laino, T. Extraction of organic chemistry grammar from unsupervised learning of chemical reactions. _Science Advances_ 2021, _7_ , eabe4166
  * Dosovitskiy et al. 2020 Dosovitskiy, A.; Beyer, L.; Kolesnikov, A.; Weissenborn, D.; Zhai, X.; Unterthiner, T.; Dehghani, M.; Minderer, M.; Heigold, G.; Gelly, S., et al. An image is worth 16x16 words: Transformers for image recognition at scale. _arXiv preprint arXiv:2010.11929_ 2020,
  * Magar et al. 2022 Magar, R.; Wang, Y.; Farimani, A. B. Crystal Twins: Self-supervised Learning for Crystalline Material Property Prediction. _arXiv preprint arXiv:2205.01893_ 2022,
  * Zbontar et al. 2021 Zbontar, J.; Jing, L.; Misra, I.; LeCun, Y.; Deny, S. Barlow twins: Self-supervised learning via redundancy reduction. International Conference on Machine Learning. 2021; pp 12310–12320
  * Chen and He 2021 Chen, X.; He, K. Exploring simple siamese representation learning. Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. 2021; pp 15750–15758
  * Bartók et al. 2013 Bartók, A. P.; Kondor, R.; Csányi, G. On representing chemical environments. _Physical Review B_ 2013, _87_ , 184115
  * Bartók et al. 2017 Bartók, A. P.; De, S.; Poelking, C.; Bernstein, N.; Kermode, J. R.; Csányi, G.; Ceriotti, M. Machine learning unifies the modeling of materials and molecules. _Science advances_ 2017, _3_ , e1701816
  * Himanen et al. 2020 Himanen, L.; Jäger, M. O. J.; Morooka, E. V.; Federici Canova, F.; Ranawat, Y. S.; Gao, D. Z.; Rinke, P.; Foster, A. S. DScribe: Library of descriptors for machine learning in materials science. _Computer Physics Communications_ 2020, _247_ , 106949
  * Meredig et al. 2014 Meredig, B.; Agrawal, A.; Kirklin, S.; Saal, J. E.; Doak, J. W.; Thompson, A.; Zhang, K.; Choudhary, A.; Wolverton, C. Combinatorial screening for new materials in unconstrained composition space with machine learning. _Physical Review B_ 2014, _89_ , 094104
  * Chen and Guestrin 2016 Chen, T.; Guestrin, C. Xgboost: A scalable tree boosting system. Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining. 2016; pp 785–794
  * Xie et al. 2020 Xie, L. S.; Skorupskii, G.; Dincă, M. Electrically conductive metal–organic frameworks. _Chemical reviews_ 2020, _120_ , 8536–8580
  * Sheberla et al. 2017 Sheberla, D.; Bachman, J. C.; Elias, J. S.; Sun, C.-J.; Shao-Horn, Y.; Dincă, M. Conductive MOF electrodes for stable supercapacitors with high areal capacitance. _Nature materials_ 2017, _16_ , 220–224
  * Dubbeldam et al. 2018 Dubbeldam, D.; Calero, S.; Vlugt, T. J. iRASPA: GPU-accelerated visualization software for materials scientists. _Molecular Simulation_ 2018, _44_ , 653–676
  * Van der Maaten and Hinton 2008 Van der Maaten, L.; Hinton, G. Visualizing data using t-SNE. _Journal of machine learning research_ 2008, _9_
  * Vig 2019 Vig, J. A multiscale visualization of attention in the transformer model. _arXiv preprint arXiv:1906.05714_ 2019,

##  6 Supporting Information

##  7 Self-attention and the Transformer model

In this work, we use the encoder part of the Transformer model as the base for
the MOFormer framework. Transformer32 model was proposed as a language model
for sequence-to-sequence language translation. The encoder part of it learns a
latent representation of the input sequence using the self-attention
mechanism. Self-attention mechanism is previous used along with sequential
models (such as recurrent neural network33 and LSTM34) to combat the drop in
model performance with long sequence. Transformer model drop the recurrent
architecture in previous sequential models, and rely solely on the self-
attention mechanism to learn the representation of the input sequence. The
advantage of such model architecture is threefold. First, self-attention has
lower per-layer computational complexity than recurrent architectures. Second,
models can be highly parallelizable with only self-attention (as self-
attention can be calculated using matrix multiplication) thus can be trained
much faster. Third, long-range dependencies in the input sequence can be
easily learned using self-attention mechanism.

In the self-attention mechanism, each token in the sequence is linearly
projected into its corresponding query, key, and value vector ($Q$, $K$, and
$V$ matrix for the entire sequence, respectively). The weights that linearly
project the entire input sequence are denoted as the $W_{q}$, $W_{k}$, and
$W_{v}$ matrices. In this work, followed by Vaswani et al.32, the scaled dot-
product attention:

$Attention(Q,K,V)=\text{softmax}(\frac{QK^{\top}}{\sqrt{d_{emb}}})V$ (5)

is adopted, where $d_{emb}=512$. Multihead attention is to have multiple
linear projection matrices in parallel for the calculation of attention. The
$W_{q}$, $W_{k}$, and $W_{v}$ in each attention head is initialized
differently, thus multihead attention enables the model to learn
representation of the input sequence from different subspaces. Assuming we
have $h$ attention heads, the multihead attention is calculated as:

$\begin{split}MultiheadAttn(Q,K,V)&=\text{Concat}(head_{1}\dots
head_{h})W_{o}\\\
\text{where}\;head_{i}&=Attention(Q_{i},K_{i},V_{i})\end{split}$ (6)

$W_{o}$ is a matrix to linearly project the concatenated attention from each
head; $Q_{i}$, $K_{i}$, and $V_{i}$ is the query, key, and value matrix of the
$i$-th head, respectively. In this work, the number of head $h=8$. The
dimension of the $W_{q}^{i}$, $W_{k}^{i}$, and $W_{v}^{i}$ (linear projection
matrices in the $i$-th head) is $W_{q}^{i}\in\mathbb{R}^{d_{emb}\times
d_{k}}$, $W_{k}^{i}\in\mathbb{R}^{d_{emb}\times d_{k}}$,
$W_{v}^{i}\in\mathbb{R}^{d_{emb}\times d_{v}}$, where $d_{emb}=512$ and
$d_{v}=d_{k}=d_{emb}/h=64$. Moreover, $W_{o}\in\mathbb{R}^{d_{emb}\times
d_{emb}}$. The Transformer encoder in this work contains 6 self-attention
layers (Figure 1b in the manuscript), in which the MLP has a dimension of 512.

##  8 Training Details

In this section, we describe in details the hyperparameters that we used for
training the models. Our framework consists of MOFormer (transformer), CGGNN
and the SSL pretraining done to improve the performance of both the GNN and
the transformer model. During the finetuning stage, a MLP regression head is
attached to the encoder part (MOFormer or CGCNN) to make prediction. The
regression head has 4 layers, with 512, 256, 128 and 64 neurons in each layer,
respectively.

###  8.1 MOFormer

The MOFormerFinetuning hyperparameters are shown in Table S1. Different
learning rate (LR) are used for the encoder part (Transformer encoder) and the
MLP regression head. We also train the MOFormerScratch to evaluate the
improvement in performance after SSL pre-training.

Hyperparameters | Batch Size | LREncoder | LRMLP | Optimizer | Epochs | Embed. Size | Train/Val | Weight Decay
---|---|---|---|---|---|---|---|---
MOFormerScratch | 64 | 0.00005 | 0.01 | Adam | 200 | 512 | 0.7/0.15/0.15 | $10^{-6}$
MOFormerFinetuning | 64 | 0.00005 | 0.01 | Adam | 200 | 512 | 0.7/0.15/0.15 | $10^{-6}$

Table S1: Finetuning hyperparameters for the MOFormer model. LR indicate the
learing rate

###  8.2 CGCNN

The CGCNN model is structure based model that takes the “.cif” file as input.
The MOFormer is a structure agnostic method, however, information like
neighbourhood of atoms and geometry of the crystal are often lost in the text
representation. To ensure we also have a structure based benchmark to predict
MOF properties we also evaluate the performance of the CGCNN model28. The
CGCNN model with its information rich features often outperforms the structure
agnostic MOFormer. We also pretrain the MOFormer and the CGCNN together using
the self-supervised learning. The performance of pretrained CGCNNPretrain
model is found to be higher than CGCNNScratch model trained from Scratch
without pretraining. The hyperparamters we used for the model are shown in
Table S2.

Hyperparameters | Batch Size | Number of layers | LR | Optimizer | Epochs | Embed. Size | Train/Val/Test | Weight Decay
---|---|---|---|---|---|---|---|---
CGCNNScratch | 128 | 3 | 0.01 | Adam | 200 | 512 | 0.7/0.15/0.15 | $10^{-6}$
CGCNNFinetuning | 128 | 3 | 0.002 | Adam | 200 | 512 | 0.7/0.15/0.15 | $10^{-6}$

Table S2: Finetuning hyperparameters for the CGCNN model. LR denotes the
learning rate for the models.

###  8.3 SSL - Pretraining

Structure agnostic models often lack the geometry and atom neighborhood
information leading to an information bottleneck for these models. To mitigate
this issue, we develop the SSL pre-training framework. In the SSL pre-training
framework we have the MOFormer branch and the CGCNN branch, the goal of the
SSL pre-training is to ensure that the representations from both the MOFormer
and the CGCNN are similar. The pre-training helps both the MOFormer model and
CGCNN in a way that it tries to leverage the embeddings learnt by the other
model. The model benefit mutually from each other and the performance
enhancements can be seen for both models during the finetuning stage (Figure
S1). The hyperparamters we used during pretraining the model are shown in
Table S3

Hyperparameters | Batch Size | Learning Rate | Optimizer | Epochs | Embed. Size | Train/Val
---|---|---|---|---|---|---
CGCNNPretrain | 32 | 0.00001 | Adam | 15 | 512 | 0.95/0.05
MOFormerPretrain | 32 | 0.00001 | Adam | 15 | 512 | 0.95/0.05

Table S3: Hyperparameters for Pretraining the model. Figure S1: The SSL
pretraining and finetuning framework. In this work, we use self supervised
learning to pre-train both the CGCNN model and the MOFormer using self-
supervised learning. The pre-trained weights are then shared during finetuning
for down stream tasks. Using such a technique help us improve the perfromance
of both the CGCNN and MOFormer.

##  9 Distribution of the benchmark datasets

Since only part of MOF in each QMOF26 and hMOF5 dataset have available MOFid,
we plot the distribution of the whole datasets and the subsets with MOFid.
Figure S2 shows that the difference between the average band gap value between
the whole dataset and the subset is as small as 0.135. The distribution
between the data are also similar.

Figure S3 shows that the difference between the average CO_2/CH_4 adsorption
value between the whole dataset and the subset is very small. The distribution
between the data are also similar. Therefore the benchmark results comparisons
between MOFormer and other models on both datasets are fair.

Figure S2: Histogram showing the distribution of QMOF dataset. Blue bars are
all MOFs in the QMOF dataset. Orange bars are MOFs with MOFid available. Mean
$\mu$ and standard deviation $\sigma$ are reported. Figure S3: Histogram
showing the distribution of hMOF dataset with the label of (a) CO_2 adsoprtion
and (b) CH_4 adsorption at 0.05, 0.5, and 2.5 bar of pressure. Blue bars are
all MOFs in the hMOF dataset. Orange bars are MOFs with MOFid available. Mean
$\mu$ and standard deviation $\sigma$ are reported.

##  10 Parameters to create SOAP feature vector

The SOAP51 feature vector is created using the DScribe package53. Table S4
records the parameters for creating SOAP feature vector, which are chosen
arbitrarily based on previous works18, 26.

rcut | $\sigma$ | nmax | lmax | average method | crossover
---|---|---|---|---|---
8 | 0.2 | 8 | 6 | Inner average | False
Table S4: Parameters used for creating SOAP feature vector

The length of the SOAP feature vector varies greatly based on the number of
different atom elements in the dataset. When the MOFs in a dataset includes a
large variation of elements (such as QMOF, 79 elements), the SOAP feature
vector will be very long and sparse. However, when the MOFs in the dataset
consist of only a few different elements (such as hMOF, 11 elements), the SOAP
feature vector will be relatively short and less sparse. For example, SOAP
feature vector created using the 11 elements in hMOF has a length of 2772. The
SOAP feature vector created using the 79 elements in QMOF has a length of
19908, given that nmax and lmax parameters are kept the same. The high
accuracy of SOAP feature vector on hMOF is a direct result of the small
variation of element type in the dataset. Assuming that we create the SOAP
feature vectors for hMOFs using double the original number of elements (22
instead of the original 11), the length of SOAP feature vector will be doubled
to 5544. The performance of XGBoost using the new SOAP feature is
significantly reduced (Table S5). Therefore, we can conclude that SOAP feature
is not a suitable choice when exploring the chemical space for MOFs where
elements are diverse.

Model |  $\\#$ elements | SOAP vector length |  hMOF 2.5 bar CH_4 adsorption
---|---|---|---
SOAP | 11 | 2772 |  0.239$\pm$0.002
SOAP | 22 | 5544 |  0.288$\pm$0.002
CGCNNpretrain | - | - |  0.258$\pm$0.008
Table S5: Parameters used for creating SOAP feature vector

##  11 Comparison between models’ performance in QMOF dataset

Figure S4 shows the kernel density estimation (KDE) plot of the predicted band
gap values of MOFs from different models. The MOFs are in a randomly selected
test set from the QMOF dataset. The ground truth is the DFT calculated band
gap, which is shown by the blue curve. SOAP and Stoichiometric-120 features
with XGBoost model underpredict low bandgap values more frequently than CGCNN
and MOFormer.

Figure S4: (a) The kernel density estimation (KDE) plot of the predicted band
gap values of MOFs from different models. (b) Zoom-in view of the distribution
of MOFs with band gap $\leq$ 1eV