Mathias Rapacioli
email: mathias.rapacioli@irsamc.ups-tlse.fr

Fernand Spiegelman

Anthony Scemama

André Mirtschink

Modeling charge resonance in ionized molecular clusters: combining DFT-Tight Binding with configuration interaction

published or not. The documents may come

Introduction

The description of molecular clusters requires to take into account the various contributions of the intermolecular energy, including the Pauli repulsion, the Coulomb interaction, and the London dispersion. The description of the electronic structure of singly ionized molecular clusters requires to consider two further essential contributions, the first one being charge resonance which may cause the charge to be partially or totally delocalized over the molecular units, the second one being the polarization contribution due to the influence of the charge. Both of them yield a stabilization of the charged species as compared to the analogous neutrals. A proper description requires correct balance between charge delocalization and polarization forces.

While Density Functional Theory (DFT) is an appealing method to describe the electronic properties of clusters up to a few tenths, maybe a few hundreds atoms, at least in single point calculations, the use of the most common functionals is known to fail to describe properly dispersion forces. This is a first handicap to deal with molecular clusters. Search for new functionals accounting for dispersion [2][3][4][5][6][7][8][9][10] is a very active field, while semi-empirical corrections to standard DFT calculations are also used [11][12][13][14][15][16][17][18] . Description of charge resonance in molecular clusters is another serious problem in standard Density Functional approaches. Using canonical Kohn-Sham orbitals, one arising problem is the self-interaction of the delocalized charge. This aspect is particularly prevalent at dissociation, for instance in a cationic molecular complex involving two identical units, and which should dissociate into one molecular cation and one neutral, whereas in restricted scheme DFT the charge is asymptotically equally shared by the two units, breaking the energy additivity and further introducing a spurious Coulomb interaction between the two moieties. Although such an artifact is essential at the dissociation, it is also expected to play a role all over the potential energy surface, including the equilibrium geometries.

A correct description of dissociation is in principle easily obtained if one uses a multiconfigurational nature of the wave function. This can be achieved by high level methods like Configuration Interaction (CI) based methods (Multi-Configurational Self-Consistent Field MCSCF, Multi-Rereference Configuration Interaction MRCI) or Coupled Cluster approaches but at a high computational price. Such calculations may provide benchmarks on reasonably small systems (essentially dimers) but rapidly exceed today's possibilities as soon as the molecular units exceed a few tens of atoms.

One of the tracks to circumvent the drawbacks of present state DFT in an ab initio framework is to combine CI for describing long range (lr) electron-electron interactions and DFT for the short range electron-electron interactions (sr). This gave rise to the lr-sr formalism following Savin's formulation, [START_REF] Savin | Recent developments and applications of modern Density Functional Theory[END_REF][20][21][22] which make possible combinations of MP2, Coupled Cluster and/or CI approaches with Density Functional Theory, via a partition of the electron-electron Coulomb interaction. This formulation is quite attractive, nevertheless its numerical cost is significantly larger than that of a standard DFT calculation.

On the other hand, charge resonance (or excitation resonance) appears quite simply in valence bond-like approaches [23][START_REF] Pauling | The nature of the Chemical Bond[END_REF][START_REF] Murrel | The Chemical Bond[END_REF][START_REF] Shaik | A Chemist's Guide to Valence Bond Theory[END_REF] explicitely taking into account the multiconfigurational nature of the wavefunction via the definition of a basis in which charge (or excitation) is localized on a given fragment of the system. This is the essence of the excitonic models originating from solid states, but also used in molecular materials and even biological systems. An application to molecular cluster cations and in particular polycyclic aromatic hydrocarbon (PAH) cluster cations was published by Bouvier et al. [START_REF] Bouvier | [END_REF] defining a resonance charge model based on frozen molecules, and parametrized from ab initio CI calculations of dimers. Diatomics-In-Molecule modeling of singly ionized rare gas clusters can also be expressed in valence bond picture in a basis of atom-localized hole configurations with no internal geometrical structure, ? ? paving the way for extensive simulations of the electronic and dynamical properties of ionized rare gas clusters (see for instance ? ?

.

More recently, the concept of using a valence bond configuration description combined with a DFT framework was proposed by the group of Van Voorhis et al. 1,[28][29][30] to investigate charge delocalization in mixed valence compounds exhibiting possible bi-stability with the perspective of controlling charge transfer. They developed a method combining constrained DFT with a small CIlike scheme (CDFT-CI) to deal with charge delocalization in extended systems. This is extremely appealing from the computational point of view, since the size of the CI increases linearly with the number of fragments, and not as a power of the total number of active electrons in the whole system, which is a bottleneck. The Self-Consistent-Charge Density Functional based Tight Binding (SCC-DFTB) [31][32][33][34] is an alternative to DFT in the quest of addressing large systems. It is derived from DFT through several approximations allowing the use of parametrized tables to avoid the explicit calculation of overlap and interaction integrals. As SCC-DFTB is derived from DFT it also inherits its lack in describing charge resonance. Recently, we presented a preliminary transcription of the CDFT-CI method in SCC-DFTB framework using also approximations to determine the CI couplings. We provided applications to coronene clusters with constrained geometries, because of lack of the gradient 35 .

In the present paper, we present the general adaptation of the CDFT+CI method to the SCC-DFTB framework, with the aim of future investigations of charge resonance in molecular clusters with either large units or large size. This method is called the DFTB-VBCI (Valence Bond CI). In order to perform geometry optimization, we also derive analytical expressions for the energy gradient with respect to the nuclear coordinates, which we use to achieve full structural optimization for the benzene dimer and water dimer cations, respectively. Section 2 is devoted to the presentation of the general methodology and the DFTB-VBCI approach and the derivation of analytical expressions for the nuclear forces. In section Section 3, we benchmark the method on benzene and water ionic dimers on the basis of comparisons with high level calculations. Outlines and perspectives are given in section Section 4.

Methodology

The DFTB-VBCI method is an adaptation of the DFT+CI approach 1,28-30 to the SCC-DFTB scheme with the focus on treating charge resonance in ionized molecular clusters. In that approach, the wave function of the system Ψ is decomposed on a basis {Φ I } of configurations where the charge is localized on different fragments of the system. The intuitive decomposition of a molecular cluster leads to identify each of the N frag monomers to a fragment and the wave function becomes

Ψ = N frag ∑ I b I Φ I (1)
where Φ I is the configuration where the charge is fully carried by fragment I. Each charge localized configuration Φ I is a single Slater determinant, built from the molecular orbitals (MO) {φ I i } resulting from a constrained SCC-DFTB calculation. These VB like configurations then interact within a small CI-like scheme giving their coefficients b I in the wave function and the ground state energy.

In this methodological part, we first briefly recall the SCC-DFTB scheme basics (Section 2.1) before explaining the derivation of the charge localized configurations Φ I using the constrained SCC-DFTB (Section 2.2) and the CI-like scheme calculation (Section 2.3). We present then analytical expressions for the nuclear forces (Section 2.4) and some further approximations to accelerate the approach (Section 2.5). We adopt different font conventions to distinguish between matrices expressed on different basis sets. For instance, the Hamiltonian matrix is written as H in the atomic orbital (AO) basis set; H in the molecular orbital (MO) basis set and H in the determinant basis set (the basis of the charge-localized configurations.

DFT and SCC-DFTB

Several reviews on the DFTB and SCC-DFTB method can be found in the litterature. [31][32][33][34] The SCC-DFTB differs from the DFT expressed on a local basis set by the following approximations:

(i) the DFT energy is expanded up to the second order with respect to charge density fluctuations around a given reference density (ii) all three center interaction integrals are neglected (iii) the MO are expressed in a minimal atomic basis set

φ i = ∑ µ c iµ ϕ µ (2)
(iv) the short distance repulsive potential is expressed as a function of two body interactions (v)

the second order term in the DFT energy expansion is expressed as a function of atomic Mulliken charges and a Γ matrix. With these approximations, the total SCC-DFTB energy reads

E SCC-DFTB = atoms ∑ α,β =α E rep αβ + ∑ i n i φ i | Ĥ0 |φ i + 1 2 atoms ∑ α,β Γ αβ q α q β ( 3
)
where Ĥ0 is the Kohn-Sham operator at the reference density and

E
∑ ν c iν (H µν -ε i S µν ) = 0 ∀µ, i (4)
S is the atomic basis overlap matrix and the Hamiltonian matrix reads H = H 0 + H 1 with

H 1 µ∈α;ν∈β = 1 2 S µν atoms ∑ ξ (Γ αξ + Γ ξ β )q ξ (5)
where µ ∈ α means that the atomic orbital µ belongs to atom α.

Constrained SCC-DFTB

Similarly to the constrained DFT 28 , the MOs {φ I i }, used to build the configuration Φ I , are obtained from a minimization of the SCC-DFTB energy with the constraints that the charge is carried by fragment I and that the orbitals are orthonormalized. The corresponding Lagrangian is

L = E SCC-DFTB ({φ I i }) + ∑ i j Λ I i j φ I i |φ I j -δ i j +V I ∑ i n i φ I i | PI |φ I i -N I (6)
where V I is the Lagrange multiplier ensuring the charge localization constraint, PI is the projector of the density on fragment I, N I is the number of electrons on fragment I which defines the charge localization, Λ i j are the Lagrange multipliers ensuring the orbitals orthonormality constraint. Wu and Van VoorhisWu and Van Voorhis 28 discussed the effect of several localizations, based on different charge definitions (Mulliken, Boys, Löwdin), on the constrained energy and finally adopting the Löwdin approach. We used for the constrained SCC-DFTB the Mulliken charge definition because (i) the defects of Mulliken charges are less crucial than in DFT due to the use of a minimal atomic basis set (no diffuse functions) and (ii) in the mostly used version, SCC-DFTB is a Mulliken charge based approach and all the matrix elements have been parametrized for this charge defintion. This choice leads to the expression for the constraint

∑ iν µ n i c I iν c I iµ P I ν µ = N I ( 7
)
P I being the projection matrix expressed as 28

P I µν =            0 if µ / ∈ I and ν / ∈ I S µν if both µ ∈ I and ν ∈ I 1 2 S µν for other cases (µ ∈ I or ν ∈ I)
The secular equation (Eq. ( 4)) should now be solved with the H matrix modified as follows

H µν = H 0 µν + H 1 µν +V I P I µν
Similarly to the constrained DFT, Eq. ( 4) must now be solved self-consistently over the atomic charges and contains an unknown Lagrange multiplier V I . To overcome some convergence problems, we have implemented three ways of solving this equation which can be used alternatively until one of them converges.

-i-The first one (similar to that of 36 ) consists in solving the secular equation with an inner loop to calulate the V I Lagrange multiplier satisfying the constraint keeping fixed charges q α in the

Hamiltonian. An external loop over the Mulliken charges is performed to reach a self-consistent solution.

-ii-The second approach consists in inverting the two previous loops, i.e. the inner loop ensuring the self-consistence over the Mulliken charges and the external loop allowing the determination of the Lagrange multiplier V I .

-iii-The third approach is somewhat different and consists in three steps. First, a MO guess is generated for isolated fragments. The full set of MOs is is orthonormalized with a Löwdin procedure. These MOs do not correspond to an energy minimum and do not satisfy the charge localization constraint. In the second step, the MOs evolve to change the number of electrons on fragment I with the iterative procedure

φ I i (n + 1) = φ I i (n) + α p I φ I i (n) + ∑ j φ I j (n)Λ i j ∀i ( 8
)
where n is the iteration step. The third term ensures the orthonormalization constraint. Transposing this equation in the atomic basis set gives the evolution of the MOs

C I (n + 1) = C I (n) + α(S -1 P I C I (n) + XC I (n)) (9)
where X = αS -1 Λ. At each step, the α coefficient is adapted to increase or decrease the number of electrons on fragment I and the X matrix is calculated solving a second order equation equivalent to the Rickaert 37 algorithm already implemented for SCC-DFTB Car Parrinello molecular dynamics 38 . Once a solution satisfying the density constraint is achieved, the last step consists in relaxing the MO to minimize the energy, conserving the charge localization and orthonormality constraints

φ I i (n + 1) = φ I i (n) + α dE dφ I * i + ∑ j φ I j Λ i j +V I P I φ I i ( 10
)
giving the evolution of the coefficients

C I (n + 1) = C I (n) + α(S -1 HC I (n) + XC I (n) +V I S -1 P I C I (n)) (11)
This step requires both the calculation of X and V I . Starting from a given V I (the one of the previous step if n > 1), C I (n + 1) is determined, calculating X with the Rickaert algorithm. The charge carried by fragment I with these new coefficients is calculated. If this charge is too large (resp. too small), V I is decreased (resp. increased). The process is repeated until the molecular charge satisfies the constraint. Finally, the MOs converge to the the charge-localized solution.

The Configuration Interaction-like scheme

The set of MO {φ I i } obtained from a constrained SCC-DFTB is used to build the charge-localized configurations Φ I as single Slater determinants. The coefficients b I of these configurations in the total wave function Ψ (see Eq. ( 1)) are obtained by solving the CI-like scheme

       
H n1 . .. H nn                   b 1 b 2 : b n          = E          S 1 S 12 .. S 1n
S 21 S 2 .. S 2n : : :: :

S n1 . .. S nn                   b 1 b 2 : b n          (12)
where S IJ is the two configurations overlap Φ I |Φ J and H II is the energy of the configuration Φ I already calculated with the constrained SCC-DFTB. Following the approach of Wu et al. 1,29 , the coupling elements H IJ are calculated by

H IJ = 1 2 (H II + H JJ + N I V I + N J V J )S IJ - 1 2 (V I Φ I | PI |Φ J +V J Φ I | PJ |Φ J ) (13)
In the case of degenerate systems, one can also include more than one configuration to represent the charge localization on a given fragment, as it will be shown in the applications of Section 2.4.

Solving Eq. ( 12) provides both the ground state of the system and some excited states generated via charge resonance. Although these excited states are of interest, for instance in spectroscopy, we focus in this work on the ground state which corresponds to the lowest eigenvalue E g .

Analytical gradient

Calculating derivatives of the energy upon atomic nuclei is required to perform molecular dynamics or geometry optimization. Their numerical calculation is possible by finite differences but the number of energy calculations (2 × 3N atoms ) turns out to be quite large even for small systems.

In SCC-DFTB, it is convenient to use the derivatives of the matrix elements (H 0 , S, Γ) which are known and tabulated. Differentiating Eq. ( 12) with respect to the nuclear coordinate R a of atom a leads to the force expression

∇ a E g = ∑ IJ b I b J ∇ a H IJ -E g ∇ a S IJ ( 14
)
We now present the calculation of the derivatives for the diagonal and off-diagonal elements separately. The differentiation operator ∇ a is replaced by the symbol ∂ a to simplify the expressions.

Derivative of the diagonal element

The diagonal element H II is the energy of the configuration Φ I . Differentiating (Eq. Eq. ( 3)) and using the eigenvalue equation (Eq. ( 4)), the molecular charge conservation (Eq. ( 7)) and orthonormality constraints leads to the analytical

∂ a H II = ∂ a E rep + ∑ i n i ∑ µν c iν c iµ ∂ a H 0 µν +V I ∂ a P I µν + ( H 1 µν S µν -ε i )∂ a S µν + q a ∑ b ∂ a Γ ab q b (15)

Derivative of the off-diagonal elements

The differentiation of the off-diagonal elements is obtained by differentiating Eq. ( 13) to give

∂ a H IJ;I =J = 1 2 A IJ + 1 2 A JI with A IJ = (∂ a H II + N I ∂ a V I )S IJ + (H II + N I V I )∂ a S IJ -Φ I |P I |Φ J ∂ a V I -V I ∂ a Φ I |P I |Φ J ( 16
)
We must now express the derivatives of the Lagrange multipliers ∂ a V I and ∂ a V J , as well as the overlaps (real overlap and through the projectors) between Φ I and Φ J . As there is no relationship between the MO of the two configurations, there is no Hellman-Feynman type simplification for the derivatives of their overlaps. The derivatives of the orbital coefficients and of the Lagrange multipliers must be explicitely calculated. The calculation of the derivatives of the coefficients has already been expressed for DFT (see for instance 39 ) and its expression only differs here through the term containing the constraint.

For a given configuration Φ I , the derivative of the coefficients of the orbitals {φ I i } can be related to the orbitals themselves by a u matrix

∂ a c I iµ = ∑ k c I kµ u ik ( 17
)
The conservation of normalized MO already imposes the form of the diagonal term of the u matrix

u ii = - 1 2 ∑ µν c I iµ c I iν ∂ a S µν . (18)
For the off-diagonal elements u i j;i = j , differentiating Eq. ( 4) leads to

u i j = ∂ a H i j -ε j ∂ a S i j ε j -ε i ( 19
)
where ∂ a S and ∂ a H are the derivatives of the SCC-DFTB overlap and Hamiltonian matrices expressed in the molecular orbital basis set

∂ a S i j = ∑ µν c I iµ c I jν ∂ a S µν ∂ a H i j = ∑ µν c I iµ c I jν ∂ a H µν (20)
In the constrained SCC-DFTB, the Hamiltonian matrix derivative depends -i-on the derivatives of matrices H 0 , S, Γ and P -ii-on the derivatives of the coefficients and -iii-on the derivatives of the Lagrange multipliers. These three contributions are now explicitely separated

∂ a H i j = ∂ a F i j + ∑ kl A i j,kl u kl + ∂ a V I P i j ( 21
)
where ∂ a F i j contains the first contribution

∂ a F i j = ∑ µν c I iµ c I jν ∂ a F µν (22)
with

∂ a F µ∈α,ν∈β = ∂ a H 0 µν +V I ∂ a P µν + ∂ a S µν H 1 µν S µν + 1 2 S µν ∑ ξ (∂ a Γ αξ + ∂ a Γ ξ β )q ξ + ∑ i n i ∑ l ∑ ω∈ξ (Γ αξ + Γ ξ β )c iω c il ∂ a S ωl
The second term in Eq. ( 21) accounts for the Hamiltonian dependences on the orbital coefficients with

A i j,kl = ∑ µν c I iµ c I jν ∑ ω ∂ H µν ∂ c I kω c I lω ( 23
)
and

∂ H µν ∂ c I kω = 1 2 n k S µν ∑ ξ ∑ λ ∈ξ S λ ω (Γ αγ + Γ β γ + Γ αξ + Γ β ξ )c I kλ ( 24
)
where µ ∈ α; ν ∈ β ; ω ∈ γ. In the last term of Eq. ( 21), P i j accounts for the Hamiltonian differ-entiation upon the Lagrange multiplier

P i j = ∑ µν c I iµ c I jν P µν (25)
Compacting the i j indices in a single m indice and the kl indices in a single n indice, we now define

u m = u i j v m = ∂ a F i j -ε j ∂ a S i j ε j -ε i B mn = A i j,kl ε j -ε i w m = P i j ε j -ε i
and rewrite Eq. ( 19)

u = Bu + v + ∂ a V I w = (1 -B) -1 v + ∂ a V I (1 -B) -1 w = u 0 + ∂ a V I u ′ (26)
with

u 0 = (1 -B) -1 v and u ′ = (1 -B) -1 w
We now determine ∂ a V I using the fact that N I remains constant. Differentiating Eq. ( 7) leads to

∂ a N I = ∑ i n i ∑ µν c iµ c iν ∂ a P µν + 2∂ a c iµ c iν P µν = 0 (27)
which can be expressed with the u matrix

∂ a N I = ∑ i ∑ µν n i c iµ c iν ∂ a P µν + 2 ∑ i j n i u i j P i j = 0 ( 28
)
Using the previous expression for u leads to the following expression for the derivatives of the Lagrange multiplier

∂ a V I = - ∑ i ∑ µν n i c iµ c iν ∂ a P µν + 2 ∑ i j n i u 0 i j P i j 2 ∑ i j n i u ′ i j P i j (29)
H can now be calculated from Eq. ( 21), as well as the u matrix from Eq. ( 19) and finally the derivatives of the coefficients from Eq. ( 17). The S IJ derivatives are computed from this MO coefficients derivatives and the derivatives of the AO overlap matrix. For efficiency, the determinant expansions appearing in the calculation of the derivatives were calculated using the Sherman-Morrison formula. ? A similar approach is applied to the derivatives of the projected overlap matrix

Φ I |P I |Φ J .

Variants of the DFTB-VBCI: the HOMO approximation

We will consider the following approximation to the DFTB-VBCI approach: we assume that, in a molecular cluster, the MOs of the different charge localized configurations mostly differ through their Highest Occupied Molecular Orbital (HOMO). The overlaps and projected overlaps between two configurations can be simplified in

S IJ = Φ I |Φ J ≃ φ I HOMO |φ J HOMO ( 30
)
Φ I |P I |Φ J ≃ N I φ I HOMO |φ J HOMO + φ I HOMO |P I |φ J

HOMO

The off-diagonal CI matrix element becomes

H IJ ≃ 1 2 (H I + H J ) φ I HOMO |φ J HOMO (31) - 1 2 (V I φ I HOMO |P I |φ J HOMO +V J φ I HOMO |P J |φ J HOMO )
The advantage of this approach is to avoid any Slater determinant overlap calculation and only the derivatives of the HOMO coefficients need to be calculated. Following 41 , we call π a g and π o g the degenerate MO in the neutral benzene molecule. In the ionized monomer, these two levels are degenerate at the neutral geometry but undergo Jahn-Teller distortion leading to an acute angle configuration (ionization from the π a g orbital) or an obtuse angle configuration (ionization from the π o g orbital). We use the qualitative molecular description of the Dimer Molecular Orbitals Linear Combination from the Fragment Molecular Orbitals (DMO-LCFMO 41 ), to describe a benzene dimer, composed of two fragments labelled A and B. In this framework, the constrained state A + B can be obtained either from removing one electron from the π a g or from the π o g orbitals of A. Consequently, we need two configurations to describe the constrained form A + B. In the case of the symmetric D 6h sandwich stacked dimer, these two A + B configurations are degenerate and are built as follows: we use the constrained SCC-DFTB to minimize the electronic energy with an occupation of 1.5 for the two highest occupied MOs (HOMO and HOMO-1) orbitals and 2 for the lower electron energy orbitals. The two A + B configurations are then built from the obtained MOs as shown on Figure 2. The same procedure is applied to obtain two AB + configurations and the CI matrix to be diagonalized is a 4 × 4 matrix.

In the other isomers (displaced sandwiches and T-shaped), the π a g and π o g orbitals of each fragments are no longer degenerate and one could in principle calculate the energy of these configura-Figure 2: Two electronic configurations (right) obtained from constrained SCC-DFTB calculation with non integer occupation numbers (left).

tions without using fractional occupation numbers. However, we could not obtain a self-consistent solution of the A + B state (respectively AB + ) with the constrained SCC-DFTB since the π a g and π o g orbitals on fragment A (respectively on B), although not degenerate remain close in energy.

We therefore decided to keep the procedure used for the D 6h stacked sandwich isomer filling the HOMO and HOMO-1 orbitals with 1.5 electrons. Although the filling of MO is fixed, these MO relax anyway and are no longer degenerate in the final results, due to the coupling with geometry relaxation. We first discuss the results obtained for the stacked sandwich in the D 6h geometry. Figure 3 represents the energy of the dimer corresponding to dissociation along the z-axis, orthogonal to the planes of the monomers. For this example, the fragments are frozen at the monomer neutral geometry). The zero energy reference corresponds to the sum of the separated fragments calculated independently namely E(C 6 H + 6 ) + E(C 6 H 6 ). SCC-DFTB dissociation curve reminds the dis- sociation curve of H + 2 calculated with DFT and can be explained as follows. First, the SCC-DFTB energy does not converge to the sum of the energies of the fragments at the dissociation limit. At infinite distance, the charge is equally distributed over the two fragments. As the evolution of the self-interaction error with the number of electrons on a fragment is unfortunately not constant nor linear, we have 2 × E SCC-DFT B (C 6 H 0.5 6 ) = E SCC-DFTB (C 6 H + 6 ) + E SCC-DFTB (C 6 H 6 ). At shorter distances, the energy increases and a barrier is even observed before reaching the minimum which is here a metastable minimum. The responsible repulsive contribution has a 1/R behavior and can be attributed to the artificial repulsion of two half charged fragment, which is a different case of self-interaction than the previous one. Finally, the minimum is much too low in energy (29.5 kcal/mol) as compared to the reference calculations (see Table 1).

Sandwich isomer

As can be seen from Figure 3, the DFTB-VBCI method does not present the wrong behavior pattern of the SCC-DFTB curve. At the dissociation limit, the energy converges to the sum of the energies of the fragments. Actually, in Eq. ( 12) the overlaps and coupling terms vanish and the energies of the localized configurations are degenerate. These energies are calculated with the electronic density corresponding to one charged and one neutral monomer and not that of two half charged fragments. The Coulombic self-interaction 1/R repulsion also disappears with this approach as well as the corresponding barrier. Finally, the binding energy for the stacked sandwich (17.70 kcal/mol) is strongly smaller than the SCC-DFTB one (29.53 kcal/mol) and gives a much better agreement with that of EOM-IP-CCSD calculation at 18.34 kcal/mol.

The dissociation curve obtained applying the HOMO approximation detailed in Section 2.5 is also plotted on Figure 3. It is almost identical to the DFTB-VBCI curve with a binding energy of 17.91 kcal/mol vs 17.70 kcal/mol for the DFTB-VBCI.

T-shaped and displaced sandwich

The T-shaped and displaced sandwiches have been optimized without any geometrical constraint.

The binding energies are reported in Table 1. The T-shaped isomer has a binding energy of 9.23 kcal/mol which is slightly smaller than the EOM-IP-CCSD one. Another difference concerns the charge localization. With EOM-IP-CCSD, the charge is mostly localized on the stem fragment (88 %) whereas this localization drops to 56 % with the DFTB-VBCI. A possible explanation for this charge localization discrepancy could be related to lack of stabilization by polarization.

In DFTB-VBCI, the benzene π system can be polarized in the direction parallel to the benzene ring. However, due to the reduced basis set used, the polarization of the π system perpendicular to the benzene ring cannot be accounted for. This lack of polarization could be at the origin of the destabilization of the configuration where the charge is carried by the stem fragment, leading to an over-sharing of the charge and an underestimation of the binding energy. Another possibility is the definition of charges of the charge analysis in the ab initio calculation (NBO).

At the SCC-DFTB level, the x-and y-displaced dimers are overstabilized as compared to reference calculations. Similarly to what is observed for the stacked sandwich dimer the DFTB-VBCI approach gives considerably improved binding energies (20.90 and 21.26 kcal/mol) in very good agreements with those of EOM-IP-CCSD (20.14 and 20.18 kcal/mol). These two structures are almost degenerate, the y-displaced being slightly more stable. The energy difference is however probably much smaller than the expected accuracy of our approach. These are clearly cases in which quantum vibrational effects should be considered.

Applying the HOMO approximation to the DFTB-VBCI leads to very similar results. The most stable structures are the x-and y-parallel displaced ones with binding energies which differ by less than 2.5 % from those of DFTB-VBCI. The T-shaped is found to be less stable than the previous isomers but its binding energy is overestimated compared to reference calculations and DFTB-VBCI, suggesting that the differences between the two configurations can not be reduced to a change in the HOMO.

Water dimer cation

The potential energy surface of cationic water dimers has been investigated at at hight level of theories CCSD(T) ) and EOM-IP-CCSDT [52][53][54][55] . The stable structures belong to two families. In the first one, the two water monomers are superimposed in an antisymmetric pattern and the charge is equally distributed over the two units. The second one results from a proton transfer leading to two non-symmetric units [H 3 O -OH] + , in which the charge is mostly localized on the H 3 O fragment.

The structures, found to be minima by Cheng et al. 54 , have been optimized with DFTB-VBCI method within a conjugated gradient scheme. In this system, the DFTB-VBCI considers the interaction between the configurations where the positive charge is localized either on the H 3 O or on the OH subfragments. We started the The weights of the two configurations in the CI approach indicate that the charge is mostly localized (99.9 %) on the H 3 O fragment. In CCSD(T) calculations, the charge is localized also strongly localized on this fragment but only at %?????????????????????.

We notice that if the DFTB-VBCI minimum (C s -trans) is different from the C1 minimum obtained with CCSD(T). In CCSD(T), the C s -trans isomer corresponds to a transition state 0.04 kcal/mol higher in energy than the global C 1 minimum 54 (< 0.2 kcal mol -1 for the EOM-IP-CCSD fitted surface of Kamarchik et al. 59 , and 0.03 kcal/mol at the SCF+MCPF level 57 ). Such a small energy difference is far beyond the expected accuracy of the DFTB-VBCI method. The binding energy of C s -trans isomer is close (42.7 kcal/mol) to that obtained with a simple SCC-DFTB calculation (46.3 kcal/mol). This is due to the fact that the charge is not significantly delocalized between the two fragments and that the SCC-DFTB calculation already attributes most of the charge to the H 3 O fragment. The artificial stabilization by the self-interaction error is therefore less crucial. This also explains why most of the SCC-DFTB functionals give reasonable results for this structure.

Antisymmetric structure

It can be seen from Section 3.2 that, at the CCSD(T) level, the antisymmetric isomer is less sta-

Conclusion

An extended method combining a VBCI-like scheme with the SCC-DFTB has been developed.

The method has been implemented together with its analytical gradient to make possible complete optimization, including namely the intra-molecular and inter-molecular degrees of freedom. We have benchmarked this approach on ionized dimers of benzene and water.

In the benzene dimer cation, the self-interaction error is at the origin of an unphysical behavior of the SCC-DFTB dissociation energy curve. It turns to be fully corrected with the DFTB-VBCI as detailed for the stacked sandwich. This method will allow to perform explorations of potential energy surfaces by molecular dynamics or Monte Carlo sampling with the focus of studying ionized dimers dissociation. The binding energies obtained for different isomers with the DFTB-VBCI compare well with those of high-level calculations as well as experimental data, these standard energies being strongly overestimated with the SCC-DFTB. We however notice that the main error for the DFTB-VBCI binding energy concerns the T-shaped structure which is understabilized by 3 kcal/mol. This may be due to the use of point charges and a possible mistreatment of the multipolar nature of the benzene π system interacting with that of the charged top benzene. Further improvement of the DFTB-VBCI could include such multipolar description of the π system (see for instance 60 ) in order to account for this effect but at the price of a larger computational effort to derive the energy gradient.

The second benchmark system is the ionized water dimer. The two main isomers strongly differ the most stable one. This effect is corrected with the DFTB-VBCI approach which gives a binding energy close to that of high level of calculations.

Taking advantage of the SCC-DFTB in terms of computational efficiency, the DFTB-VBCI will allow to deal with systems larger than dimers. In our previous study, DFTB-VBCI has been used to characterize binding energies, ionization potentials as well as charge localization in stacked coronene clusters with frozen intramolecular geometries and equal spacings between the planes of the units 35 . This work was however performed before the development of analytical gradients and it will be of interest to characterize the effects of intra-and inter-molecular relaxation in these clusters. It is known that in rare gas clusters He + n , Ne + n , Ar + n , Kr + n and Xe + n , the hole tends to delocalize on few units (from 2 to 4, depending the rare gas) and that the other atoms tend to organize in crowns around a linear core. It will be intersesting to understand how the internal degrees of freedom, the molecular extension and shape influence the size of the core unit and the organization of the whole cluster.

• Pour l'eau : Decrire notre choix de remplissage electronique H3O-OH

• Pour l'eau : dans la litterature la valeur de la charge sur H3O donnee en DOHF et excitation seulement ?

• Pour le sandwich stacked la courbe ne correspond pas aux valeurs du tableau !!!! refaire avec la meme dispersion

• parler des comparaison avec la SIC DFT (pienazcek)

• pour le water H3O-OH charge localization Cheng donne les valeurs ROHF ....

• faire le calcul H3O-OH avec le 0.5 pr le OH

• mettre 1 ref ancienne pour le coupled perurbed equations

• formule du homo a approx

• figure 3 ecrire SCC-DFTB

• ecrire les gradient dans la conclusion

• gradient complet + colle a l'experience

Figure 1 :

 1 Figure 1: Benzene dimer cations optimized at the DFTB-VBCI level. (a) T-shaped, (b) sandwich stacked, (c) x-displaced and (d) y-displaced isomers.

Figure 3 :

 3 Figure 3: Dissociation potential energy curves of the cationic benzene dimer in the stacked sandwich configuration calculated with SCC-DFTB and DFTB-VBCI approaches.

Figure 4 :

 4 Figure 4: Water dimer cations optimized at the DFTB-VBCI level. (a) [H 3 O-OH] 2 isomer and (b) antisymmetric or [H 2 O-H 2 O] + isomer.

  by a proton transfer. The binding energy of [H 3 O -OH] + isomers calculated at the DFT level with several functionals are in good agreement with reference calculations. This is also the case with SCC-DFTB and DFTB-VBCI due the localization of the charge on the H 3 O subfragment reducing the multiconfigurational nature of the wavefunction and the self-interaction error in standard DFT based calculations. On the opposite, in the [H 2 O -H 2 O] + isomer the charge is equally carried by the two fragments and the binding energies obtained at DFT level strongly differ depending on the choice of the functional. With SCC-DFTB, this structure is overstabilized and becomes artificially

Table 1 :

 1 Binding energies of cationic benzene dimer obtained at different levels of theory. The stacked sandwich structure correspond to constrained D 6h optimization whereas the other isomers are fully optimized with its respective method. Refs: a Pieniazek et al.Pieniazek et al.41 ; b Pieniazek et al.42 ; c Ibrahim et al.44 ; d Itagaki et al.43 ; e Kryachko 45 ; f Miyoshi et al.40

		DFTB-VBCI	HOMO	SCC-DFTB EOM-IP-CCSD	DFT	CASSCF+MCC
			Approximation SCC-DFTB EOM-IP-CCSD	DFT	CASSCF+MCC
	Stacked	17.70	17.91	29.53	18.34 a	17.4 c -19.1 d
	Sandwich
	x-displaced	20.90	20.43	29.01	19.58 b	-	12.3 f
	y-displaced	21.26	20.79	29.21	19.81 b	16.57 e	10.9 f
	T-shaped	9.23	16.90	24.68	12.41 b	15.7 c

Table 2 :

 2 Binding energies of cationic water dimer obtained at different levels of theory. a Lee and KimLee and Kim 55 ; b Cheng et al.Cheng et al. 54 ; c Barnett and LandmanBarnett and Landman 56 ; d Sodupe et al. 57 MCPF = SCF + electron correlation included with size extensive Modified-Coupled-Pair Functional; e Gill and Radom 58 optimization from the C 1 global CCSD(T) minimum. The absence of symmetry of this minimum leads to non-degenerate HOMOs. For the two charge-localized configurations. The DFTB-VBCI optimization leads however to a Cs structure. Due to the planar symmetry, degeneracies must

		AntiSymmetric H 3 O-OH (Cs) H 3 O-OH
	DFTB-VBCI	35.44	42.31	-
	HOMO Approx.	35.74	42.33	-
	SCC-DFTB	68	46.3	-
	GGC	53.73 c	48.66 c	-
	BLYP	58.4 a	-	49.3 a
	B3LYP	51.5 a	-	49.8 a
	MPW1K	42.9 a	-	49.9 a
	BH&HLYP	41.4 a	-	49.9 a
	MP2	40.48 a / 43.5 e	50.9 e	46.47 a
	MP4	41.1 e	49.9 e	-
	CCSD(T)	39.53 a / 39.59 b	46.64 b	46.70 a / 46.68 b
	MCPF d	36.1	45.9	45.93 d

be considered in the detailler .... ??????????????

  ble by 7 kcal/mol than the [H 3 O -OH] + isomer. At the DFT level, the binding energy strongly depends on the choice of the functional. For instance, the [H 3 O -OH] + structure is more stable than [H 3 O -OH] + with MPW1K, BH&H and BH&LYP functionals but it is the opposite with the BLYP, BPW91, HCTH407 and B3LYP functionals. At the SCC-DFTB level, the binding energy of the OO isomer is strongly overestimated (68 versus 39 for CCSD(T)) making this isomer 22 kcal/mol more stable than the [H 3 O -OH] + isomer. The DFTB-VBCI has then a drastic effect, reducing the binding energy to 38 kcal/mol, a value close to the CCSD(T) ones (35-39 kcal/mol).

In this isomer, the charge is equally distributed between the two equivalent fragment. The overstabilization observed at the SCC-DFTB level is attributed to the self-interaction error due to this this strong delocalization and corrected by the DFTB-VBCI approach. Finally, we notice that for both the [H 3 O -OH] + and [H 2 O -H 2 O ] + isomers, the binding energies obtained with the HOMO approximation are very close to that obtained with the full DFTB-VBCI method.

Applications

We will now apply the DFTB-VBCI method to the treatment of two prototype cationic molecular clusters, namely the benzene dimer and the water dimer.

Benzene Dimer

Several authors have investigated benzene dimer clusters at high levels of theories. Let us cite for instance MultiReference Configuration Interaction including Single and Double excitations (MR-CISD 40 ) and Equation-Of-Motion Coupled-Cluster model with Single and Double substitutions for Ionized systems (EOM-IP-CCSD 41,42 ) discussing the relative energies of characteristic isomers, namely the sandwiches (stacked, parallel x-and y-displaced) and T-shaped configurations (see Figure 1). Both calculations gave the sandwich parallel displaced isomers to be the most stable structures, about 7-8 kcal mol -1 lower than the T-shaped. Both approaches can describe in principle correctly the charge resonance states. However MR-CISD may suffer from the lack of size consistence which could explain that the binding energies (12.3 kcal mol -1 for the displaced sandwich and 5.4 kcal mol -1 for the T-shaped) are significantly smaller than those obtained with EOM-IP-CCSD (20.2 kcal mol -1 and 12.4 kcal mol -1 ) whereas the energy difference between the two isomers is in good agreement (about 7-8 kcal mol -1 respectively) ??WARNING cest du MRCPA size consistent??. The D 6h symmetric stacking is another structure of interest which is slightly less stable (about 2 kcal mol -1 ) than the two sandwich displaced isomers. DFT calculations [43][44][45] performed with the B3LYP functional give reasonable binding energies for the sandwhich structure (17-19 kcal mol -1 ) but underestimate the energy difference between the two structures due to an over-stabilization of the T-shaped structure. In the following, we will use the EOM-IP-CCSD 41,42 results as references to benchmark our method as they are the most recent ab initio calculations and the corresponding binding energies for the most stable structures around 20 kcal/mol (see Table 1 are in good agreement with the experimental studies providing values in the 15-20 kcal/mol range 44,46-51 .


A good example of points 3 and 4 is connecting flights at the airport: Interaction between public and private transports, getting a lot of (local) informations, managing critical situations, etc..

HPC and Data centers (clouds)

2 Fog -small clusters (hundreds/thousands cores)

3 Edge -laptops/smart phones 4 IoT -sensors A good example of points 3 and 4 is connecting flights at the airport: Interaction between public and private transports, getting a lot of (local) informations, managing critical situations, etc.. Manage efficiently such complex infrastructures composed of multiple and heterogeneous digital/computing components.

Large amount of polymorphic data issued from multiple sources.

Resources appear and disappear.

New storage capabilities.

Heterogeneity characteristics of the different architectures, including OS kernels, compilers, etc.

Network of different types and instability (different latencies).

WHERE running each workflow component? Local scheduler (model)

Two sets of tasks:

S 1 global tasks. off-line (sequential) tasks with their release dates (corresponding to data transfers).

S 2 local tasks. Sequential tasks arriving on-line. Not too many (10 to 20 %), usually small ones.

Qarnot platform is rather homogeneous (uniform machines).

Notice here that rejection (interrupt a task and re-execute it elsewhere) makes sense.

Two directions are under investigation.

Single objective.

Minimize max flow time for S 2 under the capping on sum flow for S 1 .

Bi-objective.

Minimize average flow time or tardiness for S 1 and minimize max flow time for S 2

  Report a 3 years project (french nat. Research Agency) with the Qarnot Company Describe the problems Show and discuss preliminary scheduling solutions Simulation tools 1 1 Ideally, build a digital twin of the platform From Heaters to Computing Resources: "turn IT waste heat into a viable heating solution for buildings." The actual Qarnot platform (Bordeaux and Paris): ∼1,000 distributed QRads embedding ∼3,000 diskless computing units ∼20 local servers (QBoxes) with memory disks 1 global server (QNode) with a centralized storage server Main goal: design, implement and test different placement and scheduling policies at both QNode-and QBox-levels. sec.) from QBox to QNode with: Number of resources available for each task priority Amount of free space on disk QNode-scheduling: Sort QBoxes by least available resources first (no temperature knowledge) Sort tasks by highest priority first For each task, dispatch as many instances as possible QBox-scheduling: Retrieve data-set dependencies Schedule high priority instances on coolest QRads Schedule low priority instances on warmest QRads

https://github.com/simgrid/simgrid

https://gitlab.inria.fr/batsim/batsim/tree/temperature

https://gitlab.inria.fr/batsim/pybatsim/tree/temperature

Challenges synthesis

Some issues/questions: Compute as close as possible to the sensors/digital components to avoid data transfers (and also for privacy).

Need of new data storage abstraction: object versus file systems data reduction/lossy compression. Task flow versus data flow: data streaming. For the moment, no local tasks are considered.

Two types of computing requests

Main problem: testing on the production platform is not conceivable (and it would take time...).

→ Simulation is needed!

Remark: the theoretical analysis is difficult.

Main problem: testing on the production platform is not conceivable (and it would take time...).

→ Simulation is needed!

Remark: the theoretical analysis is difficult.

Simulation Framework

SimGrid 3 : Large-scale distributed system simulator with execution and communication models. → Used to simulate platform (including temperature) and tasks execution.

Batsim 4 : Infrastructure simulator for jobs and I/O scheduling. → Used to drive the simulation, submit tasks and communicate with the decision process.

Pybatsim 5 : Batsim's Python API exposing methods to easily communicate with the Batsim process. → Used to implement both QNode and QBox schedulers.

Implementation

Dedicated module for convert computations to temperature.

solutions have been developed:

Learning on actual logs (problems: old generation of Qrads and missing contextual informations).

Analytic expression (given by physicians).

Remark: better than what is available today at Qarnot (empirical law which causes a lot of instabilities).

Simulations

Run 4 times the logs on a week (between may 3 and june 1, 2019).

The simulated platform was composed of 3390 QMobos of 669 QRads, managed by 20 QBoxes.

Two metrics studied: mean waiting time bounded slowdown (normalized average stretch)

Some results

Qarnot platform.

Preliminary results assuming periodic arrivals in S 2 : Study the problem for one machine, then extend.

Interface between global and local schedulers: add a module which gives the pressure of local tasks.

Learning algorithms have been developed (get a better idea of the leaving habits).

Extend to more general edge platforms.

Conclusion

General message: It is not so easy to build a bridge between theory and practical tools.

Model of the actual Qarnot/edge platform.

Simulator (almost complete).

We proposed several scheduling problems . Some are still under investigation Theoretical analysis remains to be done A preliminary solution has been implemented and assessed on real execution traces.


Delphine Fradin

Clémence Mille

Nadia Naoui

Chris Groves

Mark I Mccarthy

Samy Hadjadj

François Paeou

Pierre Bougnères

Allele Dependent DNA Methylation of the Human Insulin Gene Promoter in T2D patients and Controls

Introduc)on

Numerous SNPs located at various gene loci related to insulin secreVon have been idenVfied by GWAS of T2D or of insulin--related traits (review in Bonnefondd et al. 2010). GeneVcs however seems to explain yet a limited part of these phenotypes (Wheeler et al. 2011). Non geneVc factors are certainly of utmost importance, and possibly include epigeneVc factors at the interface of geneVcs and environment (Chen et al. 2011).

Genotype dependent variaVons of neighboring epigenotype have been reported at several other loci across the genome and called ASM (Allele Specific DNA MethylaVon). Shoemaker et al. showed that 23% to 37% of heterozygous SNPs are associated with epigeneVc variaVon in several cell line (Shoemaker et al. 2010) while Zhang et al. found that ASM is likely to affect 10% of all human genes (Zhang et al. 2009). The finding of a link between geneVc variants and epigenomic marks may have major consequences for a beEer understanding of the "missing heritability".

Several studies have established a significant associaVon between the class III VNTR, in complete linkage disequilibrium with the allele T of rs689, and T2D (On et al. 1999, PermuE et al. 1990). Previous studies have found that INS gene expression is regulated by epigeneVc mechanisms (ChakrabarV et al. 2003, Mutskov et al. 2007et 2009) such as DNA methylaVon (Kuorda et al. 2009, Yang et al. 2010). We also found an associaVon between DNA methylaVon of INS promoter and Type 1 Diabetes (Fradin et al. 2012).

The INS proximal promoter encompasses 7 CpGs located at posiVon --234, --206, --180, --135, --69 and --19bp relaVve to the TSS, too sparse to be considered a CpG island, thus the INS gene promoter is considered as a low CpG promoter. Three CpG residues are important to the binding of regulatory transcripVon factors: the --234 CpG site is part of the E2 box, --206 and --180 CpG sites are parts of cAMP responsive element (CRE) 1 and 2 respecVvely (Figure 1).

We performed geneVc and epigeneVc analysis of the INS promoter in whole blood cell (WBC) from two T2D case--control cohorts to study the relaVonship between DNA methylaVon levels and rs689 genotypes.

Results

Influence of cis genotype upon methyla4on.

DNA methylaVon at CpG --69, --102, --180, --206 and --234 were found to be associated with rs689 genotypes in T2D paVents (N=491) and controls (N=361) from France (Figure 2). We replicated this finding in a second T2D case--control cohort from the UK where methylaVon at CpGs --69, --180 and --206 showed a comparable associaVon with rs689 genotypes.

Materials and Methods

Pa4ents

Isola4on of genomic DNA and bisulfite genomic conversion and pyrosequencing. Nucleic acids were extracted from WBC or Vssue using

phenol/chloroform procedure. Genomic DNA was treated with EZ--96 DNA MethylaVon--Gold Kit, according to manufacturer's protocol (Zymo Research CorporaVon). We PCR--amplified the bisulfite treated genomic DNA using unbiased nested primers (sequences on request) and performed quanVtaVve pyrosequencing (PyroMark Q96 MD, Qiagen). The percentage of methylaVon for each CpG site was determined for each subject using the Q--CpG methylaVon sooware (Qiagen) according to the standards in Tost and Gut (2001).

Rs689 genotyping. Rs689 was genotyped by restricVon fragment length polymorphism as previously described (Le Stunff et al. 2001).

Sta4s4cal analysis. Differences in DNA methylaVon of the insulin promoter between rs689 genotypes were analyzed using non--parametric Kruskal--Wallis test. CorrelaVons were calculated as adjusted R square that measures the proporVon of the variaVon in the dependent variable accounted for by the explanatory variables. Differences in DNA methylaVon of the insulin promoter between T2D paVents and non--diabeVc controls were analyzed using non--parametric Wilcoxon rank sum test. All staVsVcal analysis were conducted using R 2.10.1.

Conclusion

!""# $"# %"# &"# '"# "# ()*#+!,# # #+%, # ########+!"'# # ##+!-.# # ##########+!$"# # #####+'"%# # # # +'-&# AA AT TT
//#####/0#######00# //#####/0#######00# //#####/0#######00# //#####/0#######00# //#####/0#######00# //#####/0#######00# //#####/0#######00# !""# $"# %"# &"# '"# "# 12/#3456789:;<#=>?# )@$A!" +!% # )@-A!" +!! # )@'A!" +!! # )@-A!" +!! # )@'A!" +& #
)@.A!" +% # )@!A!" +. # )@.A!" +& # )@-A!" +!" # )@%A!" +. # BC%$,#D4<;57)4C# !"

#!" $!" %!" &!" '!!" '#!" !"#$%&'()*+,-.+ !""# $"# %"# &"# '"# "# #()*#+!,#########+%,############+!"'#############+!-. ####+!$"## #+'"% #####+'-&# /01

Methyla4on has a comparable paHern in whole blood cells and other 4ssues

Comparable methylaVon levels and paEern were found in WBC, liver, skin, peritoneum and intesVne (Figure 3). In contrast, islets of Langerhans, while following a comparable paEern of methylaVon had lower levels of methylaVon than WBC at all studied CpGs. !"#$%&'$ *)$$$$.)$$$$$()$$$$-)$ *)$$$$.)$$$$$()$$$$-)$ *)$$$$.)$$$$$()$$$$-)$ *)$$$$.)$$$$$()$$$$-)$ *)$$$$.)$$$$$()$$$$-)$ *)$$$$.)$$$$$()$$$$-)$

French cohort T2D cases

$ $%(' $ $$$$$$$$%&)*$ $ $$%&+,$ $ $$$$$$$$$$%&-)$ $ $$$$$%*)($ $ $ $ %*+.$ *)$$$$.)$$$$$()$$$$-)$ /01$2314567$ 58)9&,$ :8&9&) %' $ 58)9&;$ "8*9&) %&+ $ 58)9),$ "8)9)&$ &))$ -)$ ()$ .)$ *)$ )$ <=/$>1?@3A4BCD$2E7$ 58)9&*$ "8,9&) %; $ 20

Intergroup varia4on of methyla4on.

We found no T2D--related global direcVonal change in methylaVon level that would affect all CpGs, nor did neighbouring CpGs adopt systemaVc changes in methylaVon in T2D cases or controls. Consistently, CpG --19 showed a decreased DNA methylaVon in T2D cases compared with non--diabeVc controls (p=2.10 --16 ) in French and UK cohorts (Table 1).

  . Non--immortalized WBC samples were taken from parVcipants of the Type 2 diabetes cohorts from PoiVers (T2D cases) and from parVcipants of the MONA LISA Lille study (controls). A sample of 132 non obese T2D paVents aged less than 65 years, studied in mean 7 yrs aoer T2D diagnosis, as well as 186 age--and BMI--matched controls was provided by M. McCarthy. Tissues. Human liver (n=14), peritoneum (n=26), skin (n=4), intesVne (n=2) and islets of Langerhans (n=4) were provided by F. PaEou. Each Vssue comes from a different person.

Figure 1 .

 1 Figure 1. Epigne0cs and gene0c of INS locus.

Figure 2 .

 2 Figure 2. Boxplot of methyla0on values of the studied CpGs in the INS promoter with respect to the rs689 genotype in T2D pa0ents (top) and controls (boKom), Kruskal Wallis test.

Figure 4 .

 4 Figure 4. Correla0on between age and INS promoter methyla0on in our two T2D case--control cohorts.

Table 1 .

 1 CpG methyla0on in the INS promoter in our two case--control cohorts. Results are expressed as mean ± sd. P--values are calculated by Wilcoxon rank sum test.


Thanh-An Pham

Emmanuel Soubies

Ahmed Ayoub

Demetri Psaltis

Michael Unser

ADAPTIVE REGULARIZATION FOR THREE-DIMENSIONAL OPTICAL DIFFRACTION TOMOGRAPHY

Keywords: plug-and-play, nonlinear inverse problems, dictionary learning, computational imaging

Optical diffraction tomography (ODT) allows one to quantitatively measure the distribution of the refractive index of the sample. It relies on the resolution of an inverse scattering problem. Due to the limited range of views as well as optical aberrations and speckle noise, the quality of ODT reconstructions is usually better in lateral planes than in the axial direction. In this work, we propose an adaptive regularization to mitigate this issue. We first learn a dictionary from the lateral planes of an initial reconstruction that is obtained with a total-variation regularization. This dictionary is then used to enhance both the lateral and axial planes within a final reconstruction step. The proposed pipeline is validated on real data using an accurate nonlinear forward model. Comparisons with standard reconstructions are provided to show the benefit of the proposed framework.

INTRODUCTION

Optical diffraction tomography (ODT) is a microscopic technique that provides a three-dimensional map of the refractive index (RI) of the sample [START_REF] Wolf | Three-dimensional structure determination of semitransparent objects from holographic data[END_REF]. It proceeds by measuring the complex fields that are produced when the sample is illuminated with different tilted plane waves. These measurements are related to the RI of the sample through the wave equation that governs the light-scattering phenomenon. This allows for the deployment of numerical methods to recover the RI.

The reconstruction of the RI consists in solving an inverse scattering problem which is particularly challenging for several reasons. First, wave scattering is a nonlinear phenomenon and the most faithful forward model for ODT is iterative [START_REF] Liu | Seagle: Sparsity-driven image reconstruction under multiple scattering[END_REF][START_REF] Zhang | 2 z3 Rytov z1 = -0.49µm z2 = 0µm z3 =[END_REF][4][START_REF] Pham | Three-dimensional optical diffraction tomography with Lippmann-Schwinger model[END_REF]. Second, the wavevector of the incident wave is limited to a small cone around the optical axis (Figure 1). This leads to the well-known missing-cone problem [START_REF] Lauer | New approach to optical diffraction tomography yielding a vector equation of diffraction tomography and a novel tomographic microscope[END_REF][START_REF] Lim | Comparative study of iterative reconstruction algorithms for missing cone problems in optical diffraction tomography[END_REF] whose effect is an elongation of the reconstructed object along the optical axis. Put simply, the ODT lateral resolution is better than its axial resolution. Finally, the measured fields contain optical aberrations or speckle noise [START_REF] Goodman | Speckle Phenomena in Optics: Theory and Applications[END_REF][START_REF] Che | Reduction of speckle noise in digital holography by combination of averaging several reconstructed images and modified nonlocal means filtering[END_REF] which are usually unknown and constitute a source of model mismatch. To mitigate these problems, former works added prior knowledge in the reconstruction process using regularization [4,[START_REF] Lim | Comparative study of iterative reconstruction algorithms for missing cone problems in optical diffraction tomography[END_REF][START_REF] Sung | Optical diffraction tomography for high resolution live cell imaging[END_REF][START_REF] Kamilov | Optical tomographic image reconstruction based on beam propagation and sparse regularization[END_REF]. These include, for instance, nonnegativity constraints and total-variation (TV) regularization.

In this work, we propose to go one step further and take into account the anisotropic resolution in ODT while improving the quality q } Q q=1 ∈ R 3 of the Q incident plane waves {u in q } Q q=1 are limited to a cone around the optical axis. of the reconstructions. Our motivation is to learn highly resolved features from lateral planes and use them to enhance the quality in the axial direction. Hence, inspired by the strategy proposed by Soulez [START_REF] Soulez | A learn 2D, apply 3D method for 3D deconvolution microscopy[END_REF] for the deconvolution of fluorescent microscopic images, we deploy a dictionary-based regularizer that is learnt from the lateral planes of an initial reconstructed volume.

PHYSICAL MODEL

We denote the RI of the sample by n : Ω → R and we assume that its support is included in the region of interest Ω ⊂ R 3 (Figure 1). The RI of the surrounding medium is denoted by n b . Let the sample be illuminated with the monochromatic incident field u in : Ω → C of wavelength λ. The resulting total field u = u in + u sc ,where u sc : Ω → C is the scattered field, follows the Lippmann-Schwinger (LiSc) model

u(x) = u in (x) + Ω g(x -z)f (z)u(z) dz, (1)
where

f (x) = k 2 b n(x) 2 /n 2 b
-1 is the scattering potential and k b = 2πn b /λ is the wavenumber in the surrounding medium. Finally, under Sommerfeld' radiation condition, the free-space Green function g : Ω → C is defined by [START_REF] Schmalz | On the derivation of the Green's function for the Helmholtz equation using generalized functions[END_REF]

g(x) = exp (jk b x ) 4π x . (2)
Following [START_REF] Liu | Seagle: Sparsity-driven image reconstruction under multiple scattering[END_REF]4], we consider a two-step ODT forward model. We first discretize Ω into N = n 3 voxels. Then, the scattered field y sc ∈ C M at the camera plane Γ is computed as

u = (I -G diag(f )) -1 u in (3)
y sc = G diag(f )u, (4)
where f ∈ R N , u in ∈ C N , and u ∈ C N are sampled version of f , u in , and u within Ω, respectively. The discrete counterpart of the convolution with the Green function in ( 1) is denoted by G :

C N → C N .
Finally, given u within Ω, G : C N → C M provides the scattered field y sc . We refer the reader to [START_REF] Pham | Three-dimensional optical diffraction tomography with Lippmann-Schwinger model[END_REF] for details concerning the implementation of G and G.

THREE-STEP RECONSTRUCTION

In the spirit of [START_REF] Soulez | A learn 2D, apply 3D method for 3D deconvolution microscopy[END_REF], we designed a three-step reconstruction scheme.

1. TV-regularized reconstruction.

2. Dictionary learning (DL) from lateral planes.

3. Final reconstruction using the learnt dictionary.

TV-Regularized Reconstruction

The first step consists in solving the nonnegative TV-regularized problem

f TV ∈ arg min f ∈R N ≥0 Q q=1 Hq(f ) -y sc q 2 2 + τTV f TV , (5)
where Hq : C N → C M denotes the two-step forward model described by (3) and (4) for the qth incident wave u in q , • TV ∇ • 2,1 is the TV seminorm, and τTV > 0 is the regularization parameter which balances between the data-fidelity term and the TV term. As in [4], we deploy a forward-backward splitting (FBS), implemented using the GlobalBioIm library 1 [START_REF] Soubies | Pocket guide to solve inverse problems with GlobalBioIm[END_REF], to obtain f TV .

DL from Lateral Planes

Given the nonnegative TV solution f TV , our goal is to learn a dictionary D ∈ R n×K formed out of K atoms of size n < N such that f TV can be represented as

f TV = P XY p=1 (R XY p ) T Dαp, (6)
where R XY p : R N → R n is an operator that extracts a 2D patch of size n centered on the pth element of the input vector (its adjoint inserts the patch at the pth position), PXY denotes the number of patches, and {αp ∈ R K } P XY p=1 are sparse vectors. The superscript XY in R XY denotes the fact that the operation only extracts 2D patches from lateral planes.

We formulate the DL problem as { D, αp} ∈ arg min

D∈R n×K αp∈R K P XY p=1 (R XY p ) T Dαp -f TV 2 2
+ τDL

P XY p=1 αp 1 , (7)
1 http://bigwww.epfl.ch/algorithms/globalbioim/ where τDL > 0 controls the sparsity level.

Our formulation is fundamentally different from the pioneering approaches [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF][START_REF] Mairal | Online dictionary learning for sparse coding[END_REF] where the solution was such that each extracted patch had a sparse representation in D. Our representation f = P XY p=1 (R XY p ) T Dαp is related to convolutional dictionary learning (CDL) [START_REF] Zeiler | Deconvolutional networks[END_REF][START_REF] Mørup | Transformation invariant sparse coding[END_REF][START_REF] Garcia-Cardona | Convolutional dictionary learning: A comparative review and new algorithms[END_REF], as shown by Papyan et al. [START_REF] Papyan | Convolutional dictionary learning via local processing[END_REF][START_REF] Papyan | Working locally thinking globally: Theoretical guarantees for convolutional sparse coding[END_REF]. As opposed to traditional DL, CDL accounts for global information in the image such as shift invariance. Hence, (8) not only enjoys the global sparse representation of CDL but also benefits from the local (patch-based) processing of DL [START_REF] Papyan | Convolutional dictionary learning via local processing[END_REF][START_REF] Papyan | Working locally thinking globally: Theoretical guarantees for convolutional sparse coding[END_REF].

By introducing the auxiliary variable sp = Dαp in (7), we can deploy the alternating-direction method of multipliers (ADMM) to minimize the augmented-Lagrangian functional

L( D, αp, sp, wp) = P XY p=1 (R XY p ) T sp -f TV 2 2 + P XY p=1 ρ 2 sp -Dαp + wp ρ 2 2 + τDL αp 1, (8)
where {wp ∈ R n } P XY p=1 are the dual variables and ρ > 0 is the Lagrangian multiplier. Using the CDL terminology [START_REF] Papyan | Working locally thinking globally: Theoretical guarantees for convolutional sparse coding[END_REF], the auxiliary variable sp is referred to as the pth slice. The ADMM is implemented using the SPAM toolbox2 [START_REF] Mairal | Sparse modeling for image and vision processing[END_REF].

Final Reconstruction Using the Learnt Dictionary

Equipped with the dictionary D ∈ R n×K learnt from lateral planes in Section 3.2, we now consider the optimization problem

               f = R T s * , s * ∈ arg min s∈R nP Q q=1 Hq(R T s) -y sc q 2 2 +Rpos(R T s) + R D (s) , (9)
where s = [s T 1 • • • s T P ] T ∈ R nP is the concatenation of all the slices, P = PXY + PXZ + PYZ is the total number of slices, and

R = [R T XY R T XZ R T YZ ] T ∈ R nP ×N
with RXZ (RYZ) the counterpart of RXY for the XZ (YZ, respectively) sections of the volume. We use the differentiable functional Rpos : R N → R

Rpos(f ) = N n=1 λpos µ log (exp(-µfn) + 1) (10)
to favor nonnegative solutions. Here, λpos > 0 is a weighting factor and µ > 0 shapes the tolerance to negative values. Finally, the functional R : R nP → R in ( 9) is a regularization term designed to enforce the slices sp to have a sparse representation in D. Denoting

F = Q q=1 Hq(R T • ) -y sc q 2 2 + Rpos(R T • )
, we can deploy the FBS algorithm whose iterates are given by

s k+1 = prox γR D s k -γ∇F(s k ) , (11)
where γ > 0 is a descent parameter and prox γR D denotes the proximity operator of the functional R D . Here, we follow the plug-andplay prior philosophy [START_REF] Venkatakrishnan | Plug-and-play priors for model based reconstruction[END_REF][START_REF] Sreehari | Plugand-play priors for bright field electron tomography and sparse interpolation[END_REF][START_REF] Kamilov | A plug-andplay priors approach for solving nonlinear imaging inverse problems[END_REF] and replace prox γR D in (11) by the Algorithm 1 Proposed algorithm to solve ( 9)

Require: y ∈ C M Q , x 0 ∈ C N , D ∈ R n×K , γ > 0,τsc > 0 1: Define F = Q q=1 Hq(R T • ) -y sc q 2 2 + Rpos(R T • ) 2: s 0 = 1 n 2 Rx 0 3: k = 1 4: while (not converged) do 5: z k = s k -γ∇F(s k
where τsc > 0, and

α ∈ arg min α∈R K 1 2 s -Dα 2 2 + τsc α 1 . (13)
Numerous solvers exist to solve (13) [START_REF] Mairal | Online dictionary learning for sparse coding[END_REF][START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. Again, we use the GlobalBioIm library together with the SPAMS toolbox for this step. Finally, we summarize the complete reconstruction scheme of this section in Algorithm 1.

RESULTS

We validated our method on real data acquired with the ODT setup described in [START_REF] Ayoub | A method for assessing the fidelity of optical diffraction tomography reconstruction methods using structured illumination[END_REF]. We acquired 61 views with illumination of wavelength λ = 532nm within a cone of 35

•

. Each view has (300 × 300) (complex) measurements focused at the center of the sample. The sample is a yeast cell (Figure 1) immersed in water (n b = 1.3388). The size of the reconstructed volume is (96 × 96 × 96) with a sampling step of 99.3nm (9.53µm in each dimension).

TV-Regularized Reconstruction

We first reconstructed the sample using the method described in [4] by minimizing [START_REF] Pham | Three-dimensional optical diffraction tomography with Lippmann-Schwinger model[END_REF]. The initial guess was the solution provided by the Fig. 3. Dictionary D learnt from lateral planes of the TV regularized solution (f TV , dotted rectangle in Figure 2).

Rytov model [START_REF] Devaney | Inverse-scattering theory within the Rytov approximation[END_REF]. We used diverse regularization parameters τTV for TV (Figure 2). When the regularization is weak, artifacts due to model mismatch are hindering the quality of reconstruction. On the contrary, over-regularization results in cartoon-like solutions.

Learning the Dictionary

We learned the dictionary by minimizing [START_REF] Goodman | Speckle Phenomena in Optics: Theory and Applications[END_REF]. We used patches of size (8 × 8) (n = 64) and K = 64 atoms. We set τDL = 1/ √ 8 and ρ = 0.5 max(f TV )n 2 . The learnt atoms of the dictionary D are shown in Figure 3.

Final Reconstruction

We solved the optimization problem ( 9) and encouraged the nonnegativity of the solution with λpos = 1/96 3 and µ = 5 in (10). The initial guess was f TV . The denoising operator ( 12)-( 13) was used with regularization parameter τsc = 10 -4 .

We observe that the Rytov-based solution suffers from the missing-cone problem whereas the regularized solutions (i.e., TV and the proposed one) mitigate its effect. In addition, some features are enhanced with the proposed solution in comparison to the TV solution (Figure 4). Finally, the proposed method is able to recover features in deeper axial position whereas the TV-regularized solution is over-regularized (Figure 4, right column).

CONCLUSION

We designed an adaptive regularization that allowed us to improve the quality of optical diffraction tomography (ODT) reconstructions in the axial direction using features learned in lateral planes. The proposed regularization relies on a dictionary that is learnt from the lateral planes of an initial total-variation reconstruction. This dictionary is then used in a final step to enhance the quality of the reconstruction in all XY, XZ, and YZ sections. We applied this strategy to the reconstruction of real ODT measurements. Our results show the superior performance of the proposed pipeline over conventional regularizations.

Fig. 1 .

 1 Fig. 1. Principle of optical diffraction tomography. The wave vectors {k inq } Q q=1 ∈ R 3 of the Q incident plane waves {u in q } Q q=1 are limited to a cone around the optical axis.

  )

=Fig. 2 .

 2 Fig. 2. TV-regularized solutions obtained with regularization parameters τTV = {0.15, 0.3, 0.5}/96 3 (from left to right, respectively).

http://spams-devel.gforge.inria.fr/


Our objectives

Main goal: design, implement and test different placement and scheduling policies at both QNode-and QBox-level.

Main problem: testing on a

From Heater to Computing Resources

Main problem: testing on a production platform is not conceivable and takes time.

⇒ Simulation is what you need! SimGrid1 : Large-scale distributed system simulator with execution and communication models. → Used to simulate platform and tasks execution.

Batsim2 : Infrascructure simulator for jobs and I/O scheduling. → Used to drive the simulation, submit tasks and communicate with the decision process.

Pybatsim3 : Batsim's Python API exposing methods to easily communicate with the Batsim process. → Used to implement the QNode and QBox schedulers.

Simulation Overview

"

  A disruptive solution to turn IT waste heat into a viable heating solution for buildings." The Qarnot platform: ∼1,000 distributed QRads embedding ∼3,000 diskless computing units (QMobos) ∼20 local servers (QBoxes) with disks 1 global server (QNode) with a centralized storage server (CEPH) disappear over time: the inhabitants decide! Available resources when heating is required (QRad is On) Unavailable resources when ambient air is too warm (QRad is Off) Make global decisions at QNode-level: Decide where to dispatch (groups of) instances Ensure global load-balancing Make local decisions at QBox-level: Schedule instances on QRads Regulate room temperature (via DVFS) Ensure heating needs are satisfied at QNode-level: Decide where to dispatch (groups of) instances Ensure global load-balancing Make local decisions at QBox-level: Schedule instances on QRads Regulate room temperature (via DVFS) Ensure heating needs are satisfied Need to reach 100% of platform usage.⇒ An idle QRad is a lack of heating design, implement and test different placement and scheduling policies at both QNode-and QBox-level.

June 27th, 2019

https://github.com/simgrid/simgrid

https://gitlab.inria.fr/batsim/batsim/tree/temperature

https://gitlab.inria.fr/batsim/pybatsim/tree/temperature Clément Mommessin Scheduling at the Edge June 27th, 2019

Computing Requests

Multiple Objectives

Temperature Regulation

Frequency and temperature regulator in each QRad:

In general: DVFS on QMobos to adapt power consumption When too hot: instances are killed and re-submitted to the QNode When too cold (lack of heating): "background" compute-intensive instances are generated (best-effort blockchain mining )

Developed Extensions

Temperature model: to compute the temperature of the QRad and ambient air, based on thermodynamics formulae.

External events injector: to replay a machine failure or a temperature change.

Storage controller: to manage storage entities and data movements.

Clément Mommessin Scheduling at the Edge June 27th, 2019

Variants of the QNode dispatcher:

FullReplicate (instantaneous transfers)

Standard Qarnot scheduler at QBox-level.

1-week simulation inputs from real logs of Qarnot.

Clément Mommessin Scheduling at the Edge June 27th, 2019

List of n instances from Cloud or IoT, with for each instance j:

An estimation of work

Power diagram

Clément Mommessin Scheduling at the Edge June 27th, 2019

THE questions to answer:

Where to execute an instance?

When to execute an instance?

At which speed?

⇒ Similar to power capping problems Objectives:

Minimize ∑ j w j C j (where C j is the completion time)

Minimize power consumption distance to the power diagrams