On the Stability of Cu5 Catalysts in Air Using Multireference Perturbation Theory.

An ab initio study of the interaction of O2, the most abundant radical and oxidant species in the atmosphere, with a Cu5 cluster, a new generation atomic metal catalyst, is presented. The open-shell nature of the reactant species is properly accounted for by using the multireference perturbation theory, allowing the experimentally confirmed resistivity of Cu5 clusters toward oxidation to be investigated. Approximate reaction pathways for the transition from physisorption to chemisorption are calculated for the interaction of O2 with quasi-iso-energetic trapezoidal planar and trigonal bipyramidal structures. Within the multireference approach, the transition barrier for O2 activation can be interpreted as an avoided crossing between adiabatic states (neutral and ionic), which provides new insights into the charge-transfer process and gives better estimates for this hard to localize and therefore often neglected first intermediate state. For Cu5 arranged in a bipyramidal structure, the O-O bond cleavage is confirmed as the rate-determining step. However, for planar Cu5, the high energy barrier for O2 activation, related to a very pronounced avoided crossing when going from physisorption to chemisorption, determines the reactivity in this case.

Introduction

Metal clusters of subnanometer size, composed of a few atoms only, have emerged as a new generation of catalysts[1−3] and photocatalysts[4] with appealing properties arising from their molecule-like electronic structures. The synthesis of such quantum clusters made of just a few metal atoms has been achieved by kinetic control using an electrochemical technique without employing any surfactants or capping agents.[5,6] This technique has opened the possibility for numerous applications ranging from cancer-therapeutic drugs[5,7] to efficient hydrogen photoproduction.[8]

Triggered by the recent interest in subnanometric metal clusters as heterogeneous catalysts, we address the question of their stability under oxygen atmosphere via a high-level ab initio theory. For this purpose, we have selected the particular case of the Cu5 cluster. This choice is motivated by a recent work[9] demonstrating via ab initio modeling and selected experiments that the deposition of a single monolayer of highly stable Cu5 clusters onto a TiO2 surface makes it an innovative visible-light photo-active material. Much more energy can be harvested from sunlight, and the coated titanium dioxide stores this energy temporarily in the form of charge pairs—electrons and holes—which is a perfect prerequisite for follow-up chemistry. Moreover, experimental measurements[10] have indicated a particularly high stability of Cu5 clusters with respect to oxidation up to a temperature of 423 K. The experiment by Corma’s group[10] was realized for copper clusters, with and without water, considering N-doped graphene as the support in XPS spectroscopic measurements. Our previous work has shown that a Cu5 cluster is minimally perturbed when supported on graphene due to the dispersion-dominated nature of the Cu5–graphene interaction. Therefore, the experimental results[10] can be considered as representative of unsupported Cu5 clusters. Our selection of the Cu5 cluster is also motivated by very recent theoretical results[11] indicating that TiO2-supported Cu5 clusters might allow for CO2 activation through sunlight, as well as a spontaneous decomposition, leading to CO desorption.

The first step in the oxidation process is interpreted here in the picture of an avoided crossing between two electronic states. In general, charge-transfer states can often be considered as precursors for irreversible oxidation reactions. In the current case, charge-transfer is interpreted as switching from a diabatic state which is asymptotically correlated to neutral Cu5 and O2 fragments, to another diabatic state which leads to a cationic Cu5+ and an anionic O2– species.

The occurrence of an efficient hopping between neutral and ionic states has been demonstrated recently for the Cs2–C60 reaction in helium droplets used as cryogenic matrices.[12−14] In the latter case, which is paradigmatic for harpoon-type reactions, the crossing takes place at a large distance, and the energy difference, ΔE between the asymptotes of neutral and ionic reactant species [ΔE = Eionic(∞) – Eneutral(∞)] is very small (ca. 0.7 eV). Assuming that, at large distances, the energy of the neutral (covalent) state is constant, and that the energy of the ionic state has a Coulomb-like behavior, with e being the charge of an electron, the distance R at which the crossing occurs can be estimated viaif atomic units are used for the energies and the distance.

However, with increasing energy difference at the asymptotic limit ΔE, the position R at which the crossing occurs, shifts to shorter distances. For instance, in the case of O2 interacting with a reduced TiO2 surface (ΔE > 2 eV), high level ab initio theory has shown that the crossing is located close to the physisorption minimum.[15,16] Eventually, as opposed to a harpoon scenario occurring in the Cs2–C60 reaction, if ΔE is large enough, the crossing becomes located at the repulsive region of the electronic state, which correlates to the neutral species[17] (see also Section S1 of the Supporting Information). With an energy difference of ΔE = 7.4 eV as obtained in the calculations described below, this situation can be anticipated for the O2–Cu5 interaction. As a result, an energetic barrier appears for the hopping process between the two relevant electronic states, which are asymptotically correlated to either neutral or ionic reactants. An important question in this context is how much the probability for electron hopping, and hence for the onset of full oxidation, is influenced by the temperature. On one hand, it is clear that the fraction of suitable O2 molecules with energies above the barrier increases with higher temperatures. On the other hand, higher relative velocities reduce the hopping probability because less time is spent in the crossing region. Therefore, the final outcome is determined by a subtle balance of these two opposite effects as a function of temperature. In order to get estimates of the probability of switching between the two relevant electronic states as a function of temperature, we apply the Landau–Zener (LZ) model,[18,19] a first-order approximation suitable for small[12] as well as large[20] values of electronic couplings.

Computationally, we approach this problem as follows. Typically, reaction pathways are evaluated, for example, via the nudged elastic band algorithm[21] or a string method[22] on the high-dimensional potential energy surface of the system at hand. In the current case, it is necessary to determine a path leading from the asymptotic or “reactant” region of that surface, that is, the O2 (X3Σg–) and Cu5 (X~2A1) species, to a short-range region where the activated complex can form. Unfortunately, as both reactants are open-shell species in their corresponding ground electronic state, a costly multireference treatment is required to adequately capture their interaction. Note that in such cases, the dominant electronic configuration of the activated complex is inherently different from those of the reactant species. Moreover, the combination between a doublet (Cu5) and a triplet (O2) state gives rise to doublet and quartet states, which become degenerate in the asymptotic region. Hence, both doublet and quartet manifolds of the complete system need to be considered. The computational expense of a wavefunction-based multireference treatment forces us to approximate the actual reaction pathway by a partial freezing of less-involved nuclear degrees of freedom in our study. In other words, we provide convenient cuts through the actual PES and restrict our investigation to structural configurations of higher symmetry. Both a planar trapezoidal and a trigonal bipyramidal structure of the bare Cu5 cluster have to be taken into consideration due to their very similar energies. Preliminary explorations of the energy landscape have shown that the most favorable orientations for the stabilization of the ionic state encompass high symmetry (C2) approaches for both planar and bipyramidal Cu5 structures. This is expected for the interaction between two symmetric and homonuclear species such as O2 (D2) and Cu5 (D3/C2) and the high directional nature of the main orbitals responsible for the O2–Cu5 binding, having allowed us to exploit the C2 symmetry in accelerating the computations.

Single-reference methods based on density functional theory (DFT) have been previously applied to systematic studies of the O2–Cu interaction (n ≤ 38).[10,23] These studies have provided valuable insights into the increased stability of copper clusters in O2 when the copper cluster size decreases. Besides the probability of formation of O2–Cu5 complexes in the charge-transfer chemisorption states, the barrier for O2 dissociation into two individual O atoms attached to the copper cluster was also considered and found to be the rate limiting state for full oxidation. It is clear, however, that the applicability of single-reference methods, such as DFT, to the charge-transfer reactions between two open-shell systems such as a copper cluster and oxygen molecule, remains rather questionable, especially for low-spin states. As one illustrative example, a previous study of the O2 photodesorption from a reduced TiO2– surface highlighted the high multiconfigurational character of the wavefunction for the O2–TiO2– interaction, with TiO2– having being modeled as a cluster in a doublet spin state.[15,16] The analysis of the open-shell Cu5 cluster in this work extends these previous DFT studies[10,23] by considering the multireference nature of the O2–Cu5 interaction. Similarities, but also qualitative differences will be discussed in the following sections.

Methods

The size of the system still permits a multireference treatment of the O2 (X3Σg–) and Cu5 (X̃2A1) interaction in the two relevant electronic states. We focus on the adsorption of molecular oxygen at a bridging site of the Cu5 cluster in the planar trapezoidal as well as the bipyramidal structure. In both cases, the bridging site corresponds to the global energy minimum for the physisorption of O2. For details of the optimization procedure, please see the Supporting Information. Imposing C2 symmetry, the geometry of the Cu5 cluster in the planar structure only depends on four variables. As shown in Figure , these parameters are rCu1, the distance between the two Cu atoms facing the oxygen molecule; rCu2, the distance between the two Cu atoms on each side of the trapezoidal structure; θ, the angle between rCu1 and rCu2; and finally rh, the distance between the geometric center of rCu1 and the central atom of the planar Cu5 cluster. The geometry of the trigonal bipyramidal Cu5 cluster is determined by rCu1, the distance between two Cu atoms in equatorial positions, and rCu2, the distance between a Cu atom in equatorial and in axial position.

Figure 1

Coordinates of the O2–Cu5 system for planar trapezoidal (left-hand panel) and trigonal bipyramidal (right-hand panel) Cu5 structures.

The intermolecular distance between the O2 center-of-mass and the geometric center of rCu1 (denoted as d) can be interpreted as a geometrically simplified approximation to the reaction pathway at the onset of the oxidation process.

First, the geometries of the O2 and Cu5 fragments are optimized separately. Then, these geometries are kept frozen and potential energy curves (PECs) are calculated as a function of the O2–Cu5 distance d (see Figure ), preserving the C2 symmetry.

We use the polarized correlation-consistent triple-ζ basis of Dunning and collaborators[24] (cc-pVTZ) for oxygen atoms, whereas a double-ζ basis set was employed for copper,[25] including a small-core (10-valence-electron) relativistic pseudopotential. These electronic structure calculations are performed with the MOLPRO program package.[26] Additional calculations (presented in the Supporting Information) were performed with the ORCA[27] suite of programs (version 4.0.1.2).

Multireference Treatment

We perform multireference calculations in order to obtain the PECs of the lowest-energy states. Because of the open-shell character of the electronic wavefunctions such calculations are not straightforward. In particular, when describing the long-range region, the standard orbital optimization algorithms typically fail to retrieve physically meaningful orbital occupations. To solve this issue and to ensure correct orbital occupation numbers, we apply the following strategy: the orbitals of the Cu5 cluster are first optimized while the O2 fragment is treated as a dummy partner using ghost orbitals. Then, the orbitals of O2 are optimized while considering the Cu5 cluster as the dummy partner instead. After merging these orbitals, the Hartree–Fock electronic wavefunction is calculated in the asymptotic region with the correct orbital occupation numbers. These orbitals are then stored and used as an initial guess in the follow-up calculations using the state-averaged complete-active-space self-consistent (CASSCF) treatment. In order to keep the calculations computationally feasible, we chose the minimum active CAS space able to describe the charge transfer between O2 and Cu5 for both planar and trapezoidal Cu5 structures, considering only the geometries of C2 symmetry. For the case of the planar trapezoidal cluster, the orbitals labeled as 22a1, 9b1, 18b2, and 8a2 were included in the active space, whereas the 21a1, 9b1, 16b2, and 7a2 orbitals were kept doubly occupied. For the trigonal bipyramidal cluster, the 22a1, 13b1, 14b2, and 8a2 orbitals were active, whereas the 21a1, 13b1, 12b2, and 7a2 orbitals were kept doubly occupied. All labels refer to irreducible representations in the C2 molecular point group. Both choices of active spaces yield ten configuration state functions (CFSs) from eight doublet and two quartet spin states. For both cases, only electronic B2 states were considered in the state average CAS calculations, with B2 also corresponding to the symmetry of the ground state of O2. As the combination of doublet (Cu5) and triplet (O2) states results in doublet and quartet states, these two spin states have been considered for both planar and bipyramidal isomers, exploring several binding sites (bridge and on-top) and various orientations of the O2 reactive species.

To capture the dynamical correlation effects, we further employ the multistate complete-active-space second-order perturbation theory (MS-CASPT2). For this purpose, the reference CFSs and molecular orbitals obtained in the previous CASSCF calculations are used. In a second step, the rCu1, rCu2, θ, rh, and rO–O variables are optimized at each value of the O2–Cu5 distance, d. The calculations performed for the doublet spin states show that neutral and ionic states have different spacial symmetries and, thus, provide no probability for the hopping process.[28] Hence, our study can be restricted to the quartet spin states.

LZ Model

To further analyze hopping probabilities between neutral and ionic states as a function of temperature, we apply the LZ model,[18,19] a well-known, one-dimensional semiclassical model, which provides reasonable estimates of probabilities PLZ for non-adiabatic transitions via the approximationwith v as the relative velocity of the fragments and F12 as the difference between the two slopes, F1 and F2 of the diabatic PECs at the intersection between neutral and ionic states. The probability of hopping is defined as 1 – PLZ. In the above expression, H12 is the off-diagonal matrix element of the electronic Hamiltonian. Although the accuracy of the LZ formula depends on the values of H12, F12, and v, its validity under rather general conditions has been demonstrated (see e.g., ref (29)). As the purpose of the present work is to provide a qualitatively correct description of the molecular oxidation process, the employment of the LZ model can be considered sufficient.

Within the LZ model, half of the minimum energy splitting between the adiabatic states is identified with the value of H12 at the crossing region. In contrast to harpoon-type reactions as mentioned in the introduction, the location of the avoided crossing at the repulsive region is directly related to a large difference of the slopes of the neutral and ionic PECs at their intersection, which indicates reduced probabilities for a non-adiabatic transition. Assuming a Maxwell–Boltzmann (MB) distribution for the relative velocities of the reactants, the electron hopping probabilities can be written as a function of temperature as follows. We integrate over the hopping probabilities from the LZ model, PLZ, expressed as a function of the velocity, v in the reaction coordinate, and weighted with a Boltzmann factor, PMB,where PMB(v) denotes the MB distribution of relative velocities in one direction,with kB as the Boltzmann constant. Assuming that the Cu5 cluster is supported on a bulk surface of infinite mass, the reduced mass, μ is simply the mass of molecular oxygen. MB distributions were also assumed for calculating the fraction of O2–Cu5 pairs (referred to as f) with kinetic energies above the adiabatic energy barrier between neutral and ionic PECs. The global probability for the hopping process to occur can thus be estimated as Ptot ≈ (1 – PLZ) × f.

Results and Discussion

Planar Structure of Cu5

For all sites investigated in this work, the O2–Cu5 interaction is repulsive for the doublet spin state. This result indicates that the O2 molecule in its triplet state lacks a favorable interaction with the unpaired electron of Cu5. For the quartet state, the O2–Cu5 interaction is more favorable so that bonding occurs, with a partial charge transfer from Cu5 to O2. The lowest barrier to chemisorption is associated to a bridge site with the internuclear axis of O2 parallel to the plane of the copper cluster, as represented in Figure . The most favorable approach has a C2 symmetry because this orientation maximizes the overlap between the π* antibonding orbitals of O2 and the inner shell of d-orbitals of the Cu5 cluster, which are still well localized due to the subnanometer size of the cluster (see Section S6 of the Supporting Information). This favorable overlap greatly stabilizes the ionic state. When the energy of the ionic state starts approaching that of the ground state, the efficiency of the coupling increases, producing the energy barrier toward chemisorption and allowing a partial charge transfer from Cu5 to O2. At the transition state, the wavefunction presents a strong multireference character with weights of 80 and 20% for neutral and ionic states, respectively. The ionic character of the wavefunction increases progressively as O2 gets closer to Cu5. At the chemisorption minimum, the character of the wavefunction becomes essentially ionic. It is worth recalling here the importance of using a multireference method to properly describe the O2–Cu5 interaction and, in particular, to provide a correct representation of the wavefunction at the energy barrier. In the upper left panel of Figure , we analyze one-dimensional cuts through the potential energy surfaces of the O2–Cu5 reaction in the case of a planar cluster structure. The states are plotted as a function of the distance d between both fragments as defined in Figure , whereas all remaining nuclear degrees of freedom are kept frozen.

Figure 2

Left-hand panels: O2–Cu5 interaction energies in neutral and ionic states with Cu5 arranged in a planar trapezoidal (upper panel) or trigonal bipyramidal (bottom panel) structure. The reaction pathway is approximated by keeping the geometries of the reactant species frozen for each O2–Cu5 distance, d as defined in Figure . Right-hand panels: Hopping probability from neutral to ionic states at their crossing as a function of temperature, plotted together with the fraction of O2–Cu5 pairs with kinetic energy above the value of the adiabatic energy barrier. The highest temperature for which resistivity to oxidation has been reported from experimental measurements[10] (423 K) is shown as a gray dashed line.

Within the adiabatic representation (see Figure ), the X̃4B1 ground and the 2̃4B1 excited state are asymptotically correlated to a separation in neutral fragments or charged species Cu5+ and O2–, respectively. It can be observed that the two states feature a very pronounced avoided crossing at the repulsive region of the diabatic neutral state (dashed blue line), at about 2 Å. For larger distances, the reactant species conserve their neutral character and their interaction is dominated by van der Waals forces. It can be characterized as a molecular physisorption with a well-depth of −0.05 eV at a distance of 4.20 Å (see Table ). From Table , it can also be observed that zero-point energy (ZPE) contributions have a very modest influence (below 0.02 eV).

Table 1

Characteristics of the O2–Cu5 Reaction Along the Intermolecular Distance, d, with the O–O Bond Length and Cu5 Internal Coordinates Shown in Figure Fully Relaxeda

planar trapezoidal Cu5				trigonal bipyramidal Cu5
interaction		dmin, Å	Emin, eV	interaction		dmin, Å	Emin, eV
O2–Cu5	(neutral)	4.20	–0.05	O2–Cu5	(neutral)	3.60	–0.09
O2–Cu5	(ionic)	1.88	–1.40	O2–Cu5	(mixed ionic–covalent)	1.95	–0.12
			(0.02)				(0.01)
O2–Cu5	(barrier)	2.12	0.43	O2–Cu5	(barrier)	2.48	0.09
		[2.10]	[0.45]			[2.38]	[0.12]

Variables dmin and Emin denote the position and energy at both the potential minima and the barrier from the neutral to the ionic state (see also Figures and 3). Values in square brackets correspond to the frozen approach (see Figure ). ZPE corrections at the global minima are indicated in parentheses.

Figure 3

O2–Cu5 reaction energy pathway in the adiabatic ground state (X 4B1), with Cu5 arranged either in a planar (left-hand panel) or a bipyramidal structure (right-hand panel). The geometries of the reactant species have been relaxed at each O2–Cu5 distance d as defined in Figure .

After passing the avoided crossing (see Figure ), the interacting species access a region characterized by a strong Coulomb attractive interaction (−1.40 eV, see Table ) between the positively charged copper cluster and the negatively charged O2 species. This region is much better characterized when the geometries of the reactant species are allowed to relax. This can be seen in the left-hand panel of Figure (see also Table ), where the geometries of the O2 and Cu5 reactants are allowed to relax at each O2–Cu5 distance, d. In this case, the O–O bond is elongated until reaching a value close to the equilibrium distance for the superoxo O2– radical (about 1.31 Å). Notice also that the Cu–Cu bond closest to O2– is stretched by about 0.2 Å from the value for the bare Cu5 cluster. Once the O2–Cu5 system enters the well of the ionic state, it is highly improbable to turn back: it would be necessary to overcome an energy penalty of about 1.8 eV and to redistribute this energy in nuclear vibrational modes so that the geometries of the neutral fragments are recovered. On the other hand, the chemisorption geometry is just an intermediate state which is followed by the actual splitting of the O–O bond (i.e., the dissociation of O2 into two individual O atoms attached to the Cu5 cluster). However, the energy barrier for O2 splitting (0.12 eV, see Section S6 of the Supporting Information) in this case turns out to be lower than the barrier for state-to-state crossing (0.43 eV, see Table ). Therefore, it is clear that the transition to the precursor charge-transfer state as illustrated in Figure is the limiting step for an irreversible oxidation through the breaking of the O–O bond and the formation of new O–Cu bonds.

An estimate of the ratio τ between the reaction rates from left to right and right to left (see Figure ) can be obtained using the expression based on the Arrhenius equation,where Eminionic and Eminneutral stand for the energy at the potential minima of ionic and neutral states, respectively. Values of 2 × 10–23 and 8 × 10–17 are obtained for τ at 300 and 423 K, respectively. By definition, these ratios are equivalent to the equilibrium constants of the corresponding reaction. Their extremely small value indicates that oxidation is definitively irreversible at the experimentally relevant temperatures (see ref (10)). However, as will be shown below, despite being thermodynamically favored, oxidation cannot take place due to the high reaction barrier, which prevents the system from relaxation into a thermodynamic equilibrium. In other words, O2 molecules in the gas phase do not have enough kinetic energy in the reaction coordinate to overcome the necessary barrier, which renders the oxidation reaction thermodynamically allowed but kinetically forbidden.[30]

We further estimate the probability for surface hopping between neutral and ionic states via the LZ approximation. Note that the PECs determined by freezing all but one degree of freedom of molecular motion for the reactant species yields one-dimensional cuts of the full potential energy surface. Most importantly, only the partial freezing in a convenient higher symmetry makes a costly multireference approach actually feasible. Furthermore, it provides smooth PECs and allows for a better decoupling of the kinetic energy terms. The curves can be easily diabatized, allowing for an application of the one-dimensional LZ model. It is important, however, to calculate also the relaxed PECs in order to determine more precise values of the barrier heights and estimate how much the formed complex can be further stabilized. It can be seen from Table (by comparing unrelaxed and relaxed barrier heights) and from Figures and 3 that the effect of relaxation is rather small in the zone of the avoided crossing, and therefore barely affecting the hopping probabilities. It is more pronounced at shorter distances where the cluster relaxation is more significant.

As can be seen in the upper right panel of Figure , the LZ probabilities are very large for temperatures within the experimental range (between 293 and 423 K, see ref (10)). However, assuming a MB distribution as a function of temperature, the fraction of O2–Cu5 collisional pairs, f with energies above the threshold value which is necessary to overcome the barrier between neutral and ionic states, is negligible. As a result, access to the ionic state is not possible. The barely noticeable decrease of hopping probabilities with increasing temperature is a direct consequence of the very large value of the coupling term H12, which reflects the large energy gap between ground and excited adiabatic states in the avoided crossing region. Only at very high temperatures, the barrier at the avoided crossing is naturally overcome.

With this in mind, it is clear that the O2–Cu5 interaction is well characterized as an intermolecular physisorption at 423 K. However, considering a physisorption well-depth of 0.05 eV (see Table ), about 70% of physisorbed O2 molecules would have enough energy to escape from the physisorption minimum already at T = 373 K, if their translational and rotational degrees of freedom are accounted for (see Section S4 of the Supporting Information).

It is also worth recalling that a very different picture is obtained when applying the single-reference DFT-based theory (see Section S7 of the Supporting Information). The O2 splitting into two O atoms bound to the cluster is found as the rate-limiting step for oxidation instead.

Bipyramidal Structure of Cu5

Our results have revealed that the planar trapezoidal structure of bare Cu5 is only slightly energetically favored over the trigonal bipyramidal structure shown in Figure (see Section S2 of the Supporting Information). This suggests that the actual shape of the nanoparticle will be strongly dependent on the actual support and might also fluctuate significantly at higher temperatures. Hence, a more realistic study of Cu5 oxidation must consider both structures.

As for the planar trapezoidal Cu5 structure, the O2–Cu5 (bipyramidal) interaction is repulsive for the doublet spin state. The most favorable reaction pathway is found when the approaching molecular oxygen stays on the equatorial plane of the bipyramidal Cu5 geometry, as represented in Figure . In fact, this binding site provides an optimal overlap of the d-orbitals of Cu5 with the π* antibonding orbitals of O2 (see Section S6 of the Supporting Information). As for the trapezoidal Cu5 structure, the wavefunction at the region of the potential barrier presents a strong multireference character. Interestingly, in contrast to the case of the trapezoidal cluster, the strongly mixed ionic/neutral character of the wavefunction (35%/65%) is kept also at the chemisorption minimum. Keeping the O–O distance and the internal degrees of freedom of Cu5 frozen, we obtain curves as shown in Figure (bottom left-hand panel). The interaction energies are plotted as a function of the distance d as defined in Figure . By relaxing the geometries of the O2 and Cu5 reactants at each intermolecular O2–Cu5 distance, we obtain the values for the potential minima and for the energy barrier shown in Table . Note that the asymptote of the ionic state is significantly lower than that corresponding to the planar Cu5 structure. The avoided crossing is still located at the repulsive region of the physisorption potential, but the energy at its intersection with the ionic state is much smaller. As a result, the adiabatic energy barrier is much lower (about 0.09 eV, see Table ). Also, the avoided crossing is associated with a very large value of the electronic coupling H12 so that the probability for a non-adiabatic transition is negligible and the hopping process probability is mostly determined by the fraction of O2–Cu5 complexes with kinetic energy in the reaction coordinate above the adiabatic energy barrier (see lower right panel of Figure ).

When the reactant species are allowed to relax (see right-hand panel of Figure ), it can be observed that the potential minimum of the electronic state asymptotically correlating with the charged fragments is much less deep for the bipyramidal geometry than for the planar structure (−0.12 vs −1.40 eV, see Table ). A Mulliken analysis of the charge population indicates that the electron transfer from Cu5 to O2 is only partial in this case (about −0.2 |e|), so that the Cu–O bonding can be characterized as intermediate between ionic and covalent. Also, in contrast to the clear ionic nature of the O2–Cu5 (planar) interaction at the global minimum, the Cu5 structure is barely different from that of the isolated Cu5 cluster. The potential energy profile characterizes a reversible transition from physisorption to chemisorption. At 423 K, about 18% of the O2–Cu5 pairs would have kinetic energy above the adiabatic energy barrier in the collisional process (see bottom right-hand panel of Figure ). If the O2–Cu5 system decays on the potential minimum of the mixed ionic–covalent state, a barrier of 0.22 eV would have to be overcome in order to return to the physisorption state. In contrast to the case of planar O2–Cu5, there is no need to redistribute the kinetic energy in various vibrational modes because the compact Cu5 structure has not been distorted at the potential minimum of the mixed ionic–covalent state.

Another interesting difference between the planar and bipyramidal case is the multiconfigurational character of the wavefunction at the chemisorption minima, which is much more pronounced for the bipyramidal O2–Cu5 complex. However, we notice that in both cases, the mixing of ionic and neutral electronic configurations is significant at the barrier, which underlines the importance of multireference treatments.

Using eq to get an estimate of the ratio between reaction rates from left to right and right to left, values of ca. 0.33 and 0.45 are obtained at 300 and 423 K, respectively. The last value is ca. 6 × 1015 larger than for the planar Cu5 structure.

This similarity of chemisorbed and physisorbed geometries has an obvious consequence: for the bipyramidal structure, the O–O splitting remains the rate determining step (see Section S6 of the Supporting Information). This direct comparison highlights the very different nature of the oxidation reaction pathway for bipyramidal and planar structures at experimentally relevant temperatures.

Finally, we would like to mention that many other reaction pathways apart from those shown in Figure are possible, but the pathways discussed are those exhibiting the lowest activation energies for chemisorption, for both planar and bipyramidal Cu5 structures. Therefore, they are expected to be more relevant for the oxidation chemistry.

Conclusions

The onset of oxidation in Cu5 atomic clusters has been investigated computationally via multireference perturbation theory. Two isomers of similar energy have been taken into consideration. When the Cu5 cluster adopts a planar trapezoidal structure, a high-energy barrier (0.43 eV) arising from an avoided crossing between the electronic states correlating with neutral and ionic (O2– + Cu5+) fragments prevents the access to the precursor charge-transfer state for the oxidation of Cu5. Even at T = 1500 K, the probability to access the potential region where irreversible oxidation takes place is less than 10%. This region is characterized by a very deep potential minimum (−1.40 eV) due to the strong Coulomb attraction between ionic O2– and Cu5+ species. In this region, both O2 and Cu5 geometries become significantly distorted from those of separated Cu5 and O2 species. Within the experimentally relevant temperatures,[10] the O2–Cu5 (planar) interaction is almost completely characterized as a physisorption-type interaction with a well-depth of −0.05 eV.

A slightly deeper physisorption well was found when the Cu5 cluster adopts a trigonal bipyramidal structure. The energy barrier from a physisorption-type to a mixed ionic–covalent state is very low for this structure (0.09 eV), but so is the actual energy difference between the corresponding minima (0.03 eV). The necessary charge transfer from Cu5 to O2 is very small (−0.2 |e|), and so is the actual geometric distortion of the Cu5 structure. In other words, the activation of the oxygen molecule in the chemisorbed state is very modest, rendering this very stable cluster configuration rather nonreactive. In this case the O–O splitting remains the costly, rate-determining step for irreversible oxidation. These findings underline the very strong impact of the cluster geometry at the subnanometer scale.

Altogether, this work provides new insights into the nature of the bonding occurring between atomic open-shell metal clusters and molecular oxygen. A mixed covalent–ionic bond is formed, which is characterized by both hybridization of d-shell Cu5 orbitals with the antibonding π* orbitals of O2 and a partial charge transfer. Characterizing the nature of the O2–Cu5 binding mechanism is the first step toward the understanding of the catalytic properties of Cu5 clusters in air when supported on technological relevant materials. Work is currently in progress to analyze if the high stability of Cu5 in O2 is preserved when the clusters are supported on surfaces of TiO2.

From a methodological perspective, our study clearly points out the multireference character of the air-stability and oxidation process of novel subnanometer-sized materials such as open-shell Cu5 clusters, and more generally, of the activation barriers in relevant heterogeneous catalytic reactions on atomic metal clusters. Efforts in this direction have also included the development of multireference DFT-based approaches[31,32] (see ref (33) for a very recent review on multireference-based modeling in heterogeneous catalysis). From the perspective of method development, subnanometer catalysts can be considered as highly challenging benchmark systems due to the very recent availability of cutting-edge experimental measurements of their structure and properties.