Theory as a driving force to understand reactions on nanoparticles: general discussion

Carlo Lamberti opened discussion of the introductory lecture by Bruce Gates: As you have shown, EXAFS plays a crucial role in determining the structure and the nuclearity of nanoparticles (NPs). For each shell, the accuracy of this determination depends on the error bar associated to the coordination number, that strongly correlates with the corresponding Debye–Waller (DW) parameter. This becomes even more important when in situ operando experiments are performed at reaction temperature. Based on your experience, what suggestions can you give to reduce this correlation and increase the potentiality of the technique? Do you believe it is possible to fix or to determine, in a reliable way, DW parameters from independent experimental or computational works? Do you believe that in temperature-dependent experiments it is reliable to adopt the Debye or the Einstein model^1,2 to parametrize the evolution of DW parameters?

1 G. Dalba, P. Fornasini, R. Grisenti and J. Purans, Phys. Rev. Lett., 1999, 82, 4240–4243.

2 S. Øien, G. Agostini, S. Svelle, E. Borfecchia, K. A. Lomachenko, L. Mino, E. Gallo, S. Bordiga, U. Olsbye, K. P. Lillerud and C. Lamberti, Chem. Mater., 2015, 27, 1042–1056.

Bruce Gates answered: You raise good points, and all I would like to state is that it is valuable to have corroborating evidence from other techniques to determine the coordination number. For example, triosmium clusters on a support, if synthesized with high precision, can be characterized by STEM to determine the cluster nuclearity (hence the Os–Os coordination number) for comparison with the EXAFS value. Furthermore, IR spectra of triosmium carbonyls provide evidence of the structure (including the cluster nuclearity). Insofar as such comparisons have been made, for various osmium clusters on supports, the data confirm the Os–Os coordination number determined by EXAFS spectroscopy. See, for example, ref. 1–3. Perhaps samples such as these can be used in experimentation to address the questions you have raised.

1 N. L. Okamoto, B. W. Reed, S. Mehraeen, A. Kulkarni, D. G. Morgan, B. C. Gates and N. D. Browning, Determination of Nanocluster Sizes from Dark-Field Scanning Transmission Electron Microscopy Images, J. Phys. Chem. C, 2008, 112, 1759.

2 A. Kulkarni, S. Mehraeen, B. W. Reed, N. L. Okamoto, N. D. Browning and B. C. Gates, Nearly Uniform Decaosmium Clusters Supported on MgO: Characterization by X-ray Absorption Spectroscopy and Scanning Transmission Electron Microscopy, J. Phys. Chem. C, 2009, 113, 13377.

3 S. Mehraeen, A. Kulkarni, M. Chi, B. W. Reed, N. L. Okamoto, N. D. Browning and B. C. Gates, Triosmium Clusters on a Support: Determination of Structure by X-Ray Absorption Spectroscopy and High-Resolution Microscopy, Chem.– Eur. J., 2011, 17, 1000.

Rutger van Santen said: You mentioned the unique stability of cubic Ir nanoparticles of particular size. Do the total number of atoms in these particles agree with the ideal structures? Can one exclude the possibility that the presence of particular surface edge features or surface reconstruction, that is only stable when a particular surface size is reached, are the explanations?

Bruce Gates responded: The TEM images show a range of sizes of Ir species, ranging from the single-atom complexes to the nanoparticles that are all about 1 nm in diameter or less. The images indicate various nanoparticle morphologies, and one that was emphasized in the presentation is evidently cubic. Whether the distribution of nanoparticles evolves to cubic nanoparticles after long times is not determined by the data. Neither the experimental results nor the theory of Pawlik et al. (ref. 51 in the paper) exclude the possibilities you have suggested.

Justin Hargreaves asked: Taking the analogy with organometallic chemistry further, to what extent is it possible to produce an empirical ranking of supports in terms of some parameter akin to the Tolman electronic parameter?

Bruce Gates replied: The expectations one would have on the basis of organometallic catalysis in solution, in my view, extend seamlessly to supported metal complex (and, presumably, metal cluster) catalysts when they have a high degree of uniformity. For example, ref. 1 reported correlations of the activities of supported mononuclear iridium complexes (measured as turnover frequencies) for ethylene hydrogenation and for ethylene dimerization with the carbonyl stretching frequencies of the iridium complexes in the catalysts after they were exposed to CO to form anchored iridium gem-dicarbonyls. These frequencies are a measure of the electron-donor tendency of the supports, which are ligands. Thus, the correlations provide the kind of empirical ranking that you are referring to, and they represent a family of supports, some being electron donating and some being electron withdrawing, and account for orders of magnitude ranges in the catalytic activities. It is important in this context that the supported catalysts are not highly non-uniform on the supports and thus nearly unique (and essentially molecular).

1 D. Yang, S. O. Odoh, T. C. Wang, O. K. Farha, J. T. Hupp, C. J. Cramer, L. Gagliardi and B. C. Gates, J. Am. Chem. Soc., 2015, 137, 7391–7396.

Kassim Badmus commented: Why does a zeolite not give consistent result in its characterization and when it is used as a support? Can you tell us the factors that must be considered before we choose a support for a nanoparticle? Can the pore size be responsible for the inconsistency in characterization of zeolite systems?

Bruce Gates replied: For a given zeolite sample, our data show good reproducibility. But zeolite syntheses give samples with variable compositions (distributions of Si and Al sites) and the initially formed crystals of a zeolite generally don’t match the ones formed later in a batch synthesis; furthermore, synthesis of many zeolites is challenging to reproduce. Some syntheses give more than one zeolite, and some samples of zeolites incorporate amorphous material. In general, in finding porous supports (zeolites or others) for metal nanoparticle catalysts, one must consider the support surface chemistry, because it influences the synthesis of the supported species, and the pore size distribution, because mass transfer of reactants and products in the pores can affect rates of catalytic reactions (and blocking of small pores by the nanoparticles can occur).

Maurits Boeije asked: Based on your previous work,¹ can you draw the general conclusion that partially encapsulated nanoparticles in a matrix can improve stability of a catalyst by restricting nanoparticle motion and preventing coalescence?

1 J. Zhang, L. Wang, Y. Shao, Y. Wang, B. C. Gates and F.-S. Xiao, Angew. Chem., Int. Ed., 2017, 56, 9628.

Bruce Gates answered: With proper preparation, zeolite-encapsulated nanoparticles of various metals can indeed be stabilized against sintering and coke formation; details are to be published soon; see ref. 48 of the paper.

Carlo Lamberti asked: Based on your experience, do you believe that EXAFS is able to discriminate among metal–carbon, metal–oxygen and metal–nitrogen bonds? Do you believe that the recent experimental and theoretical progress of X-ray Emission spectroscopy (XES)^1,2 will promote the technique as a standard characterization tool in the near future? Recently XES has been able to discriminate between first-shell Cu–O and Cu–N bonds in Cu–CHA catalyst during NH₃-assisted selective catalytic reduction of NO_x.^3,4 In this regard, XES even succeeded in the discrimination between Al and P in the second shell environment of Ti-AlPO-5 catalyst.⁵ An exhaustive understanding of the XES spectra however requires the theoretical support of DFT calculations.^6,7

1 P. Glatzel and U. Bergmann, Coord. Chem. Rev., 2005, 249, 65–95.

2 J. Singh, C. Lamberti and J. A. van Bokhoven, Chem. Soc. Rev., 2010, 39, 4754–4766.

3 F. Giordanino, E. Borfecchia, K. A. Lomachenko, A. Lazzarini, G. Agostini, E. Gallo, A. V. Soldatov, P. Beato, S. Bordiga and C. Lamberti, J. Phys. Chem. Lett., 2014, 5, 1552–1559.

4 K. A. Lomachenko, E. Borfecchia, C. Negri, G. Berlier, C. Lamberti, P. Beato, H. Falsig and S. Bordiga, J. Am. Chem. Soc., 2016, 138, 12025–12028.

5 E. Gallo, A. Piovano, C. Marini, O. Mathon, S. Pascarelli, P. Glatzel, C. Lamberti and G. Berlier, J. Phys. Chem. C, 2014, 118, 11745–11751.

6 E. Gallo, C. Lamberti and P. Glatzel, Phys. Chem. Chem. Phys., 2011, 13, 19409–19419.

7 E. Borfecchia, K. A. Lomachenko, F. Giordanino, H. Falsig, P. Beato, A. V. Soldatov, S. Bordiga and C. Lamberti, Chem. Sci., 2015, 6, 548–563.

Bruce Gates responded: Thanks for this timely question. I am basically in agreement with your views stated in your question. I believe that structures determined by EXAFS spectroscopy have substantially more value when they are bolstered by results determined by complementary characterization methods, such as vibrational spectroscopies, TEM imaging, and theory. Structures inferred from EXAFS spectra, in my view, are far better justified when they are based on data characterizing structurally nearly uniform samples – and most solid catalysts incorporate surface structures that are non-uniform, with the catalytically relevant species often being minorities, sometimes too sparse to even characterize by EXAFS spectroscopy. I believe your point about XES is pertinent and that the value of this technique in catalysis will become even more evident when XES is applied to structurally well-defined samples such as the Ti-AlPO-5 you mention and catalysts like the ones I mentioned in my talk.

Graham Hutchings remarked: You have shown some very elegant microscopy, in which there are Ir dimers in one example and Os trimers in another example; what happens if you put a second metal e.g. Ir into the osmium system or vice versa, would you still observe separate dimers and trimers?

Bruce Gates answered: With today’s aberration-corrected STEM capabilities, the experiments you have suggested are quite challenging, because a good structure determination requires that the two metals in a bimetallic be readily distinguished from each other, which requires that they have significantly different atomic numbers. Os and Ir are too close in atomic number. In principle, one could distinguish, say, Rh and Ir, although Rh is so light that there are still only a few examples in the literature showing Rh atoms with atomic resolution on a support. A further limitation of investigating bimetallics on supports by aberration-corrected STEM is that the metals on the support need to be quite different in atomic number from the atoms of the support; thus, for example, Ir atoms or Os atoms on MgO yield excellent images; it helps that the MgO (powder), if properly prepared, is highly crystalline and allows identification of various MgO crystal faces.

Annette Trunschke commented: Thank you very much for your interesting lecture. I was particularly impressed by your results that clearly show the structural changes of metal clusters depending on the reaction conditions. In this regard, I am interested in your opinion on general approaches in catalyst characterization in the future. Is it worth or necessary to investigate the fresh catalyst with high precision or should we concentrate our efforts on operando experiments?

Bruce Gates replied: No doubt careful characterization of fresh catalysts is valuable, especially insofar as it helps us to understand what is going on in catalyst synthesis. But I agree that, at least in prospect, in operando investigations provide the most valuable catalyst characterizations. This is easy to say and not always easy to do; in part because it can be challenging (especially when the relevant conditions involve complex feedstocks, high pressures, and high temperatures) to apply the methods to ensure that they all provide characterization of the same catalyst. In principle, the more such methods can be applied the better, and it is advantageous when all the measurements are done with the same apparatus – optimally, in my view, this prospect would allow measurement of vibrational spectra, X-ray absorption spectra, images, and more, along with catalyst performance data. This is in my opinion an essential question and can help motivate advances in the characterization methods.

Philip Davies opened the discussion of the paper by Rutger van Santen: In your very interesting paper you highlight the difficulties of building accurate kinetic models of catalytic systems. What, in your view, are the experimental advances needed to provide the data necessary to make calculations more realistic?

Rutger van Santen responded: Reliability of the microkinetics simulated predictions is improved by validation of the supporting quantum-chemically calculated elementary reaction rates by experiments that focus on a comparison of such rates measured at a molecular level rather than comparison with macroscopic kinetics. Agreement with the latter is never a guarantee that the predicted mechanism of the reaction is actually correct, because kinetics will lump the molecular information together.

Richard Catlow commented: One general issue that needs to be considered when discussing the interplay between theory and experiment is often the state of the catalysts is not well defined making detailed comparisons difficult. This problem needs to be addressed by a joint computational–experimental approach.

Cynthia Friend responded: I wholeheartedly agree.

Bruce Gates said: I concur and would emphasize that comparisons of experiment with theoretical predictions of supported catalysts can be especially fruitful when the supported species are synthesized precisely to give samples that have a high degree of uniformity, and there are now some examples showing good agreement between theory and experiment and opportunities for further work in this direction. For some recent examples, see ref. 1.

1 B. C. Gates, M. Flytzani-Stephanopoulos, D. A. Dixon and A. Katz, Atomically dispersed supported metal catalysts: perspectives and suggestions for future research, Catal. Sci. Technol., 2017, 7, 4259

Katerina Soulantica stated: I strongly agree with this comment and I am persuaded that synthetic protocols which reproducibly afford a variety of real-catalysts with the specific characteristics predicted to be necessary from modeling and mechanistic studies on model systems, are a prerequisite. The standard existing procedures for real-catalyst preparation are not well adapted for such a high degree of control. The best catalyst configuration may correspond to a real synthetic challenge, but several procedures of well-defined nanoparticles are already available and this is a good starting point. I believe that in the future we could think about creating a library of benchmark synthetic procedures in parallel to benchmark modeling techniques in order to be able to make the optimal calculated catalyst as predicted from theory and in situ measurements a reality.

Rutger van Santen commented: Yes, I completely agree. Methods are available to study surface reconstruction, using molecular dynamics or equilibrium approaches that establish the state of a surface in equilibrium with a reactant medium. It may even be necessary to consider the transition between different surface states when reactions are oscillatory. For not too complex reactions (oxidation of CO by transition metals or oxides) kinetic Monte Carlo simulations are available that demonstrate this computationally. However the timescales of surface reconstruction (sometimes activated and slow) will usually not match the timescale of the catalytic reaction cycle and then the two have to be simulated independently. There is usually very limited information on the state of the working catalyst under practical conditions. The latter is essential because it often sensitively depends on conditions (the pressure, temperature gap). Such measurements will help to validate kinetics simulations.

Keith Whiston asked: Is it possible using your microkinetic modelling approach to incorporate aspects of zeolite geometry as predictors within the model? Either by using them to represent the Brønsted acidic properties of the catalysts and also to predict the effects of diffusion on product distribution and deactivation rate?

Rutger van Santen answered: One property that is very sensitive to zeolite pore geometry is the adsorption isotherm of hydrocarbons. There is a very approximate relation between matching molecular size and cavity shape. But there are methods, such as the Configurational-Bias Monte Carlo (CB-MC) method, that provide the relevant numbers directly.

The case of diffusion is more complicated. Knudsen diffusion does not apply to zeolites, since diffusion in zeolites is of the ballistic type. Also for this case excellent molecular dynamics approaches are available to estimate the corresponding diffusion constants. If one would like to define predictors, they relate again to a match between molecule shape and volume with that of the microcavity.

A Monte Carlo approach would enable the inclusion of diffusion in the kinetics modelling. This implies definition of size of crystallite and explicit consideration of zeolite nanopore topology. Since the diffusion of the small reactant molecules we consider is fast compared to reaction rate, we did not include diffusion explicitly in our microkinetics modelling.

The effect of varying proton acidity, assuming that zeolite structure remains the same, can be readily incorporated in the microkinetics simulations by varying the activation energies of the elementary reaction steps using the BEP linear activation energy–reaction energy relationships, that are valid as long as the structure of the reaction intermediates does not vary. A probe of the proton reactivity is its bond strength, which can be measured in several ways.

There is no general rule to describe activation free energies as a function of curvature. Generally one expects them to increase when curvature increases, because it inhibits approach of reactant to proton site.

Andrea Russell commented: When reading Prof. van Santen’s paper and then also listening to Prof. Gates’ opening address to this conference, a recent opinion piece published in Chemistry World¹ came to mind. In that article the author speculated that developments in machine learning would make synthetic chemists, especially organic chemists, redundant when it came to the discovery of new reactions and reaction schemes. Do you think that the same can be said for catalysis and, if so, what information needs to be provided in experimental reports to facilitate such an advance?

1 Derek Lowe, Will robots make you redundant?, Chem. World, 29 March 2018.

Bruce Gates responded: This is a provocative thought, but my sense is that robots will not in the foreseeable future take many jobs away from scientists working to find better catalysts, because of the complexity and subtlety of catalysis and the complexity of the structures of the surfaces of solid catalysts. Even if robots could predict optimal catalyst structures, they would be challenged (as we are) to synthesize them and find ways to stabilize them. Nonetheless, the idea seems to be an extension of the technology of rapid-throughput testing in catalyst discovery, and its value is, in my judgement, well demonstrated (if not well documented), but mostly for indications of material compositions offering tantalizing initial catalytic activities and selectivities. Professor van Santen has written thoughtful books about the future of technology, and his thoughts about this matter will have much more substance than what I have stated here.

Rutger van Santen remarked: Machine learning requires training of systems on correlations between existing data. This process is not model based. Empirically it would be useful to have ready access to such a database, but it can be hardly be considered to be predictive. It would be a poor expert system, that so far for catalysis has been of little use. Essentially because we still have no ultimate predictive understanding.

With machine learning, data are not used to construct a mechanistic model of the relationship between the performance of a reaction and the catalyst structure and composition. For many reactions such mechanistic understanding, including information on the structure of the catalyst during a reaction, is absent, so the machine does not have the information to be trained on. Whereas such mechanistic models are necessary to be predictive. The theoretical catalysis programme that has been of increasing relevance the past twenty years has as its very aim to provide performance–structure relationships based on mechanistic models of the reactions and catalyst site reorganisation. It seems that machine learning techniques when applied to a large data set that contain substantial errors in accuracy are useful to reduce the error margin of the actual numbers to be used. Then it can be a useful tool to reduce the accuracy of calculations using approximate methods only applicable to very large and complex systems.

Francesca Baletto asked: Could you provide more information on the importance of site reconstruction and comment whether the lattice mobility should be include into a microkinetic model?

Rutger van Santen answered: Especially in transition metal catalysis the phenomenon of site reconstruction when the catalyst is exposed to a reactive medium is quite general. In order to actually determine the structure and composition of the surface overlayer that forms, displacement of the lattice atoms has to be taken into account. However the timescale of a catalytic reaction cycle is such that on that timescale the reaction can usually be assumed to take place on a surface where the transition atom mobility is slow. Then the determination of the surface structure and composition can be done independently once an overlayer concentration of adatoms has been established.

This however is not generally the case. Exceptions are surface reactions that self-organize, as the Ertl-related systems and systems where a liquid overlayer forms, as is most likely the case in oxychloride systems.

Carlo Lamberti said: Concerning the machine learning (ML) approach, it is worth mentioning the recent work of Frenkel and co-workers,¹ which has shown that the size, shape and morphology of Pt nanoparticles (NPs) can be obtained from XANES spectroscopy supported by a ML approach. The ML method was trained with ab initio XANES simulations on a huge library of clusters. Consistent results were obtained simulating the spectra with both FEFF-9² and FDMNES³ codes, resulting in the correct 3D reconstruction of the NPs. On the other hand, a library of XANES spectra for Pd hydrate and Pd carbide phases has been created, on DFT-optimized geometries, changing the Pd–Pd distance and the x and y stoichiometries of the PdH_x or PdC_y phases (ref. 4 and DOI: 10.1039/c7fd00211d) This allowed the authors to determine both structure and composition of Pd NPs under hydrogenation reactions. These kinds of studies demonstrated that XANES spectroscopy can be applied for high-throughput, time-dependent, studies typical of operando investigation of a catalytic system.

1 J. Timoshenko, D. Lu, Y. Lin and A. I. Frenkel, J. Phys. Chem. Lett., 2017, 8, 5091–5098.

2 J. J. Rehr, J. J. Kas, F. D. Vila, M. P. Prange and K. Jorissen, Phys. Chem. Chem. Phys., 2010, 12, 5503–5513.

3 S. A. Guda, A. A. Guda, M. A. Soldatov, K. A. Lomachenko, A. L. Bugaev, C. Lamberti, W. Gawelda, C. Bressler, G. Smolentsev, A. V. Soldatov and Y. Joly, J. Chem. Theory Comput., 2015, 11, 4512–4521.

4 A. L. Bugaev, O. A. Usoltsev, A. A. Guda, K. A. Lomachenko, I. A. Pankin, Y. V. Rusalev, H. Emerich, E. Groppo, R. Pellegrini, A. V. Soldatov, J. A. van Bokhoven and C. Lamberti, J. Phys. Chem. C, 2018, 122, 12029–12037.

Bruce Gates replied: I would add the thought that further progress in this direction might be facilitated by work with structurally well-defined and nearly uniform supported species and not just samples such as those of Frenkel et al. that consist of a smear of structures.

Yaroslav Odarchenko said: Thank you very much for your talk. Our group is also studying the deactivation mechanism of the FTS catalyst. In your paper you discuss only the promotional effect of water. However, in our experimental work on Co-based catalyst¹ we have observed that metal nanoparticles oxidize most probably due to the presence of water. What are your thoughts about this?

1 P. Senecal et al., ACS Catal., 2017, 7, 2284–2293.

Rutger van Santen replied: This is a well known phenomenon that is observed for many transition metal particles when on the nanoscale. It will be a strong function of particle size. In Fischer–Tropsch catalysis it has been a long open question whether the deactivation of the catalyst for Co particles less than 6 nm in size is due to oxidation or intrinsic. The consensus now is that this deactivation is an intrinsic property of the small particles.

Cynthia Friend opened the discussion of the paper by Roy Johnston: Because your motivation for studying these systems was hydrogenation reactions, the titania support will be partially reduced under catalytic conditions. Prior work has shown that this leads to overgrowth of metal nanoparticles, including Rh (so-called strong metal support interactions (SMSI)). Did you consider such changes? If not, how would you approach this using theory? What methodology is required?

Roy Johnston answered: No we have not considered partial reduction of the titania support, though I agree this will be important in a future study of hydrogenation on AuRh catalysts. There should not be too much of a problem actually carrying out these calculations, though we would probably have to use a larger surface cell, so they will be more computationally expensive. However, generating configurations with overgrowth of partially reduced titania maybe more of a problem. In the absence of reliable empirical potentials to describe all of the required interactions, it may be necessary to carry out short DFT molecular dynamics simulations.

Aram Bugaev asked: Have you considered the effect of the thickness of the support on the bonding energies and geometries?

Roy Johnston answered: For the titania support we tested the convergence of the energy with the slab thickness. Reasonable convergence (when balanced against computational cost) was found for a slab of three TiO₂ (110) layers, corresponding to 9 layers of atoms. When calculating the surface binding energies of the AuRh and PdIr nanoparticles, the bottom TiO₂ layer was kept fixed, to mimic bulk TiO₂, and the top two layers were allowed to relax. However, we have not investigated how the cluster surface binding energy or geometry depends on the slab thickness. This is a good idea for a future study.

Hans-Joachim Freund queried: Have you considered charge transfer as a function of the distribution of Au and Rh with respect to the distance to the surface and in particular how that would be influenced by defects on the surface of TiO₂ or even below the surface.

Roy Johnston answered: No we have not investigated these aspects of charge transfer, but I agree that this would be a useful future study.

Philip Davies remarked: In a real catalytic system there is always a solvent present and the surface of the support will change accordingly. In particular, there is a lot of experimental evidence for the important role played by hydroxyls on the surface, see for example the review by Davis.¹ Have you considered the role of these species on the stability of the nanoparticles you have studied?

1 M. S. Ide and R. J. Davis, Acc. Chem. Res., 2014, 47, 825–833.

Roy Johnston responded: I agree that the explicit inclusion of solvent and solvent-induced changes to the substrate will be important when modelling catalysts under realistic operating conditions. So far, we have studied idealised nanoparticles and substrates, in order to establish the effect of chemical order (in the bimetallic nanoparticles) on cluster–substrate and cluster–adsorbate binding. In our future work, we plan to include solvent effects and surface modification.

Arunabhiram Chutia said: Your study on the interaction of Au–Rh and Pd–Ir on TiO₂ is really very interesting. Could you please comment on the distribution of charges on these nanoalloys due to their interaction with the TiO₂ surface, and secondly have you seen any electron transfer phenomenon from these nanolloys to the TiO₂ surface and vice versa?

Roy Johnston replied: In our calculations on titania-supported AuRh nanoparticles (ref. 24 in the paper), we observe electron transfer from the nanoparticle to the TiO₂ surface, which is greater when Rh (rather than Au) is bound to the surface. In the free AuRh nanoparticles, there is Rh to Au electron transfer. Due to the effect of the surface (mentioned above), for the supported AuRh nanoparticles, the Au–Rh charge separation increases significantly when Rh is in contact with the TiO₂ surface (for Janus-Rh and Au@Rh configurations), but there is little change when the Au is in contact with the surface (Janus-Au and Rh@Au). So far, we have not performed this analysis for the supported Pd–Ir nanoparticles.

Michele Carosso remarked: In your paper you suggest a strategy for stabilizing supported metal nanoparticles by adding a small amount of a second element, in order to increase the strength of the metal–support interaction. However, you studied this effect on a reducible, strongly interacting support such as TiO₂. Do you think that it could be possible to take advantage of this strategy also with a less-interacting support, such as, for example, an activated carbon? Do you expect that in these cases your nanoalloys will display a behavior intermediate between the unsupported and the TiO₂-supported ones?

Roy Johnston answered: Yes, I believe this is definitely possible. If the nanoparticle–substrate interaction is weaker then the effect on the nanoalloys (in terms of stabilising Janus or ball–cup structures) is indeed likely to be less than for the titania surface.

Valerii Bukhtiyarov commented: The aim of your paper is the understanding of elemental composition for different types of alloy catalysts, including surface composition. So I would like to ask did you analyse the ratio between metals on the catalyst surface in reaction conditions? Indeed, in the conditions of preferential CO oxidation, both CO and oxygen adsorbed on the catalysts can change the surface composition due to selective segregation of one of the elements.

Roy Johnston answered: Although we have not yet investigated thermal effects, for the 38-atom AuRh particles we found that the greater adsorption energies of CO and O₂ to Rh (as compared to Au) means that (in the presence of these molecules) configurations with some degree of Rh migration to the surface are lower in energy than the Rh@Au core–shell structure (ref. 23 in the paper).

Graham Hutchings remarked: The method of preparation you have used for your supported bimetallic nanoparticles is an impregnation technique and will produce atoms, clusters and nanoparticles which are evident in your micrographs. Do the clusters contain both metals? Or is there a minimum cluster size at which the second metal can be included or become stable? For some reaction such bimetallic clusters could be very effective and so accessing such structures could be useful if such structures can be readily made.

Roy Johnston replied: I believe that most of the clusters contain both elements. After heat treatment (700 ^oC) of the AuRh system, some pure Au and Rh NPs are indeed observed, in addition to Janus AuRh particles, with predominantly Rh at the nanoparticle–TiO₂ interface.¹ Since there is some degree of Au overgrowth at the edges of the Rh sub-clusters, these can be described as “ball–cup” configurations.

1 L. Piccolo et al., Sci. Rep., 2016, 6, 35226.

Chris-Kriton Skylaris asked: Have you included or examined including thermal effects in your calculations of binding of metal nanoparticles to titania? As they stand, binding energies are computed at 0 K while we know that in real applications we have finite temperatures and the binding is determined by the free energy. What are your thoughts about possible ways of including thermal and entropic contributions to your binding energies?

Roy Johnston responded: We haven't included any thermal effects in our calculations yet. I agree that, going forward, it will be important to calculate free energies to enable us to get closer to the experimental studies, which of course are at finite temperatures. The simplest approach would be to include vibrational contributions (at least in the harmonic approximation) to both the energy and the entropy, as these will lead to quantitative changes in the surface binding and adsorption energies, and (probably more importantly) may cause qualitative differences as regards to adsorption site preferences and kinetically preferred reaction pathways.

Stephen Shevlin asked: Is there a significant difference in the binding (i.e. atomic charge distribution) of your Janus nanoparticles depending on which face binds to the surface? Would this also have implications for catalytic properties?

Roy Johnston replied: We have calculated the overall metal-to-metal and cluster-to-support charge transfer for the Janus and core shell AuRh clusters (ref. 22 in the paper). For the Au₁₉Rh₁₉ Janus cluster, the Bader charges indicate that a total of 1.31 electrons are transferred from the Rh to the Au half of the cluster. For the Janus-Rh supported cluster (i.e. where Rh is in contact with the TiO₂ support), overall 2.50 electrons are transferred to the support, with the resulting total Rh and Au charges being +3.35 and −0.85, respectively. For the Janus-Au supported cluster, only 1.78 electrons are transferred to the support and the Au half of the cluster (which is now in contact with the support) has a small (+0.34) positive charge, with a higher charge (+1.44) on the Rh half. We have not yet analysed the charges on a facet-by-facet basis, but I believe that the surface charge distribution should indeed have an important influence on molecular adsorption and catalytic properties.

Lucas Garcia Verga remarked: I found this paper extremely interesting, especially the analysis about how the support affects different nanoalloys. Your results show that the interactions between metal nanoalloys and the support induce changes in the metal–metal bond lengths for metallic facets close to and far from the support. In the literature, these effects are usually followed by changes in the d-band centres and widths. These are useful electronic descriptors for the binding energies of reactive species such as O and CO in the metallic surface.

I understand that the focus of the work was to assess the stability of the isolated and supported nanoalloys; however, I was wondering if you calculated the shifts of the d-band centres for metallic facets close to and far from the support? If yes, do you see a trend between the shifts of the d-band centres and the binding strengths between nanoalloy and support?

Roy Johnston answered: We have carried out an analysis of the relationships between d-band centres and adsorption energies for CO and O₂ on free 38-atom Au–Rh clusters (ref. 23 in the paper) and for CO on free 38-atom Pd–Ir clusters (ref. 25 in the paper). In all cases, there is no simple correlation due to the importance of elastic (strain) effects in addition to the electronic effects. We have also investigated the relationship between d-band centre and adsorption energy for CO and O₂ on TiO₂-supported 38-atom Au–Rh clusters (ref. 24 in the paper). For the supported clusters, the elastic effects are reduced relative to the free clusters and there is a better correlation between d-band centre and adsorption energy. Going from the free to the supported clusters, there is a downward shift in the d-band centre, which is accompanied by a slight reduction in CO and O₂ adsorption energies. This is consistent with a net transfer of electron density from the nanoalloy to the support (as measured by Bader charges). Finally, I should say that we have calculated overall d-band shifts, not for specific facets, though this would be a good idea for future work.

Caetano Rodrigues Miranda asked: Can you rationalize your findings based on the electronic properties of the systems studied (difference charge densities, bands, size effects, etc.)?

Roy Johnston answered: In the Faraday Discussions paper and our other papers cited therein, we have performed analyses based on d-band filling, charge transfer and elastic (strain) effects and their contribution to the relative stabilities of free and supported pure and bimetallic nanoparticles and to the adsorption energies of small molecules. In some cases, we have found that these effects support each other, but sometimes they act against each other. In particular, surface binding and adsorption energies trends seem to be clearer for Au–Rh than for Pd–Ir.

Nia Richards asked: In your conclusion you use the term “Janus-like structure”. Does this imply that the structures you see are intermediate structures and not fully Janus structures?

Roy Johnston replied: Yes. For example, some of the AuRh nanoparticles are like Janus structures but often with Au overgrowth at the sides of the Rh part (but leaving Rh atoms in contact with the TiO₂ surface), giving rise to “ball–cup” structures, which are intermediate between core–shell and true Janus particles.

Cynthia Friend asked a general question to Rutger van Santen, Roy Johnston and David Willock: In modeling reactions, entropy is important to include. For complex reaction networks, low frequency modes need to be included and we need to go beyond the harmonic approximation. Can you all comment on what advances are needed to accurately include entropy?

David Willock answered: I agree, entropy is an important factor in chemical processes and is largely ignored in most theoretical approaches based on a transition state theory interpretation of potential energy surfaces. As you mention a common way to talk about reaction “Free energy” is to take these minima structures and transition states, perform a frequency calculation and the use the harmonic approximation to extract entropy changes. This is very approximate as low frequency modes have closely spaced energy levels which contribute significantly to the entropy. These are also the modes that have the greatest effect from non-harmonic effects. Entropy due to changes in translational and rotational degrees of freedom are also added in a general way based just on the mass and moments of inertia of the molecules involved in the reaction and the temperature.

What is needed are techniques that sample the immediate region around the key minima and transition states on the potential energy surface so that the number of states that are thermally accessible around each point can be estimated and so the entropy extracted directly. There are many more advanced methods that do this and that have been around for a number of years; Umbrella sampling,¹ transition path sampling,² metadynamics³ among many others. Each of these techniques use some form of molecular dynamics or Monte Carlo simulation to carry out the required sampling. They have been widely applied in biological systems and enzyme catalysis^4,5 with the use of QM/MM methods to speed up the sampling calculations.

These methods have also been applied to some homogeneous⁶ and heterogeneous catalysis reactions for example in ZnO catalysed methanol synthesis.⁷ These methods do require additional computational time and investment by the researchers to interpret the data in terms of reaction rates. Even so I would expect these methods to become increasingly important in the future.

1 S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen and P. A. Kollman, J. Comput. Chem., 1992, 13, 1011–1021.

2 P. G. Bolhuis, D. Chandler, C. Dellago and P. L. Geissler, Annu. Rev. Phys. Chem., 2002, 53, 291–318.

3 A. Barducci, G. Bussi and M. Parrinello, Phys. Rev. Lett., 2008, 100, 020603.

4 J. L. Gao and D. G. Truhlar, Annu. Rev. Phys. Chem., 2002, 53, 467–505.

5 K. Świderek, I. Tuñón, I. H. Williams and V. Moliner, J. Am. Chem. Soc., 2018, 140, 4327–4334.

6 A. Urakawa, M. Iannuzzi, J. Hutter and A. Baiker, Chem.– Eur. J., 2007, 13, 6828–6840.

7 J. Kiss, J. Frenzel, N. N. Nair, B. Meyer and D. Marx, J. Chem. Phys., 2011, 134, 064710.

Roy Johnston answered: Manzhos et al.¹ have reported the calculation of anharmonic vibrational frequencies and couplings for water on Pt(111) based on DFT calculations and solving the vibrational Schrödinger equation using a neural network. Perhaps such an approach could be used to fit parameters for cubic or quartic vibrational energy functions for cluster-adsorbed molecules.

1 S. Manzhos, T. Carrington, K. Yamashita, Surf. Sci., 2011, 605, 616–622.

Rutger van Santen answered: Since elementary reaction rates depend on activation free energies, it is essential to include properly calculated or estimated activation free energies. Partition functions can be used to calculate those. However the harmonic approximation can only be used for vibrational energies that are large compared to kT. Frustrated rotations are typical examples of modes where non harmonic corrections apply. Especially in zeolite catalysis this is a critical issue, since the intermediate carbenium ions are often nearly free moving. Molecular dynamics-related approaches, quantum-mechanical or quasi-classical have been fruitfully applied. Metadynamics is a quickly developing tool to address this issue.

Richard Catlow remarked: As I commented elsewhere, the landmark paper from Sauer’s group¹ calculated entropies and free energies including the contributions of anharmonic terms to the former.

1 Piccini et al., Angew. Chem., Int. Ed., 2016, 55, 5235.

Carlo Lamberti opened discussion of the paper by David Willock: Our experimental evidence on oxides,¹ zeolites² and MOFs³ support your theoretical prediction (DOI: 10.1039/C8FD00005K) that you need to have a reduced Cu(I) site to efficiently bond CO. It will be very interesting if you could extend your theoretical study on the overall redox cycle for CO oxychlorination in order to have a deeper understanding of the structure of the oxychloride phase, the formation of which has been foreseen after interaction of the reduced form of the catalyst with oxygen.^1,4–8 Finally, as your catalyst contains 10% CuCl₂ and also 8% KCl, if you want to have a realistic picture of its redox property, you should include potassium in your model, because it is known that its presence strongly favors the oxidized state of copper chloride.^4,6,7

1 G. Leofanti, A. Marsella, B. Cremaschi, M. Garilli, A. Zecchina, G. Spoto, S. Bordiga, P. Fisicaro, G. Berlier, C. Prestipino, G. Casali and C. Lamberti, J. Catal., 2001, 202, 279–295.

2 F. Giordanino, P. N. R. Vennestrøm, L. F. Lundegaard, F. N. Stappen, S. Mossin, P. Beato, S. Bordiga and C. Lamberti, Dalton Trans., 2013, 42, 12741–12761.

3 L. Braglia, E. Borfecchia, K. A. Lomachenko, A. L. Bugaev, A. A. Guda, A. V. Soldatov, B. T. L. Bleken, S. Øien-Ødegaard, U. Olsbye, K. P. Lillerud, S. Bordiga, G. Agostini, M. Manzoli and C. Lamberti, Faraday Discuss., 2017, 201, 265–286.

4 C. Lamberti, C. Prestipino, F. Bonino, L. Capello, S. Bordiga, G. Spoto, A. Zecchina, S. D. Moreno, B. Cremaschi, M. Garilli, A. Marsella, D. Carmello, S. Vidotto and G. Leofanti, Angew. Chem., Int. Ed., 2002, 41, 2341–2344.

5 G. Leofanti, A. Marsella, B. Cremaschi, M. Garilli, A. Zecchina, G. Spoto, S. Bordiga, P. Fisicaro, C. Prestipino, F. Villain and C. Lamberti, J. Catal., 2002, 205, 375–381.

6 N. B. Muddada, U. Olsbye, L. Caccialupi, F. Cavani, G. Leofanti, D. Gianolio, S. Bordiga and C. Lamberti, Phys. Chem. Chem. Phys., 2010, 12, 5605–5618.

7 N. B. Muddada, U. Olsbye, G. Leofanti, D. Gianolio, F. Bonino, S. Bordiga, T. Fuglerud, S. Vidotto, A. Marsella and C. Lamberti, Dalton Trans., 2010, 39, 8437–8449.

8 T. Zhang, C. Troll, B. Rieger, J. Kintrup, O. F.-K. Schlüter and R. Weber, J. Catal., 2010, 270, 76–85.

David Willock responded: It is good to hear that the experimental evidence from your experimental work also shows this requirement to have Cu(I) present to adsorb CO. We are currently working on models of the chlorination process itself, i.e. the transfer of Cl from the lattice to the adsorbed CO. We have also created some higher index planes which require termination with water or OH groups so that the competition between phosgene synthesis and oxidation to CO₂ can be modelled. The introduction of KCl would be a separate study; we have no structural model for the location of the KCl and the way the two chlorides are mixed. This means that constructing a reliable model with KCl present is difficult. A starting point may be to simply dope the CuCl₂ lattice with K⁺ and see the effect on the defect formation energies for the Cl defects that we have presented here.

Paul Sermon remarked: Your paper mentions CuCl₂/alumina ethene oxychlorination catalysts. You characterize your CuCl₂/attapulgite catalyst by CO conversion to Cl₂C=O at 633 K. Your results reminded me of the ethene (1kPa)/He temperature-programmed titration (from 298–773 K at 5 K min⁻¹ (i.e. below the melting point of bulk CuCl₂)) of 100 mg CuCl₂/alumina, silica and titania and PdCl₂–CuCl₂/titania catalysts by Keith Rollins.¹ This revealed different peaks (T_max) of maximum rates of 1,2-ethylene dichloride (EDC) and vinyl chloride (VCM) production, along with integrated numbers of EDC molecules produced overall, that varied with the support and the addition of Pd. Might one be able to titrate your catalysts with CO and see maximum rates of phosgene and CO₂ production (kinetically limited at low temperature and thermodynamically limited above 473 K)? Does the CO₂ come from a shift reaction?

1 K. Rollins and P. A. Sermon, J. Chem. Soc., Chem. Commun., 1986, 1171–1173.

David Willock responded: This is a very good point. In the current paper we have concentrated on the structure of the material used and the removal of Cl from the CuCl₂ lattice during phosgene production. We are also working now on a publication which covers our reactor work in more detail. As you suggest we can titrate the Cl active site with CO monitoring the reduction of phosgene production as a function of time. CO₂ activity follows a similar trend although we have not checked if we can quantify the number of sites involved with CO₂ production and so rule out a water gas shift reaction. Even so the levels of water in the gas feeds are kept as low as possible during these reactions; as we monitor the products with FTIR we can confirm that there are no significant bands in the water stretching region of the spectra.

Graham Hutchings said: You have carried out a very interesting combined experimental and theoretical study and your conclusion is that a Cu(I)–Cu(II) redox cycle is active. Copper(II) chloride at 370 °C will have an appreciable vapour pressure; it was the catalyst in the original Deacon process and is known to deactivate rapidly. Can you predict or suggest alternative metals to copper that could be more stable? For example, ruthenium oxide is a more stable Deacon catalyst.

David Willock replied: We would point out that the reaction temperature used here is lower than that for the Deacon process which operates at 400–450 °C. The commercial catalyst used in our experiments also contains KCl which is thought to act as a stabiliser for the CuCl₂ supported on the clay. Even so we would agree that stability would need to be carefully tested for long term application of the catalyst for phosgene production. It would be interesting for us to carry out calculations on the ruthenium chloride system too for comparison. In our latest calculations we are examining the elementary steps that lead to CO₂ or CCl₂O over a higher index termination of the CuCl₂ structure terminated with a mixture of OH and Cl. This will allow the selectivity of the catalyst to be examined. It would be interesting to find a chloride that was more selective to phosgene over carbon dioxide.

Andrea Russell remarked: As an electrochemist, I recognise that Cu⁺ is not a stable species, with the reaction 2Cu⁺ → Cu + Cu²⁺ being spontaneous. In your paper you present a conundrum in that you observe the Cu species being oxidised upon the addition of CO, which is more normally thought of as a reducing agent. Do you think that it is really the spontaneous reaction between two Cu⁺ ions that is occurring, which becomes possible as the Cl is consumed in the reaction? It appears you rejected this idea in your paper as you did not observe much Cu⁰, but it looks to me that the post edge features in the XANES may be indicative of this species. I don’t think that a Cu foil is necessarily the best reference for Cu⁰ in this case, as the local coordination environment also influences the XANES features.

David Willock responded: The stability of Cu⁺ will depend on reaction conditions; we agree that in the aqueous chemistry of an electrochemical cell this disproportionation will take place. Indeed we thought about this when trying to understand the apparent Cu oxidation on introduction of CO. However we saw no evidence of Cu⁰ in our data and in the overlayered XANES spectra of Fig.5a in the paper we note an isobestic point which would suggest direct interconversion of Cu⁺ and Cu²⁺ without the generation of any Cu⁰. This is a high temperature gas/solid reaction and we know that in the solid state Cu⁺ can be stabilised, for example in the synthesis of Cu₂O.

Keith Whiston asked: What is the effect of the clay support on the reducibility and performance of the CuCl₂ catalyst? Does the KCl modifier used in the commercial catalyst inhibit copper reduction or otherwise improve the lifetime or performance of the reaction?

David Willock responded: We have not studied the role of the support and KCl modifier ourselves as we have only presented results for the commercial catalyst. Experiments making a direct comparison of alumina supported CuCl₂ with and without KCl modifiers have been reported.¹ These show that KCl acts to stabilise the higher Cu oxidation state chloride which allows higher temperature operation of the catalysts.

1 C. Lamberti, C. Prestipino, F. Bonino, L. Capello, S. Bordiga, G. Spoto, A. Zecchina, S. D. Moreno, B. Cremaschi, M. Garilli, A. Marsella, D. Carmello, S. Vidotto and G. Leofanti, Angew. Chem., Int. Ed., 2002, 41, 2341.

Francesca Baletto remarked: In your paper you have shown calculations using two different versions of U, namely 4 and 7, which provide considerably different band gaps of 0.3 and 0.9 eV, respectively. Could you comment on the U effects on surface defect formation, binding energy of CO and charge transfer?

David Willock answered: We have chosen to look at U = 4 as a lower limit of the parameter at which a clear band gap appears in the calculation and U = 7 which seems to be more widely used in the literature for calculations on Cu salts. The U = 7 parameter tends to give defect formation energies around 0.1 eV higher than U = 4. However, the trends on comparing the surface defect formation with second layer defects and comparing the small and large supercell results are the same irrespective of the choice of U. So both sets of data show the second layer defect formation energy around 0.1 eV higher than that for the surface layer and the larger supercell giving lower defect formation energies (by up to 0.06 eV). We rationalised this by thinking about the electronic character of the defect. When a Cl⁻ ion is removed as 1/2 Cl₂ the electron remaining will reduce one Cu centre. The U parameter ensures that this electron is localised on one of the Cu centres neighbouring the defect, while a calculation without a U correction would tend to delocalise the electron in the conduction band. It appears that once the U parameter is large enough to introduce a band gap the electron localisation at a Cu centre is ensured and so the results are only relatively weakly affected by the choice of U. For the adsorption of CO the calculated energies also seem to be only weakly affected by the choice of U.

Carlo Lamberti said: I would like to add two comments here. First, in all our studies on the ethylene oxychlorination reaction,^1–9 we never observed evidence of a measurable fraction of Cu(0) species. Second, there is evidence suggesting that CuCl₂ and CuCl₂/CuCl supported catalysts should be in the form of a molten salt under oxychlorination reactions. This holds for both ethylene and CO oxychlorination, that are performed around 200 and 370 °C, respectively. This is obviously very difficult to prove on a structural point of view. Indeed, in our experience the 10 wt% CuCl₂ loaded catalyst on γ-alumina, even at room temperature, exhibits a CuCl₂ phase that is highly dispersed and probably of amorphous nature, as it has been well detected by EXAFS, being however XRD silent. Probably an in situ PDF study would be required to fully understand this point.

1 G. Leofanti, M. Padovan, M. Garilli, D. Carmello, G. L. Marra, A. Zecchina, G. Spoto, S. Bordiga and C. Lamberti, J. Catal., 2000, 189, 105–116.

2 G. Leofanti, M. Padovan, M. Garilli, D. Carmello, A. Zecchina, G. Spoto, S. Bordiga, G. T. Palomino and C. Lamberti, J. Catal., 2000, 189, 91–104.

3 G. Leofanti, A. Marsella, B. Cremaschi, M. Garilli, A. Zecchina, G. Spoto, S. Bordiga, P. Fisicaro, G. Berlier, C. Prestipino, G. Casali and C. Lamberti, J. Catal., 2001, 202, 279–295.

4 G. Leofanti, A. Marsella, B. Cremaschi, M. Garilli, A. Zecchina, G. Spoto, S. Bordiga, P. Fisicaro, C. Prestipino, F. Villain and C. Lamberti, J. Catal., 2002, 205, 375–381.

5 C. Lamberti, C. Prestipino, F. Bonino, L. Capello, S. Bordiga, G. Spoto, A. Zecchina, S. D. Moreno, B. Cremaschi, M. Garilli, A. Marsella, D. Carmello, S. Vidotto and G. Leofanti, Angew. Chem., Int. Ed., 2002, 41, 2341–2344.

7 C. Prestipino, S. Bordiga, C. Lamberti, S. Vidotto, M. Garilli, B. Cremaschi, A. Marsella, G. Leofanti, P. Fisicaro, G. Spoto and A. Zecchina, J. Phys. Chem. B, 2003, 107, 5022–5030.

8 N. B. Muddada, U. Olsbye, L. Caccialupi, F. Cavani, G. Leofanti, D. Gianolio, S. Bordiga and C. Lamberti, Phys. Chem. Chem. Phys., 2010, 12, 5605–5618.

9 N. B. Muddada, U. Olsbye, G. Leofanti, D. Gianolio, F. Bonino, S. Bordiga, T. Fuglerud, S. Vidotto, A. Marsella and C. Lamberti, Dalton Trans., 2010, 39, 8437–8449.

Nia Richards asked: Have you investigated the effect of chlorine concentration in the pre-treatment stream, and how does this effect phosgene selectivity?

David Willock replied: In the laboratory experiments for testing phosgene production (see Fig. 1 in the paper), we have tried different concentrations of Cl₂ pre-treatment. The more and the longer Cl₂ is passed, the more phosgene (and less CO₂) is observed. Eventually, the saturation of Cl in the clay was reached and the amount of phosgene produced became constant. In the samples prepared for XANES analysis this saturation level of Cl₂ was used.

Mzamo Shozi asked: Could electron spin resonance be used to detect formation of chlorine radicals during pre-treatment of the catalysts with Cl₂ gas?

David Willock responded: Electron spin resonance would be useful in this area to show how chlorine is stored in the material. We see the formation of CuCl₂via X-ray diffraction and the XANES data of the chloride catalyst. However, our observations of the oxidation state of Cu as CO flow is introduced may suggest that there are other Cl species on the catalyst and the way that Cl₂ is taken up by this material would be an interesting study in its own right.

Richard Catlow continued the discussion of the paper by Rutger van Santen: How far can we simulate full reaction cycles? There are many examples in the literature and there is no doubt that the field has made great progress in recent years, but can you comment on how reliable the quantitative aspects of the results are?

Rutger van Santen replied: Once the mechanism of the catalytic reaction has been formulated and on this basis elementary reaction rate constants have been computed and the catalytic reaction cycle has been closed, the ordinary differential equations can be formulated that enable us to calculate kinetics of the reaction. There is no automatic procedure to establish the reaction mechanism. One needs to use available experimental or additional computational information to make a proposal and several alternatives have sometimes to be included in the calculations. They may be operated in parallel. It is essential then to, as we do, solve the corresponding microkinetics equations without making an assumption on a rate controlling step. This should come automatically from the simulation. Implicit to this procedure is the assumption of the mean field approximation. In case reaction intermediates show an inhomogeneous surface distribution, the ODE’s may have to be replaced by correponding Monte Carlo simulations and surface diffusion has to be explicitly taken into account.

When direct comparison is made with experiment, a decision has to be made whether simulations should be done at differential or integral conditions. In the latter case, convection as in reaction engineering simulations has to be at least taken into account. Also, in relation to experiment, one has to be aware that dependent on the evolution of the reaction, the surface state may be close to the initial state, may have converged to a steady state configuration or composition, or is a in a state of deactivation. Theory is available that is able to predict the state of a surface in equilibrium with a reactive medium. It may be required to establish this separately from the kinetic simulations because the timescale of surface equilibration or reconstruction may be long compared to that of the catalytic cycle. When simulations are completely “ab initio” reliability of the calculations will depend on the accuracy of the used molecular information that is contained in the elementary reaction rate constants.

Generally when DFT simulations are used, activation energies have at least an error of 10 kJ mol⁻¹ and pre-exponents may have substantial errors if only calculated within the harmonic approximation. Care has to be taken that elementary reaction equilibria are correct. Systematic energy errors usually cancel out (except between gas phase and surface), but this will not be the case for reaction entropies. So the use of proper pre-exponents of the elementary reaction rates is essential.

To predict an overall reaction rate properly the temperature of reaction has to be predicted right. This very often depends on the equilibrium of a molecule between the gas phase and the surface. When based on DFT this equilibrium has to be usually adjusted to experiment. In the case of the zeolite simulations we discussed in our paper, the adsorption isotherms of propylene and isobutane based on experimental values took care of this. An additional issue is the question whether the interaction between adsorbates or reaction centers can be considered ideal and lateral interactions can be ignored. In addition to concentration dependent correction terms to reaction energies this may also lead to surface reconstruction effects and island formation. This may be a strong function of conditions.

Clearly absolute catalyst performace prediction by full kinetics simulations have to be considered with care. However, when used to predict trends as a function of surface reactivity, as when Sabatier volcano’s are constructed, due to cancellation of systematic errors such simulations may provide considerable insight into the microscopic interactions that determine activity or selectivity differences.

Caetano Rodrigues Miranda asked: How do the variation and errors at DFT level (accuracy and functional dependency) affect and propagate within the microkinetic calculations? Are there studies in this direction on how sensitive or robust the microkinetics model results are regarding the variation of DFT ones?

Rutger van Santen answered: See also my reply to the previous question. The general comment is that for surface reactions that can be considered quasi-equilibrated, systematic errors will largely cancel. This is not the case for equilibria between the gas phase and the surface. Critical also are the free energy values used for those reaction steps that are rate controlling. When adsorption equilibria between gas and surface determine largely surface vacancies, such as for low temperature reactions or with reactants that strongly adsorb, the temperature of reaction will strongly depend on adsorption free energies. Then it is advisable to use experimental data or computed data of higher accuracy than DFT calculations can provide.

Bruce Gates commented: You mentioned the importance of curvatures or shapes of zeolite pores in the modeling of your hydrocarbon reactions. In a recent report from the U.S. DOE (Basic Research Needs for Catalysis), a prominent recommendation was for research about the environments immediately surrounding catalytic sites. Please let us know your thoughts about how important this issue is and whether you have some suggestions about how to formulate questions about it.

Rutger van Santen replied: In zeolite catalysis the size and shape of zeolite cavities play an important role, because they determine the strength of the van der Waals interaction with occluded molecules. The adsorption isotherms of reactant and product molecules are a sensitive function. Steric matches of reactant structure and size and zeolite nanopore dimensions are relevant. This affects rates of diffusion as well as reaction. The curvature of the zeolite cavity inhibits sterically extended intermediates to become close to catalyst reaction centers. In acid catalysis this will strongly affect activation energies for protonation and deprotonation and relative stability of carbenium ions versus alcoxy species.

Medium effects due to the presence of a high concentration of molecules in the zeolite micropore may be important and act quite differently from comparable effects in solvents or solution. Medium effects in apolar media are well understood in high pressure hydrocarbon conversion catalysis¹ where it has been demonstrated in hydrocracking catalysis that packing of hydrocarbon fragments in the zeolite, that equilibrate, determine the selectivity of the reaction. The interplay between polar solvent molecules as water and proton catalysis is physico-chemically complex. It relates for instance to the difficulty to predict computational prediction of the pH of an acid.

It is known from enzyme catalysis that the presence of a few water molecules in a hydrophobic environment will have a dramatic effect on rates of proton transfer reactions. Proton channeling through the water proton bridges has been demonstrated to play an important role in redox reactions as the Wacker reaction or in the selective oxidation reaction of glucose and related molecules in zeolites.

Computational complexity arises due to the need to combine through molecular dynamics the mobile adjustments of solvent molecules around the reacting complex with calculations of (partially) ionic transition states or reaction intermediates.

1 B. Smit et.al., Chem. Rev., 2008, 108, 4125.

Katharina Brinkert asked: Could you comment on the limitations of DFT calculations/simulations for catalysis? Which theoretical tools do we need to describe catalytic systems in a better way? Which information would you like to obtain from experimentalists in this respect?

Rutger van Santen replied: There is a great variety in DFT computational methods in the way they deal with functionals or basis functions. The advantage of the technique is that many now allow for calculation on complex systems close to those in practice of reaction intermediates as well as transition states in combination with molecular dynamics simulations around local energy minima or maxima. A drawback is the still substantial error in energy accuracy, that is typically still at least 10 kJ mol⁻¹ and that it is not trivial to obtain calculated pre-exponents of elementary reaction rates beyond the harmonic approximation. Those are needed for systems with frustrated bending or rotational frequencies or other low frequencies related to the reaction coordinate.

The aim should be to predict properly the spectroscopic properties of adsorption intermediates (less challenging and often doable with spectroscopic accuracy, unless one deals with highly electron correlated systems as the oxides or sulfides. Here there is a need for detailed information of electron structure. High quality first principle calculations embedded in larger matrices are probably the solution) or elementary reaction rates or adsorption energies (for the latter two, experimental numbers of high quality are highly needed.

Richard Catlow opened a general discussion of the papers by Rutger van Santen, Roy Johnston and David Willock: It is important to stress the progress that has been made in the application of modelling techniques in catalytic science. A notable recent development was the landmark paper from Sauer’s group¹ which calculated rate constants of catalytic reactions within zeolites with chemical accuracy. These techniques are far from routine, but they illustrate what can now be achieved.

1 Piccini et al., Angew. Chem., Int. Ed., 2016, 55, 5235.

Roy Johnston commented: At present, it is feasible to perform computational studies at a higher level of theory (e.g. coupled cluster and Time Dependent DFT) on small clusters in the gas phase, which can be compared with experimental UV-vis and IR spectroscopy (generally involving photodepletion coupled with mass spectrometry), and magnetic or electrostatic deflection measurements. With the development of faster computers and more efficient computer codes (e.g. linear scaling DFT), it will be possible to extend these methods to larger clusters and nanoparticles.

Carlo Lamberti asked: Concerning the theoretical work of Sautet et al.,¹ predicting the coverage-dependent reshaping of a 13-atoms Pt cluster supported on γ-Al₂O₃ in the presence of different numbers of adsorbed H atoms, it is worth mentioning that at the Operando VI conference Prof. E. Groppo presented work² where she showed synchronous IR (in DRIFT mode), XANES/EXAFS (in transmission mode) and MS data supporting on an experimental ground the theoretical predictions of Prof. Sautet.¹

1 C. Mager-Maury, G. Bonnard, C. Chizallet, P. Sautet and P. Raybaud, ChemCatChem, 2011, 3, 200–207.

2 E. Groppo, Dynamics of reactive species and reactant-induced reconstruction of Pt clusters in Pt/Al₂O₃catalysts, presented at Operando VI conference (Estepona, Spain, April 15–19, 2018).

Rutger van Santen replied: This is very nice confirmation of agreement between state of the art computational prediction of the state of small transition metal nanoparticles at ambient conditions and experiment.

Rosa Arrigo commented: With reference to the challenge of electrocatalyst prediction and in particular to the case of the carbon dioxide electrochemical reduction, some of the computational studies present in the literature^1,2 use the binding energy of for instance carbon monoxide to the surface as a reactivity descriptor. Whilst this static description explains the reactivity towards the formation of C1 products such as carbon monoxide and methane, to explain the formation of C3 molecules using this model one intuitively would invoke a more complex structure of the active sites. To complicate things further, in the case of metals such as Fe³ which is able to dissolve C and form carbides (as opposed to Cu), it could well be that metastable subsurface carbide could be involved in the reaction mechanism. Would it be possible by means of the computational tools available nowadays to model such surface/subsurface dynamics, which are possibly driven by kinetic factors rather than thermodynamics and link this to the products evolved?

1 K. P. Kuhl, T. Hatsukade, E. R. Cave, D. N. Abram, J. Kibsgaard and T. F. Jaramillo, J. Am. Chem. Soc., 2014, 136, 14107–14113.

2 A. Bagger, W. Ju, A. S. Varela, P. Strasser and J. Rossmeisl, ChemPhysChem, 2017, 18, 3266–3273.

3 R. Arrigo, M. E. Schuster, S. Wrabetz, F. Girgsdies, J.-P. Tessonnier, G. Centi, S. Perathoner, D. S. Su and R. Schloegl, ChemSusChem, 2012, 5, 577–586.

Rutger van Santen responded: For Fischer–Tropsch catalysis there is no correlation between chain growth selectivity and CO adsorption energy. The electrochemically active system that is most selective with respect to hydrogen evolution is the Cu electrode. Also in this case there is no correlation with the CO adsorption strength. The electrocatalytic mechanism that produces longer hydrocarbons is substantially different from that of the Fischer–Tropsch reaction.

As I discuss in the paper, molecular dynamics simulations based on fitted forcefields can be done to study the dynamics and surface reconstruction that result from high adatom coverage and subsurface adatom incorporation.

Graham Hutchings remarked: Returning to the complexity of supported catalysts when made by deposition precipitation it is apparent that atoms, clusters and nanoparticles are all present. The support will have a different influence on each of these as one can envisage a single atom being affected more so than a nanoparticle. Is theory now at a level that it can comment on the reactivity of these different species so we can refine what catalysts we prepare?

David Willock responded: Computer simulation can build models of the structures that we think are important in the catalysis. We have seen in this conference examples of single metal atoms on oxide supports, isolated clusters of metal atoms and discussed the chemical composition of particles in response to their environment such as the oxidation of Pt particles. One advantage of modelling is that the structure that is used to represent the “active” site is well defined and we usually only work on one particular type of structure at a time. But by building different models and testing the energetics of catalysed reactions over each model a comparison can be made. This could be used to narrow down which is the more active, the single atoms, the nanoparticles or even nanoparticle/support interface sites.

Roy Johnston commented: I believe that this is possible provided the various-sized systems can all be calculated at the same level of theory and using the same functionals, basis sets etc. For DFT calculations, this is a realistic goal due to the development of linear scaling DFT (ONETEP).

Rutger van Santen replied: For model systems that are well defined, such as ideal surfaces of non reducible oxides like zeolites, but also reducible oxides, such as CeO₂ or TiO₂, to calculate the state and relative energies of adsorbed metal atoms, small clusters and even small nano particles are feasable. Also for non ideal model surfaces, when hydroxylated or containing vacancies this is doable .

The real issue I believe is to predict the state of the surface or particles at particular stages of catalyst preparation. Often a complex solvent is present etc. This is not only computationally a challenge but also relates to a proper understanding of the physical chemistry and inorganic chemistry of these complex systems and critical conditions for particular transformations.

Hans-Joachim Freund said: We need to be careful in making statements about predictability of structures based purely on calculations of the ideal system because the state of the support is complex.

David Willock answered: I do agree that the support materials used in catalysis can be complex with different degrees of defects present, stepped surfaces and reactions with the environment that change the surface chemically, for example hydroxylation. This gives us a drive to work with experimentalists to use characterisation to understand the likely chemical state of the surface and probe the structure of the surface. What calculations then do is to link the structural characterisation to the observed reactivity and to attempt to understand what the important features of the surface are that lead to catalytic activity.

Wilke Dononelli commented: Related to the question that was asked about the accuracy of theoretical calculations and the use of different levels of theory, I want to add that the theory that has to be chosen depends on the investigated problem. For example DFT might predict very good structures compared to experiment. It could also give good energies compared to experiment. But in some cases it could also fail to predict accurate energetics. In order to predict such precise energies it might be a good choice to go beyond DFT to higher levels of theory.

Richard Catlow replied: I fully agree with your comment including the need to go beyond DFT to higher levels of theory; and perhaps the coupled cluster approach offers an opportunity in this respect.

Cynthia Friend asked: There has been a focus on studying materials structure using theory more so than studying reactions in the first session. While this is important, catalysis depends on predicting reactivity and kinetics, meaning that complex reaction schemes and knowledge of elementary steps are required. As noted in my previous question, entropy is also very important. This was demonstrated in Prof. van Santen’s paper. He had a network with 140 steps. Confinement in pores of the zeolite was important in the transition state. How can experimental work help constrain this problem so it is tractable? What theoretical advances are necessary to better address these issues?

Rutger van Santen responded: See my reply to your earlier question.

Julien Marbaix opened a general discussion of the paper by Nora de Leeuw: As we know that the cathode in SOFC can be divided into two layers, conduction and reaction, did you try to integrate the incoming oxygen flow to understand its influence on the conductive layer performance?

Nora de Leeuw answered: We did not study the oxygen migration within the YSZ material as we were interested mainly in the geometric and electronic structures of the Ni clusters on top of the oxide. Here, we focused on the behaviour of the metal atoms: whether they prefer to aggregate or wet the surface, their mobility on the surface, and their interaction with the YSZ surface. However, there are other investigations where the authors have performed a theoretical study of the fuel cell, including the migration of oxygen through the electrolyte.^1,2

1 X. Wang, K. C. Lau, C. H. Turner and B. I. Dunlap, J. Power Sources, 2010, 195, 4177–4184.

2 S. C. Ammal and A. Heyden, J. Phys. Chem. Lett., 2012, 3, 2767–2772.

Rutger van Santen commented: Concerning the Monte Carlo simulations of nickel cluster formation, I have a question about whether the elementary rate constants used are reversible. Is there then a relationship between agglomerisation time and the overall thermodynamics of the process?

Nora de Leeuw responded: We have performed the kinetics simulation without any restriction in the sense that all the elementary reactions were reversible and individual rate constants were calculated independently for each process. The rate constants are derived from the thermodynamics and they are linked.

Bruce Gates asked: What distances do you find between Ni and O atoms on the zirconia support surface, and what happens if the zirconia is hydroxylated?

Nora de Leeuw responded: The average distance between Ni clusters and the support surface, in Ni_n/ZrO₂(111) and Ni_n/YSZ(111), is 1.9 Å. We have not considered a hydroxylated surface, which might affect the Ni binding to the surface. This could be taken up in a future study, but here we had to start from a simpler case in order to gain initial insight into the Ni clustering on the zirconia and YSZ surfaces.

Parag Shah queried: As the surface is heated up during the simulations, does the Y segregate or move within the surface?

Nora de Leeuw responded: We have not considered segregation effects in the support. We have considered different positions for the Y atoms and noted that, in the most stable configuration, Y is positioned at the top of the surface and as the next nearest neighbour of the oxygen vacancy. We have provided more details about the Y position in YSZ in our previous paper.¹

1 A. Cadi-Essadek, A. Roldan and N. H. de Leeuw, J. Phys. Chem. C, 2015, 119, 6581–6591.

Andrea Russell asked: In your paper you showed that the rate of sintering of the Ni particles was independent of coverage. Do you think that this observation can be accounted for by the fact that you’ve only considered a high coverage regime and perhaps the rate is only pseudo zero order and may become dependent on coverage at lower overall Ni coverage?

Nora de Leeuw answered: It is correct that we have neglected the lateral interactions of non-bonded Ni, which is a fair approximation as two nickel atoms form a bond from a relatively large distance. We have considered low initial coverage (5%) for Ni₁₀ (Fig. 10(c) in the paper) and found that the sintering rate is similar to an initial coverage of 10% (Fig. 10(d) in the paper). For a low atomic coverage, a single Ni atom might be too far to bind to another structure. We could have included random movements for further evaluation but that was outside the scope of this paper.

Parasuraman Selvam asked: What will happen to the stability and mobility of supported nickel clusters if we start with regular shaped clusters such as tetrahedron (Ni₄), octahedron (Ni₆), icosahedron (Ni₁₂) or cuboctahedron (Ni₁₂) or anticuboctahedron (Ni₁₂) as the starting geometry rather than 1–10 atoms? In fact, in the realistic situation, we get regular shaped clusters/particles, viz., spherical, cubic or isosahedron/cuboctahedron/anticuboctahedron with 2–4 nm sizes onto the supported system. Is it not appropriate to start with such clusters for a rational understanding of the metal–support interaction?

Nora de Leeuw responded: We built our Ni clusters’ shape by cutting the Ni(111) surface; a similar approach to previous theoretical studies.¹ Thus, the most stable surface of the cluster is facing the gas phase. We could have also tried other geometries such as tetrahedron (Ni₄) but the result would have been the same, i.e. that a 3D shape will be more stable than a flat configuration. This is the main outcome of our study: Ni atoms prefer to aggregate rather than wetting the surface.

1 M. Shishkin and T. Ziegler, J. Phys. Chem. C, 2009, 113(52), 21667–21678.

David Willock remarked: In the transition states you show in Fig. 7 of your paper nickel atoms are moving to join a cluster. It looks like the barrier is to do with diffusion of the Ni atom over the surface rather than a barrier to it joining the cluster. Did you look at different directions of approach for the Ni atom across the surface of the oxide? If so, are there easy directions of travel for Ni atoms that will make the formation of large clusters along certain crystallographic directions easier than others?

Nora de Leeuw responded: A single Ni atom will move across the surface following random movements before finding another Ni atom. The symmetry of the ZrO₂ structure would lead to isotropic paths, i.e. same energies. We have considered this fast Ni diffusion on the surface, but have not noted asymmetric growth in certain directions; perhaps because the clusters were as yet too small to show such a phenomenon.

Hans-Joachim Freund asked: Have you looked at oxygen vacancies on the clean ZrO₂(111) or more reactive surfaces (with respect to recent calculations^1,2 by G. Pacchioni on nano-ZrO₂) in comparison to Y stabilized ZrO₂?

1 A. R. Puigdollers et al., J. Phys. Chem. C, 2016, 120, 15329–15337.

2 A. R. Puigdollers et al., J. Phys. Chem. C, 2016, 120, 4392–4402.

Nora de Leeuw answered: We did study various oxygen vacancies in ZrO₂ and YSZ,¹ however, we have not considered them in this work, which main goal was to understand the Ni adsorption on YSZ(111) surface.

1 A. Cadi-Essadek, A. Roldan and N. H. de Leeuw, J. Phys. Chem. C, 2015, 119, 6581–6591.

Michele Carosso commented: You have shown that the Ni_n clusters supported on both ZrO₂ and Y-modified ZrO₂ present a pyramidal shape where the base is characterized by a partial positive charge while the apex is either charge-neutral or partially negative. This of course is of great interest for catalysis because adsorption at the cluster surface of both electrophilic and nucleophilic substances could be possible at the same time. Is this effect stronger for one of the two supports? Do you expect that a similar effect is present also in combination with other supports?

Nora de Leeuw replied: The effect is similar for both ZrO₂ and YSZ supports: the average Bader charge of the Ni atoms at the base is +0.1 e while the charge of the Ni atoms at the apex is −0.1 e or nil. Considering other combined supports, as long as the adsorption of the Ni clusters is not favourable, i.e. aggregation and formation of a Ni pyramid, we should observe a similar charge distribution to the one observed in ZrO₂ and YSZ. However, if the Ni cluster adsorption is favourable, i.e. wetting of Ni atoms over the surface, we should observe stronger charge transfer from the metal atoms to the surface. In the latter case, all the adsorbed Ni atoms would have a positive charge.

Wilke Dononelli commented: Your results are very interesting, especially the different structures you found for different sizes of nanoparticles on the supports. Giordano et al. found a similar tetrahedral structure for Ni₄ nanoparticles on MgO.¹ In addition they found a flat configuration. For gas phase Pt₄ clusters, Demirogulo et al. found a slightly bent rhombus configuration as the global minimum configuration and a planar rectangular configuration for Ru₄.² Do you think that there might be other dominant configurations on your investigated supports, depending on the adsorption position? It would additionally be very interesting to see how the structure of the nanoparticles change when they are used for catalysis and reactants are adsorbed. Did you, for example, look at the change in structure of the tetragonal Ni₄ or pyramidal Ni₁₀ nanoparticles, when these are covered with oxygen?

1 L. Giordano et al., Surf. Sci., 2001, 473, 213–226.

2 I. Demiroglu, K. Yao, H. A. Hussein and R. L. Johnston, J. Phys. Chem. C, 2017, 121, 10773–10780.

Nora de Leeuw answered: In our first investigation¹ we scanned all the different adsorption sites of Ni on top of both ZrO₂(111) and YSZ(111). We noted that on ZrO₂(111), Ni prefers to adsorb on top of O_d which is the oxygen belonging to the 3rd atomic layer (Fig. 1(a) in the current paper), slightly off the perpendicular. On YSZ(111), the preferential adsorption site is on top of the vacancy and away from the Y atoms. We therefore decided to build Ni_n clusters around the oxygen vacancy. Then, we considered different cluster shapes. We agree that it would have been interesting to evaluate other adsorption sites to see how the cluster shape would be affected but we chose to focus on the most stable geometries, as the next step of our study will be an investigation of the reactivity at the interface between the cluster and the YSZ surface.

1 A. Cadi-Essadek, A. Roldan, and N. H. de Leeuw, J. Phys. Chem. C, 2015, 119, 6581–6591.

Carlo Lamberti opened discussion of the paper by Wilke Dononelli: There are several examples where the adsorption of carbon monoxide on a surface site exhibits an equilibrium between C-end and O-end adducts: M⋯CO ⇄ M⋯OC, see e.g. the examples reviewed in Table 9 of ref. 1. Did you try to calculate also the adsorption of the CO molecule from the oxygen side?

1 S. Bordiga, C. Lamberti, F. Bonino, A. Travert and F. Thibault-Starzyk, Chem. Soc. Rev., 2015, 44, 7262–7341.

Wilke Dononelli replied: We know from studies of CO adsorption on rutile(110) that CO adsorbed with the oxygen atom at a 5-fold coordinated Ti atom can be a local minimum in the potential energy landscape.¹ We did not consider this configuration in our underlying study of CO adsorption energies on coinage metal nanoparticles. It might be possible to find the configuration with O bound to the coinage metal surface, but it is not reported in the literature. This configuration should be energetically less favourable.

1 H. Spieker and T. Klüner, Phys. Chem. Chem. Phys., 2014, 16, 18743–18748.

Philip Davies commented: In your paper, you mention separate work where you have shown that water catalyses the dissociation of oxygen. Does that statement refer to oxygen dissociation on all 3 of the metals studied here and does this result help explain the well known experimental observation^1,2 that water can enhance the rate of CO oxidation?

1 M. Haruta et al., J. Catal., 2001, 201, 221–224.

2 D. A. H. Cunningham, W. Vogel and M. Haruta, Catal. Lett., 1999, 63, 43–47.

Wilke Dononelli replied: Until now, we have just studied the role of water towards oxygen activation on Au(321) with and without silver impurities and on Au(310).¹ Here, we found that water can be used to activate molecular oxygen in a catalytic cycle with activation barriers of maximum 0.4 eV. Others also found that water might be a possible key in activating oxygen at gold surfaces.² Our calculations indicate that water can enhance the activation of O₂ on gold. On Cu(321) we found an activation barrier of 0.6 eV for the reaction of CO and atomic oxygen. On copper this reaction step seems to be the rate limiting step and not the dissociation of O₂. Here, an associative mechanism might be even more favourable, so the role of water in this context was not determined by our calculations and we cannot give any suggestion about the role of water towards CO oxidation on copper.

1 G. Tomaschun, W. Dononelli, Y. Li, M. Bäumer, T. Klüner and L. V. Moskaleva, J. Catal., 2018, 364, 216–227.

2 F. Xu, I. Fampiou, C. R. O’Connor, S. Karakalos, F. Hiebel, E. Kaxiras, R. J. Madix and C. M. Friend, Phys. Chem. Chem. Phys., 2018, 20, 2196–2204.

Francesca Baletto commented: A couple of years ago, we studied the adsorption of CO on various monometallic clusters.¹ We have reported the effect of the addition of Grimme’s, DFT-D2 and DFT-D3 corrections, and the optPBE vdW-DF on the site preference of CO. Our study shows clearly the importance of studying the adsorption on the various sites, but it was difficult to say what is the best DFT-functional. It would be quite important to do a close comparison between CCSD(T) and those DFT studies to indicate the best ab initio strategy and to quantify clearly the importance of various adsorption sites.

1 J. B. A. Davis et al., J. Phys. Chem. A, 2015, 119, 9703–9709.

Wilke Dononelli answered: In our study we focused on the atop adsorption at the 6-fold low coordinated metal atoms. In your work, you focused on several adsorption positions.¹ The results you present in your work¹ are very interesting. You show that not just quantitative values like adsorption energies but even qualitatively the preferred adsorption site changes depending on the used dispersion correction scheme. As you stated before, it would be highly interesting to make a comparison between your DFT energies and CCSD(T), which we have not focused on, yet. This should be part of a follow-up study.

If other adsorption sites are investigated using the QM/QM′ embedding used in our paper, more atoms have to be described in the high level region, as the coordination number of the other surface atoms are higher. This will of course result in higher calculation time. Additionally we will have to check carefully how many atoms have to be considered in the high level region in order to converge the energies of the embedded system to CCSD(T) results.

1 J. B. A. Davis et al., J. Phys. Chem. A, 2015, 119, 9703–9709.

Francesca Baletto remarked: What are the differences in the adsorption energy for various adsorption sites?

Wilke Dononelli replied: On a M₅₅ nanoparticle two different atop positions and additional positions between two and three metal atoms may be possible adsorption sites. We just focused on the atop adsorption on the 6-fold low coordinated atoms.

Richard Catlow asked: Can you discuss how far we can use these coupled cluster calculations as a benchmark for DFT and other calculations on sorption and reactivity in catalytic science?

Wilke Dononelli replied: First let me start by pointing out that experiments should always be a good choice to be the benchmark for our calculations. When DFT is used, depending whether an oxide surface or a metal catalyst is investigated, either a hybrid functional or a pure GGA functional usually gives the best results with respect to experiment. A problem of DFT in this context is the lack of systematic improvement. If a functional is more expensive from a computational point of view it is not always a “better” functional. For example, Janesko et al. showed that some hybrid functionals tend to incorporate unphysical features when used for describing reactions at metal catalysts.¹ In contrast, high level ab initio methods have the advantage that there is a systematic hierarchy in accuracy. Starting from MP2, over CCSD to CCSD(T) the results should improve. For this reason coupled cluster could serve as a benchmark, especially in the case of relative energies like sorption or activation energies. One of the challenges of coupled cluster theory is the high computational cost.

Nevertheless, being “the gold standard” of quantum chemistry, CCSD(T) calculations definitely can be used as a benchmark for more approximate DFT.

1 B. G. Janesko, T. M. Henderson and G. E. Scuseria, Screened hybrid density functionals for solid-state chemistry and physics, Phys. Chem. Chem. Phys., 2009, 11, 443–454.

Aram Bugaev remarked: When you apply CCSD you face some computational limitations and have to, for example, reduce the basis set, which can have an even more negative effect on your energies than the application of the less precise DFT approach. Have you tried to add not one, but two or maybe even six molecules in order to preserve the initial symmetry of the nanoparticles? How much can it reduce the computational cost?

Wilke Dononelli answered: In the case of the M₁₃ nanoparticles (M = Au, Ag or Cu), we performed CCSD(T) calculations for the entire system. By adding other CO molecules the number of basis functions will increase rapidly, which should be seen as an increase of the calculation time. Of course, the use of symmetry can reduce the calculation time, but we did not consider comparing the calculation time of a symmetric model turning on symmetry operations of the program packages and turning them off again. However, when the system size has to be increased in order to achieve symmetry, the time loss due to enlargement should be greater than the time gain due to symmetry.

In case of M₅₅ nanoparticles (M = Au, Ag or Cu) we used a QM/QM′ embedding scheme in order to get CCSD(T)/PBE results for CO adsorption energies. In this embedding scheme, the adsorbate, the adsorption centre, and the nearest neighbours were considered in the high level region. If another CO molecule was adsorbed on the other side of the nanoparticle within a perfectly symmetric structure, the number of atoms of the high level region had to be twice as high. With a formal scaling of N⁷ the calculation time would be 2⁷ times as long. Even if one naively suggests that an inversion centre or mirror axis would enhance the speed of the calculation by a factor of two, then no speed-up can be achieved by doubling the size of the subsystem. For these systems with 55 metal atoms in our embedding scheme, adding another adsorbate would always result in much higher calculation times.

David Willock commented: In your CCSD(T) calculations the effect of dispersion interactions between the adsorbate and the metal cluster will be taken into account. In the DFT calculations it is now common practice to introduce dispersion using an additional parameterised term (e.g. D2 or D3 corrections). Can you make a comparison of your DFT and CCSD(T) results to make some comment on the accuracy and appropriateness of such dispersion corrections? This seems particularly important for metal nanoparticles for which an atom-by-atom parameterised calculation of dispersion seems difficult to justify.

Wilke Dononelli answered: Thank you for this question. There is an ongoing debate about using semi-empirical corrections like D2 or D3, where in most cases people use additional DFT functionals in order to give a statement about whether dispersion corrections have to be used or not. We calculated the adsorption energy of CO on Au₁₃ using D3. Without dispersion correction we found E_ads = −1.21 eV at the PBE level of theory. Using D3 we found −1.30 eV. At the CCSD(T) level of theory we found −0.88 eV as shown in our paper. On Au₅₅ we found E_ads = −0.84 eV for PBE, E_ads = −0.86 eV for CCSD(T)/PBE and E_ads = −1.00 for PBE-D3. If we claim that CCSD(T) gives good results, then PBE-D3 seems to overestimate the binding strength of CO. An indication for difficulties in the parametrisation you are mentioning may be seen in the differences of the total energies. The total energy of CO in the gas phase using pure PBE or PBE-D3 is virtually identical, whereas the absolute value of the total energy of Au₅₅ is 17.53 eV (∼11%) higher for PBE-D3 compared with pure PBE. This should not be a statement that dispersion corrections are incorrect in a general sense. For example, others showed that in the case of alcohols or alcoxy species, van der Waals interactions might be favourably described by using semi-empirical dispersion corrections.¹

1 Y. Xu, W. Chen, E. Kaxiras, C. M. Friend and R. J. Madix, J. Phys. Chem. B, 2018, 122, 555–560.

Alexander Genest remarked: The results of your plane-wave based methods and your atom-centered basis set methods matched perfectly, even though they should not due to the basis set superposition error (BSSE), which affects only the latter approach. Did you correct the results of your atom-centered results for this error? Do you have an estimate for the magnitude of the error?

Wilke Dononelli responded: In our atom-centred based calculations we did not consider the BSSE. As a first step we tried to find a basis set where we find adsorption energies at the PBE level of theory that were in good agreement with calculations using the PAW method. The adsorption energies calculated using the two different methods only show the same results by chance. Using the counterpoise correction (CP correction) by Boys and Bernadi, we find an adsorption energy of −1.01 eV at the PBE level of theory for CO on Au₁₃ (−1.21 eV without counterpoise correction). Here the BSSE is small. On the other hand the basis set incompleteness error (BSIE) has to be considered as well. Using a bigger def2-qzvp basis set at the PBE level of theory we find −1.25 eV without CP correction and −1.29 eV using CP correction. Unfortunately, we were just able to calculate adsorption energies for bigger def2-qzvp basis sets at the DLPNO-CCSD(T) level of theory. Here we found CO adsorption energies of −1.02 eV on Au₁₃, −0.49 eV on Ag₁₃ and −1.15 eV on Cu₁₃. But at the moment, we are not able to claim whether this deviation results from the BSSE or is a fact of the DLPNO Ansatz used for the bigger basis sets. Nevertheless, all results show that PBE tends to overestimate the binding energy of CO on M₁₃ nanoparticles (M = Au, Ag or Cu).

Rutger van Santen commented: The result you report on the independence of the CO bond strength on nanoparticle size for the Au₁₃ and Au₅₅ particles using CCSD(T) calculations are remarkable. In my experience the 13 atom nanoparticle has exceptionally strong bonds along the surface because of low coordination numbers, compared with that of the 55 particle. For this reason, for atop chemisorption the interaction strengths do not necessarily favour bonding to the 13 atom particle. If this reasoning is correct, the main importance of the CCSD(T) calculations is to describe metal bonding in the cluster correctly. One would be able to deduce this by comparing the bandwidths, or differences between HOMO and LUMO levels in the two systems.

An important point is that the attractive part of the chemical bonding within the metal particles is dominated by the s-p valence electrons and the repulsive part by the doubly occupied d-orbitals. It is critical how the contributions of the (occupied) d-valence electrons are accounted for. They give a repulsive contribution to chemical bonding, that most likely is sensitive to correlation. If there is a change in the extension of the d-atomic orbitals, this will affect also the (repulsive) part of the interaction with CO. Did you study different ways to include the d-valence electrons in your study? Also have relativistic effects been accounted for? The vibrational frequency of the CO stretch mode should show an upwards shift compared to the gas phase (strengthening of CO sigma bonds). It most likely correlates with ad-molecule bond strength in your case. Do you actually find this?

Wilke Dononelli responded: Thank you for your comment. The first point you raised is that you expect differences in the bandwidth for the two cluster sizes we investigated. Your suggestion is completely correct. In order to estimate the bandwidth in an accurate way, we performed time dependent DFT calculations using the PBE functional (TD-PBE) and evaluated the excitation energies of the nanoparticles. To verify the TD-PBE results we performed equation of motion CCSD (EOM-CCSD) calculations for the smaller M₁₃ nanoparticles. All excitation energies are summarized in Table 1.

Table 1 Excitation energies of M₁₃ and M₅₅ nanoparticles (M = Au, Ag or Cu). Energies in eV

Nanoparticle	E* in eV
Au13 (EOM-CCSD)	0.29
Au13 (TD-PBE)	0.20
Au55 (TD-PBE)	0.07
Ag13 (EOM-CCSD)	0.20
Ag13 (TD-PBE)	0.15
Ag55 (TD-PBE)	0.07
Cu13 (EOM-CCSD)	0.12
Cu13 (TD-PBE)	0.13
Cu55 (TD-PBE)	0.07

For Au and Ag the EOM-CCSD excitation energies are a little bit higher than TD-PBE but still in good agreement. For the smaller M₁₃ nanoparticles the excitation energy ranges from 0.13 eV to 0.20 eV at TD-PBE level of theory whereas the energy is 0.07 eV for the bigger Au₅₅, Ag₅₅ and Cu₅₅ clusters, respectively. These results indicate that the bigger M₅₅ nanoparticles show a metal-like character, whereas the smaller nanoparticles exhibit a small band gap.

In the second part of your question you are giving good insight into the understanding of the influence of d-orbitals or d-bands to chemical bonding in metal–adsorbate interactions. In our study, we compared VASP calculations using the PAW approach to Gaussian calculations where we used a Hay and Wadt basis set with pseudo potentials. The basis set was chosen in order to give good results with respect to the adsorption energies calculated with VASP at the PBE level of theory. Within these two different approaches, two different ways of including the d-valence electrons were used. Both methods result in similar adsorption energies of CO on the same nanoparticle (using PBE). Relativistic effects have been taken into account by using pseudo potentials. We think that your comment about the repulsive character of doubly occupied d-orbitals of the metal nanoparticles could be correct. In order to investigate the role of electron correlation, a possible way could be to freeze the electrons of the MOs constructed from the d-AOs during the correlation part of the calculation. The resulting adsorption strength should be stronger. The difficulty will be to decide which molecular orbitals have to be considered. Your argument is based on a one-particle picture, which might be difficult to transfer into a many particle perspective. In addition, binding should be a local phenomenon, whereas orbitals in metal nanoparticles exhibit a delocalised character. Despite these difficulties, we will try to examine your comment in future work.

To answer your last question, we summarised the wavenumbers of CO stretching vibrations in the gas phase and on the different Au nanoparticles calculated at the PBE level of theory in Table 2. Comparable to metal carbonyl species a red shift of the CO stretching frequency (compared to CO in the gas phase) was observed for all systems in our investigations. However, due to the different nature of the electronic structure of different sized nanoparticles no simple correlation between the bond strength and the shift in frequency can be made.

Table 2 CO stretching frequencies on different Au surfaces and nanoparticles calculated at the PBE level of theory. Wavenumbers are given in cm⁻¹

Surface/nanoparticle	Wavenumber in cm^−1
CO-Au(321)	2056.9
CO-Au55	2058.5
CO-Au13	2076.7
COgas	2122.8

Valerii Bukhtiyarov asked: When you study oxygen adsorption on Group 11 metals, you try to compare the activation energy of adsorption on different clusters. Are there any differences in oxygen adsorption on Ag₁₃ and Ag₅₅ clusters? I mean the activation energy of adsorption or subsequent reactivity of adsorbed oxygen species.

Wilke Dononelli answered: In the underlying study, we compared the dissociation energy of O₂ on M₅₅ (M = Au, Ag or Cu) nanoparticles to dissociation barriers on periodic M(321) surfaces. Additionally, we focused on the adsorption energy of CO on M₁₃, M₅₅ and M(321). We did not focus on reaction barriers of oxygen and other molecules on M₁₃ nanoparticles or adsorption energies of oxygen on M₁₃ or M₅₅ nanoparticles. Nevertheless, we can use the geometries of the initial state of the O₂ dissociation reaction on the M₅₅ nanoparticles and the M(321) surfaces to calculate the adsorption energies of these geometries as shown in our paper. Here, E_ads is calculated as:

E _ads = E(O₂,adsorbed) − E(O₂,gas phase) − E(surface)

The adsorption energies are listed in Table 3.

Table 3 PBE adsorption energies for O₂ on M₅₅ and M(321) surfaces (M = Au, Ag or Cu). Energies are given in eV

Surface/nanoparticle	E ads in eV
Au(321)	−0.15
Au55	−0.25
Ag(321)	−0.43
Ag55	−0.42
Cu(321)	−1.24
Cu55	−1.56

Bruce Gates opened discussion of the paper by Arunabhiram Chutia: What happens when the isolated positively charged/negatively charged gold species on the support surface are probed with CO?

Arunabhiram Chutia replied: Based on whether the CO probe molecule is adsorbed on the isolated positively or negatively charged Au species it may display different electronic properties, i.e., if CO is adsorbed on a negatively charged Au species then the antibonding π* orbitals of CO will be populated weakening the C=O bonding. On the other hand, if CO is adsorbed on a positively charged species it may give an opposite effect. Therefore, a calculation of vibrational frequency for CO adsorbed on negatively charged Au may give a lower frequency compared with CO adsorbed on a positively charged Au species. This is certainly something we will be considering as we extend our work further.

Carlo Lamberti commented: I found very interesting the in-depth knowledge of the electronic density of states (DOS), both occupied and non-occupied, of the gold on ceria system that you report in your study, and how it can be significantly influenced by the adsorption site. It would be very interesting to perform some experiments (XPS, XANES, XES) to confirm your calculations. In this regard, are you already in contact with some experimental group?

Arunabhiram Chutia replied: Thanks for your comment. In one of our previous studies on the interaction of Cu and CuO clusters with CeO₂(110) surface, we reported our theoretical studies in conjunction with XAFS experiments¹ but for this study we have not yet contacted any experimental groups. However, we would be open to such a collaboration so that a direct comparison between experiments and our calculations on Au/Au₂ clusters on CeO₂ surfaces could be done.

1 A. Chutia, E. K. Gibson, M. R. Farrow, P. P. Wells, D. O. Scanlon, N. Dimitratos, D. J. Willock and C. R. A. Catlow, Phys. Chem. Chem. Phys., 2017, 19, 27191–27203.

Parasuraman Selvam asked: Have you considered the surface coordination for the interaction? If so, what is the influence on the reduction of cerium ions as compared to other low-index planes, viz., (100) and (111)?

In the case of defect formation, two electrons from the oxide ions are transferred to two cerium ions neighbouring the vacancy site, so that the cerium ions are reduced from tetravalent to trivalent state. After formation of the vacancies, reactive sites that are present can interact with gold atom/atoms allowing reduction or a partial negative charge. Alternatively, it may also cause partial reoxidation of the ceria surface. How do you view the partial positive charge on gold atom/atoms?

Arunabhiram Chutia answered: This study was performed on Au adatoms adsorbed on CeO₂(110) surface only. Previously however, the adsorption of Au on the other low index surfaces were reported and in those cases similar reduction of the Ce ions has been observed. In our study we have seen Au⁺, Au⁻ and Au^δ− species when an Au atom is adsorbed on CeO₂(110) with and without O-defects. Our observations of these species are based on our analyses on partial density of states for Au s orbitals and Bader charges, which clearly showed that the Au adatom on pristine CeO₂ transferred its s electrons to reduce a surface Ce ion on the surface giving rise to Au⁺ species. On the other hand, when Au atoms were adsorbed on an O-vacancy, we show two cases: first, the Au adatom is partially reduced giving rise to a Au^δ− species, and secondly an Au− species when it was fully reduced. However, in our analysis no Au^δ+ species were observed for the Au adatom on CeO₂ with and without O-defects. If there were any Au^δ+ species , then we would have seen a small fraction of the occupied s orbital signature above the Fermi energy.

Wilke Dononelli commented: In our own work we looked at the Mulliken atomic charges of the gold atoms in Au₁₃ and Au₅₅ clusters. Here we found highly positive charges of +5.7 e and +6.1 e at the central atoms of Au₁₃ and Au₅₅, respectively. All surface atoms are negatively charged. What would you suggest? How does this charge distribution change when the nanoparticles are adsorbed on an oxide surface like CeO₂? Could you please additionally comment on what is the explanation for this charge transfer you saw in your examples?

Arunabhiram Chutia answered: In our studies on the interaction of Au adatoms and Au₂ clusters with the CeO₂(110) surface we analysed Bader charges and we saw that, in the Au₂ case, the Au atom close to the CeO₂(110) surface has a positive charge (+0.274 e) and a negative charge (−0.308 e) for the non-interacting Au atom away from the surface leading to a Au^δ+–Au^δ− like species. From this analysis we concluded that the Au atom closer to the surface may share its electron partially with the other Au atom and simultaneously partially reduce a Ce ion on the surface. However, for bigger clusters such as those of Au₁₃ and Au₅₅ interacting with CeO₂ the distribution of charge may display a different behaviour. Previously, Tereschchuk et al. used the DFT+U method and analysed Bader charges to report the interaction of 13 atom clusters of Au, Pd, Ag and Pt with CeO₂ and they concluded that the topmost layer of these pyramidal TM₁₃ clusters interacting with the CeO₂(111) surface had a slightly negative charge (−0.16 e for Au), the middle layer is almost zero (+0.03 e for Au) displaying bulk like behaviour, and the layer in direct contact with the O atoms of the CeO₂(111) surface was positively charged (+2.54 e for Au) indicating charge transfer.¹ We looked into the electron transfer phenomenon in Au (or Au₂) interacting with the CeO₂(110) surface and found that, based on the initial geometry, we can observe interesting Au species such as Au⁺, Au⁻, Au^δ− and Au^δ+–Au^δ− due to electron transfer, which could be clearly understood in terms of the electronic configuration of Au atom ([Xe]4f¹⁴5d¹⁰6s¹). Since Au has one electron in its valence s-orbital, it can either donate or gain an electron, which gives rise to Au⁺ or Au⁻ species, respectively. There is also a possibility that the Au atom is partially reduced due to its interaction with the CeO₂(110) surface with an O-vacancy and in such a case we see an Au^δ− species. We confirmed these observations by our analysis of Bader charges and partial density of states.

It is always difficult to link the calculated charges from any method to a formal oxidation state of the Au atom unless good model compounds have been considered. In an earlier work for Au₁₀ on Fe₂O₃(001) Willock and co-workers used AuCl and AuCl₃ molecular species calculated at a consistent level of theory to estimate the Bader charge of Au(I) and Au(III), respectively. This gave Bader charges of 0.33 e for Au(1+) and 0.81 e for Au(3+) so that the calculated charges always seem to underestimate the formal oxidation state, presumably as the bond is partially covalent.² With this in mind, the charges you quote seem rather large and would imply a formal oxidation state for the Au atoms at the centre of your cluster which would be outside of the normal range for Au.

1 P. Tereschchuk, R. L. H. Freire, C. G. Ungureanu, Y. Seminovski, A. Kiejna and J. L. F. Da Silva, Phys. Chem. Chem. Phys., 2015, 17, 13520–13530.

2 Kara L. Howard and David J. Willock, Faraday Discuss. 2011, 152, 135–151.

Aram Bugaev asked: You have obtained a very interesting result that the binding energy of gold on ceria surface is higher on top of the oxygen vacancy. However, it would be even more interesting to analyse this result from the point of catalytic applications. For example, it has been shown that addition of platinum nanoparticles affects the reducibility of the ceria.¹ In our recent work, we have also shown that there is a critical temperature above which the metal/ceria interface does not play a dominant role in the catalytic oxidation of carbon monoxide.²

From the results, that you have shown for gold it can be concluded that reduction of ceria should be favorable at the gold/ceria interface. In addition, it would be interesting to see also the energy barriers of removing an oxygen atom and correlate them with the experimental data.

1 G. N Vayssilov, Y. Lykhach, A. Migani, T. Staudt, G. P Petrova, N. Tsud, T. Skála, A. Bruix, F. Illas, K. C. Prince, V. Matolín, K. M. Neyman and J. Libuda, Nat. Mater., 2011, 10, 310.

2 A. A. Guda, A. L. Bugaev, R. Kopelent, L. Braglia, A. V. Soldatov, M. Nachtegaal, O. V. Safonova and G. Smolentsev, J. Synchrotron Radiat., 2018, 25, 989–997.

Arunabhiram Chutia responded: Thanks for your comments. We can see from our calculations that the adsorption energy of Au on top of an O-vacancy is stronger (−1.992 eV) compared with the most stable site on CeO₂(110) surface without O-defects (−1.132 eV), which, as you mentioned, also means that the reduction of the CeO₂ surface can be significantly enhanced by Au atom adsorption. So far as the energy barriers for removing an O-atom is concerned, Hernandez et al. previously proposed a detailed mechanism on this.¹ They found formation energies of O-vacancies are between 0.53–1.07 eV (depending on the initial adsorption site) and since these values are lower than the formation energy for the vacancy on the clean surface (2.56 eV) it can be expected that the presence of Au will result in a higher concentration of vacancies. Willock and co-workers also looked at this effect for Au₁₀ on Fe₂O₃(0001) and found lower O defect formation energies with the Au cluster present than for the clean surface. They analyzed charges to show that for defects formed near the metal/support interface, the electrons that are left behind can be partially accommodated on the metal cluster and it is likely that this will be a general observation when metal nanoparticles are deposited on reducible oxides.²

1 N. C. Hernández, R. Grau-Crespo, N. H. de Leeuw and J. F. Sanz, Phys. Chem. Chem. Phys., 2009, 11, 5246–5252.

2 S. W. Hoh, L. Thomas, G. Jones and D. J. Willock, Res. Chem. Intermed., 2015, 41(12), 9587–9601.

Cynthia Friend opened a general discussion of the papers by Nora de Leeuw, Wilke Dononelli and Arunabhiram Chutia: Should there be an “industry standard” by which we benchmark calculations? The nice study on gold clusters comparing various DFT methods to the CCSD(T) is an example. However, as there are so many different “flavors” of DFT and also different details in the mode, can we reliably compare results across different groups? What is a good compromise, given that the methodology is rapidly evolving?

Nora de Leeuw answered: The use of computational tools to describe coalescence phenomena is relatively new and a standard is missing. Furthermore, the gap between experimental particle size and the sizes that can be reached computationally makes benchmarking still rather difficult. The transferability of results is the same as for the comparison of any other computational studies. However, a concerted effort by groups working in the field to identify general trends and define a "standard" method would be very welcome.

Arunabhiram Chutia replied: As you mentioned, methodologies are rapidly evolving and there are attempts to attain a very high level of chemical accuracy, which of course comes at a considerable computational cost. Therefore, a good compromise is to compare our results obtained from our calculations at different levels of theory with reliable experimental data and then choose the appropriate theory for our studies. Perhaps in this direction there is also a very strong need for both theoreticians and experimentalists to work together to design benchmarking experiments. It may also be very helpful if experiments can guide us in modelling. For example, we recently reported studies in which Inelastic Neutron Spectroscopy (INS) and XAFS experiments guided us in predicting reliable models and performing DFT based theoretical calculations to study the geometrical and electronic properties of catalytically interesting species on metal, metal oxide surfaces and in zeolites.^1–3

1 A. Chutia, E. K. Gibson, M. R. Farrow, P. P. Wells, D. O. Scanlon, N. Dimitratos, D. J. Willock and C. R. A. Catlow, Phys. Chem. Chem. Phys., 2017, 19, 27191–27203.

2 A. Chutia, I. P Silverwood, M. R. Farrow, D. O. Scanlon, P. P. Wells, M. Bowker, S. F. Parker and C. R. A. Catlow, Surf. Sci., 2016, 653, 45–54.

3 A. J. O’Malley, S. F. Parker, A. Chutia, M. R. Farrow, I. P. Silverwood, V. García-Sakai and C. R. A. Catlow, Chem. Commun., 2016, 52, 2897–2900.

Wilke Dononelli responded: Let me start again by saying that the standard as a benchmark for theoretical calculations should be experiments. But sometimes experiments are not available or cannot be defined for “perfect” systems, like the ones that are used for most theoretical investigations. It will not be easy to find a standard or benchmark in the landscape of DFT. Some functionals are very good in describing reactions on oxide surfaces, but these might fail in describing reactions at metal nanoparticles. Coupled cluster on the other hand might be a good choice as a benchmark method. This does not mean that coupled cluster will always give results that are more accurate with respect to experiment than DFT. But there is an intrinsic accuracy within coupled cluster theory, like I have stated before. By considering higher cluster amplitudes, the result should always become more accurate, because it will converge towards the limit of the full configuration interaction, which gives the exact energy of a system (in the basis set limit). A problem of coupled cluster theory is the high computational cost. For example, CCSD(T) scales formally as N⁷, where N is the system size. On the other hand, standard LDA or GGA (DFT) implementations scale like N³-N⁴. The high computational cost of coupled cluster theory makes it unavailable at the moment for most standard applications, where more than 20 heavy transition metal atoms and reactants or adsorbates have to be described, but workarounds are available. These could be embedding techniques as shown in our work, expansion techniques like the method of increments as used by Paulus¹ and Staemmler² for oxide surfaces or Voloshina³ for metal surfaces, and local variations like LCCSD(T) by Werner⁴ or DLPNO-CCSD(T) by Neese.⁵ But these workarounds still have to be checked carefully for each individual system of interest.

1 C. Müller, B. Herschend, K. Hermansson and B. Paulus, J. Chem. Phys., 2008 128, 214701.

2 V. Staemmler, J. Phys. Chem. A, 2011, 115(25), 7153–7160.

3 E. Voloshina, Phys. Rev. B, 2012, 85, 045444.

4 H.-J. Werner and M. Schütz, J. Chem. Phys., 2011 135, 144116.

5 M. Saitow, U. Becker, C. Riplinger, E. F. Valeev and F. Neese, J. Chem. Phys., 2017, 146, 164105.

Francesca Baletto asked: In the near future of first-principle calculations applied to investigate metallic nanoparticle properties, what is your suggestion about the use of DFT+U in comparison with other schemes, such as CCSD(T), TD-DFT and GW?

Arunabhiram Chutia responded: The DFT+U methodology takes into account strong on-site Coulomb interaction of localised electrons and has been widely used to explore the electronic structure of materials with f-electrons and metal oxides. A recent study by Beridze et al. suggested that use of DFT+U, with the Hubbard U parameter derived by ab initio methods could be a good choice over CCSD(T); the reason being that even though the CCSD(T) method gives a very accurate electron correlation energy, this method is however computationally very expensive and can be used to study small clusters.¹ However, as we heard from our previous speaker (Dononelli et al.), with a combination of CCSD(T) and DFT methods employing the QM/QM scheme, we can perform more accurate calculations and in the near future this could be perhaps routinely done. On the other hand, if we are interested in studying the excited states of these materials then we can employ the TDDFT+U method. In this regard it is also worth mentioning that recently, Tancogne-Dejean et al. reported the implementation of self-consistent DFT+U and TDDFT+U methods.² So far as the GW method is concerned, similar to other Green’s function approaches, it replaces the unknown XC-potential by a self-energy and in recent years it has been able to calculate transport properties of single molecules, and correctly predicted the band gap of a large number of semiconductors. In terms of the computational cost, it is also expensive, however, in due course we may be able to use this method more routinely. Finally, the choice of methods greatly depend on which properties of materials we are looking at.

1 G. Beridze et al., J. Phys. Chem. A, 2014, 118, 11797–11810.

2 N. Tancogne-Dejean et al., Phys Rev. B, 2017, 96, 245133.

Nora de Leeuw responded: DFT+U and GW are usually used to improve the calculated band gap of semiconductors and insulators, in comparison with the experimental results. CCSD(T) and TD-DFT are employed to study excited states. Neither of those properties are in the scope of our paper as we are interested in the electronic and geometric structures, and the mobility of the Ni clusters on the oxide surfaces. Additionally, the next step is to study the reverse water gas shift reaction at the interface between the cluster and the oxide surfaces: this does not require any of the latter techniques that are more time consuming than DFT.

Cynthia Friend commented: In order to evaluate whether water produced in situ can promote O₂ dissociation, the short lifetimes of both O₂ and water need to be considered. Dioxygen specifically is very weakly bound and will have a short surface lifetime under reaction conditions. Kinetic modeling that accounts for such a short lifetime under specific conditions needs to be considered.

Wilke Dononelli replied: You are completely right. Especially the adsorption strength of O₂ is very weak. By co-adsorption of water and O₂ the adsorption strength can be increased. But like you have stated, kinetic factors are important to understand this reaction.

Laura Torrente-Murciano continued the discussion of the paper by Arunabhiram Chutia: Your paper nicely discusses how the presence of Au atoms/NP on the surface of ceria modifies the oxidation state of the cerium. My question is related to the potential effect of gold atoms/clusters/NPs when they are inside the crystal structure of the ceria (e.g. a few atomic layers below the surface) and your opinion about its effect on the ceria surface (e.g. increased amount of oxygen vacancies due to higher concentration of Ce³⁺)?

Arunabhiram Chutia responded: Even though we have investigated the effect of creating O-vacancies in the bulk of CeO₂ on the reduction of Ce(IV) ions to Ce(III) ions and its influence on the adsorption of Au adatoms, we have not yet investigated the potential effects of Au atoms/clusters in the bulk of CeO₂. But there is an interesting work by Kehoe et al. where they investigated the interaction of a range of divalent dopants with CeO₂ using the DFT+U method.¹ They found that these dopants adopted the coordination of their own oxide and they also reported that different coordination environments can create weakly or under-coordinated oxygen ions, which could be more easily removed than in pure CeO₂. Therefore, if we have Au atoms in the bulk we may see Au(III) oxide like structures.

1 A. B. Kehoe et al., Chem. Mater., 2011, 23, 4464–4468.

Conflicts of interest

There are no conflicts to declare.

---