PhD Romain Vayron

Hieu Vu

PhD Nguyen

Benoît Lecuelle

PhD Guillaume Haiat
email: guillaume.haiat@univ-paris-est.fr

Guillaume Haïat

Evaluation of dental implant stability in bone phantoms: comparison between a quantitative ultrasound technique and resonance frequency analysis Running Head (abbreviated

Keywords: Quantitative ultrasound, Resonance frequency analysis, Dental implant stability, Bone phantom, Bone quality

Background

Resonance frequency analyses and quantitative ultrasound methods have been suggested to assess dental implant primary stability.

Purpose

The purpose of this study was to compare the results obtained using these two techniques applied to the same dental implants inserted in various bone phantoms.

Materials and Methods

Different values of trabecular bone density and cortical thickness were considered to assess the effect of bone quality on the respective indicators (UI and ISQ). The effect of the implant insertion depth and of the final drill diameter was also investigated.

Results

ISQ values increase and UI values decrease as a function of trabecular density, cortical thickness and the screwing of the implant. When the implant diameter varies, the UI values are significantly different for all final drill diameters (except for two), while the ISQ values are similar for all final drill diameters lower than 3.2 mm and higher than 3.3 mm. The error on the estimation of parameters with the QUS device is between 4 and 8 times lower compared to that made with the RFA technique.

Conclusions

The results show that ultrasound technique provides a better estimation of different parameters related to the implant stability compared to the RFA technique.

Introduction

Dental implant stability, which is determinant for the surgical success [START_REF] Haiat | Effects of biomechanical properties of the bone-implant interface on dental implant stability: from in silico approaches to the patient's mouth[END_REF] , is determined by the quantity and biomechanical quality of bone tissue around the implant [START_REF] Franchi | Influence of different implant surfaces on peri-implant osteogenesis: Histomorphometric analysis in sheep[END_REF] . Two kinds of implant stability may be distinguished. The primary stability occurs at the moment of implant surgical insertion within bone tissue. Dental implant primary stability should be sufficiently high in order to avoid micromotion at the bone-implant interface after surgery, but should not be too high to avoid bone necrosis due to overloading of bone tissue [START_REF] Mathieu | Biomechanical determinants of the stability of dental implants: influence of the bone-implant interface properties[END_REF] . Secondary stability is obtained through osseointegration phenomena, a complex phenomenon of a multi-time and multiscale nature [START_REF] Mathieu | Biomechanical determinants of the stability of dental implants: influence of the bone-implant interface properties[END_REF] , which strongly depends on the implant primary stability.

Dental implant stability remains difficult to be assessed clinically because it depends on the implant properties (geometry, surface properties...), on the patient bone quality, as well as on the surgical protocol. In particular, there is a lack of standardization of the surgical procedures used in dental implantology, for example in the choice of the duration between implant insertion and loading, which may vary from 0 up to 6 months [START_REF] Raghavendra | Early wound healing around endosseous implants: a review of the literature[END_REF] . A compromise should be found in a patient specific manner between i) an early (or even immediate [START_REF] Buser | Long-term stability of osseointegrated implants in augmented bone: A 5-year prospective study in partially edentulous patients[END_REF] ) implant loading in order to stimulate osseointegration phenomena and ii) a late implant loading in order to avoid degradation of the consolidating bone-implant interface in early postsurgical stages [START_REF] Serra | Sequential bone healing of immediately loaded mini-implants[END_REF] .

Meanwhile, shortening the time of implant loading has become a challenge in recent implant developments to i) minimize the time of social disfigurement and ii) avoid gum loss. As a consequence, accurate measurements of implant biomechanical stability are of interest since they could be used to improve the surgical strategy by adapting the choice of the healing period in a patient-specific manner.

Assessing the implant stability is a difficult multiscale problem because of the complex heterogeneous nature of periprosthetic bone tissue and to remodeling phenomena [START_REF] Frost | Bone's mechanostat: A 2003 update[END_REF][START_REF] Wolff | The law of bone remodeling[END_REF] . Different approaches have been used to assess the implant stability in vivo. So far, most surgeons still rely on their proprioception because it remains difficult to monitor bone healing in vivo [START_REF] Serra | Sequential bone healing of immediately loaded mini-implants[END_REF] .

Accurate quantitative methods capable of assessing implant stability are required to guide the surgeons and eventually reduce the risk of implant failure. Magnetic resonance imaging [START_REF] Hecht | Magnetic resonance imaging susceptibility artifacts due to metallic foreign bodies[END_REF] and X-ray based [START_REF] Albrektsson | Osseointegrated oral implants. A Swedish multicenter study of 8139 consecutively inserted Nobelpharma implants[END_REF] techniques remain of limited interest to measure implant stability because of diffraction phenomena occurring at the bone-implant interface due to the presence of metal. Therefore, biomechanical methods have been developed, their main advantage consisting in the absence of ionizing radiation, inexpensiveness, portability and noninvasiveness. The measurement of the insertion torque to assess dental implant primary stability has been evoked but such approach remains limited 11 because the result is not only related to the properties of the bone-implant interface and because it cannot be used for secondary stability assessment. The Periotest (Bensheim, Germany) is a percussion test methods [START_REF] Schulte | Periotest-A new measurement process for periodontal function[END_REF][START_REF] Van Scotter | The Periotest method for determining implant success[END_REF] . Its sensitivity to striking height and handpiece angulation complicates the clinical examination [START_REF] Meredith | Relationship between contact time measurements and PTV values when using the Periotest to measure implant stability[END_REF] and limits the reproducibility of the measurements. The most commonly used biomechanical technique is the resonance frequency analysis (RFA) [START_REF] Valderrama | Evaluation of two different resonance frequency devices to detect implant stability: A clinical trial[END_REF] , which consists in measuring [START_REF] Meredith | Quantitative determination of the stability of the implanttissue interface using resonance frequency analysis[END_REF] the first bending resonance frequency of a small rod attached to the implant. The RFA technique allows to assess the implant anchorage depth into bone [START_REF] Meredith | The application of resonance frequency measurements to study the stability of titanium implants during healing in the rabbit tibia[END_REF] , marginal bone level [START_REF] Friberg | One-year prospective three-center study comparing the outcome of a ''soft bone implant'' (prototypeMkIV) and the standard Branemark implant[END_REF] and the stiffness of the bone-implant structure [START_REF] Ersanli | Resonance frequency analysis of one-stage dental implant stability during the osseointegration period[END_REF][START_REF] Pattijn | Resonance frequency analysis of implants in the guinea pig model: Influence of boundary conditions and orientation of the transducer[END_REF] . However, RFA cannot be used to identify directly the bone-implant interface characteristics [START_REF] Aparicio | Validity and clinical significance of biomechanical testing of implant/bone interface[END_REF] . No correlation between the implant stability quotient (ISQ) and bone implant contact (BIC) nor between ISQ and cortical thickness has been evidenced so far [START_REF] Seong | Correlations between physical properties of jawbone and dental implant initial stability[END_REF] . Finite element numerical simulation tools showed that the orientation and fixation of the transducer have an important effect on ISQ values 20 obtained with the older Osstell version with an L-shaped, wired transducer.

An alternative method has been developed by our group and consists in using a quantitative ultrasound (QUS) method [START_REF] De Almeida | Proposal for an ultrasonic tool to monitor the osseointegration of dental implants[END_REF] to investigate the properties of the bone-implant interface. The principle of the measurement relies on the dependence of ultrasonic propagation within the implant on the boundary conditions given by the biomechanical properties of the bone-implant interface [START_REF] Mathieu | Ultrasonic evaluation of dental implant biomechanical stability: an in vitro study[END_REF] . An in vitro preliminary study was carried out using a prototype titanium cylinder shaped implants inserted in bone tissue, showing the sensitivity of the ultrasound response of the implant to the quantity of bone in contact with the implant [START_REF] Mathieu | Ultrasonic evaluation of dental implant biomechanical stability: an in vitro study[END_REF] . The principle of the measurement was then validated experimentally by showing the sensitivity of the echographic response of a planar bone-implant interface to healing time using coin-shaped implant models [START_REF] Mathieu | Influence of healing time on the ultrasonic response of the bone-implant interface[END_REF] . Significant variations of the ultrasonic response of dental implants embedded in a bone substitute biomaterial (a tricalcium silicate-based cement) was shown to occur when the implants are subjected to fatigue loading [START_REF] Vayron | Variation of the ultrasonic response of a dental implant embedded in tricalcium silicate-based cement under cyclic loading[END_REF] . More recently, another in vitro study proved the potentiality of QUS methods to assess dental implant primary stability [START_REF] Vayron | Assessment of the in vitro dental implant primary stability using an ultrasonic method[END_REF] . A preclinical validation of the device in rabbits was carried out [START_REF] Mathieu | Influence of healing time on the ultrasonic response of the bone-implant interface[END_REF][START_REF] Vayron | Ultrasonic evaluation of dental implant osseointegration[END_REF] and showed that i) the measurement was sensitive to healing time and ii) a significant correlation of the measurement with the boneimplant contact (BIC) ratio measured with histology. Moreover, finite difference [START_REF] Mathieu | Numerical simulation of ultrasonic wave propagation for the evaluation of dental implant biomechanical stability[END_REF] and finite element numerical simulations [START_REF] Vayron | Finite element simulation of ultrasonic wave propagation in a dental implant for biomechanical stability assessment[END_REF][START_REF] Vayron | Assessment of dental implant biomechanical stability using quantitative ultrasound: a 3-D finite element study[END_REF] were carried out to understand the interaction between an ultrasonic wave and the bone-implant system, leading to a better performance of the device.

The comparison of the experimental results obtained in a controlled configuration using the QUS and RFA techniques would be of interest in order to better assess the performance of the different approaches. The aim of the present study is to compare these the RFA and QUS methods which have been evoked to assess dental implant stability. To do so, our strategy consists in using dental implants inserted in bone mimicking phantoms made of polyurethane because it allows to work under standardized and reproducible conditions. Different parameters related to the implant stability (such as the density of the bone phantom, the thickness of cortical bone, the insertion depth and the drill diameter) were investigated and the related variations of different quantities such as i) the RFA response and ii) the ultrasonic response were investigated.

II) Materials and Methods

II)1 Bone mimicking phantoms and dental implant

The bone mimicking phantoms used in the present study were made of rigid polyurethane foam (Orthobones®, 3B Scientific, Hamburg, Germany) with different values of bone density and of cortical thickness (1 and 2 mm). Cortical bone was modeled by the material type #40 PCF with a mass density equal to 0.55 g/cm³. Three types of trabecular mimicking phantoms were considered (#10, #20, #30 PCF) with mass density values equal to 0.16 g/cm³, 0.32 g/cm³ and 0.48 g/cm³, respectively, according to the manufacturer.

All dental implants used in this study were commercial implants manufactured by Zimmer Biomet (Winterthur, Switzerland) under the reference TSVT4B10. The geometrical characteristics are a length of 10 mm and an external diameter of 4.1 mm.

A conical cavity was then created in the block bone test, using color-coded, 10-mm-length surgical drills preconized by Zimmer Biomet and used as a reference protocol by dentists. Different values were considered for the final drill diameter, as described in details in subsection II)3.

II)2 Measurement methods of dental implant stability

Two different and complementary methods were used in order to measure the implant stability of each implant in each configuration, which are described in subsection II)3.

a) Resonance frequency analysis

For each implant and each configuration, the RFA response of the implant was measured in ISQ units (on a scale from 1 to 100) using the Osstell device (Osstell, Göteborg, Sweden).

Figure 1a shows the configuration of the measurements which were realized using a smart peg placed in the implant, as recommended by the manufacturer. Each measurement was performed 5 times (in order to assess the reproducibility of the measurements) in two perpendicular directions denoted 0° and 90° relatively to an arbitrary axis chosen for each block. The average and standard deviation values ISQm and ISQstd of the ten different ISQ values were determined for each implant and each configuration.

b) Quantitative ultrasound device

Figure 1b shows the dedicated QUS device, which consists in a 5 mm diameter planar ultrasonic contact transducer (Sonaxis, Besançon, France) generating a 10 MHz broadband (the frequency bandwidth is approximately equal to 6-14 MHz) ultrasonic pulse propagating perpendicularly to its active surface (monoelement transducer used in echographic mode). The probe was attached rigidly to a titanium alloy dental healing abutment which can be screwed into the implant so that the measurements are not influenced by positioning problems of the probe relatively to the healing abutment, as it was the case in previous studies [START_REF] Vayron | Variation of the ultrasonic response of a dental implant embedded in tricalcium silicate-based cement under cyclic loading[END_REF] . The QUS device was screwed into the implant in order to realize each measurement for each implant and each configuration, as shown in Fig. 1b. The ultrasonic probe was linked to a pulser-receiver via a standard coaxial cable. A transient recorder was used to record the radiofrequency (rf) signal with a sampling frequency equal to 100 MHz.

For each measurement, the transducer was screwed in the implant with a torque of 0.035N.m, which is around 10 times lower than torque values recommended by implant manufacturers [START_REF] Kanawati | Measurement of clinicians' ability to hand torque dental implant components[END_REF] and the echographic measurement was made instantaneously. The transducer was then unscrewed and the same measurement was carried ten times in order to assess the reproducibility of the measurements.

The same signal processing technique as the one used in [START_REF] Vayron | Variation of the ultrasonic response of a dental implant embedded in tricalcium silicate-based cement under cyclic loading[END_REF] was used to derive a quantitative ultrasonic indicator UI which had been shown to be related with dental implant stability. The indicator UI was devised to quantitatively estimate the average amplitude of the signal between 10 and 120 µs. Therefore, the envelop S(t) of the rf signal s(t) was determined and the indicator UI is given by:

𝑈𝐼 = ∑ 𝑆(𝑖𝑇 0 ) 12000 𝑖=1000 , (1)
where T0=0.01 µs corresponds to the sampling period. The time window chosen to determine the indicator UI was chosen as follows. The upper bound of the time window equal to 120 µs was chosen, which corresponds to a compromise between a sufficient duration to obtain pertinent information and the requirement of a sufficient signal to noise ratio for all rf signals.

Moreover, the time window used to compute the indicator UI starts at a time of 10 µs because the amplitude of the envelop of the rf signals before 10 µs is approximately constant due to a saturation of the amplitude.

The average and standard deviation values UIm and UIstd of the indicator UI obtained for the ten measurements were determined for each implant and each configuration.

II)3 Experimental protocol

The effect of varying the different parameters described below on the variation of the ISQ and the indicator UI was investigated in this study, as summarized in Figure 2. For each configuration described below, the values of the averaged and standard deviation values ISQm and ISQstd (respectively UIm and UIstd) of the ISQ (respectively the ultrasonic indicator UI) were measured.

a) Variation of trabecular bone density

The effect of modifying trabecular bone properties was first investigated by comparing the results obtained with the two methods described in subsection II)2 for implants fully inserted in different bone mimicking phantoms having the same cortical thickness equal to 2 mm, which corresponds approximately to the clinical situation [START_REF] Tabassum | Influence of surgical technique and surface roughness on the primary stability of an implant in artificial bone with different cortical thickness: a laboratory study[END_REF] . Three different materials were used to mimic trabecular bone (#10, #20 and #30 PCF). Three implants were inserted with a final drill of 3.4 mm of diameter in each type of bone mimicking phantom, resulting in a total of 9 configurations.

b) Variation of cortical bone thickness

The effect of changing the cortical bone thickness was then investigated because cortical thickness is known to be an important determinant of implant stability 11 . To do so, two bone mimicking phantoms with two different cortical thicknesses (1 and 2 mm) and the same material to mimic trabecular bone (#10 PCF) were considered. Again, three implants were fully inserted with a final drill of 3.4 mm of diameter in each type of bone mimicking phantom, resulting in a total of 6 configurations.

c) Variation of the diameter of the final drill

The effect of the diameter of the final drill used to machine the cavity in two given bone mimicking phantom was investigated. The two phantoms used had the same cortical thickness (1 mm) and trabecular bone was mimicked by #10 and #30 PCF. Ten values of diameter of the final drill comprised between 2.7 and 3.6 mm with a step of 0.1 mm were considered. Three implants were considered for each diameter tested, leading to a total number equal to 30 configurations.

d) Variation of the implant insertion depth

Eventually, the effect of the implant insertion depth was investigated because it corresponds to a simple way of imposing a variation of dental implant stability, similarly as what has been done in [START_REF] Vayron | Ultrasonic evaluation of dental implant osseointegration[END_REF] . To do so, three cavities were realized in a given bone mimicking phantom (cortical thickness of 1 mm and #30 PCF for trabecular bone) with a final drill diameter of 3.4 mm.

Then, an implant was screwed in each cavity so that half of the implant was inserted into the cavity. The average and standard deviation values of the ISQ and of the ultrasonic indicator UI were measured. The implant was then screwed by π rad in order to increase the implant surface area in contact with bone mimicking phantom. The two measurements described above (ISQ and UI) were realized again. The procedure was repeated until the implant was completely inserted in bone mimicking phantom, this screwing level is denoted "0". Then, the implant was then unscrewed by π rad in order to reduce the surface area of the implant in contact with bone block. The same two measurements (ISQ and UI) were again carried out. The procedure was reproduced until the implant was completely detached from the test block.

II)4 Statistical analyses

A one-way analysis of variance (ANOVA) and Tukey-Kramer tests were performed to evaluate the significance of variations of the ultrasonic indicator UI and the ISQ as a function of the different configurations. Statistical differences were defined at a 95% confidence level. However, the results obtained with RFA measurements in #20 PCF and #30 PCF are not significantly different (p-value = 0.93).

III) Results

III)1 Effect of bone density

III)2 Effect of cortical thickness

5.

ANOVA shows a significant effect of final drilling diameter on the values of the UI (p < 10 -12 )

and on the values of ISQ (p < 10 -3 ). The stars in Fig. 5 indicate the configurations for which the results obtained for the UI and for the ISQ are statistically similar, which was obtained using a Tukey-Kramer test. The results obtained with Tukey-Kramer test indicate that the UI obtained for all final drill diameters are significantly different to the other ones, except for the results obtained for the final drill diameters equal to 2.8 and 2.9. However, the results obtained for ISQ with all drill diameters lower than 3.3 (respectively higher that 3.2) are statistically similar.

III) 4 Effect of the insertion depth

Figure 6 shows the average and standard deviations of the ISQ and of the UI for the same implant with different values of the penetration depth. The arrows in Fig. 5 indicate the evolution of the results when the implant is screwed (from a bone level of -5π rad to 0 rad) and then unscrewed (from a bone level of 0 rad to -4π rad). The UI value is shown to decrease during the insertion of the dental implant in the bone block while the UI value increases when the implant in unscrewed from the test block. The ISQ values increase during the insertion stages and then decrease when the implant is unscrewed. In order to clarify the figure, we have chosen to present only one curve (instead of three) because the results obtained using the other test block are qualitatively similar to those shown in Fig. 6.

IV) Discussion

The originality of the present study was to compare two different techniques (RFA and ultrasound measurements) in order to investigate the primary stability of dental implants inserted in artificial bone blocks (Orthobones®). Different stability conditions were considered by varying different parameters such as the type of bone block (trabecular bone density, cortical thickness), the final diameter drill and the insertion depth, in order to simulate different situations mimicking variations of dental implant primary stability.

IV) 1 Physical interpretation and comparison with the literature

The values of the ISQ are shown to increase when i) trabecular density increases (see Fig. 3),

ii) cortical thickness increases (see Fig. 4), and iii) the insertion depth increases (see Fig. 6).

These results are in agreement with previous studies [START_REF] Bayarchimeg | Evaluation of the correlation between insertion torque and primary stability of dental implants using a block bone test[END_REF][START_REF] Meredith | Relationship between contact time measurements and PTV values when using the Periotest to measure implant stability[END_REF][START_REF] Meredith | The application of resonance frequency measurements to study the stability of titanium implants during healing in the rabbit tibia[END_REF][START_REF] Andreotti | Relationship Between Implant Stability Measurements Obtained by Two Different Devices: A Systematic Review[END_REF] and can be explained by the higher rigidity of the block-implant system induced by an increase of trabecular density, of cortical thickness and of implant insertion depth.

The values of the UI are shown to decrease when i) trabecular density increases (see Fig. 3), ii) cortical thickness increases (see Fig. 4), iii) the final diameter drill decreases (see Fig. 5) and iv) the insertion depth increases (see Fig. 6). These results are in agreement with previous papers by our group that showed that i) the UI also increases when a dental implant is unscrewed in a bone sample in ex vivo [START_REF] Vayron | Assessment of the in vitro dental implant primary stability using an ultrasonic method[END_REF] as well as in silico [START_REF] Vayron | Assessment of dental implant biomechanical stability using quantitative ultrasound: a 3-D finite element study[END_REF] studies, ii) the UI decreases when the Bone-Implant contact ratio increases for in vivo measurements [START_REF] Vayron | Ultrasonic evaluation of dental implant osseointegration[END_REF] and iii) the UI decreases when trabecular bone quality increases and when cortical thickness increases for in silico measurements [START_REF] Mathieu | Numerical simulation of ultrasonic wave propagation for the evaluation of dental implant biomechanical stability[END_REF][START_REF] Vayron | Finite element simulation of ultrasonic wave propagation in a dental implant for biomechanical stability assessment[END_REF][START_REF] Vayron | Assessment of dental implant biomechanical stability using quantitative ultrasound: a 3-D finite element study[END_REF] . The aforementioned results can be explained by the fact that the UI is related to the amplitude of the echographic response of the implant, which depends on the boundary conditions applied to the implant external surface. The implant acts as a wave guide in which the acoustical energy is trapped [START_REF] Mathieu | Numerical simulation of ultrasonic wave propagation for the evaluation of dental implant biomechanical stability[END_REF][START_REF] Vayron | Finite element simulation of ultrasonic wave propagation in a dental implant for biomechanical stability assessment[END_REF][START_REF] Vayron | Assessment of dental implant biomechanical stability using quantitative ultrasound: a 3-D finite element study[END_REF] . When bone quantity and quality increase around the implant surface, the gap of mechanical properties between the implant and its environment decreases, which induces an increase of the transmission coefficient. Therefore, the ultrasound energy decreases faster since ultrasound may propagate in the surrounding medium, resulting in a faster decrease of the implant echographic response. This physical interpretation may qualitatively explain the results corresponding to the variation of the UI obtained herein.

Figure 6 shows that the values of the UI (respectively ISQ) are lower (respectively higher) for the same insertion value during the screwing phase compared to the unscrewing phase. These results may be explained by damage and wear phenomena occurring during the implant insertion, which lead to a higher bone-implant contact ratio during the screwing phase compared to the unscrewing phase, thus explaining the lower (respectively higher) values of UI (respectively ISQ) during the insertion phase compared to the unscrewing phase.

The results obtained in the present study depend on the displacement of the implant surface generated by the ultrasound transducer, which is difficult to measure using an experimental approach because accessing the surface of an implant embedded in a bone mimicking phantom remains a difficult task. However, the analysis of the interaction between an ultrasonic wave and a dental implant in a situation similar to the one described herein was carried out in previous in silico studies by our group [START_REF] Mathieu | Numerical simulation of ultrasonic wave propagation for the evaluation of dental implant biomechanical stability[END_REF][START_REF] Vayron | Finite element simulation of ultrasonic wave propagation in a dental implant for biomechanical stability assessment[END_REF][START_REF] Vayron | Assessment of dental implant biomechanical stability using quantitative ultrasound: a 3-D finite element study[END_REF] . Briefly, the importance of the lateral wave was evidenced in preliminary studies considering simple cylindrical shaped implants [START_REF] Mathieu | Numerical simulation of ultrasonic wave propagation for the evaluation of dental implant biomechanical stability[END_REF][START_REF] Vayron | Finite element simulation of ultrasonic wave propagation in a dental implant for biomechanical stability assessment[END_REF] . In a more recent study considering the geometry of an actual dental implant [START_REF] Vayron | Assessment of dental implant biomechanical stability using quantitative ultrasound: a 3-D finite element study[END_REF] , the sensitivity of the ultrasonic indicator to changes occurring around 30 µm around the implant surface was evidenced.

IV) 2 Comparison between the RFA and QUS methods

The results obtained in this study allow to compare quantitatively the ultrasound and RFA approaches to assess dental implant primary stability in controlled configurations.

a) Drill diameter

As shown in Fig. 5, the values obtained for the UI are significantly different when all final drill diameters vary between 2.7 mm and 3.6 mm (except for 2.8 and 2.9 mm) whereas the values obtained for the ISQ are similar when the final drill diameter varies between 2.7 mm and 3.2 and when it varies between 3.3 and 3.6. This results show that the ultrasound device can discriminate implants inserted in test blocks with a wider range of final diameter size compared to the RFA measurements.

b) Sensitivity analysis

The results can be analyzed in order to determine the sensitivity of QUS and RFA methods to the different parameters related to the implant stability. To do so, a simple two step method is described in what follows (refer to [START_REF] Vayron | Assessment of the in vitro dental implant primary stability using an ultrasonic method[END_REF] for further details on this method) and constitutes a simple approach in order to estimate the sensitivity of both methods to a variation of a given parameter X related to the implant primary stability. Note that X may correspond to trabecular bone density, to cortical thickness or to the implant insertion depth. The sensitivity of the ISQ to variations of the final drill diameters was not investigated herein because the results obtained with the RFA techniques are weakly sensitive to the drill diameter (see subsection a) above).

The first step is to perform a linear regression analysis of the average value of the ISQ (respectively of the UI) as a function of the parameter X by analyzing the results shown in Fig. 3456, which leads to the following approximated relations:

𝐼𝑆𝑄 ̃= 𝑎 𝐼𝑆𝑄 𝑋 + 𝑏 𝐼𝑆𝑄 (2)
𝑈𝐼 ̃= 𝑎 𝑈𝐼 𝑋 + 𝑏 𝑈𝐼 , ( 3
)
where X is the investigated parameter, 𝐼𝑆𝑄 ̃ and 𝑈𝐼 ̃ are the approximated values of the ISQ and the UI respectively and 𝑎 𝐼𝑆𝑄 and 𝑏 𝐼𝑆𝑄 (respectively 𝑎 𝑈𝐼 and 𝑏 𝑈𝐼 ) are the coefficients found by applying a linear regression analysis to the variation of the ISQ (respectively the UI) as a function of X.

The second step of the method consists in using the averaged reproducibility error corresponding to a given configuration in combination with the linear regression analyses corresponding to Eqs. 2 and 3, in order to assess the error realized on the estimation of the parameter X, noted in what follows XISQ (respectively XUI) for the error realized using the RFA (respectively ultrasound) technique:

XISQ = ISQstd / aISQ (4) XUI = UIstd / aUI (5)
As expected, the error on the estimation of the parameter X increases when the reproducibility error (given by ISQstd and UIstd) increases and when the sensitivity of the method (given by 𝑎 𝐼𝑆𝑄 and 𝑎 𝑈𝐼 ) decreases.

Table 1 shows the values of the error made on the determination of the trabecular density, of the cortical thickness and of the insertion depth obtained by applying the aforementioned procedure on the data shown in Fig. 3, 4 and 6 respectively. As shown in Table 1, the error made on the estimation of the trabecular density using the ultrasound method is more than four times lower compared to that using the RFA analysis, which was obtained by analyzing the results obtained in Fig. 3. Note that the values of ISQ obtained for #20 and #30 PCF test blocks were not significantly different, while the results obtained for the UI for #20 and #30 PCF test blocks were significantly different. Therefore, the quantitative ultrasound method is more sensitive than the RFA method to retrieve the trabecular density.

Table 1 also shows the results corresponding to the error realized on the cortical thickness estimation using the data shown in Fig. 4. The results indicates that the error realized on the estimation of the cortical thickness using the ultrasound method is approximately eight times lower compared to that using the RFA analysis. Note that the values of ISQ obtained with cortical thicknesses equal to 1 and 2 mm were not significantly different while the results obtained for the UI obtained with cortical thicknesses equal to 1 and 2 mm were significantly different. Therefore, the quantitative ultrasound method is more sensitive the RFA analysis to retrieve the cortical thickness.

Table 1 shows the results corresponding to the error realized on the insertion level estimation using the data shown in Fig. 5. The results indicates that the error realized on the estimation of the insertion level using the ultrasound method is approximately four times lower compared to that using the RFA analysis. Therefore, the quantitative ultrasound method is more sensitive the RFA analysis to retrieve the cortical thickness.

As shown in Table 1, the performances of the quantitative ultrasound device to retrieve dental implant primary stability are better than the performances of RFA technique in the simple configurations considered herein. These results can be explained by the fact the ultrasound device is sensitive to bone properties at the intimate contact with the implant, which is precisely the important parameter determining the implant stability [START_REF] Haiat | Effects of biomechanical properties of the bone-implant interface on dental implant stability: from in silico approaches to the patient's mouth[END_REF][START_REF] Huja | Microdamage adjacent toendosseousimplants[END_REF][START_REF] Luo | The effect of surface roughness on the stress adaptation of trabecular architecture around a cylindrical implant[END_REF] . More specifically, it has been shown in a previous in silico study [START_REF] Vayron | Finite element simulation of ultrasonic wave propagation in a dental implant for biomechanical stability assessment[END_REF] that the ultrasound technique is sensitive to the bone properties at a distance lower than around 30 µm. In the case of the RFA technique, the frequency response of the implant is sensitive to the biomechanical properties of the boneimplant interface as well as to the surrounding medium since the entire test block is likely to be excited by the RFA device. Although the amplitude of the displacement field induced by the RFA method is small compared to the implant size, the excitation mode corresponds to the first resonance of the bending mode, which indicates that a vibration is induced by the transducer.

This vibration is applied to the bone-implant system as a whole (and not only to the implant interface) through the implant surface. However, in the case of the QUS device, the excitation induced by the transducer differs from the RFA device because propagation phenomena should be considered [START_REF] Vayron | Assessment of dental implant biomechanical stability using quantitative ultrasound: a 3-D finite element study[END_REF] , instead of resonance effects. The ultrasound wave propagating within the implant [START_REF] Vayron | Assessment of dental implant biomechanical stability using quantitative ultrasound: a 3-D finite element study[END_REF] is therefore sensitive to its interface only, which has been quantified in [START_REF] Vayron | Assessment of dental implant biomechanical stability using quantitative ultrasound: a 3-D finite element study[END_REF] .

Moreover, the results obtained with the RFA device depend on the direction of solicitation, which is not the case of the ultrasound device. In the case of the RFA method, the dependence of the results on the orientation of the transducer has been explained by the anisotropic nature of bone tissue but it remains difficult to extract useful information from such variation [START_REF] Simunek | Development of implant stability during early healing of immediately loaded implants[END_REF][START_REF] Sim | Factors influencing resonance frequency analysis assessed by Osstell mentor during implant tissue integration: I. Instrument positioning, bone structure, implant length[END_REF] .

IV) 3 Limitations of the study

The main limitation of this study is that bone mimicking phantoms were considered, which allows to work under standardized and reproducible conditions. Note that both techniques had been validated previously independently in vivo [START_REF] Vayron | Ultrasonic evaluation of dental implant osseointegration[END_REF][START_REF] Andreotti | Relationship Between Implant Stability Measurements Obtained by Two Different Devices: A Systematic Review[END_REF] . However, the artificial bone blocks are adapted to mimic the healthy jaw bone because the mean bone mineral density was 0.31 g/cm³ for the posterior maxilla and 0.55 g/cm³ for the anterior maxilla [START_REF] Devlin | A comparison of maxillary and mandibular bone mineral densities[END_REF] . For healthy jaw bone, the mean cortical thickness for the mandible is 2.22±0.47 mm and the mean cortical bone thickness for the maxilla is 1.49±0.34 mm [START_REF] Tabassum | Influence of surgical technique and surface roughness on the primary stability of an implant in artificial bone with different cortical thickness: a laboratory study[END_REF] . Further in vivo validations should be carried out clinical studies should also be made.

Another limitation consists in the fact that only one implant type was used in this study, because the goal was to investigate the effect of variations of the implant stability on the UI and on the ISQ. It would be of interest to carry out the same study with other implant types. However, previous studies [START_REF] Vayron | Variation of the ultrasonic response of a dental implant embedded in tricalcium silicate-based cement under cyclic loading[END_REF][START_REF] Vayron | Ultrasonic evaluation of dental implant osseointegration[END_REF] have shown the validation of the ultrasound device using other implant types.

The present study does not address the important clinical question corresponding to the

Conclusion

This study allows to compare the results obtained with two different approaches aiming at estimating primary dental implant stability, which are realized with the same implants under various configurations. Namely, we considered a variation of the following parameters: i) trabecular bone density, ii) cortical thickness, final drill diameter and iv) penetration depth. All results are consistent and can be explained by physical analyses of the biomechanical phenomena occurring around the implant. Moreover, we found that ultrasound technique provides a better estimation of different parameters related to the implant stability compared to RFA techniques. Therefore, the present study paves the way for the development of an ultrasonic device to estimate dental implant stability that could be used in the clinic provided further in vivo investigations. 1. Results obtained for the linear regression analysis (a and b), the average standard deviation of each indicator and for the error realized on the estimation of each parameter X.

Figure 3

 3 Figure3shows the results obtained for the nine configurations corresponding to three models

Figure 4 III) 3 Figure 5

 435 Figure 4 shows the variation of the ISQ and of the UI for different implants inserted in test

  determination of the extent to which the RFA and QUS methods can actually accurately measure an important information from a clinical point of view. The question of the relationship between the accuracy of each method versus the precision needed in the context of implant stability remains open and is not in the scope of the present study. Further clinical studies are needed to determine whether either method (RFA or QUS) can accurately measure a clinically relevant information.

Figure 1 :

 1 Figure 1: Description of the experimental protocol. Measurement realized using a) the

Figure 2 :

 2 Figure 2: Summary of the experimental protocol describing the parameters investigated in this

Figure 3 :

 3 Figure 3: Variation of the values of the ISQ and of the UI for implants inserted in test blocks

Figure 4 .

 4 Figure 4. Variation of the values of the ISQ and of the UI for 3 implants inserted in a test block

Figure 5 .

 5 Figure 5. Variation of the values of the ISQ and of the UI for implants inserted in cavities

Figure 6 .

 6 Figure 6. Variation of the ultrasonic indicator and ISQ values for different screwing levels

Acknowledgements

This work was supported by the French National Research Agency (ANR) through the PRTS program (project OsseoWave n°ANR-13-PRTS-0015-02). This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 682001, project ERC Consolidator Grant 2015 BoneImplant).


Christophe Raynaud

Iker Del

Franck Jolibois

Laurent Maron

Romuald Poteau

Iker Del Rosal

Multicentered effective group potentials: ligand-field effects in organometallic clusters and dynamical study of chemical reactivity

Introduction

Organometallic nanoparticles (NPs) elicit a great interest due to their physical and chemical properties intermediate between small molecular compounds and the bulk material, resulting from surface or quantum size effects [START_REF] Schmid | Nanoparticles. From theory to application[END_REF] . Quantum chemists face a challenge, owing to the size of these species and the difficulty to determine their wavefunction. Three complementary approaches may bring different insights: (i) the bottom-up route, which consists in the accurate description of small organometallic clusters and in an extrapolation to NPs, (ii) the top-down route, which consists in the description of ligands adsorbed at the surface of slab models by means of periodic quantum chemistry methods, (iii) the actual description of NPs, which may be achieved by density functional theory (DFT) methods, but which are rather the scope of semi-empirical approaches [START_REF] Calvo | [END_REF] . Recently, structural and spectroscopic properties of scale models of ruthenium NPs have been investigated: several [Ru 4 ] and [Ru 6 ] organometallic clusters, previously synthesized or non-existent, have provided useful results in the context of organometallic ruthenium NPs 3 . We have more particularly focused on the 1 H NMR properties of these compounds and established a link between electron-deficiency and low-field resonance of hydrides adsorbed at the surface of Ru clusters or NPs. From a more general point of view, 1 H NMR and 2 H MAS-NMR techniques are of high practical importance for characterizing the presence of hydrides (or deuterides) in ruthenium clusters and nanoparticles [4][5][6] . According to the work done in our group in the specific context of Ru clusters and mononuclear complexes, DFT calculations of chemical shielding tensors (σ) as well as electric field gradient (EFG) tensors seem accurate 3,7,8 . However, one should pay attention that owing to the NMR timescale, 1 H and 2 H NMR properties are highly dependent of the motion of hydrides in organometallic complexes or at the surface of Ru NPs. EFG calculations in Ref. 8 were mainly compared to low-temperature MAS-NMR data which were considered close to the so-called rigid limit, but such experiments are not so easy to achieve 9 . At sufficiently high temperature, proton spectra of fluxional compounds consist in the superposition of several bonding situations. A relevant, but computationally expensive, way to take into account the motion of protons consists in performing molecular dynamics (MD) simulations at the ab initio or DFT level of theory and to average the theoretical spectroscopic data [10][11][12] . We have developed in our group a Car-Parrinello MD (CPMD) code using atom-centered basis functions, which is a tool of choice for taking into account the role of temperature on the motion of nuclei 13,14 , and to calculate several properties with quantum chemistry methods more familiar to theoretical chemists. The importance of the Berry pseudorotation process in the interpretation of the 19 F NMR spectrum of phosphorus pentafluoride (PF 5 ) has been shown in Ref. 11. Ab initio MD (AIMD) simulations have been performed [a] Université de Toulouse ; INSA, UPS ; LPCNO, IRSAMC ; 135 avenue de Rangueil, F-31077 Toulouse, France and CNRS ; UMR 5215 (IRSAMC) ; F-31077 Toulouse, France [b] Institut Charles Gerhardt, CNRS 5253, Université Montpellier 2, CC 1501 Place Eugène Bataillon, 34095 Montpellier, France e-mail: romuald.poteau@univ-tlse3.fr to generate a large number of geometries considered for NMR parameter computations at the DFT level. It has been shown that averaged 19 F chemical shifts and spin-spin coupling constants calculated for the five fluorine atoms converge towards the NMR equivalence of the five atoms when the Berry pseudorotation mechanism is accounted for. Although AIMD can also be applied to larger systems, there is a need for lowering the cost of quantum chemistry calculations. It is the purpose of effective group potentials (EGPs), which aim at simplifying molecular ab initio calculations for large systems involving bulky ligands as long as these ligands can be supposed to play the role of spectator groups. In an EGP, a group of atoms is replaced by one or several fictitious atoms associated with energy levels and orbitals necessary for representing the active/spectator bond(s), in conjunction with a mono-or multi-centered pseudopotential. The principle of molecular pseudopotentials was briefly reported by Durand and Malrieu in 1987 [START_REF] Durand | Advances in Chemical Physics : ab initio methods in quantum chemistry[END_REF] . Primary developments and promising tests were carried out in the Toulouse group in the late eighties, but were not published since the transferability of this new tool was not obvious. This was the prolegomena to EGPs in Toulouse. EGPs witnessed a renaissance in the late nineties and have been extensively developped and used in their present formulation. They are hereafter considered as the first generation of EGPs. Their definition and the way to derive their parameters have been in particular given in Ref. 16. The accuracy and limitations of this tool has been discussed in several papers, mainly devoted to organometallic compounds [17][18][19][20] and its usefulness in terms of solution to the link atom problem has also been considered 21 . A summary of the cons and pros of EGPs has temporarily ended this successful attempt to reduce the computational time by reducing both the number of electrons and nuclei explicitly treated in ab initio calculations 22 . The reader interested with an overview of the miscellaneous propositions of molecular pseudopotentials is suggested to read the introduction of the article of Slavíček and Martínez 23 , which introduces a multicentered valence electron effective potential description of functional groups which also succeeds in the context of excited electronic states. Note also the interesting proposition to modify carbon ECPs for substituted cyclopentadienide ring carbons for use in QM/MM calculations 24 .

The main purposes of this article is firstly to provide new directions for the derivation of multicentered EGPs and secondly to show that AIMD simulations can take advantage of the reduction of the cost involved by EGPs, owing to some limitations, which should be overcome in the next generation of these molecular pseudopotentials. New results for a tetrahedral ruthenium cluster, namely Ru 4 Ar # 4 H 4 (where Ar # stands for a pseudo η 6 -C 6 H 6 ligand), are also provided as well as a joint EGP/ab initio MD study aiming at finding free energies of activation in a bis-cyclopentadienyl lanthanide complex.

2 Effective Group Potentials 2.1 The case of η 6 -C 6 H 6 ligands
The construction of an EGP first requires to identify the characteristics of the target functional group from simple MO theory. An η 6 -C 6 H 6 ligand is considered as a 6 π-electrons system. Metal-ligand interactions will essentially involve π (shown in figure 1a), and possibly π * , MOs of benzene and the symmetry-relevant AOs of the metal. EGPs considered in this section will be extracted in that context. Although σ and π MOs are not well separated in energy (actually the highest occupied σ MOs lie above the lowest occupied π MOs), σ MOs are not considered due to a weak overlap with the metal AOs [START_REF] Albright | Orbital interactions in chemistry[END_REF] .

The basic methodology underlying the derivation of EGPs, rather simple, has been described in Ref. 16 and recalled in Refs. 21 and 22. Although it will not be explained in details in the present paper, some guidelines are necessary in order to present new ideas which could serve as a basis for deriving new multicentered molecular pseudopotentials. The role of EGPs is to remove some of the electrons and nuclei of functional groups which are considered as inactive, or spectators. This is done in two indivisible steps, achieved in the context of the Hartree-Fock (HF) method and based on reference calculations performed on benzene, in a large basis set {χp} of Gaussian-Type orbitals (GTOs, namely a polarized double-zeta basis set in the present case) :

1. The definition of a fictitious chemical group, with a reduced number of electrons and nuclei. An η 6 -C 6 H 6 ligand is substituted by an hexagonal system with six pseudo-carbon atoms C # which occupy the position of the actual carbon atoms in benzene. Each C # atom brings a single electron in a p atomic orbital (AO), in order to describe the three occupied π MOs φ i of benzene, which define the minimal set of active orbitals of this functional group (N act = 3, see figure 1a). We have already checked that this minimal truncated basis set {fq } is sufficient for a proper description of the ligand-metal bonding. The exponent of the truncated GTO basis set {fq } is done in order to obtain the best molecular valence pseudo-orbitals (MVPOs), that is pseudo-MOs φ ps i built in this truncated basis set which are as close as possible to the N act active molecular orbitals of the real system [START_REF] Poteau | [END_REF] . This step, analogous to the definition of valence pseudo-orbitals in the shape consistent version in the effective core potentials (ECPs) determination [START_REF] Durand | [END_REF]27 , involves de facto a σ/π separation. Note that the optimized exponent for benzene (ζp : 0.332) is very close to the exponent optimized 19 for η 5 -C 5 H 5 (ζp : 0.312), thus suggesting that a unique exponent could be considered for carbon atoms involved in aromatic compounds coordinated by their π electrons. The pitfall of the drastic simplification introduced in this first step lies in the fact that such calculations are done in the context of Self-Consistent-Field (SCF) methods and not extended-Hückel or tight-binding calculations. As a consequence, although the actual ligand and its C # 6 counterpart are expected to be isolobal analogues according to the definition given by Hoffmann 28 , the SCF process is ruled by the optimal stabilization of the C # 6 system. It is shown in figure 1b that the energies and shapes of the three actual HF MOs of C # 6 are far from the A 2u and Eg π MOs of benzene. The lowest one consists in a bonding symmetric σ MO built from radial p AOs. The second one is the expected totally bonding A 2u π MO, but it lies 18 eV above its C 6 H 6 counterpart. The third one is a σ bonding MO which exhibits two nodal planes. They are far from the target π MOs of an η 6 -benzene ligand, both in terms of shape and energy. It it thus mandatory to introduce an operator in the hamiltonian in order to obtain at the end of the SCF process the relevant active MOs of the functional group. If one defines F ref and F # as the Fock operators of the reference and model systems, respectively, the operator which should be added to F # is the exact FEGP operator. This operator is obviously molecule-dependent and not transferable to other systems. It is the reason why we hereafter introduce a transferable operator, namely ŴEGP (see refs. 16 or 21 for more details). 2. The definition of an operator associated to the fictitious chemical group, herafter called ŴEGP , which will compel the SCF process to converge on three π MOs which must be as close as possible to those of the real chemical group in terms of shape and energy. The main part of this molecular pseudopotential is formulated with a non-local form ( P P αnm|gn >< gm|), which presents some numerical advantages 29 and which ensures its transferability [START_REF] Poteau | [END_REF] :

ŴEGP = X n X m αnm|gn >< gm| - n # X i=1 Z eff i r + n # X i=1 Ŵps(i) (1)
where n # is the total number of fictitious atoms, the functions |gn > are gaussians functions with predetermined exponents, and Ŵps(i) is the ECP of the real atom. The advantage of this ECP is to already take into account, if necessary, the repulsion of core electrons, although repulsion could have been included as a criterion to be fulfilled by the non-local part of ŴEGP (in that case, ŴEGP would be reduced to the sum of the non-local part and the Coulomb term). Thus the calculation of the EGP operator in a GTO basis set is obtained by simple products of gaussian overlaps, in addition to the integrals of the Coulomb and ECP terms which are available in standard quantum chemistry packages. The exponents of the gn functions are chosen so that two conditions are simultaneously fulfilled: (i) the isolobal analogy between the active MOs of the real functional group and the SCF MOs of the fictitious group; (ii) ŴEGP must be a short-range operator since it must only act on the electrons of the spectator chemical group which are kept. Similarly to the cyclopentadienyl case, the operator basis set has been firstly defined as an even-tempered basis set of p-type gaussian functions, with the exponents 3, 1.5, 0.75 and 0.375. The parameters αnm are obtained by a least-squares fit method proposed in 1980 by Nicolas and Durand 30 (so far other EGPs will be considered in the following, the non-local part of ŴEGP fitted with this basis set will herafter be called EGPa). The mathematical form of the non-local part of ŴEGP is very general, and there is no a priori introduction of σ or π parts in the multinuclear core. In principle, that is in the favorable case of a large number of informations defining the reference system, in conjunction with a wide basis set {gn}, the least-squares fit method of Nicolas and Durand should automatically provide the most important |gn >< gm| terms, i.e. with the largest αnm parameters. We will see later that EGPs are rather generated in the context of a small number of informations, which may involves numerical instability.

typo: + sign

The truncated GTO basis set, the effective charge Z eff = 1 for each C # atom, the operator basis set |gn > and the parameters αnm of this operator completely define the pseudo η 6 -C 6 H 6 , which will be herafter called Ar #a . The three SCF π MOs of Ar #a , shown in figure 1c, shed light on the good behavior of the ŴEGP operator. However, one should pay attention that systematic tests must be done in order to check the accuracy and transferability of any new EGP, owing to some numerical instability of the least-squares fit method. The limitation of this method will be commented on in the next section. The EGP routines have been implemented in the Gaussian03 package 31 . During geometry optimizations, the distance between the chemical group replaced by an EGP and the active part is kept fixed at a standard value.

Discussion about the pseudopotential and new prospects

As already commented in Ref. 21, the least-squares fit method used for the determination of optimized αnm parameters may yield chemically irrelevant results albeit the norm is close to zero. This is related to the rather poor sample of informations considered in the least-squares fit process, i.e. a N act × N act Fock operator, (so far N act = 3 for η 6 -C 6 H 6 ). A solution to avoid this artefactual procedure would be to increase the number of informations to be fit. Note that in the shape-consistent extraction of atomic pseudopotentials, the condition imposed to fit the parameters of the pseudopotential is Fpsϕ

ps i = ε i ϕ ps i
, where ε i is the energy of the real valence AO ϕ i and ϕ ps i is the valence pseudo-orbital which fits ϕ i in the valence region. The specification of ϕ ps i (r) for any value of r provides enough informations to fix Ŵps. In the energy-consistent approach 32 , Fpsϕ ps i = ε i ϕ ps i is solved not only for the electronic ground state, but also for electronic excited states, Ŵps being determined in order to fit the energy difference between these states. Note that extraction of EGPs could take advantage of the energy-consistent approach used for defining ECPs. Unfortunately, the method proposed by Nicolas and Durand is not free of numerical instability. We shall now propose a new way to perform the second step of the EGP extraction process. For that, we shall first try to answer the following question: what is the shape of an EGP ? While ECPs on one hand or quantum capping potentials proposed as link atoms in QM/MM methods 33,34 on the other hand have a repulsive shape, easily shown by the plot of Ŵps(r), it is less straightforward to visualize the nature of ŴEGP . A way to more deeply understand the role of this molecular pseudopotential is to formulate the non-local part of the operator ( Ŵ NL EGP ) in a diagonal form:

ŴEGP = NEGP X I=1 Λ I |G I >< G I | - n # X i=1 Z eff i r + n # X i=1 Ŵps(i) (2)
= Ŵ NL EGP - n # X i=1 Z eff i r + n # X i=1 Ŵps(i) (3)
where Λ I and |G I > are the eigenvalues and eigenvectors of the αnm matrix. Each eigenvector |G I > is a linear combination of the N EGP gaussian functions which define the basis set of ŴEGP :

|G I >= NEGP X n=1 c In |gn > (4)
ŴEGP is thus a finite sum of projectors |G I >< G I | with weights Λ I , in addition to Coulomb and ECP operators.

In the case of EGPa, the dimension of the operator basis set is N EGP = 72, with four p-type gaussian functions per atom. Among the 72 eigenvectors of the αnm matrix, only 24 actually play a role in the total hamiltonian since the eigenvalues Λ I of the 48 other eigenvectors cancel. One can however raise doubts on the interest of 21 of the 24 remaining eigenvectors, since the absolute value of the associated eigenvalues are smaller than 0.1 hartree, whereas the three main components of the non-local part of ŴEGP are close to -0.7 hartree. These three eigenvectors, namely |G 1 >, |G 2 > and |G 3 >, are plotted in figure 3. The striking similarity between the main components of Ŵ NL EGP and the three active π MOs of benzene sheds a new light on the role of EGPs. Ŵ NL EGP is nothing else than a level-shift operator which aims at stabilizing the active MOs. In the present case the weights Λ 1-3 are of the order of magnitude of the energy difference between the real MOs φ i and the SCF MOs of C # 6 calculated in the truncated basis set. Nevertheless, Λ I parameters strongly depend on the spatial extension of the operator basis set, namely {gn}. The less diffuse the basis set, the larger the absolute values of the weights Λ I . Although in the present development Ŵ NL EGP only acts in order to stabilize the active occupied MOs in the SCF process, the active set of MOs could easily be extended, for example to virtual MOs in order to suitably fix back-donation effects, such as for the carbonyl ligand 20 . This analysis suggests that the least-squares fit of Nicolas and Durand could be replaced by a new method for optimizing the parameters of Ŵ NL EGP , by directly expressing this operator as a finite sum of projectors. It involves to tune N act Λ I weights, where N act is the number of active MOs (6 or 3 in the case of benzene, whether π * MOs are taken into account or not), instead of optimizing N EGP parameters αnm (72 in the case of EGPa). The lack of informations at the origin of the numerical instability of the Nicolas and Durand fit could thus be solved. Step 2 would now be replaced by the following approach:

1. fix predetermined exponents for the basis set {gn} which defines Ŵ NL EGP typo: + sign typo: + sign (a) calculate the elements of the S matrix:

Snm = Nact X i=1 < gn|φ i >< φ i |gm > n, m = 1, N EGP (5)
(b) diagonalize S and select the N act eigenvectors with eigenvalues closest to 1.

3. optimize the weights Λ I at the same time so that the energies and shapes of the active MOs are well reproduced This proposition will be evaluated in a forthcoming paper.

2.3

The η 6 -C 6 H 6 EGP. A useful tool in the context of organometallic nanoparticles ?

The accurate calculation by means of molecular quantum chemistry methods of the geometric and spectroscopic properties of organometallic NPs is desirable, but for the time being it requires appropriate simplifications. One of them could be supplied by EGPs, in the context of molecular quantum chemistry. Small organometallic clusters are useful in order to validate this suggestion, such as the ruthenium clusters recently investigated in our group 3 . We will consider now benzene ruthenium clusters, which also exhibit hydrido ligands. The water-soluble tetranuclear ruthenium cluster [(H) 4 Ru 4 (C 6 H 6 ) 4 ] 2+ , synthesized for the first time in 1994 35 , is an electron-deficient species (58e) according to electron-counting rules [START_REF] Albright | Orbital interactions in chemistry[END_REF][START_REF] Mingos | Introduction to cluster chemistry[END_REF][START_REF] Fox | [END_REF] the triangular faces of the Ru core, whereas the benzene ligands are coordinated to the ruthenium centers via their π electrons. The Ar # -substituted compound is shown in figure 3.

For a 60e tetrahedral cluster such as (H) 4 Ru 4 (C 6 H 6 ) 4 , the electronic ground state is defined by the occupation of ligands and metal-ligands MOs, four Ru-H MOs, six skeletal Ru-Ru MOs and eight metal-metal MOs which have mainly non-bonding character. We have considered several EGPs, which differ by the basis set which define Ŵ NL EGP : EGPa, EGP b and EGPc are all built according to four p-type gaussian functions, EGPa (ζ: 3.0, 1.5, 0.75, 0.375) being more diffuse than EGPc (ζ: 12.0, 6.0, 3.0, 1.5), EGP b being in-between (ζ: 6.0, 3.0, 1.5, 0.75). EGPa is extracted in the same basis set than the cyclopentadienyl EGP, which turned out to be accurate 19 . It does not exhibit similar qualities in the case of η 6 -benzene ligands, since the DFT calculation applied on the optimal geometry of the reference structure calculated in ref.

3 did not converge. EGP b and EGPc did slightly better, but the SCF process failed to converge at the second step of the optimization process. EGPc, a short range operator, performs better, but it yields poorly accurate geometries. Although the metal-hydrogen distance is close to the reference value (1.82 Å vs. 1.83 Å), the metal-metal distances are smaller than the reference calculation by 0.14 Å (2.63 Å vs. 2.77 Å). The last pseudopotential, namely EGP d , is built in a larger basis set which spreads from 12.0 to 0.375. It also failed at the second step of the optimization process.

This new multicentered EGP is disappointing, whereas according to the previous work on a pseudo-C 5 H 5 it was expected to fairly behave. The main difference between cyclopentadienyl and benzene is the lowest HOMO-LUMO energy gap in the case of C 6 H 6 , which involves a more important π * back-donation in metal-benzene bonding. Four additional EGPs have been extracted, defined by the same basis sets than EGP a-d . A simultaneous fit of the shape and energies of the three highest π MOs and of the three lowest π * MOs has been imposed. EGP a * -d * all accurately reproduce the π * MOs energies (ε * =0.1353×2, 0.3482) and shapes. As can be seen in table 1, no convergence failure did occur in the case of (H) 4 Ru 4 Ar # 4 , thus showing the necessity to suitably describe back-bonding effects in such compounds. The most accurate Ru-Ru bond length is found with EGPa * , the benzene counterpart of the EGP developed for cyclopentadienyl in Ref. 19, with special attention paid to virtual π * MOs in the case of C 6 H 6 . However, the Ru-H bond length is slightly overestimated by 0.044 Å. EGPc * , the short-range operator, is the less accurate since the Ru-Ru bond length is underestimated by 0.104Å, although on the other hand the Ru-H bond length is incidentally identical to the reference value. EGP d * , whose operator spreads on the widest range, performs fairly well, with the Ru-Ru distance smallest than the reference value by 0.019 Å, whereas the Ru-H distance is overestimated by 0.029 Å. Energy diagrams of the occupied MOs close to the Fermi level are reported in figure 4. On the contrary to EGP a-d , the t 1 and e skeletal Ru-Ru MOs lie below the most stable t 2 non-bonding metal MO, in agreement with the reference calculation. The energy differences between all MOs of (H) 4 Ru 4 Ar # 4 obtained with EGP a * -d * are in excellent agreement with those of (H) 4 Ru 4 (C 6 H 6 ) 4 . The main apparent discrepancy is a systematic shift between the energies of the reference and pseudo molecular systems, probably related to the insufficient shielding of the ŴEGP operator. Note that it could be possible to add repulsion by an appropriate parametrization of conventional effective core potentials 33,34 .

The shape of some selected MOs, namely the six skeletal (or metal-metal) MOs and the four metal-hydrogen MOs are shown in figure 5. As already observed in the case of benzene, there is a remarkable similarity between the reference MOs of (H) 2+ compound is a more tricky system, with the lowest electronic triplet state T 1 which lies only 3.8 kcal/mol above the singlet ground state S 0 . As can be seen in table 1, all EGPs extracted from occupied and virtual π MOs of benzene, namely EGP a * -d * , perform well, both in terms of geometry and singlet-triplet energy difference (∆E ST ). EGPa * and EGP d * are the most accurate, whereas ∆E ST as well as Ru-Ru bond lengths are underestimated with EGPc * . As already mentioned, that short-range EGP is the less reliable. However, it behaves better than a fifth EGP, which spans a more diffuse space (EGPe * , ζ: 3.0, 1.5, 0.75, 0.375, 0.1875, 0.09375). The diffuse functions are probably responsible of the attractive character of this EGP, which manifests itself by the optimal geometry of [(H) 4 Ru 4 Ar #e * 4 ] 2+ , with two edge-bridging H and two terminal-H instead of four face-capping H atoms. Each hydride in top position is located very close to an Ar # fragment. ∆E ST , found to be 216 kcal/mol (!), is also an indication of the unreliability of this EGP. It confirms that long-range Ŵ NL EGP operators act both on the fictitious part and unexpectedly on the active part close to the active/spectator boundary.

A joint EGP/Molecular Dynamics study

Ab initio Molecular Dynamics

From a quantum chemical point of view, the characteristic values of a physico-chemical reaction are obtained from the knowledge of specific point on the potential energy surface associated to the process. In most cases, kinetic or thermodynamic data are determined using calculated Gibbs free energies for which the entropic contributions are estimated using the harmonic approximation [START_REF] Mcquarrie | Statistical Mechanics[END_REF] . In some cases, this harmonic approximation leads to results in poor agreement with experimental knowledge, since anharmonic contributions must be taken into account. This can be performed as a correction of harmonic approximation and is implemented in most standard quantum chemical codes [START_REF] Barone | [END_REF] . Another way to explore potential energy surface and to estimate the temperature effect added to electronic energies is based on a dynamical approach. In this specific case, temperature is explicitly controlled by simulating the canonical ensemble by means of thermostats 41 in a molecular dynamics simulation scheme. The thermodynamic and kinetic data can then be determined by employing a constraint 42 that must correspond to the reaction coordinate.

Molecular dynamic simulations were initially based on empirical interaction potentials determined from experimental data. Since the appearance of the Car-Parrinello molecular dynamics technique 43 , more general molecular dynamics approaches have known a tremendous development both from an application and development point of view. Rather than using empirical interaction potentials, the electronic energy and the forces on nuclei that are necessary to resolve Newton classical equations of propagation are computed by wavefunction-based or density-based methods. Mainly two types of molecular dynamics schemes can be highlighted: The so-called Born-Oppenheimer [START_REF] Bolton | Modern Methods for Multidimensional Dynamics Computation in Chemistry, chapter Direct Dynamics of Reactive Systems[END_REF][START_REF] Millam | [END_REF][46][START_REF] Marx | Modern Methods and Algorithms of Quantum Chemistry, chapter Ab initio molecular dynamics: Theory and implementation[END_REF] and Car-Parrinello molecular dynamics techniques. In the first case, the electronic energy and its associated nuclear derivatives are calculated after an energy minimization procedure that must be fulfilled at each step of the molecular dynamic simulation. This procedure is necessary in order to obtain stable behavior of the molecular dynamics simulation relative to the molecular system total energy for long time propagation. In the case of Car-Parrinello molecular dynamics, a new degree of freedom is introduced, namely the wavefunction. Consequently, the nuclear degrees of freedom and the wavefunction evolve according to the Newton law and the SCF convergence procedure can be avoided.

It must be noted that most of the Car-Parrinello algorithms are based on the density functional theory (DFT) approach coupled with a plane wave basis set description of the electronic structure. While this approach presents many advantages, the use of plane waves may be ill-suited for study of non-periodic chemical systems to describe cusps and high inhomogeneities. Consequently, following the initial proposition of E. Carter [START_REF] Hartke | [END_REF]49 , we have developed a few years ago an algorithm based on Car-Parrinello Lagrangian but using atom centered basis functions (more precisely Gaussian type basis sets) for the electronic structure calculation part of the procedure 13 . This algorithm has been coupled with several tools in order to improve our way to explore chemical reaction and physico-chemical processes. For instance, Nosé-Hoover chain of thermostats 41,50 have been implemented in our molecular dynamics code in order to perform molecular dynamic simulations at finite temperature. This specific tool has allowed to explore molecular dynamic treatment of Infra-Red [51][52][53] , NMR 11 and UV 12 spectroscopic properties. Moreover, in order to explore chemical reactivity, it has been necessary to introduce a statistical mechanical procedure that samples the phase space along a defined reaction coordinate and allows the estimation of relative free energies (i.e. free energy barriers of a chemical process). In this context, two types of chemical reaction have been studied 14,54 . These studies have permitted to demonstrate the failure of the "static" approach based on the harmonic approximation used for Gibbs free energy calculations in the case of flat potential energy surfaces. . This reaction is a prototype of reactions in lanthanide chemistry, extensively studied by quantum chemistry methods 56,57 . Theoretical studies have been extended for the activation of other bonds (C-H, C-F, Si-H, etc.), leading to successfull experimental and computational works [58][59][60] . Quantum chemistry can indeed provide very relevant insights on multi-step reactions. However such studies are generally restricted to the search of points potential energy surfaces. Although MD simulations could bring additional informations, the study of chemical reactions by AIMD is still a challenge. The CPMD method has however been applied to study the hydrogen exchange catalyzed by a lanthanocene complex 14 . In this study both cyclopentadienyl rings have been substituted by chloride ligands in order to reduce the computational cost, taking advantage of the fair isolobality of Cp and Cl ligands. Interestingly, it has been shown that anharmonicity effects are of great importance for the evaluation of the activation barrier of this reaction and that the harmonic approximation which is usually used to estimate vibrational entropic contributions may fail.

EGPs could be very appropriate tools to efficiently model a spectator ligand and to drastically reduce the associated computational cost. It has already been recalled in section 2.1 that the pseudo η 5 -C 5 H 5 EGP (Cp # ) has successfully been used to mimic a real η 5 -C 5 H 5 cyclopentadienyl ligand 61 . For these previous reasons, we focus on an AIMD study of the hydrogen exchange reaction catalyzed by a Cp 2 LaH complex, where two Cp # EGPs fill in for the real Cp rings:

Cp # 2 La -H * + H -H -→ Cp # 2 La -H + H * -H
Evaluation of free energy activation barriers with AIMD. In the case of free (i.e. non constrained) ab initio molecular dynamics simulations, two important thermodynamic quantities for a chemical process, namely entropy and free energy, cannot in general be derived from a statistical average. These global properties depend on the extent of the phase space accessible to the molecular system during the simulation [START_REF] Mcquarrie | Statistical Mechanics[END_REF][START_REF] Frenkel | Understanding Molecular Simulation -From Algorithms to Applications[END_REF] . The probability of finding the system in a transition state region is so small that, therefore, no relative free energies of a chemical process can be calculated by natural molecular dynamics simulation. Several statistical mechanical procedures circumvent this problem. For instance, the original "umbrella-sampling" method [START_REF] Torrie | [END_REF]64 , which adds a coordinate-dependent potential, leads to the system to sample a specific region of the phase space. Other approaches which are based on similar ideas also allow the estimation of free energies for chemical processes, for instance the recent metadynamic methodology [START_REF] Laio | Proc. Natl. Acad. Sci[END_REF][START_REF] Iannuzzi | [END_REF] or the "blue-moon" method 67,68 . The latter has been used to estimate free energy differences by sampling the phase space along a defined action pathway. By means of thermodynamic integration 62 , free energy difference is obtained as:

F(ξ 2 ) -F(ξ 1 ) = ξ2 Z ξ1 dξ ′ fi ∂H ∂ξ fl ξ ′ ( 6
)
where H is the Hamiltonian of the system, ξ(r) is a chosen reaction coordinate and . . . ξ ′ is an ensemble average at ξ(r) = ξ ′ . This conditional average could be estimated by a time average over a constrained trajectory with the chosen reaction coordinate fixed at a specified value. Recently, generally applicable expressions for the average force of the constraint f ξ = -∇ ξ F, which is the opposite of the integrand in eq. ( 6), have been outlined 42,69 . It has been shown that the bias introduced by the constraint can be removed by considering the following expression:

f ξ = D Z -1/2 (-λ + k B T G) E ξ ˙Z-1/2 ¸ξ (7)
where k B is the Boltzmann constant, T is the temperature, Z and G are respectively weight and correction factors, λ is the Lagrange multiplier of the associated holonomic constraint. Therefore, the obtained free energy differences might not be dependent on the choice of the constraint and accordingly to the Jarzynski inequality it leads to an upper bound for the activation barrier.

Computational details. Potential energy calculations. Calculations were carried out at the DFT level of theory using the hybrid B3PW91 70,71 functional. A large core relativistic effective potential (ECP) optimized by the Stuttgart group and the corresponding basis sets have been used for the lanthanum center 72 , whereas hydrogen atoms have been described with a 6-31G(d,p) double-ζ quality basis set 73 . The effect of Cp has thus been mimicked by the previously published EGP 19 , five electrons are allowed to occupy the space spanned by p atomic orbitals (α = 0.312) centered at the edges of the pentagonal ring. The computations of potential energy and the nuclear gradients have been carried out with the Gaussian 98 suite of programs modified to incorporate EGPs 74 . Free energy differences at room temperature have been estimated by means of the "Blue-Moon" approach. CPMD simulations with Gaussian type orbitals have been performed with our own code 13 . Fictitious electronic mass was set to 170 a.u. and equations of motion have been integrated with a time step of 0.25 fs by means of a velocity Verlet scheme 75 . These simulations were performed in the canonical ensemble with Nosé-Hoover chains of thermostats 41 and holonomic constraints associated to the reaction coordinate have been applied considering the method of undetermined parameters. For each trajectory, a thermalization procedure has been performed for at least 6 ps, followed by a production simulation achieved for additional 5 to 6 ps. The property of interest, namely the force of the constraint, has been averaged as a function of time during the production step and its convergence has been checked at the end of the whole process.

Description of geometrical constraints. The choice of the constraint which represents the reaction coordinate can be a difficult task. Such constraint can be described by several geometrical parameters. Moreover, an incorrect choice of the type of constraint can induce a bias for free energy calculation. For the considered reaction, it is not obvious to find a function that will constrain the transition state geometry which can be mainly described by the concerted breaking and forming of four bonds. For this study, we have employed a well-tried type of constraint: it corresponds to the projection of the H a H b vector on the H a H c one (fig. 6). This type of constraint has already been successfully applied for the hydrogen exchange catalyzed by Cl 2 LaH 54 .

La

Cp # Cp # H c H a H b
Fig. 6 Atom definition for the geometrical constraint AIMD linked with the EGPs. Ab initio molecular dynamics simulations can take advantage of the reduction of the cost involved by EGPs. On the other hand, the current formulation of EGPs (eq. 1) does not include an explicit repulsive term in charge to mimic the screening effect of the inactive MOs of the real chemical group. As a consequence, the calculation of the gradient of the potential energy with respect to the position of each center which define a given ŴEGP operator would be meaningless. This limitation binds the EGPs skeleton to be frozen during the molecular dynamics simulations. Moreover for the same reasons, atoms directly "bonded" to EGPs cannot be relaxed during a geometry optimization procedure, or allowed to move during a MD trajectory. Consequently, only the positions of hydrogen atoms were allowed to evolve during the MD simulations and the geometrical structure of the Cp # 2 La fragment has been kept frozen. This leads to a modification for the number of degrees of freedom. Thus, the Cp # 2 La fragment could be seen as two fixed "rugby balls" linked by the top (fig. 7). The H - 3 fragment is allowed to evolve in the field created by these associated rugby balls. Conservation of the total energy. In order to illustrate the implementation of the approach discussed above, we focus our attention on the quality of the potential energy surface explored by the dynamics. In particular, the conservation of the total energy is carefully examined. One may note that the fictitious kinetic energy associated to the classical propagation of the Kohn-Sham orbitals in CPMD is not taken into account for the total energy. The latter includes the potential energy, the nuclear kinetic energy and the contributions of the Nosé-Hoover chain of thermostats. One can suspect that the classical propagation of the Kohn-Sham orbitals by means of the Car-Parrinello procedure could be approximate, in particular in the transition state region where the concerted breaking and forming of four bonds occurs. The evolution of the relative total energy as a function of time for the trajectory which corresponds to the exploration of the transition state region (ξ = 0.5) is shown on figure 8. After 5 ps of simulation time, the total energy deviation is lower than 2 • 10 -5 a.u. The range of CP dynamics' error for this trajectory remains below a very satisfactory limit. Moreover considering the employed integrator -a velocity Verlet scheme -and the integration time step -0.25 f s -this range of error is very close to the one expected for a reference Born-Oppenheimer MD simulation. The similar range of error has been measured for the trajectories performed for different values of the constraint ξ, which represents the reaction coordinate. Finally this satisfactory behavior and in particular the conservation of total energy leads us to be confident in the presented scheme, combining CPMD with GTO basis sets and EGPs.

La

Estimation of the free energy activation barrier. Using the previously discussed "blue-moon" method, the activation barrier has been estimated with a set of trajectories where the constraint ξ has been fixed to different values in a range of 0.05 -0.5. Two properties are plotted on figure 9 as a function of the constraint ξ: (i) the opposite of the force of the constraint (eq. 7) and (ii) the free energy difference computed after integration of the force (eq. 6). The calculated free energy barrier is found to be 7.1 kcal • mol -1 , whereas the corresponding barrier computed considering the Cl 2 LaH catalyst with the same AIMD methodology 54 is 6.2 kcal • mol -1 . These close values are in agreement with the fact that the potential energy activation barriers ∆rE # are also close (1.43 kcal • mol -1 and 1.93 kcal • mol -1 with Cp # 2 LaH and Cl 2 LaH, respectively). These kinetic data lead us to affirm that anharmonicity effects are of great importance for this reaction. As a matter of fact, they are far away from the corresponding free energy activation barrier computed in the case of the chloride ligand and with the harmonic approximation, ∆rF # = 11.2 kcal • mol -1 . This discrepancy has already been attributed to anharmonicity effects and the failure of the harmonic approximation usually employed to compute vibrational entropic contributions for free energy. This is due to the flatness of the potential energy surface, a significant number of normal modes of low frequencies (below 500 cm -1 ) being calculated for the transition state 54 .

These anharmonicity effects lead to deviations for some geometrical parameters. Table 2 gathers some selected distances for the transition state. Significant differences occur between the optimized geometry on the potential energy surface and the corresponding values averaged on the ξ = 0.5 MD trajectory. These discrepancies strengthen our doubts concerning the relevance of the harmonic approximation for the free energy activation barrier estimation for this hydrogen exchange reaction.

Conclusion

Multicentered EGPs have proven their ability to generate a relevant ligand-field in organometallic clusters. It has more specifically been demonstrated in the case of η 6 -benzene, which has been added to the EGP library. The remarkably accurate description of the electronic structure, geometry and energy of [Ru 4 (H) 4 (C 6 H 6 ) 4 ] n+ is very promising. It opens the road to the use of such molecular pseudopotentials in the context of organometallic NPs, where new methodological developments are required in order to reduce the computational cost of DFT calculations. Considering that dynamical effects occur at the surface of organometallic NPs, we have also checked the behavior of EGPs in the context of ab initio molecular dynamics. A prototype system has been considered, namely an hydrogen exchange reaction catalyzed by a Cp 2 LaH complex. Energy conservation during 5 ps of this joint EGP/MD is not violated, thus confirming that EGPs do not introduce spurious trajectories in MD simulations. In addition, the barrier height of this reaction, evaluated within the so-called blue-moon algorithm, is in close agreement with reference data. In other words, the part of the potential energy surface associated to the reaction is accurately obtained with EGPs, even the anharmonic part. These results confirm that EGPs could be useful tools in the context of large organometallic clusters. An extension of the EGP library would however require to perform secure extractions of the molecular pseudopotential. Some guidelines are also provided in order to build a new generation of EGPs in a more routine way.

Fig. 1

 1 Fig. 1 Design of an η 6 -C 6 H 6 EGP (a) Reference HF π MOs in the 6-31G(d,p) basis set ; (b) HF MOs obtained in the minimal truncated basis (i.e. one p GTO / C # atom, see section 2.1), with carbon ECPs in order to take into account the repulsion of the 1s core, but without the non-local part of ŴEGP (i.e. without P P αnm|gn >< gm|); (c) HF π MOs obtained with EGPa. The MO energies are given in hartree.

Fig. 2

 2 Fig. 2 Graphical representation of the main terms |G I > of Ŵ NL EGP operators. EGPa and EGPa * are respectively defined by |G 1-3 > and |G 1-6 >. The weights Λ I (in hartree) are given below for each EGP.

  2. find normalized functions G I as close as possible to the N act reference MOs φ i . It can be done in the same spirit than the determination of the MVPOs, which are as close as possible to the reference MOs, but which are expressed in the truncated basis set {fq }. The G I would be obtained in two steps :

Fig. 3

 3 Fig. 3 Tetrahedral Ru 4 Ar # 4 (µ 3 -H) 4 cluster.

Fig. 5 4 ,

 54 Fig. 5 Comparison of the skeletal (T 1 , E and A 1 symmetry) and Ru-H (A 1 and T 2 symmetry) MOs of (H) 4 Ru 4 (C 6 H 6 ) 4 and (H) 4 Ru 4 Ar #d * 4 , with Ar #d * = 1p per C # (ζ = 0.332) + EGP d * .

3. 2

 2 Reactivity of a lanthanocene complex: a test case Motivations. Lanthanide chemistry has witnessed an increasing interest over the last three decades. In particular, Watson and Parshall have shown that Cp * 2 LuH (Cp * = C 5 Me 5 ) or Cp * 2 LuMe can react rapidly with H 2 55

Fig. 7

 7 Fig. 7 Equivalence between the Cp # 2 La fragment and two rugby balls.

Fig. 8

 8 Fig.8Evolution of the relative total energy for a trajectory as a function of time. The constraint is fixed to ξ = 0.5, which corresponds to the exploration of the transition state region.

1 Fig. 9

 19 Fig.9Evolution of the opposite of the force -∇ ξ F (dotted line) of the constraint and the free energy difference ∆rF (solid line) obtained after integration.

Table 1

 1 Comparison of the properties of [(H) 4 Ru 4 (C 6 H 6 ) 4 ] n+ and [(H) 4 Ru 4 Ar # 4 ] n+ (n = 0 or 2) for EGP defined according to occupied a and lowest unoccupied π MOs (EGP a * -d * ). Distances in Å, energy differences in kcal/mol. [(H) 4 Ru 4 (C 6 H 6 ) 4 ] n+ [(H) 4 Ru 4 Ar # 4 ] n+

				EGPa *	EGP b *	EGPc *	EGP d *
				n = 0, ground state (T d )
	R RuRu	2.77	2.76	2.72	2.66	2.75
	R RuH	1.83	1.87	1.84	1.83	1.86
				n = 2, singlet state S 0 (Cs)
	R RuRu	2.84×4	2.85×4	2.79×4	2.74×4	2.83×4
			2.59×2	2.60×2	2.55×2	2.50×2	2.58×2
	R RuH	b	1.84×2	1.89×2	1.87×2	1.87×2	1.88×2
			1.89	1.90	1.88	1.84	1.89
			n = 2, triplet state T 1 (∼C 3v )
	R RuRu	2.75×6	2.77×6	2.71×6	2.65×6	2.75×6
	R RuH	1.85×3	1.89×3	1.87×3	1.85×3	1.88×3
				n = 2, ∆E ST
			3.8	4.8	2.8	1.0	3.8

a : severe numerical instabilities have been observed for most EGPs built by taking into account occupied π MOs only. EGPa: SCF not converged, EGP b and EGP d : geometry optimization not converged, EGPc: R RuRu 2.63 and R RuH 1.82 Å. b : averaged values

  Energy-level diagrams of the occupied MOs of (H) 4 Ru 4 (C 6 H 6 ) 4 and (H) 4 Ru 4 Ar # 4 close to the Fermi level, calculated for the optimal geometry of (H) 4 Ru 4 (C 6 H 6 ) 4 . : non-bonding or weakly bonding Ru-Ru MOs (E and T 2 symmetry); q: skeletal Ru 4 MOs (T 1 , E and A 1 symmetry); v: C 6 H 6 or Ar # MOs; ♦: C 6 H 6 -Ru or Ar # -Ru MOs; : Ru-H MOs (A 1 and T 2 symmetry).

		-0.15							-0.15
		-0.2							-0.2
	Energy / hartree	-0.3 -0.25							-0.3 -0.25
		-0.35							-0.35
		-0.4							-0.4
		-0.45							-0.45
		Ref	EGP a*	EGP b	EGP b*	EGP c	EGP c*	EGP d	EGP d*
	Fig. 4
									-levels

4

  Ru 4 (C 6 H 6 ) 4 and the MOs of (H) 4 Ru 4 Ar #d * 4 , with the exception of the triply degenerate Ru-H MOs. While the totally symmetric Ru-H MO of (H) 4 Ru 4 Ar # 4 is almost identical to the corresponding MO in (H) 4 Ru 4 (C 6 H 6 ) 4 , the three other MOs of the model system do not match the reference MOs due do a slight mixing of the Ru 4 H 4 fragment MOs with σ MOs of the (C 6 H 6 ) 4 fragment. It is not possible to obtain such mixing, since σ(C 6 H 6 ) MOs cannot be described by an Ar # EGP. The energy difference between the Ru-H MOs in (H) 4 Ru 4 Ar # 4 is nevertheless analogous to (H) 4 Ru 4 (C 6 H 6 ) 4 , owing to the small σ(C 6 H 6 ) contribution. The energies and shapes of the MOs as well as the geometrical properties of the tetrahedral Ru 4 cluster show that most of the η 6 -C 6 H 6 EGPs are able to create an accurate ligand field around the Ru 4 H 4 moiety. The [(H) 4 Ru 4 (C 6 H 6 ) 4 ]

Table 2

 2 Selected distances (in Å) for the hydrogen exchange transition state Cp # 2 La(H 3 ), see fig.6for atom definition. a First column: bond length from a "static" TS geometry optimization; b second column: average values and standard deviations from the ξ = 0.5 MD trajectory.

		Static a	Dynamic b
	r(La -H a )	2.284	2.350 ± 0.069
	r(La -H b )	2.257	2.313 ± 0.054
	r(H a -H b )	1.015	1.048 ± 0.027

5

, which corresponds to the exploration of the transition state region.

Acknowledgements We thank the CALcul en MIdi-Pyrénées (CALMIP) for generous allocations of computer time. Laurent Maron is member of the Institut Universitaire de France (IUF). Financial support by the CNRS and the ANR (SIDERUS project, ANR-08-BLAN-0010-01 and ANR-08-BLAN-0010-02) is gratefully acknowledged. This work is dedicated to the memory of Dr Jean-Pierre Daudey, who was affectionately called by us "Ze Guru". EGP extractions and applications for benzene, as well as the corresponding sections of this article were essentially achieved during a wonderful week in the département of Ardèche, in April 2009. Je me souviens 76,77 d'une fameuse fondue au LPQ autour de la table à saucissons. Je me souviens des rires et des doutes existentiels et scientifiques. Je me souviens du xylophone à roulettes d'Emilie que tu faisais tintinabuler pour nous convier à prendre le thé. Je me souviens de la trousse Babar qui t'accompagnait à chaque réunion. Je me souviens des impulsions scientifiques décisives que tu as données à ces travaux, après tant d'autres. Je me souviens de ta mélancolie.