Ahmad W Huran

Nadia Ben Amor

Stefano Evangelisti
email: stefano.evangelisti@irsamc.ups-tlse.fr

Sophie Hoyau

Thierry Leininger

Véronique Brumas

The Electronic Structure of Beryllium Chains

We present an ab initio theoretical study of quasi one-dimensional beryllium chains, Be N , from an electronic structure perspective for N=3, 4,..., 12. In particular, linear and cyclic systems have been compared by using high-quality Coupled-Cluster formalism. Both linear and cyclic species have been found to be local minima on the corresponding Potential-Energy Surface, for all the considered values of N . The linear geometry is the most stable one only in the case of Be 4 . Several indicators (energy gap, position spread tensor, locality of the molecular orbitals) clearly show that both linear and cyclic 1D structures, unlike 3D bulk beryllium, have a covalent insulating nature.

Introduction

The interest in atomically-thin linear surface-deposited structures has experienced huge growth in the last years. This growth was driven, for the most part, by the possibility of synthesizing one-dimensional electron systems at stepped surfaces. [START_REF] Emberly | Electron standing-wave formation in atomic wires[END_REF][START_REF] Bahn | Chain formation of metal atoms[END_REF][START_REF] Himpsel | One-dimensional electronic states at surfaces[END_REF][START_REF] Amorim | Short linear atomic chains in copper nanowires[END_REF][START_REF] Crain | End states in one-dimensional atom chains[END_REF] These systems exhibit a common characteristic, namely the presence of edge orbitals, which correspond to what Tamm and Shockley introduced long ago as "surface states" in a one-electron picture. [START_REF] Davison | Basic Theory of Surface States[END_REF] The edge effects were predicted in graphene-based nanoribbons and nanodots, [START_REF] Hod | Edge effects in finite elongated graphene nanoribbons[END_REF][START_REF] Hod | Half-metallic graphene nanodots: A comprehensive first-principles theoretical study[END_REF] that can give rise to a wide spectrum of interesting electric and magnetic phenomena. [START_REF] Son | Half-metallic graphene nanoribbons[END_REF] These states, in one-dimensional setting, are referred to as "end states", the presence of which has been observed, for e.g., in one-dimensional chains of self-assembled gold atoms on silicon surfaces. [START_REF] Crain | End states in one-dimensional atom chains[END_REF] Ab initio theoretical treatment of linear beryllium chains, close to equilibrium configurations, predicted the existence of two edge orbitals localized at the extremities of the chains. [START_REF] Vetere | A theoretical study of linear beryllium chains: full configuration interaction[END_REF] These edge orbitals are partially filled, and lead to two low-lying quasi-degenerate states, namely, 1 Σ g and 3 Σ u . The energy difference between the singlet ground state and the triplet state presents an exponential decay with increased length of the atomic chain.

Recently, a multitude of interesting effects and features exhibited by beryllium in different chemical environments have been theoretically predicted by Yáñez and coworkers.

The formation of the so-called beryllium bonds between a Lewis base and a beryllium Lewis acid [START_REF] Yáñez | Beryllium bonds, do they exist?[END_REF] was found to act as a modulator of intra-molecular interactions, [START_REF] Mó | Modulating the strength of hydrogen bonds through beryllium bonds[END_REF] and as the driver of spontaneous radical formation. [START_REF] Brea | Exergonic and spontaneous production of radicals through beryllium bonds[END_REF] Be-Be bonds on the other hand were proposed as possible moiety to trap anions utilizing the electron deficiency of the bond in what is known as beryllium-based anion sponges. [START_REF] Brea | Beryllium-based anion sponges: close relatives of proton sponges[END_REF][START_REF] Montero-Campillo | Beryllium-based fluorenes as efficient anion sponges[END_REF] A new kind of Be-Be bonding was also predicted, namely a one electron Be-Be bond of a hybrid character dominated by the contribution of s orbitals. [START_REF] Brea | On the existence of intramolecular one-electron Be-Be bonds[END_REF] Versatility of Be-based species extends well beyond small molecular structures and has spurred interest in the context of linear surface-deposited atomic chains, as it serves as a prototypical subject for studying such chains. Linear beryllium chains have been extensively investigated by our group. [START_REF] Vetere | A theoretical study of linear beryllium chains: full configuration interaction[END_REF][START_REF] Monari | Theoretical study of Be N linear chains: optimized geometries and harmonic Frequencies[END_REF][START_REF] Helal | Electronic bistability in linear beryllium chains[END_REF][START_REF] Pastore | A theoretical study of Be N linear chains: variational and perturbative approaches[END_REF][START_REF] Pastore | Mixed valence character of anionic linear beryllium chains: a CAS-SCF and MR-CI study[END_REF][START_REF] Evangelisti | Beryllium chains interacting with graphene nanoislands: from anti-ferromagnetic to ferromagnetic ground state[END_REF] Because of the presence of unpaired electrons, unsaturated linear chains, would be certainly highly reactive. However, the terminal edge orbitals in a Be N chain can be saturated by means of hydrogen atoms, finally obtaining a linear HBe N H chain. Indeed, HBe 2 H clusters have been observed in beryllium vapors by Tagues and Andrews. [START_REF] Tague | Reactions of beryllium atoms with hydrogen. Matrix infrared spectra of novel product molecules[END_REF] Another possibility of getting a saturated species is by wrapping the chain into a ring, thus creating a topologically closed structure. [START_REF] Fertitta | Investigation of metal-insulator-like transition through the ab initio density matrix renormalization group approach[END_REF][START_REF] Fertitta | On the calculation of complete dissociation curves of closed-shell pseudo-onedimensional systems via the complete active space method of increments[END_REF][START_REF] Koch | Calculation of the static and dynamical correlation energy of pseudo-one-dimensional beryllium systems via a many-body expansion[END_REF] Notice that, although a ring is, strictly speaking, a 2-Dimensional (2D) system, its local topology has a close similarity to the 1-Dimensional (1D) chains. For this reason, both structures are considered as quasi-1D ones.

In this work, we studied the electronic structure of open (linear) and closed (cyclic) Be N clusters, where N in any integer number ranging from 3 to 12. We optimized the geometry of both linear and cyclic chains by using a high-quality Coupled-Cluster (CC) formalism. on the Potential Energy Surface (PES) for each structure. An important point that is addressed in this work concerns the nature of the Be-Be bond at equilibrium distance. Bulk 3D beryllium is a metal, and it is often said that the chains are quasi-1D metallic systems.

Such an opinion can even be found, for example, in an article of one of the authors of the present paper. [START_REF] Bendazzoli | Full configurationinteraction study of the metal-insulator transition in model systems[END_REF] In order to investigate this point, we studied the behavior of the Total-Position Spread (TPS) tensor as a function of the system size. Indeed, it has been shown by Resta and co-workers that the value of the TPS per-electron, or "Localization Tensor", [START_REF] Resta | Electron localization in the insulating state[END_REF][START_REF] Resta | Why are insulators insulating and metals conducting?[END_REF][START_REF] Resta | Kohns theory of the insulating state: A quantum-chemistry viewpoint[END_REF][START_REF] Resta | Polarization Fluctuations in Insulators and Metals: New and Old Theories Merge[END_REF] tends toward a finite quantity if the system is an insulator, while it diverges in the case of metals. For this reason, we computed the eigenvalues of the spread tensor for both types of chains.

The body of the article is structured as follows. Firstly, we consider beryllium chains from an electronic structure perspective. The case of linear and cyclic geometries will be compared. Then the theoretical framework of the investigation is given, and the details concerning the used computational techniques are described. The results are presented and discussed thereafter: the symmetry and the canonical and localized orbitals of these systems are considered, and the Highest Occupied Molecular Orbital-Lowest Unoccupied Molecular Orbital (HOMO-LUMO) energy gaps discussed. Subsequently, results on the optimized geometries, harmonic frequencies, total energies and their extrapolation in the limit of infinitelength chains as well as the behavior of the TPS tensor are presented. Finally, we close by summarizing our main results and give concluding remarks.

Electronic Structure

Beryllium has four electrons in the isolated atom, and a closed-shell electronic configuration, 1s 2 2s 2 . This leads to an electronic structure somewhat resembling that of a rare gas.

Beryllium exhibits an Ionization Potential (IP) that is larger than the IP of its period successor, a relatively infrequent behavior in the Periodic Table, that beryllium shares with other Alkaline-Earth metals.

With only two valence electrons, beryllium is found in molecular moieties of co-linear geometrical configurations, for e.g., the well known case of beryllium hydride (BeH 2 ), and the less studied, Be N H 2 linear structures. [START_REF] Tague | Reactions of beryllium atoms with hydrogen. Matrix infrared spectra of novel product molecules[END_REF] Consider a linear arrangement of N beryllium atoms that, for simplicity, defines the z axis of our reference Cartesian coordinate system. We assume a D ∞h as the symmetry point group of the chain, to be verified, and label the orbitals by g and u according to their character under the inversion symmetry operation.

For a qualitative description of the electronic structure of the systems at hand it is enough to consider σ valence orbitals, as the role of π valence orbitals can be neglected due to their relatively high energies. For a chain made of N Be atoms, the total number of valence electrons is 2N . 2N -2 of these electrons are accommodated in the N -1 σ-bonds formed from the inner sp z hybrids along the chain. The outer sp z hybrids combine to form two edge orbitals a σ g , and a σ u hosting the remaining two valence electrons. In terms of these two symmetry-adapted orbitals, the lowest energy configurations of the two edge electrons are 1 √ 2 (σ g σg -σ u σu ), and 1 √ 2 (σ g σu + σ u σg ), yielding a singlet 1 Σ g and a triplet 3 Σ u state, respectively. Due to the large separation between the edge sp z hybrids, the aforementioned states are quasi-degenerate.

We consider now the total energy of both linear and cyclic systems. In the case of linear systems, we have shown [START_REF] Monari | Theoretical study of Be N linear chains: optimized geometries and harmonic Frequencies[END_REF][START_REF] Helal | Electronic bistability in linear beryllium chains[END_REF][START_REF] Pastore | A theoretical study of Be N linear chains: variational and perturbative approaches[END_REF][START_REF] Pastore | Mixed valence character of anionic linear beryllium chains: a CAS-SCF and MR-CI study[END_REF] that the total energy can be expressed as:

E linear (N ) = E ST (N ) = (N -1)E Be-Be + E Be• (1)
where E ST indicates the mean singlet-triplet energy, E Be-Be is the energy of a Be atom having two Be-Be bond, while E Be• is the energy associated to the terminal edge electrons.

In practice, for large systems, the singlet and the triplet are degenerate, so one can extract the total energy from single-reference calculations on the triplet. Notice that by using this definition the energies associated to the inner core electrons are absorbed in the bond or edge energy definition. The total energy per Be atom can be rewritten as

E linear (N )/N = E Be-Be + 1 N (E Be• -E Be-Be ) (2)
This means that the energy per atom is a linear function of 1/N , N being the number of beryllium atoms. The limit for N → ∞ gives the energy per atom of the infinite chain, while the (formal) extrapolation of the curve to the value N = 1 gives the value of the E Be• parameter.

In the case of cyclic systems, on the other hand, no terminal effects are present. In turn, in this case, the presence of a strain energy must be considered. This is due to the fact that the two bonds on each beryllium atom, -Be-, form an angle smaller than π, which is the ideal one for the sp hybridization of beryllium atoms. In the case of a regular ring containing N atoms, the difference between the ideal angle π and the actual one is given by θ = 2π/N .

For small deformations from the ideal local geometry, the strain energy is expected to be a quadratic function of the deformation. In other words, for small values of θ (which means in the limit of long chains), the strain energy is expected to be a quadratic function of θ.

Therefore, the total energy per beryllium atom will become

E cyclic (N )/N = E Be-Be + 1 N 2 E Be (3)
In this expression, E Be-Be is the Be-Be bond energy, as in equation (2), while E Be represents the strain energy. Therefore, for cyclic chains, the total energy per atom is expected to be a quadratic function of 1/N . As for the case of linear systems, limit for N → ∞ gives the energy per atom of the infinite chain, while the line slope gives the E Be parameter. If the present analysis is correct, the parameter E Be-Be should be exactly the same in the two energy expressions, eq.s (2) and (3). Notice that in these equations the parameters E Be-Be and E Be• include also the large core energy associated to the 1s electrons.

In this work we are particularly interested in comparing linear and cyclic systems up to relatively long chains (N =12). This means that the use of high-quality size-consistent methods is particularly appropriate. With this respect, the Coupled-Cluster (CC) formalism is certainly the best tool nowadays available in Quantum Chemistry. However, CC formalism is presently restricted essentially to the treatment of closed-shell systems. This is not a problem for the treatment of cyclic chains, due to the single-reference nature of their ground state. Linear chains, on the other hand, have an open-shell singlet ground state, whose description is highly problematic at CC level. In the course of our previous investigations, however, we have shown that the single-reference first excited triplet is quasi degenerate with the ground state. Their energy split, in fact, decreases exponentially as a function of the number of atoms in the chain, and can be safely neglected in energetic considerations except for the shortest chains for which the triplet is preferred. In the present investigation, in order to apply a high-quality Coupled-Cluster formalism to the study of the systems, we made the choice of studying the singlet ground state of cyclic chains and the triplet first-excited state in the case of linear chains.

Computational Details

Quasi-1D systems are the linear and cyclic chains. The linear chains have a D ∞h symmetry, while cyclic ones have D nh , where N is the number of Be atoms. These symmetries have been checked a posteriori by means of a frequency calculation, as described in the following section. However, since most ab initio packages can only handle abelian groups, the actual calculations were carried out considering the largest abelian subgroup corresponding to the system at hand. This is the D 2h subgroup for linear and even cyclic chains, while odd cyclic chains were studied in the C 2v subgroup.

Crystalline metallic beryllium has a Hexagonal Close-Packing structure. In order to confirm the non-metallic nature of the Be-Be bonds in 1D structures, it is instructive to consider also three-dimensional (3D) clusters. Since Hexagonal Close-Packing (HCP) arrangements of weakly interacting atoms are often quasi-degenerate with Face-Centered Cubic (FCC) structures, we considered highly symmetric clusters as fragments of these crystal structures.

The structures we chose are the tetrahedral Be 4 structure (T d symmetry), and two isomers of Be 13 : the icosahedral structure (I h symmetry, which can be seen as a fragment of a FCC crystal), and a D 3h structure which is the seed of a HCP crystal. At HF level, Be 4 has a closed-shell ground state, and do not create particular difficulties. The two forms of Be 13 , on the other hand, have open-shell ground states that require some caution. In particular FCC Be 13 , has a triply-degenerate HOMO manifold hosting four electrons, and it is not possible to have a neutral non broken-symmetry wavefunction. For this reason, we also considered the two ionic undistorted structures that are closest to the neutral one, i.e. the quartet Be + 13 and the singlet Be 2- 13 . FCC Be 13 , on the other hand, has a doubly degenerate HOMO pair containing two electrons, and we considered therefore the triplet lowest state.

The ANO-type basis set optimized by Widmark and co-workers was chosen for the atomic orbitals, [START_REF] Widmark | Density matrix averaged atomic natural orbital (ANO) basis sets for correlated molecular wave functions: I. First row atoms[END_REF] adopting the largest contraction recommended by the authors. We performed Restricted Hartree-Fock Self Consistent Field (RHF-SCF) and Complete Active Space SCF (CAS-SCF) calculations for the singlet ground state of cyclic system, and the triplet lowest state of diradical linear chains. Geometries have also been optimized at Coupled-Cluster level restricted to Single and Double excitations (CCSD). Finally, the contribution of triple excitations was taken into account at a perturbative level by performing a CCSD(T) calculation at the CCSD optimized geometry. The harmonic frequencies were evaluated numerically at the SCF and CCSD levels. All CC calculations were performed with doubly occupied RHF-frozen 1s orbitals.

Finally, we note that all the calculations reported in this work were carried out using the MOLPRO package. [START_REF] Werner | version 2015.1, a package of ab initio programs[END_REF] In particular, we performed RHF-SCF, CAS-SCF, CCSD, [START_REF] Hampel | A comparison of the efficiency and accuracy of the quadratic configuration interaction (QCISD), coupled cluster (CCSD), and Brueckner coupled cluster (BCCD) methods[END_REF] and CCSD(T) [START_REF] Watts | Coupled-cluster methods with noniterative triple excitations for restricted open-shell Hartree-Fock and other general single determinant reference functions. Energies and analytical gradients[END_REF] calculations with the corresponding MOLPRO sections. Geometry optimizations were obtained through the rational-function approach and the geometry DIIS algorithm. [START_REF] Eckert | Ab initio geometry optimization for large molecules[END_REF]

Results and Discussion

We discuss here the behavior of the singlet ground state of cyclic Be N chains and the first excited triplet state of linear chains containing N Be atoms for values of N ranging from 3 to 12. Geometries were optimized at CCSD level (see Tables S1 andS2 in the Supplementary Information), and single-point calculations at the CCSD optimized geometries were performed at CCSD(T) level (Tables S3, and S4).

Linear-Chain Symmetry

The symmetry group of the linear chains Be N is D ∞h . The different equilibrium-bond-lengths of linear structures coincide up to few percents. The symmetry of the valence occupied orbitals for even chains (N = 2m) are, in ascending energy order, σ g , σ u , σ g ,... The two singly occupied molecular orbital (SOMO) have σ u and σ g symmetry and become quasi-degenerate orbitals for long chains. At CAS-SCF level, the two orbitals have very close occupations in the ground state singlet, 1 two adjacent atoms only. As usual, we assume a β value < 0 for the hopping integral between two adjacent orbitals (its actual value is arbitrary), while the monocentric α integral does not play any role, and can be put equal to zero. All the resulting orbitals will be doubly occupied, so this model is not able to describe any excitation. Nevertheless, it can give some useful information on the Highest Occupied Molecular Orbital (HOMO) structure. In fact, in the case of even-N cyclic chains, the HOMO is a non degenerate orbital. In the Hückel model, its energy has a constant value of -2β. The HOMO in odd cyclic chains is given by a pair of degenerate orbitals, whose Hückel energy is given by -2β cos(π/N ). The two different behaviors are illustrated in Figure 1, where the orbital energies for the two cases of C 5 H 5 and C 6 H 6 are reported.

Canonical and local orbitals

In Figure 2a, some of the ROHF valence orbitals are shown, for a linear chain of intermediate length (Be 6 ) at the equilibrium geometry. The five doubly occupied canonical orbitals are different bonding combinations of Be-Be sp z hybrid orbitals. The singly occupied edge orbitals are hybrid sp z orbitals pointing towards the outer part of the chain.

Local orbitals have been obtained by a unitary transformation performed on the ROHF canonical orbitals for the case of N = 6, with the DOLO code. [START_REF] Maynau | Direct generation of local orbitals for multireference treatment and subsequent uses for the calculation of the correlation energy[END_REF] In Figure 3, the valence local orbitals of linear Be 6 are shown: there are five doubly occupied valence bonding orbitals, together with the corresponding empty antibonding orbitals.

Between these two sets, two singly occupied edge orbitals are found. They come from the linear combinations of the σ g and σ u quasi-degenerate orbitals that are located at the Fermi level of the system.

In Figure 2b, RHF valence orbitals are shown for the cyclic (Be 6 ) chain at its equilibrium geometry. In Figure 4, the corresponding local orbitals are illustrated. In this case, there are N sets of orbitals that are equivalent because of symmetry reasons. Therefore, only four types of different valence orbitals are found: bonding and antibonding, both having a σ character, as well as in-plane (σ) and out-of-plane (π) non-bonding orbitals. For this reason, only one doubly occupied bonding orbitals and its corresponding empty antibonding orbital have been reported. The in-plane non-bonding orbital is also illustrated in the figure.

In Figure 2, the orbital energies of the valence doubly occupied or singly occupied edge orbitals are also reported. For both linear and cyclic chains, the highest occupied orbitals have negative energies. In the linear Be 6 chain, the SOMO energies are closed to the highest doubly occupied orbital one.

HOMO-LUMO Energy Gap

Although the HF calculations lack the correlation contribution, the HOMO-LUMO energy gap computed at HF level is an indicator of the metallic character of an extended system.

For linear chains, HOMO correspond to the highest SOMO (the two SOMO become quasidegenerate for N > 8). In Figure 5 (see also Table 1), the HOMO-LUMO gap is illustrated as a function of N for the two types of geometries. In both cases, the gap is roughly constant as a function of the system size, and no sign of metallicity is shown. The linear gap is almost constant as a function of the size, and converges toward a limit close to 0.26 hartree for long chains. As discussed previously, the cyclic gap shows an alternating behavior between even and odd structures. Because of this oscillating trend, its large-size limit is more difficult to evaluate, but should be close to 0.23 hartree. The difference between the two limit values are clearly due to the presence of the edge orbitals in the linear systems. In both cases, these results are consistent with an infinite limit of the gap well different from zero. For a comparison, the gap is 0.344943 hartree for the isolated Be atom, and becomes 0.260674 hartree for Be 4 . In the case of Be 13 , the gap goes down to 0.136318 (0.170492) hartree for the HCP (FCC) symmetry structure of the closed-shell neutral (dianion) compound. Although the comparison between gaps of systems having different charges should be taken with some caution, these values show a clear descending trend as a function of the system size, as one would expect from a metallic system.

Optimized Geometries and Harmonic Frequencies

In Tables S1 andS2, the optimized bond lengths are reported for the linear and cyclic structures, respectively. All bond distances of chain structures are systematically shorter than the value in the bulk, which is about 4.195 bohr. [START_REF] Kittel | Introduction to solid state physics[END_REF] This fact can be understood if we consider that each electron is shared among a much large number of "bonds" in the bulk than in the chains. From these tables, it appears also that the bond lengths are weakly correlated with the position of the bond, or the chain length. The terminal bonds of open chains are an exception to this behavior, since their length are considerably longer than all the other bonds. It is interesting to see that the same type of phenomenon has been described in crystalline beryllium, where the inter-planar spacing between the first and the second planes at the (0001) crystal surface has been found to be expanded by about 0.2 bohr. [START_REF] Davis | Anomalous interplanar expansion at the (0001) surface of Be[END_REF] In cyclic Be N , all the bonds are equivalent. The optimized bond length are reported in Table S2 and Figure 6. For small systems, the values of even cyclic chains are significantly shorter than those corresponding to odd chains. The bond lengths decreases when N increases for odd values of N and, on the contrary, increases for even N in the case of the short chains (up to N = 7). This can be related to the different behavior of even and odd closed chains, as discussed in the previous section. For larger values of N , on the other hand, the two curves are barely distinguishable, and they tend to the common limit of an infinite system for large values of N . The asymptotic value can be estimated to be about 4.04 bohr.

In Table S1 and Figure 7 the distances for linear chains are reported. In our notation the numbering of the bonds begins from the center and ends at the edge of the chain. Figure 7 gives a clear depiction of the general trend in bond lengths, where each curve describes an N membered chain (N = 3, 4, ..., 12). The external bond ((N -1)/2 bond number for odd systems and N/2 for even ones) has a value of about 4.087 bohr. In each chain, this bond is the longest one, being about 0.056 bohr longer than the previous one ((N -3)/2 for odd chains and N/2 -1 for even ones). We note that the outer Be atoms are surrounded by three valence electrons compared to four electrons for the inner atoms, explaining the relatively long external bonds. The preceding ((N -3)/2 for odd chains and N/2 -1 for even ones) bond is the shortest one in each chain, reaching around 4.031 bohr. Then, after slight increments in bond lengths, saturation is reached at a value ∼ 4.044 bohr. We note that bond lengths in chains of different sizes are reduced for the shorter chains, with the exception of Be 3 . Finally, we notice that consecutive odd (N + 1) and even (N ) chains have very similar bond lengths when N increases.

It was found that the linear geometry is the most stable one only in the case of Be 4 .

Concerning the linear chains, we computed harmonic frequencies at CAS(2/2) level, for the singlet and triplet states for short chains. It turns out that singlet and triplet have extremely close frequencies in the case of long chains, since these two states are almost degenerate. In other words, this means that the singlet-triplet quasi degeneracy is not removed by small bending (or stretching) of the system. For this reason, we report in this work only the frequencies for the triplet state computed at the RHF level (Table S5). The harmonic frequencies for the cyclic chains have been calculated for the RHF singlet state (Table S6).

All the second derivatives where found to be positive (i.e., all the harmonic frequencies are real and positive). This means that, in principle, it should be possible to produce these molecular species, although their stability is also related to dynamical and thermodynamical aspects. The lowest bending and stretching frequencies vanish asymptotically with the number of atoms N , whereas the highest frequencies approach finite values. CCSD harmonic frequencies have also been calculated for both linear and cyclic chains, until N=8.

The results, presented in Table S7 and S8, do not change significantly compared to the ones obtained at the RHF level . Indeed, the difference is no more than ±35 cm -1 except for the Be 4 cyclic chain for which one frequency is found 129 cm -1 lower at the correlated level (344.73 cm -1 compared to 481.05 cm -1 ).

Extrapolated Energies

The total CCSD(T) energies can be found in Tables S3 andS4 while the per-atom CCSD(T) energies are plotted in Figure 8 as a function of 1/N (linear clusters), and 1/N 2 (cyclic clusters). For sake of comparison, the energies for Be, Be 4 , Be 13 are given in Table 2. In As already anticipated, the odd cyclic chains have a behavior that is slightly different from the even chain ones. In particular, the points lies on a straight line to a worse approximation (this is illustrated in Figure 8). Although this is hardly surprising in the case N = 3, since the strain angle is in this case very large, the behavior is much more surprising in the case N = 5. In fact, this point is less well described by a linear fitting in 1/n 2 than the even chain N = 4, which in principle has a much larger strain. We should also mention the fact that linear chains studied in previous works [START_REF] Monari | Theoretical study of Be N linear chains: optimized geometries and harmonic Frequencies[END_REF] showed effects somehow reminding the present one: in particular, we remarked a parity effect in both the ST splitting in neutral chains and in the interaction in cationic or anionic bistable chains.

This alternating behavior is found also by looking at the equilibrium bond lengths of the cyclic forms, reported in Table S2 and Figure 6 and already discussed.

The Total-Position Spread Tensor

The TPS tensor gives useful information about the metallicity of the system. The isolated atomic value is Λ = 3.9491 bohr 2 , while the tetrahedral Be 4 has the three degenerate components equal to Λ = 24.1284 bohr 2 , which gives a per-atom value of Λ = 6.0321 bohr 2 . This trend is consistent with the well known metallic nature of bulk beryllium.

The chain values are reported in Tables 3 (cyclic) and 4 (linear), and illustrated in Fig- indicates the doubly degenerate components that lie on the chain plane, while Λ ⊥ is the component in the direction orthogonal to this plane. Both singlet and triplet values were reported for the linear case. Singlet and triplet spreads, however, are almost coincident but for the shortest chains, since the corresponding wavefunctions are almost identical, except for the spin coupling of the electrons in the singly occupied orbitals.

In all cases, the spin-summed values become quickly proportional to the number of atoms, as shown in Figures 9a and10a. Moreover, from Figures 9b and10b it appears that Λ/N converges to a constant value from above, thus excluding the possibility of a slow divergence of this quantity. This is a clear indication of the non-metallic nature of the chains. As expected, the cyclic chains saturate more rapidly than the linear ones, because of the border effects that are present in the latter. In particular, the orthogonal components Λ ⊥ /N for both cyclic and linear systems converge toward a common limit that is close to 4.2 bohr 2 .

Remarkably, the parallel components are even smaller than the perpendicular ones, showing a very little mobility of the electrons in the directions of the chains. These components have obviously different behaviors for the linear and cyclic systems, since the spread is purely longitudinal in the former case, and a mixture of longitudinal and transversal spread in the latter case. However, in the case of a local character of the spread and because of its quadratic nature as a function of the coordinates, one would expect that Λ /N for the cyclic systems is close to the average between Λ /N and Λ ⊥ /N for the linear case. This is indeed the case, with Λ /N = 3.9 bohr 2 for the cycles, while Λ ⊥ /N = 4.2 bohr 2 and Λ /N = 3.5 bohr 2 for the linear geometries (average value: 3.9 bohr 2 ). This is slightly above the atomic spread (3.9491), but well below the per-atom Be 4 value for the cluster with T d symmetry (6.0321).

The per electron eigenvalues of the spread tensor for both types of chains quickly saturate for large structures (actually, they even show a tiny reduction of their magnitude). These facts clearly indicate, in our opinion, that quasi-1D beryllium systems do not have a metallic nature. On the contrary, we believe that they can be seen as prototypes of 1D covalent chains.

Spread data for FCC Be 13 , admit a less clear interpretation, because of the ionic character of the wave function. The spread for the Be + 13 cluster is obviously much less than the one for Be 2- 13 , partly because of the slightly larger bond lengths and the presence of three more electrons in the latter, but in particular for the more diffuse nature of the anionic wave function with respect to the cationic one.

Conclusions

A theoretical Coupled-Cluster study on the equilibrium geometry and harmonic frequencies of quasi-1D beryllium structures, having either linear or cyclic geometries, was reported. It was shown that the effect of dynamical correlation do not significantly change the HF results.

In conclusion, quasi-1D Be N systems are likely to be metastable clusters.

The equilibrium bond lengths of the studied systems are weakly dependent on the size of the chain or the position of the bond along the chain. Moreover, harmonic frequencies, computed at CCSD level, are in good agreement with previous CAS-SCF computations.

The reported calculations give information about the nature of the ground states of beryllium linear structures. All the indicators point toward a non-metallic, covalent type of bond between neighbor atoms. In particular, we found a HOMO-LUMO gap larger than 6 eV for the biggest clusters, there is no tendency to shrink as the size of the system increases.

The valence Molecular Orbitals of the chains can be well localized, and they correspond to bond orbitals located on pairs of neighbor atoms, contrary to what happens in the case of metallic structures. Finally, what is probably the strongest evidence of a non-metallic nature of the systems is the fact that the TPS tensor has a linear dependence on the system size.

According to our results, both linear and cyclic structures are local minima on the corresponding PES. For a given value of N, the cyclic structures are more stable than the linear ones, except in the case of Be 4 . Nevertheless, more compact (i.e., 3D) clusters are much more stable, with total energies considerably lower. The possible existence of these clusters in quasi-1D geometries depend on the barrier toward more stable arrangements. Transition states connecting these minima should be investigated to conclude.

• Table S1: CCSD equilibrium bond distances as a function of the chain length N for the linear chains.

• Table S2: CCSD equilibrium bond distances as a function of the chain length N for the cyclic chains.

• Table S3: Total energies for cyclic Be N clusters (singlet ground state).

• Table S4: Total energies for linear Be N clusters (lowest triplet state).

• Table S5: Harmonic frequencies of the normal modes for linear Be N , for N = 3, 4, ..., 12.

• Table S6: RHF Harmonic frequencies of the normal modes for cyclic Be N , for N = 3, 4, ..., 12.

• Table S7: CCSD Harmonic frequencies of the normal modes for linear Be N , for N = 3, 4, ..., 8.

• Table S8: CCSD Harmonic frequencies of the normal modes for cyclic Be N , for N = 3, 4, ..., 8.

This material is available free of charge via the Internet at http://pubs.acs.org/.

N 3 Σ u Λ 3 Σ u Λ ⊥ 1 Σ g Λ 1 Σ g Λ ⊥ 1 Σ g Λ αα+ββ 1 Σ g Λ αβ+βα 3

  the case of cyclic structures, the odd values on N show a behavior that is slightly different from the even values, as discussed earlier, the curves shown in the figure have a remarkably linear character for the largest values of N . By interpolating the energies corresponding to the last seven values (N =6,...,12) we obtain a behavior given by E lin (N )/N = -14.658 + 0the abscissa equal to zero of these straight lines corresponds to the limit for N → ∞ of the energy per monomer. The fact that the values obtained from the two different types of systems are coincident within a mhartree confirms the consistency of the model. It is also a further indication of the strong covalent nature of the bonds, since the absence of size effects (except N = 3 and N = 4 for the cyclic case) indicates the absence of any remarkable sign of electron delocalization.

ures 9 ,

 9 10, and 11. We computed both the spin-summed quantities and the spin-partitioned αα + ββ and αβ + βα ones. For a single Slater determinant, the αβ + βα term is zero for the cyclic case, and therefore only the spin-summed component is shown. For linear chains, Λ indicates the component in the chain direction, while Λ ⊥ are the doubly degenerate components in the directions orthogonal to the chain. For cyclic chains, on the other hand, Λ

Figure 1 :

 1 Figure 1: The Hückel π Molecular Orbitals of (a) benzene C 6 H 6 and (b) cyclopentadienyl anion C 5 H 5 . Orbital energies are given as a function of the Hückel parameters α and β.

Figure 2 :

 2 Figure 2: Canonical Occupied Molecular Orbitals (MO) of (a) Linear Be 6 (D ∞h ) computed at ROHF level and (b) Cyclic Be 6 (D N h ) computed at RHF level. The MO symmetries are specified and the orbital energies are given (hartree).

Figure 3 :

 3 Figure 3: The 12 valence orbitals having σ character of the Be 6 chain, local picture: the five occupied bonding orbitals (bottom); the two singly occupied edge orbitals (middle); the five empty virtual orbitals (top). The π orbitals are not shown.

Figure 4 :

 4 Figure 4: Three out of the 18 valence orbitals having σ character of cyclic Be 6 , local picture: a bonding occupied orbital (bottom); an empty non-bonding orbital (middle); an empty virtual orbital (top). The remaining 15 σ orbitals are equivalent by symmetry to these three ones by rotation around the C 6 axis. The π orbitals are not shown.

Figure 5 :

 5 Figure 5: The HOMO-LUMO energy gap obtained by ROHF/RHF calculations for both type of chains, as a function of the number of beryllium atoms N . Energies in hartree.

Figure 6 :

 6 Figure 6: The optimized CCSD bond lengths (bohr) as a function of the number of Be atoms (N), for cyclic Be N , N=3,...,12.

  0

Figure 8 :

 8 Figure 8: The CCSD(T) energies per Be atom (hartree), as a function of 1/N (linear chains) and 1/N 2 (cyclic chains).

  The spin-summed TPS for cyclic chains: planar (degenerate) components (Λ ) and out-of-plane component (Λ ⊥ ). The isolated Be atom value is 3.949061, while the Be 4 value is 24.128468. All values in bohr 2 .

Λ

  || /N Λ /N (b) The spin-summed TPS for cyclic chains, divided by the number of beryllium atoms N : planar (degenerate) components (Λ /N ) and out-of-plane component (Λ ⊥ /N ).

Figure 9 : 2 ) 1 Σ

 921 Figure 9: The spin-summed TPS for cyclic chains 34

1 Σ

 1 g Λ || /N 1 Σ g Λ /N 3 Σ u Λ || /N 3 Σ u Λ /N (b)The spin-summed TPS for linear chains, divided by the number of beryllium atoms N : perpendicular (degenerate) components and parallel component.

Figure 10 :Figure 11 :

 1011 Figure 10: The spin-summed TPS for linear chains 35

  The per-atom values are comprised between 5.5561 for Be + 13 and 7.1731 for Be 2- 13 . HCP Be 13 has values 5.9775 for the axial component and 6.0953 for the two degenerate equatorial components. It is clear from these results that the longitudinal components of the TPS in quasi-1D structures do not show any signature of a metallic nature of the wave function, contrary to what happens for the 3D case.

Table 1 :

 1 HOMO-LUMO energy gap obtained by ROHF/RHF calculations for the linear and cyclic structures. Energy values are in hartree.

	N	linear	cyclic
	3 0.260950 0.250322
	4 0.260679 0.212153
	5 0.260644 0.235867
	6 0.260686 0.217744
	7 0.260904 0.233666
	8 0.261172 0.222158
	9 0.261512 0.232881
	10 0.261825 0.225006
	11 0.262146 0.232277
	12 0.262413 0.226369

Table 2 :

 2 Total energies for the singlet ground state of the isolated Be atom and three 3-D clusters (Be 4 , Be 13 FCC and Be 13 HCP). For Be 13 FCC two undistorted clusters have been studied: the high-spin (quartet) Be + 13 and the closed-shell Be 2- 13 . For Be 13 HCP the triplet with an "ideal" geometry has been considered. Energies are in hartree.

	system	RCCSD(T)	RCCSD	RHF-SCF
	Be	-14.61896573 -14.61896573 -14.57298670
	Be 4	-58.60931559 -58.58345127 -58.35870626
	FCC Be + 13 FCC Be 2-13	-190.62705587 -190.52910207 -189.76368163 -190.93162621 -190.81969066 -189.96562397
	HCP Be 13 -190.76001277 -190.65777865 -189.85624311

Table 3 :

 3 TPS tensor for cyclic Be N clusters. The isolated Be atom value is 3.949061 bohr 2 , while the Be 4 value for the cluster with T d symmetry is 24.128468. For the FCC clusters with I h symmetry, the TPS value is 72.229108 for Be + 13 and 93.249672 for Be 2- 13 , while the HCP Be 13 has values 77.705190 (the axial component) and 79.392777 (the two degenerate equatorial components). All TPS values are in bohr 2 .

	N	Λ	Λ ⊥
	3 22.40055271 14.61175820
	4 20.10512254 17.33366066
	5 24.19144676 21.57474198
	6 26.26399922 25.42874598
	7 29.76517495 29.56328841
	8 32.95341909 33.63912437
	9 36.48514744 37.77095934
	10 39.98501422 41.89523029
	11 43.58183491 46.03325189
	12 47.19360128 50.17274976

Table 4 :

 4 TPS tensor for linear Be N clusters.

  11.36190259 13.24719585 12.25496116 13.36032816 66.51043287 -54.25547170 4 15.00274761 17.45137526 15.24779377 17.47028013 122.41199734 -107.16420357 5 18.51186839 21.62742504 18.56716211 21.63077807 193.68876105 -175.12159893 6 21.96918759 25.79452283 21.98197352 25.79502058 281.08133953 -259.09936600 7 25.40448964 29.95848755 25.40676696 29.95858478 384.79795293 -359.39118597 8 28.82989461 34.12105801 28.83050954 34.12106795 504.88414711 -476.05363756 9 32.25083491 38.28292822 32.25088095 38.28293151 641.34745876 -609.09657781 10 35.66977175 42.44445216 35.66981185 42.44445205 794.18957727 -758.51976542 11 39.08783634 46.60574585 39.08783014 46.60574593 963.40895434 -924.32112419 12 42.50563131 50.76689007 42.50563579 50.76688997 1149.00329929 -1106.49766350

Acknowledgement

This work was partly supported by the French "Centre National de la Recherche Scientifique" (CNRS, also under the PICS action 4263). It has received fundings from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement n o 642294. This work was also supported by the Programme Investissements d'Avenir under the program ANR-11-IDEX-0002-02, reference ANR-10-LABX-0037-NEXT. The calculations of this work have been partly performed by using the resources of the HPC center CALMIP under the grant 2016-p1048. One of us (AWH) acknowledges the support of the "Theoretical Chemistry and Computational Modeling" (TCCM) Erasmus-Plus Master program. Finally, we would like to thank Ulf Saalman (Max-Planck-Institut für Physik komplexer Systeme, Dresden), for interesting discussion and suggestions.

Supporting Information Available

The document SI-Be.pdf contains:


Houda Nouasse
email: houda.nouasse@mines-douai.fr

Arnaud Doniec

Guillaume Lozenguez

Eric Duviella

Pascale Chiron

Bernard Archimède

Constraint satisfaction problem based on

Keywords: Resilience, Constraint satisfaction problems, Graphs, Navigation networks, Water management, Global change

INTRODUCTION

Inland navigation transport is part of the multimodal transport that is promoted by the Trans-European network program (TEN-T 1 ). This program aims at developing new transport infrastructure policy to close the gaps between Member States' transport networks and to guaranty seamless transport chains for passenger and freight. To well integrating the inland navigation transport in this framework, an efficient water management strategy is required. It consists in guaranteeing the navigation conditions even if an increase of navigation demand and extreme events due to climate change is expected [START_REF] Bates | Climate change and water[END_REF]. The intergovernmental panels as the IPCC (Intergovernmental Panel on Climate Change) has defined RCP scenarios (Representative Concentration Pathways) on which future forecasts on temperature and rain can be generated IPCC (2014). Based on these scenarios, the flood and drought events will be more frequent and mode intensive in close future Boe ét al. (2009); [START_REF] Dechter | Climate change impacts on water resources and hydrological extremes in northern france[END_REF]; [START_REF] Wanders | Human and climate impacts on the 21 st century hydrological drought[END_REF]; [START_REF] Li | Assessing vulnerability and adaptive capacity to potential drought for winter-wheat under the rcp 8.5 scenario in the huang-huai-hai plain[END_REF]; [START_REF] Park | Drought prediction till 2100 under {RCP} 8.5 climate change scenarios for korea[END_REF]. Thus, constraints on water resource management for navigation will be bigger.

To deal with the navigation demand increase and the climate change effects, adaptive water management strategies have to be designed. An adaptive and predictive control architecture was proposed in [START_REF] Duviella | Adaptive and predictive control architecture of inland navigation networks in a global change context: application to the cuinchy-fontinettes reach[END_REF]. It is based on the multiscale modeling approach proposed in [START_REF] Duviella | Multi-scale modeling approaches of inland navigation networks for their management in a global change context[END_REF] to

⋆
This work is a contribution to GEPET-Eau project which is granted by the French ministry MEDDE-GICC, ORNERC and the DGTIM. http:// gepeteau.wordpress.com/enversion/. 1 http://ec.europa.eu/transport/themes/ infrastructure/ten-t-guidelines/index_en.htm reproduce the dynamics of inland navigation networks during flood and drought events. This architecture is improved in this paper to consider events that can impact large areas on larger periods. It consists in dispatching volumes of water through the network to guarantee enough water in each part of the network. Thus, the designed tools aim at determining the resilience of the inland navigation networks. They are based on an integrated model of inland navigation network that allows locating the navigation reaches, the locks and gates and identifying the main water intakes. Then a flow graph is proposed to determine the possible paths between the navigation reaches and the main constraints on water volume exchanges. Network flow problems have been widely used for the computation of maximum flow or minimum cost flow in several areas such as transportation Silver and de Weck (2007), telecommunication [START_REF] Fekete | The maximum energy-constrained dynamic flow problem[END_REF], job scheduling or flood attenuation [START_REF] Nouasse | A flood lamination strategy based on transportation network with time delay[END_REF]. Finally, a constraint satisfaction problem is defined to determine the water volumes that have to be exchanged between each navigation reach. Mathematical programming in general [START_REF] Passchyn | Mathematical programming models for lock scheduling with an emission objective[END_REF] and constraint satisfaction in particular [START_REF] Sun | Combining csp and mpc for the operational control of water networks: Application to the richmond case study[END_REF] is a major tool to address the problems of transport.

The management objectives of waterways are given in Section 2. Section 3 allows formulating the problem of this management during drought and flood periods. The integrated model is detailed. Then, the constraint satisfaction problem is proposed. All the designed tools allows considering inland navigation networks that are composed of several confluents and diffluents. In Section 4, an academical example of inland navigation networks composed of five interconnected navigation reaches is considered. Its characteristics have been determined We can assume that the resilience study is an approach based on the fact to propose a set or rather a system of rules for maintaining "the proper" functioning or sometimes simply the functioning of an inland navigation network. The main objective would be to adapt to extreme conditions such as floods or drought. The needs of each of these situations are unlike those of the other that is why we need to establish a stable and adaptive resilience "system" (as we suppose that is a system of rules). Stable in the sense to resist change and adaptive in the sense of accompanying that change. When one fails to offer a solution, the system does not fail but reveals its limits (it can also be considered a form of adaptation). An understanding of the functioning and needs of each component of the network, in particular NR, is an important step in determining strategies that will be adapted on a microscopic level to each of them and on a macroscopic point of view to their interaction with the entire network and climate change. To address this issue, authors in [START_REF] Duviella | Adaptive and predictive control architecture of inland navigation networks in a global change context: application to the cuinchy-fontinettes reach[END_REF] proposed an adaptive and predictive control architecture.

PROBLEM MODELING

Inland navigation reach model

The integrated model is proposed to model several configurations of inland navigation networks by considering two elementary configurations: a confluence and a diffluence (see Figure 1.a). Networks are composed with a finite number η of interconnected NR. NR are numbered and denoted N R i , with i ∈ 1 to η. The N R i is modeled as a tank that contains a volume of water, denoted V i (t). According to the boundaries of the navigation rectangle, to the NNL and to the geometrical characteristics of each N R i , it is possible to determine the volumes that corresponds to the NNL, the HNL and the LNL such as

V LN L i ≤ V N N L i ≤ V HN L i . The management objective is V i (t) = V N N L i and at least V LN L i ≤ V i (t) ≤ V HN L i
. If this condition is broken, the navigation has to stop.

A N R i is supplied and is emptied by controlled and uncontrolled water volumes (see Figure 1.b). Controlled water volumes gather the water that is coming from controlled gates and from the lock operations. Uncontrolled water volumes are all the withdrawals and supplies from water intakes located along the N R i . It is also possible to consider the water exchanges with groundwater. Thus, the set of controlled water volumes is composed of:

• controlled volumes from the upstream NR that supply the N R i , denoted V s,c i (s: supply, c: controlled), • controlled volumes from the N R i that empty the N R i , denoted V e,c
i (e: empty), • controlled volumes from water intakes that can supply or empty the N R i , denoted V c i . These volumes are signed; positive if the N R i is supplied, negative otherwise.

The set of uncontrolled water volumes is composed of:

• uncontrolled volumes from natural rivers, rainfall-runoff, Human uses, denoted V u i (u: uncontrolled). These volumes are signed depending of their contribution to the volume V i (t) in the N R i .

• uncontrolled volumes from exchanges with groundwater, denoted V g,u i (g: groundwater). These volumes are also signed.

Based on the definition of the water volumes that contribute to the volume contained in the N R i , it is possible to model its dynamics by:

V i (t) = V i (t -1) + V s,c i (t) -V e,c i (t) + V c i (t) + V u i (t) + V g,u i (t). (1)
The dynamics of the N R i have to take into account the configuration of the network. For a confluence, the controlled volumes coming from all the NR that are located upstream the N R i are added. For a diffluence, the controlled volumes that empty the N R i correspond to the sum of the controlled volumes that supply the downstream NR (see relation ( 2)).

       V s,c i (t) = j∈Ω i V e,c j (t), V e,c i (t) = j∈Θ i V s,c j (t), (2)
where Ω i gathers all the index of the NR that supply the N R i , and Θ i all the index of the NR that are supplied by the N R i .

The proposed integrated model requires the knowledge of all the possible water volume contributions. It is necessary to identify all the controlled volumes and to estimate all the uncontrolled volumes. This knowledge comes from the waterways managers.

In this paper, we propose to set some assumptions. Water exchanges with groundwater are nonexistent. The controlled water volumes V c i from water intakes are bounded with the known minimal and maximal daily volumes of water. The water volume that corresponds to each lock operation is known. Finally, the daily average number of ships that cross the inland navigation network is known.

Flow-based network modeling

We choose in the following to model inland navigation network as a graph. The answers related to questions about the specific requirements for operating a navigation system are made by considering a flow problem on this graph. The volume management in such a network is guaranteed by satisfying the constraints formulated for this problem.

Let consider a network, G = (G x ; G a ), defined as a connected Directed Acyclic Graph (DAG), where G x is the set of nodes and G a is the set of arcs. The set of nodes G x contains the nodes representing each reach in the navigation network N R i and two additional nodes, O and N , respectively source and sink. We denote each directed arc a = (N R i , N R j ) with nodes corresponding to navigation reaches, such as N R i is the leaving node (the tail of arc a) and N R j is the entering node (the head of arc a). To simplify notations, we will consider that nodes are labeled with the index of each N R, i.e. a = (N R i , N R j ) = (i, j) or a = (O, N R j ) = (O, j). As mentioned in (3.1), the link between two adjacent reaches N R represents a water volume amount. Thus, we consider an s-t flow routing (sourceto-sink flow) in the network described by the water volume path between different reaches. Figure 2 shows an example of a flow network model by considering the navigation network proposed in Figure 1. The arcs are directed according to the configuration of the navigation network. According to the proposed assumptions, the node O gathers all the volumes of water that supply the navigation network from natural rivers. The node N retrieves all the volumes of water from the navigation network. On every arc in G a , it is defined a flow variable φ a , a ∈ G a that can be expressed by φ i j . This flow varies with respect to capacities constraints and demands relation on each node.

One step in the management of water volumes in the inland navigation network is, for example, if one refers to this network G, to determine the flow that guaranties different navigation scenarios for a given period of time. That can be formulated as follows:

Maximize

φ (3) subject to l a ≤ φ a ≤ u a for a ∈ G a Φ ≥ 0 (4)
by considering these additional constraints

φ a + -φ a = d(i) for i ∈ G x -{O, N } , (5)
This constraints allows flow on a given node i, not to be conservative. In order to respect the navigation conditions on each reach N R i , the demand of water volumes on each node representing a N R is defined as:

d(i) = V N N L i -V i (t), i ∈ G x -{O, N } (6)
that must verify the following relation:

V N N L i -V HN L i ≤ d(i) ≤ V N N L i -V LN L i i ∈ G x -{O, N } (7)
where Φ is a vector that contains all the flow φ a , l a and u a the lower and upper bound capacities of the arc a, a + the arcs leaving the node i and a -, the arcs entering the node i, and d(i) the demand of the node i. We also assume that

φ a + (O) = d(O) and φ a (N ) = -d(O) as d(O)
is the known supply on the source node.

The design of the integrated model (see Figure 1.b) helps to build the flow network and to determine the capacities that are expressed in volume. The lower and upper bound capacities of the arc a, i.e. l a and u a , depend on the configuration and equipment of the inland navigation networks. Hence, we define Ψ the set of the index of the upstream N R that are not supplied by another N R, and Ξ the set of the index of the downstream N R that not supply another N R. By considering example in Figures 1 and2, these sets are Ψ = {i -2, i -1}, and Ξ = {i + 1, i + 2}. Thus:

• upper bound capacities for arcs between two N R, i.e. {(i -2, i), (i -1, i), (i, i + 1), (i, i + 2)} in Figure 2, are computed as the sum of the maximum available controlled water volumes (V s,c i ), • lower bound capacities for arcs between two N R, are only the required water volumes for the navigation (V s,c i ), depending on the b ∈ N the number of ships that cross each N R, • upper bound capacities for arcs between O and N R j , j ∈ Ψ, i.e. {(O, i -2), (O, i -1)}, are computed as the sum of the maximum available controlled water volumes (V s,c i ), the maximum available water volumes from water intakes (maximum positive V c i ) and the positive uncontrolled water volumes (V u i ), • lower bound capacities for arcs between O and N R j , j ∈ Ψ, are the sum of the required water volumes for the navigation (V s,c i ) and the positive uncontrolled water volumes (V u i ), • upper bound capacities for arcs between O and N R j , j ∈ Ψ, i.e. {(O, i), (O, i + 1), (O, i + 2)}, correspond to the sum of the maximum available water volumes from water intakes (maximum positive V c i ) and the positive uncontrolled water volumes (V u i ), • lower bound capacities for arcs between O and N R j , j ∈ Ψ, are only the positive uncontrolled water volumes (V u i ), • upper bound capacities for arcs between N R j , j ∈ Ξ, and N , i.e. {(i -2, N ), (i -1, N ), (i, N )}, correspond to the sum of the maximum water volumes that can empty the N R (minimum negative V c i ) and the negative uncontrolled water volumes (V u i ), • lower bound capacities for arcs between N R j , j ∈ Ξ, and N , are only the negative uncontrolled water volumes (V u i ), • upper bound capacities for arcs between N and the N R l , l ∈ Ξ, i.e. {(i + 1, N ), (i + 2, N )}, are the sum of the maximum available controlled water volumes (V e,c i ), the maximum water volumes that can empty the N R (minimum negative V c i ) and the negative uncontrolled water volumes (V u i ), • lower bound capacities for arcs between N and the N R l , l ∈ Ξ, are the sum of the required water volumes for the navigation (V e,c i ) and the negative uncontrolled water volumes (V u i ).

CSP & Relaxation

In the following, we propose a CSP approach to determine the resilience of our system in the flow network G. Indeed, by definition CSP allows to evaluate in a finite domain a number of possible solutions. Unlike classical methods of solving a flow problem where the solution consists of finding a s-t flow augmenting path, CSP solving is based on search algorithms that aim on domain exploration for values of consistent variables. For example one may want to "augment" the water volume between two reaches but in that case we can not maximize flow between these two reaches so the arc connecting them is excluded from the optimal path. While in a CSP perspective, we will try to "best" satisfy this constraint of water volume required between the two reaches and if it is not satisfied it will be considered as an event to address in our resilience study. Hence, the approach we consider, for volume management in the navigation network, is to first solve a constraint satisfaction problem that can check as a priority that the level of water in the reaches is as close as possible to the NNL, and releases and river inputs are near zero. Such a problem is over-constrained and given the conflictual aspect of the problem: in order to maintain the water level around the NNL between two reaches while one may be a reach water lack while the other is in surplus of water; we then need to consider more flexible constraints such tolerating releases and river inputs. Which generally comes back to consider a Max-CSP problem [START_REF] Dechter | Climate change impacts on water resources and hydrological extremes in northern france[END_REF] when it comes to satisfy a maximum number of constraints. Starting from an over-constrained problem we relax to ensure the satisfaction of new and more permissive constraints. In this formulation as variables, we have:

• φ a , flow variable on each arc a ∈ G a , the domain D(φ a )
of φ a expressed as a discrete interval whose bounds are calculated according to the number of ships b. • d(i) demand variable for every reach i in the network. The domain D(d(i)) is set according to navigation conditions (see 3.2).

In addition to flow conservation constraint on each node, we consider that the two following constraints are possible:

|φ a -min(D(φ a ))| = 0, (8) when a ∈ {(i -2, N ), (i -1, N ), (i, N )} and a ∈ {(O, i), (O, i + 1), (O, i + 2)}.
With min(D) and max(D) respectively the minimal and maximal values of variables domain. This constraint helps reduce the relative contributions to releases or external resources compared to the ensembles of contributions.

|d(i) -V N N L i | = 0, i ∈ G x -{O, N } (9)
This constraint helps to keep water volume on each reach around NNL. We consider a relaxation of the problem in order to enlarge the set of solutions. For this, during the resolution, the constraints above ( 8) and ( 9) will be relaxed one after the other which comes in some way to progressively identify the responsible arcs of problem infeasibility.

| a φ a - a max(D(φ a ))| = 0, (10)
a ∈ {(i -2, N ), (i -1, N ), (i, N )} and a ∈ {(O, i), (O, i + 1), (O, i + 2)} | i d(i) -z| = 0, i ∈ G x -{O, N } (11) such as, i min(D(d(i))) ≤ z ≤ i max(D(d(i))), (12)
The following algorithm show how constraint programming is integrated to calculate the maximum flow.

SIMULATION RESULTS

We apply our method on a network model composed of five reaches respecting the configuration described in Figure 3.a. This configuration is characterized by the presence of a junction reach. Each reach is supplied with a daily volume V s,c i . For example, N R 1 is supplied by operations of the lock L 1 with a

N R N R 1 N R 2 N R 3 N R 4 N R 5 V N N L i 11.
V s,c,L i 5b 1 3b 2 3.7b 1 + 3.5b 2 6b 4 3b 5 V s,c,G i 0 [0, 0.5T ] [0, 3T ] [0, 2T ] [0, 1T ] V e,c,G i [0, 2T ] [0, 1T ] [0, 2.5T ] [-0.5T, 0] [0, 3T ] V e,c,L i 3.7b 1 3.5b 2 6b 4 + 3b 5 25b 4 2b 5 V c,O i [0, 1T ] [0, 0.5T ] [0, 1T ] [0, 3T ] 0 V c,N i [-1T, 0] 0 [-2T, 0] [-5T, 0] [-2T, 0] V u i 1T 0.5T 0.75T 0.5T -0.25T
T is period of management strategy, T in hours. These volumes were estimated based on average discharges. Negative volume values are used to distinguish an emptying situation. These data represent a normal navigation period given a management horizon T expressed in hours. We considered two management scenarios during drought and flood. The parameters of the problem described in 3, were calculated in Table 1 on a 24hour management horizon. The algorithm 1 was implemented in Java using the Choco Open Source Java library [START_REF] Jussien | Choco: an open source java constraint programming library[END_REF] for computing the solution of the CSP proposed. The results presented in Figures 4 and5, refer respectively to scenarios of drought and flood. We assumed that the number of ships navigating in the network is known a priori, denoted b, with b = b j , j = 1, . . . , 5. In each of these situations, we represent the network state when the maximum number of ships is reached and when the volume of water in reaches is solicited.

In Figure 4, we studied three scenarios with drought case considering that volumes of input water resources compared to the normal period (see Table 1) were decreased by 25% (see Figure 4.a), 50% (see Figure 4.b) and finally 75% (see Figure 4.c). The impact of drought manifests itself in the network, so that, from a number of ships; downstream reaches (N R 3 and N R 4 ) are undersupply because the demand for water consumption becomes important (values in red (see Figure 4)) in this case study one of the upstream reaches (N R 2 ) is in surplus. Resources and water requirements in the network are not being distributed in a uniform basis, so that each part of the network respond to this lack of water differently. In addition to this navigation limitations are to expect in this network in a drought case, i.e., navigation can not be tolerated for a number of ships greater than 36/day when a decrease of 25% of water resources is indicated and not to a number greater than 28/day when there is a decrease of 75%. Downstream reaches emptied as soon as we have 9 ships that circulate in the network when the percentage of drought is 75% and 21 ships when there is a percentage of 25%. In the Figure 5, two flood scenarios are In this case, we notice that, at the contrary of the drought case, upstream reaches (N R 1 and N R 2 ) are both in surplus. N R 3 as a junction reach, it is rather on surplus or undersupply which is the case of the downstream reach N R 4 . Water volume changes in reaches are apparent even for one ship navigating in the network in the case of flood percentage of 200% (See Figure 5(b) left). The difference is substantial compared to a 100% increase when the reaches were filled of water or emptied as from a number of ships equal to 28 in the network (See Figure 5(a) left). On the other hand, the maximum number of ships allowed to navigate in the case 100% percentage of flood is 43 whereas it is 45 in the case of 200%. This means that the way in which the volume of water is distributed in the network, depending on the interactions between the reaches, is important in determining the maximum capacity of the network in terms of navigation. In this paper, a constraint satisfaction problem based on flow transport graph is proposed to determine the resilience of inland navigation networks against increase of the navigation demand and extreme events as drought and flood periods. The proposed tools are based on an integrated model of inland navigation networks. The constraint satisfaction problem allows determining the optimal management of water volumes during extreme events. The main objective is to keep the navigation conditions on each section of the navigation networks. If the navigation conditions are guaranteed, navigation transport could take part of the multimodal transport as a competitive transport mode. The proposed tools are tested using a network composed of five interconnected navigation reaches. The results show the efficiency of the method to evaluate the ability of the network to satisfy navigation conditions in the case of extreme climatic events. In future works, uncertainty on navigation demand and on available water resource will be considered.

Fig. 1 .

 1 Fig. 1. (a) Inland navigation network, (b) its integrated model.

Fig. 2 .

 2 Fig. 2. Flow Network Model.

Fig. 3 .

 3 Fig. 3. (a) Inland navigation network composed of 5 NRs, (b) its integrated model. daily volume V s,c,L 1, and by the gate G 1 with a daily volume V s,c,G 1 . The characteristics of the NR are summarized in Table1, with volumes from gate V s,c,G

Fig. 4 .

 4 Fig. 4. (a) Drought 25% left b = 21 and right b = 36, (b) Drought 50% left b = 18 and right b = 32, (c) Drought 75% left b = 9 and right b = 28 studied in the network. Figure 5(a) describes the distribution of water volumes in the network when a flood estimated at an increase of 100% of water resources is indicated. The second scenario (see Figure 5(b)) represents an increase of about 200% of water resources.In this case, we notice that, at the contrary of the drought case, upstream reaches (N R 1 and N R 2 ) are both in surplus. N R 3 as a junction reach, it is rather on surplus or undersupply which is the case of the downstream reach N R 4 . Water volume changes in reaches are apparent even for one ship navigating in the network in the case of flood percentage of 200% (See Figure5(b) left). The difference is substantial compared to a 100% increase when the reaches were filled of water or emptied as from a number of ships equal to 28 in the network (See Figure5(a) left). On the other hand, the maximum number of ships allowed to navigate in the case 100% percentage of flood is 43 whereas it is 45 in the case of 200%. This means that the way in which the volume of water is distributed in the network, depending on the interactions

Fig. 5 .

 5 Fig. 5. (a) Flood 100% left b = 28 and right b = 43, (b) Flood 200% left b = 1 and right b = 45

  H. Nouasse * A. Doniec * G. Lozenguez * E. Duviella * P. Chiron * * B. Archimède * * K. Chuquet * * * Constraint satisfaction problem based on flow graph to study the resilience of inland navigation networks in a climate change context ⋆considering real navigation networks. This case study allows detailling the design step of the proposed tools. Finally, drought and flood events are simulated to highlight the performances of the designed tools.

	2. WATERWAYS MANAGEMENT
	Waterways are equipped and opened to passenger transport,
	cargo and boating. Their different uses have an economic bene-
	fit Mihic et al. (2011); Mallidis et al. (2012); Brand et al. (2012).
	It is observed that inland navigation network is constituted by
	a number of hydraulic structures, including locks. A part of
	an inland navigation network between two locks is defined
	as a Navigation Reach (NR). It is assumed that in general
	navigation requirements are the same for each NR. Navigation
	is maintained in a reach respecting the so called navigation
	rectangle. The boundaries of the navigation rectangle are the
	High Navigation Level (HNL) and the Low Navigation Level
	(LNL). The main management objective consists in keeping the
	water level in each NR is inside the navigation rectangle and
	close to the Normal Navigation Level (NNL). The waterways
	have to be supplied with natural rivers. Climate change impacts
	severely on the availability of water resources, more accurately
	during floods and drought. Flood periods are as problematic
	as droughts. For example, in extreme climate scenarios natural
	water reserves may reach their ecological limits in the absence
	of management of excess volumes of water. This case repre-
	sents a deficit situation in water resources during drought. Thus,
	the resilience study of inland navigation networks is necessary.

Table 1 ,

 1 with volumes from gate V s,c,G

	i	and from lock V s,c,L i	. Where

Table 1 .

 1 Characteristics of the NR, with volumes in thousand of m 3 , volumes as V s,c,G i computed according to available discharges that are expressed in m 3 /h, b j the operation number of lock L j