Fluorine spillover for ceria- vs silica-supported palladium nanoparticles: A MD study using machine learning potentials

It has long been speculated that catalyst supports often participate actively in heterogeneous catalysis rather than acting as inert bystanders, merely providing a high surface area material on which catalysts are dispersed. One of the ways that support can participate in the reaction is through the spillover of adsorbate species from the metal catalyst to the support.^1–3 Theoretical studies of spillover are currently being pursued, but mostly for carbon supports.^4–6 To date, studies of metal oxide supports are more limited.^5,7 This is due to the inherent complexity of the systems, where the structure of the oxide support, the interface between the catalyst and support, and the shape and structure of the catalyst nanoparticles all provide challenges for modeling. In addition, it is necessary to realistically describe the adsorption and diffusion of chemisorbed reactant species in these complex systems.

Metallic nanoparticles supported on oxides, as well as adsorption and reaction on these nanoparticles, provide a distinctive challenge for computational studies.⁸ Their sizes are often too large for direct ab initio calculations but too small for more phenomenologically coarse-grained continuum modeling. In addition, the need to treat a combination of transition metals, oxides, and chemical reactions makes the development of effective empirical potentials difficult. Machine learning naturally provides an enticing new approach to the problem, and there are some new developments in applying machine learning to the general problem of nanoparticles, supported or unsupported.^9–11 It should, however, be noted that other phenomena such as morphological evolution for larger metal nanoparticles on time scales up to minutes or longer, are best addressed via stochastic modeling, ideally incorporating a realistic description of surface diffusion kinetics.^12–14 

Recent advances in the development of machine learning potentials^15–19 make molecular dynamics (MD) simulations feasible for up to thousands of atoms and for up to a few nanoseconds. In this contribution, we will demonstrate that simulations at a size and time scale achievable through machine learning potential can facilitate answering many questions for heterogeneous catalysis on supported metal nanoparticles. Using fluorine adsorption on ceria and silica supported Pd nanoparticles as an example, we develop simulation methods that can be used to study the structure of the oxide support (Sec. III), size and shapes of the supported catalyst nanoparticles (Sec. IV), and adsorption and migration of reactants on supported nanoparticles as well as spillover to the supports (Sec. V). The general methodology for using DeePMD to train ML potential from DFT data is given in Sec. II, and the conclusion is given in Sec. VI.

Some aspects of the methodology and code development have already been reported.^20,21 Plane-wave density-functional theory (DFT) calculations were performed using the VASP code²² with Perdew–Burke–Ernzerhof (PBE)²³ functionals and the recommended PAW potentials released with VASP 5.4. The energy cutoff is 400 eV for all elements. The extraction of energy and forces and training of an artificial neural network with the se_e2_a descriptor were performed using the DeePMD-kit (v2.1.5) software.²⁴ MD simulations and data analysis were performed with Python, heavily relying on the atomistic simulation environment (ASE).²⁵ MD simulations are mostly performed using constant NVT (canonical ensemble), and the temperature is controlled by a fluctuating Langevin force and equilibrated by a friction force with a coefficient of 0.002. The time step is 2 fs for systems where the lightest element is O and 5 fs for systems where the lightest element is F.

Through trial and error, we have formulated some guidelines for the efficient development of machine learning potentials from ab initio data for systems involving nanoclusters and surfaces. We train individual potentials for all combinations of relevant elements. For example, for fluorine adsorption on ceria supported palladium particles, potentials for Ce, F, O, Pd, CeF, CeO, CePd, FO, … , and CeFOPd are all developed instead of just a single potential for CeFOPd. This is because it is much more efficient to perform calculations using the potential tailored to individual components, and they are generally more accurate with the same amount of training. Although, strictly speaking, potentials developed for individual components do not give exactly the same result as potentials for more complete systems, they are generally close enough not to cause any physical problems. For example, structure minimized using the potential for smaller systems will also be stable using the potential for the larger system. This is also related to the feature that the library of DFT data is shared between all components, so that DFT results for Pd are also used for training potential for PdO, CePd, CeOPd, etc.

The generation of DFT data is an iterative process. The initial configuration for training is generally a nanocluster consisting of the relevant elements. In cases where there is a common and stable bulk phase for those elements, such as fcc for the single element Pd, we choose it as the starting point. Otherwise, we search for the densest phase in the Material Project.²⁶ Nanoclusters with side lengths around 1.3 nm are placed in a cubic supercell with twice the volume so that they are allowed to expand and/or melt at 3000 K, simulated for 2000 time steps using the VASP code. Note that we use a single Γ k-point and so tend to choose supercells of at least 1.5 nm in side length. The results are then used for training the initial ML potential. Using Pd as an example, we first create a cubic Pd nanocluster with 108 atoms and place it in a cubic supercell with a 1.5 nm cell length. We perform DFT MD simulations at 3000 K for 2000 time steps, or 10 ps. The initial ML potential is obtained by training these DFT data for 10⁵ training steps. We then create another system of 256 Pd atoms with the fcc structure in a supercell with periodic boundary condition and an equilibrium volume of 15.8 nm³. The system is then simulated using the ML potential at 4000 K, while gradually increasing the supercell size from 1.58 to 1.6 nm for another 2000 time steps (10 ps) (this slight increase in volume is optional; it is used here for ease of systematically adjusting the volume when the need for additional densities arises). The resulting configuration is then used as the starting point for another DFT simulation at the same temperature. The new DFT data, combined with the old data, is used to train a new ML potential. We use this ML potential to simulate the 256 Pd atom systems at 4000 K but this time compress it to a supercell with a length of 1.4 nm, so that density is about 1.43 times the equilibrium value. We then simulate this denser structure with DFT and combine the new data to train a new ML potential. We find that it is important to vary the density of the system to obtain a ML potential with general applicability, not just for the equilibrium bulk system. Even though, as a whole system, the palladium system will not reach a density as high as 1.4 times the equilibrium volume under normal conditions, the ability to model correctly much shorter bond lengths than the equilibrium value is important for describing defects at step edges and kink sites.

In general, after three iterations, a model with no obvious pathologies can be obtained. Under conditions of normal temperature and pressure, the ML potential does not give a catastrophically wrong answer. The results can be considered generally acceptable for instructive purposes. However, for more serious applications, further refinements and additional validation are often necessary. For a specific application of the ML potential, we perform DFT calculations first on a system of reasonable size achievable by DFT, generally starting with a configuration obtained by simulations with ML potentials. We then compared the DFT results with the ML potential. For example, simulated annealing of Pd nanoparticles of 192 atoms was performed using the ML potentials, and configurations below and near the melting temperature (1000 and 1500 K) are used for the initial configuration of further DFT simulations. If the rms errors for energy comparing DFT and ML potentials is too big (e.g., larger than 0.05 eV/atom), further training of the ML potential are conducted using the new DFT data. In the case of Pd, we find the rms error is only 0.01 eV, even without further training. Most systems, however, do require additional training and further iterations.

The methodology for systematic validation of the ML potentials is still in development. At this stage, much of the validation relies on our own physical and chemical knowledge and intuition. We inspect manually the results of MD simulations using ML potentials and compare them with existing knowledge of the system. When a system behaves “strangely,” it usually means that the ML potential is exploring a part of the PES where it has not been trained. Additional DFT calculations are then performed, the results of which are used for further refinement of the ML potential.

For a more quantitative assessment, we need to compare ML potentials with those of DFT calculations. Here, we use a two-pronged approach. The first is an application specific approach. We use ML potentials to create systems that are expected to be encountered during the specific application, such as pure Pd, ceria, and silica nanoparticles, ceria and silica supported Pd, and fluorine adsorbed on those nanoparticles. They are of sufficiently small sizes so that DFT calculations can be performed, and the results are compared with the predictions of the ML potentials.

Table I shows error estimates of two ML potentials used in this work for these selected systems. The RMSE for the total energy (per atom) and forces on each atom is obtained from 30 randomly selected configurations for each DFT MD simulation.

The accuracy of the ML potentials for these scenarios is the most important factor specific to our applications. However, it may not be an accurate gauge of their general applicability. In fact, the choice of the training dataset inevitably reflects the desired application due to the nature of our iterative procedure, so the training is biased toward these applications. To provide a more unbiased assessment, we conduct more systematic tests. First, we randomly fill up a simulation box with atoms of each element. The system is equilibrated at 1600 K for 2000 time steps between each addition of atoms by performing MD simulations using the ML potentials. The resulting configurations are then used as the starting point for DFT MD simulations for another 200–2000 time steps. The DFT trajectories are then compared with ML potentials. The equilibration using the ML potentials vs purely randomly filling the cell with atoms serves two important considerations: first, it allows us to validate the most energetically favorable part of the ML PES; second, assuming the ML potential is reliable, it allows us to focus on the physically meaningful part of the PES. There is no need to model the unphysical part of the PES very accurately.

Table II summarizes the RMSE in energy/atoms and forces for various randomly synthesized systems using the above described procedures. As expected, the errors for random systems are generally larger than those for application-specific systems. There are some variances, such as the fact that agreement for pure Pd systems is much better than for Ce, perhaps reflecting the fact that there is more extensive DFT data for Pd systems than for Ce. Fortunately, it is in line with our application since we expect to see pure Pd nanoparticles, but not pure Ce nanoparticles in our system. For systems that do involve pure Ce, naturally, more DFT data for Ce should be generated.

A more physical illustration of the validity of the ML potential can be represented by the comparison of the compositional formation energy between DFT and the ML potential. The results for selected pairs of elements are shown in the supplementary material.

Ceria (CeO₂) and silica (SiO₂) are two common catalyst supports. The crystal structures of ceria and silica are well-studied and can be used to design the model surfaces for catalyst support. For example, it is known from DFT calculations that the CeO₂(111) surface is the most stable.^27,28 Thus, it is natural to study, e.g., catalytic properties, using the ideal CeO₂(111) surface as the basis but potentially also introducing defects such as oxygen vacancies using energy derived from DFT calculations²⁹ or experiments.³⁰ 

Here, we develop a methodology for preparing metal oxide surfaces using MD simulations with ML potentials that could be widely applicable to a variety of systems. For ceria, we start with a nanocube of Ce₆₀₀O₁₂₀₀ and raise the temperature to 3000 K at a rate of 25 K/ps using the ML potential, so the precise starting configuration is not so relevant. We then put the nanoparticle in a 3.8 × 3.8 × 7.6 nm³ supercell and gradually increased the supercell to a size of 4.8 × 4.8 × 9.6 nm³ at 3000 K in 10 ps. The nanoparticle dimension is sufficiently close to 3.8 nm that at this temperature and time scale, there is enough mobility for the ceria nanoparticle to bridge the simulation cell and, given the periodic boundary conditions, convert it into a near flat slab. We then reduce the system temperature at various cooling rates to 800 K while maintaining the same supercell.

Figure 1 shows the results of such a simulated annealing process for ceria. With a faster cooling rate (above 1 K/ps), the system is not very ordered, with only a short range hexagonal structure present. The surface consists of many single-atom-wide ridges. With a slower cooling rate, the system becomes more ordered, with smooth hexagonal surfaces. At the slowest cooling rate, a step structure is also present on the surface. The energy vs temperature plot has a clear signature of a first-order transition corresponding to crystallization around 1200 K with a slower cooling rate. The fact that this is much lower than the usually cited 2700–2900 K bulk melting temperature for ceria is expected to be associated with effects of both small sizes and hysteresis in simulations and does not necessarily indicate any deficiency of the ML potentials. For comparison, MD simulations of a homogeneous ceria system (Ce₁₀₈O₂₁₆ in a cube with 16.40 Å side length and periodic boundary conditions) yield a melting temperature of around 2800 K at a heating rate of 2.5 K/ps.

Silica of the type used as catalyst support is less amenable to standard DFT analysis since the material is generally amorphous. Using a similar method as for ceria, Fig. 2 shows results for simulated annealing of SiO₂. A slower cooling rate leads to a more stable structure. However, even at 0.5 K/ps, no ordered structure is found. There are signatures for latent heat in the E − cT vs T plots, but the transition is much further spread out and does not become obviously sharper as the cooling rate increases. In addition, the energy at lower temperatures is not as sensitive to the cooling rate, at least for the range simulated. Those are indicative of the formation of amorphous silica, which is consistent with the observations that amorphous and glassy silica are common products of solidifying silica melt or from the sol-gel process, while crystallization requires more carefully controlled cooling. MD simulations have long been used to study amorphous silica, but mostly using empirical potentials.³¹ However, recently, machine learning potential has started to compete with empirical potential.³² 

After preparing model ceria and silica surfaces, we can now model supported Pd nanoparticles. One could consider the following two approaches to describe this process: The simpler one is to directly place a preformed Pd nanoparticle on the surface of the support, then study the evolution of its shape and structure due to interaction with the support using MD simulations. The other approach is to simulate impregnation of the support with Pd precursors via incipient wetness followed by calcination and reduction, as often utilized in experiments.²¹ 

While the second approach is beyond our current capabilities, we have developed a methodology that mimics the formation of metal nanoparticles (albeit on much shorter time scales) and can, thus, be used to systematically explore this process. We assume individual free Pd atoms spontaneously emerge through reduction during the formation process. For Pd/CeO₂, we use the model system obtained in Sec. III (also see Fig. 1) by cooling to 800 K with the slowest cooling rate. Pd atoms are randomly placed in the 4.8 × 4.8 × 9.6 nm³ simulation cell, with the condition that each new atom is at least 0.3 nm away from any existing atoms. One of the parameters that can be adjusted is the rate at which Pd atoms are added to the simulations. In principle, this should be a stochastic process that occurs randomly throughout a simulation. However, for simplicity, in this work, we just add a fixed number of Pd atoms for every 4 ps of simulations.

Figure 3 shows the results of simulations with different Pd deposition rates after 1.024 ns of simulations at 800 K. Most of the Pd islands have a semi-spherical shape, but some of them are triangular or trapezoidal. Those more exotic shapes are due to the merging of smaller islands and represent transient configurations during a reshaping process.

Figure 4 shows the results for the formation of Pd islands on a silica support. A smaller 3.8 × 3.8 × 7.6 nm³ cell is used here, containing approximately the same number of atoms as the ceria. We observe a higher mobility of Pd nanoparticles compared with ceria. This feature is consistent with the observation that the shape of Pd nanoparticles is more spherical. The smaller contact area with the support indicates weaker binding between Pd and silica than for ceria, which correlates with higher mobility.

Before we present the results of fluorine adsorption on the supported Pd catalyst nanoparticles, we briefly discuss the results of fluorine adsorption separately on each of the components (on Pd, on ceria, and silica). We perform MD simulations using the same method as in Sec. IV for Pd, but now instead using F atoms. This is not intended to simulate any particular reaction but rather to provide a generic picture of F adsorption.

The prepared model Pd slab has a (111) surface orientation, with a Pd island on one side of the slab. At low F dosages, F adatoms occupy three-fold sites. At the elevated temperature of 800 K, no preference for step edges is observed. As the coverage of F increases, the surface consists of a large semi-ordered Pd–F structure and smaller Pd–F complexes that are scattered on the lower terrace [see Fig. 5(a)]. At even higher dosages, palladium fluoride grows in height and forms 3D islands.

For fluorine adsorption on ceria, we observe that F preferentially adsorbs at the step edge at low coverage. Some of the Ce atoms at kink sites can have two F atoms attached. As F coverage increases, more and more Ce atoms detach from the step edges and have multiple F atoms attached to them. Those CeF_x units can also form small Ce–F complexes, but no large, ordered surface cerium fluoride is observed. At much larger F coverage, CeF_x molecules start to detach from the surface (at 800 K).

Fluorine adsorption on the model amorphous silica is much harder to characterize, with no obvious stages. It can be generally described as the competition of fluorine with oxygen to form bonds with silicon atoms. Unlike Si and O, Si and F do not form complex networks. Finally, at very high F coverage, SiF_x molecules also start to detach from the surface.

Some quantitative analysis can be performed on the MD simulation results. By measuring the change in potential energy as the number of F atoms increases, we can deduce the heat of adsorption and use it as the adsorption energy for fluorine adatoms. Using the energy (per atom) of a F₂ dimer in the gas phase as the reference, the fluorine adsorption energy on a clean surface is estimated to be −2.4 eV for Pd, −2.5 eV for ceria, and −1.5 eV for silica. A higher level of negative energy means a stronger preference for F to attach to the surface. Note that on a flat CeO₂(111) surface, the adsorption energy is only −1.6 eV from our MD simulations, so defect sites are crucial for understanding the catalytic properties of ceria supports. On the other hand, we also studied fluorine adsorption on a α-quartz(0001) surface with no defects and no dangling bonds. A similar adsorption energy of −1.5 eV is obtained for low F coverage as for the model silica slab analyzed above. The F adsorption energy on a flat Pd(111) surface is −1.9 eV from MD simulations.

Using a ceria-supported Pd catalyst obtained with the methodology described in Sec. IV, we simulate F adsorption as in Sec. V A. Figure 6 shows the heat of adsorption as a function of the number of F atoms as they are incrementally added to the simulation. As the ceria support is generally free of defects, initially most of the F adatoms are on the Pd nanoparticles. In addition, note that the initial heat of adsorption has a value indicating much stronger adsorption for F than either Pd or ceria surfaces obtained in the section above. This is likely due to the presence of some very small Pd particles on the ceria support. As the number of F atoms increases, one starts to see some Ce atoms that are free of oxygen being lifted out of the surface around the Pd-ceria interface. As more of the Ce atoms are reduced, the oxygen density on the surface increases. Eventually, these oxygen atoms start to migrate toward the Pd nanoparticles and are ultimately chemisorbed on the Pd surface.

The hill-and-valley feature of the differential heat of the adsorption curve is sensitive to the F deposition rate used in the simulations. The peak tends to shift to the left and flatten out as the deposition rate decreases. Thus, the feature indicates some activated or slow processes that delay the system from finding the more energetically favorable configurations. Here, we tentatively identify the rate limiting process for fluorine spillover from the Pd nanoparticle to the ceria support as being the migration of oxygen from the ceria to the Pd. We are currently developing methodologies to derive activation energy from MD simulations, which should further clarify this issue.

The results for F deposition on a silica-supported Pd catalyst model is shown are Fig. 7. The differential heat of adsorption initially assumes the same value as for F/Pd and becomes less negative as the F coverage increases due to repulsion between chemisorbed F adatoms, eventually reaching the value for F/silica obtained in Sec. V A. There is no big drop in the curve as in the case of Pd/ceria (cf. Fig. 6), consistent with the observation that no significant migration of oxygen from silica to Pd is observed. It is worth mentioning that when the number of F atoms is sufficiently high, F atoms do adsorb onto the silica surface, but only after Pd is fully covered by F. Thus, we conclude from this MD study that the spillover of F from Pd to silica is not significant.

Molecular dynamic simulations using machine learning potential trained by DFT data enable us to provide a truly dynamic picture of the migration of reactants during heterogeneous catalysis on supported nanoparticles, an important phenomenon that is difficult to study theoretically with traditional methods. There are several novel applications of machine learning potential here: (a) preparing model surfaces using simulated annealing as catalyst supports; (b) mimicking the formation of catalytic metal nanoparticles on model supports; and (c) adsorption on catalyst surfaces and its study using the differential heat of adsorption. The latter constitutes an extension of an earlier study of S adsorption on coinage metal surfaces.²⁰ 

Part of the motivation for this study is a previous experimental study²¹ that shows that a ceria-supported Pd catalyst has an exceptional ability to remove fluorine from fluorophenols (compared to Pd catalysts on other supports). Fluorine spillover from Pd to ceria support was conjectured from some preliminary MD studies to be one of the factors enabling this enhanced ability. The presence of fluorine in ceria has been verified experimentally.²¹ In this study, we use a more realistic crystalline ceria substrate rather than the amorphous nanocluster used in the previous study, thus more realistically modeling the experiments. We suggest that the rate limiting process for spillover is likely to be the reverse migration of oxygen from the ceria support to the Pd nanoparticles.

The existence of potential complex synergetic effects involving supported metal nanoparticles, as well as the supports themselves, in catalytic processes means that traditional methods of separating complex systems into different components for analysis might not capture key aspects of behavior. In the study of fluorine adsorption separately on Pd and ceria, no significant differences were found. Worse yet, using CeO₂(111) as the substrate, even weaker F adsorption is found than on Pd. Only through MD studies using a ML potential, which can access time scales on the order of nanoseconds, is it revealed that the migration of fluorine in this system is a result of the interplay between multiple elements. The ability of Pd to partially reduce ceria is crucial, especially at the later stages of fluorine adsorption. By systematically comparing ceria and silica for their roles in fluorine adsorption, we also understand why this is a unique feature of ceria, despite the fact that fluorine is known to also reduce silica.

Insights obtained from studies using the ML potential also identify the limitations and pitfalls of relying only on DFT calculations for modeling catalysis. Since the CeO₂(111) surface is shown to be the most stable surface and is also observed on ceria nanocubes in experiments, it is natural to use a CeO₂(111) slab as a model surface for DFT studies. As mentioned earlier, the F adsorption energy is much weaker on a flat Ce₂(111) surface than on a Pd(111) surface, so a standard DFT study will not reveal any F spillover on ceria supported Pd. One can systematically introduce defects such as steps, kinks, and even some amorphousness to the ceria model. However, one also needs to account for the energetics of the defects in addition to the many different ways that fluorine can adsorb those defects. Such a study is beyond the capability of typical DFT calculations, while machine learning potential is most suitable for this type of study.

A recent experimental study³³ of Ir clusters on CeO₂(111) thin films grown on a Cu(111) substrate is consistent with many of the findings of this work, including the stability of nanoparticles and reverse oxygen spillover effects. It provides further indication that the use of ML potentials derived from DFT calculations can be a reliable method to study heterogeneous catalysis on supported nanoparticles.

See the supplementary material for lists of systems used for training the potentials, comparisons between DFT and the ML potentials for select pairs of elements (Ce–O, Si–O, and F–Pd), potential files for Ce–F–O–Pd and F–O–Pd–Si systems, and an input file for training.

This work was supported by the US Department of Energy, Office of Science, Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biological Sciences, and was performed in the Computational and Theoretical Chemistry (CTC) project at Ames National Laboratory, which is operated by Iowa State University under Contract No. DE-AC02-07CH11358.

The authors have no conflicts to disclose.

Da-Jiang Liu: Conceptualization (equal); Formal analysis (equal); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). James W. Evans: Conceptualization (equal); Funding acquisition (lead); Investigation (supporting); Methodology (supporting); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (supporting).

The data that support the findings of this study are available within the article and its supplementary material.