Aggregation and support effects in the oxidation of fluxional atomic metal clusters. The paradigmatic Cu₅ case

Abstract

The recent development of new synthesis techniques has allowed the production of monodisperse metal clusters composed of a few atoms. Follow-up experimental spectroscopic characterization has indicated the stability of these atomic metal clusters (AMCs). Despite the common assumption that the occurrence of an irreversible oxidation becomes more likely as the cluster size decreases, its quenching and reversible nature has been experimentally identified in the particular case of Cu₅ clusters, making them paradigmatic. This work aims to address the influence of aggregation and the effects of a chemically inert carbon-based support on the oxidation of AMCs, considering the case of Cu₅ as a model system. For this purpose, we present an extended first-principles study of the oxidation of Cu₅–Cu₅ and circumpyrene-supported Cu₅, comparing it with that of unsupported Cu₅, and combine dispersion-corrected density-functionals, first principles thermochemistry, and ab initio molecular dynamics (AIMD) simulations within an adiabatic approach. Our results indicate that a molecular chemisorption/desorption model is sensible upon consideration of aggregation and support effects in such a way that the predicted (p–T)-phase diagrams do not differ significantly from those obtained for unsupported Cu₅. We also provide insights into the decoupling of the Cu₅–Cu₅ dimer into Cu₅ sub-units through activated fluxional rotational motion, upon heating, as well as the adsorption of multiple O₂ molecules at high oxygen gas pressures. Furthermore, numerical evidence shows the likelihood of a support-mediated mechanism leading to the dissociation of chemisorbed peroxo O₂²⁻ species, delivering states with very similar energies to those characterized by molecular chemisorption. A Boltzmann-weighted average of the free energies of formation is computed as well, coming up with a diagram of the dominant copper oxidation states as a function of temperature and oxygen gas pressure.

---

1 Introduction

This work is motivated by the impact that aggregation and the effects of a chemically inert support might have on the oxidation of atomic metal clusters (AMCs) composed of a few atoms. The recent development of new synthesis techniques has allowed the production of monodisperse AMCs such as, for instance, Cu₅ clusters.¹ Previous works have shown their catalytic applications in, e.g., CO oxidation,^2,3 CO₂ reduction,⁴ and methanation,⁵ as well as their ability to modify the properties of a relevant semiconductor (TiO₂), making it a visible photo-active material⁶ (see, e.g., ref. 7 and 8 for recent reviews). One major concern has been on the stability of Cu₅ clusters in an O₂-rich atmosphere as it is generally assumed that nanoparticles suffer from an increasing susceptibility to oxidation with decreasing size. In contrast to that assumption, recent works have indicated that Cu₅ clusters experience reversible oxidation,^7,9–13 making them paradigmatic.

The possibility of reversible oxidation of Cu₅ clusters has been pointed out through ex situ experimental measurements at temperatures below 423 K¹⁴ as well as theoretical studies considering the chemisorption of either single^9,11,14 or multiple¹² oxygen molecules into Cu₅ clusters in the two structures identified by first-principles theory: planar trapezoidal (2D) and trigonal bipyramidal (3D). Thus, for a planar structure, a multireference ab initio treatment¹¹ has indicated that a high energy barrier prevents the Cu₅–O₂ collisional pair to access the precursor charge-transfer state for irreversible oxidation of Cu₅. In contrast, for a bipyramidal structure, the potential energy landscape has revealed a reversible transition from physisorption to chemisorption states, with the energetically costly O–O splitting remaining as the rate determining step for irreversible oxidation to occur. The most recent experimental evidence has shown a reversible oxidation of surface-supported Cu₅ clusters at thermodynamic equilibrium using in situ X-ray absorption spectroscopy and X-ray photoelectron spectroscopy.¹³ Highly oriented pyrolytic graphite (HOPG) was considered as a chemically inert support, at different oxygen gas pressures, and in the temperature range going from room-temperature (RT) up to 773 K.¹³

In general, the interaction of molecular oxygen with metal and metal oxides such as TiO₂ rutile (110) has been a long-standing problem over the last few decades as extensively reviewed by Yates's^15–17 and Henderson's^18,19 groups from an experimental perspective. As shown in the earliest first-principles studies,^20–24 molecular oxygen can be adsorbed onto TiO₂ rutile (110) as either a superoxo O₂⁻ radical or a peroxo O₂²⁻ species, both in molecular and dissociative forms, existing in several (photo-)desorption channels as well. The existence of multiple stable charged species and reaction pathways becomes even more challenging in the case of AMCs due to their fluxional nature. In fact, when the size of a metal cluster is reduced to a very small number of atoms, the d-band of the metal splits into a network of molecular orbitals with the interconnections between different atoms having the length of a chemical bond (1–2 Å).⁷ The ‘floppy’ character of the resulting ACM geometries invokes the concept of structural fluxionality. This concept has been recently advocated to explain the catalytic activities of AMCs (see, e.g., ref. 4, 7, 25 and 26). The occurrence of fluxional dynamics has also been suggested from recent spectroscopic observations.²⁷ This characteristic might lead to the existence of many local minima, having very similar energies, when AMCs interact with multiple O₂ molecules. Such an energy landscape is implicitly associated with fluxional dynamics and wide amplitude AMC–AMC and AMC–O₂ motions.

Despite the attention that AMCs have recently attracted, little is known about the effects that their aggregation might have in the most fundamental process (i.e., their oxidation). Interestingly, previous work on TiO₂-supported atomic copper clusters has shown that oxidation might quench their aggregation.²⁸ Besides complementing previous research on the interaction of O₂ molecules with the Cu₅ monomer,^7,9–13 we have selected the unsupported Cu₅–Cu₅ dimer as a model system to explore the aggregation effects in the oxidation of AMCs. Next, the interaction of the circumpyrene-supported Cu₅ monomer with multiple O₂ molecules has been analyzed as well with the main purpose of providing insight into the role of a chemically inert support in the oxidation of a model AMC. As reported in ref. 6, graphene has very little impact on the electronic structure of supported Cu₅ due to the dispersion nature of the Cu₅–graphene interaction. This way, the net charge donation from the copper cluster to the support was found to be insignificant (less than 0.02 |e⁻|). In this work, the consideration of circumpyrene as a molecular model of a carbon-based surface has aided in the analysis of its influence on the interaction of supported Cu₅ with the environmental gas-phase O₂ molecules.

This article is structured as follows. In Section 2, we provide a detailed description of the computational methods. The results are presented in Section 3 as follows: first, we compare the stability of (Cu₅)–(O₂)_n and (Cu₅–Cu₅)–(O₂)_2n (n < 9) clusters as a function of temperature and oxygen gas pressure. Second, we analyze the stability of circumpyrene-supported Cu₅–(O₂)_n (n < 9) complexes, comparing it with that of unsupported Cu₅−_n (n < 9) complexes. Finally, our findings are summarized and conclusions are presented in Section 4.

2 Methods

2.1 Electronic structure calculations and ab initio molecular dynamics simulations

2.1.1 Structural models. In all calculations, a trigonal bipyramidal (3D) structure of the Cu₅ cluster has been considered for both unsupported and surface-supported Cu₅ and Cu₅–Cu₅ complexes. The main reason for this selection is that the 3D Cu₅ structure has been the one identified in recent experimental measurements.¹³ A carbon-based support has been modeled with the polycyclic aromatic hydrocarbon (PAH) circumpyrene (C₄₂H₁₆).²⁹ Using the 3 × 3 graphene supercell model reported in ref. 30, the Cu₅–graphene interaction energy differed by less than 10% from that calculated for the Cu₅–circumpyrene complex (ca. 1.2 eV). These periodic electronic structure calculations were carried out using the Vienna Ab initio Simulation Package (VASP 5.4.4).^31,32 In particular, electron–ion interactions were described by the projector augmented-wave method,^32,33 using Perdew–Burke–Ernzerhof (PBE) pseudo-potentials as implemented in the program. The electrons of the C(2s, 2p) and Cu(3d, 4s) orbitals were treated explicitly as valence electrons using plane wave basis sets with a kinetic energy cutoff large enough to obtain converged interaction energies. For the sake of comparison, a previous work has shown that the PAH circumcoronene (C₅₄H₁₈) is large enough to model the graphene–Au₁₂ interaction.³⁴

2.1.2 Density functional theory and ab initio calculations. Dispersion-corrected density functional theory (DFT) has been applied to study the adsorption of up to 8 O₂ molecules on bare Cu₅ and Cu₅–Cu₅ systems as well as circumpyrene-supported Cu₅ clusters. Given the well-proven performance of the DFT-D3(BJ) scheme^35,36 in describing the adsorption of small silver and copper clusters on titanium dioxide (TiO₂) surfaces,^6,37,38 we have applied the same ansatz, using the Perdew–Burke–Ernzerhof (PBE) density functional,³⁹ but replacing the D3(BJ) dispersion correction with the most recent D4's Grimes parameterization.⁴⁰ As mentioned in ref. 40, the D4 dispersion correction outperforms the D3 ansatz for reference metal-containing compounds. However, in the case of the Cu₅–(O₂)₃ complex (see ref. 12), the substitution of D3(BJ) by the D4 dispersion correction modified its optimized structure (e.g., the Cu–Cu bond length) very little (by less than 0.001 Å).

Structural optimizations and frequency calculations have been carried out using the atom-centered def2-TZVPP⁴¹ basis set for copper and oxygen atoms, while the def2-SVP basis set has been used for carbon and hydrogen atoms of circumpyrene. Binding energies E_bind of n and 2n O₂ molecules binding into the bare Cu₅ monomer and the Cu₅–Cu₅ dimer are defined as

E_bind = E_Cu5–(O2)n − E_Cu5 − n·E_O2 (1)

E_bind = E_{(Cu5–Cu5)–(O2)2n} − E_Cu5–Cu5 − 2n·E_O2, (2)

while, for the case of engagement of n O₂ molecules into circumpyrene-supported Cu₅, we have used the expression

E_bind = E_{(Cu5/support)–(O2)n} − E_Cu5/support − n·E_O2, (3)

These binding energies have been calculated using the def2-QZVPP basis set,⁴¹ including the Boys-Bernardi counterpoise correction.⁴² The reduction states of the adsorbed O₂ molecules (neutral O₂, O₂⁻, and O₂²⁻) as well as each copper metal atom [Cu(0), Cu(I), and Cu(II)] have been assigned from a correlated analysis of Mulliken charges,⁴³ atomic spin populations of the molecular orbitals using the Hirshfeld method,^44,45 and the values of the optimized O–O distances, comparing the latter with those corresponding to gas-phase O₂, O₂⁻, and O₂²⁻ species. The adequacy of our analysis has been assessed for the case of the Cu₅–(O₂)₃ complex by comparing it with that reported in ref. 12 through an analysis of the Mulliken charges and atomic spin populations of the natural orbitals obtained at a high-level of ab initio theory. All dispersion-corrected DFT-based calculations have been performed using the ORCA suite of programs^46–48 (version 5.0.1).

For assessment purposes, density fitting single-state multi-configurational self-consistent-field (DF-CASSCF) calculations have also been carried out, as implemented in the most updated version of the MOLPRO code.⁴⁹ Since the CASSCF wave-function is built on a basis of configuration state functions and separate calculations have been carried out for different spin states, no spin contamination has been identified. The MOLPRO implementation of the internally contracted RS2C method⁵⁰ has been applied on top of DF-CASSCF calculations to cover dynamical correlation effects. We have used the polarized correlation-consistent triple-ζ basis set of Dunning and collaborators⁵¹ (cc-pVTZ) for oxygen atoms and the cc-pVTZ-PP basis set for copper atoms,⁵² including a small (10-valence-electron) relativistic pseudopotential. For density fitting, the associated MP2FIT and JKFIT bases have been used in DF-CASSCF and R2SC calculations, respectively. For the sake of comparison, the domain-based pair natural orbital coupled-cluster approach DLPNO-CCSD(T)⁵³ has been applied as well, as implemented in the ORCA code (version 4.0.1.2).

2.1.3 Ab Initio molecular dynamics simulations. Ab Initio molecular dynamics (AIMD) simulations have been carried out using the ORCA MD module written by Brehm and collaborators.^54–56 Newton's equations of motion have been solved through numerical integration using the velocity Verlet algorithm within the canonical ensemble with time steps of 0.5 and 1.0 fs. The simulations have been accomplished using the Berendsen thermostat at temperatures of 673, 473, and 573 K for bare Cu₅–Cu₅, circumpyrene-supported Cu₅–(O₂)₄, and circumpyrene-supported Cu₅–Cu₅, respectively.

2.2 Equilibrium thermodynamics: thermodynamic potential

At a given temperature (T) and partial oxygen gas pressure (p), considering that we are dealing with an open system with an arbitrary number of gas-phase O₂ molecules, the relative stability of Cu₅–(O₂)_n (or (Cu₅–Cu₅)–(O₂)_2n) complexes is determined by calculating the thermodynamic potential ω^57–60 (referred to as the (Gibbs) free energy of formation in, e.g., ref. 58 and 59).

ω(T,μ_O2,n) = ΔE_F,corr(T) − T·S_Cu5–(O2)n(T) + T·S_Cu5(T) − n·[small mu, Greek, tilde]_O2(p,T) (4)

It should be stressed that the Cu₅ clusters are thus treated as immobilized on the support, therefore not contributing to the enthalpy and are coupled to a heat bath of temperature T and an infinite reservoir of O₂ gas-phase molecules at pressure p. Under these idealized conditions, the ω potential becomes minimal at thermodynamic equilibrium. ΔE_F,corr, the first term on the right-hand side of eqn (4), corresponds to the formation energy of the Cu₅–(O₂)_n complex and is defined as

ΔE_F,corr(T) = E_Cu5–(O2)n − E_Cu5 − nE_O2 + E_corr(T,n), (5)

where E_Cu5–(O2)n and E_Cu5 denote the electronic energies of the oxygen-covered and pure (surface-supported) cluster, respectively, and E_O2 is the electronic energy of molecular oxygen. E_corr introduces corrections to the internal energy such as zero-point energies and thermal vibrational contributions (i.e., arising from the population of excited vibrational states at a given temperature). The next term on the right-hand side of eqn (4) stands for a correction with respect to the entropy s_Cu5–(O2)n of the cluster. The last term on the right side of eqn (4) is the temperature and pressure-dependent part of the chemical potential of molecular oxygen expressed as:where the pressure enters through the ratio p/p₀, with the reference oxygen pressure p₀ set to 1 atm (ca. 1033 mbar). Note that μ_O2 = E_O2 + [small mu, Greek, tilde]_O2. This way, the T, p-independent contribution to the O₂ chemical potential, E_O2, is shifted to the ΔE_F,corr(T) term (eqn (5)). The change of enthalpy is given by Δh_O2 = h(p₀,T) − h(p₀,T = 0 K). For maximum accuracy, values for h_O2 and s_O2 from the NIST database have been used.^61,62 Next, the (p,T)-phase diagram has been created by determining the number n of adsorbed O₂ molecules which minimizes ω at a specified temperature and a given oxygen gas pressure. The thermodynamic potential ω can also be expressed as the difference between the Helmholtz free energies of the Cu₅–(O₂)_n complex and the Cu₅ cluster [F_Cu5–(O2)n–F_Cu5] minus n × μ_O2 as reported in ref. 58 and 59.

We emphasize that our methodological protocol to analyze the stability of supported copper clusters holding multiple adsorbed O₂ molecules can be extended to other relevant processes such as hydrogen storage onto AMCs as a function of temperature and hydrogen gas pressure.⁶³ For this purpose, we have provided a general code, delivering phase diagrams using the thermochemistry output of the ORCA package (see the ESI†).

3 Results and discussion

3.1 Oxidation of unsupported Cu₅ and Cu₅–Cu₅

Following previous work on individual Cu₅ clusters,^12,13 we start by assuming that the Cu₅ monomer and the Cu₅–Cu₅ dimer can form complexes with multiple gas-phase O₂ molecules and calculate the number of O₂ molecules that can be adsorbed depending on the oxygen gas pressure and temperature. As theoretically shown in ref. 11 and 12, we also consider that the energy barriers from physisorption to chemisorption states are very low (ca. 0.1 in ref. 11 and 12), while those associated with O₂ bond breaking are very high (of a few eV for the unsupported Cu₅–O₂ and Cu₅–(O₂)₃ complexes^11,12). This assumption holds true for the reaction pathway to chemisorption states exhibiting the lowest activation energies for O₂ adsoprtion into both planar and bipyramidal Cu₅ structures.¹¹

Table 1 and Fig. 1 (upper left-hand panel) show the values of the binding energies as a function of the number of the adsorbed O₂ molecules onto the Cu₅–Cu₅ dimer, considering singlet and triplet spin states. For the sake of comparison, they are also shown for the Cu₅ monomer in the lowest-energy doublet spin state (see also ref. 13). The bottom panels of Fig. 1 present the free energies of formation at ambient pressure (1 atm), while the most likely adsorption states as a function of oxygen gas pressure and temperature appear in the phase diagrams (upper panels of Fig. 2). The optimized structures for bare Cu₅ and Cu₅–Cu₅ along with the most stable Cu₅–(O₂)_n and (Cu₅–Cu₅)–(O₂)_n complexes are shown in the bottom panels of Fig. 2 instead.

Table 1 Binding energies as a function of the number of O₂ molecules adsorbed onto Cu₅ and Cu₅–Cu₅. The number of adsorbed neutral [n_O2], superoxo , and peroxo species is also indicated. The ‘*’ quotation indicates the dissociation of one peroxo species in two O²⁻ ions

	Cu5–(O2)n (doublet)		(Cu5–Cu5)–(O2)2n (singlet)		(Cu5–Cu5)–(O2)2n (triplet)
		E bind/eV		E bind/eV		E bind/eV
n = 1	0/1/0	−1.4	0/2/0	−0.4	2/0/0	−1.2
n = 2	0/2/0	−2.5	1/2/1	−4.1	1/2/1	−5.1
n = 3	0/3/0	−4.9	0/2/4	−8.1	0/2/4	−9.1
n = 4	0/4/0	−5.8	0/4/4	−9.2	0/4/4	−10.3
n = 5	2/3/0	−5.5	2/4/3	−9.6	4/3/3	−10.7
n = 6	1/5/0	−7.4	2/8/2	−11.0	1/7/4*	−13.8
n = 7	0/5/2	−8.9	0/10/4	−15.0	0/9/5	−16.2
n = 8	1/5/2	−9.0	2/8/6	−12.7	—	—

---

Fig. 1 Free energies of formation of Cu₅–(O₂)_n and (Cu₅–Cu₅)–(O₂)_n complexes. Upper panels: electronic binding energies of the clusters as a function of the number of O₂ molecules attached to Cu₅ (panel a) and (Cu₅–Cu₅) (panel b) complexes in doublet and singlet spin states, respectively. Bottom panels: free energies of formation (w potential) comparison at 1 atm for Cu₅–(O₂)_n (panel c) and (Cu₅–Cu₅)–(O₂)_2n (panel d) complexes. The color coding in solid lines follows dot colors in the predicted phase diagrams (see Fig. 2). The energy values corresponding to less stable complexes are plotted with dashed lines. The predicted phases can be retrieved from the convex hull of the lowest-energy curves.

Fig. 2 Phase diagrams of Cu₅–(O₂)_n and (Cu₅–Cu₅)–(O₂)_n complexes. Upper panels (a–c): (a) phase diagrams of Cu₅–(O₂)_n (doublet spin state), (b) (Cu₅–Cu₅)–(O₂)_2n (singlet spin state) and (c) (Cu₅–Cu₅)–(O₂)_2n (triplet spin state) complexes. The phase diagrams show the n values corresponding to the most stable adsorption states as a function of the oxygen gas pressure and temperature. The color coding is used to indicate the different (p − T) regions where the given Cu₅–(O₂)_n (or (Cu₅–Cu₅)–(O₂)_2n) complex is the most stable, with the corresponding n value having been specified (e.g., n = 7 in the green-colored region). Bottom panels (d–f): structures corresponding to the most stable adsorption states indicated in the phase diagrams for Cu₅–(O₂)_n (doublet spin state, panel d), (Cu₅–Cu₅)–(O₂)_2n (singlet spin state, panel e) and (c) (Cu₅–Cu₅)–(O₂)_2n (triplet spin state, panel f) complexes.

Let us first compare how differently the Cu₅–Cu₅ dimer and the Cu₅ monomer interact with the surrounding O₂ molecules. As can be observed in Table 1, binding energies per Cu₅ cluster are smaller for the Cu₅–Cu₅ dimer. However, the most stable complexes shown up in the phase diagram (see Fig. 2) have 3 and 7 adsorbed O₂ molecules per Cu₅ cluster for either the monomer or the dimer.

The optimized structures for both Cu₅–(O₂)_n and (Cu₅–Cu₅)–(O₂)_2n complexes show enhanced stability when the O₂ molecules attach to bridge Cu₅ positions. Depending on their number, the O₂ molecules are absorbed as neutral (O₂) or charged superoxo (O₂⁻), or peroxo (O₂²⁻) species (see Table 1). The copper clusters thus donate electron charge to the O₂ molecules but the entire system remains neutral. Upon stretching of the Cu–Cu distances, Cu₅ and Cu₅–Cu₅ adapt their shape to accommodate the charged O₂ species, featuring larger O–O bonds, at its bridge sites.

A clear difference between the Cu₅ and Cu₅–Cu₅ cases is on the major contribution of adsorbed peroxo O₂²⁻ species over superoxo O₂⁻ radicals for the latter. This outcome can be simply explained by considering that the O₂ molecule adsorbed in between the two Cu₅ clusters tends to elongate its O–O distance to benefit from the attractive interaction with both copper clusters. This enlargement might cause its dissociation even at 0 K, as shown in Table 1 for the (Cu₅–Cu₅)–(O₂)₁₂ complex in the singlet spin state. For the triplet state, the binding energies are slightly larger but hardly compensate for the singlet-triplet energy splitting on the bare Cu₅–Cu₅ dimer (see below).

3.1.1 On the Cu₅–Cu₅ dimer. As can be observed in Fig. 2, the optimized structure of the Cu₅–Cu₅ dimer in the lowest-energy state is characterized by two coupled bipyramidal structures, considering both singlet and triplet spin states. Local minima have also been found for the structure shown in the left-hand panel of Fig. 3, being 0.3 eV above the lowest-energy state. Further AIMD simulations reveal the rotational motion of the two bipyramids upon heating (see the right-hand panel of Fig. 3).

Fig. 3 Snapshots from AIMD simulations showing the rotational decoupling of a Cu₁₀ cluster into two Cu₅ bypiramidal 3D subunits upon increasing the temperature (up to 673 K in the given snapshot after 6 ps of simulation).

For assessment purposes, DF-CASSCF and R2SC calculations have been carried out for the structures corresponding to both the global and the local minimum identified at the DFT-D4 level. Using CAS active spaces of increasing size [from (4,4), with four electrons in four orbitals, to (10,10), with ten electrons in ten orbitals], it is found that the lowest-energy structure is that predicted via the DFT-D4 ansatz. Its re-optimization using the DF-CASSCF method modifies the Cu–Cu bond length by less than 0.1 Å on average.

The interaction energy between the two Cu₅ monomers is very large (−4.9, −5.2, and −5.4 eV at LCPNO-CCSD(T), R2SC, and DFT-D4 levels, respectively). This attractive interaction is mainly determined by the dynamical correlation contribution (by about 70%), as covered by the R2SC approach. The large values of Mulliken charges on each Cu atom (larger than 0.1 |e⁻|) suggest that strong electrostatic contributions favor the stabilization of the Cu₅–Cu₅ dimer. With values above 0.1 |e⁻|, the analysis of Mulliken charges on copper atoms reveals that the interaction is favored by electrostatic contributions. A detailed investigation into the reorganization of the Cu₅ electronic charge upon dimerization is a topic of future prospect.

At DFT-D4 and R2SC levels, the triplet state is 0.8 and 1.0 eV above the singlet spin state, respectively. The optimized structures are very similar (see Fig. 2), with the average Cu–Cu bond length being just 0.02 Å larger for the triplet spin state. This similarity can be understood by considering that the dynamical correlation contribution to the attractive Cu₅–Cu₅ interaction differs very little (by only 0.2 eV) in singlet and triplet spin states. The energy difference in the interaction energy arises from non-dynamical correlation effects and can be explained as follows: the single-electron occupied s-type orbitals of each Cu₅ cluster lead to the formation of bonding and anti-bonding orbitals when the two Cu₅ clusters come together. Two electrons with opposite spins occupy the resulting bonding orbital in the singlet spin state of the Cu₅–Cu₅ dimer. Contrarily, two electrons with the same spin cannot occupy the same orbital. This way, both bonding and anti-bonding orbitals become single occupied in the triplet state of the Cu₅–Cu₅ dimer, destabilizing the system.

Despite the strongly attractive nature of the Cu₅–Cu₅ interaction, our results indicate that a quenching of the Cu₅–Cu₅ aggregation might occur in an O₂-rich environment. As mentioned in the Introduction section, previous work on TiO₂-supported atomic copper clusters has shown that ACMs’ oxidation is preferred over their sintering, with the AMC–surface interaction playing a key role.²⁸ These outcomes suggest strategies to individualize Cu₅ clusters. This way, a clear spatial separation between the Cu₅ monomers happen when the energy gained by the adsorption of O₂ molecules is larger than the Cu₅–Cu₅ interaction energy itself (e.g., for n = 7 in Fig. 2). The decoupling of the two Cu₅ monomers through fluxional rotational motion is favored by an increase in the temperature as well (see Fig. 3).

3.1.2 Comparison with recent experimental measurements. Besides the quenching of aggregation due to adsorption of multiple O₂ molecules, the prediction of very similar phase diagrams for the three considered systems (see Fig. 2) is worth mentioning. Dimerization causes an overall decrease in the binding energies of O₂ molecules to the Cu₅ clusters. As a result, the inter-phase boundaries occur at higher pressures and lower temperatures for the Cu₅–Cu₅ dimer. Interestingly, this overall shifting improves the agreement with recent experimental measurements on fractions of copper oxidation states at low concentrations of copper clusters.¹³ Under ambient conditions (1 atm and RT), the experiment indicates a dominance of the Cu(II) oxidation state. This is consistent with the appearance of complexes with 7 O₂ molecules per cluster in the phase diagrams at 1 atm and RT (phases shown in green in the upper panels of Fig. 2). Upon lowering the pressure to 0.15 mbar, from RT up to 673 K, the dominance of the Cu(I) oxidation state has been experimentally resolved instead.¹³ This observation agrees with the occurrence of complexes with 3 O₂ per Cu₅ cluster as the most stable at 0.15 mbar from 350 K for the dimer in the singlet spin state (phase shown in yellow in the middle upper panel of Fig. 2). Under high-vacuum conditions (down to 10⁻⁷ mbar), from RT to 673 K, the dominance of the Cu(0) oxidation state has been experimentally identified.¹³ This outcome is also even with the predicted Cu₅ phase from 400 K at 10⁻⁷ mbar for the dimer (phase shown in the dark-red color in the middle upper panel of Fig. 2).

3.2 Oxidation of circumpyrene-supported Cu₅ clusters

In order to consider the effect of a carbon-based support, the PAH circumpyrene has been used (see Fig. 4). The dispersion D4 correction accounts for 82% of the attractive Cu₅-circumpyrene interaction at the optimized geometry. For the sake of comparison, a dispersion-dominated nature of the ACM–PAH interaction has been also found for the Au₁₂ cluster but is less pronounced (45% vs. 82%).

Fig. 4 Panel (a): phase diagram of carbon-supported Cu₅–(O₂)_n complexes within the molecular chemisorption model. The color coding is used to indicate the different (p-T) regions where the given Cu₅–(O₂)_n complex is the most stable, with the corresponding n value having been specified. Panels (b–d): structures corresponding to supported Cu₅–(O₂)₇ [panel (b)], Cu₅ Cu₅–(O₂)₇ [panel (c)], and Cu₅–(O₂)₃ [panel (d)].

Let us now analyze the contribution of a carbon-based support to the interaction of molecular oxygen with Cu₅. Considering the molecular chemisorption model, Fig. 5 shows the binding energies as a function of the number of adsorbed O₂ molecules onto circumpyrene-supported Cu₅. The values of the binding energies can be compared with those obtained for bare Cu₅ in Table 2. The left-hand panel of Fig. 5 presents the corresponding free energies of formation at ambient pressure (1 atm), while the most likely adsorption states as a function of oxygen pressure and temperature are depicted in the phase diagram of Fig. 4 (upper left-hand panel). The optimized structures of circumpyrene-supported Cu₅, as well as the most stable (Cu₅/circumpyrene)–(O₂)_n complexes are shown in Fig. 4.

Fig. 5 Free energies of formation of carbon-supported Cu₅–(O₂)_n complexes. Right-hand panel: electronic binding energies of the clusters as a function of the number of O₂ molecules attached. Left-hand panel: free energies of formation (w potential) comparison at 1 atm for circumpyrene-supported and unsupported Cu₅–(O₂)_n complexes. The color coding in solid lines follows dot colors in the predicted phase diagram. The energy values corresponding to less stable complexes are plotted with dashed lines. The predicted phases can be retrieved from the convex hull of the lowest-energy curves.

Table 2 Binding energies (second column) as a function of the number of O₂ molecules (n, first column) adsorbed onto circumpyrene-supported Cu₅. For the sake of comparison, the binding energies of O₂ molecules onto unsupported Cu₅ are also indicated (values in parentheses). The number of adsorbed neutral [n_O2], superoxo , and peroxo species is also specified. The ‘*’ and ‘**’ quotations indicate the dissociation of a single and all the adsorbed peroxo species, respectively. These species dissociate in three-fold coordinated O²⁻ ions and carbon-supported ad-atoms (or ad-ions)

	Cu5–(O2)n
		E bind/eV
n = 1	0/1/0 (0/1/0)	−1.4 (−1.4)
	0/0/1*	−1.7
n = 2	0/2/0 (0/2/0)	−3.4 (−2.5)
	0/0/2*	−3.9
	0/0/2**	−4.8
n = 3	0/3/0 (0/3/0)	−5.6 (−4.9)
	0/1/2*	−4.7
n = 4	0/4/0 (0/4/0)	−6.9 (−5.8)
	0/1/3*	−7.1
	0/1/3**	−7.2
n = 5	0/5/0 (2/3/0)	−7.8 (−5.8)
	0/3/2*	−7.6
	0/1/4**	−8.4
n = 6	0/6/0 (1/5/0)	−8.0 (−7.4)
	0/5/1*	−8.6
n = 7	0/7/0 (0/5/2)	−9.1 (−8.9)

---

Except for single O₂ adsorption, it can be observed from Table 2 that the support stabilizes the Cu₅–(O₂)_n complexes so that the binding energies become 2–26% larger. As a result, additional phases do appear with n = 4 and 5, and the inter-phase boundaries for n = 7 and 3 become shifted to lower pressures and higher temperatures as compared with unsupported Cu₅. As expected from the dispersion-dominated nature of the Cu₅-circumpyrene interaction, the Cu₅ trigonal bipyramidal structure is slightly modified by the support, with one side face becoming almost parallel to the circumpyrene plane. The same holds true for the adsorption of three O₂ molecules (see Fig. 2). When the number of adsorbed O₂ increases further up to 7, however, the Cu₅ pyramidal structure becomes much elongated along one of their vertices due to the stretching of Cu–Cu distances to host additional O₂ molecules. Concurrently, the binding of the complex to the support becomes weaker.

To gain insights into the fluxional nature of the complexes, we have carried out an AIMD simulation of the Cu₅–(O₂)₄ motion at 473 K (see Fig. 6). As indicated in Table 2, the four adsorbed O₂ molecules become charged as superoxo O₂⁻ ions. The conversion of the four O₂ molecules in superoxo radicals causes the enlargement of the O–O bond lengths (ca. 1.3 Å) upon that corresponding to the neutral O₂ molecule (ca. 1.2 Å). During the 12 ps of the AIMD simulation, no desorption of neutral O₂ molecules is observed. Accordingly, the O–O distance (shown in blue) changes very little. In contrast, the fluctuations of the Cu–O distances (to within 0.2 Å) are noticeable. They are associated with Cu–O bond breaking/formation which is concurrent to concerted contractions/elongations of Cu–Cu distances. The wide amplitude nature of the AMC motifs is also reflected in the large fluctuations (within 0.5 Å) of the Cu–Cu distances (shown in yellow and red in panel a of Fig. 6).

Fig. 6 Upper panel: figure showing the evolution of average Cu–Cu, Cu–O, and O–O distances during the first 8000 fs of an AIMD simulation of circumpyrene-supported Cu₅–(O₂)₄ at 473 K. Middle and bottom panels: snapshots from AIMD simulations: (1) the snapshot taken at 73.5 fs (middle panel) shows the Cu–O bond breaking and contraction of the Cu–Cu distance from the closest Cu atoms; (2) the snapshot taken at 101 fs shows the Cu–O bond formation and elongation of the Cu–Cu distance from the nearest Cu atoms.

3.2.1 Oxidation including support-mediated dissociative chemisorption. Besides the molecular chemisorption model, we have considered a support-mediated dissociative mechanism of peroxo species as well. As illustrated in Fig. 7, it consists in the dissociation of chemisorbed peroxo species into one O²⁻ ion and one O ad-atom (or ad-ion) along with the desorption of one superoxo O₂⁻ radical as a neutral O₂ molecule. The donated charge from the desorbing superoxo radical is transferred to one co-adsorbed superoxo molecule which thus becomes a peroxo species. Through dissociation, one O²⁻ ion is formed, becoming three-fold coordinated to the closest Cu cations of the copper cluster. The second product of the dissociative channel is one O atom (or O⁻ ion), which becomes adsorbed on top of two carbon atoms of the support (see Fig. 7).

Fig. 7 Picture illustrating the support-mediated dissociation of one peroxo O₂²⁻ species along with the desorption of one neutral O₂ molecule.

As can be observed in Fig. 8, the states characterized by the adsorption of O₂ molecules as superoxo radicals are very close to those arising from the occurrence of support-mediated dissociation. The lowest energy state (see Table 2) corresponds to a complex that can be characterized as Cu₅O₅. As typically found for fluxional atomic metal clusters, the energy differences between global and local minima are very small (to within 0.3 eV). Therefore, it seems sensible that the inter-conversion between them be reversible at experimentally relevant ranges of temperatures and oxygen gas pressures as discussed for single O₂ chemisorption on bare Cu₅ in ref. 11. In order to include the likelihood of inter-conversion between these energetically close-lying adsorption states (see Fig. 8), we have computed the Boltzmann-weighted average of free energies of formation for each circumpyrene-supported Cu₅–(O₂)_n complex as well as their associated distributions of copper oxidation states [Cu(0), Cu(I), Cu(II)].

Fig. 8 Picture showing the lowest-energy states of circumpyrene-supported Cu₅ upon adsorption of n = 4 O₂ molecules. The values of the binding energies are also indicated. Panel (a): the Cu₅ cluster with four O₂⁻ superoxo radicals attached to it. Panel (b): the Cu₅ cluster with one peroxo species having been dissociated into one three-fold coordinated O²⁻ anion and one O adatom onto the support. Panel (c): the Cu₅O₅ complex bearing three dissociated peroxo species and one superoxo O₂⁻ radical.

3.2.2 Comparison with recent experimental measurements. To compare with the experimental measurements,¹³ we have calculated a phase diagram showing the dominant oxidation states of circumpyrene-supported Cu₅ as a function of temperature and oxygen gas pressure. As can be observed in Fig. 9, under ambient conditions, the Cu(II) oxidation state is dominant, in agreement with the experiment. Also consistently with the experiment, upon lowering the oxygen pressure down to ca. 0.1 mbar at RT, the Cu(I) oxidation state becomes dominant. This feature is preserved within the experimentally relevant temperature range (up to 773 K in ref. 13). At 500 K, upon decreasing the oxygen pressure down to ca. 10⁻⁸ mbar, the Cu₅ cluster recovers its Cu(0) metallic state. In contrast, the Cu(0) state becomes predominant already at RT and under high-vacuum conditions (between 10⁻⁷ and 10⁻³ mbar) in the experimental measurements.¹³ This deviation can be explained as follows:

Fig. 9 Diagram showing the dominant oxidation states of copper atoms of circumpyrene-supported Cu₅ as a function of temperature and oxygen gas pressure.

(1) The aggregation of two Cu₅ clusters into the Cu₅–Cu₅ dimer induces a shift of the inter-phase boundaries to higher pressures and lower temperatures. This outcome might suggest that a portion of Cu₅ clusters might experience an aggregation process upon surface deposition.

(2) Strong non-adiabatic effects featured by AMCs^7,12 induce a ‘trapping’ of O₂ molecules approaching Cu₅ in an electronic state correlating with the neutral O₂ and Cu₅ fragments at the asymptotic limit. These effects prevent the system to access the charge-transfer chemisorption state, thus favoring the dominance of the Cu(0) metallic state over that predicted through the consideration of chemisorption states only.

3.2.3 The circumpyrene-supported Cu₅–Cu₅ dimer. Aimed to analyze temperature effects on circumpyrene-supported Cu₅–Cu₅, we have carried out AIMD simulations at a temperature of 573 K. As can be seen in the left-hand panel of Fig. 10, the lowest-energy structure at 0 K suffers a modification upon heating through the rotation of the bipyramidal structures forming it. As for unsupported Cu₅–Cu₅, the identity of the Cu₅ sub-units composing the dimer shows up due to the ‘floppy’ nature of the Cu₅–Cu₅ rotational motion. Temperature-dependent experimental studies, using deposition of the Cu₅ clusters in an ultrahigh vacuum, could provide very useful insights into their coupling to form the Cu₅–Cu₅ dimer and larger aggregates. Our results suggest that their subsequent exposition of the sample in an O₂ atmosphere could be an useful tool to individualize them for further spectroscopic characterization as well.

Fig. 10 Snapshots from AIMD simulations showing the decoupling of the Cu₅–Cu₅ dimer through the rotation of Cu₅ bypiramidal 3D structures upon increasing the temperature (up to 573 K in the given snapshot after 1.5 ps of simulation). For the sake of clarity, the Cu–Cu bonds joining the two Cu₅ subunits are not shown.

4 Concluding remarks and future perspectives

We have identified aggregation and the role of a chemically inert support in the interaction of gas-phase O₂ molecules with Cu₅ as a paradigmatic case of a fluxional AMC. Our results can be summarized as follows:

(1) The aggregation of two Cu₅ clusters into a Cu₅–Cu₅ dimer formed by coupled bipyramidal 3D structures is predicted to be energetically much favored, with the attractive Cu₅–Cu₅ interaction being about −5.2 eV. However, the decoupling of the Cu₅–Cu₅ dimer into two Cu₅ clusters has been apparent in the corresponding (p, T)-phase diagram upon increasing the oxygen gas pressure. Molecules of O₂ adsorb in between Cu₅ clusters, causing an enlargement of their O–O bonds, the formation of peroxo species, and even their dissociation. AIMD simulations have illustrated the rotational decoupling of the Cu₅–Cu₅ dimer into two Cu₅ sub-units, with and without including support effects, upon heating.

(2) Without including support effects, Cu₅–(O₂)_n and (Cu₅–Cu₅)–(O₂)_2n complexes with n = 3 and 7 are the most stable. The phase diagrams of Cu₅ and Cu₅–Cu₅ are very similar, differing just in an overall shifting of the inter-phase boundaries to higher pressures and lower temperatures for the Cu₅–Cu₅ dimer.

(3) Considering circumpyrene as a molecular model of a carbon-based surface, the role of an inert support has been manifested in the overall stabilization of Cu₅–(O₂)_n complexes in such a way that additional (p,T)-phases do appear. Moreover, lower pressures and higher temperatures are necessary for the desorption of neutral O₂ molecules when the copper clusters are supported.

(4) An AIMD simulation has illustrated the fluxionality and wide amplitude nature of the Cu–Cu and Cu–O motions in circumpyrene-supported Cu₅–(O₂)₄. We have identified the contraction of Cu–Cu distances upon Cu–O bond breaking and their enlargement upon Cu–O bond formation.

(5) A support-mediated mechanism for dissociation of peroxo species adsorbed onto Cu₅ has been revealed, with the energies of the dissociative states being very close to those characterized by molecular chemisorption. Considering both molecular and dissociative reaction channels, a diagram has been provided to show the dominant copper oxidation states as a function of temperature and oxygen gas pressure, being consistent with recent experimental measurements.¹³

Through first-principles modelling, we have shown that gas-phase O₂ molecules can be adsorbed as either superoxo or peroxo species onto atomic metal clusters (AMCs) and that both molecular and dissociative chemisorption can occur. The energetically favored aggregation of AMCs and occurrence of multiple energetically closely lying chemisorption states upon their interaction with O₂ molecules have the potential of impacting both their stability and catalytic properties, either under aggregation and/or oxidation conditions. On the one hand, the structural fluxionality of ACMs allows their rotational decoupling from formed aggregates, such as the Cu₅ dimer both upon heating, and their spatial separation at exposition into an O₂-rich environment. On the other hand, the same properties might be the one responsible for the appearance of many reaction channels, such as support-mediated dissociative chemisorption, with the resulting adsorption states having very similar energies to those characterized by, e.g., molecular chemisorption. This characteristic might favor a reversible inter-conversion between different adsorption states at experimentally relevant ranges of temperatures and oxygen gas pressures. Considering circumpyrene-supported Cu₅ as a model system, we have shown that a Boltzmann average of a representative ensemble of close-lying states and their corresponding copper oxidation states allow reaching a sensible agreement with recent in situ experimental measurements on HOGP-supported Cu₅ clusters.¹³ In order to gain further insights into the oxidation of Cu₅, it is, however, necessary to consider additional channels for chemisorption of O₂ both in molecular and dissociative forms, particularly those involving peroxo species. The achievement of a complete picture of the aggregation process for two Cu₅ clusters, with and without O₂ molecules, requires extensive AIMD simulations as a function of temperature. Future work is planned to address them, using periodic models of graphene, properly addressing the diffusion of Cu₅ clusters as well. An investigation of the influence of the ACMs’ aggregation in the electronic structures and geometries of individual ACMs is also a future perspective. Despite the need for further research, we consider that our study provides a solid support to previous works addressing the reversible nature of the oxidation of the smallest copper clusters,^7,9–13 including in situ experimental measurements.¹³

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank Andreas W. Hauser for having provided Python and MATLAB codes that have served to write a general-purpose code, delivering phase diagrams. We are also thankful to Alexandre Zanchet, Lyudmila Moskaleva, Berta Fernández, and Alexander O. Mitrushchenkov for very fruitful discussions, and Mehrdad R. Osanloo for the proof-reading of the revised manuscript. This work was supported by the Spanish Agencia Estatal de Investigación (AEI) under Grant No. PID2020-117605GB-I00. This publication is also based upon work of COST Action CA21101 “Confined molecular systems: from a new generation of materials to the stars” (COSY) supported by COST (European Cooperation in Science and Technology). J. A.-G. acknowledges the ‘JAE-Intro 2021’ program for having provided a scholarship at CSIC. The CESGA super-computer center (Spain) and the CTI-CSIC are acknowledged for having provided computational resources. We also acknowledge the support of the publication fee by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI). The EU project of reference HORIZON-MSCA-2021-DN-01 is acknowledged as well.

Notes and references

1. S. Huseyinova, J. Blanco, F. G. Requejo, J. M. Ramallo-López, M. C. Blanco, D. Buceta and M. A. López-Quintela, J. Phys. Chem. C, 2016, 120, 15902–15908.
2. A. Halder, L. A. Curtiss, A. Fortunelli and S. Vajda, J. Chem. Phys., 2018, 148, 110901.
3. S. Hirabayashi and M. Ichihashi, Phys. Chem. Chem. Phys., 2014, 16, 26500–26505.
4. P. López-Caballero, A. W. Hauser and M. P. de Lara-Castells, J. Phys. Chem. C, 2019, 123, 23064–23074.
5. A. Halder, C. Lenardi, J. Timoshenko, A. Mravak, B. Yang, L. K. Kolipaka, C. Piazzoni, S. Seifert, V. Bonačić-Koutecký, A. I. Frenkel, P. Milani and S. Vajda, ACS Catal., 2021, 11, 6210–6224.
6. M. P. de Lara-Castells, A. W. Hauser, J. M. Ramallo-López, D. Buceta, L. J. Giovanetti, M. A. López-Quintela and F. G. Requejo, J. Mater. Chem. A, 2019, 7, 7489–7500.
7. M. P. de Lara-Castells, J. Colloid Interface Sci., 2022, 612, 737–759.
8. J. Juraj, A. Fortunelli and S. Vajda, Phys. Chem. Chem. Phys., 2022, 120, 12083–12115.
9. E. Fernández, M. Boronat and A. Corma, J. Phys. Chem. C, 2015, 119, 19832–19846.
10. P. Concepción, et al. , ACS Catal., 2017, 7, 3560–3568.
11. A. Zanchet, P. López-Caballero, A. O. Mitrushchenkov, D. Buceta, M. A. López-Quintela, A. W. Hauser and M. P. de Lara-Castells, J. Phys. Chem. C, 2019, 123, 27064–27072.
12. A. O. Mitrushchenkov, A. Zanchet, A. W. Hauser and M. P. de Lara-Castells, J. Phys. Chem. A, 2021, 125, 9143–9150.
13. D. Buceta, S. Huseyinova, M. Cuerva, H. Lozano, L. J. Giovanetti, J. M. Ramallo-López, P. Lóopez-Caballero, A. Zanchet, A. O. Mitrushchenkov, A. W. Hauser, G. Barone, C. Huck-Iriart, C. Escudero, J. C. Hernández-Garrido, J. Calvino, M. Lopez-Haro, M. P. de Lara-Castells, F. G. Requejo and M. A. López-Quintela, ChemRxiv, 2021, 1–36,  DOI:10.26434/chemrxiv.13661081.v1.
14. P. Concepción, M. Boronat, S. García-García, E. Fernández and A. Corma, ACS Catal., 2017, 7, 3560–3568.
15. T. L. Thomson and J. T. Yates, Jr., Chem. Rev., 2006, 106, 4428–4453.
16. J. T. Yates Jr., Surf. Sci., 2009, 603, 1605–1612.
17. Z. Zhang and J. T. Yates Jr., Chem. Rev., 2012, 112, 5520–5551.
18. M. A. Henderson, Surf. Sci. Rep., 2011, 66, 185–197.
19. M. A. Henderson, M. M. Shen, Z. T. Wang and I. Lyubinetsky, J. Phys. Chem. C, 2013, 117, 5774–5784.
20. C. Shu, N. Sukumar and C. P. Ursenbach, J. Chem. Phys., 1999, 110, 10539–10544.
21. M. P. de Lara-Castells and J. L. Krause, J. Chem. Phys., 2001, 115, 4798–4810.
22. M. P. de Lara-Castells and J. L. Krause, Chem. Phys. Lett., 2002, 354, 483–490.
23. M. P. de Lara-Castells and J. L. Krause, J. Chem. Phys., 2003, 118, 5098–5105.
24. M. P. de Lara-Castells, A. O. Mitrushchenkov, O. Roncero and J. L. Krause, Isr. J. Chem., 2005, 45, 59–76.
25. H. Zhai and A. N. Alexandrova, ACS Catal., 2017, 7, 1905–1911.
26. Q.-Y. Fan, Y. Wang and J. Cheng, J. Phys. Chem. Lett., 2021, 12, 3891–3897.
27. O. V. Lushchikova, H. Tahmasbi, S. Reijmer, R. Platte, J. Meyer and J. M. Bakker, J. Phys. Chem. A, 2021, 125, 2836–2848.
28. S. K. Iyemperumal, T. G. Fenton, S. L. Gillingham, A. D. Carl, R. L. Grimm, G. Li and N. A. Deskins, J. Chem. Phys., 2019, 151, 054702.
29. Q. Chen, D. Schollmeyer, K. Müllen and A. Narita, J. Am. Chem. Soc., 2019, 141, 19994–19999.
30. M. P. de Lara-Castells, A. O. Mitrushchenkov and H. Stoll, J. Chem. Phys., 2015, 143, 102804.
31. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
32. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
33. P. E. Blöch, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
34. A. Muñoz-Castro, T. Gomez, D. M. Carey, S. Miranda-Rojas, F. Mendizabal, J. H. Zagal and R. Arratia-Perez, J. Phys. Chem. C, 2016, 120, 7358–7364.
35. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
36. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
37. M. P. de Lara-Castells, C. Cabrillo, D. A. Micha, A. O. Mitrushchenkov and T. Vazhappilly, Phys. Chem. Chem. Phys., 2018, 20, 19110–19119.
38. P. López-Caballero, J. M. Ramallo-López, L. J. Giovanetti, D. Buceta, S. Miret-Artés, M. A. López-Quintela, F. G. Requejo and M. P. de Lara-Castells, J. Mater. Chem. A, 2020, 8, 6842–6853.
39. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
40. E. Caldeweyher, S. Ehlert, A. Hansen, H. Neugebauer, S. Spicher, C. Bannwarth and S. Grimme, J. Chem. Phys., 2019, 150, 154122.
41. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
42. S. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553–566.
43. R. S. Mulliken, J. Chem. Phys., 1962, 36, 3428–3439.
44. P. Bultinck, C. Van Alsenoy, P. W. Ayers and R. Carbó-Dorca, J. Chem. Phys., 2007, 126, 144111.
45. F. L. Hirshfeld, Theor. Chim. Acta, 1977, 44, 129–138.
46. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 73–78.
47. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1327.
48. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2022, 12(5), e1606.
49. H.-J. Werner, P. J. Knowles, F. R. Manby, J. A. Black, K. Doll, A. Heßelmann, D. Kats, A. Köhn, T. Korona, D. A. Kreplin, Q. Ma, T. F. Miller, A. Mitrushchenkov, K. A. Peterson, I. Polyak, G. Rauhut and M. Sibaev, J. Chem. Phys., 2020, 152, 144107.
50. P. Celani and H.-J. Werner, J. Chem. Phys., 2000, 112, 5546–5557.
51. D. E. Woon and T. H. Dunning, Jr., J. Chem. Phys., 1994, 100, 2975–2988.
52. D. Figgen, G. Rauhut, M. Dolg and H. Stoll, Chem. Phys., 2005, 311, 227–244.
53. A. Kubas, D. Berger, H. Oberhofer, D. Maganas, K. Reuter and F. Neese, J. Phys. Chem. Lett., 2016, 7, 4207–4212.
54. M. Mukherjee, D. Tripathi, M. Brehm, C. Riplinger and A. K. Dutta, J. Chem. Theory Comput., 2021, 17(1), 105–116.
55. M. Weiβ and M. Brehm, Molecules, 2020, 25(24), 5861.
56. C. Slawik, C. Rickmeyer, M. Brehm, A. Böhme and G. Schüürmann, Environ. Sci. Technol., 2017, 51(7), 4018–4026.
57. A. W. Hauser, J. Gomes, M. Bajdich, M. Head-Gordon and A. T. Bell, Phys. Chem. Chem. Phys., 2013, 15, 20727–20734.
58. X. Yu, A. R. Oganov, Q. Zhu, F. Qi and G. Qian, Phys. Chem. Chem. Phys., 2018, 20, 30437–30444.
59. S. Bhattacharya, S. V. Levchenko, L. M. Ghiringhelli and M. Scheffler, Phys. Rev. Lett., 2013, 111, 135501.
60. Y. Xu, W. A. Shelton and W. F. Schneider, J. Phys. Chem. B, 2006, 110, 16591–16599.
61. NIST Chemistry WebBook, NIST Standard Reference Database Number 69, ed. P. J. Linstrom and W. G. Mallard, National Institute of Standards and Technology, Gaithersburg MD, 2019, p. 20899 DOI:10.18434/T4D303.
62. M. W. Chase, Jr., NIST-JANAF Themochemical Tables, Fourth Edition, J. Phys. Chem. Ref. Data, Monograph 9, 1998, 1–1951.
63. P. L. Rodríguez-Kessler, A. R. Rodríguez-Domínguez, D. MacLeod-Carey and A. Muñoz-Castro, Adv. Theory Simul., 2021, 4, 2100043.

---

Footnote

† Electronic supplementary information (ESI) available: Complementary code written in Octave/MATLAB to extract (p–T)-phase diagrams from the thermochemistry output of the ORCA code. See DOI: https://doi.org/10.1039/d2cp02169b

---