HAL is a

Working closely with experimentalists on the comprehension of the surface properties of catalytically active organometallic nanoparticles (NPs) requires to develop several computational strategies which significantly differ from the cluster domain where a precise knowledge of their optimal geometry is a mandatory prerequisite to computational modeling. Theoretical simulations can address several properties of organometallic nanoparticles: the morphology of the metal core, the surface composition under realistic thermodynamic conditions, the relationship between adsorption energies and predictive descriptors for reactivity. It is in such context that an integrated package has been developed or adapted in our group: (i) one tool aims at building a wide variety of the typical shapes exhibited by nanoparticles. Using Reverse Monte Carlo modeling, a given shape can be optimized in order to fit pair distribution function data obtained from X-Ray diffraction measurements; (ii) trends in density functional theory (DFT) adsorption energies of surface species can be rationalized and predicted by making use of simple descriptors. This is why we have proposed an extension of the d-band center model, that leads to the formulation of a generalized ligand-field theory. A comparison between cobalt and ruthenium is proposed in the case of a 55-atoms nanocluster. The accuracy of the generalized coordination number [Angew. Chem.

Int. Ed. (2014) 53, 8316], a very simple coordination-activity criterion, is also assessed; (iii) the builder package is completed by the steric-driven grafting of ligands on the surface of metal NPs. It easily generates structures with adjustable surface composition values and coordination modes; (iv) after a local optimization at the DFT level of theory, DFT energies and normal modes of vibration can feed a general tool based on the ab initio thermodynamics method. This method aims at easily calculating an optimal surface composition under realistic temperature and pressure conditions. On top of that, we also show to which extent the knowledge of the density of states (DOS) and of the crystal overlap Hamilton population (COHP), both projected from a plane-wave basis set to a local basis set, sheds light on the metal core -ligand chemical bonding.

NPs. It usually consists in the reduction of a metal salt precursor in solution, the presence of a stabilizing agent preventing aggregation of the formed NPs. Nanocatalysis, which combines the design of heterogeneous catalysts (monodispersity, shape control, presence of ligands) with the monitoring of their surface activity using new microscopic and spectroscopic methods, has recently emerged as a major new field for the rational design of improved catalysts. Provided that understanding at the atomic level reactive processes that occur on their surface allows the fine-tuning of nanocatalysts, first-principle calculations can guide their conception, both in terms of activity and selectivity [START_REF] Gerber | Nanocatalysis with Transition Metal Particles: Where are we Now[END_REF][START_REF] Carchini | [END_REF] , although it remains challenging. The reason lies in the difficulty to accurately describe such complex systems, which exhibit a nanoscale metallic core with ligands on its surface, often partially mobile, and not so easy to characterize experimentally. Yet, prior to the calculation of data with a fairly accurate computational method, a valuable theoretical rationalization requires to design a relevant chemical model. Let us focus on ruthenium NPs, thoroughly investigated in our group these last years. As good catalytic agents, they can decompose ancillary ligands which are introduced in the chemical medium to prevent coalescence. 4 Several features are now firmly established by various experimental and theoretical approaches and can be used as a guide to define a theoretical model for dressed metal NPs: (i) the RuNPs we are concerned with have a very small metal core (typically between 1 nm 5 and 4 nm [START_REF] Philippot | Organometallic Ruthenium Nanoparticles and Catalysis[END_REF][START_REF] Amiens | [END_REF] ), which characterization by X-ray diffraction shows a crystalline character with the expected hcp structure of bulk ruthenium; 8 (ii) they can accommodate at least 1.3 hydride per surface ruthenium atom, a value probably slightly underestimated owing to the method employed; 9 (iii) a combination of [START_REF]From theory to application[END_REF] H and 2 H NMR has shown that hydrides are very mobile at the surface of RuNPs and that hydride-deuteride exchange between adsorbed hydrides and an isotopic D 2 gas source occurs under mild conditions; 10 (iv) a joint DFT/NMR study has proven the coexistence of ontop, edge-bridging, and face-capping deuterides at NPs surfaces, in contrast to what is observed and calculated on bare Ru(0001) surfaces; 11 (v) in the context of the Fischer-Tropsch reaction, ab initio thermodynamics and gas-phase chromatography measurements have shown that surface hydrides are replaced by CO ligands under a syngas flow. 12 This remarkable and broad knowledge of the metal core, of the surface composition and of the role of protecting ligands both on the size, on the shape, and on the catalytic properties has been collected during several years of experimental and theoretical studies performed on ruthenium NPs. On the theoretical side, there is a need to speed up the rationalization and optimization processes by integrating several tools and models into a rational framework. This involves to have a reliable structural description of the metal core and to describe the species coordinated on the metal surface both in terms of steric and electronic effects. Owing to the electronic states degeneracy at the Fermi level, periodic-DFT methods are considered as a tool of choice to tackle first principles calculations on large metal clusters. [START_REF] Gerber | Nanocatalysis with Transition Metal Particles: Where are we Now[END_REF] Inorganic NPs can exhibit a wide variety of morphologies which are rarely available in molecular builders. A versatile tool aiming at providing the coordinates of several crystalline NPs, polyhedra, is presented in section 2. It also includes a parallelized Reverse Monte Carlo module in order to optimize a geometry with respect to X-ray diffraction data. 13,[START_REF] Neder | Diffuse Scattering and Defect Structure Simulations. A cook book using the program DIS-CUS[END_REF] The electronic structure of metal NPs is analyzed in section 3 in terms of a variant of the d-band center model of Hammer and Nørskov. 15 This conceptual DFT approach, in line with the Sabatier principle, 16 the Br{\o}nsted-Evans-Polanyi relationship, 17,18 and Balandin's volcano curves, [START_REF] Balandin | Advances in Catalysis[END_REF] provides theoretical descriptors for adsorption strength. Such monoelectronic descriptors could be useful guides to design efficient nanocatalysts with sites having a specific activity. [START_REF] Del Rosal | [END_REF] They are depicted as color maps, which give a straightforward point of view of possible reactivity spots. A comparison between cobalt and ruthenium nanoclusters will be shown. Based on the standard d-band center model, we recently proposed the definition of effective d atomic orbitals for each surface atom, which energy depends on the field generated by the other metal atoms and by the surface ligands. [START_REF] Del Rosal | [END_REF] The standard d-band center model and its by-product will then be compared to the generalized coordination number, which is a simple alternative to electronic structure analysis for predicting adsorption properties on pure metal NPs. [21][22][23] We also show both in sections 3 and 4 that the Crystal Orbital Hamilton Population (COHP) is a complementary analysis method of the electronic structure that provides insights on the chemical bonding. [24][25][26] Colloidal NPs are not bare metal crystals, they accommodate several chemical species at their surface, and possibly underneath. Surface species may not be innocent with respect to a catalytic process, they will modulate the electronic properties of the catalyst and will alter the generalized coordination number of surface metal atoms. To understand the catalytic properties of NPs from a mechanistic point of view, it is therefore crucial to consider decorated metal surfaces, with realistic coverage values. In order to easily generate metal NPs dressed by atoms or ligands, we have developed a tool (dressNPs) which accounts for the steric hindrance at the surface by a CPK model (section 4). The steric clash is minimized using the Monte-Carlo simulated annealing method. dressNPs facilitates the automatic generation of ligand-covered NPs, prior to DFT optimization. Following the ab initio thermodynamics methodology [START_REF] Scheffler | Handbook of Surface Science[END_REF] implemented in our aithermo utility, the resulting energies as well as the normal modes of vibrations can then be used to determine stability phase diagrams. Overall, this integrated package mainly based on firstprinciples calculations can contribute to design realistic models of nanocatalysts. Owing to the prohibitive cost of hybrid functionals applied to such very large compounds, [START_REF] Paier | [END_REF] periodic DFT studies were systematically achieved with the widely used GGA-PBE functional. However, nothing prevents to take into consideration in our packages results obtained within the hybrid functional class of approximation to DFT.

Morphology of the metal core of colloidal NPs

Introduction

The chemical and physical properties of compounds are determined by their atomic structure, so the determination of experimental or theoretical data related to the structure are essential in chemical science. Controlling the morphology of nanoparticles is of key importance for most applications. For example, in magnetic devices, a good control of the morphology is very important in order to have well-defined magnetization axes to store or to process information 29 . Different shapes can be obtained by varying the synthesis conditions like the type of metallic precursors and stabilizers and the surrounding atmosphere. [30][31][32] NPs can be classified as single-crystalline polyhedral (cubes, octahedra, tetrahedra, truncated octahedra, ...) or nonpolyhedral (plates, rods, ...) shapes, or as multiply twinned shapes (decahedra, Marks decahedra, icosahedra, ... See Nepijko et al. 33 ).

Growth and nucleation reactions at the basis of colloidal synthesis are governed by intricated thermodynamic and kinetic parameters. 34 Among these, surface energies and surface area-to-volume ratio are one of the key parameters in understanding the morphology of colloidal NPs. However, little is known at the atomistic level on the growth mechanism nor on the relationship between the final structure and the experimental conditions. On the quantum chemistry side, the potential energy surface (PES) is, within the Born-Oppenheimer approximation, the solution of the Schrödinger equation in the framework of frozen nuclei. Yet, in the case of large clusters, the dimensionality of the configuration space required to describe them becomes very large. [START_REF] Wales | Energy Landscapes: Applications to Clusters, Biomolecules and Glasses[END_REF] Finding in silico the optimal structure of a nanoparticle in its environment is accordingly out of reach of current computational algorithms and facilities, even when the energy is calculated with empirical potentials. Although Wulff-von Laue constructions [START_REF] Wulff | [END_REF][37][38] based on first-principles surface energies 39,40 of bare or passivated surfaces can easily be sketched, the kinetic parameters that underlie the growth process are not understood. On top of that, its application to alloys 41 or to twinned NPs 42 is not straightforward and it is also unable to bring details on surface defects, which are expected to play an important role in catalytic processes. To some extent, computer simulations can provide additional information concerning the different parameters having an influence on the morphology of formed NPs but also concerning the intrinsic link between shape and catalytic, magnetic or optical properties of NPs. However, this high wealth in terms of shapes, the weak experimental knowledge regarding the surface composition and the diversity of chemical events that may occur on the surface are theoretically challenging. In order to overcome a first issue, i.e., the building of realistic metal NPs, we have developed a software, polyhedra, that generates a wide range of nanostructures.

Generation of atomic structures by a molecular builder

Let us first consider the case where size, shape and crystal packing are experimentally well-defined. It is possible to design a 3D model that accounts for experimental NPs. They can be regarded as arising from two different approaches. In the "bottomup" approach, NPs are considered as aggregates of individual atoms leading to well defined clusters. The "top-down" approach is characterized by the size reduction of bulk material towards the obtainment of nanometric material. Both exist experimentally, [43][44][45] even though "top-down" synthesis often leads to less control over the size and/or morphology of NPs and are therefore less common in chemistry. [46][47][48] The polyhedra code provides different NPs shapes by either a "top-down" or a "bottom-up" approach. The first kind of structure that polyhedra can design is model clusters. These highly symmetrical particles are built via a "bottom-up" approach. The final particles, with the wanted size, are obtained by adding several layers to the smaller cluster of the desired shape. Eleven "bottom-up" cluster are implemented in polyhedra (Figure 1) : tetrahedron (a) and double tetrahedron (b); pentagonal decahedron (c), Ino's decahedron (d) and Mark's decahedron (e); octahedron (f) and its derived truncated polyhedra, the truncated octahedron (g) and the cuboctahedron (h); the icosahedron (i), the truncated icosahedron (j) and the icosidodecahedron (k). Cuboctahedra and icosidodecahedra correspond to a strong truncation of the edges ( 1 /2) of their respective parent structure, whereas truncated octahedra and icosahedra are obtained by a weak truncation of the edges ( 1 /3). NPs can also be

1-20 | 3
shaped from bulk materials, as in the "top-down" approach. So far, the available lattices in polyhedra are : simple cubic, facecentered cubic, body-centered cubic, hexagonal close packed and β manganese phases (Figure 2a). Those bulk materials can be shaped into nanocrystals by (hkl) crystallographic planes at a given distance from origin. This is at the basis of the Wulff-von Laue method if the distance is proportional to the surface energy of the plane. It can yield a wide variety of nanostructures, such as nanorods, cubes,.... Spheres or ellipsoids with tunable radii can easily be obtained, leading to the sphere and to the oblate and prolate shapes shown in Figure 2b. In the same way, cluster with preferential growth directions can be designed, such as "cross"shaped nanoparticles (Figure 2c), by the reproduction of the bulk unit cell in a specific directions. In this case, branches length and thickness can be tuned. The polyhedra software can also model twining particles [START_REF] Hofmeister | Encyclopedia of Nanoscience and Nanotechnology[END_REF] . Using a set of user-defined twinning planes, polyhedra can alter the growth pattern of bulk materials leading to a newly ordered layer structure. That way, instead of stacking layers in the same order, a symmetry is obtained along the plane generating a new stacking pattern. For example, considering the fcc bulk structure ABCABC, an appropriate twinning plane leads to a CBABC stacking. This feature can be used on all the structures listed hereinabove. This is not the only pattern alteration implemented in polyhedra. As the hcp and fcc stacking match respectively in the (001) and (111) direction a compact stacking with random apparition of AB (hcp) and ABC (fcc) patterns can also be designed as bulk material that can be cut (Figure 2d).

X-ray or neutron powder diffraction and Reverse Monte-Carlo modeling

The global search for relevant minima on a first-principles PES is computationally prohibitive. Despite the limitations of pure theoretical methods, a joint theoretical and experimental approach, such as Reverse Monte Carlo (RMC) modeling of atomic pair distribution function (PDF) analysis of X-ray or neutron powder diffraction, is a powerful and relevant way to assess the structure of nanoclusters, since the PDF does not require periodicity. Actually, RMC modeling is a general method that provides atomic structures based on experimental data. 13,[START_REF] Mcgreevy | [END_REF] It was originally developed for applications to the structure of liquids and glasses but has now also been applied to crystalline and magnetic structures, 51 as well as to NPs. [52][53][54][55][56] The normalized scattering intensity S(Q) provided by X-Ray or neutron diffraction is the cornerstone of RMC modeling. The sine Fourier transformation of S(Q) yields the PDF, G(r), also called the radial distribution function (RDF):

G exp (r) = 4πr 2 [ρ(r) -ρ 0 (r)] = 1 2π ∞ 0 Q [(S(Q) -1] sin(Qr)dr (1)
where ρ(r) is the microscopy pair density and ρ 0 is the aver- age number density of the material (a rigorous definition of the negatively sloping baseline, -4πrρ 0 , has been given for finite-size objects such as NPs by 57 ). G exp (r) gives the probability of finding pairs of atoms at distances r. A RDF function can also be calculated from atomic coordinates, using:

G calc (r) = A ∑ i ∑ j b i b j < b > 2 δ (r -r i j ) -4πr 2 ρ 0 (2)
where r i j is the interatomic distance between two atoms i and j belonging to the model crystal, b i is the scattering power of atom i, < b > is the average scattering power of the sample and A is a fitting parameter for the amplitude of the signal. Alternatively to this so-called chemists definition, 58 another formula is often used, where G ′ (r) = G(r)/r. The delta function, δ (rr i j ) is replaced by a Gaussian distribution function of the form:

G(r -r i j ) = 1 √ 2πσ (r i j ) exp - 1 2 r -r i j σ (r i j ) 2
While in its simplest formulation, the peak width σ (r i j ) is set up to a constant value σ 0 , an r-dependence can be introduced to account for correlated motion and for the Q-resolution of the diffractometer. [START_REF] Neder | Diffuse Scattering and Defect Structure Simulations. A cook book using the program DIS-CUS[END_REF] Experimentally, RuNPs are usually observed to be hcp spherical nanocrystals. 4,59 Recently, although the fcc structure does not exist in the bulk Ru phase diagram, pure fcc RuNPs were obtained by the standard reduction method of various organometallic precursors. 60 Fivefold-symmetry twinned NPs 49 having a decahedral structure were directly observed in that work with high resolution transmission electron microscopy (HRTEM), and their crystal structure was investigated by powder X-ray diffraction (XRD) analysis. The calculated RDF profiles for ∼2.5 nm hcp and fcc RuNPs are shown in Figure 3. They both exhibit a first sharp and intense peak at 2.7 Å, which corresponds to the interatomic distance between nearest neighbors. Both profiles allow to discriminate between fcc and hcp at pair separations larger than ca. 5.5 Å, i.e. beyond second nearest-neighbors. The direct comparison between experimental RDF profiles and profiles calculated from structural models such as those provided by polyhedra (Figures 1 and2) allows discrimination between competing structural models. But, thanks to RMC modeling, it is possible to learn more from atomic pair distribution function analysis, in particular the level of local or global structural disorder. For example, it has been shown by Bedford et al. 52 that a non-crystallographic model is the only way to interpret PDF data of 2 nm thiol-capped RuNPs.

Provided that an experimental profile, G exp (r) has been determined, a structure can be produced by minimizing a goodness-offit function by a MC process:

χ 2 (r) = ∑ N k=1 G exp (r k ) -G calc (r k ) 2 ∑ N k=1 G exp (r k ) 2 (3)
This is similar to the Metropolis Monte-Carlo method, 61 but the main advantage is that an energy potential is not required in the standard formulation of the RMC method. Whereas a configuration space is explored by sampling energies in Metropolis MC simulations, in RMC modeling the configuration space is explored by minimizing χ 2 , in a manner reminiscent of the Metropolis MC method. As in any MC simulation, a global exploration of the cost function landscape can be set up by tuning the temperature parameter and the random changes of configurations.

Basin-Hopping 62 and simulated annealing 63 algorithms are implemented in polyhedra, in a similar way as proposed by Henkelman and co-workers 54 in the same context of RMC modeling. The RDF function is calculated as convolved gaussian functions with full width at half maximum σ 0 = 0.2 Å (i.e. the resolution of the experimental data). Owing to the finite-size of a NP, a baseline term consisting of this function convolved with a Gaussian function of width 1 Å is subtracted in eq. 2 instead of the function -4πrρ 0 . This term makes the G calc (r) signal vary around zero on this length scale and its limit is 0 beyond the diameter of the NP. Yet, nothing can prevent such minimization process to yield spurious structures. This is why structural constraints are imposed in RMC modeling, for example by fixing atomic coordination or neighbor constraints. 51,64 Another possibility to introduce constraint consists in using empirical potentials that provide the total energy E(r) of a given configuration r. The cost function to be minimized thus becomes

F(r) = ω E E(r) + ω χ 2 χ 2 (r)
. This is the strategy adopted in poly- hedra, in which Sutton-Chen 65,66 and Lennard-Jones potentials are implemented. Other constraints can be easily added, such as keeping the moments of inertia (see their definition in the Computational Details section): ω I ∑ 3 q=1 [I q -I q (ref)] 2 . Several geometry modification operators are implemented in polyhedra. One operator is randomly selected at each MC step. Some of them are identical to those used by Rondina and da Silva 67 : the cartesian displacement of atoms on a random direction, with a random magnitude; the random angular displacement atoms with the constraint that surface atoms are allowed to have greater displacement than core atoms; the displacement of an atom on a sphere centered on the atom of a radius equal to the distance between the atom and the center of the structure; the same operator but for surface atoms only. The same operator for a total of 7.5 % of the whole structure and the inclusion of a surface atom into the cluster core was added. We also implemented an exchange between two planes of atoms, in order to generate stacking defaults in the structure. Another new operator consists in moving under-coordinated atoms, in order to regain compactness. A "Rubik's cube" operator is also available. It consists in twisting one moiety of the NP and may lead to helical structures.

1-20 | 5
We are now going to demonstrate on a test case the capabilities of this method. Let us take a spherical 147-atom hcp RuNP as the reference. It is possible to retrieve this morphology starting from a cuboctahedral Ru 147 cluster (Ru 147 -CB) or from an icosahedral Ru 147 cluster (Ru 147 -CB) using the RMC method, despite the high stability of these structures, especially the icosahedron. As a matter of fact, it is worth noting that setting up ω χ 2 to zero involves a standard global search for the most stable Sutton-Chen isomer. Starting from any 147-atom structure, it almost systematically leads to the icosahedron. Figure 4a shows the initial RDF profile of the cuboctahedron and the one obtained at the end one of the RMC simulations, after setting up ω χ 2 such as the goodness- of-function represents approximately 1/3 of the total cost function F(r). The same conditions were applied starting from an icosahedral cluster (figure 4b). The comparison between the reference RDF profile and the RMC-optimized structure profile is very good. For the cuboctahedron, the characteristic peaks of an fcc stacking, as the one found at 6 Å, decay to let the hcp signature profile emerge. The geometry of the resulting nanoparticle can be compared to the reference one, plotted in figure 4c. They are not exactly the same, but it is safe to assert that the hcp stacking is found and that a spherical shape is obtained. This is even truer when the icosahedron is considered as a starting point as this time an almost perfectly spherical morphology is obtained and the agreement between the reference RDF profile and the final structure RDF profile is even better than in the case of Ru 147 -CB. Note that the quite good overlap between the Ru 147 -IC and reference RDF profiles makes the global optimization search more challenging than using Ru 147 -CB as a starting point. According to the qualitative concept of Paul Sabatier, catalytic properties will be hindered if the reactants adsorb too strongly, whereas no reaction will occur if the interaction is too weak. On the other hand, the Brønsted-Evans-Polanyi (BEP) equation 17,18 relates the change in activation energy of a reaction to the change in the reaction energy:

δ E act = αδ E r .
As a consequence of the Sabatier principle and of the BEP relationship, the catalytic activity of heterogeneous catalysts within the periodic table is well described by Balandin's volcano relations between reaction rates and adsorption energies. 16,[START_REF] Balandin | Advances in Catalysis[END_REF][68][START_REF] Swiegers | Mechanical Catalysis: Methods of Enzymatic, Homogeneous, and Heterogeneous Catalysis[END_REF][START_REF] Tedsree | [END_REF] Volcano plots have been extensively and successfully used in the analysis of surface elementary steps of heterogeneous catalytic processes. [71][72][73] Even if calculating adsorption strengths on the active catalytic sites of NPs is less time-consuming and tedious than finding transition states, evaluating adsorption processes on all the active sites of a metal NP is a huge task. The DFT-based design of catalytic NPs thus needs to define a relevant and easyto-calculate index for adsorption strength, at low computational and human cost. Whereas coordination numbers of the metal atoms calculated from a 3D model are probably the simplest and cheapest reactivity indexes, as we will see in section 4.2, in a previous paper, [START_REF] Del Rosal | [END_REF] we chose to derive a variant of the d-band center. It requires to calculate the DOS and to project it on a local atomic basis set. We showed that the adsorption strength for a given coordination site µ k can be assessed by calculating an index, εd (µ k ).

It can be nicely depicted as a color map, which provides a first overview of the overall d-accessibility of a µ k site. It is calcu- lated as the normalized, energy-weighted integral of the density of states (DOS), projected onto all d atomic orbitals (AOs) of the surface atoms which characterize the µ k coordination site. The maps given in this paper are obtained as a superposition of the εd values calculated for atop, edge-bridging, face-and squarecapping sites. We have shown in our seminal paper that this model can be refined by considering the PDOS for individual d AOs. [START_REF] Del Rosal | [END_REF] However, the εd (µ k ) index already provides a good estimation of the adsorption strength of ligands on the surface, and it will be the only considered index in the present paper, which only aims at roughly comparing the electronic fingerprint of various naked and dressed nanoclusters. It will also later be compared to generalized coordination numbers [21][22][23] (CN) in the case of Ru 55 .

Comparison between Co 55 and Ru 55 . For a sake of a rigorous comparison, we have considered in the present paper the same Marks-decahedron-like structure than in our seminal work. Thus, a Marks-decahedron-like cobalt nanoparticle (Co 55 -MD) has been optimized and its d-band center has been computed. As can be seen on Figure 5, even if the two nanoparticles are equally shaped, they have a very different color map. This color map gives an immediate view of the possible spots for strong (dark red), normal (white) and weak (dark blue) adsorption with respect to Ru(0001) and Co(111), the origin of the ε d energy scale being chosen as the value calculated for surface atoms of respectively slabs. Firstly, the ruthenium displays a wider range of ε d values : centered on the Ru(0001) hcp surface ε d value of -2.54 eV, the d-band center spreads from -3.17 eV to -1.91 eV. For ruthenium it is noteworthy that the ε d color map highlights the pres- ence of different special sites at the surface of this cluster, such as a butterfly like pattern. The cobalt nanoparticle is very different. The ε d color map depicts a situation where no highly reactive sites can be found. The principle for the scaling is the same as for the ruthenium, except the scale is made narrower because of the low d-band center variation of the nanoparticle : the mapping is centered on the ε d value of -1.88 eV obtained for the Co(111) fcc surface. For Co 55 -MD, the ε d values range from -2.17 eV to -1.59 eV. Some dark blue spots are visible evidencing sites on which the adsorption energy should be lower than at the Co(111) surface. White sites, with an adsorption energy equivalent to the Co(111) surface are prevalent. No dark red sites can be spotted out on the map. Only some light pink sites, with adsorption energies barely higher than on Co(111), can be located on the surface of Co 55 -MD. The correlation between εd (µ k ) and adsorption energies of ligands will be the scope of a separate paper.

The molecular orbital language

The molecular orbitals (MOs) model is very successful in explaining various chemical properties and for providing trends. Other models based on a MO description or on the direct analysis of the electronic density also arose (such as NBO 74,75 , QTAIM, 76 ELF, 77,78 ...), some of them aiming at giving an information localized on few sets of atoms (e.g. involved in a chemical bond).

Such methods are of course relevant, but MO analysis already provides valuable local informations that can explain local properties of a molecule, such as the regioslectivity of a catalyst. It has been promoted as a cornerstone of molecular chemistry owing to its ability to rationalize a broad range of molecular properties. [START_REF] Fukui | Theory of orientation and stereoselection[END_REF][START_REF] Hoffmann | [END_REF] Moreover, Hoffmann's view of bonding in solid state chemistry 81 fills the gap between MO theory and the language of band theory in solid state physics. Canadell et al. [START_REF] Canadell | Orbital Approach to the Electronic Structure of Solids[END_REF] have also significantly contributed to demystify crystal orbitals (CO), and to strengthen the bridge between molecular chemistry and solid state chemistry. Whereas in both cases, the underlying Hamiltonian is of extended Hückel or tight-binding type, Kohn-Sham MOs or COs can be analyzed as well.

A lot of first principles calculations of crystals, surface models or NPs, i.e. calculations involving periodic boundary conditions, are performed within plane-wave basis sets. Combined with pseudopotentials and thanks to Fast Fourier Transform such basis sets are very efficient because they allow to work in the reciprocal space. Unfortunately, the quantum chemical information contained in the total wavefunction is not easy to handle. This is why Hoffmann proposed the crystal orbital overlap population (COOP) concept, [START_REF] Hughbanks | [END_REF] which provides an overview of the bond strength in crystals and of the bonding, non-bonding or anti-bonding character of the states. The crystal orbital Hamilton population (COHP) was then introduced in the framework of DFT by Dronskowski and Blöchl. 24 COHP is a partitioning of the band-structure energy in terms of orbital-pair contributions, and it is therefore based on a local basis set (eq. 4).

COHP µν (ε) = ∑ nk c * µ c ν H µν δ (ε -ε n (k)) (4)
The Lobster package which has been used in this work 25,26 allows to calculate COHP curves projected in a local atomic basis set (pCOHP), and also reliable atom-projected density of states (pDOS), both directly based on plane-wave wavefunction calculated with the VASP package. Whereas in the case of small molecules COOP or COHP profiles (calling them MOOP and MOHP would be more rigorous) do not bring much with respect to the close examination of MOs, it becomes an interesting approach for large systems. We have recently demonstrated the usefulness of this approach in the case of RuNPs. 12 We are now going to confirm the accuracy of the pDOS and pCOHP analysis provided by Lobster on Ru 4 (µ 3 -H) 4 (η 6 -C 6 H 6 ) 4 , which belongs to the cluster realm, thoroughly investigated for decades [84][85][86][87] . This 60-electron Ru 4 H 4 L m tetrahedral cluster, stable according to the Mingos-Wade rules, 88,[START_REF] Mingos | Introduction to cluster chemistry[END_REF] exhibits six skeletal MOs as well as four Ru-H MOs. An energy-level diagram calculated with Gaus-sian09 at the DFT-PBE level of theory is reported in Figure 6, together with the pDOS and pCOHP profiles calculated with VASP and Lobster at the same level of theory. We shall now focus on some remarkable MOs. According to the MOs plot in Figure 6a, which analysis may be slightly biased by the cutoff value for the isosurface, the two degenerate HOMOs are slightly bonding and span ruthenium and carbon atoms. The six so-called skeletal MOs lie below two sets of triply degenerate MOs which exhibit an antibonding character between Ru atoms. The first totally symmetric Ru-H MO, which lies approximately at -8 eV, also has a strong component on the π MOs of the η 6 -benzene rings. The same analysis performed from the pDOS and pCOHP profiles (Figure 6b) is less straightforward, but it has the added benefit to provide a well-founded identification of the bonding vs. antibonding character between a pair of atoms in complex systems, i.e. it is a way to get rid of the bias introduced by the isosurface cutoff value in MO plots. The pDOS profiles show that the MOs that lie between the Fermi energy and -6 eV essentially span the Ru (green curve) and C (red curve) atoms, as observed on the isosurface plots. Among them, the skeletal bonding MOs can easily be identified from the pCOHP Ru-Ru profile (green curve). A black peak rises on the bonding side of the pCOHP plot at ca. -8.5 eV, together with Ru-Ru and C-C bonding interactions. It indicates the location of the Ru-H MO. Beyond these specific cases, a careful comparison between the MOs and associated energy levels on one hand and the pCOHP(ε) and pDOS(ε) profiles on the other hand does not reveal a wrong projection of the wavefunction on the local basis set. Owing to the relevant and accurate description of the bonding provided by the Lobster package on this test case, it can be safely applied to more realistic models of NPs, as we will be exemplified in section 4.2 (see Cusinato et al. 12 ).

Surface species

Motivation

Ligands spread over the surface of colloidal NPs according to electronic as well as to steric effects. As a consequence of the electrondeficient nature of surface atoms in the case of middle transition metals, the optimal number of surface ligands is to some extent ruled by steric hindrance. Let us consider now the case of the Fischer-Tropsch reaction catalyzed by ruthenium NPs. [START_REF] Martínez-Prieto | [END_REF] Fischer-Tropsch synthesis (FTS) is a catalytic process that converts a mixture of carbon monoxide and dihydrogen in the gas phase into mainly linear hydrocarbons and water (eq. 5). It is known for around ninety years now 91 and it is still of high interest because of its ability to produce fuels with lower environmental impact, but higher financial costs, than regular production pathways. 92,93 In a recent paper, 12 we have shown that according to DFT energies only a very low number of hydrides or carbon monoxide ligands should adsorb on the surface of the Ru 55 -hcp cluster, in contradiction with titration experiments 9,94 . It was also demonstrated that the weakening of the NP-ligand bond strength as the coverage value increases is closely related to the stabilization of the d-band center with respect to the Fermi energy. Let us now consider the bis(diphenylphosphino)butane (dppb) which is one of the stabilizers used by Mart\'{i}nez-Prieto et al. [START_REF] Martínez-Prieto | [END_REF] This bidentate ligand is not totally innocent with respect to the FTS process, as it does not only act as a capping and protecting agent but it also stays on the surface of RuNPs during reaction and it improves their activity and selectivity toward C2-C4 hydrocarbons. On a bare 55-atoms RuNP, the adsorption of dppb is largely exothermic by ca. -70 to -125 kcal.mol -1 . The most stable isomer is shown in Figure 7a (E ads : -124.6 kcal.mol -1 ). It is characterized by a gauche conformation of the butane arm, by a geminal grafting of the diphenylphosphino ligands on an apex reinforced by a π-bonding of two phenyl rings with the metal surface. The adsorption energy on an hydrogenated cluster (Ru 55 H 64 , H-coverage value: 1.5, Figure 7b) becomes significantly less exothermic (E ads : -51.4 kcal.mol -1 ), both owing to the steric hindrance of the hydrogen atoms which prevent the πbonding and to the lowering of the Ru-P bond strength. An anti conformation of the butane arm is also possible, for example on an edge, as shown in Figures 7c and7d. The adsorption energy on the bare RuNP can be estimated to be 37 kcal.mol -1 per Ru-P bond, whereas it is slightly lowered in its hydrogenated counterpart (-35 kcal.mol -1 per Ru-P bond).

(2n + 1)H 2 + nCO -→ C n H 2n+2 + nH 2 O (5)

Electronic structure

Effect of ligands on the electronic signature of the metal core. In contrast to the use of ligands to decorate metal nanoparticles and modulate their surface chemistry, the interaction of metal complexes with the surface of nanoparticles is not well developed. The modes of coordination of tin species on nanoparticles were recently studied by several spectroscopic techniques and by DFT as well, whereas a simple catalytic test, styrene hydrogena- tion, was also achieved. 4 The d-band center map of 55-atom hcp nanocluster is shown in Figure 8a. It also exhibits sites with an expected strong adsorption ability. On this basis, let us now consider a Ru 55 (SnMe) 9 nanocluster (Figure 8b). With the exception of the ruthenium atoms bound to the SnMe ligand, the d-band center map looks very similar to the bare nanocluster. It involves a local influence only of the nine methyltin complexes on the electronic structure of the ruthenium core. Let us now consider the same Ru 55 -HCP cluster, but now fully saturated with 1.5 CO per surface ruthenium atom, a value that correspond to a realistic coverage at r.t. and under a 1 bar pressure of CO. 12 The average d-band center of the surface atoms is lowered by ca. 0.5 eV (-3.05 eV vs. -2.58 eV), which, on an electronic basis only, makes ruthenium surface atoms weaker grafting sites than naked surfaces. On top of that, surface species may as well involve steric effects that will further impede the coordination of additional ligands. pDOS and pCOHP profiles. The electronic structure of another Ru 55 H 64 (dppb) isomer is reported in Figure 9 in terms of pDOS and pCOHP profiles. The geminal coordination of the two phosphorus groups and the π interaction of one phenyl ring are responsible for a relatively strong coordination of the dppb ligand (adsorption energy: -56.4 kcal.mol -1 ). Let's start tracing the surface-ligands interactions. The metal d-band dominates the states close to the Fermi energy. According to pCOHP, the highest occupied ruthenium states are anti-bonding, as already observed in the Ru 4 H 4 (C 6 H 6 ) 4 cluster (see Figure 6). The high number of surface hydrides make the pCOHP(Ru-H) and pDOS(H) peaks almost as high as the pCOHP(Ru-Ru) and pDOS(Ru) peaks, whereas the Ru-P and Ru-C curves barely arise, owing to the small number of these interactions with respect to Ru-Ru or Ru-H. The bond strength can be assessed by integrating the pCOHP function up to the Fermi energy (IpCOHP = E F E min pCOHP(ε)dε). It provides a qualitative information of the average bond strength in covalent systems. Such one-electron derived energy must not be compared to bond dissociation energies (BDE) but it gives a valuable bond strength index which variation is expected to follow BDE variations. It has been shown by Cusinato et al. 12 that in the case of carbonylated clusters it accounts nicely for the CO stretching frequencies. Moreover it is obtained with a single calculation, whereas the determination of a BDE require to perform three calculations. Whereas IpCOHP is found to be 39 kcal.mol -1 per each Ru-Ru interaction, the Ru-P turns out to be much stronger, as expected (IpCOHP: 111 kcal.mol -1 ) and the Ru-phenyl interaction is characterized by a large IpCOHP: 232 kcal.mol -1 , in line with the energy differences between the dppb complexes shown in Figures 7a and7c.

d-band center and generalized coordination number.

Whereas the d-band center model requires to compute the density of states of the nanoparticle, other less computationally expensive descriptors have been developed so as to be able to link the adsorption strength to the structure of the nanoparticles. One of them is the generalized coordination number CN, computed for an atom i as the weighted sum of the neighboring j atoms coordination number cn, cn max being the maximal coordination number in the bulk (eq. 6). 21,22

CN(i) = n ∑ j=1 cn( j) cn max ( 6
)
In that way, the generalized coordination number is a structural descriptor that takes into account the spatial environment of the considered site in order to refine the value traditionally obtained for coordination number. It was used to determine the optimal geometry for an active site in the case of hydrogen reduction, showing that cavities dug into fcc (111) platinum surface exhibit a generalized coordination number close to the optimal value expected for this reaction (CN=8.3). 22 In another paper, CN was shown to perform better than d-band center when applied to the correlation between adsorption energies of O, O 2 , OH, OOH, H 2 O and H 2 O 2 on platinum surfaces. We report here an example in which the d-band center model performs better than generalized coordination numbers, for the case of the adsorption of H 2 or CO on the Ru 55 -MD nanocluster previously investigated in Ref. [START_REF] Del Rosal | [END_REF] . Dissociative adsorption energies of H 2 (in black) and CO (in violet) are reported as a function of ε d in figure 10a and as a function of of CN in figure 10b. The correlation between dissociative adsorption energies of H 2 and d-band centers is very good, with a regression coefficient of 0.96. The generalized coordination number method performs fairly well, with a 0.80 regression coefficient. This shows that, in that case, an electronic descriptor give a more suitable description than a structural one. The CO case is a more stringent test because, as previously described by del Rosal et al., [START_REF] Del Rosal | [END_REF] there is a competition between σ donation and π back-donation that differs according to the coordination site. In some cases the averaged dband center value should be replaced by the projection onto the In that way, the electronic descriptor of the interaction between CO and the metal atoms of the nanoparticle leads to a very good correlation between adsorption energies and ε d value. Such a distinction is not possible with the generalized coordination number. Even if it takes into account the neighboring spatial environment of the coordination sites, it does not consider the way the ligand is coordinated nor the influence of the local electronic structure, which is of importance in the coordination of ligands, especially for sites that are not flat or perfectly shaped, like the ones found at the surface of the Ru 55 MD nanocluster considered here. As a result, the d-band center index performs better than generalized coordination numbers in the case of CO adsorption, as the regression coefficient for CN is only r = 0.68, slightly lower than the averaged d-band center, and significantly less good than our refined d-band center model. Color maps are also plot for each index. Only red sites appear on the CN map, which means that this criterion somewhat overestimates the adsorption strength since it systematically predicts an increased adsorption energy with respect to Ru(0001), in contradiction on some sites with the calculated energies. Overall, whereas generalized coordination numbers performed well in the case of platinum nanocatalysts, they do not appear to be the best adsorption index for this RuNP. Moreover, it is not clear how this generalized coordination number could be extended to surfaces covered by ligands.

Steric-driven distribution of ligands on the surface

It can be assumed that surface ligands are uniformly distributed on the surface, provided that the grafting sites are uniformly distributed as well. This can be related to the general mathematical problem of positioning a large number of points uniformly over the surface of a sphere. 95 Such points can be generated via the minimization of a relevant energy. It is exactly the same as locating on a sphere identical point charges that are in equilibrium with respect to the Coulomb potential. Such points minimize the electrostatic repulsion energy E = ∑ 1≤i< j≤N x i -x j -1 . This also the underlying scheme of the popular VSEPR method in chemistry, 96 which is used to explain the optimal shape of simple molecules in terms of electrostatic repulsion between electron pairs.

More generally, it is also possible to consider points that minimize the s-energy

E s = ∑ 1≤i< j≤N x i -x j -s , s > 0 (7)
As s → ∞, the s-energy mainly depends on short-range interactions and the problem of finding a uniform distribution of points tends to the so-called best packing of hard spheres problem. 95 There is an abundant mathematical literature dealing with spherical point optimal arrangement, and several algorithms have been proposed. But finding the absolute optimal configuration of ligands on the surface can be considered as meaningless, in particular because of their dynamical motion on the surface and the permanent adsorption/desorption processes under experimental conditions. Less refined algorithms can be used, such as a Monte-Carlo simulated annealing, 63,[START_REF] Binder | Monte-Carlo simulation in statistical physics[END_REF] which is known to lead to a very good approximation of the optimal configuration in most cases. Although in the general case NPs do not have a strictly spherical shape, some of them even being far from spheres, the mathematical background of this problem, i.e. minimizing the s-energy, can be used in order to design an efficient algorithm aiming at finding relevant geometrical configurations and a steric-driven titration of ligands grafted at the surface of metal NPs. Such information can be directly compared to the experimental titration of ligands on the surface, usually achieved by performing well-known chemical reactions (e.g. hydrogenation of olefins) involving the ligands (e.g. hydrides).and to measure the amount of products (e.g. alkanes) by chromatography or any other analytic technique. 9 Thus, in addition to the polyhedra program, the dressNPs homemade utility has been developed in our group, as a part of a multiscale strategy aiming at decorating the surface of metal NPs with ligands. The algorithm implemented in dressNPs is related to a coarse-grained method, the metal core and the ligands being kept frozen during this process. DressNPs first identifies surface atoms, The points correspond to the individual values of adsorption energies whereas the line is obtained by linear regression from these points. r, the regression coefficient, is given in each case. The colors of the ε d and CN 3D activity maps cannot be directly compared owing to the arbitrariness of the scales. In each case, the Ru(0001) surface is chosen as the reference value, in white. prior to an analysis of the surface in order to identify apexes, edges, µ 3 sites,. . . The user also has to define which ligands will be grafted on the surface, and their coordination mode. The optimal geometry of the ligands in interaction with a metal surface must have been previously optimized, for example at the DFT level of theory. We chose to make a library of ligands whose geometry were optimized with the DFT-B3PW91 functional implemented in Gaussian09, [START_REF] Frisch | Gaussian 09 Revision D.01[END_REF] using a 6-31G(d,p) basis set. The typical metal-ligand distance must also be set up. The steric hindrance index E s is then minimized by Monte Carlo simulated annealing, in order to find the optimal distribution of a given number of ligands, N. Finally, the optimal number of ligands, N * , that could be grafted on the surface without being sterically discomforted can be found that running the algorithm for various values of N. We have recently applied this algorithm in the case of tin-decorated ruthenium NPs, and found N * values close to experiments. 4 An example of a typical decoration of a 1.9 nm hcp RuNP with B 5 sites 99 is shown in Figure 11. It is also possible to add terms to the s-energy (eq. 7) for example in order to uniformly distribute a given ligand and to force the segregation of another type of ligands. It can be achieved by increasing the repulsive term between two different surface species.

1-20 | 11
It becomes easy to graft surface species with various coverage values and on different coordination sites that can then be optimized at the DFT level of theory. Such high-throughput generation has facilitated the calculation of the ab initio phase diagrams presented in the following subsection.

From microscopic to macroscopic scale: ab initio thermodynamics 4.4.1 Methodological background

It is possible to calculate the free energy and other thermodynamic functions of solids and liquids, using first principles methods. Such method has successfully been applied to explain or predict thermodynamic material properties, and in particular surface properties at the solid-gas interface. [100][101][102][103][104] In practice, it extends the T = 0K and p = 0Pa ab initio energies of surfaces covered by ligands coming from the surrounding medium to realistic environmental conditions in terms of temperature, pressure and composition of the gas phase. Let us now consider the usual field of application, the adsorption of organic species on clean crystallographic planes. The adsorption Gibbs free energy per unit area of n ligands on the bare (hkl) surface, ∆ a G hkl (nL, T, p), requires to calculate the Gibbs free energy for the reaction M(hkl) + nL = nL * (hkl), with the bare surface chosen as the reference:

∆ a G hkl (nL, T, p) = µ(nL * , T, p) -µ(M hkl , T, p) -nµ(L, T, p) 2A hkl (8)
∆ a G is normalized per unit area by dividing through the surface area 2A hkl . The temperature and pressure effect for ligand (L) adsorption on an (hkl) crystallographic plane for a given metal (M) will be introduced by the appropriate definition of the chemical potentials. The surface energy of (hkl) facets dressed by n ligands L can be expressed as: 102,105

γ hkl = γ (0) hkl + ∆ a G hkl (nL, T, p) (9)
While chemical potentials for solids weakly depend on pressure, the influence of temperature should be taken into account. Nevertheless, it is assumed that for condensed phases the difference µ(nL * ) -µ(M hkl ) can be approximated by the difference between the calculated electronic energy including (∆G • , eq. 10) or not (∆E, eq. 11) thermal variations of internal energies, as well as pV terms and entropy contributions:

∆ a G = [∆G • -nµ(L, T, p)] /(2A hkl ) (10)
∆ a G = [∆E -nµ(L, T, p)] /(2A hkl ) (11)
As the surrounding medium acts as a reservoir of ligands, the chemical potential for the ligand L can then be calculated from the standard chemical potential and the activity of the ligand: µ(L, T, p) = µ ⊖ (L, T, p • ) + kT ln a(L)

At the gas-phase -solid-surface interface, the p-dependence is introduced in the L chemical potential through an ideal gas law:

a(L) = p(L) p • (13)
µ ⊖ (L, T, p • ) can usually be calculated from H • T and S • T values given in thermodynamic tables 106

(µ ⊖ (L, T, p • ) = H • T (L) - T S • T (L))
or computed from first-principles calculations done at 0K. In this case the standard chemical potential, i.e. G • T (L), is given by: 107

µ ⊖ (L, T, p • ) = E DFT (L) + H • T (L) -T S • T (L) (14)
where H • T (L) and S • T (L) are usually calculated from the L partition function, i.e. by a straightforward application of the statistical thermodynamic equations. 107 In the case of colloidal nanoparticles with low concentration of ligands in the solvent medium, L ligands are considered as a solute and their activity can be approximated to:

a(L) = c(L) c • (15)
The chemical potential is defined with respect to the ideal infinitely dilute reference state: µ ⊖ (L, T, p • ) = µ • ∞ (L, T ). This reference chemical potential can as well be found in reference thermodynamic tables or can also be determined with a good accuracy by first-principle calculations performed on the ligand L embedded in the so-called polarizable continuum medium (PCM). [108][109][110] In the next examples, the vibrational contribution to the free Gibbs energy is systematically taken into account:

∆G • = E DFT (nL * ) + F vib (nL * ) + F ′ vib (M hkl ) (16
)

-E DFT (M hkl ) + F vib (M hkl ) (17)
where F vib (M hkl ) and F ′ vib (M hkl ) stand for the vibrational contribution of the metal surface atoms either without and with adsorbed ligands. By considering these two quantities to be almost equal, eq. 16 becomes:

∆G • = ∆E + F vib (nL * ) = ∆E + H vib (nL * ) -T S vib (nL * ) (18)
Let us now consider a co-adsorption process of two species, L 1 and L 2 , which is the starting point of the Langmuir-Hinshelwood mechanism in heterogeneous catalysis. It can be summarized as

M(hkl) + n 1 L 1 + n 2 L 2 = n 1 L * 1 (hkl) + n 2 L * 2 (hkl)
, and the Gibbs free energy for this reaction is calculated as:

∆ a G = [∆G • -n 1 µ(L 1 , T, p) -n 2 µ(L 2 , T, p)]/(2A hkl ) (19)
with:

∆G • = ∆E + F vib (n 1 L * 1 ) + F vib (n 2 L * 2 ) (20)
Surface free energies can be plotted as a function of the two chemical potentials µ(L 1 ) and µ(L 2 ) for various compositions (n 1 , n 2 ) and plausible geometries and can then be summarized as a surface phase diagram in the [µ(L 1 ), µ(L 2 )] space. It has for example been applied to determine the DFT phase diagram of surface structures for the model catalyst RuO 2 (110) in contact with a gas environment of O 2 and CO. 111,112 As underlined by the authors of these studies, knowing the various stable surface phases in equilibrium with the given gas phase is an important prerequisite.

Application on an hydrogenated 55-atoms ruthenium nanocluster

It is thus possible to combine DFT energy calculations with the chemical potential in the gas phase or in solution in order to

1-20 | 13
relate computed energies to realistic experimental conditions in terms of temperature and pressure. We reported earlier an application of this ab initio thermodynamics approach to the dissociative adsorption of H 2 on a bare Ru(0001) surface. 104 The main conclusion arising from the computed phase diagrams was that it is unlikely to observe a H coverage greater than 1 H / Ru surface . Yet, it was recalled in the introduction that according to titration experiments, small RuNPs stabilized by polyvinylpyrrolidone (RuNP@PVP) can be covered by at least 1.3 hydride per surface ruthenium atom. 9 These two different values must probably not be considered as a discrepancy between theory and experiments. They suggest that it is not relevant to extrapolate results obtained on a compact surface to the case of Ruthenium NPs for which the presence of numerous defects may certainly have a direct influence on the hydrogen coverage. This is what we are now going to examine, following the same computational methodology as in Refs. 104 and 12 . The variation of the Gibbs free energies of adsorption, ∆ a G (eqs. 10 and 16), at different temperatures and for different coverage values is given in Figure 12. They are compared to the adsorption on a Ru(0001) slab. 104 Below 520 K the compact Ru(0001) surface can at best accommodate 1H per surface ruthenium atom. The surface of a 1 nm RuNP without ancillary ligands, i.e. typically a RuNP@PVP 6 , is saturated by 1.6 H / Ru surface below 350K, and still covered by 1.2 H / Ru surface up to 550 K. This optimal coverage value diminishes up to 780 K, whereas according to these calculations, a bare Ru(0001) surface should already be obtained at 640 K. This comparative study confirms that this methodology provides thermodynamical data in agreement with experiments, with an optimal coverage value comprised between 1.2 and 1.6 H / Ru surface . The saturation threshold on small RuNPs being significantly larger than on a ruthenium surface, it involves that the use of finite size models is to be preferred if one wants to reach a reliable description of the properties of small metal nanoparticles.

55-atoms Ru cluster stabilized by carbonyl compounds

The variation of the Gibbs free energies of adsorption, ∆ a G (eqs. 10 and 16), under p=1 bar but at different temperatures and for different coverage values is given in Figure 13a. The variation of ∆ a G as a function of the CO pressure is plotted for T = 400 K in Figure 13b. Those are two different way of representing the evolution of the coverage as a function of a thermodynamical parameter and, in particular it allows to see the evolution of ∆ a G for all the considered structure and does not just define domains for the minimal ones like in a traditional (T, p) phase diagram. In that way, it can be seen that the structure with 1.74 ML, in light brown on the diagrams, which is not the most stable at any of the considered points, is quite close to be minimal at 200 K for p CO = 1 bar and could be in competition with 1.53 ML, whereas the second plot shows that such competition does not exist at T = 400 K, even at high pressures. The straightforward result that, un- der standard thermodynamical conditions, the optimal coverage is 1.53 ML is also visible on those diagrams. 12,94

Co-adsorption of hydrides and carbonyl compounds on Ru 55

Considering the adsorption of a single kind of ligand in equilibrium with a nanoparticle is generally not enough to describe a catalytic process, not to mention the formation of transient species which may accumulate in the chemical medium. For example, in the Fischer-Tropsch synthesis two reactants are needed in order to produce the resulting hydrocarbons. In the experimental set-up, the metallic nanoparticles are placed in equilibrium with syngas, generally 3 bar (1:1 molar mixture of H 2 and CO). It is therefore necessary, in order to define a theoretical model reproducing the expected nanoparticles obtained under these experimental conditions, to take into account the nanoparticle in equilibrium with both H 2 and CO. Phase diagrams depending on the pressure of H 2 (p H 2 ) and the pressure of CO (p CO ) at a given temperature for the co-adsorption of H 2 and CO on small ruthenium nanoparticles (∼ 1 nm) have shown that it was hard to reach co-adsorption of the two species under FTS standard pressures and temperature (∼ 450 K) and that CO only covered nanoparticles are expected under those conditions. 12 If the co-adsorption of H 2 and CO is not thermodynamically favorable at p H 2 = 1 bar, p CO =1 bar and T = 450 K, one can try to find which thermodynamical conditions allow to reach this situation while remaining experimentally real- istic. Considering a phase diagram depending on the temperature and on the pressure of one ligand, the pressure of the other one being kept constant, is one way to figure it out. The figure 14 shows such a phase diagram, depending on (T, p CO ) for p H 2 = 1 bar. From this diagram, the maximal pressure of carbon monoxide ensuring H 2 /CO co-adsorption can be read for temperature ranging from 0 to 1000 K, and especially from 350 to 450 K, range in which FTS occurs on small ruthenium nanoparticle. At 350 K the pressure of CO must be kept under 10 -14 bar in order to reach a co-adsorption domain, namely 0.25 H / 0.75 CO, whereas at 450 K p CO can be up to 10 -10 bar in order to reach that same area of the phase diagram. Those are very low pressures, similar to high vacuum ones, and therefore not reachable in a FTS process.

Conclusion

Inorganic metal nanoparticles are more complicated to accurately characterize than inorganic metal clusters. For small metallic clusters such as Ru 4 H 4 (C 6 H 6 ) 4 , there is a unique thermodynamically optimal isomer. It can be experimentally characterized by liquid-state NMR and infra-red spectroscopies. Since it crystallizes, it can be as well characterized by X-ray experiments and by solid-state NMR. It is this specific geometry that has to be modeled as a starting point for reactivity studies or electronic structure analysis. In the case of NPs, the focus is less on being able to mimic an hypothetical unique and exact experimental geometry than on being able to extract information on the general morphological features of the synthesized NPs and to reproduce them in silico. In that way, the rationalization of the morphology, of the electronic properties and coverage of NPs under experimentally realistic conditions turns out to be a crucial point in order to design relevant models for theoretical studies. We have presented here a multiscale framework (polyhedra, dressNPs and aithermo) aiming at helping such a rationalization. A new way to build a wide variety of the typical shapes exhibited by the nanoparticles completed by a steric-driven grafting of different ligands on its surface has been presented in this paper. A

1-20 | 15
RMC procedure that allows to refine geometries to match welldefined experimental RDF profiles is proposed in addition to this systematic generator. By using it, it has been possible to find a structure that perfectly matches the spherical Ru 147 -HCP starting from an icosahedral geometry. We have also compared two 55-atoms nanoclusters, Ru 55 -MD and Co 55 -MD, in light of the d-band center model. The local adsorption strength at the surface of both clusters has been analyzed and interpreted by an on-site d-band center descriptor depicted as a color map. According to those color maps, that give a straightforward point of view of the expected adsorption strength at the surface of the clusters with respect to a reference value (in this case the close-packed (001) plane), the reactivity of Co 55 -MD is expected to be very different from its ruthenium counterpart, even though the structures are very similar. Comparison between d-band center and generalized coordination number shows that the lack of directionality of the latter causes it to perform worse than d-band center, regarding the correlation with adsorption energies. A chemical bond analysis tool, the COHP index, is also applied to NPs so as to fill the gap between chemical knowledge of molecular orbitals and band structures obtained via calculations on metallic NPs. A Monte Carlo process for grafting ligands on the surface on NPs by taking into account steric repulsions between surface species is described, leading to a wide range of grafted NPs. The relative stability of those ligand-covered NPs is then studied via ab initio thermodynamics, exemplified in the experimental context of FTS by evidencing the stability of CO-only covered NPs. The combined use of those structural, electronic and thermodynamic approaches is expected to play a significant role in advancing our knowledge of complex materials.

Computational Details

Regarding DFT calculations on large metal clusters, i.e. species that exhibit a HOMO-LUMO degeneracy at the Fermi level, periodic DFT calculations have been achieved, using a GGA functional (see below). It is usually considered as a good compromise between accuracy and computational efficiency, so it has extensively been employed in particular when the metal cluster size is large. [START_REF] Gerber | Nanocatalysis with Transition Metal Particles: Where are we Now[END_REF][START_REF] Carchini | [END_REF] It has also been stated by other authors 113 that no semilocal functional is capable of describing properly material properties (adsorption and atomization energies, d-band center, equilibrium volume of the bulk) and including non-local exchange also only improves some but worsens other properties.

Periodic DFT calculations of metal clusters. Software: Vienna ab initio simulation package, VASP.

εd (µ k ) = ∑ α∈µ k ∑ m E F E min εn d m (α, ε)dε ∑ α∈µ k ∑ m E F E min n d m (α, ε)dε (21)
where m runs over the five d AOs and n d m (α, ε) is the atomprojected density of states on the d m AO of atom α; µ k is reminis- cent of the symbol which designates bridging ligands in coordination chemistry and the bar sign above ε d means that it is averaged over all d AOs; E min is readily set to the bottom of the occupied dband; DOS integrated up to the Fermi level E F (see the discussion in ref. [START_REF] Del Rosal | [END_REF] ); all calculated values are plotted as colored maps: hot spots for adsorption in red, blue for potentially weak interactions, white shows where intermediate adsorption processes should occur (white is defined by the electronic feature of the Ru(0001) or Co(111) slab. 197.65 J.mol -1 .K -1 ); the chemical potential for H is calculated as 1 2 µ(H 2 , T, p), since it is the dissociative adsorption of H 2 that is considered. The surface area of the Ru 55 cluster was calculated by dressNPs, after identification of core and surface atoms (A = 258 Å 2 ).

Moments of inertia.

They are the eigenvalues I n (n = 1 -3) of the inertia tensor: 123

I qs = 1 N N ∑ i=1 (q i -q G )(s i -s G ) (22)
where q, s = x, y, z and G refers to the center of mass of the NP. The eigenvalues I n are related to the lengths a, b and c of the semiprincipal axes of the ellipsoid that contains the NP by: a

= 2 √ 5I 1 , b = 2 √ 5I 2 , c = 2 √ 5I 3 .

  Shape, electronic structure and steric effects of organometallic nanocatalysts: relevant tools to improve the synergy between theory and experiments Lucy Cusinato, a Iker del Rosal a and Romuald Poteau * a

Fig. 1

 1 Fig. 1 NPs shapes generated by polyhedra following a bottom-up approach: (a) tetrahedron; (b) double-tetrahedron; (c) pentagonal decahedron; (d) Ino's decahedron; (e) Marks's decahedron; (f) octahedron; (g) truncated octahedron; (h) cuboctahedron; (i) icosahedron; (j) truncated icosahedron; (k) icosidodecahedron.

Fig. 2

 2 Fig. 2 (a) Stacking implemented in polyhedra; (b) Examples of morphologies from those bulks, here for a hcp bulk; (c) helical nanorod and fcc cross-like nanoparticle; (d) hcp stacking with randomly distributed defective fcc areas.

Fig. 3

 3 Fig.3Comparison of calculated RDF profiles G calc (r) for 2.5 nm metal NPs and a metal-metal distance R of 2.7 Å. Black (plain line): 666-atoms spherical hcp NP; red (dashed line): 428-atoms decahedron (fcc). The RDF profiles are displayed up to 1.6 nm only, for the sake of clarity.

3

  Electronic structure of NPs and relation with their surface chemistry 3.1 d-band center model, Sabatier principle, volcano plots and the BEP relationship Introduction.

Fig. 4

 4 Fig. 4 Initial and final geometries and RDF profiles (in red) for a RMC procedure starting from (a) a cuboctahedron Ru 147 -CB or (b) an icosahedron Ru 147 -IC. In each case, the reference RDF profile is plotted in black. (c) Reference Ru 147 hcp geometry.

Fig. 5 ε

 5 Fig. 5 ε d mapping for (a) Ru 55 -MD (b) Co 55 -MD.

Fig. 6

 6 Fig. 6 Electronic structure of the Ru 4 H 4 (C 6 H 6 ) 4 tetrahedral cluster. (a) MO diagram calculated with Gaussian09. The energy levels for the [Ru 4 ] skeletal MOs are plot in green; (b) pDOS(ε) and pCOHP(ε) calculated from VASP.

Fig. 7

 7 Fig. 7 Adsorption of dppb on a bare Ru 55 cluster and on an hydrogenated Ru 55 H 64 cluster. Adsorption energies, given in kcal.mol -1 , are calculated as E ads (dppb) = E(dppb * ) -E(Ru 55 ) -E(dppb[anti]).

Fig. 8 d

 8 Fig. 8 d-band center coordination map for (a) a bare Ru 55 -hcp nanocluster; (b) the metallic core of Ru 55 (SnMe) 9 (the nine SnMe ligands are shown in white); (c) the metallic core of Ru 55 (CO) 66 (compound shown in Figure13).

Fig. 9

 9 Fig. 9 pDOS(ε) and COHP(ε) profiles for the Ru 55 H 64 (dppb) isomer shown in the inset.

  d σ , d π or d δ orbitals, in line with the ligand-field theory. The correlation between adsorption energies and the averaged d-band center value is only 0.79, whereas the one reported in figure 10a is 0.91. The difference lies in the fact that the error bars plot on the graph come from the decomposition of ε d values into d σ , d π and d δ contributions, allowing a refinement of the d-band center value according to the most important contribution in the CO coordination for each site.

Fig. 10

 10 Fig. 10 Correlation between adsorption energies of H atoms (in black) and CO ligands (in violet) on Ru 55 -MD and (a) ε d ; (b) CN.The points correspond to the individual values of adsorption energies whereas the line is obtained by linear regression from these points. r, the regression coefficient, is given in each case. The colors of the ε d and CN 3D activity maps cannot be directly compared owing to the arbitrariness of the scales. In each case, the Ru(0001) surface is chosen as the reference value, in white.

Fig. 11

 11 Fig. 11 Ru 288 -(µ 3 -SnBu) 31 H 62 , an hydrogenated ruthenium nanoparticle stabilized by tin ligands (0.2SnBu:Ru surface and 0.4H:Ru surface ). See Bonnefille et al. 4 for more details.

Fig. 12

 12 Fig. 12 Variations of the Gibbs free energy of adsorption ∆ a G (in kcal.mol -1 .nm -2 , eq. 10) on the Ru(0001) surface (dashed line, see our previous work 104 ) and on a spherical hcp Ru 55 NP (solid lines) as a function of the temperature at different coverage (all optimized geometries were first built with dressNPs). The H 2 chemical potential, µ(H 2 ,T, p) is related to temperature for p H 2 = 3 atm. The Ru 55 H n clusters corresponding to the most stable phases are shown below (1.6 ML: n = 70; 1.2 ML: n = 53; 0.8 ML: n = 35; 0.4 ML: n = 17).

Fig. 13

 13 Fig. 13 Variations of the Gibbs free energy of adsorption ∆ a G (in kcal.mol -1 .nm -2 , eq. 10) on a spherical hcp Ru 55 NP (a) as a function of the temperature at different coverage (all optimized geometries were first built with dressNPs). The CO chemical potential, µ(CO,T, p) is related to temperature for p CO = 1 atm; (b) as a function of pressure at T = 400 K. The most stable Ru 55 (CO) n clusters in each phase are shown below (1.53 ML: n = 66; 1.28 ML: n = 55; 1.00 ML: n = 43; 0.77 ML: n = 33; 0.51 ML: n = 22; 0.26 ML: n = 11).

Fig. 14

 14 Fig. 14 Phase diagram of the most stable H n CO m *Ru 55 structures in equilibrium with H 2 and CO as a function of T (in K) and p CO (in bar) for p H 2 = 1 bar.

  pDOS and pCOHP calculations. pbeVASPfit basis set of Lobster. Ru: {4p, 4d, 5s, 5p}; H: {1s}; C, O, P: {ns, np}. At least 12n + m + 8k bands are calculated for a Ru 55 H m (X) k . The charge spilling, a criterion that assesses the quality of the projection, is systematically lower than 0.7%. RMC process. ω E = 1, ω χ 2 = 10000, T = 1000K, number of MC steps = 10000 (although the optimal solution is usually found within 1000 steps). Sutton-Chen parameters for Ru taken from Popoola's work. 121 ab initio thermodynamics. H • T (H 2 , g), S • T (H 2 ,g), H • T (CO, g) and S • T (CO,g) are taken from the JANAF tables 122 (for example, H • 298.15 (H 2 , g) = 8467 J.mol -1 , S • 298.15 (H 2 , g) = 130.68 J.mol -1 .K -1 , H • 298.15 (CO), g = 8671 J.mol -1 and S • 298.15 (CO), g =

  Ruand Co respectively); kinetic energy cutoff: 500 eV; 104 Γ-centered calculations;119 Gaussian smearing of 0.02 eV width; geometry optimization threshold: residual forces on any direction less than 0.02 eV/Å; supercell size set to ensure a vacuum space of ca. 16 Å between periodic images of metal clusters (for example, 25×25×26 Å for Ru 55 ). Software: VASP; surfaces modeled by a periodic six-layer metal slab, six being kept frozen; ligands were adsorbed on one side of the slab; computational details are the same as in ref.104 , i.e. (2×2) cell; geometry optimizations: (7×7×1) Γ-centered k-points grid with a Gaussian smearing of 0.2 eV and 0.02eV/Å forces threshold; energy calculations on optimized geometries: (10×10×1) k-points grid and tetrahedron method proposed by Blöchl 120 for the k-space integration; kinetic energy cutoff: 500 eV. Coordination and averaged d-band centers were calculated with our homemade tools4vasp suite of utilities, which uses the DOS projected on a local basis set by the Lobster package (pDOS and pCOHP calculations, vide infra). Coordination d-band centers εd (µ k ) are calculated as:

	Periodic DFT calculations of surfaces.

114,115

; spin polarized DFT; exchange-correlation potential approximated by the generalized gradient approach proposed by Perdew, Burke, and Ernzerhof (PBE); 116 projector augmented waves (PAW) full-potential reconstruction; 117,118 PAW data sets for metal atoms treating the (n-1)p, (n-1)d and ns states (14 and 15 valence electrons for Calculation of coordination d-band centers.

| 1-20

| 1-20

[START_REF]From theory to application[END_REF][START_REF] Gerber | Nanocatalysis with Transition Metal Particles: Where are we Now[END_REF][START_REF] Carchini | [END_REF][4][5][START_REF] Philippot | Organometallic Ruthenium Nanoparticles and Catalysis[END_REF][START_REF] Amiens | [END_REF][8][9][10][11][12][13][START_REF] Neder | Diffuse Scattering and Defect Structure Simulations. A cook book using the program DIS-CUS[END_REF][START_REF] Hammer | [END_REF][16][17][18][START_REF] Balandin | Advances in Catalysis[END_REF][START_REF] Del Rosal | [END_REF]

| 1-20

| 1-20

Acknowledgments

We thank Drs K. Philippot and B. Chaudret for all our fruitful discussions these ten last years. We acknowledge financial support from the ANR-DFG (MOCANANO 2011-INTB-1011-1) project as well as the HPCs CALcul en MIdi-Pyrénées (CALMIP-HYPERION and CALMIP-EOS, grant P0611) and the Grand Equipement National de Calcul Intensif (GENCI-TGCC, grant 6211) for generous allocations of computer time.