Migration of zeolite-encapsulated subnanometre platinum clusters via reactive neural network potentials

Abstract

The migration of atoms and small clusters is an important process in sub-nanometre scale heterogeneous catalysis, affecting activity, accessibility and deactivation through sintering. Control of migration can be partially achieved via encapsulation of sub-nanometre metal particles into porous media such as zeolites. However, a general understanding of the migration mechanisms and their sensitivity to particle size and framework environment is lacking. Here, we extend the time-scale and sampling of atomistic simulations of platinum cluster diffusion in siliceous zeolite frameworks, by introducing a reactive neural network potential of density functional quality. We observe that Pt atoms migrate in a qualitatively different manner from clusters, occupying the dense region of the framework and avoiding the free pore space. We also find that for cage-like zeolite CHA there exists a maximum in self diffusivity for the Pt dimer beyond which, confinement effects hinder intercage migration. By extending the quality of sampling, NNP-based methods allow for the discovery of novel dynamical processes at the atomistic scale, bringing modelling closer to operando experimental characterization of catalytic materials.

---

1. Introduction

Incorporation of noble and coinage metal atoms and clusters, such as Ir, Pd, Pt and Cu into zeolites has become a common strategy for generation of stable, monodisperse active sites for heterogeneous catalysis.^1–5 Electrostatic stabilisation may be achieved via pinning of partially oxidised atoms and clusters to the framework via the charge imbalance generated by substitution of silicon atoms with heteroatom sites (e.g. Al, B, Ga).^6–8 Alternatively, physically trapping metallic clusters within the cages of zeolites has been proposed to stabilize particles against sintering.⁹ Nevertheless, for reduced atoms and clusters, sintering is energetically preferred, and the stability of encapsulated particles against sintering is often low, with rapid loss of internal metal loading over repeated catalytic cycles or at elevated temperatures.¹⁰

For catalytic metal species: atoms, clusters and nanoparticles, the loss of active site density through thermodynamically favourable sintering of particles is the primary deactivation route. Diffusion also plays a role in the kinetics of reactions – as access to active sites is affected by the migration of particles. This can even manifest during the course of the catalysed reaction itself, as has been proposed for the zeolite-catalysed selective catalytic reduction of NO_x in automotive emissions abatement.¹¹ Hence, the maximisation of sintering resistance, through the suppression of particle migration is an important goal for the optimization of nanocatalyst activity and longevity.¹²

Characterization of zeolite-encapsulated sub-nanometre clusters is experimentally challenging, owing to the size resolution limitations of most imaging methods, the sensitivity of zeolites to beam damage in electron microscopes and the lack of sensitivity of spectroscopic methods towards the low loadings of metal (often <0.1 wt%).^13–15 High resolution TEM imaging remains a preferred method for characterising metal atoms and clusters^2,7 and has been successfully employed for detecting and determining cluster size, coordination environment and oxidation state, in cooperation with EXAFS and XPS.⁸ Nevertheless, ex situ characterization of as-prepared or spent catalyst is restricted to a less than fully realistic depiction of the system under operando conditions and is unable to directly interrogate dynamic processes such as diffusion.

The diffusive behaviour of noble and late transition metal clusters supported on metal oxide surfaces has been explored extensively via computational methods in the last two decades. Several reports have noted that the self-diffusivity of clusters is non-monotonic in cluster size.^16,17 These findings establish that the most diffusive particles may be intermediate cluster sizes, and thus the mechanism and kinetics of Ostwald ripening and particle migration processes may be controlled by small clusters, rather than single atoms.

A thorough, atomistic understanding of the migration processes that control sintering is necessary. However, in theory, diffusion modelling has been limited by both the precision of the method and the short duration of the simulations. For electronic structure calculations, such as ab initio molecular dynamics, the simulations are generally limited to the tens of picoseconds timescale, while classical dynamics is usually performed with simplified, non-reactive models. The understanding of metal–zeolite interactions is still incomplete, but it has been shown by several studies that methods of quantum chemical accuracy are required to capture the physical nature of these interactions.¹⁸ For example, reversible bond breaking events have been shown to occur spontaneously in several zeolites when interrogated with small metal clusters.^8,14

The expense of quantum chemical modelling has so far limited the direct study of diffusive processes to short timescales, with the exception of recent studies on molecular diffusion in zeolites.¹⁸ For metal particles, which bind strongly to the zeolite framework, direct ab initio simulations of sufficiently long duration to observe rare events such as intercage migration are not feasible. However, a promising recent approach to bypass this problem is the development of artificial neural networks, which generate forcefields that are trained upon quantum chemical data, thus maintaining a high-quality description of the underlying potential energy surface, while reducing the cost of simulations by several orders of magnitude.^19,20 This approach has been applied, for example, by some of the current authors toward reactive transformations in silicate and aluminosilicate zeolites,²¹ by Yu et al. towards small metal clusters in MOFs²² and very recently by both Ma et al. for binary metal oxide clusters in zeolites for catalysis of propane dehydrogenation²³ and Millan et al., for the diffusion of Cu⁺ ions in the presence of ammonia during the selective catalytic reduction process in zeolite CHA.²⁴

In this work, we present an investigative study into the diffusive behaviours of ultrasmall metal clusters inside siliceous zeolite frameworks. We develop a reactive neural network forcefield which models metal–zeolite interactions at the (generalized gradient + dispersion) DFT level and we apply it to the study of sub-nanometre Pt in zeolites CHA, TON, MWW and MFI. By extending the simulation timescale by several orders of magnitude, we are able to observe hitherto unreported migration mechanisms of the Pt atom, and clear differences between atoms and clusters that will have consequences on their accessibility and activity as confined nanocatalysts.

2. Results and discussion

2.1 Development of NNP

Training of accurate NNPs first requires the generation of a structurally diverse DFT database for robust interpolation of the Pt-silica potential energy surface (PES). The starting point of the NNP development was the DFT database (PBE+D3 level) from a previous work²¹ containing about 33k structures covering the silica PES from low-density, siliceous zeolites to high-density silica polymorphs. We extended the DFT database using structures from AIMD and global optimization (GO) runs of Pt_n (n = 1–5) in siliceous LTA and the gas phase (taken from ref. 8). Additional AIMD simulations (10 ps at 1500 K) were performed for gas phase models of Pt₁₉ and Pt₃₈ as well as for bulk Pt to further diversify the DFT database. All AIMD and GO trajectories were subsampled by farthest point sampling (FPS) using the smooth overlap of atomic positions (SOAP) descriptor as similarity metric of the Pt atomic environments. This selection procedure (see ref. 21 for details) allowed us to collect a subset of 23k structures with highly diverse Pt environments for spin-polarized DFT single-point calculations.

The created database was used to train an initial ensemble of six NNPs, enabling active learning to further extend the database with configurations of sparsely interpolated or yet unseen parts of the PES.^20,21,25 This is achieved by monitoring the deviation of the six energy predictions from the average ensemble prediction which acts as estimator of the prediction accuracy during NNP simulations. We performed 70 NNP level MD runs (1 ns at 1000 K and 2000 K) for Pt_n (n = 1–5) in TON, MOR, MWW, LTA, CHA, MFI, and FAU (primitive unit cell). These MD runs sampled large parts of the PES of Pt in industrially relevant zeolites, including highly energetic reactive events. To also cover lower energy parts of the PES, local optimizations were applied to 1000 structures extracted from each MD trajectory. All configurations were collected if the energy predictions deviated by more than 10 meV per atom from the ensemble average. The resulting structure set was then subsampled by FPS followed by DFT single-point calculations to extend the initial database.

The final DFT database contains about 73k structures with 33k silica configurations, 28k Pt containing zeolites, 11k gas phase Pt clusters and 500 bulk Pt structures. The SchNet NNPs trained on this database learned a latent representation of atomic environments encoding structure, chemical composition, and energetics of the training set.^19,26 We averaged the atomic representation vectors of each structure in the database and applied T-distributed stochastic neighbor embedding (t-SNE) to visualize the NNPs’ capability to learn chemical and structural features from DFT data. The 2D t-SNE components show not only a clear correlation with DFT calculated energies (Fig. 1a) but also reasonably group systems with similar chemical composition and distinguish between bulk and gas phase Pt models (Fig. 1b). Both t-SNE maps demonstrate that the trained NNPs successfully learned a latent representation describing large parts of the Pt PES in the gas phase and in siliceous zeolites.

Fig. 1 t-SNE maps of the DFT database using the SchNet representation vectors averaged over all atoms of a structure (each point corresponds to one configuration). The colour codes represent (a) energies and (b) purely siliceous structures (black), Pt containing structures (Pt_n@zeo, red shades), Pt_n gas phase models (blue shades) and bulk Pt (purple).

2.2 Performance of NNP

To quantify the NNP accuracy, Fig. 2 shows the NNP energy error distributions with respect to their DFT reference for different test cases. Fig. S1† shows the error analysis in diagonal form for energies and forces. Table 1 summarizes the corresponding root mean square errors (RMSE) of energies and forces (see Table S1† for mean absolute errors). The database test set consists of about 7300 configurations randomly chosen from the generated DFT database. These structures were excluded from the NNP training procedure and show an overall RMSE of 7.5 meV per atom and 141 meV Å⁻¹ for energy and forces, respectively. We split the test set into three categories: purely siliceous structures, Pt containing zeolites (Pt@zeolite) and pure Pt configurations (bulk and gas phase). The NNPs have approximately the same quality for siliceous and Pt containing zeolites and the RMSEs are about the same as for the previously published pure silica NNPs,²¹i.e., the NNP extension did not result in any accuracy loss.

Fig. 2 Energy error ΔE distributions for (a) the test set taken from the DFT database and (b) structures taken from MD trajectories of small Pt cluster in CHA and on silicatene, e.g., Pt₆@silicatene shown in (c).

Table 1 Root mean square errors (RMSE) of energies and forces for different systems taken from the DFT database and MD production runs

RMSE		Energies [meV per atom]	Forces [meV Å^−1]
Database test set	Silica	4.6	135
	Pt@zeolite	3.3	156
	Pt (gas + bulk)	16.9	114
MD test sets	Pt1,3,5@CHA	0.9	98
	Pt1@silicatene	7.4	100
	Pt3@silicatene	7.5	99
	Pt6@silicatene	8.7	115
	Pt9@silicatene	11.2	149

---

In particular, the NNPs accurately approximate both energy (3.3 meV per atom RMSE) and forces (156 meV Å⁻¹ RMSE) in the case of Pt@zeolite which is the main focus of this work. Only the pure Pt structures (bulk and gas phase clusters) show a higher energy RMSE of about 17 meV per atom, but similarly low force errors. One factor leading to the larger energy errors is the definition of the loss function which puts high weight on forces and energies of structures with a higher number of atoms (see Computational details). Another reason is that the NNP models do not consider different spin states of small Pt gas phase clusters but were trained on spin polarized DFT data. It was shown that spin multiplicities have a pronounced effect on the relative energies of small Pt₃–Pt₆ cluster isomers in the order of 30–100 meV per atom.^8,27 Therefore, the trained NNPs are an approximate closed-shell interpolator of the Pt PES with limited accuracy for energy predictions of gas phase clusters but with low errors for Pt@zeolite interactions. In addition, the NNPs accurately model forces for all test cases which is most relevant for dynamical simulations of, e.g., particle migration. Therefore, the pure Pt training data rather acts as limiting case to improve the NNP robustness for generalization to systems unseen during NNP training.

It has been found for previously generated potentials with this architecture that the NNPs are more accurate for the reproduction of dynamic energetics than statics, owing to error cancellation. Nevertheless, we tested the ability of the NNPs to reproduce geometric configurations for relevant Pt₁@CHA local minima (Fig. S2†) showing good agreement between the NNP and DFT for bond lengths and angles. Nudged elastic band calculations were also performed to determine the ability of the NNPs to reproduce reactive processes for migration paths of Pt₁: (i) through the centre of the double six-ring d6r, with a barrier of 70 kJ mol⁻¹, which corresponds to a 5.4% reduction when compared with the DFT barrier (74 kJ mol⁻¹), (ii) migration of Pt₁ through the 8-ring of the LTA zeolite, with a barrier of 142 kJ mol⁻¹, which is a 12.9% reduction from the DFT barrier (163 kJ mol⁻¹ (ref. 8) and iii) the NEB pathway of the initial process of Pt₁ migration in CHA, which contains two elementary steps. It has an effective barrier of 111 kJ mol⁻¹ at the NNP level, and 138 kJ mol⁻¹ at the DFT level. This underestimation of 19.7% is almost entirely due to a single local minimum along the pathway, with most stationary points in close agreement between NNP and DFT (Fig. S3†) and the same rate determining step – the initial transfer of Pt out of the plane of the 6-ring. Thus, we conclude that the NNPs are able to reproduce static configurations, energetics and barriers to a good precision in relevant cases, and they perform well in dynamic simulations.

To test the generalization capability of the NNPs, we performed additional MD simulations (500 ps, T = 1000 K) for Pt_n (n = 1, 3, 6, 9) on a silica bilayer (silicatene), which we note is structurally highly distinct from a three-dimensional zeolite model, and was not included in the training dataset, thereby representing a stern test. Fig. 2b compares the energy error distributions of Pt_n@silicatene and Pt_n@CHA (n = 1, 3, 5, T = 1000 K) calculated for a subset of the MD trajectories. Table 1 also summarizes the energy and force RMSEs for this MD test set. In the case of Pt_n@CHA we found low energy and force errors since such systems were part of the training dataset. In contrast, Pt_n@silicatene show slightly higher errors which increase with cluster size. The shift of the energy error distributions shows a systematic NNP error of total energies but the narrow distributions indicate a good NNP performance for relative energies (e.g., with respect to an optimized reference configuration). In addition, the energy errors in the order of 10 meV per atom are still much smaller compared to, e.g., ReaxFF for siliceous systems (about 100 meV per atom).²¹ Finally, the NNPs accurately model forces (RMSE of less than 150 meV Å⁻¹) even for Pt₉@silicatene, which is of particular importance for modelling dynamic averages of silica supported Pt from MD simulations.

The NNPs are thus shown to be robust, transferable among zeolite topologies and able to reproduce the DFT reference with high fidelity along with reasonable performance in out-of-domain tests. They appear to be a useful tool for extending the timescale and scope of computational investigations of cluster diffusion in zeolites.

2.3 Migration of Pt₁, Pt₃ and Pt₅ in CHA

As an archetypal small pore zeolite, of great importance in the catalytic industry,^28,29 we chose CHA as the first zeolite topology to consider. We investigated the nature of Pt migration dynamics inside the pore of siliceous CHA, which is comprised of four, six and eight rings, via 25 ns canonical molecular dynamics simulations with our NNP. Preliminary simulations showed that 25 ns simulations at room temperature are insufficient to induce migration events. Therefore, we consider elevated temperatures of 750 K, 1000 K and 1250 K. These temperatures are within the range in which spontaneous framework destruction and volatilization of the Pt cluster is avoided, and only extend mildly beyond common catalytic applications or calcination temperature ranges for zeolites.^30,31 We apply these temperatures to enhance sampling of Pt migration, while retaining the underlying physical nature of the system relevant at lower temperatures. Results from 750 K and 1000 K are shown in the ESI,† and reflect the observation that at 1250 K, the same qualitative behaviours are observed, with higher frequency.

Pt₁ is found to be strongly attracted to the 6-ring plane and occupies this site almost exclusively throughout the simulation (Fig. 3). The symmetry-equivalent O–Pt–O sites in each 6-ring plane of the double 6-ring (d6r) are occupied equally. There is a frequent local transition between adjacent 6-ring faces of the d6r. This transition is achieved via a reaction pathway that involves Pt insertion into a framework Si–O bond, with concomitant defect formation in the 6-ring, followed by migration of Pt into the internal space of the d6r, with healing of the framework Si–O bond. Pt then migrates to the face of the second 6-ring, inserting into another Si–O bond (see ESI,† animation of Pt₁ in d6r of CHA). Such breakage of the Si–O bond by Pt insertion has been observed to be spontaneous for Pt atoms bound to the single six rings of the siliceous LTA framework,⁸ and also in the case of the single Ir atom binding in zeolite MWW.¹⁴ This pathway was found to have a moderate barrier of 70 kJ mol⁻¹via nudged elastic band (NEB) calculations. Further analysis of this transition pathway via umbrella sampling of the free energy path at 300 K provided a barrier of 62 kJ mol⁻¹ (Fig. 4), suggesting that this is a feasible transition even at room temperature. It is worth noting, however, that this is not a productive pathway for diffusion, as Pt₁ must leave the d6r to move through the pore system. The feasible d6r motion of Pt was additionally investigated via alternative methods and is described in the ESI (Fig. S4†). A negligible occupation probability density for Pt₁ is observed outside of the d6r. This implies that in siliceous CHA, the available 8-ring sites which interconvert adjacent cages are to a large extent unoccupied. Hence, Pt₁ does not spend time in the free pore volume of the CHA cage, but rather is trapped in the denser d6r region of the framework.

Fig. 3 Isosurfaces of the average Pt occupation density in CHA at 1250 K from 25 ns MD trajectories for (a) Pt₁ and (b) Pt₃. Pt₁ is observed to occupy the broken 6-ring face of the d6r (right). Pt is in grey, silicon in yellow and oxygen in red.

Fig. 4 Helmholtz free energy profiles for inter-cage cluster migration of, Pt_x in CHA from umbrella sampling at 300 K. For x = 2–5, this migration is through the 8-ring, while for x = 1, it is through the d6r. The x-axis is given in scaled units spanning the range of the chosen CV (Fig. S11†). A reasonable free energy pathway for Pt₁ was not found through the 8-ring. The (athermal) NEB barrier for inter-cage migration is included (at NNP level) via the pathway shown in Fig. S3.† The effective barrier is given by the black dotted line.

Previous calculations for Pt clusters encapsulated into the topologically similar cage-like zeolite LTA, showed that while Pt₁ lies in the plane of the 6-ring, forming strong Pt–O bonds to the framework, platinum clusters (n > 1) interact more weakly and occupy the free volume of the pore.⁸ This difference in binding mode suggests a difference in the migration behaviour. We investigated this in CHA via molecular dynamics simulations. Mean squared displacements of the clusters are shown in Fig. S5.† Platinum clusters Pt₃ and Pt₅ are indeed found to differ in their migration behaviour from the Pt atom in CHA. At both 750 K and 1250 K, the clusters move freely within the CHA cage, adopting liquid-like configurations. No evidence of Pt–Pt dissociation was observed over the duration of the simulations. Migration of the clusters through the pore system occurs via transitions between adjacent cages through the 8-ring, which are rare events, even at elevated temperatures (see ESI,† animation of Pt₃ jump through 8-ring). Overall, the clusters are found to diffuse in a qualitatively opposite manner to the single atom. While the clusters occupy the free space in the zeolite and move through the pore via the largest available ring, which in CHA is the 8-ring, the Pt atom occupies the dense phase of the zeolite, hopping between faces of the d6r. This difference in location and migration mechanism is likely to have pronounced effects on the accessibility of the cluster towards other encapsulated species, such as gas phase adsorbates in catalytic processes.

In order to obtain the free energy barrier for the inter-pore migration processes, umbrella sampling calculations were performed for the migration of Pt₂, Pt₃, Pt₄ and Pt₅ through the 8-ring at 300 K according to a one-dimensional collective variable which varies the cluster centre of mass along a coordinate which passes through the plane of the 8-ring (Fig. 3). The larger clusters, Pt₄ and Pt₅ exhibit high barriers of 140 and 149 kJ mol⁻¹, respectively, in accordance with the limited mobility observed in the MD simulations for Pt₅.

These clusters are likely to be trapped within the cage under mild conditions and contribute little to the overall cluster diffusivity. Pt₃, which is sterically less hindered than the tetramer or pentamer when passing through the 8-ring, has a lower migration barrier of 125 kJ mol⁻¹, but this is still a highly activated rare event. Interestingly, the dimer has a substantially lower barrier of 54 kJ mol⁻¹ (Table 2 and Fig. S6, S7†). This barrier is lower than that of the short-range Pt atom migration between faces of the d6r (62 kJ mol⁻¹). This is due to a combination of two factors: a weaker interaction with the framework than the single atom, in which Si–O bonds are not disrupted by platinum, and a sufficiently small size to move through the 8-ring without significant steric clashes. Finally, we did not observe a free energy pathway for Pt₁ to migrate directly through the 8-ring, as Pt₁ is found to prefer to hop between 6-ring faces. An athermal NEB calculation of a pathway through the 8-ring was calculated instead, as an approximation of the barrier for the direct 8-ring migration process, and found to have a high barrier of 111 kJ mol⁻¹. This pathway is detailed in Fig. S3.† Diffusivity is therefore non-monotonic in cluster size inside CHA. The same trend is observed for Pt_x (x = 1–4) inside the similar zeolite, LTA (Fig. S6†).

Table 2 Helmoltz free energies for migration of Pt clusters in siliceous CHA at 300 K. For Pt₁ the barrier is for a local transition between 6-rings. The asterisk denotes a NEB barrier (athermal potential energy) for the intercage transition for Pt₁

Zeolite	Cluster	Barrier/kJ mol^−1
CHA	Pt1	62
	Pt1*	111*
	Pt2	54
	Pt3	125
	Pt4	140
	Pt5	149

---

2.4 Effect of zeolite topology on Pt migration

CHA is a cage-like zeolite, which presents a clear hierarchy of diffusion limitations. Inter-cage migration is controlled (for the clusters) by migration through the 8-ring which connects adjacent pores. A similar topological arrangement is in place in zeolite LTA (biased dynamical simulations are presented in the ESI†). However, industrially relevant zeolites adopt a range of topologies, including 1D, tubular architectures (TON, ABW), 2D layered structure (UTL, MWW) or intersecting channels (MFI). Our neural network potentials were trained to retain flexibility across silicate topologies. Hence, we consider the role of framework topology on the diffusive properties of Pt₁, Pt₃ and Pt₅, by selecting TON, MWW and MFI as representative of a broad range of zeolites. MD trajectories are available under: zenodo.org/doi/10.5281/zenodo.10371972.

In all three tested frameworks, irrespective of framework density or channel connectivity, Pt₁ is observed to preferentially occupy the dense regions of the framework, and to migrate through the denser regions of the unit cell, avoiding the open pore volume, in agreement with the findings for Pt₁@CHA.

In TON, which is comprised of a zigzag one-dimensional channel made up of distorted 10-rings, surrounded by a dense phase containing bent 5-rings and bent 6-rings, the Pt atom accesses the dense region via the bent small rings. Sequential hops between these rings allows the Pt atom to move through the framework. At 1250 K, the diffusivity is notably higher than observed in CHA, which can be observed via the more diffuse occupation densities in Fig. 5 and the mean squared displacement variations over time in Fig. S5.† This increased diffusivity is facilitated by the easier migration of Pt₁ through the bent rings of TON and the fact that 5 and 6 rings in TON are directly connected to eachother, in contrast to the perfectly planar CHA d6r, which is isolated from other d6rs.

Fig. 5 Occupation isosurfaces from MD for Pt₁ (top) and Pt₃ (bottom) in zeolite TON. Pt₁ occupies the extensive dense region of the framework, while Pt₃ diffuses along the 1D channel.

MFI consists of four, five, six and ten rings, and has a more complex pore system comprised of intersecting, perpendicular zigzag and straight channels. The dense phase consists of small cages, which may be accessed via bent 6-rings (cas, mel) and 5-rings (mfi/t-pen). The MD simulations at lower temperatures (750 K and 1000 K), showed Pt₁ to primarily occupy 6-rings which are accessible to the pores (Fig. 6 and S7†). An equilibrium exists between fully connected and partially broken 6-rings, as is the case for Pt₁ in the 6-rings of CHA and LTA. Hops between 6-rings are rare events. An additional migration mechanism is observed, which is a rare event on the timescale of the simulations, in which Pt₁ occupies the interior of the t-pen cage, which is accessed via a ring-opening mechanism that is similar to that observed in the 6-rings of CHA (Fig. 3). Pt binds to a five ring and breaks an Si–O bond, forming an O–Pt–O moiety and expanding the ring. It subsequently moves inside the cage, and the framework heals behind, trapping the Pt atom, which can form multiple bonds to nearby framework oxygen atoms in the t-pen unit and achieve coordinative saturation (see ESI,† animation of Pt₁ insertion into t-pen unit). Migration of Pt between unit cells is thereby primarily achieved by migration between t-pen cages, mediated by adjacent 6-rings. Thus, at higher temperatures the accessibility of platinum may decrease, as Pt₁ occupies t-pen cages more frequently.

Fig. 6 Occupation isosurfaces from MD trajectories of Pt₁ (top) and Pt₃ (bottom) in zeolite MFI at 1250 K. Pt₁ occupies t-pen cages, while Pt₃ diffuses through both straight and sinusoidal channels.

MWW consists of a stack of dense 2D siliceous layers comprised of mel composite building units, connected together via a d6r. Two disconnected, inequivalent 2D gallery spaces are available: one containing the d6r, in which the accessible space is bounded by 5-rings that make up the mel units, known as the in-layer gallery, and one containing a pocket, which is terminated by the face of the d6r, known as the between-layer gallery. Our simulations indicate that Pt₁ has a strong preference for the 5-ring of the mel groups (Fig. 7), to which the Pt atom binds, forming an O–Pt–O moiety, which breaks an Si–O bond of the ring (Fig. S8 and S9†). This has previously been predicted to be the stable binding mode of the Ir atom in MWW.¹⁴ Migration through the unit cell in the lateral directions takes place via hops between mel 5-rings, without significant interaction with the d6r unit. It is notable that the disconnected 2D galleries do not restrict the Pt atom from reaching the preferred mel 5Rs, as migration along the c direction, through the dense phase between layers is feasible, both at 750 K and 1250 K. Hence, Pt₁ can effectively tunnel between disconnected regions of the zeolite. This can be seen clearly by comparison of two simulations, denoted MWW and MWW* in Fig. 7, in which Pt₁ isosurfaces are similar, despite Pt originating in different layers.

Fig. 7 Occupation isosurfaces from MD for Pt₁ (top) and Pt₃ (bottom) in zeolite MWW. Initial cluster locations are in the “in-layer” gallery (left) and the “between-layer” gallery (right).

By comparison between diverse zeolite frameworks CHA, TON, MWW and MFI, we conclude that for Pt atoms in siliceous frameworks, neither the void space, nor the connectivity of the pore system appear to be the crucial factors in controlling diffusivity or migration mechanism. It is the connectivity, shape and size of rings which dictate the Pt atom location and migration pathway, in particular, the presence or absence of 5- and 6-rings. The sequential formation and breakage of Pt–O and Si–O bonds are used to migrate through the zeolite via the dense phase of the siliceous zeolite.

As in the case of Pt clusters in CHA, Pt₃ and Pt₅ inside TON, MFI and MWW are observed to behave in a qualitatively similar way to each other, but distinctly from the Pt atoms. Migration occurs entirely within the free pore volume, with no dissociation of the cluster, nor disruption of the framework Si–O network. In TON, the Pt₃ cluster moves along the one-dimensional channel, with the occupation density map describing the void space (Fig. 5 and S10†). In MFI, diffusion occurs along both the straight and sinusoidal channels, and the t-pen cages are unperturbed by the cluster (Fig. 6), which is too large to open the ring and move inside. Interestingly, no preference is found for either channel type, but there is a reduced occupation of the intersection for Pt₃ and an increased occupation of the intersection for Pt₅. This difference may be due to the fact that Pt₅ is more strongly confined in the channels than Pt₃. In MWW, the inability for Pt₃ to tunnel through the dense phase means clusters remain trapped within the 2D interlayer space into which they are placed. However, the cluster ergodically samples the available free space of a unit cell within 25 ns at 1250 K. This weak cluster-framework interaction is in agreement with the recent predictions for PtSn clusters in MFI by Ma et al.²³

Comparing the diffusivity of the clusters across the topologies (Fig. S5†) we observe that migration is fastest in MFI and slowest in CHA, with TON and MWW intermediate. For clusters, the high diffusivity in MFI is explained by the highly interconnected channel network in three dimensions, while CHA traps clusters within the cage, making long-distance migration a highly activated rare event. TON is intermediate, as it has only one-dimensional channels, but little hindrance to cluster motion along the channel. MWW has two dimensions for cluster migration, but the landscape is hindered by the presence of pillars between layers made up of d6rs. From Fig. 3–7, it is clear that Pt₁ and Pt_x (x > 1) form two distinct regimes in terms of their diffusive behaviour, with occupation densities that are near mirror images of eachother, and that this trend is independent of the zeolite topology.

3. Conclusions

We have successfully developed a reactive neural network potential for zeolite-encapsulated sub-nanometre sized Pt clusters, that extends the accessible range of dynamical simulations by several orders of magnitude, while retaining the quality of dispersion-corrected GGA DFT. This potential is transferrable to a wide range of topologically distinct zeolites, even extended two-dimensional surfaces with reasonable accuracy. By extending the simulation timescale, and further accelerating via biased dynamics, we were able to determine the preferential locations and novel migration mechanisms for several small Pt clusters. We observed that Pt atoms and clusters form two distinct classes, moving according to different pathways. For the atom, the defining feature is framework pinning, while for larger clusters, physical confinement, such as within cages controls diffusion. The structure of the dense phase controls the location, accessibility and diffusivity of the single Pt atom, implying that the paradigm of solid, impassable silicate walls that control the movement of intercalated atoms must be reconsidered, in favour of a more dynamic view, in which metal atoms interact chemically with the framework and may tunnel through dense regions. The nature of the dense region, e.g. presence of 5 and 6-rings, and their connectivity defines the diffusivity of Pt₁. Reactive, bond-breaking events are common between Pt atom and siliceous frameworks and occur commonly as part of migration routes, while these reactive events were not observed for the Pt clusters. Five and six rings can be used as access points into the dense region of zeolites, which are favoured locations for single Pt atoms. The different preferred locations and migration routes for atoms and clusters leads to non-monotonic trends in migration barriers, suggesting that the metal dimer may be more diffusive than larger clusters and single atoms in siliceous zeolite frameworks.

4. Experimental/computational details

NNPs were generated via the SchNetPack software within the atomistic simulation environment (ASE). Dynamical simulations were performed using ASE. The SchNet NNP setup¹⁹ used six interaction blocks with a feature vector size of 128, a cutoff radius of 6 Å and 60 Gaussians for expansion of pairwise distances. Training of an ensemble of six NNPs employed the Python package SchNetPack³² with random splits of the reference DFT database into training (80% of the database), validation (10%), and test sets (10%). Mini-batch gradient descent optimization with the ADAM optimizer was applied (batch size of eight structures) for minimization of same squared loss function for energy and forces as in ref. 19 along with a trade-off factor of 0.01 (high weight on forces). The learning rate (between 10⁻⁴ and 10⁻⁶) was halved if the validation loss did not decrease after 10 epochs. MD simulations at the NNP level were performed within the atomic simulation environment (ASE)³³ with 1 fs time step and Nosé–Hoover thermostat.^34,35 The NNPs and all MD trajectories are available under: zenodo.org/doi/10.5281/zenodo.10371972.

Biased dynamical simulations were performed within ASE, using the PLUMED plugin.^36–38 Umbrella sampling was performed using a minimum of 50 ps simulations across a minimum of 25 equally spaced windows along a collective variable based on the method of Cnudde et al.,¹⁸ in which the cluster moves along a vector normal to the plane of a ring (Fig. S11†). All biased dynamics simulations were performed at 300 K. The weighted histogram method was used to generate free energy pathways from umbrella sampling simulations.³⁹ Test calculations to determine the role of parameter choice and the extent of numerical variation gave an estimated maximum error of ±5 kJ mol⁻¹.

All density functional theory (DFT) calculations were performed using the VASP 5.4 code,^40,41 with the exchange–correlation functional of Perdew, Burke and Ernzerhof⁴² and an additional D3(BJ)^43,44 dispersion correction. For convergence of electronic states, Gaussian smearing was applied with a width of 0.1 eV to smooth the Fermi–Dirac distribution. Wavefunctions were described by a plane-wave basis with a kinetic energy cutoff of 450 eV. The electronic structure was calculated at the gamma point only. Local minimization of geometry was performed with the conjugate gradient method, and convergence was defined by an energy and force tolerance of 10⁻⁶ eV and 0.01 eV Å⁻¹, respectively. For ab initio molecular dynamics (AIMD) simulations, the canonical ensemble was chosen, with target temperatures set to 750, 1000 and 1250 K, which was controlled by the Nosé–Hoover thermostat.

Migration barriers were calculated by locating transition states between pairs of adjacent local minima on the potential energy surface, via the climbing image nudged elastic band method provided within the transition state tools package of the VASP software package for DFT calculations.⁴⁵ For barriers calculated at the NNP level, the NEB module of the atomistic simulation environment was used to connect local minima optimized at the NNP level.

Author contributions

CJH: conceptualization, resources, funding acquisition, investigation, project administration, writing. LG: resources, funding acquisition, methodology. AE: methodology, writing, investigation.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

CJH acknowledges the support of the Czech Science Foundation (20-26767Y). Charles University Centre of Advanced Materials (CUCAM) (OP VVV Excellent Research Teams, project number CZ.02.1.01/0.0/0.0/15_003/0000417) is acknowledged along with e-INFRA CZ (ID: 90254) from the Ministry of Education, Youth and Sports of the Czech Republic. CJH and LG acknowledge support from the Charles University Centre of Excellence award UNCE/SCI/014. LG acknowledges the Primus Research Program of Charles University (PRIMUS/20/SCI/004). The authors thank Jakub Szmitek for helpful discussions and analysis regarding free energy simulations.

References

1. N. Felvey, J. Guo, R. Rana, L. Xu, S. R. Bare and B. C. Gates, et al., Interconversion of Atomically Dispersed Platinum Cations and Platinum Clusters in Zeolite ZSM-5 and Formation of Platinum gem-Dicarbonyls, J. Am. Chem. Soc., 2022, 144(30), 13874–13887.
2. P. Serna, A. Rodríguez-Fernández, S. Yacob, C. Kliewer, M. Moliner and A. Corma, Single-Site vs. Cluster Catalysis in High Temperature Oxidations, Angew. Chem., Int. Ed., 2021, 60(29), 15954–15962.
3. L. Liu, U. Díaz, R. Arenal, G. Agostini, P. Concepción and A. Corma, Generation of subnanometric platinum with high stability during transformation of a 2D zeolite into 3D, Nat. Mater., 2017, 16(1), 132–138.
4. M. Moliner, J. E. Gabay, C. E. Kliewer, R. T. Carr, J. Guzman and G. L. Casty, et al., Reversible Transformation of Pt Nanoparticles into Single Atoms inside High-Silica Chabazite Zeolite, J. Am. Chem. Soc., 2016, 138(48), 15743–15750.
5. J. D. Kistler, N. Chotigkrai, P. Xu, B. Enderle, P. Praserthdam and C.-Y. Chen, et al., A Single-Site Platinum CO Oxidation Catalyst in Zeolite KLTL: Microscopic and Spectroscopic Determination of the Locations of the Platinum Atoms, Angew. Chem., Int. Ed., 2014, 53(34), 8904–8907.
6. J. Van der Mynsbrugge, M. Head-Gordon and A. T. Bell, Computational modeling predicts the stability of both Pd+ and Pd2+ ion-exchanged into H-CHA, J. Mater. Chem. A, 2021, 9(4), 2161–2174.
7. H. A. Aljama, M. Head-Gordon and A. T. Bell, Assessing the stability of Pd-exchanged sites in zeolites with the aid of a high throughput quantum chemistry workflow, Nat. Commun., 2022, 13(1), 2910.
8. D. Hou, L. Grajciar, P. Nachtigall and C. J. Heard, Origin of the Unusual Stability of Zeolite-Encapsulated Sub-Nanometer Platinum, ACS Catal., 2020, 10(19), 11057–11068.
9. L. Liu and A. Corma, Metal Catalysts for Heterogeneous Catalysis: From Single Atoms to Nanoclusters and Nanoparticles, Chem. Rev., 2018, 118(10), 4981–5079.
10. K. Morgan, A. Goguet and C. Hardacre, Metal Redispersion Strategies for Recycling of Supported Metal Catalysts: A Perspective, ACS Catal., 2015, 5(6), 3430–3445.
11. L. Chen, T. V. W. Janssens, P. N. R. Vennestrøm, J. Jansson, M. Skoglundh and H. Grönbeck, A Complete Multisite Reaction Mechanism for Low-Temperature NH3-SCR over Cu-CHA, ACS Catal., 2020, 10(10), 5646–5656.
12. K. Wettergren, F. F. Schweinberger, D. Deiana, C. J. Ridge, A. S. Crampton and M. D. Rötzer, et al., High Sintering Resistance of Size-Selected Platinum Cluster Catalysts by Suppressed Ostwald Ripening, Nano Lett., 2014, 14(10), 5803–5809.
13. L. Liu, M. Lopez-Haro, C. W. Lopes, C. Li, P. Concepcion and L. Simonelli, et al., Regioselective generation and reactivity control of subnanometric platinum clusters in zeolites for high-temperature catalysis, Nat. Mater., 2019, 18(8), 866–873.
14. L. Liu, M. Lopez-Haro, J. A. Perez-Omil, M. Boronat, J. J. Calvino and A. Corma, Direct assessment of confinement effect in zeolite-encapsulated subnanometric metal species, Nat. Commun., 2022, 13(1), 821.
15. L. Liu, D. N. Zakharov, R. Arenal, P. Concepcion, E. A. Stach and A. Corma, Evolution and stabilization of subnanometric metal species in confined space by in situ TEM, Nat. Commun., 2018, 9(1), 574.
16. R. Ferrando and A. Fortunelli, Diffusion of adatoms and small clusters on magnesium oxide surfaces, J. Phys.: Condens. Matter, 2009, 21(26), 264001.
17. L. Xu, G. Henkelman, C. T. Campbell and H. Jónsson, Small Pd Clusters, up to the Tetramer At Least, Are Highly Mobile on the MgO(100) Surface, Phys. Rev. Lett., 2005, 95(14), 146103.
18. P. Cnudde, E. A. Redekop, W. Dai, N. G. Porcaro, M. Waroquier and S. Bordiga, et al., Experimental and Theoretical Evidence for the Promotional Effect of Acid Sites on the Diffusion of Alkenes through Small-Pore Zeolites, Angew. Chem., Int. Ed., 2021, 60(18), 10016–10022.
19. K. T. Schütt, H. E. Sauceda, P. J. Kindermans, A. Tkatchenko and K. R. Müller, SchNet – A deep learning architecture for molecules and materials, J. Chem. Phys., 2018, 148(24), 241722.
20. J. Behler, First Principles Neural Network Potentials for Reactive Simulations of Large Molecular and Condensed Systems, Angew. Chem., Int. Ed., 2017, 56(42), 12828–12840.
21. A. Erlebach, P. Nachtigall and L. Grajciar, Accurate large-scale simulations of siliceous zeolites by neural network potentials, npj Comput. Mater., 2022, 8(1), 174.
22. Y. Yu, W. Zhang and D. Mei, Artificial Neural Network Potential for Encapsulated Platinum Clusters in MOF-808, J. Phys. Chem. C, 2022, 126(2), 1204–1214.
23. S. Ma and Z.-P. Liu, Zeolite-confined subnanometric PtSn mimicking mortise-and-tenon joinery for catalytic propane dehydrogenation, Nat. Commun., 2022, 13(1), 2716.
24. R. Millan, E. Bello-Jurado, M. Moliner, M. Boronat and R. Gomez-Bombarelli, Effect of Framework Composition and NH3 on the Diffusion of Cu+ in Cu-CHA Catalysts Predicted by Machine Learning Accelerated Molecular Dynamics, ACS Cent. Sci., 2023, 9(11), 2044–2056.
25. O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky and K. T. Schütt, et al., Machine Learning Force Fields, Chem. Rev., 2021, 121(16), 10142–10186.
26. V. Fung, J. Zhang, E. Juarez and B. G. Sumpter, Benchmarking graph neural networks for materials chemistry, npj Comput. Mater., 2021, 7(1), 84.
27. P. C. Jennings and R. L. Johnston, Structures of small Ti- and V-doped Pt clusters: A GA-DFT study, Comput. Theor. Chem., 2013, 1021, 91–100.
28. F. Gao, Y. Wang, N. M. Washton, M. Kollár, J. Szanyi and C. H. F. Peden, Effects of Alkali and Alkaline Earth Cocations on the Activity and Hydrothermal Stability of Cu/SSZ-13 NH3−SCR Catalysts, ACS Catal., 2015, 5(11), 6780–6791.
29. J. Song, Y. Wang, E. D. Walter, N. M. Washton, D. Mei and L. Kovarik, et al., Toward Rational Design of Cu/SSZ-13 Selective Catalytic Reduction Catalysts: Implications from Atomic-Level Understanding of Hydrothermal Stability, ACS Catal., 2017, 7(12), 8214–8227.
30. L. E. Burris and M. C. G. Juenger, Effect of calcination on the reactivity of natural clinoptilolite zeolites used as supplementary cementitious materials, Constr. Build. Mater., 2020, 258, 119988.
31. X. Guo, L. Qiao, S. Zong, R. Ye, Y. He and J. Cheng, et al., Effect of NaY Zeolite at Different Calcination Temperatures on the Activity in Hydroformylation of Formaldehyde, ChemistrySelect, 2022, 7(36), e202201574.
32. K. T. Schütt, P. Kessel, M. Gastegger, K. A. Nicoli, A. Tkatchenko and K. R. Müller, SchNetPack: A Deep Learning Toolbox For Atomistic Systems, J. Chem. Theory Comput., 2019, 15(1), 448–455.
33. A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen and M. Dułak, et al., The atomic simulation environment—a Python library for working with atoms, J. Phys.: Condens. Matter, 2017, 29(27), 273002.
34. S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys., 1984, 81(1), 511–519.
35. W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A, 1985, 31(3), 1695–1697.
36. M. Bonomi, D. Branduardi, G. Bussi, C. Camilloni, D. Provasi and P. Raiteri, et al., PLUMED: A portable plugin for free-energy calculations with molecular dynamics, Comput. Phys. Commun., 2009, 180(10), 1961–1972.
37. M. Bonomi, G. Bussi, C. Camilloni, G. A. Tribello, P. Banáš and A. Barducci, et al., Promoting transparency and reproducibility in enhanced molecular simulations, Nat. Methods, 2019, 16(8), 670–673.
38. G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi, PLUMED 2: New feathers for an old bird, Comput. Phys. Commun., 2014, 185(2), 604–613.
39. J. Kästner, Umbrella sampling, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1(6), 932–942.
40. G. Kresse and J. Hafner, Ab initio molecular-dynamics simulation of the liquid-metal–amorphous-semiconductor transition in germanium, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49(20), 14251–14269.
41. G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47(1), 558–561.
42. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized Gradient Approximation Made Simple [Phys. Rev. Lett. 77, 3865 (1996)], Phys. Rev. Lett., 1997, 78(7), 1396.
43. S. Grimme, S. Ehrlich and L. Goerigk, Effect of the damping function in dispersion corrected density functional theory, J. Comput. Chem., 2011, 32(7), 1456–1465.
44. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, J. Chem. Phys., 2010, 132(15), 154104.
45. G. Henkelman, B. P. Uberuaga and H. Jónsson, A climbing image nudged elastic band method for finding saddle points and minimum energy paths, J. Chem. Phys., 2000, 113(22), 9901–9904.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4nr00017j

---