Universitat Politècnica de València

# Effect of framework composition and NH3 on the diffusion of Cu+ in Cu-CHA
catalysts predicted by machine-learning accelerated molecular dynamics

Reisel Millan  Department of Materials Science and Engineering, Massachusetts
Institute of Technology, Cambridge, MA 02139  Estefanía Bello-Jurado
Instituto de Tecnología Química, Universitat Politècnica de València-Consejo
Superior de Investigaciones Científicas, Avenida de los Naranjos s/n, 46022
Valencia, Spain  Manuel Moliner  Instituto de Tecnología Química, Universitat
Politècnica de València-Consejo Superior de Investigaciones Científicas,
Avenida de los Naranjos s/n, 46022 Valencia, Spain  Mercedes Boronat
<EMAIL_ADDRESS>Instituto de Tecnología Química, Universitat Politècnica de
València-Consejo Superior de Investigaciones Científicas, Avenida de los
Naranjos s/n, 46022 Valencia, Spain  Rafael Gomez-Bombarelli<EMAIL_ADDRESS>Department of Materials Science and Engineering, Massachusetts Institute of
Technology, Cambridge, MA 02139

###### Abstract

Cu-exchanged zeolites rely on mobile solvated Cu+ cations for their catalytic
activity, but the role of framework composition on transport is not fully
understood. Ab initio molecular dynamics simulations can provide quantitative
atomistic insight but are too computationally expensive to explore large
length- and time-scales or diverse compositions. We report a machine-learning
interatomic potential that accurately reproduces ab initio results and
effectively generalizes to allow multi-nanosecond simulations of large
supercells and diverse chemical compositions. Biased and unbiased simulations
of [Cu(NH3)2]+ mobility show that aluminum pairing in eight-membered rings
accelerates local hopping, and demonstrate that increased NH3 concentration
enhances long-range diffusion. The probability of finding two [Cu(NH3)2]+
complexes in the same cage - key for SCR-NOx reaction - increases with Cu
content and Al content, but does not correlate with the long-range mobility of
Cu+. Supporting experimental evidence was obtained from reactivity tests of
Cu-CHA catalysts with controlled chemical composition.

###### keywords:

NNP, diffusion, SCR-NOx

Copper-exchanged zeolites play a crucial role as redox catalysts for some
environmentally relevant processes, like the partial methane oxidation to
methanol or the selective catalytic reduction of nitrogen oxides with ammonia
(NH3-SCR-NOx). In both cases, the small pore Cu-SSZ-13 zeolite with the CHA
structure has been reported as an efficient catalyst. 1, 2, 3, 4, 5, 6, 7, 8,
9, 10, 11

The NH3-SCR-NOx reaction is currently employed for the removal of nitrogen
oxides (NOx) from exhaust gases in diesel vehicles and stationary plants,
through a redox catalytic cycle in which Cu+ is oxidized to Cu^2+ by NO2 or
NO+O2, and then reduced to Cu+ by reaction of NH3 and NO forming harmless N2 +
H2O (Scheme Effect of framework composition and NH3 on the diffusion of Cu+ in
Cu-CHA catalysts predicted by machine-learning accelerated molecular
dynamics). 12, 13, 14, 15 This understanding of the reaction mechanism has
enabled development of optimized catalysts by tuning framework topology,
composition and copper speciation. In the as-prepared catalysts, Cu+ and Cu^2+
cations are directly coordinated to the zeolite framework forming
heterogeneous active sites, while under reaction conditions NH3 solvates the
Cu+ cations forming mobile [Cu(NH3)2]+ complexes that act as dynamic active
sites, resembling homogeneous catalysts but within the confinement of the
zeolite pores. At low temperature (T $<$ 523 K) the oxidation step involves
transient dimeric Cu+(NH3)2-O2-Cu+(NH3)2 species whose formation requires the
simultaneous presence of two [Cu(NH3)2]+ monomers in the same cha cage. The
hops between adjacent cha cages are modulated by size exclusion effects and
also by the attractive interaction between the positively charged [Cu(NH3)2]+
complexes and the negatively charged framework Al sites8, 9, 15, 16, 17. Thus,
structural properties like Al content and distribution, Cu loading or Brønsted
acid site density, as well as the interaction of the Cu active sites with the
reactants, in particular NH3, might affect the mobility of Cu cations and
consequently the NH3-SCR-NOx reaction rate. This has been evidenced by recent
studies combining catalytic activity tests with operando XAS or EPR
spectroscopy, 18, 19, 20, 21, 22 and ab initio molecular dynamics (AIMD)
simulations have been successfully applied to provide atomistic insight into
the dynamic nature of the Cu+ cations under reaction conditions9, 16, 17.

[H]  Illustration of the low temperature SCR-NOx redox cycle.

The cost of AIMD simulations limits their applicabiltiy to a few selected
systems at a time, at small length- and time- scales. The timescale limitation
may be partially bypassed with enhanced sampling methods, such as umbrella
sampling (US), which have been used to study the slow diffusion of copper
complexes in CHA.9, 17 However, the systematic exploration of parameters such
as Si/Al ratio, Al distribution, Cu/Al ratio, NH3 concentration or the
presence of Brønsted acid sites and compensating NH4+ cations has not been yet
possible.

Machine learning (ML) has demonstrated broad applicability in materials
science23, 24 and heterogeneous catalysis. 25, 26, 27, 28 Machine learning
potentials (MLPs), when trained with a sufficiently large and diverse dataset,
can match the accuracy of quantum chemistry methods at a fraction of the
computational cost29, 30, 31, 32, 33, 34, 35. This allows the study of larger
and more realistic systems and more complex scientific problems 36, 24, 37,
38, in particular those requiring the use of molecular dynamics simulations
39, 40, 41. A broad variety of MLPs based on neural networks - so called
neural network potentials (NNPs) - have been developed in the last years: ANI
42, 43, 44, 45, Deep Tensor Neural Networks 46, SchNet47, DeepPotentialNet 48,
MEGNet49, DimeNet50, OrbNet51, PaiNN52, NequIP 53, and successfully used to
study solids systems54, 37, 55, 56, ion diffusion 57, and chemical reactions
58, 59, 60.

Here, we leveraged these innovations and trained an NNP capable of describing
[Cu(NH3)2]+ species in aluminosilicate CHA with varying composition and NH3
concentration. The trained NNP proved accurate and transferable, and acquiring
all the training data was less costly than one traditional AIMD simulation.
Biased MD simulations reproduced free energy profiles from DFT, and provided
insight into transport for over a dozen combinations of Al distribution and
NH4+ presence. Unbiased MD simulations were scaled to thousands of atoms for
nanoseconds, and achieved a more realistic representation of the importance of
Al density and distribution, Cu loading, and adsorbed NH3 on the mobility of
Cu+ cations in Cu-CHA catalysts.

These results show that the activation free energy for [Cu(NH3)2]+ hops
between adjacent cages is lower for windows containing Al pairs, but also that
this is a local effect with only weak influence on long-range mobility.
[Cu(NH3)2]+ migration to remote cages requires the simultaneous displacement
of charge-compensating NH4+ which show a lower mobility that is enhanced by
excess NH3. Lastly, simulations with large supercells show that the
probability of finding two [Cu(NH3)2]+ complexes in the same cage, a pre-
requisite for the SCR-NOx reaction, increases with Cu loading, and also with
the Al content in the zeolite. We confirm these trends experimentally through
catalytic tests of Cu-CHA samples with controlled Si/Al and Cu/Al ratios.

##  1 Results

###  1.1 Neural network potential

NNPs are highly accurate but they struggle to extrapolate outside their
training data. In order to ensure robust and accurate production simulations,
our NNP was trained on data gathered through multiple generations of active
learning (AL) using a query by-committee approach 61, 62, 63, 64, 65, 66, 67,
68, 69. A committee (ensemble) of NNPs was trained on the available labeled
data at each iteration, and new data was collected based on the disagreement
(variance) of the prediction of the committee members on newly generated
geometries, as described in Computational Details Section and illustrated in
Figure 1a. The last generation of the NNP was trained on a complete dataset
containing 42k revPBE+D3 force calculations on structural models with a
diverse set of atomic local environments, ranging from 290 to 323 atoms per
supercell, summarized in Table S1 (see structural models in Figure S1). The
chemical compositions in the dataset (Figure 1b) cover a range of Si/Al ratio
from 47 to 13. All DFT calculations in the training data were electrically
neutral with negative charges arising from Al substitution being compensated
by [Cu(NH3)2]+, NH4+, or H+70, 71, 72.

Figure 1: Neural Network Potential. (a) Illustration of the active learning
cycle. (b) Distribution of the chemical composition of the dataset. (c)
Correlation between the predicted and target energies (left) and forces
(right) of the last generation potential.

The active learning strategy was capable of automatically adding new, diverse,
and informative chemical environments to the training pool at each of the pre-
selected compositions through a combination of MD and uncertainty
quantification. It generated informative training data for a number of
chemical processes that occur during the reaction, but were not present in the
initial training data. These include adsorption and protonation of NH3 on the
Brønsted acid sites to form NH4+ cations, exchange between a gas phase NH3
molecule and one of the two NH3 ligands of the [Cu(NH3)2]+ complex, and proton
transfer from NH4+ to NH3. The diffusion of [Cu(NH3)2]+ complexes through the
8R windows that connect adjacent cha cages has a higher activation barrier.
Therefore representative training data was obtained through the same enhanced
sampling approach as the production simulations (Figure S2).

This strategic combination of biased MD with uncertainty quantification
allowed efficient sampling of the relevant regions on the PES with a small and
diverse number of DFT evaluations. Figure S3 illustrates the structural
diversity in the final dataset by means of a 2D projection of the local
chemical environments around each Al atom in our data using UMAP73 on the
feature vectors learned by the NNP 74. Atoms with similar local environments
have similar feature vectors and appear close on the UMAP plot. The overlap
among the chemical compositions suggests a nearly continuous sampling of the
Al local environment.

Figure 1d shows the correlation between predicted and target energies and
forces for a held-out test set. The mean absolute error of the predicted
energies and forces are 0.98 kcal/mol and 1.2 kcal/mol/Å respectively,
indicating that the NNP is capable of predicting the energies and forces with
chemical accuracy.

###  1.2 Effect of Al distribution on [Cu(NH3)2]+ diffusion through 8R windows
from biased simulations

The favorable speed of the NNP accelerates US MD simulations by orders of
magnitude over DFT, and enabled systematic exploration of the role of Al
distribution in well converged simulations. Ten different structural models
with H12Al2Cu2N4O192Si94 composition were built (Fig 2a), each containing two
framework Al atoms compensated with two [Cu(NH3)2]+, and distributed either in
the same 8R (SR1, SR2, SR3 and SR4), in different 8Rs (DR1, DR2, DR3 and DR4),
in the same 4R (S4R) or in the same 6R (S6R).

Figure 2: Structural models and accuracy of NNP biased simulations. (a)
Representation of the different Al pair distributions used in the biased
simulations. The atoms in the ring containing the Al pair are highlighted. Si,
O and Al atoms are depicted in orange, red and light brown, respectively. (b)
Comparison of free energy profiles for [Cu(NH3)2]+ diffusion through the 8R
window of SR1 (orange), SR4 (blue) and DR2 (green) systems obtained from DFT
(left) and NNP (right) based biased simulations at 423 K.

Because the error of the predicted forces is not a sufficient metric for true
performance in production simulations 75, the NNPs were further validated by
comparing NNP and DFT free energy profiles for [Cu(NH3)2]+ diffusion through
the 8R window of DR2, SR1 and SR4. Due to the high computational cost of
producing reference DFT biased simulations, a smaller hexagonal 126-atom unit
cell was used. The free energy profiles for SR1, SR4 and DR2 (Figure 2b) and
the calculated activation free energies ($\Delta G_{\mathrm{act}}$), 5.2 and
5.3 kcal/mol for SR1, 4.8 and 4.8 kcal/mol for SR4, and 6.1 and 6.4 kcal/mol
for DR2, are in excellent agreement at both computational levels.

The smoother NNP traces are a consequence of more abundant sampling (4.8 ns
total, compared to 0.8 for DFT) given the advantageous computational cost of
NNP (20 ps of AIMD required over a week on CPU cores (Intel(R) Xeon(R) CPU
E5-2650 v3 @ 2.30GHz), compared to 20 min on a Tesla V100-32 GB GPU for NNP).

Then, the free energy profiles for [Cu(NH3)2]+ diffusion between neighboring
cages in the ten systems depicted in Figure 2a were obtained from NNP biased
simulations at 423 K using the larger H12Al2Cu2N4O192Si94 models. The profiles
are plotted in Figure 3 and the corresponding values of activation ($\Delta
G_{\mathrm{act}}$) and reaction ($\Delta G$) free energies are summarized in
Table S2. The shaded area in each profile shows the standard error calculated
from three independent simulations using three different NNPs trained on the
same dataset. The average value of the standard error, 0.1 kcal/mol in all
cases, indicates a low uncertainty in the prediction of the free energy and
well-converged statistics.

Figure 3: NNP biased simulations of [Cu(NH3)2]+ diffusion. (a) Free energy
profiles for the diffusion of one [Cu(NH3)2]+ complex from a cage A to a
neighboring cage B occupied by another [Cu(NH3)2]+ complex and (b) snapshots
of the final state at $\xi$ = 2.5 corresponding to the diffusion of
[Cu(NH3)2]+ through the 8R windows of SR1, DR4, S4R and S6R systems at 423 K.

For the systems with two Al atoms in the same 8R (orange profiles in Figure
3), $\Delta G_{\mathrm{act}}$ values range from 3.9 to 5.4 kcal/mol and the
reaction is slightly endergonic with $\Delta G$ values between 0.4 and 2.3
kcal/mol. In all other cases, $\Delta G_{\mathrm{act}}$ are higher than 6
kcal/mol and $\Delta G$ are larger than 3 kcal/mol, with the only exception of
the DR4 system for which the process is slightly exergonic. The Al
distribution in the DR4 model is the same as in SR1, but the diffusion occurs
through different 8R windows (see snapshots in Figure 3). In both cases, the
stability of the final state with the two [Cu(NH3)2]+ in the same cage is
similar to that of the initial state with the two complexes in different
cages, which suggests that this particular Al distribution might favor the
formation of the [Cu(NH3)2]+-OO-[Cu(NH3)2]+ dimers involved in the low-
temperature NH3-SCR-NOx reaction. In contrast, distributions with two Al in
the same 4R or 6R hinder the formation of such dimeric intermediates, because
only one of the two [Cu(NH3)2]+ can stay near the Al atoms (see snapshots in
Figure 3) while the second [Cu(NH3)2]+ is forced to remain in too close
contact with the first, resulting in its diffusion through a different 8R to
another empty cage. These results demonstrate that Al distribution affects the
movement of [Cu(NH3)2]+ species between cages and the stability of pairs of
[Cu(NH3)2]+ complexes in the same cage, and points to a positive effect of Al
pairs in 8R on the rate of the low-temperature NH3-SCR-NOx reaction.

###  1.3 Effect of NH4+ on the diffusion of [Cu(NH3)2]+ through 8R windows
from biased simulations

The mobility of [Cu(NH3)2]+ complexes within the zeolite microporous structure
is affected by the presence of other molecules involved in the reaction, among
which NH3 is the most abundant and the one with the largest impact on
diffusion and reactivity.16, 17, 20, 21, 76, 77 Under reaction conditions NH3
is readily protonated on the Brønsted acid sites forming NH4+ cations that
remain coordinated to the framework AlO4- units and might partly block the
diffusion of [Cu(NH3)2]+ complexes through the 8R windows. To analyze this
possibility we first performed NNP biased simulations at 423 K using the
previously described H10Al2Cu1N3O192Si94 models with two framework Al atoms
compensated now with one [Cu(NH3)2]+ complex and one NH4+ cation. The free
energy profiles obtained for the three Al distributions considered (SR1, SR2
and SR3), plotted in Figure 4 are clearly different from those depicted in
Figure 3, and the calculated activation free energies for [Cu(NH3)2]+
diffusion are also reflective of how the strong coordination of NH4+ modifies
transport.

(a) SR1 ($\xi$ = -2.5)

(b) SR1 ($\xi$ = 2.5)

(c) SR2 ($\xi$ = -4.8)

(d) SR2 ($\xi$ = 2.5)

(e) SR3 ($\xi$ = -4.8)

(f) SR3 ($\xi$ = 2.0)

Figure 4: Effect of NH4+ on the diffusion of [Cu(NH3)2]+ complex. Free energy
profiles for [Cu(NH3)2]+ diffusion between neighboring cages in the presence
of NH4+ in three different systems obtained from NNP-based biased simulations
at 423 K, and snapshots of the systems SR1, SR2 and SR3 at the initial and
final minimum states. Si, O, Al, H, Cu and N atoms are depicted in orange,
red, light brown, white, green and blue, respectively.

Whereas in the previous simulations without NH4+ the most stable minimum for
the initial state occurs at $\xi$ $\sim$ -2.5 Å, with [Cu(NH3)2]+ relatively
close to the 8R, that is the case only for SR1 in the presence of NH4+. In SR2
and SR3 models the most stable minimum lies at $\xi$ $\sim$ -4.8 Å, with the
[Cu(NH3)2]+ complex closer to the center of the cavity and farther away from
the 8R to be crossed and therefore from the negatively charged AlO4- units.
The calculated $\Delta G_{\mathrm{act}}$, 5.2, 7.6 and 13.1 kcal/mol for SR1,
SR2 and SR3, respectively, and $\Delta G$ values, 3.2, 5.0 and 11.4 kcal/mol
for SR1, SR2 and SR3, respectively, are higher than those obtained for the
corresponding systems in the absence of NH4+. The snapshots of the final state
at $\xi$ $\sim$ 2.5 Å for SR1 and SR2 in Figure 4 show that the NH4+ cation
has been displaced from its initial position in the plane of the 8R to a
position relatively close to one of the AlO4- units. In SR3, however, the NH4+
cation has been displaced to the opposite side of the cage, far from the two
AlO4- sites present in the model, which would explain the instability of the
system. The deviation across replicate profiles is wider in some regions,
which we attribute to the fact that our simulations occasionally, but not
exhaustively, sample the spontaneous reversible deprotonation of NH4+ cations
to form NH3 and a Brønsted acid site, which is typically not accessible to
traditional simulations.

Altogether, the results from the biased simulations suggest a potential
blocking effect of NH4+ cations. However, their own mobility, as either NH4+
or NH3 following proton transfer to deprotonated AlO4-, and possible migration
from the 8R towards other nearby AlO4- units that are not present in this
model might modify this conclusion. The NNPs developed here open the
possibility of running long-time unbiased MD simulations on larger systems
with more realistic chemical compositions, allowing to capture the dynamics of
[Cu(NH3)2]+ and NH4+ globally and to observe long-range diffusion of both
cationic species.

###  1.4 Long-range diffusion of NH4+ and [Cu(NH3)2]+ species from unbiased MD
simulations

A large T768O1536 supercell of dimension $\sim$37 Å was used to construct
eight models representing three Al contents, low (L), medium (M) and high (H),
two copper loadings (4 and 20 atoms per simulation box), and varying NH3, NH4+
and H+ concentrations and Al pairings as described in Methods Section, Fig 5
and Table S3. Unbiased MD simulations at 500 K were conducted for at least 5
ns on each 2000-atom system.

L4

L20

M4

M20

(a) H6R

(b) H8R

(c) HR

(d) HB

Figure 5: Structural models used in NNP unbiased MD simulations. Snapshots of
the input structures corresponding to the eight different zeolite framework
compositions considered. The first letter of the label indicates the Al
content, L (low) Si/Al=28.5, M (medium) Si/Al=14.3 and H (high) Si/Al=10.3.
The number in the top models indicates the Cu loading, 4 or 20 atoms per unit
cell. The second letter in the bottom models that contain 20 Cu atoms per unit
cell indicates the Al distribution, forming pairs in 6R or 8R, random (R) and
biased (B). Si, O, Al, H, Cu and N atoms are depicted as orange, red, light
brown, white, green and blue.

Figures 6 a and b track, respectively, the time evolution of the mean square
displacements (MSDs) of the N atoms in the NH4+ cations and of the Cu atoms in
the [Cu(NH3)2]+ complexes. While both species have a net charge of +1,
[Cu(NH3)2]+ is much more mobile while NH4+ cations are more closely attached
to the AlO4- units (Figs S4 and S5). This is due to the positive charge in
[Cu(NH3)2]+ complexes being highly shielded by the two NH3 ligands.

The similarity among the MSD profiles of N atoms suggests that the mobility of
NH4+ cations is rather independent of zeolite framework composition and Cu
content (Figure 6a), while the MSD traces for Cu (Figure 6b) suggest slightly
lower mobility of [Cu(NH3)2]+ in the systems with high Al content (Si/Al
$\sim$ 10, Figure 6b). The trends are similar in the number of distinct cages
visited by the [Cu(NH3)2]+ complex over time (Figure 6c), but they are more
stratified and more clearly show an increase in mobility with decreasing
copper content (orange and red dashed lines higher than solid in Fig 6c). In 5
ns, each [Cu(NH3)2]+ complex visits on average fewer than three different
cages in the models with high Al content, which increases to 3-4 for
intermediate Al, and reaches a maximum of over 5 different cages visited for
the L4 model, which has the lowest Al content and thus few Al pairs in 8R.
This is in apparent contradiction to the biased simulations results for small
models, that showed lower hopping free energy barrier for Al pairing in 8R. A
potential explanation is that the hopping landscape is statically and
dynamically heterogeneous, with the local chemical environment of each initial
and final cages, and transient cage occupation by other mobile molecules
influencing the mobility of copper complexes. Fig S4 shows the diversity in
length and tortuosity in example trajectories of Cu atoms inside the zeolite
microporous structure. While MSD profiles representing the average movement of
the [Cu(NH3)2]+ species are relatively similar across catalyst models
($\sim$40-60 Å2 ), the local mobility of each individual [Cu(NH3)2]+ complex
depends on its local chemical environment.

Figure 6: Mobility of NH4+ and [Cu(NH3)2]+. (a) Mean square displacement (MSD)
of N atoms in NH4+ (b) MSD of Cu atoms in the [Cu(NH3)2]+ complexes obtained
from NPP-based unbiased MD simulations at 500 K. (c) Number of distinct cages
visited by the Cu atoms. The apparent contradiction with the MSD plots is
explained by the fact that the diffusion of [Cu(NH3)2]+ does not occur along a
preferential direction. One Al pair in a 8R can accelerate the diffusion of
[Cu(NH3)2]+ between two neighboring cages in multiple forward and backward
steps, so that only two distinct cages are visited by this complex. (d)
average minimum Cu-Al distances. All plots were obtained from the same
simulations. MSD was calculated using eq 1.

Since all framework AlO4- units are compensated with either [Cu(NH3)2]+ or
NH4+, lower Cu content implies larger NH4+ concentration for a given Si/Al
ratio. Comparison of L4 with L20, and M4 with M20 in Figure 6c suggests that
NH4+ assists the long-range diffusion of [Cu(NH3)2]+. A proposed explanation
is that adsorbed NH4+ shields the attractive interaction between AlO4- anionic
sites and the [Cu(NH3)2]+ complexes. Each NH4+ forms two strong hydrogen bonds
with the AlO4- site, hence the lower mobility of NH4+ while [Cu(NH3)2]+
interacts with the zeolite through the H atoms of the coordinating NH3
molecules. This ”crowding” of the anionic sites by the harder NH4+ can be
observed statistically in the simulations. The average distances between the
Cu atoms and the closest framework Al atoms plotted in Figure 6d are $\sim$
4.5 Å in the L and M models (orange and red lines), and increase to $\sim$ 5.0
Å in the systems with higher Al content and thus a larger amount of charge-
balancing NH4+ (green lines). The HB model has a heterogeneous Al distribution
and its anomalously high Cu-Al distance is due to [Cu(NH3)2]+ complexes in Al-
poor regions.

Previous studies have suggested that [Cu(NH3)2]+ migration is fast only
between the three cages sharing a common framework Al, while long-range
diffusion to non-adjacent cages is limited to $\sim$ 9 Å due to decaying
electrostatic to the [Cu(NH3)2]+ and the anionic AlO4- site. 9, 17, 19. This
argument is rigorously true when the final cages contain no additional Al
sites, as our model systems from the first section. In real systems, however,
the long-range diffusion is easily explained through a sequence of local steps
combining the crossing of 8R windows into adjacent Al-containing cages,
followed by the exchange of NH4+ as compensating cation (Scheme S8 in the
Supporting Information). The low mobility of NH4+ revealed by the present
simulations suggests that a limiting factor for such long-range diffusion of
[Cu(NH3)2]+ complexes is the slower rate of counter-migration of the charge-
compensating NH4+. While [Cu(NH3)2]+ could act as a migrating compensating
cation, this results in no net migration of copper.

To explore this hypothesis, two additional MD simulations of 5 ns were run
using two modified M20 models, one of them containing 30 protons as
compensating cations, labeled M20-H+, and another one with 60 additional NH3
molecules added to the system, labeled M20-NH3. Protons on the Brønsted acid
sites are fairly static and localized within the four O atoms directly
attached to Al, so long-range charge-compensation diffusion of H+ should be
hindered. In contrast, additional NH3 should facilitate the movement of the
positive charges via proton transfer from NH4+ to NH3 via a Grotthuss-like
chain of proton transfers, thus allowing faster charge compensation between
separated AlO4- units without the need of physically displacing the strongly
attached NH4+ cations.

The plots in Figure 7 confirm this hypothesis. The MSD of Cu atoms in Figure
7a does not change when NH4+ cations are substituted by protons (black lines)
with a similar electrostatic shielding effect as NH4+. However, copper
mobility increases significantly in the presence of excess NH3 molecules
(brown lines), suggesting that free NH3 facilitates the charge re-equilibrium
following [Cu(NH3)2]+ crossing 8R. In keeping, the average distances between
the Cu atoms and the closest framework Al atoms increases from $\sim$ 4.5 Å to
$\sim$ 6.0 Å with the added NH3 (Figure 7b). Lastly, the number of distinct
cages visited by [Cu(NH3)2]+ (Figure 7c) indicates again a slightly lower
long-range mobility of [Cu(NH3)2]+ in the model containing Brønsted acid sites
and enhanced diffusion in the presence of an excess of NH3.

Our results provide theoretical backing to recent experimental observations of
the enhancing effect of gas phase NH3 during the solid-state ion-exchange of
copper using mixtures of copper oxides and zeolites, which allows the fast
preparation of Cu-exchanged zeolites at low temperatures (473-523 K) 76, 77.
After $\sim$ 4000 ps the three models reached similar steady states with
$\sim$ 3.5 distinct cages being visited by each complex, corresponding to the
final state of the ion-exchange process with a random distribution of
[Cu(NH3)2]+ complexes occupying the whole unit cell.

Figure 7: Influence of framework charge shielding on the mobility of
[Cu(NH3)2]+. (a) Mean square displacement (MSD) of Cu atoms in the [Cu(NH3)2]+
complexes, (b) average minimum Cu-Al distances and (c) number of distinct
cages visited obtained from NPP-based unbiased MD simulations at 500 K in
models with different charge-compensating cations and additional NH3
molecules. MSDs where calculated using eq 1.

###  1.5 Bimolecular complexes and mechanistic implications for the NH3-SCR-
NOx reaction

According to the proposed mechanism,9 the reaction rate depends directly on
the probability of finding simultaneously two [Cu(NH3)2]+ complexes in the
same cage, which we analyzed separately from the mobility calculations.

Figure 8a tracks the time evolution of the average number of [Cu(NH3)2]+ pairs
present simultaneously within the same cha cage. As expected, the probability
of having two [Cu(NH3)2]+ complexes in the same cage directly correlates with
the total amount of Cu in the system (compare L4 with L20 or M4 with M20
models). More interestingly, for a given Cu loading, the average number of
[Cu(NH3)2]+ pairs in the same cage increases with increasing the Al content in
the system (compare L4 with M4 or L20 with M20). This trend runs opposite to
our observation about the number of distinct cages visited, which was highest
for the model with lowest Al and lowest Cu content. That is, the chemical
environments that result in a higher rate of long-range diffusion of copper
complexes do not result in formation of more pairs, or by extension in higher
reactivity, .

As regards the role of Al distribution, the four models with a Si/Al ratio
$\sim$10 and different Al pairings (H6R, H8R, HR, HB) (see Methods Section,
Fig 5 and Table S3) behave similarly and the plots in Figure 8a only suggest a
lower probability of finding two [Cu(NH3)2]+ in the same cage in the H6R model
with Al pairs in the 6R windows, and perhaps a slightly higher probability of
bringing two [Cu(NH3)2]+ complexes together in the H8R model with Al pairs in
the 8R windows. The increased NH3 loading in M20-NH3 also results in increased
number of pairs (Figure S6).

Figure 8: Probability of formation of the active sites and catalytic results
for the NH3-SCR-NOx reaction. (a) Time evolution of the average number of
[Cu(NH3)2]+ pairs in the same cage for the eight systems in Figure 5. The
geometric center of every cage in the unit cell are determined and two Cu
atoms are considered to be in the same cage when both have the same cage
center as the nearest cage center. (b) Turnover frequency (TOF) values on a
per Cu-ion basis for the $\sim$1.5%wt Cu-containing zeolites with different
Si/Al molar ratios (b) and CHA07 (Si/Al$\sim$7) zeolite with two different Cu-
contents ($\sim$1.5 and $\sim$3 %wt Cu) (c).

###  1.6 Experimental validation for the low-T NH3-SCR-NOx reaction catalyzed
by Cu-CHA zeolites with controlled composition

To experimentally validate the computational predictions, three CHA samples
with different Si/Al molar ratios ranging from 7.3 to 23.3, which translates
to a broad range of 1.4 to 0.5 Al sites per cha (see Table S4 in the
Supporting Information) were synthesized as described in Methods. Then, the
same Cu loading ($\sim$1.5%wt Cu) was introduced within the three CHA
materials, which resulted in a similar amount of initial Cu atoms per cha
cage, $\sim$0.17, but different Cu/Al ratios (from 0.11 to 0.35, see Table
S5). In addition, the CHA sample with Si/Al ratio $\sim$7.3 was also loaded
with 3.0%wt Cu resulting in an additional sample with increased amount of
initial Cu atoms per cha cage (0.3).

The catalytic tests to evaluate the low-temperature SCR-NOx activity of the
different Cu-CHA materials were performed at very high space velocities
(1800000 ml/h.grcat) to assure low NO conversions ($<$20%). Turnover frequency
values (TOF) were obtained for each Cu-CHA sample at five different
temperatures (see Figure 8b). The TOF values obtained for the three catalysts
with different Si/Al molar ratio and the same Cu content ($\sim$1.5%wt Cu)
exhibit a continuous activity enhancement as the Al/cha cage ratio increases
from 0.5 to 1.4 (see Figure 8b), in agreement with the theoretical conclusion
that the probability of finding simultaneously two [Cu(NH3)2]+ complexes in
the same cha cage increases with increasing the Al content in the zeolite
system. Comparison of the TOF values obtained for samples with the same Si/Al
ratio and different Cu content (full and dotted green lines in Figure 8b) or
even with similar Cu/Al ratio and different Cu content (red and dotted green
line in Figure 8b) confirm that the catalytic activity clearly improves with
increasing the Cu/cha ratio, in good agreement with the theoretical conclusion
that the probability of forming [Cu(NH3)2]+ pairs in the same cha cage
directly correlates with the total amount of Cu in the system.

The experimental apparent activation energies (Figure 8c) are similar for all
samples irrespective of the Si/Al ratio, Cu content or catalytic activity,
ranging from 10.4 $\pm$ 0.8 to 11.9 $\pm$ 0.3 kcal/mol. The catalytic activity
measured by the TOF and normalized by Cu content, however, increases in
parallel with the calculated likelihood of two-copper encounters in the same
cage (Fig 8a). This supports the argument that the formation of [Cu(NH3)2]+
pairs in the same cha cage is responsible for the generation of the binuclear
active sites that catalyze the reaction, although copper diffusion may not
necessarily be the rate-determining step of the global process.

##  2 Conclusions

Biased and unbiased MD simulations using a newly trained NNP have achieved
high accuracy, chemical diversity, length- and time scale, allowing the
systematic investigation of the influence of catalyst composition and adsorbed
NH3 on the mobility of [Cu(NH3)2]+ cations in Cu-CHA catalysts.

Biased simulations on small systems showed that single [Cu(NH3)2]+ cation hops
between adjacent cha cages are very sensitive to the Al distribution and, in
general, Al pairs in 8R windows lower the free energy barrier for diffusion
and stabilize the product configuration with two [Cu(NH3)2]+ cations in the
same cage. This might be taken to suggest that the rate of the SCR-NOx
reaction could be accelerated by selectively positioning the Al atoms as Al
pairs. However, even though those results are well beyond the limits of
traditional AIMD, the simulation cells employed are overly simple and lack
realistic Al and NH4+ concentration.

Unbiased MD simulations using time scales of multiple nanosecods and
supercells with over 2300 framework atoms at a variety of Si/Al ratios, Cu+,
NH4+ and NH3 loadings show that [Cu(NH3)2]+ cations can visit on average 3-4
cages and diffuse as far as 30 Å in a few nanoseconds. They also show that the
long-range migration to remote cages requires the simultaneous displacement of
a charge-compensating NH4+ cation. An excess of NH3 facilitates the movement
of the positive charges via proton transfer from NH4+ to NH3, thus enhancing
the long-range diffusion of [Cu(NH3)2]+ complexes.

Regarding catalytic activity, we observed that the probability of finding two
[Cu(NH3)2]+ complexes in the same cage, which is necessary for the SCR-NOx
reaction at low temperature, correlates directly with the Cu content and on
the Al content but not so much on the Al distribution. These trends were
confirmed experimentally through testing the SCR-NOx reaction at low
temperatures using Cu-CHA zeolites with different Si/Al and Cu/Al molar
ratios, where we found increasing catalytic performance with increasing Al and
Cu loading.

These results demonstrate the power of combining high-throughput DFT
calculations, machine learning and molecular dynamics simulations for
simulating transport in nanoporous catalysts. The collection of the training
data and the training of the NNP had a lower total computational cost than a
single traditional AIMD simulation and resulted in scalable, fast and accurate
simulations.

This overall strategy is broadly applicable to other unsolved questions in
nanoporous catalysts, since it enables DFT-accuracy, nanosecond-long MD
simulations for thousands of atoms, and possibly beyond, by combining ML and
enhanced sampling techniques78.

##  3 Methods

###  3.1 Catalyst Models

The Cu-CHA catalytic system was modeled using three different supercells of
increasing size (see Figure S1). The smallest system is a 1x1x1 hexagonal
supercell of CHA framework containing 2Al, 34 Si and 72 O atoms, with lattice
parameters, a = b = 13.81 Å, c = 15.00 Å, $\alpha$ = 89.86°, $\beta$ = 89.94°,
$\gamma$ = 120.41°. Three models with different lattice distribution of the
two Al atoms were created, and the two negative charges introduced by the
presence of Al in the framework were compensated with two [Cu(NH3)2]+
complexes (Figure 2a). These models were used to perform DFT-based umbrella
sampling (US) molecular dynamics (MD) simulations and in the first generation
of the NNP.

A 2x2x2 triclinic supercell containing 96 T and 192 O atoms, with lattice
parameters a = 18.68 Å, b = 18.67 Å, c = 18.67 Å, $\alpha$ = 94.68°, $\beta$ =
94.63°, $\gamma$ = 94.72°was used for the subsequent training of the NN
potentials including additional Al distributions, Al contents, and
compensating species ([Cu(NH3)2]+, NH4+ and H+) as well as additional NH3
molecules as summarized in Table S1. For the systems with 3 and 7 Al atoms,
ten random Al distributions were generated. For the systems with two Al atoms
ten different Al distributions were generated, one with 2 Al in the same 4R,
one with two Al in the same 6R, four of them containing 2 Al in the same 8R,
and four of them with the 2 Al in different rings. Single-point DFT
calculations and NNP US-MD simulations were performed using this model.

Finally, a 4x4x4 triclinic supercell containing 768 T and 1536 O atoms, with
lattice parameters a = 37.35 Å, b = 37.38 Å, c = 37.34 Å, $\alpha$ = 94.64°,
$\beta$ = 94.59°, $\gamma$ = 94.47 was employed to run NNP unbiased MD
simulations on more realistic systems and reaction conditions. Three lattice
compositions corresponding to Si/Al $\sim$ 30 (26 Al, 742 Si and 1536 O
atoms), Si/Al $\sim$ 15 (50 Al, 718 Si and 1536 O atoms) and Si/Al $\sim$ 10
(50 Al, 718 Si and 1536 O atoms) were considered, and for each of them three
random Al distributions were generated. For the model with the highest Al
content (Si/Al $\sim$ 10), three additional Al distributions were generated in
which all Al atoms were either forming pairs in the 6R units, forming pairs in
the 8R windows, or heterogeneously distributed along the system so that Al-
rich and Al-poor regions are found in the same unit cell. The negative charges
in each model were compensated with [Cu(NH3)2]+, NH4+ and H+ as summarized in
Table S3 in the Supporting Information. In the models with low and medium Al
content, two Cu loadings were considered, 4 and 20 atoms per unit cell,
corresponding to Cu spatial densities of 0.08 and 0.4 Cu/1000 Å2 respectively
and 0.5% and 2.6% Cu/Si ratios (systems labeled L4, L20, M4 and M20, Figure
5). In the models with high Al content, the Cu loading was always 20 atoms per
unit cell, and the labels in Figure 5 indicate the four different Al
distributions considered, with Al pairs in the same 6R (H6R), with Al pairs in
the same 4R (H4R), random distribution (HR) and a biased, spatially-
heterogeneous distribution containing regions with high and low Al content
(HB).

All DFT calculations (single points evaluations as well as the molecular
dynamics simulations) were carried out at the revPBE+D3 79, 80 level of theory
with the software CP2K 81. The Gaussian and Plane Waves (GPW) method 82 was
used with the TZVP basis set for all atoms except Cu, which was described with
the DZVP-MOLOPT-SR-GTH basis set. A cutoff energy of 400 Ry was used for the
auxiliary plane waves and the core electrons were represented with GTH
pseudopotentials 83.

All NNP-MD simulations were performed using NNPs trained with the PaiNN52
architecture. Regular NNP-MD simulations were run in the NVT ensemble at 500
K, and consisted of a production run of at least 3 ns after 100 ps of
equilibration. The temperature was controlled by a Nosé–Hoover chain
thermostat 84, 85 with three beads and a time constant of 200 fs. The time
step to integrate the equations of motion was set to 0.5 fs.

The mean squared displacements (MSD) obtained for [Cu(NH3)2]+ and NH4+ were
calculated via the following equation:

$MSD(t)=\frac{1}{t_{s}-t+1}\sum_{t_{0}=0}^{t_{s}=t}\frac{1}{N}\sum_{i=1}^{N}\left[r_{i}(t_{0}+t)-r_{i}(t_{0})\right]^{2}$
(1)

where N corresponds to the number of [Cu(NH3)2]+ or NH4+ cations, $t_{s}$ is
the number of time steps in the simulation, the $t_{0}$ are the different time
origins.

Gibbs free energy profiles for the diffusion of [Cu(NH3)2]+ species through
the 8R windows of Cu-SSZ-13 zeolite were obtained from umbrella sampling (US)
simulations 86 at 423 K. The success of the US technique relies on the choice
of a collective variable (CV) that represents unambiguously the process
considered. It has been shown in the literature that the collective variable
($\xi$) that best describes the diffusion of [Cu(NH3)2]+ from one cavity A
($\xi<0$ ) to another one B ($\xi>0$ ) is the projection of the position
vector of the Cu atom on the vector normal to the average plane of the 8R (see
Figure S2). This collective variable takes the value $\xi=0$ when the Cu atom
is located in the plane of the 8R. In the biased simulations, the collective
variable $\xi$ was split into 40 (for DFT-based simulations) or 80 (for NN-
based simulations) equidistant windows in the range $\xi$ = -4 to $\xi$ = 4.
In each window, an independent biased MD simulation of 20 ps was run
restricting the collective variable to one specific value but ensuring
complete sampling of the configuration space in all other degrees of freedom.
A harmonic bias potential centered at the corresponding value of the
collective variable was applied to restrict the sampling to each window
individually and to ensure sufficient overlap between the sampling of adjacent
windows. The force constant of the harmonic potential was set to 12
kcal/mol/Å2. All DFT-based biased simulations were performed with PLUMED87
interfaced to the CP2K engine.

###  3.2 Neural network potential training

All NN potentials trained were based on the PaiNN architecture 52 which uses
equivariant message-passing for the ground truth prediction. For the training
we used three convolutions, a cutoff radius of 5Å, 32 Gaussian functions for
distance expansion and 128-size vector for the atomic features. Neural network
parameters were optimized using Adam algorithm with a batch size of 32,
initial learning rate of $1\text{\times}{10}^{-4}$, learning decay of 0.5 and
learning patience of 25.

The reference energies of Si, Al, Cu, N, O, and H were evaluated using the
revPBE+D3 energies of the lowest energy structure of each composition. These
energies were used to fit a linear regression model with atomic contribution
to the reference energy of each element as predictor variables. The reference
energy was subtracted from the revPBE+D3 energy of every geometry in the
training, test and validation sets. The forces were used without
modifications. The training coefficients for the energies and the forces were
set to 0.01 and 1.0, respectively.

The acquisition of training data was performed using active learning (AL) with
a query-by-committee approach61, 62, 63, 64, 65, 66, 67, 68 as illustrated in
Figure 1a. In this approach a committee (ensemble) of NNPs is trained on the
available labelled data and new data is collected based on the maximum
disagreement (variance) of the prediction the committee members.

The first generation of the potential was trained on a randomly collected
subset of the DFT data generated from a previous study 17 and from three
biased simulations performed with DFT at 423 K, used as reference ground
truth. In total, there were $\sim$9000 geometries in the initial dataset. This
pretrained potential was then retrained in 4 active learning loops using the
2x2x2 triclinic supercell described in previous section (see Figure S1. For
each loop, biased MD trajectories were generated with the learned interatomic
potential of the previous loop. The following temperatures were used in the
active learning loops 298 K, 423 K, 500 K and 550 K. The bias was applied
using a harmonic potential on the collective variable used in US. Each MD was
run for 100 ps and geometries were collected every 100 fs. The selection of
the new geometries from the MD trajectories was carried out using as criterion
the force uncertainty from an ensemble of three NN potentials. The variances
of the forces from the ensemble of potentials were ranked in descending order
and the first geometries were selected to increase the dataset in $\sim$10% .
The nonphysical geometries and those with low uncertainty, $<$ 2 kcal/mol were
discarded. Up to this point, the dataset contained structures where the Al
substitutions where compensated with Cu+ as diaminecopper(I) so that the
trained potential did not properly describe local environments of Al
compensated with NH4+ or H+. The acquisition of new geometries with new
compositions including NH4+ and H+ was performed using adversarial attack 71
for 6 more generations with NNP trained on the last generation of AL. The
initialization of adversarial attacks were performed by displacing each
coordinate $\sigma$ $\sim$N(0,0.01 Å) of one optimized geometry for every
chemical composition containing NH4+ and H+. The resulting attack $\sigma$ was
optimized for 100 iterations using the Adam optimizer with learning rate of
0.0001 and with the normalized temperature kT set to 20 kcal/mol.

Then, five more generations of active learning was performed, with biased MD
simulations at temperatures ranging from 600 to 1000 K. All models were
trained by randomly partitioning the available data into training (60%),
validation (20%) and test sets (20%). This splitting was perform for each
composition and then the training, test and validation sets of each
composition were combined in the final training, test and validation sets.
Thus, the final training, test and validation sets had a balanced distribution
of all chemical compositions.

###  3.3 Zeolite synthesis

####  3.3.1 CHA zeolite with a Si/Al$\sim$7 (CHA07)

5.91 g of a 25 wt% aqueous solution of
N,N,N-trimethyl-1-adamantylammoniumhydroxide (TMAda, Sachem) was added to 5.41
g of deionized H2O. Next, 0.11 g of Al(OH)3 (76.5 wt%, Thermo Fisher) and 0.85
g of a 5M sodium hydroxide solution (NaOH: 16.7 wt% NaOH in deionized water;
NaOH pellets 98 wt%, Alfa Aesar) was added to the aqueous TMAdaOH solution and
the mixture was stirred under ambient conditions for 15 minutes. Finally, 2.13
g of colloidal silica (Ludox HS40, 40 wt%, Sigma Aldrich) was added to the
mixture and stirred for 2 h under ambient conditions. The final gel
composition was: SiO2 : 0.056 Al2O3 : 0.49 TMAdaOH : 0.25 NaOH : 46.5 H2O. The
resultant gel was charged into a stainless steel autoclave with a Teflon
liner. The crystallization was then conducted at 160 °C for 6 days under
dynamic conditions. The solid product was filtered, washed with abundant
water, and dried at 100 °C. The solids were calcined at 580 °C for 5 h in air.

####  3.3.2 CHA zeolite with a Si/Al$\sim$13 (CHA13)

5.89 g of a 25 wt% aqueous solution of
N,N,N-trimethyl-1-adamantylammoniumhydroxide (TMAda, Sachem) was added to 5.42
g of deionized H2O. Next, 0.042 g of Al(OH)3 (76.5 wt%, Thermo Fisher) and
0.86 g of a 5M sodium hydroxide solution (NaOH: 16.7 wt% NaOH in deionized
water; NaOH pellets 98 wt%, Alfa Aesar) was added to the aqueous TMAdaOH
solution and the mixture was stirred under ambient conditions for 15 minutes.
Finally, 2.10 g of colloidal silica (Ludox HS40, 40 wt%, Sigma Aldrich) was
added to the mixture and stirred for 2 h under ambient conditions. The final
gel composition was: SiO2 : 0.02 Al2O3 : 0.50 TMAdaOH : 0.26 NaOH : 47H2O. The
resultant gel was charged into a stainless steel autoclave with a Teflon
liner. The crystallization was then conducted at 160 °C for 6 days under
dynamic conditions. The solid product was filtered, washed with abundant
water, and dried at 100 °C. The solids were calcined at 580 °C for 5 h in air.

####  3.3.3 CHA zeolite with a Si/Al$\sim$23 (CHA23)

1.55 g of FAU zeolite (FAU, CBV760 with Si/Al=26, Zeolyst, lot number:
76004N002648) was added to 6.74 g of a 25 wt% aqueous solution of
N,N,N-trimethyl-1-adamantylammoniumhydroxide (TMAda, Sachem). The mixture was
maintained under stirring the required time to evaporate the excess of water
until achieving the desired gel concentration. The final gel composition was:
SiO2 : 0.019 Al2O3 : 0.4 TMadaOH : 5 H2O. The resultant gel was charged into a
stainless steel autoclave with a Teflon liner. The crystallization was then
conducted at 175 °C for 12 days under static conditions. The solid product was
filtered, washed with abundant water, and dried at 100 °C. The solids were
calcined at 580 °C for 5 h in air.

####  3.3.4 Cu-exchange treatments

The Na-containing CHA calcined solids were first exchanged with a 2M aqueous
solution of ammonium nitrate (NH4Cl, Sigma-Aldrich, 99% by weight) with a
liquid/solid ratio of 10, maintaining the mixture at 80 °C for 2 hours under
agitation. Afterwards, the solids were recovered by filtration. 0.3 g of the
Na-free zeolites was introduced in 30 ml of an aqueous solution of
Cu(CH3COO)2·H2O [14.1 mg of Cu(CH3COO)2·H2O dissolved in 30 ml of water or
28.3 mg of Cu(CH3COO)2·H2O dissolved in 30 ml of water for materials
containing $\sim$1.5 or 3%wt Cu, respectively], maintaining a solid/liquid
ratio of 10 g/l at 80 °C for 24 h. Finally, the solids were filtered and
washed with distilled water, dried and calcined at 550 °C in air for 4 h.

###  3.4 Characterization

Powder X-ray diffraction (PXRD) measurements were performed with a multi
sample Philips X’Pert diffractometer equipped with a graphite monochromator,
operating at 40 kV and 35 mA, and using Cu K$\alpha$ radiation ($\lambda$ =
0.1542 nm). The PXRD patterns reveal the good crystallization of the CHA
materials (see Figure S7 in the Supporting Information), all of them
presenting particle sizes within the sub-micron scale (below 1 µm, see Figure
S7 and Table S4 in the Supporting Information). Chemical analyses were carried
out in a Varian 715-ES ICP-Optical Emission spectrometer, after solid
dissolution in HNO3/HCl/HF aqueous solution. Nitrogen adsorption isotherms at
-196 °C were measured on a Micromeritics ASAP 2020 with a manometric
adsorption analyser to determinate the textural properties of the samples. The
morphology of the samples was studied by field emission scanning electron
microscopy (FESEM) using a ZEISS Ultra-55 microscope.

###  3.5 Catalytic evaluation

The catalytic activity was evaluated for the selective catalytic reduction
(SCR) of NOx with NH3 in a fixed bed, quartz tubular reactor with $\sim$1 cm
inner diameter. 20 mg sieve fractionated catalysts (200-400 µm) were placed
between 0.6 of silicon carbide fractions (SiC, Fisher Chemical, 200-400 µm)
generating an overall bed height of $\sim$ 1.2 cm. The catalysts were
introduced in the reactor and heated up to 550 °C in a 300 mL/min flow of
nitrogen and maintained at this temperature for one hour. Afterwards, the feed
was admitted over the catalyst with an overall flow of 600 mL/min. The feed
composition for the catalytic tests performed over the Cu-containing catalysts
was 500 ppm NO, 550 ppm NH3, 7% O2 and 5% H2O, resulting in a very high space
velocity (GHSV=1800000 ml/h.grcat). The evaluated reaction temperatures were
190, 180 and 170 °C. The conversion of NO was measured under steady state
conversion at each temperature using a chemiluminescence detector (Thermo
62C). The TOF were estimated by dividing the moles of NO molecules converted
per second by the moles of Cu atoms in the catalysts, according to

$TOF=\frac{moles\ NO\ converted/s}{moles\ Cu}$ (2)

From the low temperature NO conversion results, the rate constants (k) can be
calculated using a first-order kinetic equation, as described previously in
the literature: 88

$k=-\frac{F_{0}}{[NO]_{0}W}\ln(1-X)$ (3)

where $F_{0}$ is the molar NO feed rate, $[NO]_{0}$ is the molar concentration
at the inlet, W is the catalyst amount (gr) and X is the NO conversion. The
Arrhenius equation was employed to estimate the apparent activation energies
(Ea) after its linearization as follows:

$k=Ae^{-E_{a}/RT}$ (4) $\ln{k}=\ln{A}-\frac{E_{a}}{R}\left(\frac{1}{T}\right)$
(5)

where A is the pre-exponential factor, R is the universal gas constant and T
denotes the absolute temperature associated with the reaction (in Kelvin).

##  4 Data availability

The experimental and computational data that support the findings of this
study are available from the corresponding author upon reasonable request. The
datasets generated during this study are available at
https://figshare.com/projects/Dataset_and_machine_learning_potential_Cu-
CHA/167645. The code used for this study can be downloaded from
https://github.com/learningmatter-mit/NeuralForceField.

## References

  * Sushkevich et al. 2017 Sushkevich, V. L.; Palagin, D.; Ranocchiari, M.; Van Bokhoven, J. A. Selective anaerobic oxidation of methane enables direct synthesis of methanol. _Science_ 2017, _356_ , 523–527
  * Dinh et al. 2019 Dinh, K. T.; Sullivan, M. M.; Narsimhan, K.; Serna, P.; Meyer, R. J.; Dincǎ, M.; Román-Leshkov, Y. Continuous Partial Oxidation of Methane to Methanol Catalyzed by Diffusion-Paired Copper Dimers in Copper-Exchanged Zeolites. _Journal of the American Chemical Society_ 2019, _141_ , 11641–11650
  * Dinh et al. 2021 Dinh, K. T.; Sullivan, M. M.; Serna, P.; Meyer, R. J.; Román-Leshkov, Y. Breaking the Selectivity-Conversion Limit of Partial Methane Oxidation with Tandem Heterogeneous Catalysts. _ACS Catalysis_ 2021, _11_ , 9262–9270
  * Del Campo et al. 2021 Del Campo, P.; Martínez, C.; Corma, A. Activation and conversion of alkanes in the confined space of zeolite-type materials. _Chemical Society Reviews_ 2021, _50_ , 8511–8595
  * Ohyama et al. 2022 Ohyama, J.; Tsuchimura, Y.; Hirayama, A.; Iwai, H.; Yoshida, H.; Machida, M.; Nishimura, S.; Kato, K.; Takahashi, K. Relationships among the Catalytic Performance, Redox Activity, and Structure of Cu-CHA Catalysts for the Direct Oxidation of Methane to Methanol Investigated Using In Situ XAFS and UV-Vis Spectroscopies. _ACS Catalysis_ 2022, _12_ , 2454–2462
  * Kwak et al. 2010 Kwak, J. H.; Tonkyn, R. G.; Kim, D. H.; Szanyi, J.; Peden, C. H. F. Excellent activity and selectivity of Cu-SSZ-13 in the selective catalytic reduction of NOx with NH3. _Journal of Catalysis_ 2010, _275_ , 187–190
  * Paolucci et al. 2016 Paolucci, C.; Parekh, A. A.; Khurana, I.; Di Iorio, J. R.; Li, H.; Albarracin Caballero, J. D.; Shih, A. J.; Anggara, T.; Delgass, W. N.; Miller, J. T.; Ribeiro, F. H.; Gounder, R.; Schneider, W. F. Catalysis in a cage: Condition-dependent speciation and dynamics of exchanged cu cations in ssz-13 zeolites. _Journal of the American Chemical Society_ 2016, _138_ , 6028–6048
  * Gao et al. 2017 Gao, F.; Mei, D.; Wang, Y.; Szanyi, J.; Peden, C. H. Selective Catalytic Reduction over Cu/SSZ-13: Linking Homo- and Heterogeneous Catalysis. _Journal of the American Chemical Society_ 2017, _139_ , 4935–4942
  * Paolucci et al. 2017 Paolucci, C.; Khurana, I.; Parekh, A. A.; Li, S.; Shih, A. J.; Li, H.; Di Iorio, J. R.; Albarracin-Caballero, J. D.; Yezerets, A.; Miller, J. T.; Delgass, W. N.; Ribeiro, F. H.; Schneider, W. F.; Gounder, R. Dynamic multinuclear sites formed by mobilized copper ions in NO x selective catalytic reduction. _Science_ 2017, _357_ , 898–903
  * Peden 2019 Peden, C. H. Cu/Chabazite catalysts for ‘Lean-Burn’ vehicle emission control. _Journal of Catalysis_ 2019, _373_ , 384–389
  * Vennestrøm et al. 2022 Vennestrøm, P. N.; Thøgersen, J. R.; Gabrielsson, P. L.; Van Tendeloo, L.; Schütze, F.-W.; Moliner, M. Advances and perspectives from a decade of collaborative efforts on zeolites for selective catalytic reduction of NOx. _Microporous and Mesoporous Materials_ 2022, 112336
  * Moreno-González et al. 2019 Moreno-González, M.; Millán, R.; Concepción, P.; Blasco, T.; Boronat, M. Spectroscopic Evidence and Density Functional Theory (DFT) Analysis of Low-Temperature Oxidation of Cu + to Cu 2+ NO x in Cu-CHA Catalysts: Implications for the SCR-NO x Reaction Mechanism. _ACS Catalysis_ 2019, _9_ , 2725–2738
  * Chen et al. 2020 Chen, L.; Janssens, T. V.; Vennestrøm, P. N.; Jansson, J.; Skoglundh, M.; Grönbeck, H. A Complete Multisite Reaction Mechanism for Low-Temperature NH3-SCR over Cu-CHA. _ACS Catalysis_ 2020, _10_ , 5646–5656
  * Borfecchia et al. 2018 Borfecchia, E.; Beato, P.; Svelle, S.; Olsbye, U.; Lamberti, C.; Bordiga, S. Cu-CHA – a model system for applied selective redox catalysis. _Chemical Society Reviews_ 2018, _47_ , 8097–8133
  * Signorile et al. 2022 Signorile, M.; Borfecchia, E.; Bordiga, S.; Berlier, G. Influence of ion mobility on the redox and catalytic properties of Cu ions in zeolites. _Chemical Science_ 2022, _13_ , 10238–10250
  * Millan et al. 2020 Millan, R.; Cnudde, P.; Hoffman, A. E.; Lopes, C. W.; Concepción, P.; Van Speybroeck, V.; Boronat, M. Theoretical and Spectroscopic Evidence of the Dynamic Nature of Copper Active Sites in Cu-CHA Catalysts under Selective Catalytic Reduction (NH3-SCR-NOx) Conditions. _Journal of Physical Chemistry Letters_ 2020, _11_ , 10060–10066
  * Millan et al. 2021 Millan, R.; Cnudde, P.; Speybroeck, V. v.; Boronat, M. Mobility and Reactivity of Cu+ Species in Cu-CHA Catalysts under NH3-SCR-NOx Reaction Conditions: Insights from AIMD Simulations. _JACS Au_ 2021, _1_ , 1778–1787
  * Martini et al. 2022 Martini, A.; Negri, C.; Bugarin, L.; Deplano, G.; Abasabadi, R. K.; Lomachenko, K. A.; Janssens, T. V.; Bordiga, S.; Berlier, G.; Borfecchia, E. Assessing the Influence of Zeolite Composition on Oxygen-Bridged Diamino Dicopper(II) Complexes in Cu-CHA DeNOxCatalysts by Machine Learning-Assisted X-ray Absorption Spectroscopy. _Journal of Physical Chemistry Letters_ 2022, _13_ , 6164–6170
  * Krishna et al. 2023 Krishna, S. H.; Goswami, A.; Wang, Y.; Jones, C. B.; Dean, D. P.; Miller, J. T.; Schneider, W. F.; Gounder, R.; Davidson, C. D. Influence of framework Al density in chabazite zeolites on copper ion mobility and reactivity during NOx selective catalytic reduction with NH3. _Nature Catalysis 2023 6:3_ 2023, _6_ , 276–285
  * Marberger et al. 2018 Marberger, A.; Petrov, A. W.; Steiger, P.; Elsener, M.; Kröcher, O.; Nachtegaal, M.; Ferri, D. Time-resolved copper speciation during selective catalytic reduction of NO on Cu-SSZ-13. _Nature Catalysis 2018 1:3_ 2018, _1_ , 221–227
  * Becher et al. 2020 Becher, J.; Sanchez, D. F.; Doronkin, D. E.; Zengel, D.; Meira, D. M.; Pascarelli, S.; Grunwaldt, J. D.; Sheppard, T. L. Chemical gradients in automotive Cu-SSZ-13 catalysts for NOx removal revealed by operando X-ray spectrotomography. _Nature Catalysis 2020 4:1_ 2020, _4_ , 46–53
  * Wu et al. 2023 Wu, Y.; Zhao, W.; Ahn, S. H.; Wang, Y.; Walter, E. D.; Chen, Y.; Derewinski, M. A.; Washton, N. M.; Rappé, K. G.; Wang, Y.; Mei, D.; Hong, S. B.; Gao, F. Interplay between copper redox and transfer and support acidity and topology in low temperature NH3-SCR. _Nature Communications_ 2023, _14_ , 1–13
  * Axelrod et al. 2022 Axelrod, S.; Schwalbe-Koda, D.; Mohapatra, S.; Damewood, J.; Greenman, K. P.; Gómez-Bombarelli, R. Learning Matter: Materials Design with Machine Learning and Atomistic Simulations. _Accounts of Materials Research_ 2022, _3_ , 343–357
  * Schmidt et al. 2019 Schmidt, J.; Marques, M. R. G.; Botti, S.; Marques, M. A. L. Recent advances and applications of machine learning in solid-state materials science. _npj Computational Materials_ 2019, _5_ , 83
  * Ma and Liu 2020 Ma, S.; Liu, Z.-P. Machine Learning for Atomic Simulation and Activity Prediction in Heterogeneous Catalysis: Current Status and Future. _ACS Catalysis_ 2020, 13213–13226
  * Goldsmith et al. 2018 Goldsmith, B. R.; Esterhuizen, J.; Liu, J.-X. X. J.; Bartel, C. J.; Sutton, C. Machine learning for heterogeneous catalyst design and discovery. _AIChE Journal_ 2018, _64_ , 2311–2323
  * Schlexer Lamoureux et al. 2019 Schlexer Lamoureux, P.; Winther, K. T.; Garrido Torres, J. A.; Streibel, V.; Zhao, M.; Bajdich, M.; Abild-Pedersen, F.; Bligaard, T. Machine Learning for Computational Heterogeneous Catalysis. _ChemCatChem_ 2019, _11_ , 3581–3601
  * Gu et al. 2020 Gu, G. H.; Noh, J.; Kim, S.; Back, S.; Ulissi, Z.; Jung, Y. Practical Deep-Learning Representation for Fast Heterogeneous Catalyst Screening. _The Journal of Physical Chemistry Letters_ 2020, _11_ , 3185–3191
  * Behler 2011 Behler, J. Neural network potential-energy surfaces in chemistry: a tool for large-scale simulations. _Physical Chemistry Chemical Physics_ 2011, _13_ , 17930
  * Behler 2017 Behler, J. First Principles Neural Network Potentials for Reactive Simulations of Large Molecular and Condensed Systems. _Angewandte Chemie International Edition_ 2017, _56_ , 12828–12840
  * Botu et al. 2017 Botu, V.; Batra, R.; Chapman, J.; Ramprasad, R. Machine learning force fields: Construction, validation, and outlook. _Journal of Physical Chemistry C_ 2017, _121_ , 511–522
  * Mueller et al. 2020 Mueller, T.; Hernandez, A.; Wang, C. Machine learning for interatomic potential models. _The Journal of Chemical Physics_ 2020, _152_ , 50902
  * von Lilienfeld et al. 2020 von Lilienfeld, O. A.; Müller, K.-R. R.; Tkatchenko, A. Exploring chemical compound space with quantum-based machine learning. _Nature Reviews Chemistry_ 2020, _4_ , 347–358
  * Noé et al. 2019 Noé, F.; Olsson, S.; Köhler, J.; Wu, H. Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning. _Science_ 2019, _365_ , eaaw1147
  * Jensen et al. 2019 Jensen, K. F.; Coley, C. W.; Eyke, N. S. Autonomous discovery in the chemical sciences part I: Progress. _Angewandte Chemie International Edition_ 2019, _59_ , 2858–22893
  * Butler et al. 2018 Butler, K. T.; Davies, D. W.; Cartwright, H.; Isayev, O.; Walsh, A. Machine learning for molecular and materials science. _Nature_ 2018, _559_ , 547–555
  * Erlebach et al. 2022 Erlebach, A.; Nachtigall, P.; Grajciar, L. Accurate large-scale simulations of siliceous zeolites by neural network potentials. _npj Computational Materials 2022 8:1_ 2022, _8_ , 1–12
  * Moliner et al. 2019 Moliner, M.; Román-Leshkov, Y.; Corma, A. Machine Learning Applied to Zeolite Synthesis: The Missing Link for Realizing High-Throughput Discovery. _Accounts of Chemical Research_ 2019, _52_ , 2971–2980
  * Chmiela et al. 2018 Chmiela, S.; Sauceda, H. E.; Müller, K.-R.; Tkatchenko, A. Towards exact molecular dynamics simulations with machine-learned force fields. _Nature Communications_ 2018, _9_ , 3887
  * Sosso et al. 2012 Sosso, G. C.; Miceli, G.; Caravati, S.; Behler, J.; Bernasconi, M. Neural network interatomic potential for the phase change material GeTe. _Physical Review B_ 2012, _85_ , 174103
  * Bernstein et al. 2019 Bernstein, N.; Bhattarai, B.; Csµnyi, G.; Drabold, D. A.; Elliott, S. R.; Deringer, V. L.; Bernstein, N.; Bhattarai, B.; Drabold, R. D. A.; Csµnyi, G.; Deringer, V. L.; Elliott, S. R. Quantifying Chemical Structure and Machine-Learned Atomic Energies in Amorphous and Liquid Silicon. _Angewandte Chemie International Edition_ 2019, _58_ , 7057–7061
  * Smith et al. 2017 Smith, J. S.; Isayev, O.; Roitberg, A. E. ANI-1, A data set of 20 million calculated off-equilibrium conformations for organic molecules. _Scientific Data_ 2017, _4_ , 170193
  * Smith et al. 2017 Smith, J. S.; Isayev, O.; Roitberg, A. E. ANI-1: an extensible neural network potential with DFT accuracy at force field computational cost. _Chemical Science_ 2017, _8_ , 3192–3203
  * Smith et al. 2019 Smith, J. S.; Nebgen, B. T.; Zubatyuk, R.; Lubbers, N.; Devereux, C.; Barros, K.; Tretiak, S.; Isayev, O.; Roitberg, A. E. Approaching coupled cluster accuracy with a general-purpose neural network potential through transfer learning. _Nature Communications_ 2019, _10_ , 1–8
  * Devereux et al. 2020 Devereux, C.; Smith, J. S.; Huddleston, K. K.; Barros, K.; Zubatyuk, R.; Isayev, O.; Roitberg, A. E. Extending the Applicability of the ANI Deep Learning Molecular Potential to Sulfur and Halogens. _Journal of Chemical Theory and Computation_ 2020, _16_ , 4192–4202
  * Schütt et al. 2017 Schütt, K. T.; Arbabzadah, F.; Chmiela, S.; Müller, K. R.; Tkatchenko, A. Quantum-chemical insights from deep tensor neural networks. _Nature Communications_ 2017, _8_ , 6–13
  * Schütt et al. 2019 Schütt, K. T.; Kessel, P.; Gastegger, M.; Nicoli, K. A.; Tkatchenko, A.; Müller, K. R. SchNetPack: A Deep Learning Toolbox for Atomistic Systems. _Journal of Chemical Theory and Computation_ 2019, _15_ , 448–455
  * Zhang et al. 2018 Zhang, L.; Han, J.; Wang, H.; Car, R.; Weinan, E. Deep Potential Molecular Dynamics: A Scalable Model with the Accuracy of Quantum Mechanics. _Physical Review Letters_ 2018, _120_ , 143001
  * Chen et al. 2019 Chen, C.; Ye, W.; Zuo, Y.; Zheng, C.; Ong, S. P. Graph Networks as a Universal Machine Learning Framework for Molecules and Crystals. _Chemistry of Materials_ 2019, _31_ , 3564–3572
  * Klicpera et al. 2020 Klicpera, J.; Groß, J.; Günnemann, S. Directional message passing for molecular graphs. ICLR. 2020; pp 1–13
  * Qiao et al. 2020 Qiao, Z.; Welborn, M.; Anandkumar, A.; Manby, F. R.; Miller, T. F. OrbNet: Deep learning for quantum chemistry using symmetry-adapted atomic-orbital features. _The Journal of Chemical Physics_ 2020, _153_ , 124111
  * Schütt et al. 2021 Schütt, K. T.; Unke, O. T.; Gastegger, M. Equivariant message passing for the prediction of tensorial properties and molecular spectra. _arXiv preprint arXiv:2102.03150_ 2021,
  * Batzner et al. 2021 Batzner, S.; Musaelian, A.; Sun, L.; Geiger, M.; Mailoa, J. P.; Kornbluth, M.; Molinari, N.; Smidt, T. E.; Kozinsky, B. E(3)-Equivariant Graph Neural Networks for Data-Efficient and Accurate Interatomic Potentials. _Nature Communications_ 2021, _13_ , 1–11
  * Artrith and Urban 2016 Artrith, N.; Urban, A. An implementation of artificial neural-network potentials for atomistic materials simulations: Performance for TiO2. _Computational Materials Science_ 2016, _114_ , 135–150
  * Vandenhaute et al. 2023 Vandenhaute, S.; Cools-Ceuppens, M.; DeKeyser, S.; Verstraelen, T.; Van Speybroeck, V. Machine learning potentials for metal-organic frameworks using an incremental learning approach. _npj Computational Materials_ 2023, _9_ , 1–8
  * Bocus et al. 2023 Bocus, M.; Goeminne, R.; Lamaire, A.; Cools-Ceuppens, M.; Verstraelen, T.; Van Speybroeck, V. Nuclear quantum effects on zeolite proton hopping kinetics explored with machine learning potentials and path integral molecular dynamics. _Nature Communications_ 2023, _14_ , 1–11
  * Winter and Gómez-Bombarelli 2022 Winter, G.; Gómez-Bombarelli, R. Simulations with machine learning potentials identify the ion conduction mechanism mediating non-Arrhenius behavior in LGPS. 2022; preprint, http://arxiv.org/abs/2211.05713v2
  * Gastegger and Marquetand 2015 Gastegger, M.; Marquetand, P. High-Dimensional Neural Network Potentials for Organic Reactions and an Improved Training Algorithm. _Journal of Chemical Theory and Computation_ 2015, _11_ , 2187–2198
  * Ang et al. 2021 Ang, S. J.; Wang, W.; Schwalbe-Koda, D.; Axelrod, S.; Gómez-Bombarelli, R. Active learning accelerates ab initio molecular dynamics on reactive energy surfaces. _Chem_ 2021, _7_ , 1–32
  * Du et al. 2023 Du, X.; Damewood, J. K.; Lunger, J. R.; Millan, R.; Yildiz, B.; Li, L.; Gómez-Bombarelli, R. Machine-learning-accelerated simulations enable heuristic-free surface reconstruction. 2023; preprint, https://arxiv.org/abs/2305.07251
  * Behler 2015 Behler, J. Constructing high-dimensional neural network potentials: A tutorial review. _International Journal of Quantum Chemistry_ 2015, _115_ , 1032–1050
  * Smith et al. 2018 Smith, J. S.; Nebgen, B.; Lubbers, N.; Isayev, O.; Roitberg, A. E. Less is more: Sampling chemical space with active learning. _Journal of Chemical Physics_ 2018, _148_ , 241733
  * Schran et al. 2020 Schran, C.; Brezina, K.; Marsalek, O. Committee neural network potentials control generalization errors and enable active learning. _Journal of Chemical Physics_ 2020, _153_ , 104105
  * Musil et al. 2019 Musil, F.; Willatt, M. J.; Langovoy, M. A.; Ceriotti, M. Fast and Accurate Uncertainty Estimation in Chemical Machine Learning. _Journal of Chemical Theory and Computation_ 2019, _15_ , 906–915
  * Peterson et al. 2017 Peterson, A. A.; Christensen, R.; Khorshidi, A. Addressing uncertainty in atomistic machine learning. _Physical Chemistry Chemical Physics_ 2017, _19_ , 10978–10985
  * Lookman et al. 2019 Lookman, T.; Balachandran, P. V.; Xue, D.; Yuan, R. Active learning in materials science with emphasis on adaptive sampling using uncertainties for targeted design. _npj Computational Materials_ 2019, _5_ , 21
  * Shapeev et al. 2020 Shapeev, A.; Gubaev, K.; Tsymbalov, E.; Podryabinkin, E. In _Machine Learning Meets Quantum Physics_ ; Schütt, K. T., Chmiela, S., von Lilienfeld, O. A., Tkatchenko, A., Tsuda, K., Müller, K.-R., Eds.; Springer International Publishing: Cham, 2020; pp 309–329
  * Imbalzano et al. 2021 Imbalzano, G.; Zhuang, Y.; Kapil, V.; Rossi, K.; Engel, E. A.; Grasselli, F.; Ceriotti, M. Uncertainty estimation for molecular dynamics and sampling. _The Journal of Chemical Physics_ 2021, _154_ , 74102
  * Wang et al. 2020 Wang, W.; Yang, T.; Harris, W. H.; Gómez-Bombarelli, R. Active learning and neural network potentials accelerate molecular screening of ether-based solvate ionic liquids. _Chemical Communications_ 2020, _56_ , 8920–8923
  * Wang et al. 2020 Wang, W.; Axelrod, S.; Gómez-Bombarelli, R. Differentiable Molecular Simulations for Control and Learning. 2020; preprint, http://arxiv.org/abs/2003.00868
  * Schwalbe-Koda et al. 2021 Schwalbe-Koda, D.; Tan, A. R.; Gómez-Bombarelli, R. Differentiable sampling of molecular geometries with uncertainty-based adversarial attacks. _Nature Communications_ 2021, _12_ , 1–12
  * Tan et al. 2023 Tan, A. R.; Urata, S.; Goldman, S.; Dietschreit, J. C. B.; Gómez-Bombarelli, R. Single-model uncertainty quantification in neural network potentials does not consistently outperform model ensembles. 2023; preprint, https://arxiv.org/abs/2305.01754v1
  * McInnes et al. 2018 McInnes, L.; Healy, J.; Saul, N.; Großberger, L. UMAP: Uniform Manifold Approximation and Projection. _Journal of Open Source Software_ 2018, _3_ , 861
  * 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
  * Fu et al. 2022 Fu, X.; Wu, Z.; Wang, W.; Xie, T.; Keten, S.; Gomez-Bombarelli, R.; Jaakkola, T. Forces are not Enough: Benchmark and Critical Evaluation for Machine Learning Force Fields with Molecular Simulations. 2022; preprint, https://arxiv.org/abs/2210.07237v1
  * Shwan et al. 2015 Shwan, S.; Skoglundh, M.; Lundegaard, L. F.; Tiruvalam, R. R.; Janssens, T. V.; Carlsson, A.; Vennestrøm, P. N. Solid-State Ion-Exchange of Copper into Zeolites Facilitated by Ammonia at Low Temperature. _ACS Catalysis_ 2015, _5_ , 16–19
  * Vennestrøm et al. 2019 Vennestrøm, P. N.; Lundegaard, L. F.; Tyrsted, C.; Bokarev, D. A.; Mytareva, A. I.; Baeva, G. N.; Stakheev, A. Y.; Janssens, T. V. The Role of Protons and Formation Cu(NH 3 ) 2+ During Ammonia-Assisted Solid-State Ion Exchange of Copper(I) Oxide into Zeolites. _Topics in Catalysis_ 2019, _62_ , 100–107
  * Šípka et al. 2023 Šípka, M.; Dietschreit, J. C. B.; Grajciar, L.; Gómez-Bombarelli, R. Differentiable Simulations for Enhanced Sampling of Rare Events. 2023,
  * Zhang and Yang 1998 Zhang, Y.; Yang, W. Comment on “Generalized Gradient Approximation Made Simple”. _Physical Review Letters_ 1998, _80_ , 890
  * Grimme et al. 2010 Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. _The Journal of Chemical Physics_ 2010, _132_ , 154104
  * Vandevondele et al. 2005 Vandevondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. _Computer Physics Communications_ 2005, _167_ , 103–128
  * Lippert et al. 1999 Lippert, G.; Hutter, J.; Parrinello, M. The Gaussian and augmented-plane-wave density functional method for ab initio molecular dynamics simulations. _Theoretical Chemistry Accounts 1999 103:2_ 1999, _103_ , 124–140
  * Goedecker and Teter 1996 Goedecker, S.; Teter, M. Separable dual-space Gaussian pseudopotentials. _Physical Review B_ 1996, _54_ , 1703
  * Martyna et al. 1992 Martyna, G. J.; Klein, M. L.; Tuckerman, M. Nosé-Hoover chains: The canonical ensemble via continuous dynamics. _The Journal of Chemical Physics_ 1992, _97_ , 2635–2643
  * Hoover et al. 1982 Hoover, W. G.; Ladd, A. J. C.; Moran, B. High-Strain-Rate Plastic Flow Studied via Nonequilibrium Molecular Dynamics. _Physical Review Letters_ 1982, _48_ , 1818–1820
  * Torrie and Valleau 1977 Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. _Journal of Computational Physics_ 1977, _23_ , 187–199
  * Kumar et al. 1992 Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. THE weighted histogram analysis method for free???energy calculations on biomolecules. I. The method. _Journal of Computational Chemistry_ 1992, _13_ , 1011–1021
  * Long and Yang 2002 Long, R. Q.; Yang, R. T. Selective Catalytic Reduction of NO with Ammonia over Fe3+-Exchanged Mordenite (Fe–MOR): Catalytic Performance, Characterization, and Mechanistic Study. _Journal of Catalysis_ 2002, _207_ , 274–285

R.M acknowledges the Margarita Salas grant from the Ministerio de
Universidades, Spain, funded by the European Union-Next Generation EU. The
authors are grateful for computation time allocated on the MIT SuperCloud
cluster, the MIT Engaging cluster at the Massachusetts Green High Performance
Computing Center (MGHPCC); and Summit at the Oakridge Leadership Computing
Facility through the 2021 ALCC DOE program. R.G-B. thanks the Jeffrey Cheah
Career Development Chair. R.M thanks Gavin Winter for assistance during the
training of NNP and processing of the MD simulations, and Simon Axelrod for
implementing PaiNN. M.B. and M.M. thank financial support from Spanish
Government through CEX2021-001230-S, PID2020-112590GB-C21,
PID2021-122755OB-I00 and TED2021-130739B-I00 (MCIN/AEI/FEDER, UE), and from
CSIC through the I-link+ Program (LINKA20381). E.B. acknowledges the Spanish
Government-MCIU for a FPI scholarship (PRE2019-088360). The Electron
Microscopy Service of the UPV is acknowledged for its help in sample
characterization.

##  5 Author Contributions

R.G-B. and M. B. designed the project. R. M. developed the NNP and performed
all simulations, supervised by R.G.-B. M. M. directed the zeolite synthesis
and catalytic experiments. E. B. performed the zeolite synthesis and catalytic
experiments. R.M. wrote the first manuscript draft with contributions from
R.G.-B. and M.B. All authors participated in the discussion and interpretation
of the results, and contributed to the manuscript preparation and revision.

##  6 Competing interests

The authors declare no competing interests.

[HTML]F4F4F4 Formulas |  Si/Al |  Al |  Si |  [Cu(NH3)2]+ |  NH4+ |  NH3 |  H+
---|---|---|---|---|---|---|---
H2Al2O192Si94 |  47 |  2 |  94 |  0 |  0 |  0 |  2
H8Al2N2O192Si94 |  47 |  2 |  94 |  0 |  2 |  0 |  0
H12Al3N3O192Si93 |  31 |  3 |  93 |  0 |  3 |  0 |  0
H15Al3N4O192Si93 |  31 |  3 |  93 |  0 |  3 |  1 |  0
H18Al3N5O192Si93 |  31 |  3 |  93 |  0 |  3 |  2 |  0
H28Al7N7O192Si89 |  12.7 |  7 |  89 |  0 |  7 |  0 |  0
H12Al2Cu2N4O192Si94 |  47 |  2 |  94 |  2 |  0 |  0 |  0
H10Al2Cu1N3O192Si94 |  47 |  2 |  94 |  1 |  1 |  0 |  0
H18Al2Cu2N6O192Si94 |  47 |  2 |  94 |  2 |  0 |  2 |  0
H18Al3Cu2N6O192Si93 |  31 |  3 |  93 |  2 |  1 |  0 |  0
Table S1: Chemical composition, cationic species and molecules included in the triclinic T96O192 supercell models used for active learning and adversarial attack. [HTML]F4F4F4 Al distribution |  $\Delta G_{\mathrm{act}}$ (kcal/mol) |  $\Delta G$ (kcal/mol)
---|---|---
SR1 |  3.9 |  0.4
SR2 |  4.6 |  1.7
SR3 |  5.4 |  2.3
SR4 |  4.1 |  0.5
DR1 |  6.1 |  3.2
DR2 |  6.7 |  2.9
DR3 |  6.4 |  3.2
DR4 |  5.4 |  -0.7
S4R |  7.3 |  4.5
S6R |  6.4 |  5.4
Table S2: Free energies of activation ($\Delta G_{\mathrm{act}}$) and reaction
($\Delta G$) for the diffusion of a [Cu(NH3)2]+ complex between two
neighboring cages through 8R windows with different Al distribution obtained
from NNP-based biased simulations at 423 K.

[HTML]F4F4F4 Name | Formulas | Si/Al | Al | Si |  [Cu(NH3)2]+ |  NH4+ |  NH3 |  H+
---|---|---|---|---|---|---|---|---
L20 |  H144Al26Cu20N46O1536Si742 | 28.5 | 26 | 742 | 20 | 6 | 0 | 0
L4 |  H112Al26Cu4N30O1536Si742 | 28.5 | 26 | 742 | 4 | 22 | 0 | 0
M20 |  H240Al50Cu20N70O1536Si718 | 14.3 | 50 | 718 | 20 | 30 | 0 | 0
M4 |  H208Al50Cu4N54O1536Si718 | 14.3 | 50 | 718 | 4 | 46 | 0 | 0
M20-H+ |  H150Al50Cu20N40O1536Si718 | 14.3 | 50 | 718 | 20 | 0 | 0 | 30
M20-NH3 |  H423Al50Cu20N131O1536Si718 | 14.3 | 50 | 718 | 20 | 30 | 60 | 0
H6R |  H312Al68Cu20N88O1536Si700 | 10.3 | 68 | 700 | 20 | 48 | 0 | 0
H8R |  H312Al68Cu20N88O1536Si700 | 10.3 | 68 | 700 | 20 | 48 | 0 | 0
HR |  H312Al68Cu20N88O1536Si700 | 10.3 | 68 | 700 | 20 | 48 | 0 | 0
HB |  H312Al68Cu20N88O1536Si700 | 10.3 | 68 | 700 | 20 | 48 | 0 | 0

Table S3: Chemical composition, cationic species and molecules included in the triclinic T768O1536 supercell models used for NNP-based regular molecular dynamics simulations.  [HTML]F4F4F4 Sample | Si/Al |  Al/cha cage |  Crystal size (nm) |  BET surface area (m2/g) |  Micropore area (m2/g) |  Micropore volume (cm3/g)
---|---|---|---|---|---|---
CHA07 | 7.3 |  1.4 |  $\sim$100 |  562 |  532 |  0.26
CHA13 | 12.7 |  0.9 |  $\sim$400-800 |  556 |  551 |  0.27
CHA23 | 23.3 |  0.5 |  $\sim$100 |  550 |  490 |  0.25
Table S4: Chemical composition and physico-chemical properties of the synthesized CHA zeolites. [HTML]F4F4F4 Sample |  Si/Al |  %wt Cu |  Cu/Al |  Al/cha cage
---|---|---|---|---
CHA07_1.5Cu |  7.3 |  1.58 |  0.11 |  0.16
CHA07_3.0Cu |  7.3 |  2.92 |  0.20 |  0.30
CHA13_1.5Cu |  12.7 |  1.55 |  0.18 |  0.17
CHA23_1.5Cu |  23.3 |  1.54 |  0.35 |  0.17
Table S5: Chemical composition of the Cu-exchanged CHA zeolites.

(a)

(b)

(c)

Figure S1: Representation of the three types of base lattices used for
umbrella sampling simulations (a and b) and for unbiased MD simulations (c).
The representations correspond the all-silica CHA, however, the models used in
this work include different combinations of Al number and locations starting
from this base lattices, as explained in the main text. Color code: Si and O
are depicted as orange and red respectively. Figure S2: Diffusion of
[Cu(NH3)2]+ complex from cage A ($\xi<0$) to cage B ($\xi>0$). (b)
Representation of the collective variable $\xi$ describing the [Cu(NH3)2]+
diffusion through the 8R window (transparent yellow plane) Figure S3: UMAP
plot for the feature vectors of Al atoms from $\sim$1% of the dataset. Both
axes are on the same scale. Figure S4: Scatter plots showing the regions
visited by Cu+ in the simulations run at 500 K for 5 ns. Color code: Si, O,
and Al atoms are depicted as orange, red and light brown. Each Cu+ cation is
represented with a random color. Figure S5: Scatter plots showing the regions
visited by NH4+ in the simulations run at 500 K for 5 ns. Color code: Si, O,
and Al atoms are depicted as orange, red and light brown. Each NH4+ is
represented with a random color. Figure S6: Time evolution of the average
number of [Cu(NH3)2]+ pairs in the same cage for the models M20 (red), M20-H+
(black), M20-NH3 (brown), in unbiased MD simulations at 500 K. Figure S7: PXRD
patters of the synthesized CHA zeolites Figure S8: FESEM images of the
synthesized CHA zeolites

[H]  Representation of long-range [Cu(NH3)2]+ migration assisted by NH4+. Al,
Cu+ and NH4+ are represented with yellow, green and blue balls, respectively.