Structures and stabilities of PAH clusters solvated by water aggregates: the case of the pyrene dimer H. Leboucher, 1 A. Simon, 1 and M. Rapacioli 1 Laboratoire de Chimie et Physique Quantiques LCPQ/FERMI, UMR5626, Université de Toulouse (UPS) and CNRS, 118 Route de Narbonne, F-31062 Toulouse, France (*Electronic mail: mathias.rapacioli@irsamc.ups-tlse.fr) (Dated: 17 February 2023) Whereas clusters made of polycyclic aromatic hydrocarbon and water monomers are relevant objects in both atmospheric and astrophysical science, little is known about their energetic and structural properties. In this work, we perform global explorations of the potential energy landscapes of neutral clusters made of two pyrene units and one to ten water molecules using a Density-Functional based Tight-Binding (DFTB) potential followed by local optimisations at the Density-Functional Theory level. We discuss the binding energies with respect to various dissociation channels. It shows that cohesion energies of the water clusters interacting with a pyrene dimer is larger than that of the pure water clusters, reaching for the largest clusters an asymptotic limit similar to that of pure water clusters and that, while the hexamer and octamer can be considered as magic numbers for isolated water clusters, it is not the case anymore when they are interacting with a pyrene dimer. Ionisation potentials are also computed making use of the Configuration Interaction extension of DFTB and we show that in cations, the charge is mostly carried by the pyrene molecules.

I. INTRODUCTION

The study of interactions between Polycyclic Aromatic Hydrocarbon (PAH) clusters and water aggregates is of great interest in atmospheric science, astrophysics and astrochemistry. Indeed, PAH molecules can be found in our Earth atmosphere where they can become part of atmospheric aerosols, which are known to have a significant impact on the Earth's radiation. 1 PAH clusters can be regarded as good models for such aerosols on which small atmospheric molecules such as water condense and grow. [START_REF] Vinklárek | [END_REF] More generally, such aerosols play an important role in atmospheric chemistry, as they act as surface catalysts for heterogeneous reactions. 3,4 AstroPAHs have been of significant interest since the proposal in the eighties that they were the carriers of the Aromatic Infrared Bands (AIBs), a set of mid-IR emission bands observed in many regions of the interstellar medium (ISM). 5,6 They would account for up to 20 % of the total carbon. 7,8 Although PAHs have been believed to be ubiquitous in the ISM for decades, motivating experimental and theoretical spectroscopic studies in order to identify a specific PAH molecule, 8 it is only very recently that specific PAH molecules, the two isomers of cyano-naphthalene, have been successfully detected based on their rotational spectra. 9 In the denser regions of the ISM, PAHs may freeze on the icy mantle of interstellar dust grains 10 where they can interact with other molecules like water and undergo physical and chemical processes induced by UV photons emitted by nearby stars such as oxidation reactions. 11 Besides, astronomical observations revealed the presence of nanograins inside molecular clouds as precursors of PAHs and PAH clusters were proposed as good analogs for these nanograins. [12][13][14] Small molecules such as H 2 O may condense on such grains and PAH clusters with adsorbed water molecules could then be considered as models for carbon nanograins in planet-and starforming regions. Clusters in space conditions are difficult to study experimen-tally due to low binding energy leading to the need to trap the species in low temperature conditions. Water-PAH neutral clusters were studied in cryogenic environments such as rare gas matrices. [15][16][17][18][19] They were also formed in the gas phase and characterised by rotational spectroscopy experiments complemented with ab initio calculations. [20][21][22] Large cationic clusters were also formed in the gas phase thanks to a cluster source and analysed by mass spectrometry. 23 While technical issues limit the fields of experimental studies, theoretical modelling allows the studies of a variety of astrophysical analogs in the gas phase. From the theoretical side, modelling neutral molecular clusters of PAH and/or water units requires an accurate description of weak forces such as dispersion and Coulomb interactions between the molecular electronic densities. Indeed, the fine balance between these interactions and the Pauli repulsion drives the structural properties of these clusters, which are far from being intuitive. For the smallest clusters or when the potential energy surface (PES) exploration is restricted to few single points or to local optimisations, ab initio quantum electronic methods can still be used such as MP2, [24][25][26] or DFT schemes. 26,[28][29][30][31][START_REF] Bartolomei | The Journal of Physical Chemistry C[END_REF][33][34] Indeed, when a global exploration of the PES is intended to be done, it requires so many energy and forces calculations that only semi-empirical methods can provide the required computational efficiency. Among them, the Density-Functional based Tight-Binding (DFTB) scheme [35][36][37] appears as a good compromise between computational cost and accuracy. Let us however mention that DFTB requires some specific corrections, to properly account for intermolecular interactions. [38][39][40] In the case of cations, the proper treatment of charge resonance between the different units is a challenging task, even at the DFT level due to the well-known selfinteraction error, 41,42 which can be achieved by the DFTB-Configuration Interaction (DFTB-CI) scheme. [43][44][45] Recently, the local exploration of the PES for some PAH dimers in interaction with small water clusters was achieved. The structures and energetics of water clusters confined between two PAH molecules were determined for a series of PAHs smaller than coronene (C 24 H 12 ). 34 The PES of the neutral and cationic anthracene dimer in interaction with 1 to 4 water molecules was explored locally in complement of vacuum-UV photo-dissociation experiments. 46 This study showed in particular the dependence of the structure of the most stable isomer as a function of size and charge states. This illustrates the need for a more efficient exploration scheme in order to elucidate the complex PES of larger PAH in interaction with larger water clusters. The present study can be regarded as a follow-up of previous PES explorations performed for PAH clusters and water-PAH clusters at the DFTB level. Studies devoted to pure pyrene clusters showed that the neutral clusters are made of compact assemblies of sub-blocs containing up to three units whereas cations present a charged dimer or trimer core surrounded by neutral units. [START_REF] Dontot | The Journal of Physical Chemistry A[END_REF] Previous investigations of pure water clusters and water clusters interacting with one PAH unit showed that DFTB was able to provide a correct description of water-water and PAH-water interaction. 48 It was able to reproduce the energetic order of the isomers of the water hexamer with respect to correlated ab initio methods. 49 Finite temperature simulations showed the diffusion of the water monomer and dimer onto a coronene surface 48 at low temperature while the interaction with the PAH showed the predominance of 2D structures for the water hexamer and octamer at low temperature. 49,50 The influence of the adsorption of such clusters on a single PAH on their heat capacities was also investigated. 51 It showed that the presence of a PAH support results in a strong decrease of the melting temperature in the case of the water octamer, whereas an increase or decrease can be observed in the case of the hexamer, depending on the PAH. The efficiency of the DFTB and DFTB-CI method allowed to describe the structures and energetics of PAHs from benzene to coronene adsorbed on large water clusters modelling water ice (amorphous and crystalline). It was shown that such interactions led to some small variations of the ionisation energies of the PAHs, the trend (increase or decrease) depending on the local interactions between the PAH and the ice surface. 52 Taking into account tens of possible local PAH-water ice organizations, it was possible to reproduce the shifts of the IR features of the O-H dangling bond at the surface of water ice upon interaction with the PAH 53 , which is of interest to interpret further observations by the James Webb Space Telescope launched in December 2021. 54 In this work, we will search for the most stable structures of neutral clusters made of two pyrene units and n = 1-10 water molecules, combining global explorations at the DFTB level and local optimisations at the DFT level. Details on the exploration strategy and the electronic potentials used are given in the following section. In Section III, devoted to the results analysis, we discuss the structural and energetic properties, as well as the vertical ionisation potentials, computed for the most stable and relevant secondary isomers. Finally, summary and outlooks are given in the last section.

II. METHODOLOGY

The first step of this work consisted in the global exploration of the PES of the (Py) 2 (H 2 O) n (n = 1 -10) molecular clusters, performed using a Parallel Tempering Monte Carlo (PTMC) algorithm. This procedure is detailed in section II B. For this global exploration, the electronic structure was determined using a dedicated DFTB potential as detailed in section II A 1. The most stable structures were refined at the DFT level as specified in section II A 2.

A. Electronic structure methods

DFTB

Briefly, DFTB is an approximated DFT scheme [35][36][37] that has been the object of many developments addressed in several reviews. [55][56][57][58] The present work was achieved in the framework of the Self-Consistent Charge (SCC-)DFTB formalism. 37 In this scheme, the Kohn-Sham (KS) theorem is applied and molecular orbitals are expressed in a minimal valence basis set of atomic orbitals: φ i (r) = ∑ µ c iµ χ µ (r). In addition, the exchange-correlation energy is expanded up to the second order with respect to density fluctuation around a reference density ρ( r) = ρ 0 ( r) + δ ρ( r) where ρ 0 ( r) corresponds to the superposition of frozen atom-like densities centered at each nucleus. After some developments, the SCC-DFTB energy is expressed as:

E = ∑ α>β E rep αβ + occ ∑ i ∑ µν n i c µi c νi H 0 µν + 1 2 ∑ αβ γ αβ q α q β (1)
The first term is a pair repulsive potential contribution between atoms α and β . In the second term, n i is the occupation of the molecular orbital i and H 0 µν is the matrix element of the KS operator at the reference density ρ 0 ( r). The last term arises from the second-order expansion and is expressed as a function of atomic charge (q) fluctuations. Three-center integrals are neglected, which allows to parameterise all matrices elements from DFT atomic pairs calculations. A correct description of molecular clusters necessitates a good account of long range interactions. First, an empirical dispersion contribution was added to the total energy: 39,40,59

E disp = -∑ α =β f damp (R αβ ) C 6 αβ R 6 αβ (2)
where the C 6 αβ dispersion coefficient values were taken from the paper by Goursot et al. 60 consistently with our previous works. 40,[START_REF] Dontot | The Journal of Physical Chemistry A[END_REF]51,[61][62][63] Second, long range electrostatic interactions were improved by replacing the original Mulliken charges of the SCC-DFTB scheme with the Class IV-Charge Model 3 (CM3) charges 64 in order to improve the description of the polarity of the bonds 40 . In its original version, the partial CM3 charge of an atom α is defined as:

q CM3 α = q Mull α + ∑ β =α (D t α t β B αβ +C t α t β B 2 αβ ) (3)
where B αβ is the Mayer's bond order between atoms α and β , and D t α t β and C t α t β are empirical parameters. In this work, as in our previous studies, 40,[START_REF] Dontot | The Journal of Physical Chemistry A[END_REF]51,[61][62][63] we considered only first order corrections (C t α t β = 0). In the SCC-DFTB scheme, several types of parameters can be used. In the present work, we used the "MAT" set of parameters [START_REF] Frenzel | Semi-relativistic, self-consistent charge Slater-Koster tables for densityfunctional based tight-binding (DFTB) for materials science simulations[END_REF] in conjunction with CM3 charges and dispersion corrections as it was shown to provide a fair description of the structures and energetics of pyrene clusters. [START_REF] Dontot | The Journal of Physical Chemistry A[END_REF] We used the previously determined value of D CH = 0.1 and we took D CO = 0.0. We determined a D OH value of 0.189 to reproduce with the "MAT" parameters the SCC-DFTB partial charges of a water molecule previously validated with the "BIO" set of parameters. 61 Details of the parameters and benchmark with water hexamers can be found in the SI (Tables I andII). In the following, the SCC-DFTB potential involving dispersion and charge correction will simply be referred to as DFTB.

Computing ionisation potentials requires a proper treatment of charge delocalisation in the cationic molecular clusters, which is a challenging task at the DFT level (and consequently at the DFTB level), due to spurious self-interaction errors 41,42 which tend to over-stabilise cations with respect to neutrals. 43 We have developed, over the past decade, the DFTB-CI (Configuration Interaction) scheme [43][44][45] , derived from the similar approach at the DFT level. [START_REF] Wu | [END_REF] We showed that it allows the computation of structural, energetic and spectral properties of cationic molecular clusters [START_REF] Dontot | The Journal of Physical Chemistry A[END_REF][67][START_REF] Lei | The Journal of Physical Chemistry Letters[END_REF][69][70] and in particular ionisation potentials of pyrene clusters. [START_REF] Joblin | The Journal of Physical Chemistry Letters[END_REF] In this scheme, the wavefunction is expressed on a basis of charge-localised configurations:

Ψ + = ∑ A C A Ψ + A (4)
where Ψ +

A is built from a constrained DFTB calculation 44,72 enforcing the charge to be fully localised on the molecular fragment A. The charge delocalisation is treated in a second step within a configuration interaction scheme:

∑ B (H CI AB -E 0 S CI AB )C B = 0 (5)
where H AB and S AB are the Hamiltonian and overlap matrices expressed in the basis of the charge localised configurations. Details on their computations can be found in reference. [START_REF] Wu | [END_REF]

DFT

Local DFT optimisations were performed using the M06-2X 73 functional, a hybrid functional that includes dispersion and was shown to be particularly suitable for main group energetics and chemistry, as well as description of non-covalent interactions. Given the large size of the investigated molecular clusters, we used a double-zeta basis set that includes polarisation orbitals, D95v(d,p). [START_REF] Dunning | Modern theoretical chemistry[END_REF] Such a level of theory, referred to as DFT in the rest of manuscript, was previously used for water-PAH clusters. 48

B. Global exploration of potential energy surfaces

The PES exploration consisted in the following three steps performed independently for each (Py) 2 (H 2 O) n stoichiometry (n = 1 -10). In the first step, we defined starting point geometries and performed PTMC explorations [START_REF] Sugita | [END_REF]76 of the DFTB PES, keeping the intramolecular modes frozen as performed in previous works. [START_REF] Dontot | The Journal of Physical Chemistry A[END_REF]51 In PTMC, several MC trajectories are run at the same time, each with its own temperature, and configurations can be exchanged between these trajectories at some regular steps. We used the all-exchange strategy where exchanges can be performed between all trajectories. 77,78 In the low temperatures trajectories, the bottom of the potential wells are well explored whereas in the high temperatures trajectories, energetic barriers can be crossed. Performing exchanges between them allows to get advantage of these two behaviours, strongly improving ergodicity. In practice, ten trajectories of 100,000 steps each were run. They correspond to ten temperatures ranging from 40 to 300 K following geometric progression. Exchanges were allowed every ten steps. For n = 2 to n = 7, the structures were confined in a sphere of radius 25a 0 to prevent evaporation. For n = 8 to n = 10, this radius was increased to 30a 0 . For the particular case of n = 1, simulations of 1,000,000 steps without confinement were run in order to validate the methodological procedure used for larger water clusters. Finally, this process was repeated four times (three times for n = 1) starting from different initial configurations (see Figure 1 in the SI) to improve ergodicity and to characterise the dependence on the arbitrary chosen initial geometry. In a second step, between 900 and 2,000 of the most stable configurations found during the first step were further optimised locally using a conjugated gradient scheme at the DFTB energy level relaxing all degrees of freedom. In order to identify the main families and their lowest energetic isomers, we observed all the structures (for n=1,2,3) or the 500 most stable optimised structures (n>3) following an increasing energetic order. Each time a new pattern was detected by eye, i.e. identification of the first isomer presenting a non-previously observed geometrical pattern, we defined a new family and the current isomer is assumed to be the lowest energy isomer of this family. Such an analysis was made easy due to the fact that isomers belonging to the same family usually presented small energetic variations, leading to a plateau on the energetic profile (Figures 234of the SI) In a third step, two or three of these lowest energy isomers were locally optimised at the DFT level for each stoichiometry. All PTMC and DFTB calculations in this work were carried out with the deMonNano code [START_REF] Rapacioli | deMonNano experiment[END_REF] .

III. RESULTS

The structures and energetics of the most stable isomers of (Py) 2 (H 2 O) n (n = 1 -10) obtained from the final DFT optimisations are first presented and discussed in section III A 1. In addition, structural patterns identified from the stable isomers obtained from PTMC/DFTB simulations (see Figures 2 to 5 in the SI) are discussed in section III A 2. Intermolecular binding energies and competitive dissociation pathways as a function of water cluster size n are presented in section III B. Finally, the evolution of ionisation potentials with n is presented and discussed in section III C.

A. Structures and energetics

Most stable isomers

The most stable structures for (Py) 2 (H 2 O) n optimised at the DFT level are reported in Figures 1 (n = 1 -5) and 2 (n = 6 -10). Each line corresponds to a stoichiometry and structures are labelled by letters according to their DFT energetic order. Although the structures and energies resulting from the full procedure (after the final DFT optimisation) are expected to be the most relevant ones, it is interesting to discuss the main differences with those obtained after the second step (PTMC+local optimisation at the DFTB level). Regarding the geometries, up to 5 water molecules, the final DFT optimisation from the DFTB starting point geometries mostly consists in small hydrogen-oxygen bond rotation. Slight pyrene rotation or different alignment of molecules are also observed for the larger clusters. As can be noticed in Figures 1 and2, for most cases, the lowest energy structure obtained at the DFTB level remains the most stable one or is almost degenerate with the most stable one at the DFT level. The largest difference (about 18 kJ/mol) is observed for (Py) 2 (H 2 O) 7 .

We now focus on the most stable DFT structures and energetics. Their absolute energy values can found in the SI (Table III). In the lowest energy isomers of small clusters (n = 1 -5), all water molecules are almost always located on the side of the pyrene dimer. This leads to a relative translation of the two pyrene monomers and the final arrangement of the pyrene dimer differs slightly from the global minimum of the isolated dimer. [START_REF] Dontot | The Journal of Physical Chemistry A[END_REF] For n = 1, the oxygen atom of the water molecule interacts with one hydrogen atom of the closest pyrene monomer whereas one of its hydrogen atom interacts with an external carbon atom of the other pyrene monomer, which has a negative Mulliken charge (for structure 1-a: -0.26 against -0.08 and -0.01 for central carbon atoms). For n = 2, in structure 2-b, the water dimer organizes as a bridge between the two pyrene monomers with both oxygen atoms interacting with hydrogen atoms of the closest pyrene monomer, and the hydrogen atoms of one water monomer interacting with the π system of one pyrene molecules, the other water monomer being more external. Similar interactions can be seen in structure 2-a, except the "bridge" effect. For n = 3, 4 or 5, the water molecules form a cyclic arrangement, similar to the most stable structures of isolated clusters (triangle, square and pentagon). The same kind of interactions between pyrene and water as described before are found once again. The formation of those cycles suggests that water-water interactions are favored against water-pyrene interactions. This point will be discussed in the next part devoted to energies analysis. We notice that for n = 3, the structure with the water molecules above the pyrene dimer is the most stable for DFT, contrary to n = 1, 2, 4 and 5 (for example, for n = 4, this isomer is 20 kJ/mol above the most stable one). Indeed, as n increases, a competition arises between isomers where the water clusters are located on the side of the dimer, making a bridge between the two pyrene monomers and quoted "side" isomers in the rest of the manuscript, with the isomers where the water cluster interacts with only one pyrene monomer ("face" isomers). In this case, the main interaction occurs between the hydrogen atoms of the water monomers not involved in water-water intermolecular interactions and the π cloud of the closest pyrene monomer.

When the number of water molecules increases, more compact and three-dimensional structures are observed for the water clusters. For n = 6 and 7, the lowest energy structures are similar to the most stable structures of isolated clusters, namely the prism (6-a) and book (6-b) geometries for the hexamer and the prism-like geometry for the heptamer. Cubic (8-a) and cubic-like (9-a) geometries are also found for n = 8 and 9. For the highest stoichiometries (n = 8, 9 and 10), structures where the water cluster grows apart from its most stable structure in the gas phase to open and adopt a more 2D organisation that "wraps" preferentially one pyrene monomer (structures 8-b, 9-b and 10-b on Figure 2) start competing with the lowest energy isomer. This could be interpreted as characteristic sizes for the beginning of solvatation. In order to illustrate the complexity of the PES and the density of low energy isomers, several examples of isomers for n = 6 and 8 are reported in Figure 3. All isomers reported on this figure lie within less than 0.1 eV (the absolute energy values of all isomers can be found in Table IV of the SI). For n = 6, we recognise characteristic structures, such as the prism (a), the book (b, h) or the boat (f) that are among the most stable isomers of (H 2 O) 6 . 63 We also notice smaller cycles such as 4-or 5-carbon cycles for structures (d) and (e) respectively. These structures can be seen as some cage or bag geometries. For n = 8, we also retrieve characteristic structures such as the cube (a, e, g) or the boat (k), structures with smaller carbon cycles (c) or geometries (i, where a prism is recognised), as well as planar arrangements (b, d, f, j, etc.). Interestingly, as in the case of the interaction with one single PAH, the interaction with the pyrene dimer tends to favor energetically more planar structures for the water clusters than in the gas phase, marking the beginning of a solvatation process. 49,50

Isomeric patterns as a function of energy

From the large number of structures optimised at the DFTB level, we derived some characteristic energies associated to the appearance of specific common structural patterns (see Figure 4). These patterns were identified based on the qualitative analysis of the structural features occurring for the different "plateaus" observed on the DFTB energy curves (Figures 2 to 5 in the SI). The first two patterns concern the organisation of the water cluster, the third one that of the pyrene dimer and the fourth one the relative position of the water cluster with respect to the pyrene dimer. The energies of the first occurrences of these patterns -with respect to the most stable isomer in a given stoichiometry -as a function of n are reported in Figure 5. More detailed results, i.e. the energies of the first ten structures for each pattern are reported in Table V of the SI. A first pattern consists in structures with one or several linear chains of 3 to 5 water molecules (isomer "a" in Figure 4). Such structures are found between 6 and 19 kJ/mol above the global minima, depending on the water cluster size. Globally, the larger the water cluster becomes, the less energetically favorable such conformation appears to be, a large gap occur-ring for the water octamer (see Figure 5). The difficulty to form linear chains of water could arise from the fact that this formation leads to a reduced number of stabilising hydrogen bonds and this loss is all the more important as the number of water molecules increases. A second pattern consists in structures where all the water molecules are not gathered in a subcluster (isomer "b" in Figure 4), the most stable ones consisting in a water monomer separated from the remaining water molecules (by at least 2 Å). These structures always lie more than 9-10 kJ/mol above the most stable ones, except for n = 2 (this is linked to the evaporation energy of one water molecule that will be discussed in the next part). A third pattern where the pyrene molecules adopt a T-shaped structure was also observed (isomer "c" in Figure 4). This pattern is found above 8.8 kJ/mol with respect to the global minimum for all stoichiometries except for n = 8, where the cubic geometry seems to participate in the global stability of the structure and to "hold" together the pyrene monomers. This explains the low energy found for this stoichiometry (see orange data points in Figure 5). Finally, in a fourth and slightly less energetically favorable pattern, the water molecules are confined between the two pyrene units in a sandwich-like configuration (isomer "d" in Figure 4). Before n = 6, such structures are not favored and only intermediate structures (between displaced pyrene molecules and a "real" sandwich) are found. Energetically, they lie on average 20 kJ/mol above the global minima. For all stoichiometries, the sandwich pattern always lies above the others. Interestingly, Molina et al. 46 reported the particular stability of a sandwich structure for a cluster constituted of two anthracene molecules and a water tetramer where the latter is inserted between two anthracene molecules with parallel planes and orthogonal main axis. This sandwich-like configuration was determined to lie 2 kcal/mol (8.4 kJ/mol) above the most stable one for (Ant) 2 (H 2 O) 4 , where the water tetramer interacts with one external face of the anthracene dimer. Such sandwich-like configuration, where the water tetramer adopts a linear conformation and is found in between the two pyrene molecules that have relative orientation between parallel and T-shape, appears 15.1 kJ/mol above the most stable one in the case of (Py) 2 (H 2 O) 4 . This difference of structures and energetics between (Ant) 2 (H 2 O) 4 and (Py) 2 (H 2 O) 4 may originate from the difference in the shape of the initial PAH monomer, more linear in the case of anthracene and compact in the case of pyrene. For the higher stoichiometries (n ≥ 8), it becomes more energetically favorable to separate the cluster of water or to get a T-shaped dimer of pyrene than to form a linear chain of water. For n = 4, it is easier to divide the water cluster into subclusters than form a T-shaped dimer of pyrene, contrary to the neighboring stoichiometries. This could be explained by the fact that when we separate such cluster, we retrieve in most cases a water triangle, which is also particularly stable.

B. Binding and evaporation energies

For a given isomer labelled k = a, b, c... of the (Py) 2 (H 2 O) n cluster, several types of binding energies can be computed depending on the dissociation channel taken as a reference. We define (i) E_W k n where the dissociation corresponds to the loss of one water molecule, (ii) E_P k n where the dissociation corresponds to the loss of one pyrene molecule, (iii) E_S k n where the dissociation corresponds to the separation of the cluster into two homogeneous subclusters, i.e. one pyrene dimer and one water cluster and finally (iv) E_M k n where the dissociation corresponds to isolated pyrene and water molecules i.e. the cohesion energy of the total cluster. These energies are given by the following expressions:

E_W k n = -E k Py 2 W n + E a Py 2 W n -1 + E W (6) E_P k n = -E k Py 2 W n + E a PyW n + E Py (7) E_S k n = -E k Py 2 W n + E a Py 2 + E a W n (8) E_M k n = -E k Py 2 W n + 2E Py + nE W (9)
The energy of the initial structure is noted E k Py 2 W n , with k = a corresponding to the most stable isomer. E W and E Py are the energies of the optimised water and pyrene molecules.

E a

PyWn are the energies of the most stable structures of PyW n (n = 1 -6) reported in a previous work. 52 In these structures, the water molecules interact with the π cloud of the pyrene molecules and adopt geometries similar to those in the gas phase (triangle, square, pentagon, prism). E a Py 2 and E a W n are the energies determined for the most stable structures of respectively the isolated pyrene dimer [START_REF] Dontot | The Journal of Physical Chemistry A[END_REF] and the isolated water clusters (n = 1 -10). 52 The geometries adopted by the isolated water clusters starting from n = 3 are: triangle, square, pentagon, prism, prism-like, cubic, cubic-like and prism. Finally, we highlight that these energies are adiabatic dissociation energies (only the lowest energy dissociation channel is considered) and that we do not address here the issue of the presence of some energetic barrier along the dissociation path. The evolution of the cohesion and normalised cohesion energy with respect to the number of water molecules as a function of n are reported in Figure 6. The latter corresponds to: E_M k n -E a Py 2 /n. As can be seen on this figure, the cohesion energy (E_M k n ) of (Py) 2 (H 2 O) n shows a global increase as a function of the water cluster size similar to the evolution of the same quantity for pure water clusters. Indeed, as a first approximation, E_M k n can be estimated by adding the binding energies of pure water clusters to the binding energy of an isolated pyrene dimer (E a Py 2 = 46.7 kJ/mol, computed in this work), cf. the shifted data points in Figure 6. The energy difference between the latter data points and E_M k n corresponds to the energy gain due to the specific arrangement of the water cluster on the pyrene dimer which have been discussed in the previous section. Such stabilisation amounts from 27 to 70 kJ/mol for the most stable isomers (quoted as "a"), 17 to 66 kJ/mol for the "b" isomers and 29 kJ/mol for the "c" isomer. The subtraction of the binding energy of the pyrene dimer allows us to compare the normalised cohesion energies directly with the normalised cohesion energies of pure water clusters:

(E a W n -nE W )/n.
We notice a similar evolution of both graphs: an increase until n = 7 leading to an asymptotic limit at 55.3 kJ/mol for (Py) 2 (H 2 O) 10 . In order to assess the error due to the use of a small basis set on cohesion energies, single point calculations were performed using the same functional and a larger basis set (augcc-pVTZ) for the smallest clusters, namely isomers 1-a and 2b (see Tables VI and VII of the SI). Basis set superposition error (BSSE) was also estimated using the counterpoise method as implemented in Gaussian 16. [START_REF] Boys | [END_REF]81 We show that BSSE represents less than 6% of the cohesion energy with either basis set, and that this energy is overestimated by 15-20% with the double-ζ basis set. Figure 7 presents a comparative graph for the three other evaporation energies, E_W, E_P and E_S, with respect to the number of water molecules n (all the data can be found in Table VIII of the SI). For n = 1, the evaporation energy of a single pyrene is larger than that for one water molecule. The energy required to remove one pyrene unit is 51.5 kJ/mol, which is larger than the intermolecular binding energy within the isolated pyrene dimer (46.7 kJ/mol). This can be understood as in the most stable Py 2 W structure, the attachment of one pyrene unit is reinforced by the presence of the "bridged" water molecule. The energy required to remove one water molecule is 27.2 kJ/mol, a value larger than the binding energy of a water molecule on a single pyrene (22.4 kJ/mol), as can also be inferred from the fact that the water can interact with two instead of one pyrene units. For n = 2, the three dissociation channels are almost isoenergetics, the loss of a single water molecule being slightly less favorable. From n > 2, and with the exception of n = 5 dis- cussed hereafter, removing one water molecule and reshaping the cluster is the most energetically unfavorable process, and this is the case for all the reported isomers. This can be related to the fact that in isolated water clusters the energy required to evaporate one molecule for clusters with more that 3 units is always larger than the binding energy of the pyrene dimer: E ≥ 47.6 kJ/mol for n ≥ 3 (against 46.7 kJ/mol). Turning now to the special case of n = 5, the water molecules seem to adopt a "classic" pentagonal shape (see Figure 1). However, these pentagons are in fact misshapen (cf. Figure 8), with almost all their OOO angles smaller than 107°. Thereby, removing one water molecule from these unstable pentagons and retrieving a square is very favorable, which explains the special behaviour of E_W for this size. For n = 8 -10, the cost for the loss of one water molecule remains almost constant and of the same order as the energy required to evaporate one water molecule from the isolated equivalent water cluster. On the contrary, the energy corresponding to the separation in two clusters drops by almost a factor of two with respect to the values for the smaller sizes. This has to be related to the fact that in these clusters (in particular for n = 8, 9), the water cluster is very stable by itself with cubic-like structures, which are unfavorable to obtain large interaction areas with the pyrene dimer.

For n = 2 -6, the loss of one pyrene molecule is competitive with the dissociation in one water cluster and one pyrene dimer. However, given the energetic data computed for n = 6, the former channel is expected to become favored when increasing the number of water molecules. This can be related to the fact that adding water molecules results in an increase of the stabilising interaction area between the pyrene dimer and the water clusters. One should be cautious in generalising this trend because data are missing for n > 6 and that cluster size effects remain important in this size range (for instance for n = 8).

Beyond the previous discussion on the general trends, we now discuss the effect related to the isomers' nature. For n = 6 and 8, it is more energetically favorable to remove one water molecule from the book geometry ("b" structures) than from the prism or the cube ("a" structures), which can be related to the higher stability of the latter isomers. For n = 7 and 9, we expected that, starting from the prism-like or the cubic-like structure, it would have been quite energetically favorable to retrieve the prism or the cubic geometry by the evaporation of one water molecule due to the peculiar stability of the latter isomers. For n = 7, we indeed observe a small gap between E_W and E_S (around 6 kJ/mol difference). However, for n = 9, E_W is much larger than E_S. This can be related to an increase of the cohesion energy within the water cluster when its size increases. 52 Finally, for n = 10, we can note that there is no enhanced energy increase for the loss of one water molecule (E_W) contrary to n = 4 or 6, probably because no particular stability occurs for n = 10 (and no particular instability for n = 9). In fact, in these structures, the water cluster does not adopt a symmetric geometry, like the prism that can be found for the isolated cluster, and so is easier to break. In conclusion to this part, we expect that in the majority of cases, from an energetic point of view, it is more favorable to dissociate into two homogeneous molecular clusters rather than evaporating one single molecule. From the data reported in Figure 7, we can extrapolate this trend for water cluster sizes exceeding 10 water molecules.

In order to identify particularly stable clusters whose stoichiometries are referred to as "magic numbers", numerical second derivatives of energy with respect to n were calculated:

d 2 E dn 2 = E a Py 2 W n+1 -E a Py 2 W n -E a Py 2 W n -E a Py 2 W n -1 (10)
This calculation allows us to compare the relative stability of the structures: a large value of d 2 E/dn 2 indicates a particular stability for such structure and n can be identified as a magic number. Results are reported in Figure 9 (the energies used are those of the most stable structure in DFT for each stoichiometry, hence "a" structures). With the blue data points, corresponding to isolated water clusters, we recognise the increased stability peaks for n = 4, 6 and 8, which correspond to the known particularly stable structures of square, prism and cubic geometries of water clusters. With the red data points, corresponding to (Py) 2 (H 2 O) n , one peak clearly exceeds the others at n = 4. This peak can be linked to the overall increase of energy observed in Figures 6 and7. This stoichiometry is especially favorable, compared to the others, which could be explained by the fact that the isolated water tetramer presents a particularly stable structure (square shape) and that the lat- ter is particularly adapted to maximize the number of bonds between pyrene and water molecules (see Figure 1). Interestingly, the water hexamer and octamer do not correspond anymore to magic number when they are deposited on the pyrene dimer. Indeed, the symmetric and compact structures of these clusters (prism and cube) reduce the possibility of finding structures maximizing the interaction (i.e. overlap) between the water cluster and the pyrene system. This is less the case when adding or removing one water unit to these clusters (i.e. cluster with n = 5, 7, 9).

C. Ionisation potentials

In order to compute vertical ionisation potentials (IP):

IP = E + Py 2 W n -E Py 2 W n (11)
the energies of cations were computed at the DFTB-CI level of theory for isomers optimised in their neutral charge states (Figures 1 and2). The DFTB-CI calculations were first performed without restricting the charge localised configurations basis, which means that the basis covers all the possible charge localisation patterns (i.e. one per water or pyrene molecule). The corresponding IPs are drawn on Figure 10 ("full basis" red data points) with respect to n for the most stable isomers ("a" structures). We faced convergence issues for clusters with more than 5 water molecules. From the wave function analysis, we found that the charge was localised by more than 97.8 % on the pyrene subsystem (see Table IX in the SI) and this can be related to the fact that the IP of the pyrene molecule (7.4 eV) is lower than that of the water molecule (12.6 eV). [START_REF] Linstrom | NIST Chemistry WebBook, NIST Standard Reference Database Number 69[END_REF] On the basis of this analysis, we also conducted DFTB-CI calculations with a "restricted basis", i.e. removing from the basis (Eq. 4) the configurations where the charge is localised on the water monomer (i.e. the charge can only delocalise over the pyrene dimer subsystem). The corresponding IPs are also reported in Figure 10 ("restricted basis" blue data point). The good agreement between the frac- tions of charges on the monomers of pyrene and on the water molecules computed with the "full" and "restricted" basis can be seen in Table IX in the SI. The values of the absolute energies of the "a" and "b" cationic isomers and those of the corresponding IPs are reported in Table X of the SI for all stoichiometries. The evolution of the IP values as a function of n obtained with the restricted basis for the "a" and "b" isomers are reported in Figure 6 of the SI.

As can be seen on Figure 10, the values of the IPs computed with the "restricted basis" are close to those determined with the "full basis", even identical up to 3 water molecules and with small differences of about 0.03 eV for n = 4 and 5. This suggests again that restricting the charge delocalisation to the pyrene dimer subsystem is a reasonable approximation. This is in line with our previous work where the IPs of individual PAHs adsorbed on a water ice surface were determined. 52 In this study, we showed that constraining the charge to be localised on the PAH unit was the proper way to recover IPs determined at a higher level of theory (MP2) without any charge localisation assumption. It can be seen from Figure 10 that the vertical IP is always larger than that of the isolated dimer of pyrene (black line at 7.13 eV). Let us remind that this pyrene dimer IP is smaller than that of the pyrene monomer (green line at 7.49 eV, computed in this work) due to charge resonance stabilisation, a process that is favored when the overlap between the Highest Occupied Molecular Orbitals (HO-MOs) of the two units increases. The IPs of Py 2 (H 2 O) 1 and Py 2 (H 2 O) 2 are is close to that of the isolated pyrene, which can be understood from the structures reported in Figure 1: the water molecules reduce the overlap of the pyrene molecules and consequently charge resonance stabilisation. The balance between unfavorable interactions (hydrogen of water with the charged pyrene dimer core) with more favorable interactions (oxygen of water with the charged pyrene core) is complex and changing for each structure, making the evolution of the IP very size dependent. For instance, the decrease until n = 8 is followed for n = 9 and 10 by a significant increase of the IP, resulting from the displacement of the dimer of pyrene induced by the water cluster. The case of n = 8 is again particularly interesting as the IP is almost the same as that of the dimer of pyrene. This could be explained by the fact that, for this isomer, the pyrene dimer is hardly disturbed by the water cluster, resulting in an equal share of the charge by the two units (Table SI-IX). The case of n = 2 is similar to that of n = 8: the structure of the pyrene dimer is preserved, which favors an important overlap. However, due to the interaction with the water cluster, the configuration where the charge is localised on the pyrene surrounded by the water molecules is slightly destabilised by the latter. Finally, for "face" isomers ("a" structures for n = 3, 6, 7 and 10), the charge is localised preferentially on the pyrene molecule that does not interact with the water cluster. This could be explained by the unfavorable interaction that occurs in this case (hydrogen of water with π cloud of the charged pyrene molecule). As a conclusion, we can state that we are here in the cluster regime, i.e. non-monotonous evolution of properties with size presenting significant changes with the addition/removal of a single unit.

IV. SUMMARY AND OUTLOOKS

In this work, we explored the PES of the pyrene dimer solvated by water clusters of increasing size (Py) 2 (H 2 O) n (n = 1 -10). The complexity of the PES required an efficient exploration based on PTMC/DFTB algorithms while the lowest energy isomers were further locally optimised at the DFT level of theory. Due to the subtle competition between weak intermolecular interactions, the most stable DFT conformer does not always correspond to the lowest energy DFTB conformer. Note that it is also likely that changing the functional and the basis set might induce some variations in the relative energies. We showed that for the small sizes (n ≤ 5), in the most stable configurations, the interaction with the water cluster induces a relative displacement of the pyrene monomers and the water cluster, that globally preserves its gas phase structure, interacts on the side of the dimer ("side" isomers). An exception occurs for the water trimer whose three hydrogen atoms not involved in the water-water bonds interact with the π electrons of one pyrene monomer ("face" isomer). As for the water dimer, the most stable structure found is an intermediate one between "side" and "face", and is almost degenerated with the "side" isomer. When the size of the water cluster increases, competition is increasingly observed between isomers where the water cluster preserves its gas phase structure and isomers where the pyrene dimer is solvated by a water cluster that has reorganised. Structural patterns as a function of energy were retrieved, showing in particular that "sandwich structures", where the water cluster is confined between two pyrene molecules, are energetically unfavorable. Focusing on the lowest energy isomers for each stoichiometry, we showed that the (normalised) cohesion energy of the water clusters interacting with the pyrene dimer is larger than that of the pure water cluster of the same size and increases as a function of the cluster size (n) until it reaches an asymp-totic limit (55.3 kJ/mol), close to that of pure water clusters. The stabilisation effect is then enhanced for smaller sizes. We also showed that, except from the water pentamer, dissociating the two homogeneous clusters is more energetically favorable than evaporating one water molecule. A question that remains is whether we expect this trend to be extrapolated to larger water clusters. Interestingly, the relative stability of (Py) 2 (H 2 O) 4 was shown and the magic numbers determined for pure water clusters of 6 and 8 units were not observed upon interaction with the pyrene dimer. Finally, vertical ionisation potentials were determined using the DFTB-CI method. In the cationic clusters, the charge was found to be localised on the pyrene subsystem, as expected from the relative IPs of the components. The IPs of (Py) 2 (H 2 O) n (n < 9) lie between the IPs of the pyrene dimer and pyrene monomer. However, no monotonous behavior could be retrieved, the IP being strongly sensitive to the relative orientation of the interacting monomers. The data computed in the present work constitute a step further for the characterisation of gas-grain interactions that are lacking in databases used to develop models of the chemistry of the troposphere and of the ISM. The next step of this work is to search for the most stable structures of the cationic clusters and to retrieve interesting IR features that could be a signature for the presence of PAH and water clusters interactions in the ISM. Such data could be used to interpret further astronomical spectra from the James Webb Space Telescope.

FIG. 1 .

 1 FIG.1. Stable structures found with DFT of (Py) 2 (H 2 O) n for n = 1 -5 from top to bottom. For each stoichiometry, the structures are labelled by letters according to their DFT energetic order. The relative energies (in kJ/mol) with respect to the most stable are reported (black: DFT energies and blue: DFTB energies, both without ZPE).

FIG. 2 .

 2 FIG.2. Stable structures found with DFT of (Py) 2 (H 2 O) n for n = 6 -10 from top to bottom. For each stoichiometry, the structures are labelled by letters according to their DFT energetic order. The relative energies (in kJ/mol) with respect to the most stable are reported (black: DFT energies and blue: DFTB energies, both without ZPE).

FIG. 3 . 8 FIG. 4 .FIG. 5 .

 3845 FIG.3. Examples of low-energy isomers of (Py) 2 (H 2 O) 6 (top, (a-h) with increasing energy) and (Py) 2 (H 2 O) 8 (bottom, (a-l) with increasing energy). These isomers are located within energy ranges of 8.4 kJ/mol and 8.9 kJ/mol (0.09 eV), respectively.

FIG. 6 .

 6 FIG. 6. (a) DFT cohesion energies (in kJ/mol) for the most stable structures of (Py) 2 (H 2 O) n (E_M a n , red crosses) and for pure water clusters, raw values (blue rotated crosses) and shifted by the intermolecular binding energy of the pyrene dimer (+ 46.7 kJ/mol, cyan stars). (b) Normalised DFT cohesion energies (in kJ/mol), for the most stable structures of (Py) 2 (H 2 O) n in red and for (H 2 O) n in blue.

FIG. 7 .

 7 FIG. 7.Comparative graph of E_W (in blue), E_P (in orange) and E_S (in red) of (Py) 2 (H 2 O) n calculated with DFT. Energies for the structures labelled as "a", "b" and "c" in Figures1 and 2are represented by squares, triangles and circles respectively, and the data points for the most stable ones ("a") are joined to guide the eyes.

FIG. 8 .

 8 FIG. 8. Side view of structures 5-a (left) and 5-b (right)

FIG. 9 .

 9 FIG. 9. d 2 E/dn 2 as a function of n for the most stable structures of (Py) 2 (H 2 O) n in red and of (H 2 O) n in blue

  The black and green lines represent respectively the IP of the isolated dimer of pyrene and of the monomer of pyrene.

		7.6		restricted basis
				full basis
		7.5
	Ionisation potential (eV)	7.2 7.3 7.4
		7.1
		7	1	2	3	4	5	6	7	8	9	10
							n
	FIG. 10. DFTB-CI ionisation potential for the most stable isomers
	with respect to n.

ACKNOWLEDGMENTS

The authors thank the computing mesocenter CALMIP ("CALcul en MIdi Pyrénées", UMS CNRS 3667) for generous allocation of computer resources (projects P0059 and P17002). H.L. acknowledges fundings from the European Research Council under the EU's Seventh Framework Programme (FP/2007-2013) ERC-2013-SyG, Grant Agreement no. 610256 NANOCOSMOS. This study has also been partially supported through the EUR grant NanoX n°ANR-17-EURE-0009 in the framework of the "Programme des Investissements d'Avenir".

DATA AVAILABILITY STATEMENT

A Zenodo database containing the .xyz files of the presented sctructures is available : DOI:10.5281/zenodo. 7646613.


Karim Moussa
email: karim.moussa@irt-jules-verne.fr

Eulalie Coevoet

Christian Duriez

Maxime Thieffry

Fabien Claveau

Philippe Chevrel

Stéphane Caro

FEM-based Dynamic Model for Cable-Driven Parallel Robots with Elasticity and Sagging ⋆

Keywords: Cable-Driven Parallel Robots, Dynamic Modeling, Finite Element Method, Cable Model, Beam Theory

Cable-Driven Parallel Robots (CDPRs) are a type of robot that is growing in popularity for different kinds of applications. However, the use of cables instead of rigid links makes the modelling of this robot a complex task, and therefore their trajectory planning and control are challenging. Assumptions such as inelastic, massless and non-sagging cables made when the CDPR is small are no longer valid when the robot becomes large. This paper presents a CDPR dynamic model taking into account cable elasticity and sagging, and its implementation within an open-source framework, named SOFA. Finally, the simulation results are compared to experiments conducted on a suspended CDPR.

Introduction

Cable-Driven Parallel Robots (CDPRs, see Fig. 1) are used to perform operations in a very large workspace: from the assembly of solar panels [START_REF] Pott | Large-scale assembly of solar power plants with parallel cable robots[END_REF] to flight simulators [START_REF] Miermeister | The cablerobot simulator large scale motion platform based on cable robot technology[END_REF], through the development of a large radio telescope [START_REF] Tang | Dimensional design on the six-cable driven parallel manipulator of fast[END_REF]. Simulating this kind of robots is of interest for carrying out potentially dangerous or time-consuming tests, such as the analysis and prediction of their behaviour, e.g. during their design and the synthesis of control laws, or to conduct tests at the safety limits, i.e. failure tests of actuator and cable breakage. It can also make it possible to estimate certain parameters or variables that cannot be measured on a real system. To do so, dynamic models that are both precise and computationally tractable are mandatory.

Several cable models are available in the literature as shown in Fig. 2. A first level of dynamic modeling is assuming massless, straight or non elastic cables [START_REF] Santos | Model predictive control of largedimension cable-driven parallel robots[END_REF][START_REF] Picard | A new control scheme of cable-driven parallel robot balancing between sliding mode and linear feedback[END_REF]. Some works consider elastic, but massless cables, therefore not undergoing Fig. 1: CDPR (CAROCA) located at IRT Jules Verne, from [START_REF] Picard | Pulleys and force sensors influence on payload estimation of cable-driven parallel robots[END_REF].

cable sagging [START_REF] Baklouti | Dynamic and oscillatory motions of cable-driven parallel robots based on a nonlinear cable tension model[END_REF]. This may be of interest to model synthetic cables that have low mass density. In [START_REF] Saadaoui | H∞ control of a planar 3-dof flexible-cable manipulator[END_REF], the authors used the assumed mode approach to model a three degree-of-freedom (DOF) CDPR with sagging and non-elastic cables. The Irvine model [START_REF] Merlet | Some properties of the irvine cable model and their use for the kinematic analysis of cable-driven parallel robots[END_REF] offers the possibility to model faithfully the cable geometry taking into account cable sag and elasticity, but only in the static case. Besides, it neglects the bending stiffness of the cable. It is therefore not directly suitable for dynamic control applications. Recently, a Rayleigh-Ritz cable model that takes into account the variation in cable mass and stiffness was studied in [START_REF] Godbole | Dynamic modeling and adaptive control of a single degree-of-freedom flexible cable-driven parallel robot[END_REF]. However, results are for now limited to a single DOF CDPR. It should be noted that few CDPR simulators exist in the literature based on different cable models. In [START_REF] Okoli | Cable-driven parallel robot simulation using gazebo and ros[END_REF], a massless and inextensible cable model is implemented in Gazebo and ROS to simulate CDPR dynamic behaviour. Multiple nodes of Reissner beam are used in [START_REF] Michelin | Simulation and control with xde and matlab/simulink of a cable-driven parallel robot (cogiro)[END_REF] to model cable elasticity, sag, shear and torsion in the XDE framework. However, some of the physical parameters of the cable such as the Young modulus need to be adapted in order to ensure simulation stability. The simulation time obtained in [START_REF] Michelin | Simulation and control with xde and matlab/simulink of a cable-driven parallel robot (cogiro)[END_REF] is slow for testing and tuning controllers, namely, 10 h to simulate 10 s. CoppeliaSim offers elastic, but straight and massless cables [13]. Other frameworks like AGX Dynamics [14] and MapleSim [15] make it possible to consider elasticity and cable mass with a lumped parameter discretization.

This work proposes a cable model suitable for CDPR control design. The proposed dynamic model considers both cable elasticity and sag. This dynamic model is implemented in SOFA framework [START_REF] Coevoet | Software toolkit for modeling, simulation, and control of soft robots[END_REF], which was developed for finite element modeling and simulation of soft objects. A digital twin of CAROCA, a CDPR located at IRT Jules Verne, is developed in SOFA. The results from the simulation with SOFA are compared to the results from simulators developed with Maplesoft™ MapleSim software.

The paper is organized as follows: Section 2 describes the cable model obtained based on the FEM model of beams in SOFA and its application to a suspended CDPR. Results from SOFA and MapleSim simulations are compared with experimental results in Section 3. Conclusions and future work are drawn in Section 4.

CDPR Dynamic Modeling

This section presents the dynamic model of CDPRs including both cable elasticity and sagging. The presented methodology is intended to be as generic as possible and is valid for all configurations (suspended and fully-constrained), but will be illustrated throughout the paper with a suspended CDPR. The method relies on the finite element method (FEM) implemented within the SOFA framework, a physics-based simulation platform that uses FEM to model, simulate, and control deformable objects.

FEM Model of Beams in SOFA

This method, first presented for SOFA in [START_REF] Duriez | New approaches to catheter navigation for interventional radiology simulation[END_REF], relies on a representation based on three-dimensionnal Timoschenko beam theory [START_REF] Przemieniecki | Theory of matrix structural analysis[END_REF] and a specific corotational formulation to account for large displacements [START_REF] Felippa | A unified formulation of small-strain corotational finite elements: I. theory[END_REF][START_REF] Jourdes | Visual haptic feedback for training of robotic suturing[END_REF]. It is implemented within the BeamAdapter4 plugin of the SOFA framework.

The beam is discretized into N small elements, each of them corresponding to the one depicted in Fig. 3(a). The motion of the flexible beam is then decomposed into two parts, a rigid body motion and a deformation motion, with the assumption that the deformation remains small at the level of each elements. The equation of motion of an object according to Newton's second law is given by: Mq = f e (q, q) + H(q)λ

where q ∈ R n is the vector of generalized coordinates, M ∈ R n×n is the inertia matrix of the system and f e (q, q) : R n × R n → R n gathers external forces.

Finally, H(q) : R n → R n×m is the Jacobian matrix that gives the direction of constraints forces, λ ∈ R m the vector of Lagrange multipliers corresponding to the constraints and m is the number of algebric constraints. In the case of a beam element, q = (q T a , q T b ) T where indices a and b represent the extremities of the beam element as shown in Fig. 3(a), q a and q b contain a positional vector ∈ R 3 and a rotational term ∈ SO(3) defining the pose of each extremity. M is given in [START_REF] Przemieniecki | Theory of matrix structural analysis[END_REF] and depends on the second moments of area I y and I z with respect to the axes y and z respectively, the cable cross section area A, the length of the beam element L, the second polar moment of area J, and the material density ρ. With the corotational assumption, the deformation in the cable elements is small locally and the external force can be calculated with the linear relation:

f e (q) = K(u(q) -u(q)) (2)
where u(q) is the deformed configuration of q, u(q) is the rest configuration of q, both expressed in the local frame, and K is the stiffness matrix given in Eq. ( 3).

K =                                         EA L 0 12EI z L 3 ϕ y 0 0 12EI y L 3 ϕ z 0 0 0 GJ L symmetric 0 0 - 6EI y L 2 ϕ z 0 EI y Φ pz Lϕ z 0 6EI z L 2 ϕ y 0 0 0 EI z Φ py Lϕ y - EA L 0 0 0 0 0 EA L 0 - 12EI z L 3 ϕ y 0 0 0 - 6EI z L 2 ϕ y 0 12EI z L 3 ϕ y 0 0 - 12EI y L 3 ϕ z 0 6EI y L 2 ϕ z 0 0 0 12EI y L 3 ϕ z 0 0 0 - GJ L 0 0 0 0 0 GJ L 0 0 - 6EI y L 2 ϕ z 0 EI y Φ nz Lϕ z 0 0 0 6EI y L 2 ϕ z 0 EI y Φ pz Lϕ z 0 6EI z L 2 ϕ y 0 0 0 EI z Φ ny Lϕ y 0 - 6EI z L 2 ϕ y 0 0 0 EI z Φ py Lϕ y                                         (3)
where

Φ py = 3+ϕ y , Φ ny = 3-ϕ y , ϕ y = 1+ 12EI z GA sy L 2 , Φ pz = 3+ϕ z , Φ nz = 3-ϕ z , ϕ z = 1 + 12EI y GA sz L 2 , G = E 2(1 + µ)
, E is the Young's modulus, µ is the Poisson's ration, and A sy and A sz are the actual surfaces of the cable due to shear along the y and z directions.

In CDPR use cases, the sheering can be neglected and therefore ϕ y = ϕ z = 1.

The cables have a symmetrical cross-section and therefore I y = I z = I. Once the external force f e (q) is calculated in the local frame, it is transformed to the global frame using the rotation matrix of the local frame. The model presented in [START_REF] Przemieniecki | Theory of matrix structural analysis[END_REF] is used to model beams. However, compared to beams, a cable can bend more and is unable to withstand compression force. This difference is mainly due to the second moment of area I that is smaller for a cable. The structure of the cable used for the CAROCA robot is shown in Fig. 3(b), it is a two layers concentric contra-helical cable made of steel where cable strands have opposite lay directions. As shown in [START_REF] Papailiou | On the bending stiffness of transmission line conductors[END_REF], when the cable is straight all the wires are sticking together and act as a one homogeneous body that can bend around the central axis of the cable. The quadratic moment is at its maximum I max . However, when the curvature k of the cable increases the wires start to slip and eventually at the most extreme case each wire bends around its own axis. The quadratic moment is then at its minimum I min defined as:

I min = I c + ns i=1 I smin,i cos β s,i (4)
where I c is the quadratic moment of the steel core around its own axis, I smin,i is the minimum quadratic moment of strand i around its own axis calculated with Eq. ( 4) after replacing strand with wire, β s,i is the lay angle of strand i and n s is the number of strands. Then, I max can be calculated:

I max = I min + a i=1 n i 2 A i r 2 i cos 3 β l,i (5)
where n i is the number of strands in layer i, A i is the cross section area of each strand in layer i, r i is the radius of layer i, β l,i is the lay angle of the strands in layer i and a is the number of layers. Note that I min and I max are independent from the curvature of the cable and its tension force value.

To simplify the model, a constant value for I ∈ [I min , I max ] will be chosen.

Knowing the working region of the cable (e.g. assuming that the curvature will not reach significant values as much as k = 1m -1 ), and based on the identified transition function between EI max and EI min in [START_REF] Papailiou | On the bending stiffness of transmission line conductors[END_REF] shown in Fig. 4, the parameter α = 0.65 is chosen such that:

I = I min + α(I max -I min ) (6

Application to CDPRs

In CDPRs, the cables are attached at one end to the platform using spherical joint-equivalent connectors and are rolled at the other end around the winches. The cables can be oriented by passing through pulleys to change the CDPR configuration from suspended to fully-constrained. In addition to cables, pulleys and winches have an influence on the full model of the CDPR. Their impact varies depending on their dimensions and the friction they induce. To get the full CDPR model, each cable is modeled by Eq. ( 1) where the constraints represent the connection between the cable and the platform, and its contact with the pulley and the winch. The spherical joint does not transmit torque, therefore only position constraints are applied. A 6-DOF, 2-DOF, and 1-DOF dynamic models for the platform, the pulley, and the winch, respectively are implemented. The cable/pulley, cable/winch and gearhead frictions are neglected. In the following section, the different simulations are compared with the experimental results.

Comparison with Existing Simulators

Modeled robot

For this first trial, the implemented model in SOFA described in Sec. 2.1 is compared to other state-of-the-art models ("Rope" and "Cable" described in Sec. 3.2) of the simulation software MapleSim, and to real-robot measurements from the CAROCA experimental platform (see Fig. 1). It is a suspended CDPR with 8 cables allowing 6-DOF for the platform. Steel cables (see Fig. 3(b)) are used, they can withstand significant tensions (10 kN), but can present significant sags given their mass. The working area of the robot is 7 m × 4 m × 3 m. The moving platform's dimensions are 1.5 m × 1.5 m × 1 m, and it weighs 366 kg.

MapleSim

The MapleSim [15] software provides some tools to model cable-driven robots including cable models, pulleys, and winches. This framework provides three possible models : rope, chain and cable. Rope is the simplest one where only the elasticity is considered, and cable is the most accurate where elasticity and sagging are considered. The rope model is made up of two massless nodes connected by a spring and a damper. The cable model relies on [START_REF] Mitiguy | Determination of spring constants for modeling flexible beams[END_REF], where the cable is discretized into N equidistant nodes, being connected with each other by springs and dampers. Moreover, the cable mass is spread along the nodes. These first trials proposed in this paper are dedicated to these two models, they will serve as a reference to compare with SOFA's implemented cable model's results.

Definition of scenario

The pick-and-place trajectory used to compare the simulators with respect to the experimental results is depicted in Figs. 5 and6. Using an ideal cable model (massless and rigid), an inverse geometric model (IGM) and an inverse kinematic model (IKM) allow to find the desired joint space, i.e. the winches, angular position θ d and velocity θd , as shown in Fig. 7. Then, a PD controller is used to calculate the winch torques τ w and track the reference trajectory. Controller gains in simulators are tuned slightly different than the robot as they must compensate for the friction [START_REF] Picard | A new control scheme of cable-driven parallel robot balancing between sliding mode and linear feedback[END_REF].

Simulation Results

The parameters needed to simulate the models (BeamAdapter in SOFA, rope and cable in MapleSim) and their values are shown in Table 1. C a is the cable's axial damping and d r is its bending damping ratio. As said before, the sheer is negligible in CDPR applications and so the precision of Poisson's ratio's value has no impact on the simulation. The first ten seconds of simulations have been cut on the results depicted in Figs. 8 to 11 ; they are associated to the initialisation of the robot, where the platform is kept at the initial target Cartesian position. This transitory behaviour before reaching the static equilibrium is not of interest. Then, the 30s trajectory shown in Fig. 5 is applied.

The measured cable tensions at the platform from both frameworks are compared to the experimental measures of the force sensors located between the cables and the moving-platform (MP). The resulting plots for cables 6 and 7 are shown in Fig. 8. The remaining six cables behave in a similar manner. It appears that the plots of both models in MapleSim are almost identical. This is probably due to the weight of the cable being insufficient to make a difference under the current cable tension between the massless (rope) and nonmassless (cable) models, and therefore the cables are almost straight. A small oscillation for the rope model is noticed, the PD controller gains might be high for this model. The plots of SOFA's model match those of MapleSim. However, these three models have the following divergence with the experimental results. On one side, tensions in cables 2, 4, 6, and 8 match the experiment measurements. On the other side, they converge to a value different from the measured one in cables 1, 3, 5, and 7. This causes high Relative Root Mean Squared Errors (RRMSE) of some simulated cables tensions with respect to the the experiment (> 20%) as can be seen in Fig. 9. The order of the cables and winches is given in Fig. 6. The SOFA model has slightly lower RRMSE mean value compared to MapleSim. As a first explanation, because of the redundancy of actuation, multiple set of cables tensions can result in the same platform position. A common distribution algorithm should be implemented in the three simulators and the experimental setup to handle this aspect. As for the winches, the angular position of winch 6 is shown in Fig. 10. The plots of the remaining 7 being similar. Contrary to the cables tensions, the plots of the three simulated models match the experimental measurements. This is due to the controller being implemented in the joint space and not considering the cable tension. This can be seen in Fig. 11 where all RRMSE are very small (< 5%) with respect to the the experiment, with a slight advantage in mean value for MapleSim models.

The simulation time for every model varies and is given in Table 2. The solver used in MapleSim is a variable step solver (Rosenbrock). In SOFA, a fixed step solver was used (implicit Euler) with a time step of dt equal to 0.005s.

Ultimately, the SOFA model presents a computational burden far less important than the Maplesim cable model, while it permits to simulate more internal degrees of freedom; 294 nodes in SOFA versus 7 in MapleSim.

As the mass of the cable increases (longer cables or heavier material), the applied motor torque must be increased to keep the cables as straight as possible. Nevertheless, this might not be possible because of security or power consumption limit. To verify the capacity of the proposed simulator to model the sagging, another simulation was run where the winch 1 was locked at its initial position.

The MP pose and the shape of the cables at the end of the simulation are illustrated in Fig. 12.

Fig. 12: The pose of the platform and the state of the cables (sagging) at t = 30s when winch 1 is locked.

The rolled cable on winch 1 has an important sagging. The rolled cable on winch 5 shows also non negligible sagging due to the orientation of the platform that makes the attach point closer than the planed one. These very first results demonstrate that the SOFA model makes it possible to simulate important cable sagging. Further tests will be conducted on the robot to evaluate the precision of the simulation's sag.

The attached video in SOFA shows the vibrations of the cables but not much sagging in the first simulation. The sagging appears clearly in the second simulation. It can be found online at: https://youtu.be/uUwP0CZeMqo.

Conclusion and Perspectives

This paper dealt with the modeling and simulation of cable-driven parallel robots (CDPRs). Specifically, the model sought must account for the effects of cable sag and elasticity. Although ambitious in terms of the realism targeted, the complexity of the model must remain limited. The approach followed is that of the use of three-dimensional beam theory [START_REF] Przemieniecki | Theory of matrix structural analysis[END_REF] and corotational method [START_REF] Felippa | A unified formulation of small-strain corotational finite elements: I. theory[END_REF]. The model was implemented using finite element method within a physics-based simulation platform, the SOFA framework. To evaluate the proposed model, the simulation results were compared with real data obtained on the CAROCA CDPR experimental platform, but also with other simulation results obtained from a model implemented using MapleSim software and the MapleSim Ropes and Pulleys Library. The cable tensions, as well as the angular positions of the winches were compared. The proposed model was implemented in SOFA and proved to be as accurate as the model implemented in MapleSim. It should be noted that it is faster than the MapleSim cable model, but slower than the MapleSim massless rope model. Therefore, it is more advantageous in applications where sag is greater.

The perspectives relate to the use of other cable models [START_REF] Adagolodjo | Coupling numerical deformable models in global and reduced coordinates for the simulation of the direct and the inverse kinematics of soft robots[END_REF][START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF], still within the SOFA framework, to model CDPRs. In particular, the theory of Cosserat [START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF] will be considered to model the cables. Finally, the code of the simulator will be available in open source.

4 Fig. 2 :

 42 Fig. 2: Illustration of different cable models.

  Cable structure used for the CAROCA robot.

Fig. 3 :

 3 Fig. 3: Cable discretization.

Fig. 4 :

 4 Fig. 4: Theoretical evolution of bending stiffness EI as function of the curvature k for different tension T values, based on [21].

Fig. 5 :

 5 Fig. 5: Reference Cartesian trajectory of the robot platform.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Viewing environment in SOFA with added winch numbering w i and the reference Cartesian trajectory of the platform from P 1 to P 4 .

7 Fig. 8 :

 78 Fig. 8: Comparison between simulated and measured cable tensions.

Fig. 9 :

 9 Fig.9: RRMSE (%) for each cable tension.

Fig. 10 :

 10 Fig. 10: Angular position of winch 6.

Fig. 11 :

 11 Fig. 11: RRMSE (%) for each winch angular position.

Table 1 :

 1 Tuning parameters associated to each cable model.

	Models Parameters	E (×10 11 Pa)	µ	A (×10 -5 m 2 )	ρ (kg m -3 )	Ca (×10 9 N s m -1 )	dr α
	BeamAdapter	1.022	0.28 2.1488	7850	-	-0.65
	Rope	1.022	-	2.1488	-	22.8	-	-
	Cable	1.022	-	2.1488	7850	22.8	0.12 -

Table 2 :

 2 Simulation time for different models.

		MapleSim rope MapleSim cable Sofa model
	Time	20 min	12 h	45 min

https://github.com/sofa-framework/BeamAdapter

⋆ This work is supported by both IRT Jules Verne in the framework of the PERFORM program and the EquipEx+ TIRREX project, grant ANR-21-ESRE-0015.