Multistep CO2 Activation and Dissociation Mechanisms on PdxPt4–x Clusters in the Gas Phase

Palladium, platinum, and their alloys are promising catalysts for electrochemical CO2 reduction reactions (CO2RR), leading to the design of durable and efficient catalysts for the production of useful chemicals in a more sustainable way. However, a deep understanding of the CO2RR mechanisms is still challenging because of the complexity and the factors influencing the system. The purpose of this study is to investigate at the atomic scale the first steps of the CO2RR, CO2 activation and dissociation mechanisms on PdxPt4–x clusters in the gas phase. To do it, we use Density Functional Theory (DFT)-based reaction path and ab initio molecular dynamics (AIMD) computations. Our research focuses on the description of CO2 activation and dissociation processes via the computation of multistep reaction paths, providing insights into the site and binding mode dependent reactivity. Detailed understanding of the CO2–cluster interaction mechanisms and estimating of the reaction energy barriers facilitate comprehension of why and how catalysts are poisoned and identification of the most stable activated adducts configurations. We show that increasing the platinum content induces fluxional behavior of the cluster structure and biases CO2 dissociation; in fact, our computations unveiled several dissociated CO2 isomers that are very stable and that there are various isomerization processes leading to a dissociated structure (possibly a CO poisoned state) from an intactly bound CO2 one (activated state). On the basis of the comparison of the PdxPt4–x reaction paths, we can observe the promising catalytic activity of Pd3Pt in the studied context. Not only does this cluster composition favor CO2 activation against dissociation (thereby expected to facilitate the hydrogenation reactions of CO2), the potential energy surface (PES) is very flat among activated CO2 isomers.

Figure S1: Tetranuclear starting structures for the clusters' geometry optimization in the gas phase. W, X, Y, Z = Pd/Pt atom.
of the optimized structures as where E mol (P d x P t 4−x ) is the SCF energy of a cluster composed by x Pd-atoms and 4 − x Ptatoms, x = 1, ..., 4, E(P t), E(P d) are the SCF energies of the Pd and Pt atoms, respectively.
First we show the difference of the structures' atomization energy depending on the basis set. Each of the Figures S2 -S12, show the results for one of the starting conformations in FigureS1 and spin multiplicities 1,3, or 5. From this results, we observe that the def2-SVP and LANL2DZ basis sets yield atomization energies higher than the def2-TZVP basis set.
For the further computations only the tetrahedron like structures were used, since they were the most stable ones. Even if we started the optimization from another structure, the final optimized cluster structure was often found to exihibit rhomboidal and tetrahedral.
For each starting geometry, we found that the triplet structure is the one with most negative atomization energy (i.e., the most stable structure, energetically). Figures S13-S15 show the atomization energies of the optimized structures for singlet, triplet and quintet spin multiplicities of the clusters, respectively. Even if the initial geometry was different, most of the structures converged to a tetrahedron-like shape, so the energies of the optimized structures are very similar. As stated in the main text, the structures with tetrahedron-like geometry and triplet state are the ones with lowest energy.

Pd x Pt 4-x +CO 2 structures
We generated CO 2 cluster adducts considering six binding modes, displayed in Figure S16 (intact CO 2 ) and Figure S17 (dissociated CO 2 ). Different basis-sets have been used for the bare-cluster structure optimization, as mentioned in Section Computational Methods in the main text.

Benchmarking of Pt 4 geometry and spin multiplicity
We analyzed the stabilizes of the quadrangular and tetrahedron-like structures with spin multiplicities 3 and 5 using various Density Functionals, comparing to high level CCSD(T) computations.
We optimized the structure using the def2-TZVPP, def2-QZVP and also the def2-QZVPP basis sets in conjunction with the PBE functional, including D3 dispersion correction ( Figure   S23). Subsequently, to test how the choice of the exchange.correlation functional influences the results, we re-optimized the cluster structures with the def2-QZVPP basis set and the B3LYP, TPSS, TPSSh functionals ( Figure S24). We report in the following the values of the bond lenghts and angles for the Pt 4 isomers. As last step of our analysis, we computed the CCSD(T)/def2-QZVPP single-point energies of the structures optimized by PBE-D3/def2-TZVP method. The results are also shown in Figure S24. Our results infer that the triplet multiplicity, tetrahedron-like structure has lower energy then the quintet state quadrangular structure. Furthermore, we observe that the description given by the TPSSh functional is the one that is closer to the CC energy. We utilize this knowledge for compute with TPSSh+D3/def2-QZVPP level of theory the single point energies of the single cluster and CO 2 structures, the highest energy-barrier structure and the most stable product in the main text, see Figure S31.
We use the T 1 diagnostics to estimate the possible extent of multireference character of the Pt 4 molecules. Table 1 shows the results for the computations. All the values we obtained are larger than 0.17 (quadrangular quintet). Typically, T 1 values smaller than 0.02 indicate that single reference methods perform well 3 . However, it has been reported for 4d transition metal clusters and their oxides, hydrides, and metal halides that the 0.02 value is not practical for 4d TM-containing species and propose T 1 < 0.045 4 . Yet, our calculations are far from this value. These indicates a possible multireference charachter of all the computed structures.  Figure S24: Test of Pt 4 with different functionals.
The aim of our test was to examine whether the cluster keeps the triplet ground state upon CO 2 adsorption. We tested the starting structure of the intact and dissociated CO 2 binding modes; our computations showed that for the CO 2 binding modes the triplet structures have lower energy than the singlet and quintet adducts. For the dissociated-CO 2 binding modes, we found a few exceptions, however those few structures have very high energy and we excluded them from the list potential CO 2 activation or dissociation reaction intermediates.  Figure S27: Test of the energy of the geometry-optimized adducts for all the CO 2 binding modes. The triplet structures have lower energy.

S25
Testing the accuracy of the PBE-D3/def2-TZVP method for the CO 2 binding energies We tested the accuracy of the applied PBE-D3/def2-TZVP method compared to coupledcluster reference values for the CO 2 binding energies of the first intermediate along the reaction paths. Please note that as it is described in the main text, the T 2 1 diagnostics indicates possible multi-reference characters of the cluster's electronic structure. Indeed, it has been shown earlier by multi-reference complete active space self-consistent field computations that there are 44 and 51 electronic states within 2.2eV and 1.2 eV above the ground state in the case of Pd 4 and Pt 4 , respectively 5 , and the ground and the first excited states are quasi-degenerate. Thus, apart from the conventional CCSD(T) computations (where the CC equations were solved using the Cholesky decomposition using the default settings in the Q-Chem code), we also used the CCSD(fT) method (which is equivalent to the CR − CCSD(T ) 2 ) 6 , what in several cases can describe also strongly correlated electronic structure. Nevertheless, it is important to note, that still these computations are not of benchmark quality, but are expected to provide correct tendencies. It is interesting to note that the CCSD(T) and the CCSD(fT) computations yield similar results.
The results above show that while the coupled cluster and the PBE-D3 DFT results are quantitatively different by up to 43 kJ/mol, the trends are correct, which supports the dopant-dependent binding tendencies shown in our manuscript.

Mayer Bond Orders
To comparison to the adducts' Mayer Bond Order, in figure S28 we present the Mayer bond orders of the cluster structures without CO 2 . Figure S28: Mayer bond orders of the triplet tetrahedron structures for each cluster composition.

EDA
In this section, we provide additional information from the EDA computations. We compared the results of the charge transfer from the TM cluster to the CO 2 moiety in the cluster-intact S27 CO 2 adducts.
In Figure S29b the charge transfer energy is related to the total dipole moment of the TM clusters in the gas phase. We observed that with higher TM dipole moment, the more structures we found. This means that the number of CO 2 binding sites on a cluster is related to the total dipole moment of the cluster itself.
We display in Figure S29a the mean charge-transfer energy and the energy of the HOMO, LUMO levels of the bare clusters. The values of the HOMO-LUMO gap are very similar (| − 0.2| ± 0.02 a.u.) for Pd 3 Pt, Pd 2 Pt 2 ,PdPt 3 , whereas the mean charge transfer goes from -510 to -320 kJ/mol. In Figure S30, we display the sum of the interaction and relaxation energies as function of the charge transfer. We cannot find any relationship between these two quantities, as mentioned in the main text. On the one hand, the relaxation energy is a quantity depending on the geometry of the adduct and on the geometry of the optimized free fragment, whereas charge transfer and interaction energy only depend on the final adduct.
On the other hand, the interaction energy is also depending on the preparation energy, which has a large contributions in adducts with dissociated CO 2 . As a result, we do not find the adsorption energy by summing interaction and relaxation energy in figure S30.

Direct CO 2 dissociation from the reaction paths
In Figure S31, we compare the relative energies of the highest lying transition structures, and the lowest-energy product structures to that of the separated clusters and CO 2 . We cannot find a linear scaling of the relative energy of the dissociation barrier with the Ptdoping as reported in 7 . In order to be more accurate when comparing the three energies for each cluster, the single point energies of the structures in Figure

Further AIMD simulations analysis
In section in the main text, we reported the characteristics of the processes observed in the AIMD simulations for each cluster composition. In the following, we describe two different analysis we performed on the trajectories. The first one is on the M-C distance (M: the nearest cluster atom to the carbon atom of CO 2 at the begin of the simulation), the second one a development of a Markov State Model for the evolution of the OCO angle.

M-C distance as function of time
It is difficult to find a good reaction coordinate (RC) for the AIMD simulations. In Figure   S32(a-e), we represent the distance between M-C as function of time. The distance oscillates in a range of 1 angstrom around an average value, except in the case of Pd 2 Pt 2 . The Pd 2 Pt 2 simulation shows periodically an abrupt change of the M-C distance, so it represents the umbrella motion we described in the main text.
However, for the other cluster compositions, the difference is not very useful as reaction coordinate. For instance, in the Pd 3 Pt simulation, we observe that the cluster opens into a bent-rhomboidal shape. But we do not see it in the M-C distance, because the atoms that move to open the tetrahedral shape are not the ones, we are tracking the distance of.
A simular results is also obtained for the Pt 4 or the PdPt 3 simulations, for which the M-C distance does not show that the cluster isomerizes.
Our analysis shows that the M-C distance does not represent a descriptive RC. This is first because the CO 2 molecule migrates around the cluster surface and the cluster also modify its structure, but we are tracking only the distance between the nearest metal atom at the start of the simulation to C. All these considerations made us dispense with using the M-C distance as RC for describing the AIMD simulations.

Markov State Model of the OCO angle
The variation of the OCO angle as function of time is reported in Figure S33. To analyze the variation of the OCO angle, we developed a Markov State Model (MSM) by utilizing K-Means for the partitioning the state space. The results did not show any trends for the CO 2 activation or differentiation of this as function of the metal cluster composition, as the OCO angle fluctuates around a 20-degrees range. The script has been written by the authors and the Kmeans algorithm we use is from the scikit-learn Python package 8 . We define the state space, i.e., all the possible values in our system (state is in the mathematical sense) to be the possible values that the OCO angle can take. The OCO angle as function of time is our trajectory. This is how our MSM model works: • Consider a trajectory. In this case, the trajectory is the OCO angle as function of time.
• Assign N centroids using the Kmeans algorithm to discretize the state space into Voronoi cells centered in the centroids.
• Assign each point of the trajectory to a centroid.
• Build a count matrix for counting the transitions from a Voronoi cell to another.
• Row-normalize the count matrix to obtain a transition probability matrix P See also a very clear explanation of MSM in the work of Prinz et al. 9 . Please, note that in Figure S33 the color sequence is not relevant, it is a mere coloring for distinguishing the different kmeans clusters. Time step = 5 fs. Figure S33 shows that, for each cluster composition, the OCO angle of the CO 2 molecule opens and closes for the whole simulation. The only exception is for the simulation Pd 3 Pt simulation, in which we see that the CO 2 angle is larger in the second half of the simulation.

S33
The corresponding transition probability matrices are: We thought that the partitioning of the state space (the values of the OCO angle) and the computation of the transition matrix could give some information about the CO 2 behavior.
Instead, we see a cyclic behavior between the states: from state 1 (central value), one can go to state 2 or 3 (which is, close or open the angle), and transitions from state 2 to state 3 are not allowed (OCO angle does not change of 50 degrees in 5 fs). The model shows what happens in the simulation (the OCO angle varies in a certain range and does not change abruptly), which means that the model represents well the process. However, the information we obtain is almost obvious.  Figure S33: OCO angle as function of time for each P d x P t 4−x cluster. The color coding in each figure represents the assignment of the trajectory point to one of three centroids by K-means algorithm.

Reduction of the oxide clusters
Reduction of the oxide clusters was investigated by computing the reaction energy (∆E) of The negative reaction energies indicate thermodynamically favoured reduction of the oxide clusters.