Jordan Paul

David

Estelle Lucas

Jean-Baptiste Sami

Anis, Athmane, François Boris N Noëlig

Richard F Yannick

Karol

Houssam... William David

Ainsi

Keywords: Analyse Numérique, Dynamique Moléculaire, Équations Différentielles Stochastique, Modèles multi-échelles, Modèles de réduction de complexité Numerical Analysis, Molecular Dynamics, Stochastic Differential Equations, Multiscale models, Coarse-Graining models. vii

Développement de schémas numériques d'intégration de méthodes multi-échelles Résumé: Cette thèse concerne l'analyse et le développement de schémas d'intégration numérique de la Dynamique des Particules Dissipatives. Une présentation et une analyse de convergence faible de schémas existants est présentée, suivie d'une présentation et d'une analyse similaire de deux nouveaux schémas d'intégration facilement parallélisables. Une analyse des propriétés de conservation d'énergie de tous ces schémas est effectuée suivie d'une étude comparative de leurs biais sur l'estimation des valeurs moyennes d'observables physiques pour des systèmes à l'équilibre. Les schémas sont ensuite testés sur des systèmes choqués de fluides DPDE, où l'on montre que nos deux nouveaux schémas apportent une amélioration dans la précision de la description du comportement de tels systèmes par rapport aux schémas facilement parallélisables existants. Finalement, nous présentons une tentative d'accélération d'un schéma d'intégration de référence s'appliquant aux simulations séquentielles de la DPDE.

physique sadiques que je leur ai imposé lors de ces séances Ô combien libératrices. Merci donc à Martial, Jean-Luc, Felix, Christophe, Johann, Sophie, Marilyn et Paul-Émile. Ainsi que ceux qui n'ont été que de passage ou que je n'ai vu que temporairement. Vous aurez l'occasion de vous venger, ne vous en faites pas.

Finalement, je remercie tous mes ami(e)s extérieurs au Cea ou au Cermics avec qui j'ai pu partager les joies et peines de ce doctorat: Pavel Borisevich, Bruno Premoselli, Remy Rodiac, Romeo Hatchi, Guilhem Blanchard, Florian Vienot, Benoit Costanzo, Stanislas Morbieu, Simon Bittman, Romain Prudhomme, Camille Poret, Julie Verbeke, Bruno Guillon, Vincent Belz, Chloé Leclère, Camille Mazo et bien sûr tous ceux que j'oublie. Dédicace spéciale à Pavel Borisevich et Bruno Premoselli, avec qui je me suis tellement étalé sur ma vie de doctorant qu'ils pourraient (presque) parler de mon doctorat aussi bien que moi, ou du moins écrire un livre sur son déroulement (pour ceux que ça intéresseraient...).

Cette page de remerciements ne serait pas complète si je ne mentionnais pas ma famille, et en particulier mes parents Mokhtar et Aïcha et ma soeur Sawsan, qui eux aussi ont su m'aider et m'apporter leurs conseils, certes sur des sujets moins techniques, mais tout aussi important. Ils ont vécu avec moi cette thèse, et ont (trop?) souvent dû supporter mes sauts d'humeurs et autres sentiments qui m'ont parfois accompagné lors de ce doctorat et que je ne pouvais pas montrer à d'autres personnes qu'eux. Ils l'ont fait avec une patience, un dévouement et un amour que j'aurai du mal à décrire avec des mots. Je leur dédie cette thèse, et espère un jour leur rendre ne serait-ce qu'une partie de ce qu'ils m'ont apporté.

Et enfin, merci à toi Augustina, dont le réconfort et l'amour que tu m'as apporté m'ont aidé à traverser cette dernière année de thèse, qui a été la plus difficile. Quelle chance j'ai de t'avoir eu à mes côtés durant cette épreuve! Je remercie donc très, très, très chaleureusement toutes les personnes que j'ai cité dans ces remerciements, qui sont bien trop courts en comparaison du soutien et de l'aide qu'elles m'ont donné.

Bonne lecture à tous! Invariant measure of the considered numerical scheme, unless otherwise mentioned D Set of infinitely differentiable functions from X to R, whose derivatives grow at most polynomially D 0 Set of elements of D of null average, i.e functions ϕ ∈ D with X ϕdπ = 0 D ∞ Set of infinitely differentiable functions from X to R with compact support Π Projection of functions of L 1 (π) from X to R to functions of null average, i.e Πϕ = ϕ -X ϕdπ Π ∆t Equivalent of Π for π ∆t , i.e Π ∆t ϕ = ϕ -X ϕdπ ∆t

Résumé substantiel

Cette thèse est consacrée au développement et à l'étude de nouveaux schémas d'intégration numérique d'un modèle multi-échelle nommé Dynamique des Particules Dissipatives avec conservation de l'Énergie (DPDE).

Un des objectifs du Cea/Dam est l'étude et la compréhension de systèmes choqués de matériaux complexes. Pour cela, il est nécessaire d'effectuer des simulations numériques d'ondes de choc pour de tels systèmes. Cependant, ces systèmes d'intérêt impliquent des phénomènes dont les échelles de temps et de distance sont très différentes, nécessitant ainsi des modèles multi-échelles afin de pouvoir effectuer des simulations numériques à des coûts raisonnables en termes de nombre d'opérations flottantes et de temps de simulation.

Un des modèles choisis est la Dynamique des Particules Dissipatives, dénotée DPD, dans lequel une ou plusieurs molécules sont représentées par une seule mésoparticule, interagissant avec les autres mésoparticules non seulement par des forces conservatives issues d'un potentiel effectif mais aussi par des interactions dissipatives et aléatoires, nommées interactions de fluctuation/dissipation. Ce modèle a été introduit par Hoogerbrugge et Koelman en 1992 [START_REF] Hoogerbrugge | Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics[END_REF], et ses fondations théoriques ont été établies en 1995 par Español et Warren [START_REF] Español | Statistical Mechanics of Dissipative Particle Dynamics[END_REF]. Cependant, la DPD est un modèle préservant la température du système mais pas son énergie, et n'est donc pas adaptée à la simulation de systèmes hors-équilibre. Afin de pouvoir simuler des systèmes hors-équilibre, il est donc nécessaire d'utiliser la variante conservant l'énergie de la DPD, dénotée DPDE. Cette dynamique est construite en ajoutant un degré de liberté aux particules DPD traduisant l'évolution de tous les degrés de libertés internes de la molécule représentée par une particule DPDE. L'énergie cinétique dissipée par les interactions de fluctuation/dissipation de la DPD est ensuite réinjectée dans l'ensemble de ces degrés de liberté s internes. La DPDE préserve donc par construction une énergie E définie comme la somme de ses trois composantes liées aux trois types de degrés de liberté du système: l'énergie potentielle, l'énergie cinétique et l'énergie interne, définie comme la somme des énergies internes des particules. Cette xv énergie s'écrit donc

E(q, p, ε) = U(q) + N i=1 p 2 i m i + N i=1 ε i ,
où m i représente la masse de la particule indexée par i et U représente le potentiel conservatif considéré. La Figure 1 illustre le procédé de modélisation de la DPDE en donnant un sens physique aux différentes variables.

Figure 1 | Illustration du procédé de modélisation de la DPDE: chacune des trois molécules (en violet) interagissant entre elles par des forces conservatives (lignes droites: continues pour les interactions à courte portée et en pointillées pour les interactions à longue portée) est représentée par une seule macroparticule DPDE (cercles en pointillés). Chacune de ces macroparticules interagit avec les autres via des interactions conservatives, dissipatives et aléatoires: les interactions entre particules DPDE ne sont donc plus deterministes mais stochastiques. En plus de sa position et de son moment, chaque macroparticule DPDE est décrite par une variable supplémentaire qui est l'énergie interne (plus la particule est rouge, plus cette variable est grande) représentant les degrés de libertés internes de sa molécule associée. L'équilibration des degrés de libertés internes/externes est faite par des transferts d'énergie liés aux interactions de fluctuation/dissipation (flèches rouges).

Considérons un système de N particules, de positions q ∈ R dN , de moments p ∈ R dN et d'énergies internes ε ∈ R N , où d représente la dimension de chaque particule (usuellement, d = 3). La dynamique DPDE pour un tel système s'écrit de la façon suivante :

                               dq i,t = p i,t m i dt,
dp i,t = -∇ q i U(q t ) dt + N j=1,j i -γ ij,t χ(r ij,t )v ij,t dt + σ ij,t χ(r ij,t )dW ij,t ,

dε i,t = 1 2         N j=1,j i        γ ij,t v 2 ij,t -d σ 2 ij,t 2 µ ij        χ(r ij,t ) dt -σ ij,t χ(r ij,t )v ij,t • dW ij,t         , avec µ ij = 1 m i + 1 m j .
Dans toutes les applications de cette thèse, nous allons considérer des potentiels de type Lennard-Jones s'écrivant

U(q) = 1≤i<j≤N v(r ij ), r ij = q i -q j , v(r) = 4ε LJ σ LJ r 12 - σ LJ r 6 ,
où ε LJ et σ LJ sont des paramètres du potentiel. La fonction χ représente la portée des interactions de fluctuation/dissipation. Dans cette thèse, nous choisissons

χ(r) =            1 - r r cut 2
, si r ≤ r cut , 0, sinon.

Le terme v ij,t = v i,t -v j,t désigne la vitesse relative entre les particules indexées par i et j, et les termes (W ij,t ) 1≤i<j≤N sont des mouvements Browniens vérifiant la condition d'antisymétrie W ij,t = -W ji,t . Les paramètres γ ij et σ ij pour 1 ≤ i j ≤ N représentent respectivement l'intensité des interactions de fluctuation/dissipation entre les particules i et j.

Nous pouvons définir une température interne associée à chaque énergie interne. Considérons une particule DPDE d'énergie interne ε ∈ R + . En supposant que les temps et distances caractéristiques des degrés de libertés internes de chaque molécule représentée par une particule DPDE sont bien plus petits que ceux des degrés de libertés externes, nous pouvons supposer que les degrés de libertés sont constamment à l'équilibre, et nous pouvons donc définir l'équation d'état interne de chaque particule DPDE par

ε = T (ε) 0 C v (θ) dθ. (1)
Ici, T (ε) représente la température interne de la particule et C v sa capacité calorifique à volume constant. Nous voyons que si C v est continue et strictement positive sur R + , (1) est bien posée. En réalité, C v est croissante sur R + , ce qui, combiné à ce qui précède, implique lim ε→0 + T (ε) = 0, et lim ε→+∞ T (ε) = +∞.

Une fois définie la température interne de chaque particule, nous pouvons définir l'entropie s(ε) associée aux degrés de libertés interne par

s ′ (ε) = 1 T (ε)
, où s ′ (ε) représente la dérivée de la fonction s.

Nous remarquons que la DPDE préserve deux invariants que sont le moment total et l'énergie totale:

d        N i=1 p i        = 0, et d        U(q) + 1 2 N i=1 p i m i + N i=1 ε i        = 0.
L'ensemble thermodynamique associé à la DPDE doit donc prendre en compte ces deux invariants. Nous pouvons construire de façon similaire à l'ensemble microcanonique associé à la dynamique Hamiltonienne (voir [START_REF] Lelièvre | Free Energy Computation : A Mathematical Perspective[END_REF] pour le détail de la construction de cet ensemble), l'ensemble dit "microcanonique" de la DPDE, noté µ E 0 ,P 0 . Cet ensemble statistique s'écrit µ E 0 ,P 0 (dq, dp, dε) = Z -1 E 0 ,P 0 δ E 0 (q, p, ε)δ P 0 (p)e S (ε)

k B dqdpdε,
où Z E 0 ,P 0 est une constante de normalisation et S (ε) est l'entropie interne du système, définie comme la somme des entropies internes des particules

S (ε) = N i=1 s i (ε i ).
Cependant, afin de proposer des estimateurs d'observables physiques telles que la température, il est plus simple d'utiliser l'ensemble dit "canonique" de la DPDE s'écrivant β,C v denotent respectivement les ensembles "microcanonique" et "canonique" de la DPDE associés à un système de N particules.

Le seul résultat d'ergodicité de la DPD a été prouvé dans un cas réduit unidimensionnel avec une densité assez élevée [START_REF] Shardlow | Geometric Ergodicity for Dissipative Particle Dynamics[END_REF]. En règle générale, aucun résultat d'ergodicité n'a été prouvé pour la DPDE. Cependant, il peut être montré que, sous réserve que γ ij et σ ij satisfassent

σ ij = σ, γ ij = σ 2 4k B 1 T i (ε i ) + 1 T j (ε j ) ,
les mesures de probabilité µ E 0 ,P 0 et µ β,C v sont invariantes pour la DPDE. Dans la suite, nous supposons que la DPDE est ergodique pour µ E 0 ,P 0 .

De plus amples détails sur la Physique Statistique et ce que sont les Équations Différentielles Stochastiques ergodiques sont donnés en Section 1.2, et de plus amples détails sur la DPDE et sa thermodynamique peuvent êtres trouvés en Section 2.1.

Nous considérons dans cette thèse trois schémas numériques d'intégration de la DPDE existants décrits et analysés en Section 2.2, et notés SVV, SEM et SSA. Le schéma SVV, pour Stochastic Velocity-Verlet, est l'extension à la DPDE du schéma Velocity-Verlet pour la dynamique Hamiltonienne. Les schémas SEM, pour Splitting Euler-Maruyama et SSA, pour Shardlow's Splitting Algorithm sont eux basés sur une technique dite de splitting séparant l'intégration de la partie conservative et de la partie de fluctuation/dissipation de la DPDE.

xviii La partie conservative est donnée par la dynamique Hamiltonienne et s'écrit

           dq i,t = M -1 p i,t ,
dp i,t = -∇ q i U(q i,t ), dε i,t = 0, avec M étant la matrice de masse du système (souvent, M = Diag (m 1 , .., m N )). Cette dynamique est toujours intégrée avec un schéma de Velocity-Verlet [START_REF] Verlet | Computer "experiments" on classical fluids. i. thermodynamical properties of lennard-jones molecules[END_REF]. La partie de fluctuation/dissipation s'écrit

                           dq i,t = 0, dp i,t = - N j=1,j i γ ij,t χ(r ij,t )v ij,t dt + N j=1,j i σ χ(r ij,t ) dW ij,t , dε i,t = 1 2         N j=1,j i       γ ij,t v 2 ij,t -d µ ij σ 2 2       χ(r ij,t ) dt -σ χ(r ij,t )v ij,t • dW ij,t         .
(

) 2
La dynamique (2) est discretisée par SEM avec un schéma basique d'Euler-Maruyama, alors que SSA utilise encore une fois une technique de séparation en la décomposant en une somme de dynamiques élémentaires associées à chaque couple 1 ≤ i j ≤ N de particules s'écrivant

                           dp i,t = -γ ij,t v ij,t χ(r ij,t ) dt + σ χ(r ij,t ) dW ij,t ,
dp j,t = -dp i,t ,

dε i,t = - 1 2 d        p 2 i,t 2m i + p 2 j,t 2m j        ,
dε j,t = dε i,t .

L'intégration SSA de chaque dynamique élémentaire est effectuée en deux étapes: les moments sont intégrés avec un schéma d'ordre faible infini et les énergies internes sont mises à jour en réinjectant l'énergie cinétique dissipée par l'intégration des moments. De plus amples détails sur cette intégration sont donnés en Section 2.2.3.

Nos deux nouveaux schémas, nommés SER pour Splitting with Energy Reinjection et Hybrid, décrits et analysés en Section 2.3, utilisent aussi cette séparation de la DPDE en une dynamique conservative et une dynamique de fluctuation/dissipation, et intègrent aussi la partie conservative avec un schéma de Velocity-Verlet. La discrétisation de (2) par SER est faite en intégrant les moments avec un Euler-Maruyama

p n+1 i = p n i + δp n i , δp n i = N j=1,j i δp n ij , xix avec n ∈ N et δp n ij = -γ n ij χ(r n ij )v n ij ∆t + σ χ(r n ij )G n ij √
∆t, puis en remarquant que

∆K n i = (p n+1 i ) 2 2m i - (p n i ) 2 2m i , = δp n i 2m i 2p n i + δp n i , = N j=1,j i δp n ij • v n i + δp n i 2m i ∆ j K n i ,
La mise à jour des énergies internes se fait ensuite en réinjectant de façon symétrique les variations d'énergie cinétique de chaque particule dans les énergies internes. Ceci revient à écrire

ε n+1 i = ε n i - 1 2 N j=1,j i ∆ i K n j + ∆ j K n i , = ε n i - 1 2 N j=1,j i δp n ij •       v n ij + 1 2       δp n i m i - δp n j m j             .
Le schéma Hybrid quant à lui est, comme son nom l'indique, un mélange entre SSA et SER, construit de façon à se débarasser des problèmes de parallélisation du SSA en remplacant la discrétisation de certaines interactions de fluctuation/dissipation élémentaires par une discrétisation globale SER n'impliquant que ces interactions en question (voir la Section 2.3.2 pour plus de détails).

Ces cinq schémas sont ensuite comparés dans la Section 3.1 avec des simulations de systèmes DPDE à l'équilibre. Nous étudions tout d'abord leurs propriétés de conservation d'énergie en calculant l'énergie moyenne dépendant du temps sur N traj ≫ 1 trajectoires, i.e

E t = E µ β,C v [E(x t )] .
où x t = (q t , p t , ε t ) est une solution de la DPDE. En pratique, E t est estimée par

E n ∆t,N traj = 1 N traj N traj m=1 E (x m,n ) .
où x m,n est la n-ième configuration de la m-ième trajectoire calculée par le schéma numérique considéré. Le graphe de gauche de la Figure 2 montre que les énergies moyennes pour SVV, SSA et SER dérivent linéairement en fonction du temps. Ceci nous permet de calculer les coefficients de dérive en fonction du pas de temps utilisé. Le graphe de droite de la Figure 2 montre que ces coefficients augmentent polynomialement avec le pas de temps. Nous observons donc qu'aucun des schémas ne conserve correctement l'énergie du système, en notant cependant que l'on observe avec SER et SSA des dérives d'un ordre de grandeur inférieur à celles obtenues avec SVV. Des simulations supplémentaires montrent que des dérives linéaires s'observent aussi pour SEM, avec des coefficients de dérive d'un ordre de grandeur supérieur à ceux de SER et de SSA. Ces dérives ne permettant pas d'obtenir un état stationnaire bien défini, nous corrigeons tous les schémas par une renormalisation des énergies internes à l'issue de la discrétisation complète de l'équation, comme proposé dans [START_REF] Lisal | Dissipative Particle Dynamics at Isothermal, Isobaric, Isoenergetic, and Isoenthalpic Conditions using Shardlow-like Splitting Algorithms[END_REF]. Notant ( qn , pn , εn ) la configuration obtenue par une application complète (i.e après les discrétisations des dynamiques conservatives et de fluctuation/dissipation) de la procédure de discrétisation d'un schéma donné sur la configuration de départ (q n , p n , ε n ), ceci consiste à définir (q n+1 , p n+1 , ε n+1 ) = qn , pn , εn -

1 N (E( xn ) -E 0 ) 1 N .
Nous calculons ensuite en Section 3.2 les biais des schémas ainsi corrigés sur les estimations de valeurs moyennes d'observables pour des systèmes à l'équilibre. L'hypothèse d'ergodicité de la DPDE par rapport à µ E 0 ,P 0 nous permet, pour une observable physique ϕ donnée, d'estimer sa valeur moyenne E µ E 0 ,P 0 [ϕ] par

ϕ N it = 1 N it N it n=1 ϕ(x n ),
où (x n ) n=0,..,N it est une solution approchée de la DPDE sur 0, ∆t, .., N it ∆t calculée par le schéma numérique considéré et N it est le nombre d'itérations de la simulation. La Figure 3 présente les estimations des valeurs moyennes des composantes de l'énergie totale pour les schémas SSA, SEM, SER et Hybrid corrigés. Nous observons que les biais des schémas SER et Hybrid sont respectivement équivalents à ceux des schémas SEM et SSA, et que les biais de SSA/Hybrid sont plus petits d'au moins un ordre de grandeur que ceux de SER/SEM. Nous remarquons aussi que tous les schémas sont cohérents entre eux, et convergent vers une valeur commune propre à chaque observable lorsque le pas de temps ∆t tends vers 0.

Pour finir, les schémas SER, Hybrid et SEM sont comparés au Chapitre 4 sur des simulations d'ondes de choc de fluides DPDE. Les fluides sont tout d'abord équilibrés à T = 1000K pendant t = 5 × 10 -12 s, puis la matière est mise en mouvement à une vitesse v init = -2000m s -1 vers un piston composé de particules Lennard-Jones de masse infinie se trouvant à la gauche de la boîte de simulation. La Figure 4 illustre le procédé utilisé ainsi que les profils attendus d'observables moyennées sur des tranches normales à la direction de propagation du choc. De plus amples informations concernant la description des simulations peuvent se trouver en Section 4.1.

En utilisant un pas de temps ∆t = 10 -15 s, la simulation SEM présente des phases choquées et à l'équilibre qui chauffent au fur et à mesure de la simulation, comme le montrent les graphes du haut de la Figure 5 avec les profils de températures cinétiques et interne (voir la Section 2.1.2 pour plus de détails sur les différents estimateurs de la température en DPDE). Ces phénomènes sont accrus lorsque l'on augmente le pas de temps, ce qui montre que ce sont des artefacts numériques liès à l'utilisation d'un pas de temps fini. Ces artefacts ne se retrouvent pas dans les simulations effectuées avec SER et Hybrid, comme le montrent les profils du milieu et du bas de la Figure 5. Ces résultats montrent donc que nos deux nouveaux schémas apportent une amélioration dans la qualité de la description du comportement de matériaux sous choc par rapport aux schémas parallélisables existants. Figure 4 | Haut : schéma de la procédure utilisée: après équilibration, le système est mis en mouvement vers un piston gauche immobile, créant ainsi une onde de choc partant du piston et se propageant de la gauche vers la droite du système. Bas : profils de chocs attendus, calculés en moyennant les observables sur des tranches normales à la direction de propagation du choc. Deux types de profils sont attendus: les observables à temps d'équilibration "rapides" présenteront un pic au front de choc puis s'équilibreront vers la valeur d'équilibre de choc (courbe rouge), alors que les observables à temps d'équilibration "lents" s'équilibreront directement vers la valeur d'équilibre de choc sans présenter de pic (courbe bleue). L'onde de relaxation à droite de la figure est due à l'absence de piston à droite de la boîte de simulation, et n'affecte pas le front de choc tant que les deux ondes ne se rencontrent pas.

xxii t = 0 m = ∞ x = 0 v av ≃ v init t > 0 m = ∞ x = 0 v av ≃ 0 v choc v av ≃ v init x ϕ 0 équilibration lente équilibration rapide v choc v relax ϕ choc ϕ equ

Context

Molecular Dynamics (MD) is the field of science dedicated to the study of matter at its microscopic level, using computers to model the interactions between atoms or molecules. It has many aims, and we consider in this work two of them that are, to the extent of our knowledge, the most important:

• The first interest consists in using MD as a numerical microscope to observe the behavior of matter at microscopic level, thus performing the computer equivalent of what is traditionally understood as experiments. Physical theories can thus be tested in "computer experiments", and even help their elaboration by guiding the physicist's intuition with the observation of molecular trajectories. It was, for example, the motivation of the first molecular simulation performed by Alder and al. in 1958 [AW58].

• The second major interest of MD consists in computing macroscopic quantities according to the theory of Statistical Physics, often obtained as averages of microscopic quantities. In this case, MD can be used for a quantitative description of matter, allowing to compute phase diagrams for instance.

Molecular Dynamics often considers a system of particles whose positions and momenta are solutions of a system of differential equations. The behavior of the particles is thus obtained by computing the solutions of this differential system or dynamical system. However, in most cases, analytical solutions of these dynamical systems are not known, and we must compute approximate solutions with computers. Those solutions are computed over a finite number of times (t 1 , . . . , t N ) by a discretization of the differential equations, or numerical scheme.

Computing these approximate solutions at each time t 1 , .., t N requires many floating point operations and can be very computationally expensive. The most accurate descriptions in Molecular Dynamics, called Ab-Initio Molecular Dynamics, take into account the quantum effects in order to describe the interaction forces between the particles. However, they require a tremendous amount of floating point operations, and today computational capacities do not allow us to simulate systems composed of several hundred atoms. Less accurate but less expensive descriptions use empirical parameters to describe the particle interactions. Nevertheless, even using such empirical descriptions does not allow to simulate systems of more than several hundred million particles on current supercomputers.

One might think that a simulation involving a hundred million atoms is enough to compute whatever we want to compute, but it is unfortunately not the case. Let us illustrate this by a simple example. A mole consists of around 6 × 10 23 atoms, and a glass of water contains 5 moles of water molecules. Therefore, simulating a glass of water would require to take into account approximately 5 × 10 16 times the maximal number of atoms that we can simulate nowadays. This simple example reveals that even very small systems in our scale of length and time are much too large for our computational power to simulate using traditional Molecular Dynamics. In practice, it is difficult to describe systems of length scales greater than several hundreds of nanometers because of the computational power such simulations require. This upper limit is reduced for some complex systems. The equation of Hydrodynamics can describe most fluids of length scale greater than a micrometer: below it, the granular aspect of matter is often necessary to understand their physics. We therefore see that there is a gap in the time and length scales between systems described by Molecular Dynamics and those described by the equations of Hydrodynamics, where some hybrid models are required (see Figure 6 for a graph of the length and time scales of different descriptions of matter). In addition, some systems involve phenomena of different orders of magnitude in terms of length and time scale. A striking example is the one of polymers in a solvent. Polymers are very long molecules with a filamentary look, and the solvent is usually composed of much smaller and simpler molecules. In these systems, the characteristic times for the polymers are one or two orders of magnitude longer than those of the solvent molecules. Hence, using traditional MD would require to use a time step corresponding to the smallest characteristic time, thus implying prohibitively long and costly simulations in order to observe any characteristic movement of the polymers. One solution to model such phenomena in reasonable times and costs is to use coarse-graining models, representing the smaller and shorter phenomena in an average manner [START_REF] Mori | Transport, collective motion and brownian motion[END_REF][START_REF] Zwanzig | Nonlinear generalized langevin equations[END_REF], thus speeding up the computations by allowing the use of time steps corresponding to the larger phenomena.

The Dissipative Particle Dynamics, termed DPD, is a coarse-graining model developped by Hoogerbrugge and Koelman in 1992 [START_REF] Hoogerbrugge | Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics[END_REF], and put in firm theoretical grounds by Español and Warren in 1995 [START_REF] Español | Statistical Mechanics of Dissipative Particle Dynamics[END_REF]. It represents each molecule, or group of molecules, by one single mesoscale particle interacting with the other mesoscale particles by conservative forces deriving from the potential function, but also by friction and random interactions. Dissipative Particle Dynamics has been intensely used for its computational simplicity compared to classical Molecular Dynamics and many efficient and accurate numerical discretizations of DPD exist. However, it suffers from two downsides: the first is its incapacity to control the amplitude of the coarse-graining process (i.e the number of molecules represented by the DPD particle for instance) in another way than acting on the conservative potential, and the second is its incapacity of simulating phenomena outside equilibrium, being a model preserving the temperature and not the energy.

The Cea/Dam is interested in simulating the behavior of shocked particle systems. These systems are often very large, and require the use of massively parallel simulations. Shockwave simulations of a hundred million copper and tin atoms have already been performed using Molecular Dynamics [START_REF] Durand | Large-Scale Molecular Dynamic Study of Jet Break-up and Ejecta Production from Shock-Loaded Copper with a Hybrid Method[END_REF], and codes designed for future architectures of supercomputers allow to reach similar systems of several hundreds of thousands of atoms [START_REF] Cieren | Molecular Dynamics for Exascale Supercomputers[END_REF]. However, other systems of interest are more complex to simulate and such length scales cannot be reached. Therefore, coarse-graining models adapted to nonequilibrium simulations are desirable.

The Dissipative Particle Dynamics with Energy conservation, termed DPDE, is an energy-conserving variant of the Dissipative Particle Dynamics introduced independently by Avalos and al. and Español in 1997 [START_REF] Avalos | Dissipative Particle Dynamics with Energy Conservation[END_REF][START_REF] Español | Dissipative Particle Dynamics with Energy Conservation[END_REF]. Its energy preservation allows the modelling of systems outside equilibrium, and is therefore an interesting option for large scale numerical simulations of shocked complex systems. In addition to preserving (as its name indicates) the energy of the system, DPDE takes into account the heat capacity of the coarse-grained molecules, thus allowing to control in some way the amplitude of the coarse-graining process in opposition to the isothermal DPD. Where DPD enjoys many efficient numerical discretizations, DPDE turns out to be harder to integrate: in addition to being accurate, DPDE discretizations must preserve the invariant of the dynamics, namely the total momentum and energy. While the preservation of the momentum is often trivially ensured, preserving the energy is sometimes more difficult.

As said before, in the Cea/Dam objectives, DPDE is destined to be used for large scale and massively parallel simulations of shocked systems. Therefore, DPDE numerical schemes need to be parallelizable in addition to being accurate and energy-preserving. Up to the extent of our knowledge, no scheme manages to satisfy both conditions so far, the most promising option being a scheme introduced by Shardlow in 2003 [START_REF] Shardlow | Splitting for Dissipative Particle Dynamics[END_REF], very efficient but complex to parallelize [LBM + 14].

The object of this thesis is therefore to construct straightforwardly parallelizable numerical tools to accurately integrate the DPDE model. We organize this thesis in the following manner:

• Chapter 1 is an introductive chapter presenting all the notions needed for the following chapters. Information about the theory of Molecular Dynamics, numerical schemes to integrate classical particle systems and the theory allowing to compute macroscopic quantities from microscopic information are provided. The Dissipative Particle Dynamics being a Stochastic Differential Equation, notions about SDEs, their integration and their use for computing macroscopic quantities are then given. Finally, all of this is applied on a prototypical example similar to DPD called the Langevin Dynamics.

• In Chapter 2, we present the DPD model and its energy conservation variant (DPDE).

We also give details about three existing representative numerical schemes for the DPDE integration: Shardlow's scheme (SSA), the most efficient scheme known so far but difficult to parallelize, and two other remaining options called SVV and SEM, both parallelizable. Finally, we present the most important contribution of this work: two new straightforwardly parallelizable numerical schemes integrating the DPDE. These two schemes are named SER and Hybrid, the last one being a merge of SER and SSA.

• We then proceed to Chapter 3 with a comparative study of the existing and new DPDE schemes on equilibrium systems. The energy preservation properties of all the schemes presented in the previous chapter are studied, and we show that no scheme manages to correctly preserve the energy and that energy drifts occur which affects the estimations of other observables. A correction ensuring the energy preservation is then devised, and we study the biases of the corrected schemes on the estimations of the average values of physical observables on equilibrium systems. We show that both new schemes are consistent, and that the accuracy of Hybrid is similar to the reference SSA scheme.

• In Chapter 4, we validate our new schemes on large scale, massively-parallel simulations of shocked systems. We show that both of our new schemes handle correctly such simulations and bring a substantial improvement compared to the existing parallelizable schemes by being able to correctly reproduce the shock propagation, where SEM, using the same parameters, displays non-physical phenomena.

• The last chapter of this thesis, Chapter 5, presents an attempt at accelerating DPDE sequential simulations with a variant of the SSA scheme. However, we show that this variant cannot be extended to parallel simulations, due to the apparition of artificial phenomena that do not vanish in the limit of infinitely small time steps.

Introduction to Molecular Simulation

We start this thesis by an introductive chapter presenting notions that will be used throughout the sequel. Molecular Dynamics and the tools to our disposal to perform numerical simulations of atomistic systems are presented in this chapter.

Molecular Dynamics (MD) consists in studying microscopic particle systems using computers simulations instead of experiments. In order to do that, a model of the particles interaction describing the behavior of the particles must be devised and then numerically integrated. These models often take the form of Differential Equations, thus requiring schemes to compute approximate solutions of these equations. The quality of these approximations need then to be assessed, through information about the difference between approximate and exact solutions, or information about the difference between their respective distributions.

Molecular Dynamics relies on a microscopic description of matter. This microscopic information allows to have qualitative information about the behavior of matter: in a given state of temperature and pressure, does a given crystal melt? Does a chemical reaction occur? In fact, from this microscopic description, it is also possible to extract macroscopic information about the system, especially quantitative information, using the theory of Statistical Physics. This allows to determine for instance the macroscopic pressure, temperature or heat capacity of the system.

The classical microscopic description of matter accurately describes atomistic systems for which quantum effects are negligible. We consider in this thesis coarse-grained models, where some microscopic phenomenons are represented in an average manner, sometimes by adding random fluctuations to the differential equations satisfied by the particles, thus turning Ordinary Differential Equations (ODEs, i.e deterministic) into Stochastic Differential Equations (SDEs). Similarly to ODEs, SDEs can be used both as a numerical microscope and to compute macroscopic averages. However, the error control for the integration of SDEs is more complex than for ODEs, and there are several ways of defining the error between the approximate and true solutions.

We start in Section 1.1 by a description of matter at its microscopic level and of some tools used to perform numerical simulations of its behavior. We then give in Section 1.2 elements of the theory of Statistical Physics, which links microscopic information to av-erage macroscopic quantities, and present methods using ergodic SDEs to compute these macroscopic averages otherwise impossible to compute by traditional brute-force methods. We continue by presenting some elements on the discretization of SDEs in Section 1.3. Finally, we apply in Section 1.4 the notions introduced in the previous sections to a prototypical dynamics called the Langevin dynamics. This dynamics is similar to the main dynamics considered in this thesis and is used as a paradigmatic example to introduce the notions that we use on the next chapters.

Microscopic description of matter

The fundamental ingredient for any molecular simulation is the description of the interaction between atoms. We consider in this section atomistic systems described by the position and velocities of the nucleus of the atoms, with particles interacting through conservative forces deriving from a potential function. The electron are taken into account in the potential and are not explicitely represented in the system.

We start in this section by an introduction to the classical description of matter and then describe the dynamics of these atomistic systems along with some tools to perform numerical simulations of such systems.

The classical microscopic description of matter

In the classical description of matter, particles are represented by their positions q ∈ R d , and their momenta p ∈ R d , where d denotes the physical dimension of the system (usually d = 3). These two quantities are called the degrees of freedom of the particle. The particle's mass is denoted by m ∈ R + and its velocity is obtained from the momenta and mass as v = p/m ∈ R d . In general, if we consider a system of N particles, we index each particle with an integer i ∈ [|1, N |] and denote the corresponding positions, momenta, masses and velocities by q i , p i , m i and v i respectively. The full set of variables describing a particle is the vector composed of its position and momenta, and is denoted by x i = (q i , p i ) ∈ R 2d . We consider in this section an isolated system composed of N particles. The configuration of the system is denoted by x = (q, p) ∈ R 2dN , where q = (q 1 , . . . , q N ) ∈ R dN and p = (p 1 , . . . , p N ) ∈ R dN are respectively the position and momenta vectors of the particles defining the system. The phase space, or set of configurations, denoted by X , is the ensemble of all the possible values that q and p can take. It can denote different ensembles (i.e we might have X R 2dN ) depending on how the boundaries of the system are defined (see the discussion on the boundary conditions later on). The mass matrix is denoted by M ∈ R N ×N (usually M = Diag (m 1 , . . . , m N )). The knowledge of the three vectors q, p and m is enough to completely define the configuration of the system under consideration.

The interactions between particles are accounted for by a function U that is called the potential function, or potential energy, depending on the positions q only. The kinetic energy of the particles is defined by a function K(p) = 1 2 p • M -1 • p. The total energy H(q, p) of the system is then defined as the sum of its kinetic and potential components:

H(q, p) = 1 2 p • M -1 • p + U(q). (1.1) q 1 , m 1 q 2 , m 2 q 3 , m 3 q 4 , m 4 q 5 , m 5 v 1 v 2 v 3 v 4 v 5
Figure 1.1 | Two-dimensional representations of particles in the classical description of matter: Each particle indexed by i is represented by its position q i , its mass m i and its momentum p i = m i v i (where v i is its velocity).

In the case where M = Diag (m 1 , . . . , m N ), the total energy writes

H(q, p) = N i=1 p 2 i 2m i + U(q).

Potential functions

The potential, or potential energy, U accounts for the interaction of the particles: the forces applied on each particle due to its interaction with the other particles are described by its derivatives. Ideally, it is computed with ab-initio computations [CDK + 03].

Ab-initio potentials are non-empirical potentials, i.e they do not depend on empirical parameters and are supposed to describe exactly the physics of the system. However, their computation is very expensive, and thus limits the use of ab-initio potentials to systems of small size (typically N = 100 -1000).

Because of the computational-cost of ab-initio potentials, empirical potentials are used in practice for the study of larger systems. Empirical potentials are obtained by assuming a functional form depending on parameters for the interactions. These parameters are often tuned by fitting them to ab-initio computations.

A simple example is the Lennard-Jones potential [START_REF] Jones | On the Determination of Molecular Fields. II. From the Equation of State of a Gas[END_REF], used mostly to describe the interactions of noble gases. Lennard-Jones (LJ) potentials are called pairwise potentials because the interactions only depends on the relative distance of each particle pair. The interaction energy of two Lennard-Jones particles writes

U(q) = 1≤i<j≤N v(r ij ), r ij = q i -q j , v(r) = 4ε σ r 12 - σ r 6 . (1.2)
The parameter ε represents the minimum energy of the interaction, and the parameter σ represents the characteristic distance at which the interaction energy reaches 0. The distance at which the energy reaches its minimum is 2 1/6 σ. An example of a Lennard-Jones interaction is given in Figure 1. described by the potential. For instance, when considering Argon atoms, Lennard-Jones parameters should be set to ε k B ≃ 119.8K and σ ≃ 3.405Å. LJ parameters have been fitted for many different materials (see for instance [START_REF] Wilhelm | Estimation of Lennard-Jones (6,12) Pair Potential Parameters from Gas Solubility Data[END_REF] for values of the LJ parameters for many liquids).

Because of their simplicity, Lennard-Jones potentials are very useful for testing purposes, when one wants to test a new numerical method to compute properties of an atomistic system. However, they do not describe accurately the atomistic interactions as soon as systems are more complex than noble gases. Many more elaborate potentials have been presented in the literature, like the Embedded Atom Model (EAM) potential [DB83, DB84, FBD86], which was initially elaborated to describe the interactions between metal atoms like copper. These more elaborate potentials are more accurate than basic potentials such as Lennard-Jones in certain cases, but their evaluation require more floating point operations and are thus more computationally expensive.

However, for some applications involving complex systems, common potentials like LJ or EAM are completely useless, and an accurate description of the physics requires the development of some very elaborate and extremely computationally expensive potentials with many parameters. This is the case of reactive systems like Hydrocarbons for instance, where potential like REBO [START_REF] Brenner | Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films[END_REF], VALBOND [START_REF] Root | Valence bond concepts applied to the molecular mechanics description of molecular shapes. 1. application to nonhypervalent molecules of the p-block[END_REF] or ReaxFF [START_REF] Van Duin | A Reactive Force Field for Hydrocarbons[END_REF] are used. In order to illustrate the complexity of these potentials, let us mention that ReaxFF uses 28 parameters and REBO more than 50.

Cut-off radius, shifted and splined potentials

In theory, the potential energy is computed by calculating the interaction energy of all the particle couples (i, j), which is extremely expensive. In addition, the interaction energy v usually tending fast enough towards 0 when r → +∞, most of the computed pairwise interaction energies are of negligible influence over the final value of U(q). Therefore, in practice, we introduce what is called a cut-off radius r cut limiting the range of the potential interactions, i.e

v cut (r) =        v(r), if r < r cut , 0, otherwise.
Introducing a cut-off radius alleviates the need for computing many negligible interactions, thus making the potential energy computation must faster, without great alteration of the final result. Usually, the cut-off radius is taken two or three times bigger than the characteristic length of the interaction described by v. E.g, if v is of the form (1.2), a typical value of the cut-off radius is usually taken as r cut ≃ 3.0σ.

In order to recover the continuity at r = r cut , we shift the potential by subtracting v(r cut ) to the pairwise interaction energy. Similarly, in order to also recover the continuity of the derivative of v, we add a term proportional to rr cut so that v ′ vanishes at r = r cut . We denote potentials with all the above corrections by cut, shifted and splined potentials. They read

v eff (r) = v(r) -v(r cut ) -v ′ (r cut ) (r -r cut ) , if r < r cut , 0, otherwise. (1.3)
Note that such a cut-off can be used for short-range interactions only. Coulomb interactions for instance are proportional to 1 r , which implies that the long-range interactions cannot be neglected and must be computed using different techniques that we do not consider in this thesis (see [START_REF] Tuckerman | Statistical Mechanics: Theory and Molecular Simulation[END_REF] for details about such techniques).

Boundary Conditions

Particles are not always allowed to visit the entire space and sometimes boundary conditions must be enforced, depending on what one wants to modelise. Let us denote by Q the set of all the possible positions.

• For simulations of systems in vacuum, one can define Q = R dN , and particles are allowed to visit the entire space.

• For some simulations, the particles have to be confined to a specified region of space Q R dN . The restriction can be performed by introducing a wall of fixed particles (or particles of infinite mass) on the borders of Q, thus preventing the particles inside to go out (see Figure 1.3). Note than in the situation depicted by Figure 1.3, we still have Q = R dN .

• In many molecular dynamics simulations nowadays, periodic boundary conditions are used. Periodic boundary conditions allow to avoid boundary effects, therefore allowing bulk systems to be modeled. With periodic boundary conditions, the simulation box is replicated in all directions, and particles interact not only with the particles inside the simulation box but also with their periodic image. illustrates these boundary conditions. In a simulation, only the positions of the particle of the thickened box (box 5) are computed, and the other particles are obtained by taking the associated periodic images. If L denotes the length of each simulation box, choosing L > 2r cut allows to compute short-range interactions (e.g LJ interactions) by considering only one perdiodic image for all the particles in the simulation box. In Figure 1.4 for instance, particle B only interacts with particle C and the periodic image of particle A, and not with particle A itself nor with the periodic image of particle C.

Note that periodic boundary conditions only apply for short-range interactions.

Long-range interactions are computed normally as the bulk condition of the system requires.

Dynamics of isolated systems

We have presented in Section 1.1.1 a microscopic description of matter, and given some details about how microscopic particle systems can be represented. We now consider the dynamics followed by the particles, i.e the differential equation satisfied by their positions and momenta.

Using the same notation as in Section 1.1.1, the time evolution of isolated systems follows the Hamiltonian dynamics

       q(t) = M -1 p(t), ṗ(t) = -∇U(q(t)), (1.4)
where q(t) represent the time derivative of q(t). The above equation is a rewriting of Newton's equation of motion

M • q(t) = -∇U(q(t)),
where q(t) represents the second time-derivative of q. Straightforward computations show that equation (1.4) preserves the total energy H of the system defined by (1.1):

dH dt = ṗ(t) • M -1 p(t) + q(t) • ∇U(q(t)), = -∇U(q(t)) • M -1 p(t) -M -1 p(t) • ∇U(q(t)), = 0.
The energy H is called the Hamiltonian energy, and in fact fully describes the dynamics. Indeed, (1.1) can be rewritten in the following form:

       q(t) = ∇ p H(q(t)), ṗ(t) = -∇ q H(q(t)).
In the sequel, unless otherwise stated, we denote by x(t) = (q(t), p(t)) the general solution of the Hamiltonian dynamics (1.4) and ϕ t (x) the solution of the dynamics equal to x ∈ X at time 0, where (ϕ t ) are the flow maps of (1.4).

Time discretization

Although the properties of the Hamiltonian dynamics are well-understood by now [START_REF] Arnold | Mathematical aspects of classical and celestial mechanics[END_REF], its solutions cannot be computed analytically for a general U. Therefore, we need to compute approximations of solutions using numerical schemes, also called discretizations of the dynamics. The integration of the Hamiltonian dynamics is also well-known (see [START_REF] Leimkuhler | Simulating Hamiltonian Dynamics[END_REF][START_REF] Hairer | Geometrical Numerical Integration -Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF] for instance) but we nevertheless give below some useful elements of its principle. For information about the discretization of Ordinary Differential Equations in a more general case, one can see [START_REF] Dormand | Numerical methods for differential equations: a computational approach[END_REF].

Consider the following general Ordinary Differential Equation (ODE):

ẋ(t) = f (x(t)), (1.5)
where f : R d → R d , with d ∈ N, is a locally Lipschitz function so that solutions of the above ODE are defined everywhere in R d in a time interval [0, T max ], with T max > 0. In the sequel, we denote by . any given norm on R d . Given an integration time T > 0, a numerical scheme computes an approximation of (x(t)) t∈[0,T ] at a finite number of times t 1 , . . . , t N it , where N it is the number of times when the solution is approximated. Usually, the approximation times are regularly spaced in [0, T ], and separated by a time step denoted by ∆t, i.e t i = i∆t and N it = ⌊T /∆t⌋. The accuracy of the discrete approximation of the solution depends on the time step, but every reasonable numerical scheme should recover the true analytical solution when ∆t goes to zero.

A numerical scheme is defined by a discretization procedure Φ ∆t : given a time step ∆t, it computes an approximate solution of the equation at a time t +∆t from a starting point

x(t) ∈ R d at time t, i.e Φ ∆t (x(t)) ≃ x(t + ∆t).
Let (x n ) n=1,..,N it be a sequence generated by the successive iterations of Φ ∆t :

x n+1 = Φ ∆t (x n ). (1.6)
When the sequence starts from x 0 = x, such a sequence can also be denoted by x n = Φ n ∆t (x). The quality of the approximation (x n ) can be estimated by looking at two aspects: the accuracy of the local estimation Φ ∆t compared to the analytical solution, and the error propagation caused by the successive use of the approximation Φ ∆t instead of the true analytical solution. These two aspects correspond respectively to the notions of consistency and stability.

Definition 1.1.1. [Stability] Consider two sequences (x n ) n=0,..,N it and (y n ) =0,..,N it defined by        x n+1 = Φ ∆t (x n ), y n+1 = Φ ∆t (y n ) + δ n ,
with (δ n ) n=0,..,N it an arbitrary sequence. A scheme is said to be stable if there exists C T > 0 independent of N it or ∆t alone (but which might depend on T = N it ∆t) such that max n=0,..,N it

y n -x n ≤ C T N it n=1 |δ n | .
In order to study the consistency of a numerical scheme, we have first to define what the local and global errors are. The global consistency, or convergence, of a numerical scheme is the fact that we recover the true solutions at any given finite integration time T when the time steps goes to 0, i.e for any finite integration time T > 0 and any x ∈ R d for which a solution to (1.5) exists, we have lim ∆t→0 E ∆t,T (x) = 0.

In order to obtain the convergence, we need the scheme to be stable and locally consistent.

Proposition 1.1.1. A locally consistent and stable numerical scheme is convergent.

The order of accuracy of a numerical scheme gives information on the convergence rate, i.e how fast the global error goes to 0 with the time step. It is defined through the local error: Definition 1.1.5. [Order of a numerical scheme] A numerical scheme is said to be of order ω if for any y ∈ R d for which a solution to (1.5) exists, we have

Φ ∆t (y) -ϕ ∆t (y) ≤ C y ∆t ω+1 ,
where C y > 0 does not depend on the time step ∆t. Proposition 1.1.2. Let Φ ∆t be a stable numerical scheme of order ω > 0. Then for any x ∈ R d for which a solution to (1.5) exists, there exists K T ,x > 0 independent of ∆t such that

max 0≤n≤N it ϕ n∆t (x) -Φ n ∆t (x) ≤ K T ,x ∆t ω .
Remark 1.1.1. If the prefactor K T ,x does not depend on ∆t, it can depend on the integration time T and the initial condition x. In some cases, K T ,x increases exponentially with T , thus severely limiting the range of integration times allowing for accurate estimations of x(T ).

Examples of discretizations of the Hamiltonian dynamics

The first numerical scheme that comes to mind for the integration of the Hamiltonian dynamics is the basic Explicit-Euler scheme (EE), which reads

       q n+1 = q n + ∆tM -1 p n , p n+1 = p n -∆t∇U(q n ).
The EE scheme is explicit and computationally fast, but is of order one of accuracy, and fails to conserve the energy by exhibiting energy drifts. Other schemes like Runge-Kutta schemes are of higher order but also suffer from energy drifts. Therefore, such schemes are not adapted to the integration of the Hamiltonian dynamics. Symplectic schemes are usually preferred to standard explicit schemes for general ODEs. They approximately preserve the energy up to a certain order for exponentially long times under some conditions on the trajectories and on the potential U ([HLW06], Chapters VI, VIII and IX). Indeed, possibly after some initial transient relaxation, the energy oscillates around an average value H ∆t = H 0 +K∆t ω for some K ∈ R and ω > 0, with the amplitude of the oscillations being proportional to ∆t 2 . This means that symplectic schemes do not yield energy drift, as can be seen in Figure 1.5. .5 | Energy preservation by the Velocity-Verlet scheme (symplectic, see (1.12)) for time steps ∆t * = 0.001 and ∆t * = 0.005 on an atomistic system described by a Lennard-Jones potential. After a transition time, the energy stabilizes around an average value close to its initial value. The difference between the two increases with the time step, as well as the amplitude of the oscillations around the average stabilized value.

Splitting techniques are often used to devise symplectic schemes. Indeed, the Hamiltonian flow is a symplectic mapping, and the composition of two symplectic mappings is still symplectic. In our case, the Hamiltonian dynamics can be separated in two subdynamics that are still Hamiltonian, containing respectively the time evolution of the positions and the momenta:

       q(t) = M -1 p(t), ṗ(t) = 0. ,        q(t) = 0, ṗ(t) = -∇U(q(t)).
(1.7)

Because each of the above two subdynamics is Hamiltonian, both of their flows are symplectic. In addition, we can analytically integrate these dynamics as

       q(t) = q 0 + tM -1 p 0 , p(t) = p 0 , ,        q(t) = q 0 , p(t) = p 0 -t∇U(q 0 ). (1.8)
Therefore, splitting schemes where each subscheme that corresponds to the above analytical integrations over a time ∆t are symplectic.

In a Trotter splitting, each part is integrated successively with a full time step ∆t, thus yielding two integration schemes, called Symplectic Euler algorithms. If the positions are updated first, we obtain the Symplectic-Euler scheme A, which reads

       q n+1 = q n + ∆tM -1 p n , p n+1 = p n -∆t∇U(q n+1 ). (1.9)
If the momenta are updated first, we obtain the Symplectic-Euler scheme B, which reads

       p n+1 = p n -∆t∇U(q n ), q n+1 = q n + ∆tM -1 p n+1 .
(1.10) Remark 1.1.2. Updating the momenta or the positions first yields similar results, and does not change the order of accuracy nor the computational cost of the algorithms.

Both SE schemes (1.9) and (1.10) are explicit and have an order 1 of accuracy, just like the EE scheme.

The order of the Symplectic Euler schemes is increased by symmetrizing the integration, i.e by integrating twice one equation with time step ∆t 2 , and integrating the other equation with time step ∆t in between the two half-integrations. This procedure, called a Strang splitting, yields two schemes called the Symmetric Symplectic Euler algorithms (SES), or the Stormër-Verlet algorithms. If the positions are updated first, we obtain

                   q n+ 1 2 = q n + ∆t 2 M -1 p n , p n+1 = p n -∆t∇U(q n+ 1 2 ), q n+1 = q n+ 1 2 + ∆t 2 M -1 p n+1 .
(1.11)

If the momenta are updated first, we obtain the Velocity-Verlet scheme, introduced in [Ver67], which reads

                   p n+ 1 2 = p n - ∆t 2 ∇U(q n ), q n+1 = q n + ∆tM -1 p n+ 1 2 , p n+1 = p n+ 1 2 - ∆t 2 ∇U(q n+1 ).
(1.12)

Both SSE schemes (1.11) and (1.12) are explicit, of order 2, and symplectic. The Velocity-Verlet scheme (1.12) is the most widely used scheme in the Molecular Simulation community for modeling NVE atomistic systems, because of its higher order accuracy and its long-time stability.

Going from a microscopic description to a macroscopic one

After having described the matter at the microscopic level, we explain in this section how macroscopic quantities can be computed from such a microscopic description. This is done according to the theory of Statistical Physics, which defines macroscopic quantities, also called physical observables, as expectations of microscopic quantities according to some probability measures. However, the downside of this approach is that physical observables are expressed as integrals with respect to a probability measure on spaces whose dimension is proportional to the number of particles of the system. The complexity of their computation therefore increases extremely fast with the number of particles considered, which makes standard spatial discretization methods useless. This is called the "Curse of Dimensionality". Brute-force deterministic methods being prohibitively expensive, one option is to use Stochastic Differential Equations (SDEs) to compute the observable expectations using ergodic limits. Therefore, under appropriate conditions, the many-dimensional integral corresponding to the observable expectation is transformed into a single dimensional time integral much simpler to compute.

We organize this section as follows: we start in Section 1.2.1 by giving some elements on Statistical Physics, and then describe in Section 1.2.2 how and under which conditions SDEs can be used to compute the observables expectations.

Introduction to Statistical Physics

In the Statistical Physics framework, thermodynamical observables are defined as functions of the configuration (q, p) ∈ X . ϕ : X -→ R, (q, p) -→ ϕ(q, p).

(1.13) However, the value ϕ(q, p) does not give any information on the system because it is the value of ϕ for one specific configuration (q, p) out of the infinite set of possible configurations X . Rather, the macroscopic value corresponding to the observable ϕ, denoted by ϕ , is defined in an average manner.

Definition 1.2.1. [Average value of an observable] Consider an observable ϕ defined by (1.13). The physical value of ϕ is defined as the expectation of its microscopical values for configurations distributed according a probability distribution π. I.e

ϕ π = E π [ϕ(q, p)] = ϕ(q, p)π(dq, dp), (1.14)
The probability distribution π is called the thermodynamic state of the system, also called thermodynamic ensemble, statistical ensemble, macrostate or macroscopic state of the system.

In the remainder, we denote by D a set of smooth functions, composed of infinitely derivable functions whose derivatives grow at most polynomially. Consider a positive function W : X → [1, +∞) and define the following weighted space

L ∞ W = ϕ measurable, ϕ W L ∞ < +∞ , (1.15)
and the corresponding norm

ϕ L ∞ W = ϕ W L ∞ . (1.16)
Consider now the family of functions W m being defined for all m > 0 as W m (x) = 1 + |x| m . The ensemble D is then defined as

D = ϕ ∈ C ∞ (X , R) ∀m ∈ N, ∀l ∈ N N d s.t. |l| = m, ∃k ∈ N, ∂ l ϕ L ∞ W k < +∞ , (1.17)
where

∂ l ϕ = ∂ l 1 ∂x l 1 1 . . . ∂ l N d ∂x l N d N d
and N d is the dimension of the system (e.g considering a system composed of N particles, N d = dN if x = q, or N d = 2dN if x = (q, p)). In the remainder of this thesis, observables, unless otherwise stated, are considered to be elements of D.

Remark 1.2.1. The choice of the family of functions (W m ) m∈N : x → (1 + |x| m ) in the definition of D given by (1.17) is arbitrary. Many results presented in this chapter still hold when (W m ) m∈N is any family of positive scale functions verifying W 0 ≤ W 1 ≤ . . . ≤ W m ≤ . . .

Thermodynamic ensembles

According to Definition 1.2.1, the physical value of a given observable depends on the thermodynamic ensemble of the system considered. Several thermodynamic ensembles can be considered, and we list below two common choices:

• The microcanonical ensemble (NVE): This ensemble describes a system in which the quantity of matter N , the volume V and the energy E are conserved. Note that E here denotes the Hamiltonian energy defined by (1.1). This is the ensemble in which any closed system at equilibrium is. Its corresponding distribution writes µ mc,E 0 (dq, dp) = Z -1 mc,E 0 δ H(q,p)-E 0 dqdp.

(1.18)

The distribution δ H(q,p)-E 0 can be explicitly constructed as the limit of a uniform measure on the set of configurations (q, p) ∈ X of energy slightly higher than E 0 .

Let C(E) be the subset of X of the configurations with energy E 0 , i.e

C(E 0 ) = { (q, p) ∈ X | H(q, p) = E 0 } ,
and define

N ∆E (E 0 ) = (q, p) ∈ X E 0 ≤ H(q, p) ≤ E 0 + ∆E .
The set N ∆E (E 0 ) is endowed with a uniform measure, and we define δ H(q,p)-E 0 through the expectations of observables ϕ ∈ D as

C(E) ϕ(q, p)δ H(q,p)-E 0 (dqdp) = lim ∆E→0 1 ∆E N ∆E (E 0 ) ϕ(q, p)dqdp.
The constant Z mc,E 0 is a normalization constant ensuring that µ mc,E 0 is a probability measure, i.e

Z mc,E 0 = C(E 0 )
δ H(q,p)-E 0 (dqdp).

For more details on the construction of δ H(q,p)-E 0 , we refer to Section 1.2.3 of [START_REF] Lelièvre | Free Energy Computation : A Mathematical Perspective[END_REF].

• The canonical ensemble (NVT): This ensemble describes a system at equilibrium coupled with an infinite energy reservoir. The system can exchange energy with the reservoir, which thus acts as a thermostat keeping the temperature T of the system constant. The quantity of matter N and the volume V are also conserved. The associated probability measure reads µ β (dq, dp) = Z -1 β e -βH(q,p) dqdp, (1.19) where β = 1 k B T , with T the temperature of the system and k B Boltzmann's constant. The term Z β is a normalization constant called the partition function which is defined as

Z β = X e -βH(q,p) dqdp.
The expression (1.19) of µ β is defined as the probability distribution maximizing, under constraints of fixed average energy value, the statistical entropy. Consider a probability measure with density ρ. The statistical entropy G(ρ) corresponding to this measure is defined as

G(ρ) = - X ρ(q, p)lnρ(q, p)dqdp.
It is shown in [START_REF] Balian | From Microphysics to Macrophysics[END_REF] that the expression given by (1.19) is the probability measure maximizing G under the constraint of fixed average energy value. I.e, assuming that e -βU ∈ L 1 (X ), µ β is the solution of the following optimization problem:

sup ρ∈L 1 (X ) G(ρ) ρ ≥ 0, X ρ(q, p) = 1, X H(q, p)ρ(dqdp) = E 0 .
Therefore, β depends on E 0 . More details about how to choose E 0 are given later on when the equivalence of ensembles is discussed.

There are many more thermodynamic ensembles than the two examples presented. In fact, the choice of ensemble amounts to choosing which quantities are preserved exactly, and which are preserved in average. The resulting macrostate is then obtained by maximizing the entropy under the chosen constraints, in a way similar to the derivation performed for the canonical ensemble µ β . Details of the derivation are given in Section 1.2.3.3 of [START_REF] Lelièvre | Free Energy Computation : A Mathematical Perspective[END_REF], and are not reported here.

As we saw previously, thermodynamic ensembles are constructs derived from a "physical" choice of fixed quantities. We therefore expect the average values they predict to be consistent one with one another. This is true in the limit of infinite systems (N → ∞), as long as neighborhoods of phase transitions are avoided [START_REF] Fisher | The free energy of a macroscopic system[END_REF]. Details on the method for passing between ensembles, outlined below, are given in [LL80,[START_REF] Bernardin | Thermodynamics and non-equilibrium macroscopic dynamics of chains of anharmonic oscillators[END_REF].

Consider as an illustration the microcanonical and canonical ensembles. Consider an observable ϕ ∈ D depending only on a finite number k ∈ N of variables denoted by (x 1 , . . . , x k ). It can be shown that lim

N →+∞ E µ (N ) mc,E 0 [ϕ(x 1 , . . . , x k )] -E µ (N ) β [ϕ(x 1 , . . . , x k )] = 0, provided E µ β [H] = E 0 , where µ (N )
mc,E 0 and µ (N ) β are respectively the microcanonical and canonical measure corresponding to a system with N particle, and H is the Hamiltonian energy of the system. The equality E µ β [H] = E 0 means that the average value of H in the canonical ensemble is equal to the chosen value E 0 in the microcanonical ensemble. Rigorous details about the equivalence of statistical ensembles are given in [START_REF] Ruelle | Statistical Physics: Rigourous Results[END_REF].

Example 1.2.1. Define a temperature T such that E µ β [H] = E 0 . Consider the virial pressure of the system defined in [START_REF] Allen | Computer Simulation of Liquids[END_REF] as

P W (q) = - 1 3 N i=1 q i • ∇ q i U(q).
If the potential U involves only short-range interactions, P W is expected to be an infinite sum of terms depending only on a finite number of variables. Indeed, short-range interactions are limited in range, and each particle is expected to have a finite number of neighbors interacting with it (see the cut-off radius section of Section 1.1.1). Therefore, for each term of the above sum, we have

∀i ∈ N, lim N →+∞ E µ (N ) β q i • ∇ q i U -E µ (N ) mc,E0 q i • ∇ q i U = 0.
This implies that the virial pressure has the same average value in the canonical and microcanonical ensembles for infinite systems, provided that

E µ β [H] = E.

Curse of dimensionality

We have seen that Statistical Physics allows to obtain macroscopic information from a microscopic description through the expectations of observables. However, these expec-tations are very high-dimensional integrals that are prohibitively expensive to compute. Standard discretizations procedures are therefore impossible to use. In addition, it is usually impossible to sample the different macrostates using i.i.d variables. Considering for instance the canonical measure given by (1.19), its marginal related to the momenta is a Gaussian distribution which is easy to sample, but the position-related marginal cannot be sampled directly in the general case. Two other options exist to compute observables such as defined by (1.14):

• One option is to compute Markov chains realizations and to obtain the thermodynamic averages as ergodic limits. One example is the Metropolis algorithm [MRR + 53] computing realizations of ergodic Markov chains for the position marginal of the canonical macrostate.

• One other option is to use Markov processes defined as solutions of Stochastic Differential Equations ergodic with respect to the target probability measure.

In this thesis, we only consider the second option and detail it in the sequel.

Computing physical observables with Stochastic Differential Equations

Let us consider the following general diffusion SDE:

dx t = f (x t )dt + b(x t )dW t , (1.20)
where

x t ∈ X ⊂ R N d , (W t ) t∈R is a standard m-dimensional
       f (x) -f (y) + b(x) -b(y) ≤ C x -y , ∀(x, y) ∈ X 2 , t ≥ 0, f (x) + b(x) ≤ C (1 + x ) , ∀x ∈ X , t ≥ 0.
The above conditions on f and b ensure the existence and uniqueness of the solution of (1.20) [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF][START_REF] Arnold | Stochastic Differential Equations[END_REF]. They also ensure the continuous dependence of the solutions with respect to the initial conditions in the following sense: denoting x t and x (n) t the solutions of (1.20) starting respectively at x 0 = y and x

(n) 0 = y n , lim n→∞ y n -y = 0 =⇒ lim n→∞ E sup 0≤t≤T x (n) t -x t 2 = 0.
Because this section is only an overview on SDEs and their discretizations, we consider all the observables ϕ to be smooth functions whose derivative grow at most polynomially, i.e ϕ ∈ D. We therefore give in the sequel results that are not optimal: most of them can be extended to larger function spaces. For more optimal results, we refer for instance the reader to [START_REF] Rey-Bellet | Ergodic Properties of Markov Processes[END_REF] or [START_REF] Kliemann | Recurrence and Invariant Measures for Degenerate Diffusions[END_REF] for the ergodic theory of SDEs and to [START_REF] Arnold | Stochastic Differential Equations[END_REF] for the general theory of SDEs.

The evolution operator of an SDE, denoted by S t , is the stochastic counterpart of the flow map of an Ordinary Differential Equation. It maps a function ϕ ∈ L ∞ (X ) to the function mapping x ∈ X to the expectation of ϕ on the solution x t at time t of the SDE starting from configuration x 0 = x.

Definition 1.2.3. [Evolution operator]

The evolution operator at time t ∈ R of the dynamics given by (1.20), denoted by S t , is defined as

             L ∞ (X ) → L ∞ (X ), ϕ → S t ϕ :        X → R, x → E [ϕ(x t ) | x 0 = x] ,
(1.21)

where E denotes the expectation over all realizations of the Brownian motion in (1.20).

Under Condition 1.2.2, the family of evolution operators (S t ) t≥0 satisfies, amongst others, the following properties for all ϕ ∈ D, x ∈ X , t ≥ 0 and s ≥ 0:

             S s S t ϕ = S t+s ϕ, semigroup property, ϕ ∈ D ⇒ S t ϕ ∈ D, smoothness and polynomial growth of the solution processes [GS72], lim t→0 S t ϕ(x) = ϕ(x), No evolution at time t = 0. Definition 1.2.4.
[Infinitesimal generator of an SDE] The infinitesimal generator or generator of the SDE, denoted by L, is defined as the generator of the semigroup (S t ) t≥0 . Given a function ϕ regular enough, L is defined as

Lϕ = lim t→0 S t ϕ -ϕ t . (1.22)
Its general expression is

L = f (x) • ∇ x + 1 2 b(x)b(x) T : ∇ 2 x = N d i=1 f (x) i ∂ x i + 1 2 N d i,j=1 b(x)b(x) T i,j ∂ x i ∂ x j .
The operator L is defined where the limit (1.22) exists. For functions ϕ inside its domain of definition, we have S t ϕ = e Lt ϕ for any t ≥ 0. Finally, we define the probability transition P t (x, A) of the dynamics as the probability for its solution to reach A ⊂ X , starting in t = 0 from x ∈ X :

P t (x, A) = E 1 A (x t ) x 0 = x .
where 1 A is the indicator function of the set A, equal to 1 if x ∈ A and 0 elsewhere.

Ergodic SDEs

Thermodynamic expectations such as (1.14) can be computed using ergodic SDEs. The ergodicity of an SDE for a given probability measure π ensures that the time integral of any function of its solutions converges towards the expectation of this function for π.

Definition 1.2.5. [Ergodicity] Consider a diffusion SDE given by (1.20). We say that the SDE is ergodic with respect to π if we have, for π-almost any x

0 ∈ X , ∀ϕ ∈ L 1 (π), lim t→+∞ 1 t t 0 ϕ(x t )dt = X ϕdπ a.s. (1.23) Remark 1.2.2.
In fact, the above definition has to be contrasted with with the geometric ergodicity, describing the convergence of the law of the solution processes to the ergodic probability measure π. However, the geometric ergodicity does not imply the ergodicity as defined in Definition 1.3.4.

If (1.23) holds, we only have to compute the time average, i.e the left-hand side of the above relation, which is a one-dimensional integral much faster to compute than the N d -dimensional spatial integral of the right-hand side.

In order for (1.23) to hold, we first need the probability measure π to be invariant.

Definition 1.2.6. [Invariant measure] A probability measure π is said to be invariant for the dynamics if ∀ϕ ∈ D ∞ , X Lϕdπ = 0. (1.24)
where L is the generator of the dynamics and D ∞ is the set of C ∞ functions with compact support, also called the set of test functions. Equivalently, we say that π is invariant if

L * π = 0, (1.25)
where L * is the adjoint of L in L 2 (X , R), also called the Fokker-Planck operator in physics. It is defined as

∀(ϕ, ψ) ∈ D 2 ∞ , X Lϕ(x)ψ(x)dx = X ϕ(x)L * ψ(x)dx,
If we have the existence of an invariant probability measure π with positive density for the Lebesgue measure, then we only need some regularity results on the generator L in order to obtain the ergodicity. If L is elliptic, i.e if b(x)b(x) T is positive definite for all x ∈ X , then the associated dynamics is ergodic with respect to π. However, if we do not have the ellipticity of L, we can still obtain the ergodicity by proving that L is hypoelliptic. An hypoelliptic operator is an operator satisfying the Hörmander condition: Definition 1.2.7. [Hörmander condition] Assume that there exists M ∈ N and a family of operators (X i ) i=1..M such that L can be rewritten as

L = M i=1 X * i X i + X 0 ,
and suppose that

{X i } M i=0 , [X i , X j ] M i,j=0 , [[X i , X j ], X k ] M i,j,k=0
, ... has full rank at every point x ∈ X , where [•] is the commutator operator defined by

[X, Y ] = XY -Y X. (1.26)
If L satisfies the above Hörmander condition, then L is hypoelliptic (Corollary 7.2 in [START_REF] Rey-Bellet | Ergodic Properties of Markov Processes[END_REF]). Note that all elliptic operators are also hypoelliptic. The following result was proved by Kliemann in [START_REF] Kliemann | Recurrence and Invariant Measures for Degenerate Diffusions[END_REF].

Theorem 1.2.1. Consider an SDE admitting an invariant probability measure π with a positive density for the Lebesgue measure, and whose generator L is hypoelliptic. Then the dynamics is ergodic with respect to π.

However, obtaining the existence of an invariant probability measure is nontrivial. Sometimes, dynamics are devised with the aim of being ergodic for a target thermodynamic ensemble. Therefore, they are constructed so this target distribution is invariant, which ensure de facto the existence of an invariant probability measure. The fact that is has a positive density or not can then be obtained from its definition. In the general case however, when we do not know any invariant probability measure π, it is still possible to obtain the existence of π by using Lyapunov functions.

Definition 1.2.8. [Lyapunov functions] A Lyapunov function W is a function verifying ∀x ∈ X , W (x) ≥ 1, and lim |x|→+∞ W (x) = +∞.
Consider the two following conditions:

• (C1): The dynamics is irreducible, i.e for any x ∈ X and for any open sets A ⊂ X , there exists t 0 > 0 such that we have P t 0 (x, A) > 0. In words, this means that the SDE can reach with positive probability any open set A of X starting from any initial condition x ∈ X in time t 0 .

• (C2): The transition probability function P t (x, dy) admits a density p t (x, y) which is a smooth function of (x, y), i.e p t ∈ D. In particular, this means that the evolution operator S t has a smoothing effect and maps bounded measurable functions into bounded continuous functions. We say that S t is Strong-Feller.

The following theorem (Theorem 8.3 in [START_REF] Rey-Bellet | Ergodic Properties of Markov Processes[END_REF]) ensures the existence of an invariant probability measure for the dynamics.

Theorem 1.2.2. Consider an SDE, whose generator is denoted by L, satisfying (C1) and (C2).

Assume that there exists b > 0, c > 0, a compact set K and a Lyapunov function W such that

∀x ∈ X , LW (x) ≤ -c + b1 K (x).
Then there exists an unique invariant probability measure for the dynamics. This measure admits a smooth everywhere positive density.

Conditions (C1) and (C2) ensures the uniqueness and the smooth everywhere positive density of the invariant probability measure if it exists but do not ensure its existence (Proposition 8.1 in [START_REF] Rey-Bellet | Ergodic Properties of Markov Processes[END_REF]). The Lyapunov inequality is therefore necessary to ensure its existence.

The two first conditions of Theorem 1.2.2 are not easy to obtain when the noise is too degenerate (see the DPD/DPDE dynamics in the next chapters). However, there are some tools that allow to obtain them:

• The irreducibility can be obtained with control theory. A control system of the dynamics given by (1.20) is defined as the following deterministic dynamics:

ẋ(t) = f (x(t)) + b(x(t))u(t).
The function u is piecewise constant, and is called a control. A point y ∈ X is said to be accessible from x ∈ X in time t ≥ 0 if there exists a control u such that the above deterministic dynamics has solutions (x

(u) s ) s∈[0,t] satisfying x (u) 0 = x and x (u) t = y.
The set of accessible points from x in time t is denoted by A t (x) ⊂ X . The Stroock-Varadhan theorem, quoted as Theorem 6.1 in [START_REF] Rey-Bellet | Ergodic Properties of Markov Processes[END_REF], implies that

supp (P t (x, •)) = A t (x).
Therefore, the irreducibility of the dynamics can be obtained by finding the appropriate controls for any given points (x, y) ∈ X 2 .

• The smoothness of the density of the probability transition function is obtained by proving some regularity properties on the generator L of (1.20). In particular, if L is elliptic then (C2) is verified. If L fails to be elliptic, one option is to prove that it is hypoelliptic, i.e that it satisfies the Hörmander condition (see Definition 1.2.7).

If L is hypoelliptic, then condition (C2) is satisfied.

Therefore, once the irreducibility of the dynamics, the hypoellipticity of its generator and the Lyapunov condition are proved, we can apply Theorem 1.2.2 in order to obtain the existence of a positive invariant probability measure for the dynamics. Theorem 1.2.1 then gives us the ergodicity of the dynamics for this probability measure.

Numerical methods for Stochastic Differential Equations

We have presented in Section 1.1 a microscopic description of matter with deterministic differential equations, and explained in Section 1.2 how macroscopic quantities can be computed from microscopic information using ergodic SDEs. We now consider how to compute approximate solutions of SDEs.

In addition of being tools to compute statistical physics expectations, Stochastic Differential Equations (SDEs) are often used in Coarse-Graining models as a modeling tool. Indeed, in some cases, microscopic phenomenons can be represented in an average manner by adding random fluctuations to the dynamics of the particles coupled with some dissipation terms [START_REF] Mori | Transport, collective motion and brownian motion[END_REF][START_REF] Zwanzig | Nonlinear generalized langevin equations[END_REF]. Therefore, SDEs can be also used to provide some description of the matter behavior, often for systems of length and time scales larger than those described by empirical MD, but still at the microscopic level. Note that the coarsegraining procedures are often based on intuition and therefore deterministic Hamiltonian simulations of fully atomistic systems should always serve as reference simulations to verify the correctness of the coarse-graining process.

Most of the time we do not know how to compute analytical solutions of SDEs, and we must rely on discretization procedures to generate Markov chains of finite size that approximate the analytical solutions in some time interval. Doing so we add errors to the estimations that are not related to the finite length of the Markov chains. Estimates of these errors are of paramount importance.

We start in Section 1.3.1 by giving elements on the discretization of Stochastic Differential Equations, and present some tools to obtain standard error estimates (weak/strong errors on finite time intervals). We then consider in Section 1.3.2 the case of errors on the invariant measure and give some tools to estimate it. Finally, we give in Section 1.3.3 some splitting techniques useful to devise higher order numerical schemes.

Elements on the discretization of Stochastic Differential Equations

A discretization of equation (1.20) is a procedure applied to a given initial condition x ∈ X approximating the solution at time ∆t of the SDE starting from x. The parameter ∆t is called the time step of the scheme.

Consider an integration time T > 0, a time step ∆t > 0 and define N it = ⌊T /∆t⌋. Denote by Φ ∆t (x, G) the result of a given one-step discretization of (1.20), starting from x ∈ X , using a time step ∆t and with G being a random variable. Let us denote by (

x n ) n∈[|0,N it |]
the random sequence defined iteratively by

x n+1 = Φ ∆t (x n , G n ), (1.27) where (G n ) n∈[|0,N it |] are i.i.d random variables.
The initial condition of the above sequence is taken equal to the initial condition of the analytical solution (x t ) we want to approximate, i.e

x 0 = x 0 ∈ X .
The sequence (x n ) n=0,..,N it is aimed to be an approximation of the exact solution (

x n∆t ) n∈[|0,N it |]
of (1.20). In the remainder of this thesis, we only consider one-step schemes (i.e schemes defined by (1.27)) with Gaussian increments

G n law = W (n+1)∆t -W n∆t ∆t ∼ N (0, Id), n = 0, . . . , N it .
The evolution operator, or transition operator, P ∆t associated with a numerical scheme is the discrete counterpart of the evolution operator S ∆t = e ∆tL . In words, P ∆t ϕ(x) gives the average value of one observable after one iteration of the scheme applied to configuration x, with time step ∆t.

Definition 1.3.1. [Evolution operator] Consider a discretization of the dynamics given by

(1.20), with its discretization procedure denoted by Φ ∆t . Its associated evolution operator P ∆t is defined as

             L ∞ (X ) → L ∞ (X ), ϕ → P ∆t ϕ :        X → R x → E G [ϕ (Φ ∆t (x, G))] ,
(1.28)

where E G is the expectation over all the random variables involved in the computation of

Φ ∆t (x, G), denoted by G.
An approximation (x n ) defined by (1.27) for a specific realization of the random sequence (G n ) is called a realization of the numerical scheme Φ ∆t .

The random sequence (x n ) is in fact a Markov chain whose probability transition function is given by P ∆t . A time-homogeneous Markov chain (x n ) n∈N is a sequence in which each x n+1 is sampled from x n according to a transition probability distribution P(x n , dx) and independently from the previous samples (x 0 , . . . , x n-1 ). In our case the transition probability function, abusively denoted by P ∆t (x, A) with x ∈ X and A ⊂ X a measurable subset of X , is obtained by

P ∆t (x, A) = P ∆t 1 A (x), (1.29)
where 1 A still denotes the indicator of the set A. We say that the scheme is irreducible with respect to the Lebesgue measure when any open subset of X of nonzero Lebesgue measure can be reached from any point of X , similarly to the continuous case of Section 1.2.2.

Definition 1.3.2. A discretization whose evolution operator is denoted by P ∆t is irreducible with respect to the Lebesgue measure if for all x ∈ X and for every subset A ⊂ X of nonzero Lebesgue measure, there exists N 0 ∈ N such that

P N 0 ∆t 1 A (x) > 0.
Similarly to the continuous case of Section 1.2.2, we can define invariant measures and ergodic properties for the Markov chain generated by P ∆t . The invariant probability measure π ∆t of a scheme is a probability measure that is preserved by the evolution operator P ∆t . It means that if x n is distributed according to π ∆t , then x n+1 also is.

Definition 1.3.3. Consider a discretization of the SDE given by (1.20) whose evolution operator is denoted by P ∆t . A probability measure π ∆t is said to be an invariant for

P ∆t if it satisfies ∀ϕ ∈ L ∞ (X ), X P ∆t ϕdπ ∆t = X ϕdπ ∆t .
Similarly to the continuous case, we can define the ergodicity of numerical schemes with respect to a given probability measure: Definition 1.3.4. Consider a time step ∆t, and a numerical scheme whose evolution operator is denoted by P ∆t . A numerical scheme is said to be ergodic with respect to a probability measure π ∆t if we have, for π ∆t -almost any

x 0 ∈ X , ∀ϕ ∈ L 1 (π ∆t ), lim N it →∞ 1 N it N it n=0 ϕ(x n ) = X ϕdπ ∆t a.s.
The following result has been proved in Chapter 17 in [START_REF] Meyn | Markov Chains and Stochastic Stability[END_REF], and give conditions in order for the Markov chain generated by P ∆t to be ergodic with respect to a probability measure π ∆t . Proposition 1.3.1. Consider a scheme whose evolution operator is denoted by P ∆t . Assume that it is irreducible and that it admits an invariant probability measure π ∆t . Then it is ergodic with respect to π ∆t .

However, Proposition 1.3.1 requires the existence of a invariant probability measure π ∆t which, similarly to the continuous case with SDEs, might be a nontrivial result. However, we have some tools to prove the existence of an invariant probability measure with the following theorem, also ensuring its uniqueness.

Theorem 1.3.2. Consider a discretization whose evolution operator is denoted by P ∆t . Assume the following Lyapunov condition to hold: there exists a function

W : X → [1, +∞), K ≥ 0 and 0 < γ < 1 such that ∀x ∈ X , P ∆t W (x) ≤ γW (x) + K.
Also assume the following minorization condition: there exists a probability measure ν and

0 < α < 1 such that ∀y ∈ X , inf x∈C P ∆t (x, dy) ≥ αν(dy), with C = {x ∈ X | W (x) ≤ R}, for any R > 1 + 2K/(1 -γ).
Then there exists an invariant probability measure π ∆t satisfying

X W dπ ∆t < +∞.
Finally, there exists C > 0 and 0 < r < 1 such that

∀ϕ ∈ L ∞ W , ∀n ∈ N, P n ∆t ϕ - X ϕdπ ∆t L ∞ W ≤ Cr n ϕ - X ϕdπ ∆t L ∞ W ,
where L ∞ W is defined by (1.15). The proof of this theorem can be found in [START_REF] Hairer | Yet Another Look at Harris' Ergodic Theorem for Markov Chains[END_REF] (and is discussed in Section 2.4 of [START_REF] Lelièvre | Partial differential equations and stochastic methods in molecular dynamics[END_REF]) and is not reported here.

Remark 1.3.1. Notice that Theorem 1.3.2 not only gives us the existence of an invariant probability measure, but also gives us the exponential convergence in law of P ∆t towards π ∆t . This notion of convergence is called the geometric ergodicity, already mentioned in Remark 1.2.2.

Error estimates for the discretization of SDEs

Similarly to ODEs, once a discretization procedure is devised, we want to know how "good" it is at describing the solutions of the continuous dynamics it approximates. In opposition to the deterministic case, a stochastic discretization can be accurate in more than one way: one can be interested only in the accuracy on the average behavior of the solutions, or one could also require trajectorial accuracy, i.e that each realization of the dynamics is correctly approximated by the scheme. We can therefore define two types of discretization errors:

• Weak errors: Weak errors estimate how accurate the discretization is on the average behavior of the solutions. They thus compare the expectations of the analytical solutions and their discrete approximations.

Definition 1.3.5. [Weak errors]

There exists ω > 0 such that for any test function ϕ ∈ D ∞ and for any integration time T , we have K ϕ,T > 0 and ∆t * ϕ,T > 0 that satisfy, for any time step 0 ≤ ∆t ≤ ∆t * ϕ,T ,

sup 0≤n≤⌊ T ∆t ⌋ E [ϕ(x n )] -E [ϕ(x n∆t )] ≤ K ϕ,T ∆t ω , (1.30)
where the expectation is on the realizations of the Brownian motion in (1.20) for the analytic expectation, on the realizations the random variables G n ∼ N (0, Id) involved in the computations of x n with 0 ≤ n ≤ T ∆t for the discrete one and also on the initial conditions for both.

Equation (1.30) is called a weak error estimate, and ω is called the weak order of the numerical scheme P ∆t . A scheme is weakly consistent if it has weak errors of order ω > 0.

• Strong errors in L k norm: Strong errors are related to the trajectorial accuracy of the discretization and estimate how well the behavior per-realization of the dynamics is reproduced. In order to define them, we must construct a stochastic process xt that interpolates the discrete sequence (x n ). This process can be defined as

xt = x k + t k∆t f (x k )ds + t k∆t b(x k )dW s . (1.31)
We notice that xt is a piecewise linear function with random fluctuations passing through every points of (x n ) (see Figure 1.6 for an example of a process equal in law to xt ). In practice, one can also interpolate (x n ) by a piecewise constant function.

Definition 1.3.6. [Strong errors in L k -norm]
There exists ω > 0 such that for any integration time T , there is K T > 0 and ∆t * T > 0 for which, for any time step

0 ≤ ∆t ≤ ∆t * T , sup 0≤n≤⌊ T ∆t ⌋ E | xn∆t -x n∆t | k 1 k ≤ K T ∆t ω . (1.32)
Equation (1.32) is called a strong error estimate in L k norm, and ω is called the strong order in L k norm of the numerical scheme P ∆t . A scheme is strongly L k -consistent if it has strong errors in L k norm of order ω > 0.

Strong errors are more difficult to estimate than weak errors, both theoretically and numerically. In addition, the strong consistency of order ω > 0 implies the weak consistency of order ω. Indeed, we have

E[ϕ(x n )] = E[ϕ( xn∆t )], and therefore, for any ϕ ∈ D ∞ , n ∈ N and k > 1, |E [ϕ(x n )] -E [ϕ(x n∆t )]| = |E [ϕ( xn∆t ) -ϕ(x n∆t )]| , ≤ ∇ϕ L ∞ E | xn∆t -x n∆t | , ≤ ∇ϕ L ∞ E | xn∆t -x n∆t | k 1 k , ≤ ∇ϕ L ∞ sup n=0,..,N it E | xn∆t -x n∆t | k 1 k ,
where we have applied the mean value theorem to go from the first to the second line and the Cauchy-Schwarz inequality to go from the second to the third line.

Example 1.3.1. [The Euler-Maruyama discretization] The first example of a discretization of a stochastic dynamics is the Euler-Maruyama scheme (EM). The EM scheme is the stochastic counterpart of the Explicit-Euler scheme, and approximates (1.20) by

Φ ∆t (x, G) = f (x)∆t + b(x) √ ∆tG, (1.33)
with G ∼ N (0, Id). The Euler-Maruyama scheme is of weak order one and of strong-order 1/2 in the general case [START_REF] Milstein | Stochastic Numerics for Mathematical Physics[END_REF].

-

Obtaining weak error estimates

In order to obtain weak error estimates, one can use the stochastic counterpart of Proposition 1.1.2, allowing to obtain the weak order of a numerical scheme from a local error approximation. Such a process has been introduced by by Mil'shtein in 1986 (see Theorem 2 in [START_REF] Mil'shtein | Weak approximation of solutions of systems of stochastic differential equations[END_REF] or Theorem 2.1 in [START_REF] Milstein | Stochastic Numerics for Mathematical Physics[END_REF]) and we report the result below.

Theorem 1.3.3. Consider the SDE given by (1.20) whose generator is denoted by L, and where the functions f and b satisfy Condition 1.2.2. Consider a numerical discretization of L whose evolution operator is denoted by P ∆t . Let us assume that, for all ϕ ∈ D, P ∆t can be expanded as

P ∆t ϕ = ϕ + ∆tLϕ + ∆t 2 2 L 2 ϕ + . . . + ∆t ω ω! L ω ϕ + ∆t ω+1 r ∆t,ϕ , (1.34)
where r ∆t,ϕ remains uniformly bounded for small ∆t in some weighted space

L ∞ W m for m ∈ N, i.e ∃(m, C ϕ , ∆t * ) ∈ N × (R * + ) 2 , sup 0<∆t≤∆t * r ∆t,ϕ L ∞ W m < C ϕ .
Then the numerical scheme P ∆t is of weak order at least ω. where the remainder term r∆t,ϕ is also uniformly bounded in some weighted space L ∞ W m for some m ∈ N.

As an illustration of the procedure to obtain expansions such as (1.34), we prove below that the Euler-Maruyama scheme (1.33) satisfies (1.34) at the order ω = 1. Proposition 1.3.4. The Euler-Maruyama discretization of the SDE given by (1.20) satisfies the conditions of Theorem 1.3.3 for ω = 1. This means that its evolution operator P EM ∆t verifies ∀ϕ ∈ D, P EM ∆t ϕ = ϕ + ∆tLϕ + ∆t2 r ∆t,ϕ , where r ∆t,ϕ stays uniformly bounded for small ∆t in some weighted space

L ∞ W m , i.e ∃(m, C ϕ , ∆t * ) ∈ N × (R * + ) 2 , sup 0<∆t≤∆t * r ∆t,ϕ L ∞ W m < C ϕ .
Proof. The discretization procedure of the EM scheme for the dynamics (1.20) writes

Φ EM ∆t (x, G) = x + ∆tf (x) + √ ∆tb(x)G,
where G ∼ N (0, Id). Denote by h ∆t (x, G) = ∆tf (x) + √ ∆tb(x)G, and fix ϕ ∈ D and x ∈ X . The Taylor formula gives us

ϕ(Φ EM ∆t (x, G)) = ϕ(x) + h ∆t (x, G) • ∇ϕ(x) + 1 2 h ∆t (x, G) T • ∇ 2 ϕ(x) • h ∆t (x, G) + 1 6 l∈N N d |l|=3 (h ∆t (x, G)) l 1 .. (h ∆t (x, G)) l N d ∂ 3 ϕ ∂x l (x) + 1 24 l∈N N d |l|=4 (h ∆t (x, G)) l 1 .. (h ∆t (x, G)) l N d ∂ 4 ϕ ∂x l (x) + 1 120 l∈N N d |l|=5 (h ∆t (x, G)) l 1 .. (h ∆t (x, G)) l N d 1 0 (1 -t) 4 ∂ 5 ϕ ∂x l (x + th ∆t (x, G)).
Because h ∆t (x, G) = O(∆t 1/2 ), the integral term of the above expansion is order at least 5 2 in ∆t. In the sequel, we use the notation

(x 1 ⊗ . . . ⊗ x n ) ⊙ D n ϕ = N d i 1 =1 . . . N d i n =1 ∂ n ϕ ∂x i 1 . . . ∂x i n (x) n k=1 (x k ) i k ,
where n ∈ N and (x 1 , .., x n ) ∈ X n is a family of elements of X . The above notation denotes in fact all the terms of the sum of order n of the expansion of ϕ(Φ EM ∆t (x, G)) containing occurrences of (x 1 , .., x n ). Expanding the above expression in powers of ∆t therefore yields

ϕ Φ EM ∆t (x, G) = ϕ(x) + √ ∆tb(x)G + ∆t f (x) • ∇ϕ + 1 2 b(x) T • ∇ 2 ϕ(x) • b(x)G 2 + ∆t 3 2 b(x) T • ∇ 2 ϕ(x) • f (x)G + 1 2 (b(x) ⊗ b(x) ⊗ b(x)) ⊙ D 3 ϕ(x)G 3 , + ∆t 2 1 24 (b(x) ⊗ b(x) ⊗ b(x) ⊗ b(x)) ⊙ D 4 ϕ(x)G 4 + 1 2 (b(x) ⊗ b(x) ⊗ f (x)) ⊙ D 3 ϕ(x)G 2 + 1 2 f (x) T • ∇ 2 ϕ(x) • f (x) + ∆t
where all the terms of power greater or equal than 5/2 have been regrouped in the remainder term ∆t 5/2 R∆t,ϕ .

We denote by rϕ (x, G) the coefficient of the order 3/2 term in ∆t of the above expression, and by k ϕ (x, G) the coefficient of the order 2 term. This gives

ϕ Φ EM ∆t (x, G) = ϕ(x) + √ ∆tb(x)G + ∆t f (x) • ∇ϕ + 1 2 b(x) T • ∇ 2 ϕ(x) • b(x)G 2 + ∆t 3 2 r∆t,ϕ (x, G) + ∆t 2 k ϕ (x, G) + √ ∆t R∆t,ϕ (x, G) .
We notice that both terms corresponding to non-integer powers of ∆t have null expecta-

tions, i.e E G [b(x)G] = E G [r ∆t,ϕ ] = 0.
In addition, the expectation of the first order term corresponds to Lϕ(x), where L is the generator of the SDE. Therefore,

P EM ∆t ϕ(x) = E G ϕ(Φ EM ∆t (x, G)) = ϕ(x) + ∆tLϕ(x) + ∆t 2 E G k ϕ (x, G) + √ ∆t R∆t,ϕ (x, G) r ∆t,ϕ (x)
.

It remains to ensure that the remainder term r ∆t,ϕ is uniformly bounded when ∆t goes to zero for some weighted space L ∞ W m with m ∈ N. Because (f , b, ϕ) ∈ D 3 and because k ∆t,ϕ (x, G) and R∆t,ϕ (x, G) involve only positive powers of ∆t, we have for a given ∆t * > 0 and m ∈ N that

sup 0<∆t<∆t * E G [k ϕ (•, G)] L ∞ W m + E G [ R∆t,ϕ (•, G)] L ∞ W m < +∞, which implies that E G √ ∆t R∆t,ϕ (•, G)
goes to 0 when ∆t goes to zero. Therefore,

sup 0<∆t<∆t * r ∆t,ϕ L ∞ W m < +∞,
an P EM ∆t satisfies (1.34) at the order ω = 1, which is the desired result.

Errors on the invariant measure

Weak errors quantify the accuracy of a numerical scheme on the average behaviors of the solutions, and strong errors quantify the accuracy on their trajectories. When doing Molecular Dynamics, one can also be interested in the accuracy on the estimations of the average values of physical observables as defined in Definition 1.2.1. These estimations require stronger assumptions both on the dynamics and on the numerical scheme than weak or strong error, as made precise in this section.

Consider the dynamics (1.20), and assume that it is ergodic for a macrostate π. Consider a scheme whose evolution operator is denoted by P ∆t , and an observable ϕ ∈ D. Using the same notation than in Section 1.3.1, the estimation ϕ N it ,∆t of the average value

E π [ϕ] of ϕ reads ϕ N it ,∆t = 1 N it N it n=1 ϕ(x n ), where (x n ) n∈[|1,N it |]
is a sample of the approximated solution generated by the numerical scheme. We are interested in how well P ∆t estimates the average value of ϕ, i.e the quality of the approximation of E π [ϕ] by ϕ N it ,∆t . The first task is proving the existence of an invariant probability measure for the numerical scheme, measure that we denote by π ∆t , using for instance results such as Theorem 1.3.2. In the sequel, we assume that P ∆t admits a unique invariant probability measure π ∆t and that the scheme is ergodic with respect to this measure. In addition, we assume that it is weakly consistent, of weak order 1. Denoting by

K ∆t,ϕ = E π ∆t [ϕ] -E π [ϕ],
we have that

ϕ N it ,∆t = 1 N it N it n=1 ϕ(x n ) = E π [ϕ] + E π ∆t [ϕ] -E π [ϕ] + ϕ N it ,∆t -E π ∆t [ϕ] R∆t,ϕ , = E π [ϕ] + K ∆t,ϕ + R∆t,ϕ .
The ergodicity of the numerical scheme tells us that the rightmost term of the right-hand side of the above equality, denoted by R∆t,ϕ , goes to 0 when

N it → +∞, since 1 N it N it n=1 ϕ(x n ) ------→ N it →∞ X ϕdπ ∆t = E π ∆t [ϕ], π ∆t -a.s.
In order to obtain error estimates on the estimation of E π [ϕ] by ϕ N it ,∆t , we must quantify the rate of convergence of R∆t,ϕ towards 0.

Under some conditions on P ∆t and on ϕ, the Central Limit theorem holds for the Markov chain generated by P ∆t (see Theorem 17.4.4 in [START_REF] Meyn | Markov Chains and Stochastic Stability[END_REF] or [START_REF] Lelièvre | Partial differential equations and stochastic methods in molecular dynamics[END_REF] and references therein) and we have

N it ϕ N it ,∆t -E π ∆t [ϕ] law -------→ N it →+∞ N (0, σ 2 ∆t,ϕ ), with σ 2 ϕ,∆t = E π ∆t (Π ∆t ϕ) 2 + 2 +∞ n=1 E π ∆t (Π ∆t ϕ(x n )) Π ∆t ϕ(x 0 ) .
In the above expression, we have denoted by Π ∆t the projection operator over the set of functions of null expectations, i.e

∀ϕ ∈ L 1 (π ∆t ), Π ∆t ϕ = ϕ -E π ∆t [ϕ].
It can be shown (see Section 3 in [START_REF] Lelièvre | Partial differential equations and stochastic methods in molecular dynamics[END_REF]) that, under some conditions on the decrease of e tL in the appropriate norm, the weak consistency of the scheme implies that √ ∆tσ ∆t,ϕ → σ ϕ when ∆t goes to 0, with σ 2 ϕ being the variance of the continuous dynamic expressed as

σ 2 ϕ = 2 +∞ 0 E π e tL Πϕ (Πϕ) dt, with Π being the equivalent of Π ∆t for π, i.e ∀ϕ ∈ L 1 (π), Πϕ = ϕ - X ϕdπ.
Therefore, asymptotically (i.e taking ∆t small enough and N it large enough), we have

ϕ N it ,∆t ≃ E π ∆t [ϕ] + ξ ∆t,N it ,ϕ √ N it ∆t ,
where ξ ∆t,N it ,ϕ is a random variable whose distribution converges in the limit ∆t → 0 and N it ∆t → +∞ to a Gaussian random variable of variance σ 2 ϕ . We therefore can rewrite ϕ N it ,∆t as

ϕ N it ,∆t ≃ E π [ϕ] + K ∆t,ϕ + ξ ∆t,N it ,ϕ √ N it ∆t .
Let us consider the right-hand side of the above expression with more attention.

• The first term E π [ϕ] is the exact physical value of the observable considered, i.e the value that we want to estimate.

• The second term K ϕ,∆t is called the perfect sampling bias, i.e the error arising from the fact that the invariant probability measure π ∆t is different from the unique invariant probability measure π of the SDE.

• The third term is called the finite sampling bias, or statistical error. It is the error arising from the use of a finite sample size N it , i.e the difference between the estimation ϕ N it and the target expectation E π ∆t [ϕ]. Note that it is inversely proportional to the square root of the physical time T = N it ∆t.

Therefore, if we take a physical time T large enough so that

σ ϕ √ T
≪ K ∆t,ϕ , then the statistical error can be neglected and the accuracy of the scheme on the estimation of the average value of ϕ is given by perfect sampling bias K ∆t,ϕ . From a practical point of view, we show that we have estimates of the form K ∆t,ϕ ≤ K ϕ ∆t ω for a given ω > 0, with K ϕ > 0 independent of ∆t. Therefore, in order for the statistical error to be negligible, it is required that T ≫ σ 2 ϕ ∆t 2ω . Depending on the variance σ 2 ϕ of the continuous dynamics, the time step ∆t or the order ω of the perfect sampling bias, this might mean very (sometimes prohibitively) long simulation times.

Remark 1.3.3. If π ∆t = π, i.e if the scheme's invariant probability measure is the unique invariant probability measure of the SDE, as for the Metropolis-Hastings algorithms [MRR + 53, Has70], the perfect sampling bias is null. The only difference between the estimation and the target quantity lies therefore in the statistical error.

Similarly to weak and strong error estimates, we can define estimates of the error on the invariant probability measure to quantify K ∆t,ϕ . Definition 1.3.7. [Error on the invariant measure] Consider a dynamics given by (1.20) and a numerical scheme P ∆t , and assume that they both admit ergodic probability measures denoted respectively by π and π ∆t . We say that P ∆t has errors on the invariant measure of order ω ∈ N when

∀ϕ ∈ D, E π ∆t [ϕ] -E π [ϕ] ≤ K ϕ ∆t ω , (1.35)
with K ϕ > 0 independent of ∆t. The integer ω is called the order of the error on the invariant measure.

It is important to note that error on the invariant measure estimates require the existence of an invariant probability measure π ∆t , whereas weak or strong error estimates do not. Errors on the invariant measure are closely linked with weak errors. Under ergodicity assumptions on the dynamics and on the numerical scheme under consideration, Talay and al. in [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF] have proved for instance that a weak order ω > 0 implies errors on the invariant measure of the same order for both the Euler-Maruyama and Mil'shtein scheme [START_REF] Milstein | Stochastic Numerics for Mathematical Physics[END_REF].

In some cases, the order on the invariant measure can be higher than the weak order. Indeed, assume that expansions as in the proof of Proposition 1.3.4 yield

∀ϕ ∈ D, P ∆t ϕ = ϕ + L∆t + ω i=2 ∆t i A i ϕ + ∆t ω+1 A ω+1 ϕ + ∆t ω+2 r ∆t,ϕ , (1.36) with (A i ) i=2,.
.,ω+1 being function operators and r ∆t,ϕ a remainder term uniformly bounded for small ∆t in some weighted space, i.e

∃(m, C ϕ , ∆t * ) ∈ N × (R * + ) 2 , sup 0<∆t≤∆t * r ∆t,ϕ L ∞ W m < C ϕ .
If we have A i = 1 i! L i for i = 2, .., ω, then the scheme is of weak order ω. But sometimes, one might have A i = a i L i with a i 1 i! . In such cases, the scheme is not of weak order ω, but we can still have errors on the invariant measures of order ω. In addition, under some conditions detailed below, it is possible to have the dominant term of the error on the invariant measure.

In the sequel, we define D 0 as the subset of D containing the functions of null average with respect to. π, i.e

D 0 = ϕ ∈ D X ϕdπ = 0 .
(1.37)

In opposition to D, we can define the notion of inversibility of differential operators like

L in D 0 . Indeed, considering any ϕ ∈ D such that E π [ϕ] < +∞, we have Lϕ = L(ϕ-E π [ϕ]
), which forbids L to be invertible on D, except by imposing a condition on the value of the expectations. We also define as A † the adjoint in L 2 (π) of any operator A, i.e

∀(ϕ, ψ) ∈ D 2 ∞ , X (Aϕ) ψdπ = X ϕ A † ψ dπ.
The following theorem allows to obtain estimates on the order of the error on the invariant measure and to make the leading error term precise.

Theorem 1.3.5. Consider the SDE given by (1.20), whose generator is denoted by L, and a numerical scheme whose evolution operator is denoted by P ∆t . Assume that L and P ∆t satisfy the following conditions:

1. The dynamics is ergodic for a probability measure π that integrates all the functions W m with m ≥ 0, i.e.

∀m ∈ N, X W m dπ < +∞.
In particular, D is dense in L 2 (π).

2. Both L and L † are invertible on D 0 , and L -1 and

(L † ) -1 leave D 0 invariant, i.e ∀ϕ ∈ D 0 ,        ∃!ψ 1 ∈ D 0 , Lψ 1 = ϕ, ∃!ψ 2 ∈ D 0 , L † ψ 2 = ϕ.
3. The numerical scheme is ergodic for a probability measure π ∆t that integrates all the functions W m with m ≥ 0, i.e.

∀m ∈ N, X W m dπ ∆t < +∞.
In particular, D is dense in L 2 (π ∆t ).

4. The evolution operator P ∆t can be expanded as in (1.36).

The following relation holds for all

i = 2, . . . , ω: ∀i ∈ [|2, ω|], ∀ϕ ∈ D, X A i ϕdπ = 0, (1.38)
i.e the operators A i for i = 2, .., ω preserve π (but A ω+1 might not). In addition, A i leaves D invariant for i = 2, .., ω.

There exists

g ω+1 ∈ D 0 such that ∀ϕ ∈ D, X A ω+1 ϕdπ = X g ω+1 ϕdπ.
(1.39)

In addition, ΠA ω+1 Π leaves D 0 invariant.

Then there exists f ω+1 ∈ D 0 and ∆t * > 0 such that for any 0

≤ ∆t ≤ ∆t * , X ϕdπ ∆t = X ϕdπ + ∆t ω X f ω+1 ϕdπ + ∆t ω+1 R ∆t,ϕ , (1.40)
with R ∆t,ϕ uniformly bounded for small ∆t, i.e

∃(C ϕ , ∆t * ) ∈ (R * + ) 2 , sup 0<∆t<∆t * R ∆t,ϕ < C ϕ .
In particular, the numerical scheme P ∆t is of order ω on the invariant measure.

Let us comment on the above conditions.

• Assumption 1 is obtained with results such as Theorem 1.2.1, (see Section 1.2.2). In addition, the integrability of the scale functions W m with respect to π is satisfied by many statistical ensembles. This is the case for the canonical ensemble µ β for instance. We see in the proof of Theorem 1.3.5 that it allows to have the remainder term of (1.40) uniformly bounded in L ∞ (X ).

• Assumption 2 is not trivial to obtain in the general case. It is satisfied if L is elliptic and the position space is bounded. However, it becomes difficult to prove when one of the two condition is missing. We discuss later on in Section 1.4.1 the case of two prototypical example dynamics.

• Assumption 3 is obtained with results such as Proposition 1.3.1 (see Section 1.3.1).

• Assumption 4 is obtained by performing Taylor expansions as in the proof of Proposition 1.3.4.

• Assumption 5 is trivial if the operators A i are proportional to the powers of the generator L, i.e A i ∝ L i . In addition, the operators A i typically are differential operators whose coefficients are derived from f and b of (1.20). They are thus all elements of D if f and b are also elements of D, which is assumed (see Condition 1.2.2).

• Assumption 6 is easily obtained. Indeed, we see that g ω+1 = A † ω+1 1, where 1 is the constant function equal to 1 everywhere. Therefore, the only condition for (1.39) to hold is that A ω+1 stabilizes D 0 , i.e ϕ ∈ D 0 implies A ω+1 ϕ ∈ D 0 . We can usually check by performing integration by parts that g ω+1 ∈ D. In addition, if we consider (1.39) with ϕ = 1, we have that X g ω+1 dπ = 0 if A ω+1 1 = 0, which is the case if A ω+1 is a differential operator.

Further details about such assumptions are given in [START_REF] Lelièvre | Partial differential equations and stochastic methods in molecular dynamics[END_REF].

Remark 1.3.4. Note that in many cases, a transition operator P ∆t satisfying (1.34) at the order ω also satisfies conditions 4, 5 and 6 of Theorem 1.3.5. Indeed, we see that expansions such as (1.28) satisfied at the order ω also satisfy (1.36) at the order ω for some operator A ω+1 , defined as the first term of the expansion of P ∆t in powers of ∆t that is not equal to L k k! for some k ∈ N. In addition, for many schemes, A ω+1 is a differential operator stabilizing D 0 , thus ensuring condition 6 if the functions f and b of the dynamics are smooth enough. The lower order terms being proportional to powers of L, condition 5 is satisfied. Therefore, provided that conditions 1, 2 and 3 are satisfied, the weak order ω often imply that (1.40) holds.

Proof. First notice that we can restrain the proof to functions in D 0 without loss of generality. Indeed, assuming that (1.40) holds for any ψ ∈ D 0 , and taking ϕ ∈ D, we have

X ϕ - X ϕdπ dπ ∆t = ∆t ω X f ω+1 ϕ - X ϕdπ dπ + ∆t ω+1 R ∆t,ϕ .
Because f ω+1 ∈ D 0 , we have X f ω+1 X ϕdπ dπ = 0 and therefore

X ϕdπ ∆t = X ϕdπ + ∆t ω X f ω+1 ϕdπ + ∆t ω+1 R ∆t,ϕ .
Thus, (1.40) also holds for functions in D.

We start by proving (1.40) for ϕ ∈ Ran(P ∆t -Id). Consider ψ ∈ D. We first determine the expression of f ω+1 . For f ω+1 given, we have by (1.38) and (1.39),

X P ∆t -Id ∆t ψ (1 + ∆t ω f ω+1 ) dπ = ∆t ω X (A ω+1 ψ + (Lψ)f ω+1 ) dπ + ∆t ω+1 X r ∆t,ψ + P ∆t -Id -∆tLψ ∆t 2 f ω+1 dπ.
The first term of the right-hand side of the above equality can be rewritten as

X (A ω+1 ψ + (Lψ)f ω+1 ) dπ = X g ω+1 + L † f ω+1 ψdπ.
Thus, in order for it to vanish, we define f ω+1 as the unique solution of -L † f ω+1 = g ω+1 .

By assumption, this means that

f ω+1 = -L † -1 g ω+1 = -(L † ) -1 A † ω+1 1. With such a choice, we can write ∀ψ ∈ D, X P ∆t -Id ∆t ψ (1 + ∆t ω f ω+1 ) dπ = ∆t ω+1 X r ∆t,ψ + P ∆t -Id -∆tA 1 ψ ∆t 2 f ω+1 dπ s ∆t,ϕ
.

(1.41) Because the integrand of the right-hand side of (1.41) is bounded in some weighted space, we have by Assumption 1 that s ∆t,ϕ is uniformly bounded for ∆t small enough, i.e

∃(C ϕ , ∆t * ) ∈ (R * + ) 2 , sup 0<∆t<∆t * s ∆t,ϕ < C ϕ .
Equation (1.41) implies that (1.40) holds for any ϕ ∈ Ran(P ∆t -Id).

We now want to extend this result to ϕ ∈ D 0 , and for that we construct a pseudoinverse of P ∆t -Id. It can be constructed by truncating the formal series expression of the inverse of the operator A + ∆tB, with A = L and B = ΠA 2 Π + . . . + ∆t ω-1 ΠA ω+1 Π. The formal series writes

(A + ∆tB) -1 = n (-1) n BA -1 n A -1 .
The pseudo-inverse Q ∆t is then defined by removing from the formal series all powers of ∆t greater than ω. Because L and

ΠA i Π leave D 0 invariant for i = 2, . . . , ω + 1 (due to (1.38)), Q ∆t also leaves D 0 invariant. In addition, Q ∆t satisfies ∀ψ ∈ D 0 , Π P ∆t -Id ∆t ΠQ ∆t ϕ = ϕ + ∆t ω+1 r∆t,ϕ , (1.42)
with r∆t,ϕ uniformly bounded for ∆t small enough. More details about the construction of Q ∆t are given in [START_REF] Lelièvre | Partial differential equations and stochastic methods in molecular dynamics[END_REF].

Now that we have our pseudo-inverse operator Q ∆t , we prove (1.34) for ϕ ∈ D 0 . We first notice that

(P ∆t -Id) Πϕ = (P ∆t -Id) ϕ - X ϕdπ = (P ∆t -Id) ϕ, and therefore, if ϕ ∈ D 0 , X (P ∆t -Id) Πϕdπ = X (P ∆t -Id) ϕdπ = X P ∆t ϕ.
Using the two above equality, we can write

X Π (P ∆t -Id) Πϕ (1 + ∆t ω f ω+1 ) dπ = X (P ∆t -Id) Πϕ (1 + ∆t ω f ω+1 ) dπ - X X (P ∆t -Id) Πϕdπ (1 + ∆t ω f ω+1 ) dπ, = X (P ∆t -Id) Πϕ (1 + ∆t ω f ω+1 ) dπ - X (P ∆t -Id) Πϕdπ X (1 + ∆t ω f ω+1 ) dπ, = X (P ∆t -Id) ϕ (1 + ∆t ω f ω+1 ) dπ - X P ∆t ϕdπ.
In the last line of the above equation, the second term of the right-hand side is obtained by saying that X (1 + ∆t ω f ω+1 ) dπ = 1 since f ω+1 ∈ D 0 and π is a probability measure (so

X dπ = 1)
. The first term of the above right-hand side is proportional to ∆t ω+2 according to (1.41). Thus,

∀ϕ ∈ D 0 , X Π P ∆t -Id ∆t Πϕ (1 + ∆t ω f ω+1 ) dπ = - 1 ∆t X P ∆t ϕdπ + ∆t ω+1 s ∆t,ϕ ,
where s ∆t,ϕ ∈ D is given by (1.41). We have that ϕ is in D 0 , therefore Q ∆t ϕ is also in D 0 , and the above equality also holds for Q ∆t ϕ. Thus, by using Q ∆t ϕ instead of ϕ in the above equality and (1.42), we have for any ϕ ∈ D 0 ,

X ϕ (1 + ∆t ω f ω+1 ) dπ = - 1 ∆t X P ∆t Q ∆t ϕdπ + ∆t ω+1 X s ∆t,ϕ + r∆t,ϕ (1 + ∆t ω f ω+1 ) dπ s∆t,ϕ ,
(1.43) with s∆t,ϕ uniformly bounded for small ∆t.

In order to compute the last non-remainder term in (1.43), we first prove the following relation:

∀ϕ ∈ D 0 , X Π P ∆t -Id ∆t Πϕdπ ∆t = - 1 ∆t X P ∆t ϕdπ, (1.44)
Equation (1.44) is obtained after a small computation using

∀ψ ∈ D, X dπ ∆t = 1, X (P ∆t -Id) ψdπ ∆t = 0.
The first equality holds because π ∆t is a probability measure, and the second because π ∆t is an invariant probability measure of P ∆t . Therefore, by replacing ϕ with Q ∆t ϕ, (1.44)

gives ∀ϕ ∈ D 0 , X ϕdπ ∆t + X ∆t ω+1 r∆t,ϕ dπ ∆t = - 1 ∆t X P ∆t Q ∆t ϕdπ.
Injecting the above equation into (1.43) and grouping the remainder terms gives us the result for functions ϕ ∈ D 0 .

Splitting techniques for Stochastic Differential Equations

Splitting techniques consist in decomposing the dynamics of interest into several elementary subdynamics, which are sequentially integrated. The resulting accuracy depends both on the accuracy of the subdynamics integration and the order of the splitting formula. Because splitting formulas with arbitrary order exist (see Sections III and V of [START_REF] Hairer | Geometrical Numerical Integration -Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF] for examples), splitting techniques are used for devising higher-order schemes when the subdynamics are simple, or simpler, to discretize. However, for reaching weak orders larger than two, integration with negative time steps must be resorted to, which may alleviate the stability of the scheme. In practice, splitting schemes are a systematic way of devising weak order two schemes, which is enough for many systems of interest.

Let us consider the SDE represented in its general form in (1.20), whose infinitesimal generator we denote by L. Let us assume that we can decompose L as the sum of two subgenerators, i.e L = A + B.

Let P A ∆t and P B ∆t be transition operators associated to numerical schemes discretizing the SDEs associated with A and B respectively, and ∆t the time step. We detail in the sequel how to construct a transition operator P ∆t from P A ∆t and P B ∆t using splitting techniques, namely the Trotter splitting and Strang splitting.

Trotter splitting formula

The Trotter formula approximates an evolution operator by composing successively the evolution operators of each of the reduced parts. It was introduced by H. F. Trotter in [START_REF] Trotter | On the Product of Semi-Groups of Operators[END_REF], and writes e t(A+B) = lim ∆t→0 e ∆tA e ∆tB ⌊ t ∆t ⌋ .

(1.45)

A splitting scheme using the Trotter decomposition given by (1.45) therefore writes

P ∆t = P A ∆t P B ∆t .
(1.46) Lemma 1.3.6. Assume that P A ∆t , P B ∆t e ∆tA and e ∆tB stabilize D. Additionally assume that both associated schemes are of weak order one in the sense that for all ϕ ∈ D,

       P A ∆t ϕ = e ∆tA ϕ + ∆t 2 r A ∆t,ϕ , P B ∆t ϕ = e ∆tB ϕ + ∆t 2 r B ∆t,ϕ ,
where

∃(a, b, C ϕ , ∆t * ) ∈ N 2 × (R * + ) 2 , sup 0<∆t≤∆t * r A ∆t,ϕ L ∞ W a + r B ∆t,ϕ L ∞ W b < C ϕ .
Then the scheme P ∆t defined by (1.46) satisfies for all ϕ ∈ D:

P ∆t ϕ = e ∆tL ϕ + ∆t 2 R ∆t,ϕ , with R ∆t,ϕ L ∞ W n
< K ϕ for some K ϕ > 0, n ∈ N and any ∆t small enough. In particular, P ∆t is of weak order one.

Proof. The Baker-Campbell-Hausdorff (BCH) formula (see [START_REF] Baker | Alternants and continuous groups[END_REF][START_REF] Campbell | On a law of combination of operators (second paper)[END_REF][START_REF] Hausdorff | Die symbolische exponentialformel in der gruppentheorie[END_REF] or [START_REF] Hairer | Geometrical Numerical Integration -Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF]) allows to compute a product of exponential operators as an exponential of operator. At first order, it formally writes for any operators A and B: ≤ J ϕ for some J ϕ ≥ 0, k ∈ N and for any ∆t small enough. For a proof of this result, we refer to [START_REF] Leimkuhler | The Computation of Averages from Equilibrium and Non-Equilibrium Langevin Molecular Dynamics[END_REF] or [START_REF] Hairer | Geometrical Numerical Integration -Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF].

e A e B = e A+B+ 1 2 [A,B]+ 1 12 ([A,[A,B]]-[B,[A,B]])+... , ( 1
By assumption,       

P A ∆t ϕ = e ∆tA ϕ + ∆t 2 r A ∆t,ϕ , P B ∆t ϕ = e ∆tB ϕ + ∆t 2 r B ∆t,ϕ , where ∃(a, b, C ϕ , ∆t * ) ∈ N 2 × (R * + ) 2 , sup 0<∆t≤∆t * r A ∆t,ϕ L ∞ W a + r B ∆t,ϕ L ∞ W b < C ϕ .
Because P A ∆t , P B ∆t , e ∆tA and e ∆tB stabilize D, we have that r A ∆t,ψ ∈ D and r B ∆t,ψ ∈ D for any ψ ∈ D. Therefore, we have

P ∆t ϕ = P A ∆t P B ∆t ϕ, = e ∆tA e ∆tB ϕ + ∆t 2 e ∆tA r B ∆t,ϕ + r A ∆t,e ∆tB ϕ + ∆t 4 r A ∆t,r B ∆t ϕ
.

(1.49)

Injecting (1.48) into the last line of (1.49), we obtain

P ∆t ϕ = ϕ + e ∆tL ϕ + ∆t 2 R ∆t,ϕ , (1.50)
where R ∆t,ϕ = r ∆t,ϕ + e ∆tA r B ∆t,ϕ + r A ∆t,e ∆tB ϕ + ∆t 2 r A ∆t,r B ∆t ϕ .

It remains to show that R ∆t,ϕ is uniformly bounded in some L ∞ W m space for some m ∈ N and any ∆t small enough. We have by assumption and by (1.48) that, for some ∆t * small enough and m ∈ N,

sup 0<∆t<∆t *       r ∆t,ϕ L ∞ W m + ∆t 2 r A ∆t,r B ∆t ϕ L ∞ W m + r A ∆t,e ∆tB ϕ L ∞ W m       < +∞.
Because e ∆tA stabilize D and that r B ∆t,ϕ ∈ D, we have that e ∆tA r B ∆t,ϕ ∈ D. In addition, we can write for any ψ ∈ D,

e ∆tA ψ = ψ + ∆t r∆t,ψ , ( 1

Strang splitting formula

The Strang formula, also called symmetric Trotter formula, was introduced by G. Strang in [START_REF] Strang | On the Construction and Comparison of Difference Schemes[END_REF] for linear hyperbolic PDEs. It consists in symmetrizing the Trotter formula (1.45) and performing two iterations with half of the time step around one central iteration with full time step. A splitting scheme using the Strang decomposition writes

P ∆t = P A ∆t 2 P B ∆t P A ∆t 2
.

(1.52) Lemma 1.3.7. Assume that P A ∆t , P B ∆t e ∆tA and e ∆tB stabilize D. Additionally assume that both associated schemes are of weak order two in the sense that for all ϕ ∈ D,

       P A ∆t ϕ = e ∆tA ϕ + ∆t 3 r A ∆t,ϕ , P B ∆t ϕ = e ∆tB ϕ + ∆t 3 r B ∆t,ϕ ,
where

∃(a, b, C ϕ , ∆t * ) ∈ N 2 × (R * + ) 2 , sup 0<∆t≤∆t * r A ∆t,ϕ L ∞ W a + r B ∆t,ϕ L ∞ W b < C ϕ .
Then the scheme P ∆t defined by (1.52) satisfies for all ϕ ∈ D:

P ∆t ϕ = e ∆tL ϕ + ∆t 3 R ∆t,ϕ , with R ∆t,ϕ L ∞ W n
< K ϕ for some K ϕ > 0, n ∈ N and any ∆t small enough. In particular, P ∆t is of weak order two.

Proof. The proof of the weak order two of a Strang splitting scheme uses the symmetric BCH formula formally written as where once again the neglected terms are linear combinations of the commutators of A and B of order greater than 3.

Similarly to the proof of Lemma 1.3.6, (1.53) allows the following relations to be satisfied for all ϕ ∈ D: e ∆tA/2 e ∆tB e ∆tA/2 ϕ = e ∆tL ϕ + ∆t 3 r ∆t,ϕ , with r ∆t,ϕ L ∞ W k ≤ J ϕ for some J ϕ ≥ 0, k ∈ N and for any ∆t small enough.

Fix ϕ ∈ D. By assumption,

       P A ∆t ϕ = e ∆tA ϕ + ∆t 3 r A ∆t,ϕ , P B ∆t ϕ = e ∆tB ϕ + ∆t 3 r B ∆t,ϕ ,
where

∃(a, b, C ϕ , ∆t * ) ∈ N 2 × (R * + ) 2 , sup 0<∆t≤∆t * r A ∆t,ϕ L ∞ W a + r B ∆t,ϕ L ∞ W b < C ϕ .
A small computation allows to write P ∆t ϕ as where, in passing from the second to the third line in the above computation, we have factorized the other terms than e ∆t A 2 e ∆tB e ∆t A 2 ϕ into a power three term in ∆t and a prefactor denoted by r A,B ∆t,ϕ . Note that r A,B ∆t,ϕ ∈ D because P A ∆t , P B ∆t , e ∆tA and e ∆tB stabilize D. Denoting by R ∆t,ϕ the remainder term of the above equation, we have, by a similar argument than for Lemma 1.3.6, that R ∆t,ϕ remains uniformly bounded when ∆t goes to zero for some weighted space L ∞ W n with n ∈ N.

P ∆t ϕ = P A ∆t 2 P B ∆t P A ∆t 2 ϕ, = e ∆t A

Geometric splitting

We have proved that, provided that the discretizations of the subdynamics are of weak order high enough, Trotter splitting schemes are of weak order one and Strang splitting schemes of weak order two. However, it is possible to have an arbitrary high order on the invariant measure with Trotter splitting schemes, provided that the error on the invariant measure of the subschemes are high enough. Indeed, under conditions similar to those of Theorem 1.3.5 applying to both subdynamics, the resulting Trotter splitting scheme satisfies (1.40).

Lemma 1.3.8. Assume that P A ∆t and P B ∆t stabilize D and that both associated schemes are of weak order w in the sense that there exists two differential operators A ω+1 and B ω+1 such that for all ϕ ∈ D,       

P A ∆t ϕ = e ∆tA ϕ + ∆t ω+1 A ω+1 ϕ + ∆t ω+2 r A ∆t,ϕ , P B ∆t ϕ = e ∆tB ϕ + ∆t ω+1 B ω+1 ϕ + ∆t ω+2 r B ∆t,ϕ , where ∃(a, b, C ϕ , ∆t * ) ∈ N 2 × (R * + ) 2 , sup 0<∆t≤∆t * r A ∆t,ϕ L ∞ W a + r B ∆t,ϕ L ∞ W b < C ϕ .
Also assume that there exists

(f A ω+1 , f B ω+1 ) ∈ D 2 0 such that for all ϕ ∈ D,              X A ω+1 ϕdπ = X f A ω+1 ϕdπ, X B ω+1 ϕdπ = X f B ω+1 ϕdπ.
Additionally assume that both L and L † are invertible on D 0 , that

∀m ∈ N, X W m dπ < +∞,
and that both subdynamics preserve the invariant measure π, i.e

∀ϕ ∈ D, X Aϕdπ = X Bϕdπ = 0.
Note that this means in particular that e ∆tA and e ∆tB stabilize D.

Then the scheme P ∆t defined by (1.46) satisfies equation (1.40), i.e there exists f ω+1 ∈ D 0 and ∆t * > 0 such that for any 0

≤ ∆t ≤ ∆t * , ∀ϕ ∈ D, X ϕdπ ∆t = X ϕdπ + ∆t ω X f ω+1 ϕdπ + ∆t ω+1 R ∆t,ϕ ,
with R ∆t,ϕ uniformly bounded for small ∆t, i.e R ∆t,ϕ < C ϕ for some C ϕ > 0 and for any ∆t small enough. In particular, the numerical scheme P ∆t is of order ω on the invariant measure.

Such schemes are usually called Geometric Splitting Schemes, and, provided both subschemes are accurate enough on the estimations of observables, combine the greater accuracy on observable estimations of Strang splitting schemes and the efficiency and CPU cost of Trotter splitting schemes. However, the above result is only valid for error on the invariant measure and has no counterpart for weak nor strong errors.

Proof. We only prove here the result for D 0 . The extension to D is similar to the one performed in the proof of Theorem 1.3.5.

Fix ϕ ∈ D 0 . Similarly to the proofs of Lemma 1.3.6, we show that the expansions of P A ∆t and P B ∆t allow us to write

X P ∆t ϕdπ = X P A ∆t P B ∆t ϕdπ, = X
e ∆tA e ∆tB ϕdπ + ∆t ω+1 X e ∆tA B ω+1 + A ω+1 e ∆tB ϕdπ + ∆t ω+2 X r ∆t,ϕ dπ, with r ∆t,ϕ uniformly bonded on L ∞ W m for some m ∈ N and for all ∆t small enough. Because both subdynamics preserve π, we have

X P ∆t ϕdπ = X ϕdπ + ∆t ω+1 X B ω+1 ϕdπ + X f A ω+1 e ∆tB ϕdπ + ∆t ω+2 X r ∆t,ϕ dπ, = ∆t ω+1 X f B ω+1 ϕdπ + X ϕ e ∆tB † f A ω+1 dπ + ∆t ω+2 X r ∆t,ϕ dπ, = ∆t ω+1 X f B ω+1 + e ∆tB † f A ω+1 f ω+1 ϕdπ + ∆t ω+2 X r ∆t,ϕ dπ.
Therefore, we have

- 1 ∆t X P ∆t ϕdπ = ∆t ω X (-f ω+1 )ϕdπ -∆t ω+1 X r ∆t,ϕ dπ.
In addition, B being a differential operator, we have that B1 = 0. Therefore, Because ϕ ∈ D 0 , we have that Q ∆t ϕ ∈ D 0 . Therefore, by using Q ∆t ϕ instead of ϕ, relation (1.42) turns the above equality into

X (e ∆tB ) † f A ω+1 dπ = X e ∆tB 1 f A ω+1 dπ = X f A ω+1 dπ, Because (f A ω+1 , f B ω+1 ) ∈ D
X ϕdπ ∆t = ∆t ω X (-f ω+1 )ϕdπ -∆t ω+1 X (r ∆t,ϕ + r∆t,ϕ )dπ,
where r∆t,ϕ is the remainder term of (1.42) and is uniformly bounded on L ∞ for all ∆t small enough.

The fact that R ∆t,ϕ = -X (r ∆t,ϕ + r∆t,ϕ )dπ is uniformly bounded for all ∆t small enough is given by the fact that L ∞ W m ⊂ L 1 (π) for all m ∈ N by assumption. Equation (1.40) therefore holds for any ϕ ∈ D 0 , which gives the result by extending to D as discussed in the beginning of the proof.

A prototypical example: the Langevin Dynamics

We now apply in this section the results presented in the previous sections on a prototypical example called the Langevin dynamics. The Langevin dynamics was originally developed by Paul Langevin for the description of microscopic particles in a fluid undergoing many collisions [START_REF] Langevin | Sur la théorie du Mouvement Brownien [On the theory of Brownian Motion[END_REF] (translated in English in [START_REF] Lemons | Paul Langevin's 1908 paper "on the theory of Brownian Motion" ["sur la théorie du Mouvement Brownien[END_REF]). These collisions are modeled by friction forces slowing the particles and a Brownian motion acting on the time evolution of the particles, which allows to reinject some energy into the system.

We choose here to present results for the discretization of the Langevin dynamics since it is now a well-understood dynamics, and that it is structurally very similar the dynamics which has been the object of study of this thesis. In fact, not only the dynamics but also the discretization methods considered in the next chapters are adapted from the numerical methods introduced in this section.

We start in Section 1.4.1 by presenting the Langevin dynamics and some ergodicity results, and give then in Section 1.4.2 some discretization methods using the splitting techniques introduced in Section 1.3.

Description of the dynamics

The Langevin dynamics is defined by adding to the Hamiltonian dynamics (1.4) a friction term proportional to the velocities of the particles and a fluctuation term proportional to a Brownian motion, where the amplitude of both terms is controlled by two scalar parameters γ and σ. It writes

       dq t = M -1 p t dt, dp t = -∇U(q t ) dt -γM -1 p t dt + σdW t .
(1.54)

The term γM -1 p t is a friction term dissipating energy, and the term σdW t is a fluctuation term bringing energy into the system to compensate the loss due to the friction terms, where W t is a d-dimensional Brownian motion. In the sequel, we assume that the potential U is such that ∇U satisfies Condition 1.2.2 so the existence and uniqueness of solutions of (1.54) is obtained, as well as the smoothness of the flow maps.

Remark 1.4.1. Notice however that ∇U satisfying Condition 1.2.2 implies that U should not have any singularity, which is not the case of Lennard-Jones potentials given by (1.2). Indeed, LJ potentials have singularities for any (unphysical) configuration where two particles have the same position. In practice however, particles never collide when one uses reasonable time steps due to the steep increase of the pairwise interaction potential v(r) when r goes to 0, and Lennard-Jones potentials can be reasonably safely used in Langevin simulations.

In this thesis, we mainly use the Langevin dynamics as a sampling device of the canonical measure µ β given by (1.19). In order for µ β to be invariant for the Langevin dynamics, the friction and fluctuation parameters must satisfy the fluctuation/dissipation relation

σ 2 = 2γ β > 0, (1.55)
where β = 1 k B T is the inverse equilibrium energy corresponding to the equilibrium temperature T . We prove in the following theorem that in fact, if γ and σ satisfy (1.55), then the dynamics is ergodic with respect to µ β .

∀(q 0 , p 0 ) ∈ X , ∀ϕ ∈ L 1 (µ β ), lim t→+∞ 1 t t 0 ϕ(q t , p t )dt = X ϕdµ β , a.s.
If γ and σ verify (1.55), the Langevin dynamics can be rewritten as

           dq t = M -1 p t dt, dp t = -∇U(q t ) dt -γM -1 p t dt + 2γ β dW t .
The proof of Theorem 1.4.1 is now a well-known result [Kli87, RB06, LRS10], and we detail it below.

Proof. The proof is organized in two steps:

1. We prove that the canonical probability measure µ β defined by (1.19) is invariant for the Langevin dynamics using (1.55). Note that we do not need σ > 0 for this result.

2. We prove that the generator L of the Langevin process is hypoelliptic. Here, we need σ > 0 for the result.

The ergodicity of the Langevin dynamics with respect to µ β is then given by Theorem 1.2.1.

Let us prove point 1, i.e the invariance of the canonical measure (1.19) for the Langevin dynamics. Let L be the infinitesimal generator of (1.54), which reads

L = M -1 p • ∇ q -∇U(q) • ∇ p -γM -1 p • ∇ p + σ 2 2 ∆ p , (1.56)
The adjoint of L on L 2 (X , R) writes

∀ψ ∈ D, L * ψ = -M -1 p • ∇ q ψ + ∇ q U(q) • ∇ p ψ + γ∇ p • M -1 pψ + σ 2 2 ∆ p ψ. (1.57)
The canonical measure µ β for a separable Hamiltonian satisfies

∇ p µ β (q, p) = -βM -1 pµ β (q, p), ∇ q µ β (q, p) = -β∇U(q)µ β (q, p).
Therefore, when using µ β instead of ϕ in equation (1.57), we have

L * µ β = -β             -M -1 p • ∇U(q) + ∇U(q) • M -1 p =0             µ β + γ∇ p • M -1 pµ β + 1 2 σ 2 ∆ p µ β , = γ∇ p • M -1 pµ β + 1 2 σ 2 ∇ p • -βM -1 pµ β , = γ - β 2 σ 2 ∇ p • M -1 pµ β .
Thus, if parameters γ and σ verify relation (1.55), then L * µ β = 0 and µ β is invariant for (1.54).

Let us now assume (1.55). Point 2, i.e the hypoellipticity of the Langevin dynamics, is obtained by rewriting its generator L as

L = A 0 - N i=1 d k=1 A * i,k A i,k , (1.58)
where N is the number of particle, d the dimension of each vector q i and p i ,

A 0 = M -1 p • ∇ q -∇U(q) • ∇ p A ham -γp T •M -1 ∇ p , and A i,k = γ β ∂ p i,k , (i, k) ∈ [|1, N |]×[|1, d|].
We have that

[A p i,k , A 0 ] = γ β 1 m i ∂ q i,k -γ∂ p i,k .
The fact that

γ β > 0 implies that A 0 , {A i,k } i=0,..,N ,k=1,..,d , {[A i,k , A 0 ]} i=0,.
.,N ,k=1,..,d has full rank. Therefore, L verifies Hörmander condition given by Definition 1.2.7 and is hypoelliptic.

It can be proved that the Langevin dynamics satisfy condition 2 of Theorem 1.3.5. A careful analysis of the proof presented in [START_REF] Talay | Stochastic hamiltonian systems: exponential convergence to the invariant measure, and discretization by the implicit euler scheme[END_REF], provided in [START_REF] Kopec | Weak backward error analysis for langevin process[END_REF], gives the result. Proposition 1.4.2. Both L and L † , denoting respectively the generator of the Langevin dynamics and its adjoint on L 2 (π), are invertible on D 0 .

Limits of the dynamics for infinite friction: Overdamped Langevin

The Overdamped Langevin Dynamics, also called Brownian dynamics, is a simplified version of the Langevin dynamics where only the positions are considered. It has been developed for mesoscale simulations of particles undergoing several collisions where momenta are neglected because it is assumed that they oscillate too fast. Their effect is represented by a Brownian motion acting on the positions. The Overdamped Langevin dynamics writes

dq t = -∇U(q t )dt + 2 β dW t . (1.59) Remark 1.4.2.
The Overdamped Langevin dynamics, as presented in (1.59), is dimensionally inconsistent. Indeed, dq t and ∇U(q t )dt do not have the same physical dimension and one should multiply each term of (1.59) by a parameter of the appropriate dimension equal to 1. We do not however take this into account in the sequel for clarity purposes.

Notice that equation (1.59) is a basic diffusion SDE, admitting the marginal of the canonical measure µ β on the position as an invariant probability measure. This marginal, denoted by ν β , writes ν β (dq) = Z -1 e -βU(q) dq, (1.60)

Splitting of the dynamics

By construction, the Langevin dynamics (1.54) can be decomposed into a Hamiltonian part and a fluctuation/dissipation part. But we have seen in Section 1.1.2 that the Hamiltonian dynamics can further be decomposed into a part related to the positions and a part related to the momenta. Therefore, the Langevin infinitesimal generator can be written as follows

L = A + B + γC.
(1.61)

The operators A and B are respectively the generators of the position and momenta dynamics of the Hamiltonian dynamics (1.4) (see Section 1.1.2 for the expression of their corresponding dynamics and integration). They write

A = M -1 p • ∇ q , B = -∇ q U • ∇ p .
Note that using the notation of the previous section, A ham = A + B. The operator C is the generator of the fluctuation/dissipation dynamics which writes

           dq t = 0, dp t = -M -1 p t dt + 2 β dW t .
(1.62)

As we have seen in Section 1.1.2, the elementary dynamics associated with A and B can be integrated analytically. In addition, we recognize in equation (1.62) the expression of an Ornstein-Uhlenbeck process. The action of the semigroup (e tC ) can therefore be analytically determined, as states the following lemma: Lemma 1.4.5. Consider x t a solution of the following dynamics:

dx t = -αx t dt + σdW t , (1.63)
where W t is a standard Brownian motion and (α, σ) two scalar parameters. Consider an initial condition x 0 ∈ X . Then we have the following equality in law:

x t ∼ e -αt x 0 + σ 1 -e -2αt 2α G, (1.64)
where G ∼ N (0, Id) and x t is the solution of (1.63) starting from x 0 . In particular,

∀ϕ ∈ D, e tC ϕ(x) = E G        ϕ        e -αt x 0 + σ 1 -e -2αt 2α G               .
Proof. An Itô calculus using (1.63) gives d e αt x t = (αx t dt + dx t ) e αt = σe αt dW t .

Therefore, by integrating the above equation from 0 to t, and by multiplying both sides of the equation by e -αt , we obtain the analytical expression of an Ornstein-Uhlenbeck process as

x t = e -αt x 0 + σ t 0 e α(s-t) dW s .

We have that t 0 e α(s-t) dW s ∼ N 0, t 0 e 2α(s-t) ds and t 0 e 2α(s-t) ds = 1-e 2αt 2α . We thus have

x t ∼ e -αt x 0 + σ 1 -e -2αt 2α G,
where G ∼ N (0, Id), which is the desired result.

Given x t , (1.64) gives us the law of x t+∆t and therefore allows us to devise a discretization of (1.63) of infinite weak order as

x n+1 = e -α∆t x n + σ 1 -e -2α∆t 2α G n , (1.65) with G n ∼ N (0, Id).
Transposing the result of Lemma 1.4.5 to the case of the fluctuation/dynamics described by (1.62), the law of the momenta of each particle is given by

p i,t ∼ e - γt m i p i,0 + m i β 1 -e - 2tγ m i G i ,
where G i ≃ N (0, Id) for i = 1, . . . , N . The analytical integration of the fluctuation/dissipation dynamics is therefore given by

p n+1 i = e - γ∆t m i p n i + m i β 1 -e - 2γ∆t m i G n .
(1.66) Thus, all three subdynamics of (1.61) can be integrated analytically if taken separately. The weak order of any Langevin splitting scheme is therefore given by the order of its splitting formula.

Examples of Langevin splitting schemes

Let P A ∆t , P B ∆t and P γC ∆t be the evolution operators corresponding to the analytical integration of respectively A, B and γC, and let P ∆t be a splitting scheme applied to L. The splitting schemes introduced in Section 1.3.3 adapted to the Langevin dynamics (1.54) are presented below. .

(1.68)

One can define similar schemes by permuting the order of integration of each subdynamics in the above examples. Geometric splitting schemes, called Geometric Langevin Algorithms (GLA), can also be devised for the Langevin dynamics [START_REF] Bou-Rabee | Long-Run Accuracy of Variational Integrators in the Stochastic Context[END_REF]. They consist in splitting the integration between the conservative part generated by A + B and the fluctuation/dissipation part generated by γC. The conservative part is integrated using a Velocity-Verlet scheme and the fluctuation/dissipation part is integrated by the analytical integration presented above in order to obtain two schemes of weak order two (see Proposition 2.2.3 in Section 2.2.2 for a proof of the weak order two of the Velocity-Verlet integration). By considering all the permutations in the order of integration such that each part can be integrated at the order two, we obtain the following four schemes: (1.69)

P γC,A,
For the example, we write below the GLA scheme corresponding to the above lower right evolution operator:

P B,A,B,γC ∆t :                              p n+1/2 i = p n i - ∆t 2 ∇ q i U(q n ), q n+1 i = q n i + ∆tM -1 p n+1/2 i , pn+1 i = p n+1/2 i - ∆t 2 ∇ q i U(q n+1
),

p n+1 i = pn+1 i + e - γ∆t m i pn+1 i + m i β 1 -e - 2γ∆t m i G n .
(1.70)

As seen in Section 1.3.3, Lemma 1.3.6 and 1.3.7 tell us that all the Trotter and GLA schemes are of weak order one, and that the Strang splitting schemes are of weak order two. Let us now turn to errors on the invariant measure and the conditions of Theorem 1.3.5.

We have seen that the Langevin dynamics satisfy conditions 1 and 2 of Theorem 1.3.5. Leimkuhler and al. proved in [START_REF] Leimkuhler | The Computation of Averages from Equilibrium and Non-Equilibrium Langevin Molecular Dynamics[END_REF] that any of the splitting schemes mentioned above satisfied the remaining conditions of Theorem 1.3.5. In particular, they proved that all schemes were ergodic with respect to a probability measure and that expansions such as (1.34) could be obtained. Therefore, by Theorem 1.3.5 and Lemma 1.3.8, we have that the above Trotter schemes are of order 1 in the invariant measure, while GLA and Strang splitting schemes are of order 2.

The Dissipative Particle Dynamics (DPD)

In this section we consider the Dissipative Particle Dynamics model and its energy-conserving variant. We first present some modeling principles of DPD, then introduce the SDE corresponding to the DPD, along with some consideration about the invariance of the canonical measure and the ergodicity of the dynamics. We then introduce the DPDE model and present some discussion about its ergodic properties.

Original DPD model

In the DPD framework, each particle represents a (group of) molecule(s) that are coarsegrained into one mesoscale particle centered around the center of mass of the molecule(s), described by its position q i and momenta p i , i being the index of each (cluster of ) molecule(s). The internal degrees of freedom of the molecule(s) are represented by friction and random interactions acting on the center of masses of the DPD particles, thus turning the deterministic differential equation satisfied by the atoms into a Stochastic Differential Equation satisfied by the center of masses of the DPD particles. Figure 2.1 illustrates the coarse-graining process of DPD by one example, where molecules composed of 4 atoms are coarse-grained into one single particle. Of course, the coarse-Figure 2.1 | Illustration of DPD coarse-graining process: each of the three molecules (purple atoms) interacting through conservative forces (straight lines: full for short-range interactions and dashed for long-range ones) are represented by a single meso-scale DPD particle (dashed circles) interacting through conservative, friction and random interactions (dashed oscillating lines). Note that DPD interactions are therefore stochastic and no longer deterministic. graining process of Figure 2.1 is only one example, and many other can be tried (for instance grouping several molecules into one DPD particle). DPD was put on a firm theoretical ground in [START_REF] Español | Statistical Mechanics of Dissipative Particle Dynamics[END_REF], and was applied to study the properties of various systems [MNZ09, DOS + 07, GMT08]. In the sequel, we only consider the center of masses of the mesoscale DPD particles that we denote by (q i , p i ) for each particle indexed with i ∈ N, and no longer consider the coordinates of the atoms.

Presentation of the dynamics

Consider a system of N DPD particles with periodic boundary conditions. Each pair of particles interacts through (i) conservative forces deriving from a potential energy function U(q), (ii) friction forces proportional to the relative velocity between the two particles and (iii) random fluctuation forces. The time evolution of the configuration (q i , p i ) ot the ith particle is given by the following set of equations:

                 dq i,t = p i,t m i dt, dp i,t = -∇ q i U(q t ) dt + N j=1,j i -γ ij χ(r ij,t ) v ij,t • e ij,t e ij,t dt + σ ij χ(r ij,t )e ij,t dW ij,t , (2.1)
where e ij,t = q i -q j q i,t -q j,t

is the directional vector going from particle i to particle j, and

v ij,t = v i,t -v j,t
is the relative velocity between particle i and particle j. The Wiener processes (W ij,t ) 1≤i<j≤N are independent and verify the antisymmetric condition W ij,t = -W ji,t . Note that in (2.1), the Wiener processes are unidimensional. The function χ(r) is a smoothing function limiting the range of the fluctuation/dissipation interactions. In this thesis, we choose

χ(r) =            1 - r r cut 2 , if r ≤ r cut , 0, otherwise. (2.2)
In the remainder of this thesis, and unless otherwise mentioned, we assume the potential U is such that ∇U satisfies Condition 1.2.2.

One can notice that the construction of the DPD dynamics is similar to the one of the Langevin dynamics (1.54) (see Section 1.4.1): friction and random terms are added to conservative forces deriving from the potential in the Hamiltonian dynamics. However, in the DPD framework, the friction and random terms are proportional to the relative velocities of the particles, thus verifying the Galilean invariance, in opposition to the Langevin dynamics. In addition, the total momentum P of the system is preserved by the dynamics, i.e

d        N i=1 p i        = 0. Remark 2.1.1

. [Control of the Coarse-Graining process] The DPD model regroups several particles into one mesoscale DPD particle. However, because DPD has no associated parameter to this regrouping process, the scale of the coarse-graining has to be controlled through the potential U(q). Indeed, the effective potential acting on the center of masses of the clusters changes according to the scale of the regrouping performed. Therefore, one must choose carefully its potential before performing any realistic DPD simulation.

In this thesis, we use a model introduced in [START_REF] Stoltz | A Reduced Model for Shock and Detonation Waves. Part I. The Inert Case[END_REF] where the friction and random terms are not projected along e ij,t and where the friction and random parameters are taken constant and equal for every couple, i.e

γ ij = γ, σ ij = σ, i, j = 1, .., N , i j.
The DPD equations therefore write

                 dq i,t = p i,t m i dt, dp i,t = -∇ q i U(q t ) dt + N j=1,j i -γχ(r ij,t )v ij,t dt + σ χ(r ij,t )dW ij,t .
(2.3)

Note that in this framework, the W ij are no longer scalars but d-dimensional Wiener processes.

Ergodic properties of DPD

The parameters γ and σ respectively control the friction and fluctuation strengths. They are chosen in order to ensure the correct statistical behavior of the system, i.e in order for the canonical measure µ β given by (

σ 2 = 2γ β . ( 2.4)
This result can be straightforwardly extended to the specific case of (2.3). Note that in practice we impose σ > 0, otherwise we would recover the Hamiltonian dynamics.

Proof. We prove the invariance of µ β for the general DPD equation given by (2.1). The adaptation of the proof to the reduced model given by (2.3) is trivial.

The generator of the DPD dynamics can be written as

L = L ham + 1≤i<j≤N L ij ,
where L ij is the generator of the elementary subdynamics

           dp i,t = -γ(v ij,t • e ij,t )χ(r ij,t )e ij,t dt + σ χ(r ij,t )e ij,t dW ij,t , dp j,t = γ(v ij,t • e ij,t )χ(r ij,t )e ij,t dt -σ χ(r ij,t )e ij,t dW ij,t .
We have shown in Section 1.4.1 that µ β is invariant for L ham . It thus remains to prove its invariance for the fluctuation/dissipation part, i.e prove that L * ij µ β = 0 for all i j, with (i, j) ∈ [1, .., N |] 2 . Under matrix form, the above elementary dynamics writes dx t = Γ(x t )dt + Σ(x t )dW ij,t , where x = (p i , p j ) and

Γ(x) = -γ(e ij • v ij )χ(r ij ) e ij -e ij , Σ(x) = σ χ(r ij ) e ij -e ij .
The operator L ij therefore writes

L ij = Γ(x) • ∇ + 1 2 Σ(x)Σ(x) T : ∇ 2 .
A straightforward computation shows that

L ij = χ(r ij ) -γ(e ij • v ij ) ∇ p i -∇ p j • e ij + σ 2 2 e ij • ∇ p i -∇ p j 2 .
We now consider the adjoint of L ij writing, for any ψ ∈ L 1 (X , R) regular enough,

L * ij ψ = χ(r ij ) γ e ij • ∇ p i -∇ p j (e ij • v ij )ψ + σ 2 2 e ij • ∇ p i -∇ p j 2 ψ . Applying L * ij to µ β gives L * ij µ β = χ(r ij ) γ e ij • ∇ p i -∇ p j (e ij • v ij )µ β + σ 2 2 e ij • ∇ p i -∇ p j 2 µ β , = χ(r ij ) e ij • ∇ p i -∇ p j γ(e ij • v ij )µ β - βσ 2 2 (e ij • v ij )µ β , = χ(r ij ) e ij • ∇ p i -∇ p j γ -β σ 2 2 (e ij • v ij )µ β .
We notice in the above expression that if (2.4) holds, then L * ij µ β = 0 for all (i, j) ∈ [|1, N |] 2 with i j, which concludes the proof.

Remark 2.1.2. We notice from the proof of Proposition 2.1.1 that, should the friction and fluctuation parameters depend on the particle couple (i, j), the fluctuation/dissipation relation would become

σ 2 ij = 2γ ij β , i, j = 1, .., N , i j.
So far, the ergodicity of the DPD dynamics for the canonical measure has been proved only in a one-dimensional case where the density is high enough [START_REF] Shardlow | Geometric Ergodicity for Dissipative Particle Dynamics[END_REF]. In the general case, the irreducibility of the dynamics cannot be proved. Indeed, assume that the position space is sufficiently large so there exists a configuration where all particles are isolated, i.e all the distances between particles are bounded from below by r min > r cut . Also assume that the particles have the same velocity, i.e T = 0. Because the random interactions are of the same limited range r cut than the conservative force, the time derivative of the momenta of the particles is null, and the system stays in the same configuration for ever. Such pathological cases prevent the DPD from being ergodic. Mathematically speaking, we do not have the hypoellipticity of the DPD generator and thus cannot apply Theorem 1.2.1.

In the sequel, we nevertheless assume that the ergodicity of the DPD dynamics holds for the canonical measure conditioned to the Dirac measure of the initial total momentum P 0 . This Dirac measure δ P 0 is constructed in a similar way than the microcanonical measure δ mc,E 0 of Section 1.2.1. The DPD invariant measure, denoted by µ β,DPD , there-fore writes µ DPD,β (dq, dp) = δ P 0 (p)e -βH(q,p) dqdp.

(2.5)

The Dissipative Particle Dynamics with Energy conservation (DPDE)

As we saw in the previous section, DPD intrinsically is an equilibrium model, with a prescribed temperature, fixed through a fluctuation/dissipation relation between the magnitudes of the friction and fluctuation parameters. DPD cannot be used as such to study nonequilibrium systems, and should be replaced by a dynamics where the fluctuationdissipation relation is not fixed a priori, but evolves depending on the physical events that have happened. DPD with conserved energy (DPDE) is such a model.

In the DPDE framework, mesoparticles have an additional degree of freedom, namely an internal energy denoted by ε i , which accounts for the energy of the missing degrees of freedom. This allows to define a total energy as

E(q, p, ε) = U(q) + N i=1 p 2 i m i + N i=1 ε i , (2.6)
where the last sum of the right-hand side of (2.6) is called the internal energy of the system. The dynamics on the internal energies is constructed in order for the total energy to remain constant. This energy conservation allows to use DPDE to simulate nonequilibrium phenomena like shock and detonation waves [Sto06, MSS07, MVDS11].

Presentation of the dynamics

Considering the DPD dynamics (2.3), an Itô calculus on the mechanical energy H(q, p) allows to write dH(q t , p t ) =

N i=1 ∇ q i U(q t ) • dq i,t + 1 m i         p i,t • dp i,t + N j=1,j i dσ 2 2 χ(r ij,t )dt         , = N i=1 N j=1,j i v i,t -γχ(r ij,t )v ij,t dt + σ χ(r ij,t )dW ij,t + dσ 2 2m i χ(r ij,t )dt, = 1≤i<j≤N v ij,t -γχ(r ij,t )v ij,t dt + σ χ(r ij,t )dW ij,t + dσ 2 2 χ(r ij,t ) 1 m i + 1 m j dt, = 1≤i<j≤N -γv 2 ij,t + dσ 2 2 1 m i + 1 m j χ(r ij,t )dt + σv ij,t χ(r ij,t )dW ij,t ,
where, to go from the second to the third line, we used a simple symmetry argument. We notice that the right-hand side of the above equality can be written as a sum of terms depending on the fluctuation/dissipation interaction of each particle couple ij, i.e dH(q t , p t ) =

N i=1 N j=1,j i dE ij,t , (2.7)
with

dE ij,t = 1 2 -γv 2 ij,t + dσ 2 2 1 m i + 1 m j χ(r ij,t )dt + σv ij,t χ(r ij,t )dW ij,t .
(2.8) Therefore, by defining the time evolution of the internal energy as dε i,t =j i dE ij,t , we obtain a dynamics preserving the energy E(q, p, ε) defined by (2.6). This dynamics therefore writes:

                               dq i,t = p i,t m i dt, dp i,t = -∇ q i U(q t ) dt + N j=1,j i -γ ij,t χ(r ij,t )v ij,t dt + σ ij,t χ(r ij,t )dW ij,t , dε i,t = 1 2         N j=1,j i        γ ij,t v 2 ij,t -d σ 2 ij,t 2 µ ij        χ(r ij,t ) dt -σ ij,t χ(r ij,t )v ij,t • dW ij,t         , ( 2.9)
where we have used the same notation than in Section 2.1.1 and defined

µ ij = 1 m i + 1 m j .
(2.10) Note that, because of the introduction of the internal energy variable, the statistical ensemble is changed and so must be the fluctuation/dissipation relation (2.4). Note also that the total momentum N i=1 p i is still preserved.

Remark 2.1.3. In the above model, the internal energy variations are all due to the friction terms of the momenta variations. Physically, these energy variations are interpreted as the dissipated energy by viscous heating of the system. The original DPDE model [START_REF] Avalos | Dissipative Particle Dynamics with Energy Conservation[END_REF][START_REF] Español | Dissipative Particle Dynamics with Energy Conservation[END_REF] considers the DPDE particle to also exchange energy through heat conduction, the viscous heating term being interpreted as the mechanical part of the energy variation. These exchanges are caused by the temperature difference between the particles, and write

dε cond i = N j=1,j i κ ij χ(r ij,t ) 1 T i - 1 T j dt + α ij χ(r ij,t )dW ′ ij,t ,
where W ′ ij,t are one-dimensional Wiener processes independent of the W ij,t also verifying dW ′ ij,t = -dW ′ ji,t , and (κ ij , α ij ) 1≤i j≤N the parameters controlling the amplitude of the conduction thermal exchanges. With the heat conduction taken into account, the resulting DPDE dynamics write

                                             dq i,t = p i,t m i dt, dp i,t = -∇ q i U(q t ) dt + N j=1,j i -γ ij,t χ(r ij,t )v ij,t dt + σ ij,t χ(r ij,t )dW ij,t , dε i,t = 1 2         N j=1,j i        γ ij,t v 2 ij,t -d σ 2 ij,t 2 µ ij        χ(r ij,t ) dt -σ ij,t χ(r ij,t )v ij,t • dW ij,t         N j=1,j i κ ij χ(r ij,t ) 1 T i - 1 T j dt + α ij χ(r ij,t )dW ′ ij,t .
(2.11)

Similarly to the DPDE dynamics, the above dynamics preserves the total energy.

Thermodynamics of the DPDE model

With the introduction of the internal energy ε, we have to also define the internal thermodynamics of the particles. In the DPDE model, internal degrees of freedom are assumed to relax towards equilibrium at much faster rates than the external ones. They thus are considered to always be at equilibrium, and we can define the internal equation of state of each particle as

ε = T (ε) 0 C v (θ) dθ, (2.12)
where, C v (θ) is the internal heat capacity at constant volume of the particle and T (ε) is its internal temperature. In order for T (ε) to be well-defined by (2.12), we make the following assumptions:

1. The heat capacity C v is continuous on R + .

2. The heat capacity C v is positive on R * + and satisfies lim

θ→0 + C v (θ) = C - v ∈ [0, +∞[.
The two above condition imply that T (ε) exists and is unique for every ε ≥ 0, and positive when ε > 0. In addition, we see that T (ε) is positive and continuous on R + (because C v also is) and satisfies lim ε→0 + T (ε) = 0.

Remark 2.1.4. In the classical description of matter, the heat capacity is increasing on R + , has a null limit when θ goes to 0 (i.e C - v = 0), and tends towards C lim v < +∞ when θ goes to infinity. In many physical systems, the limit C lim v is reached whenever θ is greater than several hundreds of Kelvin. This is the heat capacity at ambient temperature, for systems where quantum effects do not apply.

After having defined the internal energy and temperature of a DPDE particle, we now define its internal entropy, denoted by s(ε). This quantity depends solely on the internal energy and is defined from the relation

s ′ (ε) = 1 T (ε) . (2.13)
Again, we see that s ′ (ε) is well defined (up to a constant) for every ε > 0 and is strictly increasing and continuous on R * + (because T (ε) > 0 for ε > 0 and is continuous). In addition, the vanishing limit of T (ε) when ε goes to zero, the continuity and the strict growth of s(ε) imply lim ε→0 + s(ε) = -∞. The total internal entropy, denoted S , is defined as the sum of the internal entropies of the particles, i.e

S (ε) = N i=1 s i (ε i ).
Note that in the above formula, each particle has a different internal entropy function. In many applications, the internal heat capacity is assumed to be constant and independent of the internal temperature, i.e C v,i (θ) = C v,i = cst. In such cases, the internal temperature and energy are proportional, and the internal entropy can be easily computed from (2.13), i.e i = 1, .., N ,

ε i = C v,i T i , s i (ε i ) = C v,i log(ε i ).
We see from (2.12) and (2.13) that the knowledge of the internal heat capacity and internal energy fully determines the internal temperature and entropy. The original DPD model was designed to sample the canonical measure µ β,DPD , given by (2.5) and defined as the canonical measure µ β conditioned with the Dirac measure δ P 0 on the total momentum. For the DPDE framework, the thermodynamical ensemble should take into account the additional degrees of freedom represented by ε.

The total energy E and the total momentum P are both exactly preserved by the dynamics. Therefore, one can define the associated Dirac measures δ E 0 and δ P 0 respectively to a fixed energy E 0 and total momentum P 0 with the procedure outlined in Section 1.2.1. The DPDE "microcanonical" macrostate µ E 0 ,P 0 is then defined as the product of these two measures, i.e µ E 0 ,P 0 (dq, dp, dε) = Z -1 E 0 ,P 0 δ E 0 (q, p, ε)δ P 0 (p)e

S (ε) k B dqdpdε, (2.14)
where Z E 0 ,P 0 is a normalization constant. The entropy term e

S (ε)

k B comes from the DPDE coarse-graining process of integrating the missing degrees of freedom into one internal energy variable.

In the remainder of this thesis, we assume that the DPDE dynamics is ergodic with respect to µ E 0 ,P 0 . We give later on in Proposition 2.1.2 a necessary condition on γ ij and σ ij in order to obtain the invariance of µ E 0 ,P 0 for the dynamics.

Similarly to the microcanonical ensemble, the DPDE "microcanonical" ensemble is not easy to manipulate for the computations of thermodynamical averages. Therefore, we can construct a DPDE "canonical" ensemble, which should be equivalent to (2.14) in the thermodynamical limit, as

µ β,C v (dq, dp, dε) = δ P 0 (p)Z -1 β,C v e -βE(q,p,ε)+ S (ε) k B dqdpdε, (2.15)
where Z β,C v is a normalization constant. In this thesis, we assume that the equivalence of ensemble between µ E 0 ,P 0 and µ β,C v holds in the thermodynamical limit. According to Section 1.2.1, this means that, provided β is chosen such that

E µ β,C v [E] = E 0 , we have lim N →∞ E µ (N ) E 0 ,P 0 [ϕ(x 1 , .., x k )] -E µ (N ) β,C v [ϕ(x 1 , .., x k )] = 0,
for any observable ϕ ∈ D depending on a fixed number k ∈ N of variables (x 1 , .., x k ), and where µ (N ) E 0 ,P 0 and µ

(N )
β,C v are respectively the "microcanonical" and "canonical" DPDE statistical ensembles corresponding to a system of N particles.

Ergodic properties of DPDE

Similarly to the DPD model, there is no result for ergodicity of the DPDE model, and the same problem arise when one wants to prove the irreducibility of the dynamics. However, the parameters γ ij and σ ij can be chosen so that µ E 0 ,P 0 is invariant for the dynamics. Proposition 2.1.2. Assume that for all (i, j) ∈ [|1, N |] 2 with i j, the parameters γ ij and σ ij satisfy

σ ij = σ, γ ij = σ 2 4k B 1 T i (ε i ) + 1 T j (ε j ) , (2.16)
where σ is a constant. Then the macrostate µ E 0 ,P 0 defined by (2.14) is invariant for the DPDE dynamics described by (2.9).

Similarly to DPD, we always impose σ > 0.

Remark 2.1.5. When the heat conduction is taken into account like in (2.11), the parameters κ ij and α ij can also satisfy a fluctuation/dissipation relation in order to obtain the invariance of µ E 0 ,P 0 . This relation reads, for all (i, j)

∈ [|1, N |] 2 with i j, α ij = α, κ ij = κ = 2α 2 .
(2.17)

In practice, we also always impose α > 0.

Proof. The proof is organized as follows:

1. We first prove that µ β,C v is invariant provided that γ ij and σ ij satisfy (2.16).

2. We extend the proof to any probability measure of the form µ f ,g (dq, dp, dε) = Z -1 f ,g f E(q, p, ε) g

       N i=1 p i        exp S (ε) dq dp dε, (2.18)
for arbitrary functions f and g sufficiently smooth such that Z f ,g = X dµ f ,g < +∞.

3. We then extend the proof to µ E 0 ,P 0 by noticing that, for any test function ϕ ∈ D ∞ , we have

X Lϕdµ f ,g = +∞ -∞ +∞ -∞ X Lϕdµ E 0 ,P 0 f (E 0 )g(P 0 )dE 0 dP 0 . Because µ f ,g is invariant, we have X Lϕdµ f ,g = 0 for all test function ϕ ∈ D ∞ , which implies ∀ϕ ∈ D ∞ , X Lϕdµ E 0 ,P 0 = 0.
This proves the invariance of µ E 0 ,P 0 .

Let us prove that µ β,C v is invariant under the assumption that γ ij and σ ij satisfy (2.16). The generator of the DPDE dynamics can be written as

L = L ham + 1≤i<j≤N L ij ,
with L ij the generator of the following elementary subdynamics:

                           dp i,t = -γ ij,t v ij,t χ(r ij,t ) dt + σ ij,t χ(r ij,t ) dW ij,t , dp j,t = -dp i,t , dε i,t = 1 2               γ ij,t v 2 ij,t -d σ 2 ij,t µ ij 2        χ(r ij,t )dt -σ ij,t v ij,t • dW ij,t        , dε j,t = dε i,t .
From the computations of Section 1.4.1, we obtain that µ β,C v is invariant for L ham . It thus remains to prove its invariance for L ij . The operators L ij do not involve the positions, therefore χ(r ij ) can be considered as constant and we do not write it in the remainder of the proof for clarity purposes.

Under matrix form, the above elementary dynamics writes

dx t = Γ(x t )dt + d k=1 Σ k (x t )dW ij,k,t , where x = (p i , p j , ε i , ε j ) ∈ R 2d × R 2 + , W ij,k,t is the kth component of the d-dimensional Wiener process W ij,t and Γ(x) =                     -γ ij v ij γ ij v ij 1 2 γ ij v 2 ij -d σ 2 ij µ ij 2 1 2 γ ij v 2 ij -d σ 2 ij µ ij 2                     , Σ k (x) = σ ij                1 d,k -1 d,k - v ij,k 2 - v ij,k 2               
, where 1 d,k is the d-dimensional vector whose components are equal to 0 except the kth component equal to 1. The operator L ij therefore writes

L ij = Γ(x) • ∇ x + 1 2 d k=1 2d+2 a,b=1 (Σ k (x)) a (Σ k (x)) b ∂ a ∂ b .
A straightforward computation shows that

L ij = -γ ij v ij • ∇ p i -∇ p j - 1 2        γ ij v 2 ij -d σ 2 ij µ ij 2        ∂ ε i + ∂ ε j + σ 2 ij 2        ∇ p i -∇ p j 2 -v ij • ∇ p i -∇ p j ∂ ε i + ∂ ε j + v 2 ij 4 ∂ ε i + ∂ ε j 2        . ( 2

.19)

Let us define

A p,ij = ∇ p i -∇ p j , A ε,ij = ∂ ε i + ∂ ε j , A ij = A p,ij - v ij 2 A ε,ij .
We notice that

                 A ij (E(q, p, ε)) = 0, A ij p i = 0, A ij (S (ε)) = - v ij 2 s ′ (ε i ) + s ′ (ε j ) = - v ij 2 1 T i + 1 T j .
(2.20)

We additionally have that

A 2 ij = A 2 p,ij -A p,ij v ij 2 A ε,ij - v ij 2 A ε,ij A p,ij + v 2 ij 4 A 2 ε,ij , =        A 2 p,ij -v ij • A p,ij A ε,ij + v 2 ij 4 A 2 ε,ij        -d µ ij 2 A ε,ij .
We notice that both terms between hooks of (2.19) and the above expression correspond. This allows us to rewrite L ij as

L ij = -γ ij v ij • A ij + σ 2 ij 2 A 2 ij .
We now consider the adjoint of L ij writing, for any ψ ∈ L 1 (X , R) regular enough,

L * ij ψ = A ij • (γ ij v ij ψ) + A ij • A ij        σ 2 ij 2 ψ        , = A ij •        γ ij v ij ψ + A ij        σ 2 ij 2 ψ               , = A ij •        γ ij v ij ψ + A ij        σ 2 ij 2        ψ + σ 2 ij 2 A ij (ψ)        .

Using (2.20) yields

A ij (µ β,C v ) = - v ij 2k B 1 T i + 1 T j µ β,C v .
Therefore, we have

L * ij µ β,C v = A ij •        γ ij v ij µ β,C v + 1 2 A ij (σ 2 ij )µ β,C v - σ 2 ij v ij 4k B 1 T i + 1 T j µ β,C v        , = A ij •               γ ij - σ 2 ij 4k B 1 T i (ε i ) + 1 T j (ε j )        v ij + 1 2 A ij (σ 2 ij )        µ β,C v .
We notice in the above expression that if we assume (2.16) to hold, we have

L * ij µ β,C v = 0.
We notice that measures µ f ,g of the form (2.18) also satisfy

A ij ρ f ,g = - v ij 2k B 1 T i + 1 T j ρ f ,g
, where ρ f ,g denotes the density of µ f ,g . This gives L * ij µ f ,g = 0 provided that (2.16) holds and proves the invariance of the probability measures µ f ,g of the form (2.18). According to point 3 of the plan at the beginning of the proof, this gives the result. Remark 2.1.6. Equation (2.16) is not the only possibility for measures of the form (2.18) to be invariant for the DPDE dynamics. However, it is the simplest and most widely used relation.

In order to highlight the similarity with the standard fluctuation/dissipation relation for Langevin dynamics (see for instance Section 2.2.3 in Ref. [START_REF] Lelièvre | Free Energy Computation : A Mathematical Perspective[END_REF]), it might be convenient to rewrite the definition of γ ij and σ as

γ ij = γ β ij β , σ = 2γ β , (2.21) with β ij = 1 2k B 1 T i + 1
T j and γ a reference friction parameter.

Sampling the DPDE "canonical" measure

We have seen that the ergodicity of the DPDE dynamics cannot be proved, mainly because the irreducibility of the dynamics seems out of range. However, under some assumptions that we assume to hold without proof, we can devise a procedure to sample the DPDE "canonical" measure µ β,C v given by (2.15) by devising ergodic dynamics for its marginals.

The macrostate (2.15) can be decomposed into a marginal on the positions and momenta (q, p) and a marginal in the internal energies, i.e µ β,C v (dq, dp) = µ DPD,β (q, p)ν β,C v (ε)dqdpdε, where µ DPD,β abusively denotes the density of the DPD invariant measure given by (2.5), and ν β,C v abusively denotes the density of the probability measure defined as

ν β,C v (dε) = Z -1 ε e -β N i=1 ε i + S (ε) k B dε = N i=1 Z -1 ε,i e -βε i + s i (ε i ) k B dε i , (2.22)
where Z ε is a normalization constant defined as

Z ε = +∞ 0 e -β N i=1 ε i + S (ε) k B dε = N i=1 +∞ 0 e -βε i + s i (ε i ) k B dε i Z ε,i .
Sampling µ DPD,β can be done by the Langevin dynamics presented in Section 1.4. However, the Langevin dynamics does not preserve the total momentum. The DPD dynamics does preserve the total momentum, but we do not have any ergodicity result of DPD with respect to µ DPD,β . Therefore, we must assume that we can sample µ DPD,β with realizations of the trajectories of solutions of the Langevin dynamics, where the configurations are constantly projected on N (P 0 ), where

N (P 0 ) =        (q, p) ∈ X q,p N i=1 p i = P 0        ,
with X q,p being the configuration space of the momenta and positions, i.e X = X q,p × R * + . The projection can be done by subtracting α n = 1

N N i=1

pn i -P 0 to each pn i obtained by the application of one iteration of a given discretization of the Langevin dynamics, i.e i = 1, .., N , Φ

proj ∆t,P 0 (q, p, G) =        Φ (q) ∆t (q, p, G), Φ (p) ∆t (q, p, G) - 1 N        N i=1 Φ (p) ∆t (q, p, G) i -P 0        1 dN        ,
where

Φ ∆t = Φ (q) ∆t , Φ (p)
∆t is any given discretization of the Langevin dynamics, Φ proj ∆t its associated projected discretization and 1 dN the dN -dimensional vector whose components are all equal to 1.

Sampling ν β,C v , is done by considering the Overdamped Langevin dynamics which reads dε

i,t = -U ′ i,β (ε i,t )dt + 2 β dB i,t ,
where B i,t is a one-dimensional Brownian motion independent of the Brownian motion of the Langevin dynamics and the potential U i,β is defined as

U i,β (ε i ) = ε i - s i (ε i ) k B β ,
with s i being the internal entropy functions defined by (2.13). We have seen in Section 1.4.1 that the dynamics is ergodic with respect to the probability measure ρ of density Z -1 i e -βU i,β (Z i being a normalization constant), provided that U ′ i,β satisfies Condition 1.2.2. In our case however, the derivative of U i,β writes

U ′ i,β (ε i ) = 1 - 1 βk B T i (ε i ) ,
and we have seen that the DPDE internal temperature goes to zero when ε goes to zero, which forbids ∇U i,β to satisfy Condition 1.2.2 and might cause instabilities when discretizing the above dynamics. Nevertheless, we never encountered such problems in our simulations. By injecting the expression of U i,β into the dynamics, we obtain

dε i,t = 1 - 1 βk B T i (ε i,t ) dt + 2 β dB i,t , (2.23)
and we assume that such a dynamics admits solutions on R * + that sample ν β,C v . Therefore, under the two previous assumptions, we can sample µ β,C v by sampling independently the positions and momenta by a Langevin dynamics projected on N (P 0 ) and the internal energies by the dynamics described by (2.23).

Temperature estimators in the DPDE context

Assuming the equivalence of the canonical and microcanonical ensembles in the DPDE context, we can define, for instance, various estimators of the temperature. Each of the three estimators of the temperature presented below involves only one of the three cat-egories of degrees of freedom of the system: the positions q, the momenta p and the internal energies ε.

• Kinetic temperature: A first estimator is the famous kinetic temperature defined as

T kin = 1 k B N eff E µ β,C v        N i=1 p 2 i m i        , (2.24)
where N eff represents the effective number of external degrees of freedom of the system. It is a priori equal to dN but since we fix the total momentum to P 0 = 0, we reduce it to N eff = d(N -1). This correction is anyway unimportant for sufficiently large systems. An easy computation shows that for every particle indexed by i = 1, .., N , we have

E µ β       p 2 i m i       = z -1 i,β R p 2 m i e -β p 2 2m i dp = β,
where z i,β = R e -βp 2 /(2m i ) dpi = 2πm i /β. We therefore have T kin = T , where T is the equilibrium temperature corresponding to β. We assume that we still have T kin = T when the expectation of (2.24) is on µ β,C v instead of µ β .

• Potential temperature: Butler and al. proved in [START_REF] Butler | Configurational Temperature: Verification of Monte Carlo Simulations[END_REF] that a second estimator of the temperature could be defined in the canonical ensemble involving only the positions q. This estimator, called the configurational temperature, or potential temperature is defined as

T pot = 1 k B E µ β,C v ∇U 2 E µ β,C v [∆U] . ( 2

.25)

The proof of T pot = T is done thanks to a simple computation. Assume that U is such that ∇U(q)e -βU(q) ---------→ U(q)→+∞ 0, i.e that ∇U(q) does not increase too fast when U(q) increases. We therefore have

Q ∇ • ∇Ue -βU dq = 0, Q ∆U -β ∇U 2 e -βU dq = 0, Q ∆Ue -βU dq = β Q ∇U 2 e -βU dq, which means that E µ β,C v [∆U] = βE µ β,C v [ ∇U 2 ]
, and therefore that T pot = T .

• Internal temperature: Mackie and al. noticed in [MAN99] that for every particle indexed by i = 1, .., N , an integration by part gives us

E 1 T i µ β,C v = z -1 β,C v +∞ 0 s ′ (ε)e -βε+ s(ε) k B dε = z -1 β,C v +∞ 0 k B βe -βε+ s(ε) k B dε = 1 T ,
where z β,C v = +∞ 0 e -βε+s(ε)/k B dε. We therefore can define a third estimator of the temperature, involving only the internal degrees of freedom, as

T int =        E µ β,C v        1 N N i=1 1 T i               -1 = T .
(2.26)

The above three estimators approximate the same thermodynamical average, and therefore should converge to the same value for systems at equilibrium. By the equivalence of µ E 0 ,P 0 and µ β,C v , we therefore have three estimators of the temperature in the DPDE "microcanonical" ensemble in the thermodynamic limit.

Remark 2.1.7. The kinetic temperature T kin is also sometimes called the external temperature in the DPDE context by opposition to the internal temperature T int .

Integrating the Dissipative Particle Dynamics with Energy-conservation: state of the Art

We have presented both the Dissipative Particle Dynamics framework and its energyconserving variant in Section 2.1. We now focus on the existing numerical algorithms used for the integration of DPD.

Among the currently known schemes for DPDE, one of them turns out to enjoy very nice properties in terms of energy conservation. It is a splitting scheme inspired from Shardlow's splitting scheme for DPD [START_REF] Shardlow | Splitting for Dissipative Particle Dynamics[END_REF], and adapted to the DPDE framework [START_REF] Stoltz | A Reduced Model for Shock and Detonation Waves. Part I. The Inert Case[END_REF]. The resulting scheme, called SSA in this thesis (Shardlow Splitting Algorithm), is so far considered as the reference scheme for the numerical integration of DPDE [START_REF] Lisal | Dissipative Particle Dynamics at Isothermal, Isobaric, Isoenergetic, and Isoenthalpic Conditions using Shardlow-like Splitting Algorithms[END_REF]. The main issue with SSA is that particle pairs have to be updated sequentially, thus preventing a simple parallelization of the scheme. The parallelization of SSA is therefore not an easy task [LBM + 14].

The remaining options are whether implicit schemes, or explicit schemes with very poor stability and energy preservation properties. Implicit schemes are very expensive computationally because they require the resolution of a fixed point algorithm at each iteration, which sometimes takes up most of the simulation time. Therefore, we have discarded implicit schemes in this thesis.

Most of the schemes developed so far are splitting schemes, introduced in Section 1.3 and developed in Section 1.4. These schemes split the DPDE dynamics into a Hamiltonian evolution given by (1.4), and a fluctuation/dissipation part given by

                           dq i,t = 0, dp i,t = - N j=1,j i γ ij,t χ(r ij,t )v ij,t dt + N j=1,j i σ χ(r ij,t ) dW ij,t , dε i,t = 1 2         N j=1,j i       γ ij,t v 2 ij,t -d µ ij σ 2 2       χ(r ij,t ) dt -σ χ(r ij,t )v ij,t • dW ij,t         .
(2.27)

The discretization of the Hamiltonian evolution is in all cases performed with a standard Velocity-Verlet integration [Ver67], detailed in Section 1.1.2. Thus, all the splitting schemes presented in this section differ only in the discretization of the fluctuation/dissipation part. We present in this section three representative numerical schemes for the DPDE integration. As said in the introduction of this chapter, we limit ourselves to explicit algorithms:

• The first scheme is an adaptation to the DPD framework of the well-known Velocity-Verlet algorithm called Stochastic Velocity-Verlet (SVV). The SVV scheme is not a splitting scheme, and has poor invariant conservation properties.

• The second scheme is a basic splitting scheme called Splitting Euler-Maruyama (SEM), where no particular care is taken to make sure that the energy is preserved by the discretization of the fluctuation/dissipation part. It is presented for pedagogical purposes, since the later schemes of Section 2.3 are based upon it, and also for comparison purposes.

• The third scheme is the Shardlow's Splitting Algorithm (SSA). As previously said, SSA has good energy conservation properties, and is considered to be the best performing scheme for DPD and DPDE [START_REF] Lisal | Dissipative Particle Dynamics at Isothermal, Isobaric, Isoenergetic, and Isoenthalpic Conditions using Shardlow-like Splitting Algorithms[END_REF]. However, SSA suffers from parallelization issues, thus making it more difficult to implement for massively-parallel simulations [LBM + 14].

We then give some explanation about the parallelization issue of SSA, and give some details about its parallelization procedure presented in [LBM + 14].

In the following sections, we will use the notion of "uniformly bounded for small ∆t" for remainder terms r ∆t,ϕ depending on the time step and on some function ϕ. This means that there exists ∆t * > 0, C ϕ > 0 and m ∈ N such that r ∆t,ϕ

L ∞ W m < C ϕ for any ∆t ∈ [0, ∆t * ].
Definition 2.2.1. Consider a function ϕ : X → R, ∆t > 0 and a function r ∆t,ϕ : X → R depending on ϕ and ∆t. We say that r ∆t,ϕ is uniformly bounded for small ∆t if

∃(∆t * , C ϕ , m) ∈ (R * + ) 2 × N, s.t (0 < ∆t < ∆t * ) ⇒ r ∆t,ϕ L ∞ W m < C ϕ .

Stochastic Velocity-Verlet Algorithm (SVV)

The Stochastic Velocity-Verlet algorithm, denoted SVV, is an adaptation to the DPD setting of the well-known Velocity-Verlet algorithm integrating the Hamiltonian dynamics [START_REF] Verlet | Computer "experiments" on classical fluids. i. thermodynamical properties of lennard-jones molecules[END_REF]. It consists in using a Strang splitting of the DPD, where the positions q i are updated once with a time step ∆t, and the momenta p i are updated twice, both before and after the position update, with a time step ∆t 2 . Each part (the position and momenta evolution) is integrated with an Euler-Maruyama discretization.

However, this technique requires the use of a set of random numbers of double size, which is expensive in terms of computational costs: random numbers are generally much more expensive to generate than performing a single standard float operation. Fortunately, it can be shown that performing half-updates with a time step ∆t and using the same set (G n ij ) in both is also weakly consistent, and of the same weak order. This boils down to taking out the 1/2 term out of the square root in the random terms and using the same random numbers in both momenta updates. The SVV for DPD thus writes

                                       p n+1/2 i = p n i -∇ q i U(q n ) ∆t 2 + N j=1 j i -γχ(r n ij )v n ij ∆t 2 + σ χ(r n ij )G n ij √ ∆t 2 , q n+1 i = q n i + ∆t p n+1/2 i m i , p n+1 i = p n+1/2 i -∇ q i U(q n+1 ) ∆t 2 - N j=1 j i γχ(r n+1 ij )v n+1/2 ij ∆t 2 + σ χ(r n+1 ij )G n ij √ ∆t 2 ,
where (here and in the sequel) (G n ij ) n∈N,1≤i<j≤N is a set of identically and independently distributed standard d-dimensional Gaussian random variables, verifying G n ij = -G n ji . Proposition 2.2.1. The SVV discretization of the DPD dynamics is weakly consistent, of weak order one. In addition, its evolution operator satisfies ∀ϕ ∈ D, P SVV ∆t ϕ = e ∆tL ϕ + ∆t 2 r ∆t,ϕ , where r ∆t,ϕ is uniformly bounded for small ∆t in the sense of Definition 2.2.1 and L is the generator of the DPD dynamics.

Proof. Because the proof is performed for the SVV DPD discretization, X refers for the duration of the proof to the configuration space of the position q and the momenta p only, and no longer includes the internal energies ε. This proof consists in showing that the SVV discretization can be rewritten as a perturbation of an Euler-Maruyama discretization, where the expectation of the difference between the two discretizations is of order ∆t 2 . Denoting respectively by Φ SVV ∆t and Φ EM ∆t the results of an SVV and an Euler-Maruyama discretization of the DPD dynamics, this means writing

Φ SVV ∆t (q n , p n ) = Φ EM ∆t (q n , p n ) + η n ∆t 3/2 + r n ∆t ∆t 2 , (2.28)
where η n is a term of null average, i.e E G n [η n ] = 0, and r n ∆t is a remainder term satisfying

r n ∆t = O(1 + |G n | a ) for some a (i.e E G n [r n
∆t ] stays bounded when ∆t goes to zero). Assuming (2.28), fix ϕ ∈ D and (q, p) ∈ X . We have

E G n ϕ Φ SVV ∆t (q, p) = E G n ϕ Φ EM ∆t (q, p) + η n ∆t 3/2 + r n ∆t ∆t 2 , = E G n ϕ Φ EM ∆t (q, p) + η n • ∇ϕ(Φ EM ∆t (q, p))∆t 3 2 + r∆t,ϕ ∆t 2 , = E G n ϕ Φ EM ∆t (q, p) + r∆t,ϕ ∆t 2 ,
where r∆t,ϕ is obtained from the remainder terms of order ∆t 2 or greater in the Taylor expansion of ϕ. Both r n ∆t = O(1 + |G n | a ) and ϕ ∈ D imply that r∆t,ϕ is uniformly bounded in some weighted space for small ∆t. We thus have ∆t,ϕ is uniformly bounded in some weighted space for small ∆t. Therefore, we can expand P SVV ∆t ϕ as P SVV ∆t ϕ = e ∆tL ϕ + r ∆t,ϕ ∆t 2 , where r ∆t,ϕ comes from r∆t,ϕ and r EM ∆t,ϕ and is uniformly bounded in some weighted space for small ∆t. Theorem 1.3.3 then gives us the weak order one of SVV.

P SVV ∆t ϕ = P EM ∆t ϕ + r∆t,ϕ ∆t 2 . ( 2
It remains to prove (2.28). From now on, we do not make precise the remainder terms, and use notations like O(∆t ω ) to denote them. For details on how to handle such terms more precisely, we refer the reader to the proof of Proposition 1.3.4.

Using the expression of p n+1/2 i of the SVV discretization, we can expand q n+1 i as

q n+1 i = q n i + p n i m i ∆t +           N j=1,j i σ χ(r n ij ) 2 G n ij           η n q,i ∆t 3/2 + O(∆t 2 ). ( 2

.30)

We notice that the order 3/2 term of the above expansion has a vanishing expectation, i.e

E G n [η n q,i ] = 0
, where E G is the expectation on all the Gaussian increments (G n ij ) 1≤i j≤N .

Let us perform a similar expansion for p n+1 i . A Taylor expansion on ∇ q i U(q n+1 ) allows us to write ∇ q i U(q n+1 ) = ∇ q i U(q n ) + O(∆t).

(2.31)

In addition, (2.30) holding for any i = 1, .., N , we can expand r n+1 ij as

r n+1 ij = q n+1 i -q n+1 j = r n ij + v n ij 2 ∆t + O(∆t 3/2 ), thus yielding
χ(r n+1 ij ) = χ(r n ij ) + χ ′ (r n ij ) v n ij 2 k n χ,ij ∆t + O(∆t 3/2 ). (2.32)
The expression of p n+1/2 i given by the SVV discretization gives us

v n+1/2 ij = v n ij +           N k=1,k i σ χ(r n ik ) m i G n ik - N k=1,k j σ χ(r n jk ) m j G n jk           η n v,ij √ ∆t + O(∆t), (2.33)
where

E G [η n v,ij ] = 0 for all (i, j) ∈ [|1, N |] 2 .
Therefore, injecting (2.31), (2.32) and (2.33) into the expression of p n+1 i given by the SVV discretization, we obtain

p n+1 i = p n+1/2 i -∇ q i U(q n+1 ) ∆t 2 + N j=1,j i           -γ n ij χ(r n+1 ij )v n+1/2 ij ∆t 2 + σ χ(r n+1 ij )∆t 2 G n ij           , = p n+1/2 i -∇ q i U(q n ) + O(∆t) ∆t 2 + N j=1,j i -γ n ij χ(r n ij ) + O(∆t) v n ij + η n v,ij √ ∆t + O(∆t) ∆t 2 + σ χ(r n ij )∆t + k n χ,ij ∆t 2 + O(∆t 3 ) 2 G n ij , = p n i -∇ q i U(q n )∆t + N j=1,j i -γ n ij χ(r n ij )v n ij ∆t + σ χ(r n ij )∆tG n ij +           N j=1,j i - γ n ij χ(r n ij )η n v,ij 2 + σk n χ,ij 4 χ(r n ij ) G n ij           η n p,i ∆t 3/2 + O(∆t 2 ).
We notice that the order 3/2 term of the above expansion has a vanishing expectation, i.e

E G n [η n p,i ] = 0.
Combining the above expansions of q n+1 i and p n+1 i and using that

E G n [η n q,i ] = E G n [η n p,i ] = 0, we obtain (2.

28).

A "natural" extension of the SVV algorithm to the DPDE framework was given in 2011 by Lisal and al. [START_REF] Lisal | Dissipative Particle Dynamics at Isothermal, Isobaric, Isoenergetic, and Isoenthalpic Conditions using Shardlow-like Splitting Algorithms[END_REF]. It consists in updating the internal energies along the momenta with an Euler-Maruyama discretization. By denoting γ n ij = γ ij (ε n i , ε n j ), this extension of SVV for the DPDE, still denoted SVV, writes

                                                                       p n+1/2 i = p n i -∇ q i U(q n ) ∆t 2 + N j=1 j i -γ n ij χ(r n ij )v n ij ∆t 2 + σ χ(r n ij )G n ij √ ∆t 2 , ε n+1/2 i = ε n i + 1 2 N j=1 j i       γ n ij (v n ij ) 2 -d µ ij σ 2 2       χ(r n ij ) ∆t 2 -σv n ij χ(r n ij )G n ij √ ∆t 2 , q n+1 i = q n i + ∆t p n+1/2 i m , p n+1 i = p n+1/2 i -∇ q i U(q n+1 ) ∆t 2 - N j=1 j i γ n+1/2 ij χ(r n+1 ij )v n+1/2 ij ∆t 2 + σ χ(r n+1 ij )G n ij √ ∆t 2 , ε n+1 i = ε n+1/2 i + 1 2 N j=1 j i       γ n+1/2 ij v n+1/2 ij 2 -d µ ij σ 2 2       χ(r n+1 ij ) ∆t 2 -σv n+1/2 ij χ(r n+1 ij )G n ij √ ∆t 2 ,
(2.34) with µ ij is defined in (2.10). Note that, again, the same random number is used in the first and second updates of the momenta. The weak order one of the SVV discretization of DPDE can be obtained by adapting the proof concerning the SVV discretization of the DPD dynamics. We however omit the proof.

Proposition 2.2.2. The SVV discretization of the DPDE dynamics is weakly consistent of weak order one. In addition, its evolution operator satisfies ∀ϕ ∈ D, P SVV ∆t ϕ = e ∆tL ϕ + ∆t 2 r ∆t,ϕ , where r ∆t,ϕ is uniformly bounded for small ∆t in the sense of Definition 2.2.1 and L is the generator of the DPDE dynamics.

Splitting Euler-Maruyama (SEM)

From now on, and unless otherwise stated, L always denote the generator of the DPDE dynamics given by (2.9).

The SEM scheme integrates the fluctuation/dissipation with a simple Euler-Maruyama discretization:

P SEM,fd ∆t :                      p n+1 i = p n i - N j=1,j i γ n ij χ(r n ij )v n ij ∆t + σ χ(r n ij )G n ij √ ∆t, ε n+1 i = ε n i + 1 2 N j=1,j i       γ n ij (v n ij ) 2 -d µ ij σ 2 2       χ(r n ij )∆t -σ χ(r n ij )v n ij • G n ij √ ∆t, (2.35)
Note that no particular care is taken to make sure that the energy is preserved by the discretization of the fluctuation/dissipation part.

The global SEM discretization is thus characterized by the transition operator

P SEM ∆t = P VV ∆t P SEM,fd ∆t .
The corresponding numerical scheme reads

                                             p n+1/2 = p n - ∆t 2 ∇U(q n ), q n+1 = q n + ∆tp n+1/2 , pn+1 = p n+ 1 2 - ∆t 2 ∇U(q n+1 ), p n+1 i = pn+1 i - N j=1,j i γ n ij χ(r n+1 ij ) ṽn+1 ij ∆t + σ χ(r n+1 ij )G n ij √ ∆t, ε n+1 i = ε n i + 1 2 N j=1,j i       γ n ij ( ṽn+1 ij ) 2 -d µ ij σ 2 2       χ(r n+1 ij )∆t -σ χ(r n+1 ij ) ṽn+1 ij • G n ij √ ∆t.
Proposition 2.2.3. The SEM discretization of the DPDE dynamics (2.9) is of weak order one. In addition, its evolution operator satisfies ∀ϕ ∈ D, P SEM ∆t ϕ = e ∆tL ϕ + ∆t 2 r ∆t,ϕ , where r ∆t,ϕ stays uniformly bounded for small ∆t in the sense of Definition 2.2.1.

Proof. The SEM is defined as a Trotter splitting of an order two Velocity-Verlet integration of the conservative part and the simple order one Euler-Maruyama discretization of the fluctuation/dissipation part given by (2.35). The Euler-Maruyama scheme satisfies (1.34) at the order ω = 1 (see the proof of Proposition 1.3.4). Let us prove that the Velocity-Verlet discretization also does and then apply Lemma 1.3.6 in order to obtain the result.

In the framework of Stochastic Differential Equations, the generators of the time evolution of the positions and momenta of the Hamiltonian dynamics, denoted respectively by A and B, satisfy

A = M -1 p • ∇ q , B = -∇ q U(q) • ∇ p .
Here, A is the generator of the left dynamics of (1.7) and B the generator of the right one. Fix ψ ∈ D and (q, p) ∈ X q,p . The evolution operators of each subdynamics read

       e ∆tA ψ(q, p) = ψ • ϕ A ∆t (q, p) = ψ(q + ∆tM -1 p, p), e ∆tB ψ(q, p) = ψ • ϕ B ∆t (q, p) = ψ(q, p -∆t∇ q U(q)),
where ϕ A t and ϕ B t are respectively the flows of the dynamics generated by A and B. Consider the discretization given by (1.12) and denotes its evolution operator by P VV ∆t . We have that P VV ∆t ψ(q, p) = ψ • Φ VV ∆t (q, p) where Φ VV ∆t is the discretization procedure of the Velocity-Verlet scheme, given by (1.12). We clearly see that

P VV ∆t ψ = ψ • ϕ B ∆t/2 • ϕ A ∆t • ϕ B ∆t/2 = e ∆tB/2
e ∆tA e ∆tB/2 ψ.

The symmetric BCH formula (see Section 1.3.3) shows us that P VV ∆t satisfies (1.34) at the order ω = 2, which therefore means that it satisfies it at the order ω = 1 too.

Shardlow Splitting Algorithm (SSA)

SSA, in opposition to SEM, further decomposes the fluctuation/dissipation dynamics (2.27) into elementary pairwise fluctuation/dissipation dynamics involving only two particles i and j. In the DPD framework, these elementary dynamics write

                 dq i,t = dq j,t = 0, dp i,t = -γv ij,t χ(r ij,t ) dt + σ χ(r ij,t ) dW ij,t , dp j,t = γv ij,t χ(r ij,t ) dt -σ χ(r ij,t ) dW ij,t ,
(2.36) with γ and σ verifying (2.4).

The original DPD algorithm presented by Shardlow in [START_REF] Shardlow | Splitting for Dissipative Particle Dynamics[END_REF] integrates the elementary dynamics (2.36) using a BBK discretization [START_REF] Bruenger | Stochastic Boundary Conditions for Molecular Dynamics Simulations of ST2 Water[END_REF], consisting in one half-step of an explicit Euler discretization followed by one half-step of an implicit Euler-Maruyama discretization. In this particular case, the implicit updates can be rewritten in an explicit manner, thus yielding a fully explicit discretization of (2.36). However, as noted in Ref. [START_REF] Serrano | A Stochastic Trotter Integration Scheme for Dissipative Particle Dynamics[END_REF], the dynamics on the momenta can be analytically integrated. Indeed, by subtracting and adding the momenta variables of equation (2.36), we obtain

         d p i,t -p j,t = -2γv ij,t χ(r ij,t ) dt + 2σ χ(r ij,t ) dW ij,t ,
d p j,t + p j,t = 0.

(2.37)

In the monoatomic case, i.e m i = m for all i = 1, .., N , we can write

dp ij,t = -2 γ m χ(r ij,t )p ij,t dt + 2σ χ(r ij,t )dW ij,t ,
where p ij,t = p i,t -p j,t . We recognize the time evolution of an Ornstein-Uhlenbeck process. Therefore, following the procedure given in Section 1.4.2 and denoting a n ij,∆t = 2γχ(r n ij )/m, an analytical integration of (2.37) yields

           p n+1 i + p n+1 j = p n i + p n j , p n+1 ij = e -a n ij,∆t p n ij + 2σ m 4γ 1 -e -2a n ij,∆t G n ij .
From the above integration, we can recover the expressions of the updated momenta by adding and substracting p n+1 ij and p n+1 i

+ p n+1 j with                p n+1 i = p n i + 1 2 e -a n ij,∆t -1 p n ij + σ m γ 1 -e -2a n ij,∆t G n ij , p n+1 j = p n j - 1 2 e -a n ij,∆t -1 p n ij + σ m γ 1 -e -2a n ij,∆t G n ij .
The above computation extends to systems with particles of different masses by defining

a n ij,∆t = γµ ij χ(r n ij )∆t (where µ ij is defined in (2.10)),                      p n+1 i = p n i + e -a n ij,∆t -1 v n ij µ ij + σ 1 -e -2a n ij,∆t 2γµ ij G n ij , p n+1 j = p n j -e -a n ij,∆t -1 v n ij µ ij -σ 1 -e -2a n ij,∆t 2γµ ij G n ij .
(2.38) Indeed, in the polyatomic case, we notice that

dv ij,t = -γµ ij χ(r ij,t )v ij,t dt + σµ ij dW ij,t .
Consider the following infinite weak order discretization of v ij,t :

v n+1 ij = e -a n ij,∆t v n ij + σ µ ij 1 -e -2a n ij,∆t 2γ G n ij .
We notice that we can rewrite p i,t as

p i,t = (p i,t + p j,t ) + m j (v i,t -v j,t ) 1 + m j m i .
We still have d(p i,t + p j,t ) = 0. Therefore,

p n+1 i = (p n+1 i + p n+1 j ) + m j v n+1 ij 1 + m j m i , = (p n i + p n j ) + m j           e -a n ij,∆t v n ij + σ µ ij 1-e -2a n ij,∆t 2γ G n ij           1 + m j m i , = (p n i + p n j ) + m j v n ij 1 + m j m i + m j 1 + m j m i           e -a n ij,∆t -1 v n ij + σ µ ij 1 -e -2a n ij,∆t 2γ G n ij           .
Noticing that m j /(1 + m j /m i ) = 1/µ ij yields

p n+1 i = p n i + e -a n ij,∆t -1 v n ij µ ij + σ 1 -e -2a n ij,∆t 2γµ ij G n ij .
Similar computations for p n+1 j allow us to recover the result of (2.38). Let us now turn to the integration of DPDE. Taking advantage of the energy and momentum conservations, the DPDE elementary fluctuation/dissipation equations can be written as

                           dp i,t = -γ ij,t v ij,t χ(r ij,t ) dt + σ χ(r ij,t ) dW ij,t , dp j,t = -dp i,t , dε i,t = - 1 2 d        p 2 i,t 2m i + p 2 j,t 2m j        , dε j,t = dε i,t .
(2.39)

The extension of SSA to the DPDE framework [START_REF] Stoltz | A Reduced Model for Shock and Detonation Waves. Part I. The Inert Case[END_REF] updates the momenta first with a BBK discretization at fixed friction γ n ij and then updates internal energies in order to ensure the overall energy conservation. Instead, we use in this thesis the analytical integration (2.38) with γ = γ n ij instead. The integration procedure therefore consists in two steps:

1. We first update the momenta of the dynamics given by (2.39) by fixing the internal energies. The resulting equation is the time evolution of an Ornstein-Uhlenbeck process, and is given by (2.38) where we use γ n ij = γ(ε n i , ε n j ) instead of γ. The integration therefore reads

                     p n+1 i = p n i + e -a n ij,∆t -1 v n ij µ ij + σ 1 -e -2a n ij,∆t 2γ n ij µ ij G n ij , p n+1 j = p n j -e -a n ij,∆t -1 v n ij µ ij -σ 1 -e -2a n ij,∆t 2γ n ij µ ij G n ij ,
where

a n ij,∆t = γ n ij µ ij χ(r n ij )∆t.
2. We then update the internal energy in order to ensure the energy conservation. We introduce the kinetic energy variation of the momenta by

∆K n ij = (p n+1 i ) 2 m i + (p n+1 j ) 2 m j - (p n i ) 2 m i - (p n j ) 2 m j ,
and then update the internal energies by

           ε n+1 i = ε n i - 1 2 ∆K n ij , ε n+1 j = ε n j - 1 2 ∆K n ij .
Notice that the above procedure is the discrete counterpart of the Itô computation leading to (2.7).

Combining the two steps gives us the following discretization of (2.39):

P ij ∆t :                                        p n+1 i = p n i + e -a n ij,∆t -1 v n ij µ ij + σ 1 -e -2a n ij,∆t 2γ n ij µ ij G n ij , p n+1 j = p n j -e -a n ij,∆t -1 v n ij µ ij -σ 1 -e -2a n ij,∆t 2γ n ij µ ij G n ij , ε n+1 i = ε n i - 1 2 ∆K n ij , ε n+1 j = ε n j - 1 2 ∆K n ij .
(2.40)

We notice that, by construction, the above discretization exactly preserves the total energy.

The global SSA integration is obtained by a Trotter splitting of the DPDE dynamics into an Hamiltonian evolution and elementary subdynamics given by (2.39) corresponding to each particle couple, where each elementary dynamics is integrated using (2.40). Using the notation defined in Section 1.2.2, this means satisfies (1.34) at the order ω = 1 in order to obtain the result. In view of the symmetries dp i = -dp j and dε i = dε j , it suffices to expand the variables related to particle i. By expanding p n+1 i at order 1 in ∆t, we obtain

P SSA ∆t = P VV ∆t P SSA,fd ∆t , ( 2
p n+1 i = p n i + 1 2         e -a n ij,∆t -1 v n ij µ ij + σ 1 -e 2a n ij,∆t 2γ ij µ ij G n ij         , = p n i + σ χ(r n ij )G n ij √ ∆t -γ n ij v n ij χ(r n ij )∆t + C n ij ∆t 3/2 + O(∆t 2 ) ∆p n i , with C n ij = -1 2 σγ n ij µ ij χ(r n ij )G n ij . For ε n+1 , ε n+1 i = ε n i - 1 4 (p n+1 i ) 2 m i + (p n+1 j ) 2 m j - (p n i ) 2 m i - (p n j ) 2 m j , = ε n i - 1 4 2∆p n i v n ij + µ ij (∆p n i ) 2 , = ε n i - 1 2 σv n ij χ(r n ij )G n ij √ ∆t +        γ n ij (v n ij ) 2 2 - σ 2 µ ij 4 (G n ij ) 2        χ(r n ij )∆t + D n ij ∆t 3/2 + O(∆t 2 ), with D n ij = -1 2 C n ij v n ij -σγ n ij µ ij v n ij χ(r n ij ) 3/2 G n ij . Let Φ EM ∆t (x, G) = Φ EM,p
∆t (p, G), Φ EM,ε ∆t (ε, G) be the result of an Euler-Maruyama discretization of the fluctuation/dissipation dynamics (2.27) (i.e. the scheme 2.35). The previous computations show that SSA can be seen as a perturbation of the Euler-Maruyama scheme:

       p n+1 i = Φ EM,p ∆t (p n , G n ) i + C n ij ∆t 3/2 + O(∆t 2 ), ε n+1 i = Φ EM,ε ∆t (ε n , G n ) i + D n ij ∆t 3/2 + O(∆t 2
). The averages of both C n ij and D n ij are equal to zero, and the above expansions agree at order 1 with the result of a Euler-Maruyama discretization. Using similar procedures than for the proof of Proposition 2.2.2, we obtain that P SSA,ij ∆t satisfies (1.34) with L ij instead of L for ω = 1.

Remark 2.2.1. When linear scaling techniques are used (such as decomposing the system into cells of size r cut ), the order in which the particles are updated might change from one iteration to the other. From the above proof, we see that Proposition 2.2.4 still holds for any permutation η q of the couples (i, j). This means that for any permutation η q of the particle couples,

∀ϕ ∈ D,         P VV ∆t 0≤i<j≤N P η q (i,j) ∆t         ϕ = ϕ + ∆tLϕ + ∆t 2 r ∆t,ϕ ,
where r ∆t,ϕ is uniformly bounded for small ∆t in the sense of Definition 2.2.1. However, the change from iteration to iteration in the update order might depend on the history of the positions when linear scaling techniques are used, thus removing the Markovian property of the numerical solutions computed by SSA. Similar considerations apply when using Verlet lists. Special care must therefore be taken in order to keep the Markovian property of the SSA chains when using such techniques.

Shardlow Splitting Algorithm parallelization issue

As we saw previously, the SSA discretization of the fluctuation/dissipation dynamics consists in integrating sequentially each elementary pairwise dynamics given by (2.39). Therefore, the coordinates of each particle couple have to be updated after each elementary integration. This causes the processors to constantly modify the coordinates of the particles.

When domain decomposition techniques are used for parallel simulations, some particles interact with others allocated to different processors. Therefore, during the corresponding pairwise updates, one must be careful to prevent the processors to access the coordinates of the same particle simultaneously. This causes many delays due to idle processors waiting for memory access on a particle dealt with by another processor. In addition, when integrating a given elementary dynamics, the information on each particle of the couple must be up-to-date, thus requiring data exchange between processors. This drastically increases the quantity of information exchange compared to standard parallelizable schemes, thus slowing down the simulation. Figure 2.4 illustrates this fact by considering the particles inside the domain allocated to one processor.

A way to parallelize SSA has been presented in 2014 by Larentzos and al. in [LBM + 14]. The technique used to avoid additional delays caused by the sequentiality of the integration is to decompose the integration into 26 steps in 3D to allow sequentially for updates between a cell and its 26 neighbors, the 26 steps being determined in a manner that processors never updates couples with a common particle.

Let us explain this procedure in a two-dimensional case. Figure 2.5 displays the four steps of the integration procedure. Each box represents a space domain allocated to one single processor. The view is centered around one central domain, enclosed in continuous lines, the zone highlighted in red defining the region of space where the central processor integrates the elementary fluctuation/dissipation interactions. The same procedure is performed simultaneously by all the processors. As we can see, the integration is decomposed in four steps:

1. Stage 1: integration of the interactions where both particles are inside the same domain.

2. Stage 2: integration of the interactions where one of the particles is in the upper neighbor domain of the processor. resents the space domain allocated to a given processor, the exterior region being allocated to other processors. The dashed boxes of length r cut represent the cells of the linear scaling techniques: a particle on a given cell can only interact with particles in the neighboring cells. The dashed particles ("ghost" particles) are particles allocated to other processors. The particles at a distance r > r cut of the frontier ("internal" particles, in red ) never interact with ghost particles: the sequentiality of the integration does not affect the computation speed because no additional idle time nor information transfer between processors are needed. The particles at a distance r < r cut of the frontier ("frontier" particles, in blue) interact with ghost particles: the processor eventually has to wait before accessing any ghost particle's data, or wait for another processors to access its frontier particle's data, thus generating delays.

3. Stage 3: integration of the interactions where one of the particles is in the right neighbor domain of the processor.

4. Stage 4: integration of the interactions where one of the particles is whether in the top-right or bottom-right neighbor domain of the processor.

Of course, between each stage, the updated data of each particle must be sent to its respective processor. We can see that, if all processors follow the above procedure, they all consider interactions in distinct regions of space at a given time, thus forbidding simultaneous access to particle data and therefore suppressing the SSA parallelization issue. Using the symmetry of the integration, the above procedures ensures that all the elementary fluctuation/dissipation interactions are integrated. It was shown that this procedure allows for a scaling efficiency in the number of processors comparable to the scaling of other reference straightforwardly parallelizable schemes like SVV [LBM + 14].

However, if the above procedure allows for efficient parallel DPD/DPDE integration using SSA, it is valid only for MPI parallelization. If an architecture like OpenMP (www.openmp.org) is used or another architecture where several processors have access to some shared memory, this procedure no longer applies, and the same issues arise than for the sequential case. In addition, Larentzos and al. parallelization procedure requires work that is not trivial, and a simpler solution for DPDE parallel accurate integration would be welcome.

New parallelizable schemes for the integration of DPDE

As seen from Section 2.2, up to our knowledge, one should choose between difficultly parallelizable schemes enjoying good energy preservation properties or straightforwardly parallelizable schemes that preserve very poorly the energy. In this section, we consider two new techniques for integrating the DPDE that are both straightforwardly parallelizable.

The two new parallelizable schemes we propose in this section also rely on the splitting between the Hamiltonian and the fluctuation/dissipation parts mentioned at the beginning of Section 2.2, but with new strategies to discretize the fluctuation/dissipation part:

• The first scheme we present is called "Splitting with Energy Reinjection" (termed SER in the sequel). Its discretization of the fluctuation/dissipation part (2.27) is similar to SEM but uses a global symmetric reinjection of the kinetic energy variation into the internal energies, instead of directly discretizing the dynamics of the internal energies. This allows to automatically preserve the total energy during the fluctuation/dissipation integration.

• The second scheme is a mix between SSA and SER, and is therefore termed "Hybrid" in the sequel. As parallel simulations are performed using a spatial repartition of the simulation box between the processors, the bottleneck for the parallelization of SSA arises from particles located in different domains [LBM + 14]. Therefore, the idea of the Hybrid scheme is to integrate the elementary fluctuation/dissipations interactions involving particles located on the same processor by a pass of the SSA algorithm, while the remaining interactions are taken care of by a SER discretization.

Splitting with Energy Reinjection (SER)

The SER integration of the fluctuation/dissipation (2.27) is performed in two steps. First, momenta are integrated using a simple Euler-Maruyama discretization as

p n+1 i = p n i + δp n i , δp n i = N j=1,j i δp n ij , (2.43)
with δp n ij = -γ n ij χ(r n ij )v n ij ∆t +σ χ(r n ij )G n ij √ ∆t.
The internal energies ε i are then updated in order to compensate for the energy variation during the update of the momenta. In order to implement this idea, we need to identify in the global kinetic energy variation of each particle the contribution of every single pairwise interaction, in order to redistribute the associated elementary energy variations. In fact, a simple computation shows that

∆K n i = (p n+1 i ) 2 2m i - (p n i ) 2 2m i , = δp n i 2m i 2p n i + δp n i , = N j=1,j i δp n ij • v n i + δp n i 2m i ∆ j K n i .
The term ∆ j K n i represents the contribution of particle j to the kinetic energy variation of particle i. The internal energies are then updated by reinjecting the elementary variations ∆ j K n i in a symmetric manner:

∆ε n i = - 1 2 N j=1,j i ∆ i K n j + ∆ j K n i , = - 1 2 N j=1,j i δp n ij •       v n ij + 1 2       δp n i m i - δp n j m j             .
(2.44)

The discretization of (2.39) can therefore be summarized as From a pratical viewpoint, the presence of both δp n ij and δp n i in (2.44) requires two sweeps on the system, a first one to compute δp n ij and sum it into δp n i , and a second one to compute the products δp n ij • δp n i . This seemingly requires to store all increments δp n ij (or at least the Gaussian increments G n ij ), which is somewhat prohibitive. This can however be avoided if the Gaussian increments G n ij can be exactly regenerated to the value they had on the first pass of the algorithm when the particle pairs are revisited. The SARU pseudo-random generator [START_REF] Afshar | Exploiting seeding of random number generators for efficient domain decomposition parallelization of dissipative particle dynamics[END_REF] for instance allows such an easy recomputation. The pseudo-number generator SARU has been designed for improving computational cost for parallel DPD computations by saving the cost of transfering the random numbers between generators. It takes three integer seeds (i 1 , i 2 , i 3 ) ∈ N 3 , with i 1 < i 2 , to generate a pseudo-random integer between -2 31 and 2 31 . However, the interesting fact of SARU for DPD calculations is that this number is completely defined by the three seeds (i 1 , i 2 , i 3 ). Therefore, the random numbers G ij do not have to be transfered between processors, thus saving computational time.

P SER,fd ∆t :                      p n+1 i = p n i - N j=1,j i γ n ij χ(r n ij )v n ij ∆t + N j=1,j i σ χ(r n ij )G n ij √ ∆t, ε n+1 i = ε n i - N j=1,j i δp n ij •       v n ij + 1 2       δp n i m i - δp n j m j             . ( 2
The global SER discretization is obtained with a Trotter splitting of DPDE. Therefore, it writes P SER ∆t = P VV ∆t P SER,fd ∆t .

Proposition 2.3.1. The SER discretization of the DPDE dynamics is of weak order one. In addition, its evolution operator satisfies

∀ϕ ∈ D, P SER ∆t ϕ = e ∆tL ϕ + ∆t 2 r ∆t,ϕ ,
where r ∆t,ϕ stays uniformly bounded for small ∆t in the sense of Definition 2.2.1.

Proof. By the same reasoning as in the proof of weak order one of the previous schemes, it suffices to prove that P SER,fd ∆t satisfies (1.34) with ω = 1 in order to obtain the result. Let us prove that the evolution operator of the SER discretization given by (2.45), denoted by P SER,fd ∆t , of the fluctuation/dissipation dynamics can be expanded as ∀ϕ ∈ D, P SER,fd ∆t ϕ = ϕ + ∆tL fd ϕ + ∆t 2 r ∆t,ϕ , where r ∆t,ϕ is uniformly bounded when ∆t goes to zero.

As reported in Section 2.3.1, the momenta are updated with an Euler-Maruyama discretization writing

p n+1 i = p n i + N j=1,j i δp n ij , with δp n ij = -γ n ij χ(r n ij )v n ij ∆t + σ χ(r n ij )G n ij √ ∆t.
A simple computation shows that (2.45) implies

ε n+1 i = ε n i - 1 2 N j=1 j i δε n ij + O(∆t 2 ), (2.46)
where

δε n ij = -γ n ij (v n ij ) 2 χ(r n ij )∆t + σv n ij χ(r n ij )G n ij √ ∆t + η n ij ∆t + C n ij ∆t 3/2 .
In this expression, we introduced

                 η n ij = 1 2 R n ij •       R n i m i - R n j m j       , C n ij = - 1 2       F n ij •       R n i m i - R n j m j       + R n ij •       F n i m i - F n j m j             , with                          F n ij = γ n ij v n ij χ(r n ij ), F n i = N j=1 j i F n ij , R n ij = σ χ(r n ij )G n ij , R n i = N j=1 j i R n ij .
The random variables C n ij involve only terms containing one Gaussian random variable. Therefore, they vanish in average: E C n ij = 0. The random variables η n ij involve terms containing products of Gaussian variables, but their averages can be easily computed as

E[G n ij G n kl ] = d δ ik δ jl -δ il δ jk . A simple computation then shows that E η n ij = d µ ij σ 2
2 χ(r n ij ). The previous computations show that SER can be seen as a perturbation of the Euler-Maruyama scheme, i.e,

                 p n+1 i = Φ EM,p ∆t (p n , G n ) i , ε n+1 i = Φ EM,ε ∆t (ε n , G n ) i - 1 2 j=1 j i       η n ij -d µ ij σ 2 2 χ(r n ij )       ∆t + C n ij ∆t 3/2 + O(∆t 2 ).
where we have used the same notation than in the proof of Proposition 2.2.4. By a similar reasoning as in the proof of Proposition 2.2.2, this implies that P SER,fd ∆t satisfies the above expansion in powers of ∆t, which gives us the result by the successive applications of Lemma 1.3.6 and Theorem 1.3.3.

We have presented a version of SER where the momenta integration is done according to an Euler-Maruyama discretization, of weak order one. However, we see that the energy reinjection procedure still applies with any momenta discretization than can be written under the form

p n+1 i = p n i + N j=1,j i δp n ij .
This could be used to devise higher-order schemes.

A scheme mixing SSA and SER: Hybrid

The Hybrid scheme can be seen as a blending of SER and SSA, where the elementary fluctuation/dissipation interactions involving particles on the same processor are integrated by a pass of the SSA algorithm, while the remaining interactions are integrated with the SER scheme.

Let us describe more precisely this algorithm. By using the notation of Section 1.2.2, we denote by P Hybrid ∆t the Hybrid discrete generator and P SSA,ij ∆t the generator of the SSA discretization of the elementary dynamics given by (2.40). Note that only the components p i , p j , ε i , ε j of (q, p, ε) are changed by P SSA,ij ∆t . Let us denote by A q the set of particle couples where both particles are in the same processor, and A c q the remaining couples. The elements of A q are denoted by (i 1 , j 1 ), (i 2 , j 2 ), . . ., (i l , j l ), where l is the number of elements in A q . Finally, we denote by P SER,A c q ∆t the discrete generator of the SER discretization of all the interactions involving the particle couples in A c q . The expression of P (2.48)

Note that A q depends on the positions q and therefore might change from one iteration to the other. Hybrid is thus subject to the same considerations regarding the update order than SSA (see the discussion in Remark 2.2.1). Therefore, in order to prove the weak consistency of Hybrid, we have to assume that the update order in P SSA,A q ∆t does not depend on the position history. Proposition 2.3.2. Assume that the update order in P SSA,A q ∆t is defined by a systematic procedure that does not depend on the position history. Then the Hybrid discretization given by (2.48) is of weak order one. In addition, its evolution operator satisfies where r ∆t,ϕ is uniformly bounded for small ∆t in the sense of Definition 2.2.1.

Proof. Denote by L A c q the generator of the fluctuation/dissipation DPDE dynamics associated with the particles indexed by the elements of A c q . By the same reasoning as in the proof of weak order one of the previous schemes, it suffices to prove that P SER,A c q ∆t satisfies (1.34) with ω = 1 and L A c q used instead of L in order to obtain the result. Fixing q ∈ Q, the extension of Proposition 2.3.1 to the SER fluctuation/dissipation discretization of the dynamics generated by L A c q is straightforward because of the pairwise structure of the friction and random interactions. Therefore, for any ϕ ∈ D, P SER,A c q ∆t ϕ = e ∆tL A c q ϕ + O(∆t 2 ). This gives the result.

Energy conservation properties of the schemes

We focus in the numerical tests we present in this section on the study of the conservation of the energy by all the schemes presented in Sections 2.2 and 2.3. We start by studying their energy conservation, and quantify the possible drifts that occur. Then, we study the drifts of each one of the three components of the energy, namely the kinetic, potential and internal energies. Each one of the components of the energy involves only one type of the three categories of degrees of freedom of the system: the positions q i , the momenta p i and the internal energies ε i . This therefore allows to assess the quality of the sampling for each degree of freedom.

We consider N = 343 identical particles, i.e m i = m for i = 1, .., N , and we set d = 3. Initial conditions are obtained by melting a simple cubic crystal of density ρ = 0.5787 over a time t = 100, using a SSA scheme with a time step ∆t = 0.005. Such thermodynamical conditions represents a liquid Lennard-Jones fluid (see Figure 3.1 for a plot of the computed radial distribution function g(r) of the system). The results are given in reduced

Drifts of the total energy

In order to precisely quantify the possible energy drifts and determine in particular the influence of the time step ∆t on the drift rate, we compute the average evolution of the energy in time by performing N traj = 10 4 trajectories of time t f = 10 each. The initial configurations of each trajectories are sampled according to the canonical measure µ β,C v given by (2.15), obtained by the procedure described in Section 2.1.2. The following procedure is therefore iterated for m = 1, .., N traj :

1. We start from an initial configuration xm,0 = ( qm,0 , pm,0 , εm,0 ).

2. The internal degrees of freedom are sampled according to the DPDE internal canonical measure (2.22) by integrating the corresponding Overdamped Langevin dynamics given by (2.23) for t = 1 and with ∆t = 0.001, using an Euler-Maruyama discretization defined in Section 1.3.1. The resulting configuration is denoted by εm,0 .

3. The external degrees of freedom are sampled according to the canonical measure (1.19) by integrating the Langevin dynamics (1.54) for t = 1 and with ∆t = 0.001, using the GLA algorithm of Section 1.4.2. The resulting configuration is denoted by ( qm,0 , pm,0 ).

4. We define xm,0 = ( qm,0 , pm,0 , εm,0 ) as the DPDE configuration composed of the results of the two previous steps.

5. The internal and external degrees of freedom are then equilibrated for t = 1 with an integration using the same scheme and time step used for the simulation. The resulting configuration, denoted by x m,0 is the initial configuration of the subsequent trajectory over which averages are computed.

6. Once the trajectory is computed, we define xm+1,0 = xm,0 . Remark 3.1.1. We notice that the above procedure is not exactly the procedure described in Section 2.1.2 because the Langevin equilibration is not projected in order to ensure the conservation of the total momentum. Therefore, this procedure yields initial conditions that do not have the same total momenta. We sample in fact the measure of density e -βE+ S k B obtained when µ β,C v is not conditioned to δ P 0 . This is not problematic however because we are interested in this section in the average drift of the total energy and of its components.

We denote by ϕ t the expectation of ϕ(x t ) over all realization of the DPDE dynamics and all the initial conditions x 0 ∼ µ β,C v for any observable ϕ ∈ D, i.e

ϕ t = E µ β,C v [ϕ(x t )] . (3.1)
In practice, ϕ n∆t is approximated by

ϕ n ∆t,N traj = 1 N traj N traj m=1 ϕ (x m,n ) , ( 3.2)
where x m,n is the n-th configuration of the m-th trajectory. We denote by ϕ n ∆t the limit of ϕ n ∆t,N traj when N traj goes to infinity, i.e Figure 3.3 shows the behavior of the time-dependent energy drift averaged over trajectories as a function of time, for various schemes with ∆t = 0.006. Note that no scheme manages to keep the total energy constant since linear drifts in time are observed in all cases, SVV being substantially worse than all the other schemes in terms of energy conservation. Figure 3.4 represents the average total energy on equivalent two-dimensional simulations, with N = 100 and ∆t = 0.006, however using SEM instead of SVV as the reference scheme. We notice that SEM, similarly to SVV, is substantially worse than the other schemes. We therefore conclude that both SVV and SEM exhibit very bad energy conservation properties.

ϕ n ∆t = lim N traj →∞ 1 N traj N traj m=1 ϕ (x m,n ) . t = 1 t = t f DPDE, ∆t DPDE, ∆t x0
Energy drifts such as those displayed in Figure 3.3 and Figure 3.4 destroy all hope of sampling the microcanonical DPDE measure µ E 0 ,P 0 nor the DPDE canonical measure µ β,C v : a stable energy is required when computing thermodynamical averages in energypreserving statistical ensembles.

Quantification of the total energy drifts

In order to decide whether the drift is acceptable for thermodynamical averages estimations, we quantify in this section the rate of variation of the energy as a function of the time step ∆t. energies as

E n ∆t = E 0 (1 + α ∆t n∆t) , ( 3.3)
which allows to identify the relative energy drift rate α ∆t . Figure 3.5 presents the evolution α ∆t as a function of the time step. Unsurprisingly, the drift rates for SVV are substantially larger in absolute value than those of the other two schemes. The energy drift rates of SER are of the same order of magnitude than those of SSA for small time steps, but they increase much faster for larger ∆t. The SER scheme therefore seems to be less stable than SSA. In addition, we notice that the drift rates of SSA and SER schemes are both greater than 1 × 10 -6 for time steps ∆t ≥ 0.007. This means that, for simulation times longer than t max = 10 4 , the energy variation would be greater than 1%, which might be considered as an upper stability limit in terms of energy conservation. Such simulation times may be long enough to estimate equilibrium properties, but they may prove somewhat short to estimate transport properties.

Moreover, the results of Figure 3.5 suggest that the drift rate α ∆t has a polynomial behavior with respect to ∆t :

α ∆t = K∆t ω .
A least-square fit in a log-log diagram gives ω ≃ 2 for SVV, and ω ≃ 4 for the other schemes. This fast increase of the drift rate places a severe limitation on the use of larger time steps in DPDE simulations.

Drifts of the components of the energy

We carefully studied in Section 3.1.1 the total energy drifts. Let us now consider the behavior of other thermodynamical observables, and consider the intensive observables corresponding to the three components of the energy involving the three categories of degrees of freedom of a DPDE system: the positions q, the momenta p and the internal energies ε. This therefore allows to assess the behavior of the energy for each degree of freedom.

All the observables considered in this section are normalized by their corresponding degrees of freedom. Dealing with such intensive observables allows to discard the size of the system considered when estimating average values (provided the system is large enough so finite size effects are negligible). The first component is the kinetic energy per degree of freedom, depending only on the momenta and given by

E kin = 1 2dN N i=1 p 2 i m i µ β,C v < (3.4)
Note that we have E kin = T /2, where T is the equilibrium temperature of the system. The second component is the potential energy per degree of freedom, depending only on the positions of the particles and defined as

E pot = 1 dN U µ β,C v , ( 3.5)
where • µ β,C v refers to averages with respect to the measure µ β,C v introduced in (2.15). The third and last quantity is the internal energy per degree of freedom, depending only on the particle internal energies ε and defined as

E int = 1 C v N N i=1 ε i µ β,C v . (3.6)
Note that, because the internal degrees of freedom of each DPDE particle is given by the heat capaciy C v , we divided E int by C v N instead of dN .

The drift rates of the estimators given by (3.5), (3.4) and (3.6) as a function of the time step are plotted in Figure 3.6. As expected, the drift rates for all components are large for SVV and small for SSA. SER on the other hand exhibits a non-trivial behavior: its drift rates are quite large for individual components of the energy, although these large drifts compensate each other in order for the total energy to drift only slowly. In fact, we even observed in some situations that the drift of some energy components was not linear as a function of time. For instance, we noticed that the SER average components of the energy are not linear functions of time, as seen in Figure 3.7, in opposition to those of SSA or SVV. The very notion of drift rate for SER is therefore dubious.

Error estimations of equilibrium observables

We have seen in Section 3.1 that no scheme could manage to keep the total energy constant, and that linear drifts of the total energy E t occurred. The amplitude of these drifts depend on the scheme, and we have seen that SER's drifts were of the same order than those of SSA, and far better than those of SVV. However, these drifts require the use of very small time steps in order to have a stable energy when computing thermodynamical averages, and additionally make the thermodynamical averages dependent on the integration time as seen in Figure 3.6. A procedure ensuring the energy conservation is therefore devised in Section 3.2.1, in order to obtain a well-defined equilibrium state, followed by a discussion about the precautions needed when using such a projection and some expectations on the numerical results. We then compute in Section 3.2.2 estimations of the three components of the energy defined in Section 3.1.2. Finally, we perform in Section 3.2.3 numerical tests ensuring that the Hybrid behavior is consistent with the spatial decomposition of the simulation box when parallel simulations are performed.

Projection in order to ensure the energy preservation

Previous to any equilibrium simulation aimed at computing thermodynamical observables, one must first remove the energy drift. A natural approach, already suggested in Ref. [START_REF] Lisal | Dissipative Particle Dynamics at Isothermal, Isobaric, Isoenergetic, and Isoenthalpic Conditions using Shardlow-like Splitting Algorithms[END_REF], is to enforce the energy conservation by an appropriate projection on a surface of constant energy. There are several possible projection procedures.

One possible option is to use a Lagrange multiplier, i.e to replace a configuration x = (q, p, ε) such that E(x) E 0 by x -λ∇E(x), with λ chosen such that E q -λ∇ q E(x), p -λ∇ p E(x), ε -λ∇ ε E(x) = E 0 , or more explicitly

E q -λ∇U(q), p - λ m p, ε -λ 1 = E 0 , (3.7)
where 1 is the N -dimensional vector whose components are all equal to 1. However, the numerical computation of the parameter λ satisfying (3.7) is a computationally expensive task typically requiring several iterations of a Newton procedure, and hence several evaluations of the energy per time step. The potential energy U(q) being expensive to compute in terms of CPU cost, this procedure appears to be prohibitive. Note also that the total momentum is not preserved, so that this extra conservation law should be subsequently enforced in an appropriate manner.

It is therefore much more convenient from a practical viewpoint to only play on the internal energies to adjust the total energy -as was also already suggested in Ref. [START_REF] Lisal | Dissipative Particle Dynamics at Isothermal, Isobaric, Isoenergetic, and Isoenthalpic Conditions using Shardlow-like Splitting Algorithms[END_REF]. More precisely, at the end of one iteration of the numerical scheme, the resulting internal energies are replaced by ε n+1 i -(E n -E 0 )/N . Considering a given discretization Φ ∆t of the DPDE dynamics, the corresponding projected (or corrected) scheme therefore consists in the following:

• We start by a configuration x n ∈ X , and we assume it is of energy E(x n ) = E 0 .

• We obtain xn = Φ ∆t (x n ) by the considered DPDE discretization. Note that E( xn ) E 0 as we saw in Section 3.1.

• We rescale the internal energies proportionnally to the energy variation, i.e

ε n+1 = εn - 1 N (E( xn ) -E 0 ) 1 N ,
and we define (q n+1 , p n+1 ) = ( qn , pn ) (i.e we do not change the positions and momenta from the result of the discretization).

Note that this allows to remove energy variations due to the discretization of both the Hamiltonian and fluctuation/dissipation parts. Enforcing the total energy conservation allows us to obtain a well-defined steady state, for which average properties can be safely computed.

In practice however, one has to be careful that internal energies should remain positive. Indeed, we saw in Section 2.1.2 that the internal thermodynamics of DPDE particles is not well-defined whenever ε ≤ 0. Nevertheless, many cases can lead to a non-negligible probability that the result of the discretization at iteration n ∈ N yields a particle having a negative internal energy ε i ≤ 0:

• If the time step ∆t is large, biases are therefore also large, which leads to non-stable behaviors. This might cause negative internal energies.

• The parameter σ controls the strength of the fluctuation/dissipation dynamics compared to the conservative part. Therefore, when it is large, the oscillations of the internal energies are also large, thus increasing the probability of one particle to have a negative internal energy.

• The equilibrium distributions of the internal energies also depends on C v , as the definition (2.22) of ν β,C v indicates. The larger the C v is the smaller the fluctuations of the internal energies are proportionally to their average value, therefore reducing the chances of a particle to have a negative internal energy.

• Finally, we saw in Section 2.1.2 that the internal and external temperatures are equal at equilibrium. Considering the DPDE internal equation of state (2.12), we see that particles with low internal temperatures also have low internal energies.

Low temperature systems thus have a higher probability of yielding negative internal energies than high temperature systems.

Procedures allowing to avoid such problems have not been considered in the simulations presented in this thesis, but are discussed in the "Perspectives" chapter at the end. Finally, one should expect the results of the simulations to yield errors on the invariant measure of order one, as proved for each scheme in Chapter 2. However, we can specificate this by saying that the errors should in fact be proportional to σ 2 ∆t. Indeed, we have seen in Section 2.1.2 that the each elementary fluctuation/dissipation generator writes

L ij = -γ ij v ij • A ij + σ 2 2 A 2 ij , with A ij = ∇ p i -∇ p j - v ij 2 ∂ ε i + ∂ ε j .
Therefore, because γ ij and σ satisfy (2.16), we have

L ij = σ 2 Lij , with Lij = 1 4 1 T i + 1 T j v ij • A ij + 1 2 A 2 ij .
We see that Lij does not involve any term in σ. Therefore, the elementary DPDE fluctuation/dissipation evolution operator can be rewritten as e ∆tL ij = e σ 2 ∆t Lij , and everything behaves as if σ 2 ∆t is the effective time step. However, the interaction between the conservative integration and the projection could alter this dependence.

Biases on the average equilibrium properties for projected schemes

In this section, different simulations are carried out using the SEM, SSA, SER and Hybrid schemes presented in Section 2.2 and 2.3, corrected by the projection procedure on the internal energies. Sequential simulations are first performed, followed by parallel simulations. The Hybrid scheme can only be tested in parallel simulations, otherwise it reduces to SSA. Therefore, it is considered only for the parallel simulations.

In this section, we consider systems of higher density ρ = 0.8095, corresponding to liquid Lennard-Jones fluid (see Figure 3.8 for the plot of the computed radial distribution function of the system). The results are still given in reduced units.

Sequential simulations

The results presented in this section are obtained by computing averages over one long trajectory of physical time t f = 1000, with the heat capacity of the particles C v = 50 and the same initial condition obtained by the procedure detailed in Section 3.1. Note that we do no longer use the procedure described by Figure 3.2 to obtain several initial conditions because we only perform one trajectory per simulation. The cut-off radius is set to r cut = 3. value as ∆t → 0 for all components of the energy. The second point is that the systematic biases depend on the time step ∆t as expected, but also on the value of σ: larger values of σ lead to much larger errors. As predicted in Section 3.2.1, the biases seem to be proportional to σ 2 ∆t (see Figure 3.11). The third point is that, in opposition to the biases of SSA which remain very small even for large values of ∆t and σ, those of SEM and SER are much larger. For instance, the bias on the estimation of the kinetic energy for SER and SEM at ∆t = 0.008 and σ = 4 is around 30% of the extrapolated reference value.

Parallel simulations

Thermodynamics averages in this section are computed over a single trajectory of total time t f = 450, with a lower heat capacity C v = 10. The cut-off radius is set to r cut = 2.5. A system of N = 148, 500 particles was decomposed on N proc = 8 processors. Figure 3.12 presents the energy estimations as a function of the time step for σ = 2, compared to reference values computed by a sequential SSA simulation of a system of N = 4000 particles. We see that, similarly to the sequential simulations results, SEM and SER estimations are very similar and that SSA estimations are significantly smaller than those of SER and SEM. The Hybrid estimations however have an accuracy similar to those of SSA. Simulation results for σ = 4 (not reported here), showed greater biases than for σ = 2 for all the schemes with Hybrid and SSA being significantly more accurate than SEM, consistently with the sequential simulations. However, SER was found to yield particles with negative internal energies with very high probability even for small time steps, thus forbidding to obtain converged estimates (i.e estimates with a negligible statistical error). This instability is probably due to the SER energy reinjection procedure leading to larger oscillations of the internal energies compared to SSA or SEM. The fact that such phenomena was not observed for the sequential simulations is related to the lower value of C v used for parallel simulations, as discussed in Section 3.2.1. CPU timings are reported in the second column of Table 3.1. For the comparison to be fair, all the timings have been measured on a single core. We also compare these timings to reference timings from a sequential SSA simulation (third column of Table 3.1). The last column of Table 3.1 gives the required CPU time of a simulation of the same system as in Figure 3.12, using a time step chosen in order to obtain a given accuracy on the estimation of the kinetic energy. For each scheme, this time step, denoted ∆t lim , is taken such that the bias of the kinetic energy estimation is equal to e lim , where e lim is defined as the largest bias of the Hybrid kinetic energy estimations of Figure 3.12. The biases are computed as e ∆t = | E kin ∆t -E kin ∆t=0 |, where E kin ∆t=0 is obtained by extrapolating the results of Figure 3.12 to ∆t = 0. The timings of the last column of Table 3.1, denoted τ E kin cpu , are then obtained by multiplying the timings of the second column of Table 3.1 with the number of iterations necessary to reach t f = 100 when using the time step ∆t lim . The choice of the kinetic energy is arbitrary (as is the choice of e lim ). However, timings related to other observables show qualitatively similar results. simulations with time steps such that the kinetic energy estimation has a given bias. This bias is chosen to be the largest bias of the Hybrid kinetic energy estimations of Figure 3.12.

Table 3.1 tells us that SER is much slower than SEM, and that using SER takes more than twice as long as using SEM in order to reach a given accuracy. This sub performance of SER is to be tempered however with the fact that we consider here all schemes to be corrected with the projection of Section 3.2.1. We have seen in Section 3.1 that, on their own, SEM requires time steps much smaller than SER in order to achieve the same accuracy. As for Hybrid, we see that its increased accuracy displayed in Figure 3.12 more than compensates the extra computational cost compared to a very simple scheme such as SEM, as it can be seen in the fourth column of Table 3.1.

Validation of the Hybrid scheme: influence of the parallelization

There are no other parameters than the time step for SSA, SER and SEM schemes. This is not the case for the Hybrid scheme, where the spatial distribution of the simulation box between processors affects the computation. It is thus necessary, in order to further validate the Hybrid scheme, to study the influence of the ratio N /N proc , where N proc is the number of processors used for the simulation. Indeed, the ratio N /N proc directly determines the fraction of interactions treated with SER or SSA.

In order to study this influence, Hybrid simulations are performed on N proc = 8 processors for systems in the same thermodynamic state as those of the parallel simulations of Section 3.2.2. The time steps of the simulations range from ∆t = 5 × 10 -4 to ∆t = 2.5 × 10 -3 and the number of particle from 40 × 10 × 10 to 160 × 10 × 10. Note that the parallelization is done according to the x axis, and we only change the number of particles along this axis. Hybrid averages are presented in Figure 3.13, together with the reference SSA and SER estimations of Figure 3.12. The N /N proc ratios according to the system sizes are displayed in Table 3 We notice in Figure 3.13 that, as expected, Hybrid estimations go from those of SER towards those of SSA as N /N proc increases. In the particular situation considered here, SER and SSA estimations have biases of opposite signs that compensate in the Hybrid scheme, thus leading to biases smaller than those of SSA. Averages of the internal energy exhibit similar behaviors as in Figure 3.13. SSA and Hybrid biases on the potential energy are too small in order to conclude anything.

4

Application of the DPDE to the simulation of shock waves

After having studied the scheme behavior on equilibrium systems, we now turn in this chapter to their behavior on systems outside equilibrium, and consider the specific case of shockwaves.

Shockwaves involve phenomenons of very small time and space scales and large energies at the shock front and much larger time and space scales and lower energies in the relaxation of the shocked materials. As an example, in the following simulations, an increase of approximately 1000K is observed in the local temperature of the particle, going from 1000K to 2000K in approximately 1.2ps, where the characteristical time for a Lennard-Jones Argon atom is 2.2ps. Therefore, coarse-graining methods preserving the energy are much needed.

In the general case, nonequilibrium properties are more difficult to compute than equilibrium properties because they need to sample trajectories instead of configurations. In addition, they often require to simulate larger systems than equilibrium properties. Indeed, equilibrium systems are well represented using periodic boundary conditions, thus allowing to reduce the size of the system, provided that finite size effects do not occur. On the other hand, shockwaves, for instance, must propagate from one side of the simulations box to the other, and some of their dynamical properties are computed by averaging properties over slabs normal to the shock propagation, or by averaging properties of the shock front over all the trajectory. Both methods need large systems that often cannot be simulated with a single processor.

As said in Chapter 2 and shown in Chapter 3, we already have an efficient and accurate scheme for computing average equilibrium properties with SSA. The main interest of the development of SER and Hybrid is therefore to be able to perform straightforwardly massively-parallel DPDE simulations, generally needed for the computations of dynamical properties, outside equilibrium or not. We thus verify in this chapter that the new schemes presented in Chapter 2 are able to correctly reproduce dynamical behavior. We compare their results with those of SEM, which is, with SVV, the only straightforwardly parallelizable DPDE integration scheme known previous to this work.

An important point in this chapter is that we no longer correct the schemes by the projection of Chapter 3.2.1 since we consider a nonequilibrium system inhomogeneous in space. This makes indeed the global energy reinjection dubious. The second important point is that all numerical simulations in this chapter are parallel simulations performed on hundreds of processors. Therefore, this chapter is also about demonstrating that the new schemes presented in Chapter 2 behave correctly on massively-parallel simulations.

Description of the simulations

The system we consider is composed of N = 450, 000 DPDE particles. In order for the numerical values to be easier to interpret from a physical viewpoint, we work in this chapter in physical units. The Lennard-Jones parameters in (1.2) are those of Argon: ε LJ = 1.65710 -21 J and r LJ = 3.40510 -10 m, while the masses of the particles are set to m = 6.634 × 10 -26 kg. The initial condition is obtained by melting a face centered cubic crystal composed of 250 × 15 × 15 = 112, 500 unit cells, at density ρ = 1228 kg.m -3 , with periodic boundary conditions. We use 50 × 2 × 2 = 200 cores for the simulations, and various time steps ranging from ∆t = 10 -15 s to ∆t = 4 × 10 -15 s (the value ∆t = 10 -15 s corresponds to ∆t * = 4.5 × 10 -4 in reduced units). A time step of ∆t = 10 -15 s is similar to the one used to integrate Hamiltonian dynamics of shock waves in Argon, but is already one order of magnitude larger than for full atom simulations of molecular systems. With the notation introduced in (2.21), we set γ = 10 -13 kg.s -1 (which corresponds to a fluctuation parameter σ = 5.25 × 10 -17 kg.m.s -3/2 or σ = 7.35 in reduced units). We also set a constant heat capacity C v = 1.38 × 10 -22 J.K -1 or C * v = 10 in reduced units. We equilibrate the system at a given temperature T ref = 1000 K by superimposing to the DPDE equations a Langevin thermostat on the momenta:

                                         dq i,t = p i,t m i dt, dp i,t = -∇ q i U(q t ) dt + N j=1,j i -γ ij,t χ(r ij,t )v ij,t dt + σ ij,t χ(r ij,t )dW ij,t , -γv i,t + 2 γ β dB i,t , dε i,t = 1 2         N j=1,j i        γ ij,t v 2 ij,t -d σ 2 ij,t 2 µ ij        χ(r ij,t ) dt -σ ij,t χ(r ij,t )v ij,t • dW ij,t         , ( 4.1)
with γ a scalar parameters verifying the Langevin fluctuation/dissipation (1.55), and (B i,t ) i=1,..,N d-dimensional Wiener process independent of the (W ij,t ) 1≤i j≤N involved in the DPDE fluctuation/dissipation interactions. The above dynamics can be seen as a superposition of a Langevin dynamics (1.54) and a DPDE fluctuation/dissipation dynamics (2.27). The numerical integration is therefore performed by a Trotter splitting of (4.1) between these two parts, the Langevin dynamics being integrated by the GLA scheme given by (1.70) and the DPDE fluctuation/dissipation dynamics by the Eulermaruyama scheme. We set γ = 10 -12 kg.s -1 (which corresponds to a fluctuation parameter σ = 2 γ β = 1.66 × 10 -16 kg.m.s -3/2 or σ = 23.27 in reduced units). Once equilibration is performed, we remove the periodic boundary conditions in the x direction, and put a wall of fixed Lennard-Jones particles of infinite masses on the left side of the simulation box. We next set the system in motion towards the left wall by adding a velocity equal to v init = -2000 m.s -1 to all the particles. This leads to the apparition of a shock wave propagating from the left to the right of the simulation box, with an average null velocity in the shocked state, as illustrated in Figure 4.1. The

t = 0 m = ∞ x = 0 v av ≃ v init t > 0 m = ∞ x = 0 v av ≃ 0 v shock v av ≃ v init Figure 4
.1 | Schematic of the shocked system simulation procedure. At t = 0, we set the equilibrated system in motion toward the left wall at a velocity v init = -2000m s -1 (left drawing). The system is compressed when it reaches the wall, and a shock wave propagates at velocity v shock > v init from the left to the right of the simulation box (right drawing). The particles on the left of the shock front are hotter and more densely distributed than those on the right. absence of a confining piston at the right end of the system causes the apparition of a relaxation wave, propagating from the right to the left end of the simulation box. This relaxation wave however does not affect the shock front until it reaches it, and averages over slabs at the shock front can safely be computed during the time where both waves do not encounter.

We present in the following sections shock profiles where thermodynamical quantities are averaged over slabs normal to the shock direction. We expect the profiles to present two zones delimited by the shock front:

• The zone on the left of the shock front is the shocked zone, where v av = 0 (red zone in Figure 4.1). In this zone, the slab averages of an arbitrary thermodynamical quantities ϕ are homogeneous and equal to a shock value ϕ shock .

• The zone on the right of the shock front is the equilibrium zone, where the particles have not yet seen the shock front and therefore are still in the equilibrium state corresponding to the equilibration procedure given above (pink zone on Figure 4.1). In this zone, the slab averages of an arbitrary thermodynamical quantities ϕ are homogeneous and equal to the equilibrium value ϕ equ .

During the simulation, the front shock delimiting the two zone will move from the left to the right of the profile. The behavior of the front shock depends on the observable considered. Observables with fast response time exhibit profiles with a peak value at the shock front and an equilibration towards ϕ shock after the passage of the shock front.

Observables with slow response time do not exhibit profiles with high values at the shock front, and directly start to equilibrate towards the shock value after the passage of the shock front. The kinetic temperature estimator (2.24) for instance has a much faster response time than the internal temperature estimator (2.26) if the internal heat capacity is large enough. Indeed, C v represents the number of internal degrees of freedom, and therefore the higher it is the more stable is the internal temperature and the longer it takes to equilibrate to the shock value. Such phenomenons are illustrated in In this case, the shock propagates along the x axis. The purple line represents the behavior of an observable ϕ that has a large inertia and slow response time: it equilibrates after the passage of the shock front towards the shock value. The red line represents the behavior of an observable ϕ that has little inertia and fast response time: it reaches a peak at the shock front and then equilibrates towards the shock value. Both observable produce a relaxation wave that propagates with a velocity v relax from the right end to the left of the simulation box. This relaxation wave however does not affect the shock front until both waves collide.

Numerical Results

The comparative study of the numerical results of the shock simulations are performed by considering two dynamical properties: the shock profiles and the shock front velocities.

Study of the shock profile

Instantaneous internal temperature profiles are computed, i.e profiles corresponding to the state of the system at a given time t of the simulation. Averages are computed over slabs of size ∆ slab ≃ 1.36 × 10 -9 m normal to the shock propagation (in this case normal to th x-axis). Figure 4.3 displays simulation results obtained with the SEM scheme with a small time step ∆t = 10 -15 s. Note that no stationary profile is obtained as both the kinetic and internal temperature increases in time. This increase is an artifact of the numerical scheme: we indeed confirmed by additional simulations reported in Figure 4.4 that the increase in the temperature is more pronounced for larger time steps. On the other hand, for shock simulations performed with the Hybrid or SER schemes, no such drift is observed. Figure 4.5 presents Hybrid and SER internal and kinetic temperature profiles for various times, for a simulation performed with the same time step ∆t = 10 -15 s. We notice that Hybrid and SER yield similar results. We also notice that, as expected, the kinetic temperature profiles exhibit peaks at the shock front while the internal temperature profiles do not. The results are very similar for Hybrid and SER simulations with time steps up to ∆t = 4 × 10 -15 s, as reported in Figure 4.6.

Study of the shock velocity

After having studied the computed shock profiles by all the schemes, we consider the shock velocities of SEM, SER and Hybrid simulations. In order to define the shock velocity, we must first define the shock position.

The particles at equilibrium have an average velocity of v init = -2000m s -1 , as explained in Section 4.1. The particles in the shocked state have a null average velocity because the particles of the left piston are of infinite mass. In the shock front, particles decelerate very fast, with the average velocity going from v init to 0, as displayed in Figure 4.7. The shock front position x shock is defined as the inflection point of the instantaneous velocity along the x-axis profiles at time t. However, we see in Figure 4.7 that the transition from equilibrium to shock is performed in no more than 2 slabs, which might cause instabilities when approximating the second derivative of the velocity profiles. Therefore, we estimate the shock position as n∆ slab , where n satisfies

v (n+1) x -v (n-1) x ≥ 2 |v init | 3 , ( 4.2)
with v

(n)

x being the average on the slab indexed by n of the velocity according to the x-axis, and ∆ slab being the size of the slabs (∆ slab ≃ 1.36 × 10 -9 m as said before).

Figure 4.8 shows that the above procedure reveals that the shock front propagates linearly in time. We thus can apply a linear fit to the estimated positions and define the velocity of the shock front as the drift of the fit, i.e

x shock (t) = x 0 + v shock t.

We display in Table 4.1 the computed shockwave velocities of SEM, SER and Hybrid schemes for time steps ranging from ∆t = 10 -15 s to ∆t = 4 × 10 -15 s. We notice that SEM, in addition to not obtaining a stationary profile, accelerates the shock propagation com- pared to SER or Hybrid, whose shock velocities are similar. Additional simulations of the exact same systems but with a different random seed for the random number generator showed that the computed shock velocities did not vary of more than 0.5% from the values presented in Table 4.1. This shows that the schemes we developed for equilibrium simulations can be used to simulate dynamical evolutions out of equilibrium, with time steps which are even larger than the ones used for standard deterministic dynamics of fully atomistic systems.

SEM

Hybrid SER ∆t = 1fs 2.268km s -1 2.201km s -1 2.202km s -1 ∆t = 2fs 2.323km s -1 2.206km s -1 2.208km s -1 ∆t = 4fs 2.412km s -1 2.213km s -1 2.206km s -1 Table 4.1 | Velocities of the shock front for SEM, SER and Hybrid simulations with time steps ranging from ∆t = 10 -15 s to ∆t = 4 × 10 -15 s. Additional simulations of the exact same systems but with a different random seed for the random number generator showed that the computed shock velocities did not vary of more than 0.5% from the above values.

an appropriate rescaling of σ, depending on the fraction of neglected pairwise interactions.

We organize this chapter by starting in Section 5.1 with the presentation of the DPDEscheme, and prove that it is of weak order one provided the fluctuation σ is rescaled appropriately. We then study in Section 5.2 the influence of the discard ratio on DPDEsimulations results, and compare DPDE-to the other schemes introduced in Chapter 2.

The SSA parallelization problem, as explained in Section 2.2.4, comes from couples where each of the particles is allocated to a different processor. We could therefore extend the idea of the DPDE-from a random rejection criterion to a spatial criterion: fluctuation/dissipation between "internal" particles are integrated using a SSA discretization while fluctuation/dissipation interactions involving one "frontier" particle are neglected (see Figure 5.2 for an illustration). The friction is then renormalized accordingly to this new criterion. However, we show in Section 5.3 that this spatial rejection criterion produce numerical artifacts in the frontier between domains that do not vanish in the limit of infinitely small time steps, and therefore cannot be used.

Discarding some interactions in SSA integration: DPDE-As mentioned in the introduction of this chapter, the DPDE-scheme consists in applying a SSA discretization where some elementary fluctuation/dissipation interactions are randomly discarded. Consider the DPDE elementary fluctuation/dissipation dynamics given by (2.39) for particles i and j. Let us denote by η ∈ [0, 1) the discard probability, and r η the corresponding rejection function, which reads where x ∈ X , G ∼ N (0, Id), u ∼ U (0, 1) and Φ SSA,ij ∆t,η (x, G) the result of an SSA discetization of the same elementary dynamics with renormalized parameters γ η and σ η used instead of γ and σ. Note that only the components p i , p j , ε i , ε j of the configuration of the system are changed by these one-step iteration maps. The definition (5.1) means that, the variables (p i , p j , ε i , ε j ) are updated according to a SSA scheme with probability 1 -η, while the configuration is unchanged with probability η. where L ij,η is the generator of the elementary DPDE fluctuation/dissipation interaction with parameters γ η and σ η . We have seen in the proof of Proposition 2.1.2 that

r η (u) =        1 if u < 1 -η,
L ij,η = -γ ij,η v ij • A ij + σ 2 η 2 A 2 ij , (5.4) with A ij = ∇ p i -∇ p j - v ij,t 2 ∂ ε i + ∂ ε j .
If σ η satisfies (5.3), then we have

P SSA,ij

Numerical results

We now study the behavior of the DPDE-scheme for equilibrium sequential simulations. We start this validation by considering the influence of the rejection probability η on estimations of equilibrium observables, and thus define some empirical upper bound for it. We then compare the DPDE-estimations to the SEM, SER and SSA estimations of Chapter 3.

Influence of the rejection probability on DPDE-

In theory, the DPDE-scheme is consistent for any rejection probability η < 1. However, when η gets close to 1, numerical artifacts can appear, and it is therefore mandatory to study the DPDE-behavior for different values of η.

We consider here the same systems than in the sequential simulations of Section 3.2.2. Figure 5.3 presents estimations of the temperature using the estimators (2.24), (2.25) and (2.26) introduced in Section 2.1.2, for simulations of time t f = 10 5 with time steps ranging from ∆t = 0.005 to ∆t = 0.025 in reduced units, for DPDE-schemes with the rejection probability ranging from η = 0.2 to η = 0.8 and for SSA (i.e DPDE-with η = 0). As seen in Figure 5.3, the DPDE-estimations are very close to those of SSA and exhibit very small biases, thus requiring very long simulation times in order to obtain the convergence. We notice that all estimations increase with η. Note however that all the estimations of the DPDE-biases are similar to those of SSA, except maybe for the estimations using the kinetic temperature estimator with η = 0.8. We additionally notice that the configurational temperature estimations of all the schemes overlap. Potential energy and virial pressure estimations exhibit a similar behavior with however a weaker overlapping. This seems to indicate that the DPDE-idea very weakly affects the estimations of position-related observables. Therefore, the above results show that DPDE-still yields equilibrium observables estimations similar to those of SSA, for rejection probability up to η = 0.8.

DPDE-estimations of equilibrium properties

We compare now the DPDE-estimations of the temperature with those obtained with the SEM, SER and SSA schemes in Section 3.2.2. Similarly to Section 3.2.2, we consider the schemes to be corrected by the projection procedure of Section 3.2.1, and we use the same parameters. The DPDE-rejection probability is set to η = 0.5. Figure 5.4 shows the average temperatures obtained from the estimators (2.24), (2.25) and (2.26) for σ = 4, while results for the smaller fluctuation strength σ = 2 are reported in Figure 5.5. We notice that the difference between DPDE-and SSA estimations are negligible compared to the differences between the estimations of SSA and those of the other schemes. This shows that DPDE-is a good alternative to SSA when considering equilibrium sequential simulations.

The limit of DPDE-in parallel simulations

Randomly discarding some fluctuation/dissipation interactions does not change the invariant measure, and, when rescaling appropriately the fluctuation parameter σ, yields results similar to the ones obtained for SSA equilibrium simulations. However, the DPDEidea yields schemes that suffer from the same parallelization problem as SSA. This can be countered by turning the random discard criterion into a spatial criterion: elementary fluctuation/dissipation interactions are discarded when the particles under consideration are allocated to different processors (see Figure 5.2). Such a rejection procedure would simply remove any problem concerning the parallelization of SSA, because we simply do not consider the particle couples where the parallelization problems arise. However, such an approach can introduce spatial heterogeneities in the discretization of the dynamics: the "frontier" particles of Figure 5.2 would interact with much less particles than the internal particles in terms of fluctuation/dissipation interactions. One therefore has to take into account this heterogeneity when devising the parallel extensions of DPDE-. Let us give some details about this "parallel DPDE-" scheme. Consider again the DPDE elementary fluctuation/dissipation dynamics given by (2.39) for particles i and j. Consider a rejection function R ij (q) no longer depending on a uniform random variable u but on the positions q, defined as

R ij (q) =       
1, if (q i , q j ) are allocated to the same processor, 0, otherwise.

Consider γ ′ and σ ′ the renormalized friction and fluctuation parameters, where the renormalization procedure is made precise later on. We define the parallel DPDE-discretization similarly to the random case by

Φ DPDE-par,ij ∆t (x, G) = R ij (q)Φ SSA,ij,′ ∆t (x, G) + (1 -R ij (q))x, (5.6)
where x ∈ X , G ∼ N (0, Id) and Φ SSA,ij,′ ∆t is the SSA discretization of the DPDE elementary fluctuation/dissipation dynamics with parameters γ ′ and σ ′ instead of γ and σ. Note that, again, only the components p i , p j , ε i , ε j of the configuration of the system are changed by these one-step iteration maps.

In constrast to the random case, where the renormalization procedure naturally comes from the definition of the DPDE-scheme, we propose here two renormalization approaches:

• Global renormalization. We compute the ratio η n between the number of discarded elementary interactions and the total number of elementary interactions, and replace σ and γ with

σ n = σ 1 -η n , γ n = γ 1 -η n .
Note that the fluctuation/dissipation relation (2.16) still holds for σ n and γ n . In opposition to the random case, the ratio η n is not constant and depends on the positions of all particles. The idea is however that, in translation invariant systems with periodic boundary conditions, the probability for two particles to be on different processors should be constant in average. In fact, η n should not vary much in time for large homogeneous systems at equilibrium.

• Local renormalization. One may argue that a global rescaling may be too drastic since the random interactions will be too strong in some parts of the system, while still too weak in others. Indeed, the global renormalization does not take into account the heteregoneity of the parallel DPDE-integration. Performing a local renormalization where the fluctuation is enhanced only for particles experiencing a depleted fluctuation/dissipation environment may lead to weaker perturbations on the dynamics. In order to formalize this idea, we introduce the ratio η n i between the number of discarded elementary interactions of particle i with its neighbors and the number of neighbors of i. There are many ways to generalize (5.3) in order to renormalize the interaction between particles i and j. We chose for instance to replace σ with σ 1 2

      1 1 -η n i + 1 1 -η n j       .
Note that, from a practical viewpoint, the implementation of DPDE-requires a first quick pass on the system to compute the ratios of discarded interactions, before actually performing the updates of the momenta and energies in a second pass. Of course, interactions between particles whose distances are larger than the cut-off radius are never taken into account in the computation of the ratios used for the renormalization.

Unfortunately, both parallel version of the DPDE-scheme yield artificial discrepancies. We consider a system of N = 16, 000 DPDE particles alocated to 8 processors, with a Lennard-Jones conservative potential fitted for Argon, with ρ ≃ 0.579, C v = 50 and σ ≃ 16.16. Figure 5.6 displays the kinetic and internal temperature as well as the density profiles on the slabs of the cross-section (y,z) averaged over the whole trajectory, for DPDE-simulations using time steps from ∆t = 0.5 × 10 -3 to ∆t = 1.8 × 10 -3 in reduced units. We notice that artificial density peaks occur in the frontier between processors, and we see that these peaks do not vanish when the time steps goes to 0: they kept the same amplitude for simulations where the time step was pushed to ∆t ≃ 1.4 × 10 -4 (not reported here).

Additional simulations not reported here also showed that the amplitude of these peaks increase with the density of the system as well as the value of σ. The profiles of the kinetic and internal estimators seem however unaffected, and the disparity observed between the average values of the different kinetic temperature profiles of Figure 5.6 are due to the biases related to the time step: as reported in Section 3.2.2, these biases depend on the time step and on σ, and are much larger for a given ∆t and σ for the kinetic temperature than for the internal temperature. The same behavior is observed for the parallel DPDE-scheme with local renormalization, as displayed in Figure 5.7.

The fact that both global and local parallel DPDE-schemes yield numerical discrepancies that do not arise from the use of finite time steps but depend on the system configuration (our simulations have exhibited their dependence on the system's density and σ, but other factors might come into play) indicate that both schemes are not consistent and therefore should not be used. negative internal energy, one might think of several procedures in order to continue the simulation anyway.

• Kinetic energy can be transformed into internal energy for the i-th particle. If the kinetic energy is not sufficient, the kinetic energy of its neighbor can also be taken to compensate. However, this procedure does not preserve the total momentum.

• In order to preserve the total momentum, one can transfer some internal energies of the neighbors to the i-th particle. However, this method may introduce biases because it cools down a group of particle in order to warm up a single one. It thus might increase the probability of another particle having a negative internal energy, therefore making the procedure unstable. This procedure however ensures the momentum conservation.

• The third option is to use adaptive time steps: whenever a particle i has its internal energy turn negative, we reject the move and use a smaller time step ∆t ′ < ∆t to discretize the dynamics of the considered particle and its neighbors over [t, t + ∆t]. However, in order for the discretization to remain consistent, one should be cautious to use Gaussian numbers G ′ such that

⌊ ∆t ∆t ′ ⌋ k=1 G ′,k ij = G n ij ,
for every neighbor j affected by the time step change.

Note that these procedures have not been implemented in the simulations displayed in this work, and that their mathematical analysis is still missing.

  Physical dimension of a particle (typically d = 3) N Number of particles of the considered system N d Physical dimension of the whole system (typically N d = 2dN ) X Set of all possible configurations (position and momentum space in classical MD) N it Number of iterations of the considered numerical scheme ∆t Timestep of the numerical scheme π Invariant measure of the considered dynamics, unless otherwise mentioned π ∆t

  µ β,C v (dq, dp, dε) = δ P 0 (p)Z -1 β,C v e -βE(q,p,ε)+ S (ε) k B dqdpdε,où Z β,C v est une constante de normalisation. Dans la suite, nous supposons qu'il existe une equivalence entre ces deux ensembles. Ceci signifie que, pourvu queE µ β,C v [E] = E 0 , nous avons lim N →∞ E µ (N ) E 0 ,P 0 [ϕ(x 1 , .., x k )] -E µ (N ) β,C v [ϕ(x 1 , .., x k )] = 0,pour n'importe quelle fonction ϕ dépendant d'un nombre fini k de variables (x 1 , .., x k ), où µ

Figure 2 |

 2 Figure 2 | Étude des propriétés de préservation d'énergie des schémas SVV, SER et SSA. Gauche: énergie moyennée sur toutes les trajectoires en fonction du temps. Droite: coefficients de dérive par rapport au temps de l'énergie moyennée sur toutes les trajectoires en fonction du pas de temps.

  xxi

Figure 3 |

 3 Figure3| Estimations numériques de biais sur les moyennes des composantes de l'énergie en fonction du pas de temps. Haut : composante cinétique de l'énergie. Milieu : composante potentielle de l'énergie. Bas : composante interne de l'énergie.

Figure 5 |

 5 Figure 5 | Profils des températures cinétiques (gauche) et internes (droite) moyennées sur des tranches normales à la direction de propagation du choc pour un pas de temps ∆t = 10 -15 s. Haut : profils obtenus avec des simulations SEM. Milieu : profils obtenus avec des simulations Hybrid. Bas : profils obtenus avec des simulations SER.

Figure 6 |

 6 Figure 6 | Approximate scales of length and time of different descriptions of matter. The orange dashed box represents the domain where coarse-graining models are desirable.

2 . 1 -Figure 1 . 2 |

 2112 Figure 1.2 | Example of a Lennard-Jones potential where the energy and the distance is expressed in terms of the reference distance σ and the reference energy ε.

Figure 1 Figure 1 . 3 |

 113 Figure 1.3 | Illustration of fixed wall boundary conditions.

Figure 1 . 4 |

 14 Figure 1.4 | Example of a two-dimensional simulation box with periodic boundary conditions applied. Particle B interacts with particle C in the box 5 (main box) but also with the periodic image of particle A in box 8.

Figure 1

 1 Figure1.5 | Energy preservation by the Velocity-Verlet scheme (symplectic, see (1.12)) for time steps ∆t * = 0.001 and ∆t * = 0.005 on an atomistic system described by a Lennard-Jones potential. After a transition time, the energy stabilizes around an average value close to its initial value. The difference between the two increases with the time step, as well as the amplitude of the oscillations around the average stabilized value.

Condition 1 . 2 . 2 .

 122 Brownian motion [RB06] and f : R N d → R N d and b : R N d → R d×m satisfy the following assumption: Both f and b are functions of D, and there exists C > 0 such that:

  .47) where [•, •] defines the commutator operator given by (1.26). The neglected terms of (1.47) contain linear combinations of the commutators of A and B of order higher than 3. The power term of the right-hand side of (1.47) can be written as an infinite sum, where each term is a composition of the commutators of A and B. When considering the semigroups generated by A and B, we formally have e ∆tA e ∆tB = e ∆t(A+B)+ ∆t 2 2 [A,B]+... , for any ∆t > 0. Fixing ϕ ∈ D, this allows us to write e ∆tA e tB ϕ = e ∆t(A+B) ϕ + ∆t 2 r ∆t,ϕ , (1.48) where r ∆t,ϕ L ∞ W k

  e A/2 e B e A/2 = e A+B-1 24 [A,[A,B]]-1 12 [B,[A,B]]+... , (1.53)

2 0

 2 by assumption, this means that f ω+1 ∈ D 0 . Consider the projection operator Π on D 0 and Q ∆t the pseudo-inverse operator of P ∆t -Id ∆t introduced in the proof of Theorem 1.3.5 (see (1.42)). By injecting (1.44) into the above equation, we have X Π P ∆t -Id ∆t Πϕdπ ∆t = ∆t ω X (-f ω+1 )ϕdπ -∆t ω+1 X r ∆t,ϕ dπ.

Theorem 1 .

 1 4.1. [Ergodicity of the Langevin dynamics] Let us suppose that the parameters γ and σ satisfy (1.55). Then the Langevin dynamics is ergodic for the canonical probability measure µ β :

Figure 2 . 2

 22 builds upon the example of Figure 2.1 and illustrates the additional modifications of DPDE with the internal/external equilibration.

Figure 2 . 2 |

 22 Figure2.2 | Illustration of DPDE coarse-graining process: in addition to the DPD coarse-graining process (dashed lines and circles, see Figure2.1), the meso-scale particles have internal temperatures and energies (the redder the color the hotter the particle). Internal/external equilibration is done through energy exchanges, turning internal energy into kinetic energy and reversely.

Figure 2 . 3

 23 takes the example of Figure 2.2 and illustrates the additional conduction energy exchanges between DPDE particles.

Figure 2 . 3 |

 23 Figure 2.3 | Incorporating the heat conduction to the DPDE: in addition to the mechanical energy variations, thermal exchanges occur between particles with different internal temperatures. The intensity of the heat flows are proportional to the differences between internal temperatures and are represented by the orange arrows: the larger the arrow the more heat is transferred.

Figure 2 . 4 |

 24 Figure 2.4 | Illustration of the SSA parallelization issue: the enclosed box in the full lines rep-resents the space domain allocated to a given processor, the exterior region being allocated to other processors. The dashed boxes of length r cut represent the cells of the linear scaling techniques: a particle on a given cell can only interact with particles in the neighboring cells. The dashed particles ("ghost" particles) are particles allocated to other processors. The particles at a distance r > r cut of the frontier ("internal" particles, in red ) never interact with ghost particles: the sequentiality of the integration does not affect the computation speed because no additional idle time nor information transfer between processors are needed. The particles at a distance r < r cut of the frontier ("frontier" particles, in blue) interact with ghost particles: the processor eventually has to wait before accessing any ghost particle's data, or wait for another processors to access its frontier particle's data, thus generating delays.

Figure 2 . 5 |

 25 Figure 2.5 | Illustration of the SSA parallelization procedure of [LBM + 14] in a two-dimensional case. Each box represents a space domain allocated to a single processor. The zone highlighted in red defines the region of space where the central processor integrates the elementary fluctuation/dissipation interactions. The same procedure is used by the other processors in the neighboring domains (delimited by dashed lines). Between each stage, all the updated particle's data are sent to their respective processor.

Figure 2 . 6 |

 26 Figure 2.6 | Illustration of the idea behind the Hybrid scheme for a system decomposed in two domains associated with two processors p 1 and p 2 . Red interactions are integrated sequentially according to a SSA discretization, and the remaining blue interactions are dealt with by an SER sweep.

Figure 3 . 1 |

 31 Figure 3.1 | Computed radial distribution function g(r) of the considered fluid at equilibrium. The radial distribution function represents the probability for each atom to find a neighboring atom according to a distance r. The above plot represents a liquid, where particles are semi-organized, yielding peaks of decreasing amplitude that do not touch zero.

Figure 3 .

 3 Figure 3.2 gives a schematic of the above initial conditions sampling procedure.

  ∆t = 0.001 Langevin, ∆t = 0.001 & Overdamped, ∆t = 0.001 Langevin, ∆t = 0.001 & Overdamped, ∆t = 0.001

Figure 3 . 2 |

 32 Figure 3.2 | Procedure to sample initial configurations according to the DPDE canonical measure (2.15). The part concerning the generation of the initial conditions are highlighted in red.

Figure 3 .Figure 3 . 3 |Figure 3 . 4 |

 33334 Figure 3.3 | Time-dependent energy drift averaged over trajectories, as a function of time when ∆t = 0.006, for three-dimensional simulations.

Figure 3 . 5 |

 35 Figure 3.5 | Average energy drift rate α ∆t as a function of the time step. A fit α ∆t = K∆t ω is superimposed to the data.

Figure 3 . 6 |

 36 Figure 3.6 | Drift rates per unit time of the kinetic (top), potential (middle) and internal (bottom) energies.

Figure 3 . 7 |

 37 Figure 3.7 | Time-dependent drift for the kinetic (top), potential (middle) and internal (bottom) components of the energy, averaged over trajectories, as a function of time when ∆t = 0.006.

Figure 3 . 8 |

 38 Figure 3.8 | Computed radial distribution function g(r) of the considered fluid at equilibrium. We see that the higher density yields oscillations of greater amplitude than in Figure 3.1.

Figure 3 .Figure 3 . 9 |

 339 Figure 3.9 | Numerical estimations of the components of the energy as a function of the time step for σ = 4. Top: Kinetic energy (3.4) and potential energy (3.5). Bottom: Internal energy (3.6).

Figure 3 .

 3 Figure 3.10 | Numerical estimations of the components of the energy as a function of the time step for σ = 2. Top: Kinetic energy (3.4) and potential energy (3.5). Bottom: Internal energy (3.6).

Figure 3 .

 3 Figure 3.11 | Numerical estimations of the components of the energy as a function of σ 2 ∆t. Top: Kinetic energy (3.4) and potential energy (3.5). Bottom: Internal energy (3.6).

Figure 3 . 12 |

 312 Figure 3.12 | Numerical estimations of the equilibrium properties as a function of the time step for σ = 2 obtained with parallel simulations. Top: Kinetic energy estimates (3.4). Middle: Average potential energy (3.5). Bottom: Internal energy estimates (3.6). SSA estimates have been computed on a system of N = 4000 particles on a single processor.

  schemet cpu ratio to SSA τ

Figure 3 .

 3 Figure 3.13 | Kinetic energy estimations for Hybrid simulations with a number N x of particles along the x axis equal to 40, 100 and 160. The SSA and SER estimations displayed are those reported in Figure 3.12.

Figure 4 Figure 4 . 2 |

 442 Figure 4.2 | Example of a profile averaged over slabs normal to the shock propagation for an arbitrary observable ϕ.In this case, the shock propagates along the x axis. The purple line represents the behavior of an observable ϕ that has a large inertia and slow response time: it equilibrates after the passage of the shock front towards the shock value. The red line represents the behavior of an observable ϕ that has little inertia and fast response time: it reaches a peak at the shock front and then equilibrates towards the shock value. Both observable produce a relaxation wave that propagates with a velocity v relax from the right end to the left of the simulation box. This relaxation wave however does not affect the shock front until both waves collide.

Figure 4 . 3 |

 43 Figure 4.3 | Kinetic (left) and internal (right) temperature average on the slabs of the crosssection (y, z) for the SEM scheme integrated with a time step ∆t = 10 -15 s.

Figure 4 . 4 |

 44 Figure 4.4 | Instantaneous kinetic (left) and internal (right) temperature profiles on the slabs of the cross-section (y, z) for the SEM scheme integrated with time steps ranging from ∆t = 10 -15 s to ∆t = 4 × 10 -15 s at time t = 4 × 10 -11 s.

Figure 4 . 5 |

 45 Figure 4.5 | Instantaneous kinetic (left) and internal (right) temperature profiles on the slabs of the cross-section (y, z) for the Hybrid (top) and SER (bottom) schemes integrated with a time step ∆t = 10 -15 s.

Figure 4 . 6 |

 46 Figure 4.6 | Instantaneous kinetic (left) and internal (right) temperature profiles on the slabs of the cross-section (y, z) for Hybrid (top) and SER (bottom) integrated with time steps ranging from ∆t = 10 -15 s to ∆t = 4 × 10 -15 s at time t = 4 × 10 -11 s.

Figure 4 . 7 |

 47 Figure 4.7 | Instantaneous velocity on the x-axis profiles averaged on the slabs normal to the x-axis of an Hybrid simulation or ∆t = 10 -15 s (SER yields similar results).

Figure 4 . 8 |

 48 Figure 4.8 | Computed shock positions defined by (4.2) for SEM, SER and Hybrid simulations for ∆t = 10 -15 s. The SER and Hybrid positions are on top of each other while SEM exhibits a slightly faster shock propagation.

Figure 5 . 2 |

 52 Figure 5.2 | The enclosed box in the full lines represents the space domain allocated to processor p, the exterior region being allocated to other processors. Internal-internal or internalfrontier elementary fluctuation/dissipation interactions are discretized with a SSA interaction. Frontier-ghost elementary fluctuation/dissipation interactions are neglected.

  0 otherwise.Consider γ η and σ η the friction and fluctuation parameter renormalized according to η, the renormalization being made precise later on (the parameter γ η is defined by the alternative DPDE fluctuation/dissipation relation given by (2.21)). The result of a DPDEdiscretization of the elementary DPDE interaction (2.39) between particle i and j writesΦ DPDE-,ij ∆t,η (x, G, u) = r η (u)Φ SSA,ij,η ∆t (x, G) + (1 -r η (u))x,(5.1)

3 )

 3 Respectively denote by P DPDE-,ij ∆t,η and P SSA,ij,η ∆t the transition operators associated with Φ DPDE-,ij ∆t,η and Φ SSA,ij ∆t,η . A straightforward computation shows thatP DPDE-,ij ∆t,η ϕ(x) = ηϕ(x) + (1 -η)P SSA,ij ∆t,η ϕ(x).(5.2).Proposition 5.1.1. The DPDE-discretization given by (5.1) is weakly consistent of order 1In addition, its evolution operator satisfies ∀ϕ ∈ D, P DPDE-,ij ∆t ϕ = e ∆tL ij ϕ + ∆t 2 r ∆t,ϕ , where r ∆t,ϕ stays uniformly bounded for small ∆t in the sense of Definition 2.2.1 and L ij is the generator of the elementary DPDE fluctuation/dissipation dynamics given by (2.39). Using renormalization (5.3), we see that the DPDE fluctuation/dissipation relation (2.16) still holds, i.e γ ij,η = γ Fix ϕ ∈ D and x ∈ X . We have proved in Proposition 2.2.4 that P ∆t,η ϕ(x) = ϕ(x) + ∆t L ij,η ϕ(x) + O(∆t 2 ),

1 .

 1 ∆t,η ϕ(x) = ϕ(x) + 1 1 -η L ij ϕ(x) + O(∆t 2 ).Therefore, applying (5.2) then (5.4) gives usP DPDE-,ij ∆t,η ϕ(x) = η ϕ(x) + (1 -η)P SSA,ij ∆t,η ϕ(x), = ϕ(x) + (1 -η)∆t L ij,η ϕ(x) + O(∆t 2 ), = ϕ(x) + ∆t L ij ϕ(x) + O(∆t 2 ).which proves the result.The full DPDE-scheme uses the same splitting than SSA, and its transition operator therefore writes P DPDE- the following proposition immediatly derives from Proposition 5.1.Proposition 5.1.2. The DPDE-discretization of the DPDE dynamics is weakly consistent of order one, provided thatσ η = σ 1 -η , and γ η = γ 1 -η .In addition, its evolution operator satisfies ∀ϕ ∈ D, P DPDE- ∆t,η ϕ = e ∆tL ϕ + ∆t 2 r ∆t,ϕ , where r ∆t,ϕ stays uniformly bounded for small ∆t in the sense of Definition 2.2.1.Remark 5.1.1. [Parallelization of DPDE-] As it stands here, DPDE-obviously suffers from the same parallelization issues than SSA. Unless using the implementation procedure introduced in [LBM + 14], DPDE-cannot be used for parallelization. Extensions of DPDE-for parallel simulations using a spatial rejection criterion are discussed in Section 5.3.

Figure 5 . 3 |

 53 Figure 5.3 | Numerical estimations of the equilibrium temperature as a function of the time step for σ = 4, for DPDE-with various rejection probability and SSA. Top: Kinetic temperature (2.24) (left) and potential temperature (2.25) (right). Bottom: Internal temperature (2.26).

Figure 5 . 4 |

 54 Figure 5.4 | Numerical estimations of the equilibrium temperature as a function of the time step for σ = 4 (with η = 0.5). Top: Kinetic temperature (2.24) and potential temperature (2.25). Bottom: Internal temperature (2.26). Results of DPDE-and SSA are on top of each other.

Figure 5 . 5 |

 55 Figure 5.5 | Numerical estimations of the equilibrium temperature as a function of the time step for σ = 2 (with η = 0.5). Top: Kinetic temperature (2.24) and potential temperature (2.25). Bottom: Internal temperature (2.26). Results of DPDE-and SSA are on top of each other.

3 Figure 5 . 6 |

 356 Figure 5.6 | Averaged virial XX pressure profiles on the slabs of the cross-section (y, z) for the parallel DPDE-scheme with global renormalization integrated with time steps ∆t ≃ 0.58 × 10 -3 (violet), ∆t ≃ 1.16 × 10 -3 and ∆t = 1.74 × 10 -3 in reduced units. Top: average kinetic (left) and internal (right) temperature profiles. Bottom: average density (right) profile.

3 Figure 5 . 7 |

 357 Figure 5.7 | Averaged virial XX pressure profiles on the slabs of the cross-section (y, z) for the parallel DPDE-scheme with local renormalization integrated with time steps ∆t ≃ 0.58 × 10 -3 (violet), ∆t ≃ 1.16 × 10 -3 and ∆t = 1.74 × 10 -3 in reduced units. Top: average kinetic (left) and internal (right) temperature profiles. Bottom: average density (right) profile.

  Numerical methods for SDEs . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 1.3.1 Elements on the discretization of SDEs . . . . . . . . . . . . . . . . . 25 1.3.2 Errors on the invariant measure . . . . . . . . . . . . . . . . . . . . . 31 1.3.3 Splitting techniques for SDEs . . . . . . . . . . . . . . . . . . . . . . 39 1.4 A prototypical example: the Langevin Dynamics . . . . . . . . . . . . . . . 45 1.4.1 Description of the dynamics . . . . . . . . . . . . . . . . . . . . . . . 45 1.4.2 Numerical Integration of the Langevin dynamics . . . . . . . . . . . 49 Integrating DPDE: state of the art . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Stochastic Velocity-Verlet Algorithm (SVV) . . . . . . . . . . . . . . 2.2.2 Splitting Euler-Maruyama (SEM) . . . . . . . . . . . . . . . . . . . . 2.2.3 Shardlow Splitting Algorithm (SSA) . . . . . . . . . . . . . . . . . . 2.2.4 Shardlow Splitting Algorithm parallelization issue . . . . . . . . . . 2.3 New parallelizable schemes for the integration of DPDE . . . . . . . . . . . 2.3.1 Splitting with Energy Reinjection (SER) . . . . . . . . . . . . . . . . 2.3.2 A scheme mixing SSA and SER: Hybrid . . . . . . . . . . . . . . . . 3 Comparative study of numerical schemes 3.1 Energy conservation properties of the schemes . . . . . . . . . . . . . . . . 3.1.1 Drifts of the total energy . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Drifts of the components of the energy . . . . . . . . . . . . . . . . . 3.2 Error estimations of equilibrium observables . . . . . . . . . . . . . . . . .

	Table des matières	xxv
	Context	1
	1 Introduction to Molecular Simulation 1.1 Microscopic description of matter . 2 DPD and its numerical integration Conclusion 2.Perspectives Bibliographie	5 53
	xxvi

. . . . . . . . . . . . . . . . . . . . . . . 6 1.1.1 The classical microscopic description of matter . . . . . . . . . . . . 6 1.1.2 Dynamics of isolated systems . . . . . . . . . . . . . . . . . . . . . . 11 1.2 Going from a microscopic description to a macroscopic one . . . . . . . . . 16 1.2.1 Introduction to Statistical Physics . . . . . . . . . . . . . . . . . . . . 16 1.2.2 Computing observables with SDEs . . . . . . . . . . . . . . . . . . . 20 1.3 1 The Dissipative Particle Dynamics (DPD) . . . . . . . . . . . . . . . . . . . 54 2.1.1 Original DPD model . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 2.1.2 The Dissipative Particle Dynamics with Energy conservation (DPDE) 58 xxv 2.2 3.2.1 Projection in order to ensure the energy preservation . . . . . . . . . 3.2.2 Biases on the average equilibrium properties for projected schemes 3.2.3 Validation of the Hybrid scheme: influence of the parallelization . .

4 Application of the DPDE to shock waves 4.1 Description of the simulations . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 An attempt at accelerating SSA equilibrium simulations 5.1 Discarding some interactions in SSA integration: DPDE-. . . . . . . . . . . 114 5.2 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 The limit of DPDE-in parallel simulations . . . . . . . . . . . . . . . . . . .

  Definition 1.1.2. [Local error]The local error of a numerical scheme is the error of its discretization procedure Φ ∆t compared to the analytical solution, i.e ∀y ∈ R d , e ∆t (y) = ϕ ∆t (y) -Φ ∆t (y) ,

	where (ϕ t ) t∈[0,T max ] denotes the semigroup formed by the flow maps of (1.5).
	Definition 1.1.3. [Global error] Let T > 0 be an integration time. The global error of a
	numerical scheme is the maximum difference between the true solution and its approximation
	over a time interval [0, T ], i.e
	∀y ∈ R d , E ∆t,T (y) = max 0≤n≤N it	ϕ n∆t (y) -Φ n ∆t (y) .
	We now define the notion of local consistency needed later on to obtain the conver-
	gence of a scheme.
	Definition 1.1.4. [Local consistency] A numerical scheme is said to be locally consistent or
	consistent if, for any x ∈ R d for which a solution to (1.5) exists, we have
	lim ∆t→0	e ∆t (x) ∆t	= 0.

  The dots are realizations of the EM scheme and the line is a discrete approximation of the process xt using a very small time step.

				dx t = x t /(1+x t 2 )dt + σdW t
		8						x
		5
	x	2
		-1
		4
		0	5	10	15	20	25	30	35	40
						t

n approx Figure 1.6 | Example of the approximation xt defined by (1.31) of the Euler-Maruyama discretization of the dynamics dx t = x t 1+x 2 t dt + σdW t , with σ = 3.

  .51) with r∆t,ψ ∈ D uniformly bounded in L ∞ W n for some n ∈ N and any ∆t small enough. Therefore, R ∆t,ϕ also is and we have the result.

	Therefore,
	e ∆tA r B ∆t,ϕ = r B ∆t,ϕ + ∆t r∆t,r B ∆t,ϕ	.
	Both r B ∆t,ϕ and r∆t,r B ∆t,ϕ

are uniformly bounded in some weighted space for any ∆t small enough: r B ∆t,ϕ by assumption and r∆t,r B ∆t,ϕ by (1.51) and because r B ∆t,ϕ ∈ D.

  2 e ∆tB e ∆t A 2 ϕ +

				∆t 3 8	r A ∆t,e ∆t,ϕ
	+	∆t 6 8	r A ∆t,r B ∆t,e ∆tA ϕ	+ e ∆tA r B ∆t,r A ∆t,ϕ	+	1 8	r A ∆t,e ∆tB r A,ϕ +	∆t 9 64	r A ∆t,r B ∆t,r A ∆t,ϕ

∆tB e ∆tA ϕ + e ∆tA r B ∆t,e ∆tA ϕ + e ∆tA e ∆tB r A , = e ∆t A 2 e ∆tB e ∆t A 2 ϕ + ∆t 3 r A,B ∆t,ϕ , = e ∆tL ϕ + ∆t 3 r ∆t,ϕ + r A,B ∆t,ϕ ,

  Proposition 2.2.4. The SSA discretization of the DPDE dynamics (2.9) is of weak order one. In addition, its evolution operator satisfies ∀ϕ ∈ D, P SSA ∆t ϕ = e ∆tL ϕ + ∆t 2 r ∆t,ϕ , where r ∆t,ϕ is uniformly bounded for small ∆t in the sense of Definition 2.2.1.Proof. Similarly than for the proof of Proposition 2.2.3, we only need to show that P

	where	P SSA,fd ∆t	= P 1,2 ∆t . . . P N -1,N ∆t	.	(2.42)
					ij
					∆t
					.41)

  .45) Remark 2.3.1. Let us emphasize that, by construction, SER automatically ensures the exact conservation of the total energy during the numerical integration of the fluctuation/dissipation part (2.27) -a very nice feature only enjoyed by SSA among the previously known schemes to integrate DPDE. On the other hand, in opposition to SSA, SER does not rely on a further splitting of the fluctuation/dissipation into elementary pairwise interactions, and is therefore straightforward to parallelize.

Table 3 .

 3

1 | Column 2: average CPU time per particle per time step and per processor. Column 3: timings of column 2 divided by the SSA timing. Column 4: Required CPU time corresponding to

  .2.

		1.502 1.504				SSA SER Hybrid 40 Hybrid 100 Hybrid 160
		1.5
	kin (J)	1.498
	E	1.496
		1.494
		1.492
		1.49
		0	0.0006	0.0012	0.0018	0.0024	0.003
				∆t(s)

Table 3 .

 3

2 | N /N proc ratio according to the system size.

R∆t,ϕ (x, G),

Remerciements

with Z = X e -βU(q) dq is a normalization constant. The probability measure ν β is in fact the unique invariant probability measure of (1.59), since the Overdamped Langevin dynamics is ergodic and ν β is invariant.

Theorem 1.4.3. The Overdamped Langevin dynamics given by (1.59) is ergodic for the probability measure ν β given by (1.60), provided that the potential U is such that ∇U satisfies Condition 1.2.2.

The proof of the above theorem is done in Section 2.2 of [START_REF] Lelièvre | Free Energy Computation : A Mathematical Perspective[END_REF]. We however report it below.

Proof. For any (ϕ, ψ) ∈ (L 2 (ν β )) 2 , the generator L over of (1.59) and its adjoint write

We clearly have that (L over ) * ν β = 0 with ν β given by (1.60). In opposition to the Langevin dynamics, the generator of the Overdamped Langevin dynamics is elliptic, i.e its diffusion matrix σ = 2 β Id is positive definite. The Overdamped Langevin dynamics is therefore ergodic with respect to ν β .

Remark 1.4.3. Overdamped Langevin can be seen as the limit of the Langevin dynamics when the friction parameter γ goes to infinity, or equivalently when the particles mass m goes to infinity, hence the use of "Overdamped" in its name. The proof of the convergence is given in [START_REF] Lelièvre | Free Energy Computation : A Mathematical Perspective[END_REF], and we do not report it here.

Similarly to the Langevin dynamics, it can be proved that the Overdamped Langevin dynamics satisfy condition 2 of Theorem 1.3.5. The proof however is simpler than for the Langevin dynamics if the position space is bounded, the generator of the Overdamped Langevin dynamics being elliptic (see discussion in Section 1.3.2). The case where the position space is unbounded is treated in [START_REF] Kopec | Weak backward error analysis for overdamped langevin processes[END_REF]. Proposition 1.4.4. Both L and L † , denoting respectively the generator of the Overdamped Langevin dynamics and its adjoint on L 2 (π), are invertible D 0 .

Numerical Integration of the Langevin dynamics

We have seen in Section 1.4.1 that the Langevin dynamics is composed of a Hamiltonian part and a fluctuation/dissipation part. We show in the sequel that its generator can be decomposed into three subgenerators that can be analytically integrated. Thus, splitting schemes seem to be a convenient way to integrate the dynamics.

However, when the ergodicity of the splitting schemes cannot be proved as in the case for dynamics where the position space Q is unbounded, one might turn to implicit schemes to devise stable schemes [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally lipschitz vector fields and degenerate noise[END_REF][START_REF] Talay | Stochastic hamiltonian systems: exponential convergence to the invariant measure, and discretization by the implicit euler scheme[END_REF][START_REF] Kopec | Weak backward error analysis for langevin process[END_REF]. The case of implicit schemes is not considered in this thesis however.

2

The Dissipative Particle Dynamics and its numerical integration

The Dissipative Particle Dynamics (DPD) [HK92] is a particle-based coarse-grained model in which atoms, molecules or even groups of molecules are represented by a single mesoscale DPD particle. The time evolution of the mesoscale particles is governed by a stochastic differential equation, where dissipative and random forces allow to take into account some effect of the microscopic variability.

The original DPD model is however an isothermal model, and therefore cannot modelise systems where energy-conserving phenomena occur. An energy-conserving variant of the model was introduced independently in 1997 by Avalos and al. in [AM97] and Español in [START_REF] Español | Dissipative Particle Dynamics with Energy Conservation[END_REF], allowing for such simulations. This variant is called Dissipative Particle Dynamics with Energy conservation (DPDE). However, the efficient numerical integration of DPDE still requires some effort. While numerous efficient schemes were developed for DPD [GW97, PHF98, BVKP00, VKBP02] (see in particular [START_REF] Leimkuhler | On the numerical treatment of dissipative particle dynamics and related systems[END_REF] for a careful comparison), their adaption to DPDE leads to numerical schemes for which errors on average properties may be large even for time steps standardly used to integrate Hamiltonian dynamics [LBM + 14], and thus may require the use of extremely small time steps [START_REF] Eiyad | Natural Convection Heat Transfer Simulation using Energy Conservative Dissipative Particle Dynamics[END_REF][START_REF] Eiyad | Application of Dissipative Particle Dynamics to Natural Convection in Differentially Heated Enclosures[END_REF].

We start this chapter by presenting in Section 2.1 the DPD model and its energyconserving variant. We then present in Section 2.2 a review of the principal algorithms presented in the literature so far. Because of the high computational cost of implicit schemes, we limit ourselves in this thesis to explicit schemes. Finally, we present in Section 2.3 two new schemes for the integration of the DPDE.

Let us emphasize that the schemes we develop for DPDE may be of direct interest for other dynamics with similar structures, such as the Smoothed Dissipative Particle Dynamics (SDPD) [ER03]. This stochastic dynamics also preserves the energy of the system, and can be seen as the superposition of two elementary energy preserving dynamics, a conservative part and a fluctuation/dissipation part. It is therefore not a surprise that the numerical methods we develop here can be readily applied for the large scale simulation of SDPD on parallel architectures.

Comparative study of numerical schemes on simulations at equilibrium

In this chapter, we perform numerical experiments in order to validate the new schemes proposed in Chapter 2. A good numerical scheme integrating the DPDE should fulfill two conditions: (i) preserve the invariants of the dynamics and (ii) have biases on average properties as small as possible.

Considering the first task, while the total momentum conservation is easily ensured, the energy conservation on the other hand is much more demanding. As said in Chapter 2, the numerical schemes we consider are obtained by a composition of a Verlet scheme, which approximately preserves the energy, and a discretization of the fluctuation/dissipation dynamics, which may or may not preserve the energy. Even when the energy is exactly preserved by the integration of the fluctuation/dissipation, the interaction between these two integrators may lead to drifts in the total energy, as already observed in Ref. [START_REF] Lisal | Dissipative Particle Dynamics at Isothermal, Isobaric, Isoenergetic, and Isoenthalpic Conditions using Shardlow-like Splitting Algorithms[END_REF]. This is described in Section 3.1.

All simulations performed in this chapter consider a system of N identical particles, with pairwise conservative interactions described by a shifted, splined and truncated Lennard-Jones potential (1.2). The heat capacity is supposed constant, which implies that equation (2.12) simply reduces to ε i = C v T i and that we have s(ε) = C v log(ε), where s is the internal entropy function. In addition, we only consider equilibrium average properties: dynamical properties are considered later on in the next chapter.

We start this chapter by studying in Section 3.1 the conservation of the energy by all the schemes introduced in Chapter 2. We notice that no scheme manages to keep it constant and that energy drifts occur. Similar drifts also occur in the components of the energy, thus making the estimations of observables dependent on the simulation time chosen, and forbidding the system of entering in a well-defined equilibrium state. In order to correct this, we devise in Section 3.2 a procedure projecting the system onto a constant-energy surface after each scheme iteration and compute the biases on equilibrium properties for the projected schemes, both on sequential and parallel simulations.

5

An attempt at accelerating Shardlow's Splitting Scheme at equilibrium

We have seen in Chapter 3 that SSA is the best performing scheme for sequential DPDE equilibrium simulations, even when SER and Hybrid are taken into account. The elementary splitting used by SSA for the discretization of the fluctuation/dissipation dynamics of DPDE thus seems, to the extent of our knowledge, to be much more advantageous than any global integration scheme presented to the community so far.

In this chapter, we perform an attempt at accelerating SSA estimations of equilibrium observables. A simple computation following the lines of the derivation of the invariant measure for DPDE [START_REF] Avalos | Dissipative Particle Dynamics with Energy Conservation[END_REF][START_REF] Español | Dissipative Particle Dynamics with Energy Conservation[END_REF] shows that measures of the form (2.18) are still invariant when some fluctuation/dissipation interactions are missing. Therefore, we can devise a new scheme, that we denote DPDE-, that uses the same splitting strategy than SSA but randomly discards some elementary fluctuation/dissipation interactions in the discretization of (2.27). This procedure allows to reduce the computational cost of the integration of the fluctuation/dissipation dynamics, thus accelerating the simulations. Figure 5.1 illustrates this idea of randomly discarding some fluctuation/dissipation interactions. However, the dynamical behavior of the system is modified by this procedure:

(5) (6) intuitively, the fluctuation/dissipation part is too weak, and should be strengthened by

Perspectives

We have seen in the previous chapters that both new schemes SER and Hybrid behaved consistently on equilibrium and nonequilibrium systems. Hybrid was found to be the most accurate scheme, with an accuracy comparable to SSA's, and the best scheme so far considering the CPU cost/accuracy ratio, as displayed in Table 3.1, when considering straightforwardly parallelizable schemes. Both Hybrid and SER schemes have been found to correctly reproduce the shocked behavior of Lennard-Jones DPDE systems whereas SEM could not.

However, Hybrid is a blending of SER into SSA, and therefore is cursed with the same problem as SSA: it is not adapted to thread parallelization. Thread parallelization is a process where several processor units have access to the same shared memory, in opposition to MPI parallelization (abusively named parallelization in this thesis), where each processor unit has its own memory, and every information between processors must be sent through special routines. The future architectures of supercomputers are architectures mixing thread and MPI parallelization, where nodes consisting of several processor units sharing the same memory communicate through the MPI standard. It is therefore desirable that numerical algorithms are adapted to these new architectures. Newtonian MD has already been adapted to this new standard [START_REF] Cieren | Molecular Dynamics for Exascale Supercomputers[END_REF]. As for DPDE, SEM, SVV and SER schemes are adapted to thread parallelization (we say that such schemes are threadable), but Hybrid, while constructed to be straightforwardly MPI parallelizable, is not threadable because of the sequentiality of its SSA part.

We mentioned in Section 3.2.2 that the SER integration for σ = 4 and C v = 10 yielded negative energies even for small time steps due to the energy reinjection procedure. More generally, additional simulations not reported in this work indicate that SER was more likely to produce particle with negative internal energies than Hybrid or SSA. In addition, we saw that the biases estimations of SER simulations, corrected by the energy reinjection procedure, were an order of magnitude larger than those of SSA or Hybrid. We thus have with Hybrid an accurate DPDE scheme that is MPI-paralellizable but inadapted to thread parallelization. On the other hand, we have with SER a threadable scheme, much more accurate than the remaining parallelizable options previous to this work, but much less accurate than Hybrid or SSA, and for which particles with negative internal energies are 123 Development of new numerical integration schemes for multiscale coarse-graining methods more likely to appear.

A fix to this accuracy, stability and thread-parallelization problem could be the construction of a weak order two scheme inspired from the SER integration. The SER-like integration would be threadable, and the higher accuracy arising from the weak order two could improve both the accuracy and the stability of the simulations. In order to devise such a scheme, we could start by devising a weak order two discretization of the elementary DPDE dynamics. The pairwise structure of the fluctuation/dissipation could then allow us to extend this discretization to the full fluctuation/dissipation DPDE dynamics. The Velocity-Verlet discretization being of order two, a Strang splitting would then yield a DPDE discretization of weak order two.

The idea behind the construction of a weak order two discretization of (2.39) is to "evaluate" γ ij,t+∆t/2 before integrating the momenta from t to t + ∆t with the analytical integration (2.40) and the evaluated friction parameter. This procedure allows to take into account the dependence of γ ij with the internal energies in the discretization of (2.39), and thus achieve weak order two accuracy. Indeed, in opposition to DPD, the friction parameter γ ij is no longer constant in DPDE but depends on T i and T j , and therefore on ε i and ε j by (2.12). Not taking this into account forbids obtaining a weak order two accuracy. This is why the adaptation of SSA to the DPDE framework reported in Section 2.2.3 is of weak order one, in opposition to the original SSA made for DPD which is theoretically of weak order two (if a Strang splitting is used).

Consider a procedure that estimates γ ij,t+∆t/2 as

The analytical integration (2.40) performed using γn+1/2

(5.7) Therefore, one should look for terms C n and D n such that (5.7) is of weak order two, according to procedures detailed in Section 1.3.1.

Apart from devising a weak order two scheme, we can devise corrections that can be used for existing schemes whenever negative internal energies are encountered. Indeed, even stable scheme can have a nonzero probability of yielding DPDE particles with negative energies. Therefore, a procedure allowing the simulation to continue whenever such events occur is desirable.

In a pratical point of view, whenever some particle, indexed by i for instance, has a

Conclusion

Previous to this thesis, DPDE lacked accurate and straightforwardly parallelizable schemes for its integration. The best-performing scheme, named SSA, is accurate but its sequential nature prevents it from being simply parallelizable. The other two remaining options are schemes called SEM or SVV that have both very poor energy preservation properties and thus need prohibitively small time steps in order to correctly integrate DPDE systems. The main Cea/Dam motivation for the use of DPDE being to perform massively parallel simulations of shocked systems in decent times, accurate and parallelizable DPDE schemes were desirable. Two new numerical schemes satisfying both conditions have been introduced in this thesis. The first one is termed SER and is adapted to the future architecture of supercomputers by being threadable while the second, a merge of SER and SSA appropriately termed Hybrid, is as accurate as the reference SSA scheme. We have shown that both schemes bring a substantial improvement compared to SEM or SVV. Indeed, on the first hand, we have shown that both of them are consistent with other existing schemes and have significantly better energy conservation properties than SEM or SVV. On the other hand, we showed that both accurately describe a shock propagation in a DPDE fluid for parameters at which SEM could not obtain a stationary state. In the negative side, SER was found to be more unstable than the other schemes by yielding particles with negative energies with higher probability. As for Hybrid, its blending into SSA makes it incompatible with thread parallelization. Nevertheless, this thesis defines Hybrid to be the new best option today for parallel DPDE integration, with an accuracy similar to SSA. However, SER seems to be the most promising option, and work on a weak second order scheme adapted from its energy reinjection procedure is currently being performed.

This thesis therefore gives two consistent options for massively parallel DPDE simulations, one adapted to nowadays supercomputers and the second that will allow massively parallel accurate DPDE simulations on future supercomputers. This will help the understanding of complex shocked systems, which is one of the objectives of the Cea/Dam.