# Simulations multi-échelles en Dynamique Moléculaire

Frédéric Boussinot

► To cite this version:

Frédéric Boussinot. Simulations multi-échelles en Dynamique Moléculaire. 2023. hal-04267212

HAL Id: hal-04267212

https://minesparis-psl.hal.science/hal-04267212v1

Preprint submitted on 1 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
Simulations multi-échelles en
Dynamique Moléculaire

Frédéric Boussinot

frederic.boussinot@gmail.com

Novembre 2023
# Table des matières

<table><tbody><tr><td>1</td><td>Introduction</td><td>5</td></tr><tr><td>1.1</td><td>La dynamique moléculaire</td><td>5</td></tr><tr><td>1.2</td><td>Modélisation multi-échelle</td><td>8</td></tr><tr><td>1.3</td><td>Créations/destructions dynamiques</td><td>10</td></tr><tr><td>1.4</td><td>Plan du texte</td><td>10</td></tr><tr><td>2</td><td>Alkanes</td><td>13</td></tr><tr><td>2.1</td><td>Liaisons</td><td>14</td></tr><tr><td>2.2</td><td>Angles de valence</td><td>15</td></tr><tr><td>2.3</td><td>Angles de torsion</td><td>16</td></tr><tr><td>2.4</td><td>Forces inter-moléculaires</td><td>17</td></tr><tr><td>3</td><td>Application des forces</td><td>19</td></tr><tr><td>3.1</td><td>Liaisons</td><td>19</td></tr><tr><td>3.2</td><td>Angles de valence</td><td>20</td></tr><tr><td>3.3</td><td>Angles de torsion</td><td>22</td></tr><tr><td>3.4</td><td>Forces inter-moléculaires</td><td>25</td></tr><tr><td>4</td><td>Implémentation</td><td>29</td></tr><tr><td>4.1</td><td>Programmation Réactive</td><td>30</td></tr><tr><td>4.2</td><td>SugarCubes</td><td>31</td></tr><tr><td>4.3</td><td>Méthode de résolution</td><td>34</td></tr><tr><td>5</td><td>Simulations en AA</td><td>37</td></tr><tr><td>5.1</td><td>Stabilité</td><td>38</td></tr><tr><td>5.2</td><td>Chaos déterministe</td><td>40</td></tr><tr><td>6</td><td>Approches multi-échelles</td><td>43</td></tr><tr><td>6.1</td><td>Echelle UA</td><td>43</td></tr><tr><td>6.2</td><td>Echelle CG</td><td>44</td></tr></tbody></table>
<table><tr><td>6.3 Complexité</td><td>47</td></tr><tr><td>7 Méthode Boltzmann-inverse</td><td>51</td></tr><tr><td>7.1 Forces intra-moléculaires UA</td><td>52</td></tr><tr><td>7.2 Forces inter-moléculaires UA</td><td>55</td></tr><tr><td>8 Boltzmann-inverse CG</td><td>57</td></tr><tr><td>8.1 Liaisons CG</td><td>57</td></tr><tr><td>8.2 Angles de valence CG</td><td>59</td></tr><tr><td>9 Méthode de minimisation</td><td>61</td></tr><tr><td>9.1 Liaisons UA</td><td>62</td></tr><tr><td>9.2 Angles de valence UA</td><td>62</td></tr><tr><td>9.3 Angles de torsion UA</td><td>64</td></tr><tr><td>9.4 Forces inter-moléculaire UA</td><td>65</td></tr><tr><td>10 Minimisation CG</td><td>69</td></tr><tr><td>10.1 Liaisons CG</td><td>69</td></tr><tr><td>10.2 Angles de valence CG</td><td>70</td></tr><tr><td>10.3 Pondération</td><td>71</td></tr><tr><td>10.4 Potentiel des liaisons CG</td><td>72</td></tr><tr><td>10.5 Potentiel des angles de valence CG</td><td>73</td></tr><tr><td>10.6 Forces inter-moléculaires CG</td><td>73</td></tr><tr><td>10.7 Comparaison avec l'approche Boltzmann-inverse</td><td>76</td></tr><tr><td>11 Simulations en UA</td><td>79</td></tr><tr><td>12 Simulations en CG</td><td>83</td></tr><tr><td>12.1 Détermination du potentiel CG</td><td>83</td></tr><tr><td>12.2 Exécution CG</td><td>86</td></tr><tr><td>13 Reconstructions</td><td>91</td></tr><tr><td>13.1 Forces inter-moléculaires</td><td>94</td></tr><tr><td>13.2 Reconstructions automatiques</td><td>97</td></tr><tr><td>14 Conclusion</td><td>101</td></tr></table>
# Chapitre 1

## Introduction

La simulation numérique en physique consiste à modéliser les systèmes sous forme de programmes informatiques et à exécuter ceux-ci pour étudier les propriétés des systèmes considérés. En dynamique moléculaire (DM) [3] les systèmes sont des molécules, formées d'atomes et l'objectif est de déterminer comment celles-ci évoluent au cours du temps. Les premiers travaux concernant la simulation à l'échelle atomique datent des années 50 [15], à une époque où les moyens informatiques étaient extrêmement limités. Les premiers systèmes simulés étaient formés d'atomes indépendants (non regroupés en molécules) soumis à des chocs parfaitement élastiques (gaz parfait). Ensuite, les forces de van Der Waals ont été introduites, pour obtenir des simulations plus proches de la réalité. Puis les interactions intra-moléculaires ont été prises en compte, pour obtenir de véritables modélisations moléculaires. La DM est actuellement principalement utilisée pour l'étude des matériaux ainsi qu'en biologie moléculaire.

### 1.1 La dynamique moléculaire

À la base de la DM on trouve la physique classique (newtownienne), avec l'équation fondamentale :

$$ \vec{F} = m\vec{a} \qquad (1.1) $$

où $\vec{F}$ est la force subie par une particule de masse $m$ et $\vec{a}$ son accélération (la dérivée seconde de la variation de sa position, par rapport au temps).

Les composants élémentaires utilisés pour modéliser les molécules sont les suivants :
* Les atomes, avec 6 degrés de liberté (coordonnées *x*, *y*, *z* et vitesses *sx*, *sy*, *sz*).

* Les liaisons (*bonds*) qui lient deux atomes *a*, *b* (dits *liés*) au sein d'une molécule ; la liaison tend à maintenir constante la distance *ab*.

* Les angles de valence (*valence angles*) définis par deux atomes *a* et *c* liés à un même troisième atome *b*, au sein d'une molécule ; l'angle de valence tend à maintenir à une valeur fixe l'angle *abc*.

* Les angles de torsion (*dihedral angles*), souvent appelés *dièdres*, qui relient quatre atomes *a*, *b*, *c*, *d* au sein d'une molécule ; *a* est lié à *b*, *b* à *c* et *c* à *d* ; l'angle de torsion tend à faire coincider les plans *abc* et *bcd*.

* Les interactions de van Der Waals, qui lient deux atomes, non nécessairement appartenant à une même molécule ; ces interactions varient suivant la nature des atomes.

Les modélisations moléculaires peuvent également prendre en compte les interactions électrostatiques (loi de Coulomb) ; nous ne considérerons pas cet aspect ici.

Les molécules sont formées d'atomes liés. Une distance intra-moléculaire est définie, pour délimiter les atomes appartenant à la molécule pour lesquels on considère que les forces de van Der Waals sont prises en compte par les liaisons, angles de valence et angles de torsion, et ne doivent donc pas être considérées en tant que telles.

La DM repose sur un algorithme de résolution numérique, dont le rôle est d'appliquer les forces engendrées par les divers composants en suivant les équations de la physique newtonienne. L'algorithme de résolution appliqué à une molécule a pour paramètre le *pas-de-temps* $\Delta t$ de simulation. Celui-ci est de l'ordre de la femto-seconde ($10^{-15}$s), pour des molécules en tout-atome (voir plus loin).

La méthode de résolution généralement utilisée est appelée (*velocity*) *Verlet*, du nom de son concepteur qui l'a proposée en 1967 [23]. Elle repose sur une résolution en deux étapes : lors de la première étape, on calcule pour un pas-de-temps de $\Delta t/2$ la vitesse et la position, et lors de la seconde étape, on calcule pour le pas-de-temps $\Delta t/2$ suivant, la nouvelle vitesse. Cette méthode possède deux atouts principaux : (1) elle est très stable (les simulations peuvent être longues, dans un sens que l'on précisera plus tard) ; (2) elle préserve l'énergie. La préservation de l'énergie (c'est-à-dire le fait que la résolution n'ajoute ni ne fait disparaître de l'énergie au cours de la simulation) est un critère fondamental pour les simulations en physique. Elle correspond à une symétrie profonde de la physique classique (indépendance par rapport
au sens de l'écoulement du temps).

Les forces intra-moléculaires (liaisons, valence, dièdres) ainsi que les forces inter-moléculaires (van Der Waals) sont conservatives : le travail entre deux points est indépendant du chemin suivi. Elles peuvent donc être définies comme des dérivées de champs scalaires, appelés *potentiels*. On a alors :

$$ \vec{F}(\mathbf{r}) = -\nabla U(\mathbf{r}) \qquad (1.2) $$

où **r** désigne les coordonnées du point sur lequel la force $\vec{F}(\mathbf{r})$ s'applique et $U$ est le potentiel dont la force dérive.

Les potentiels des divers composants sont regroupés sous le nom de champ de forces. Plusieurs champs de forces existent et sont employés dans divers contextes. Le champ de force OPLS [10] est souvent utilisé dans le cadre des simulations de liquides [3]; c'est ce champ de force que nous considérerons dans le texte.

La définition d'un champ de force peut nécessiter des calculs à l'échelle quantique; nous n'aborderons pas ici cet aspect.

Un système de DM comprend donc les éléments suivants :

* Des moyens de définir les atomes et molécules.

* L'implémentation des divers potentiels du champ de force.

* L'implémentation de la méthode de résolution.

* Des moyens de définir et d'effectuer les simulations.

De plus, on trouve généralement dans les systèmes de DM :

* Des moyens pour effectuer les simulations dans le contexte d'un thermostat (contrôle de la température) ou d'un barostat (contrôle de la pression).

* Des moyens pour simuler des quantités non-bornées d'atomes (conditions périodiques).

On ne considérera pas ces deux derniers aspects, dans la mesure où on se préoccupe ici uniquement des fondements de la DM.

Plusieurs dizaines de systèmes de DM sont disponibles, dont certains sont issus du monde académique et libres d'accès. Parmi ces derniers, on peut citer DL_POLY [1], GROMACS [19] et CHARMM [22]. Ces systèmes sont écrits en FORTRAN ou en C/C++ et ils sont interfacés avec des outils de visualisation 3D, à l'aide d'outils de traduction entre formats appropriés. Le travail présenté dans ce texte a été réalisé lors du développement du système de DM[6] qui a servi à mener toutes les simulations que l'on va décrire plus loin.
## 1.2 Modélisation multi-échelle

Les simulations en DM au niveau “tout-atome” (AA) constituent un outil très puissant pour analyser les systèmes moléculaires mais elles ont deux limitations liées au nombre d’atomes pouvant être simulés, d’un côté, et d’un autre côté au pas-de-temps du schéma de résolution utilisée.

En ce qui concerne la taille des systèmes pouvant être simulés, la restriction principale se situe dans la taille de la mémoire devant être utilisée. Quelques milliers d’atomes peuvent raisonnablement être traités par une machine standard. Ce nombre peut être significativement augmenté par l’utilisation de réseaux de machines distribuées s’exécutant en parallèle et partageant une mémoire distribuée à travers le réseau. Ainsi, on s’attend à ce que des grappes de machines distribuées puissent simuler des systèmes allant jusqu’au million d’atomes.

Le pas de temps des simulations est imposé par la physique : c’est une fraction de la période de la vibration la plus courte apparaissant dans le système. Par exemple, pour ce qui concerne les molécules hydro-carbonées décrites à l’échelle AA, le pas-de-temps dépend de la vibration de la liaison CH qui est $10^{-14}$ s ($\lambda = 2860 \text{cm}^{-1}$ dans le spectre infra-rouge). Habituellement, le pas-de-temps est au maximum le dixième de cette vibration, c’est-à-dire $10^{-15}$ s (1 femto-seconde). Avec un tel pas-de-temps, un million d’étapes sont nécessaires pour prédire la dynamique du système pendant une nano-seconde seulement. Ainsi, les modèles AA sont les plus réalistes au niveau chimique mais le temps pendant lequel ils peuvent être simulés dépasse rarement quelques nano-secondes. De très nombreuses simulations, traitant des phénomènes de diffusion par exemple, nécessitent des durées nettement plus longues. Cette limitation concernant le pas-de-temps de simulation apparaît comme étant la limitation principale, par rapport à celle concernant la taille du système simulé.

Il existe deux approches pour contourner la limitation concernant le pas-de-temps, en augmentant ainsi la durée de simulation. La première approche est appelée *Hyper-molecular dynamics* ou *Accelerated-molecular dynamics* [24, 16]. Elle consiste en une modification du champs de forces, obtenu en réduisant la barrière d’énergie potentielle entre deux états correspondant à deux événements rares du système moléculaire, de façon à augmenter la probabilité de transition entre ces états.

La seconde approche, largement répandue et qui est celle que l’on adopte dans ce texte, consiste à réduire la complexité du système moléculaire en re-
groupant certains atomes entre eux de façon à former des “grains”. Cette approche générale est appelée “*Coarse-Grained*”. La plus simple des réductions consiste à regrouper dans un seul grain un atome de carbone et les atomes d’hydrogène qui lui sont liés ; on dit que l’on a alors un modèle “*Atome Uni-fié*” (UA). On considérera également dans la suite, des grains formés de deux atomes liés de carbone, avec les atomes d’hydrogène qui leur sont liés ; c’est à cette réduction que l’on réservera le nom de CG.

La réduction de complexité des systèmes moléculaires autorise des durées de simulation plus importantes. Ceci résulte de plusieurs points :

* Le pas-de-temps de la méthode de résolution peut être augmenté par rapport à l’échelle *AA* puisque la plus courte période de vibration l’est également (la masse des grains CG est plus importante que celle des atomes *AA*).

* Le nombre de degrés de liberté et le nombre de composants intra-moléculaires (liaisons et angles) est réduit par rapport à *AA*.

* Le nombre d’interactions inter-moléculaires (van Der Waals) entre grains est réduit par rapport au nombre d’interactions inter-moléculaires entre atomes. Ce point est particulièrement important car les interactions inter-moléculaires peuvent être très nombreuses et donc coûteuses à simuler.

Tous ces points étant pris en compte, le temps de simulation peut être augmenté par plus de deux ordres de grandeur, en passant de l’échelle *AA* à l’échelle CG. Cette augmentation est tellement significative que la plupart des simulations en DM sont en fait effectuées aux échelles UA ou CG.

Le choix d’utiliser l’échelle UA ou CG nécessite cependant de transporter le champ de forces *AA* à cette échelle. Ce transport est assez facile de *AA* à UA [10, 11, 12, 22]. En fait, l’intégration des atomes d’hydrogène dans les grains UA nécessite uniquement un changement des paramètres des potentiels des dièdres et des interactions de van Der Waals. Les nouveaux paramètres peuvent être assez facilement déduits analytiquement de ceux de *AA*.

Pour le transport de *AA* à CG, les paramètres sont déduits d’une analyse statistique dans les matériaux denses pour reproduire certaines fonctions fondamentales de ces matériaux (densité, énergie, distribution) valides à l’échelle *AA*. Plus précisément, le champ de forces CG est construit à partir d’un traitement Boltzmann-inverse dans les références[18, 21, 20]. Un point important est que ce type de transfert suppose que la forme des potentiels en CG est la même qu’en *AA*. C’est un point que l’on contestera dans la suite du texte, en proposant une approche différente, justifiant une construction analytique de
CG et autorisant des *reconstructions* qui sont des changements dynamiques d'échelle entre AA, UA et CG.

## 1.3 Créations/destructions dynamiques

Ne pas pouvoir simuler d'autres systèmes que les systèmes figés, dans lesquels le nombre d'atomes et leurs liaisons ne varient pas, est une contrainte qu'on peut vouloir dépasser dans au moins deux cas :

* Pour simuler une réaction chimique, dans laquelle on peut avoir l'ap- parition au cours de la simulation d'une nouvelle liaison entre deux atomes (par exemple, un atome d'oxygène se liant à un atome de fer).

* Dans le contexte d'un changement d'échelle d'une molécule *en cours de simulation*. Le changement d'échelle peut être implémenté par un ajout dynamique de la nouvelle version de la molécule *simultanément* avec la destruction de l'ancienne version. Ce point est particulièrement important dans le contexte des reconstructions.

Il semble donc intéressant qu'un système de DM offre des mécanismes permettant de détruire ou de créer dynamiquement, c'est-à-dire en cours d'exécution, des composants de la simulation. C'est un système ayant une telle capacité fondée sur une technique de programmation particulière, la programmation réactive, que l'on a construit et que l'on a utilisé dans ce texte.

## 1.4 Plan du texte

Les alkanes, qui sont des molécules hydro-carbonnées linéaires que l'on utilisera tout au long du texte, sont présentées dans le Chap. 2. Les composants de ces molécules (liaisons, angles) sont décrits et les formes de leurs potentiels définissant leur énergie potentielle sont présentées sous forme de courbes (celles-ci correspondant au champ de forces OPLS[9] qui est celui de l'échelle "tout atome" (AA) dans ce texte).

Le Chap. 3 décrit à l'aide d'algèbre vectorielle comment les forces intra-moléculaires et inter-moléculaires sont appliquées sur les atomes à l'échelle AA.

L'implémentation est considérée dans le Chap. 4 en commençant par la programmation réactive, dans laquelle les programmes sont définis par rap-
port à une notion d'instant. La librairie Java SugarCubes[7] est brièvement décrite ainsi que la méthode de résolution temporelle velocity-Verlet.

Les simulations en AA sont considérées dans le Chap.5 sous l'angle de la stabilité et du déterminisme. Le cadre général de la DM est le chaos déterministe.

Dans le Chap.6 on introduit les échelles de description UA et CG dont on cherchera à déterminer les potentiels dans la suite.

La méthode Boltzmann-inverse pour déterminer un potentiel UA est présentée dans le Chap.7. Le cas de CG est considérée dans le Chap.8. La méthode Boltzmann-inverse s'avère insatisfaisante par plusieurs aspects qui justifient de considérer une autre approche "par minimisation".

Cette méthode alternative est définie pour l'échelle UA dans le Chap.9. La méthode par minimisation repose sur un lien géométrique établi entre les échelles UA et AA et sur des minimisations de l'énergie potentielle au niveau AA. Le cas de CG est considérée dans le Chap.10. La méthode par minimisation s'avère beaucoup plus satisfaisante que celle Boltzmann-inverse.

Les simulations à l'échelle UA sont considérées dans le Chap.11. Celles à l'échelle CG, dans le Chap.12.

Les reconstructions entre les trois échelles AA, UA et CG sont décrites dans le Chap.13.

Finalement, le Chap.14 conclut le texte.
# Chapitre 2

## Alkanes

Les *alkanes* sont des chaînes linéaires d'atomes de carbone auxquelles sont rattachés des atomes d'hydrogène. Les alkanes sont désignées par des formules de la forme $C_nH_{2n+2}$ où n est le nombre d'atomes de carbone. La Fig. 2.1 montre une molécule alkane $C_6H_{14}$, composée de 6 atomes de carbone et de 14 atomes d'hydrogène.

<figure><img src="image_2.png" /><figcaption>FIG. 2.1 – Chaîne hydro-carbonnée $C_6H_{14}$ (6 atomes de carbone, 14 atomes d'hydrogène).</figcaption></figure>

En DM, les molécules sont structurées en composants qui déterminent leur structure et les forces qui s'appliquent à leurs atomes. A chacun de ces composants est associé un *potentiel* qui est une fonction décrivant l'énergie du composant.
Les potentiels sont habituellement regroupés sous le terme de “champ de forces”. Les forces inter-moléculaires s’exerçant entre les molécules font également généralement partie des champs de forces. Dans la suite, à l’échelle AA on utilise exclusivement un champ de forces particulier, appelé OPLS[9], à la base du système de DM DL_POLY[1].

En OPLS, les principaux composants des molécules sont les liaisons reliant deux atomes, les angles de valence reliant trois atomes, et les angles de torsion reliant quatre atomes. Par exemple, en OPLS, la molécule $C_6H_{14}$ possède 19 liaisons, 36 angles de valence et 45 dièdres.

Aux deux extrémités des alkanes on trouve trois atomes d’hydrogène, tandis que deux seulement sont liés aux autres atomes de carbone. Pour simplifier, en ne distinguant plus les carbones aux extrémités des autres, on sera amené à considérer des *fragments* d’alkanes qui sont de la forme $C_nH_{2n}$ c’est-à-dire pour lesquelle tous les carbones sans exception ont deux hydrogène liés. La Fig.2.2 montre un fragment $C_6H_{12}$. Le fragment $C_6H_{12}$

<figure><img src="image_3.png" /><figcaption>FIG. 2.2 – Fragment $C_6H_{12}$ (6 atomes de carbone, 12 atomes d’hydrogène).</figcaption></figure>

possède 17 liaisons, 30 angles de valence et 39 angles de torsion.

On va maintenant décrire plus en détail les constituants principaux du champ de forces OPLS.

## 2.1 Liaisons

Une liaison (*bond*) modélise un partage d’électrons entre deux atomes, engendrant une force entre ceux-ci. En OPLS, les potentiels des liaisons sont
harmoniques : un *potentiel harmonique de liaison* est un champs scalaire $\mathcal{U}$ définissant l'énergie (potentielle) de liaison entre deux atomes placés à distance $r$ comme étant :

$$ \mathcal{U}(r) = k(r - r_0)^2 \qquad (2.1) $$

où $k$ est la raideur de la liaison et $r_0$ est la distance d'équilibre (distance à laquelle la force exercée sur les deux atomes est nulle).

La Fig.2.3 montre les courbes paraboliques des potentiels OPLS de liaison entre deux carbones (CC) et entre un carbone et un hydrogène (CH). Les

<figure><img src="image_4.png" /><figcaption>FIG. 2.3 – Potentiels OPLS (échelle AA) de liaisons entre deux atomes : CC deux carbones, CH un carbone et un hydrogène.</figcaption></figure>

énergies sont données dans les unités internes du système de DM (cf. Chap.5).

## 2.2 Angles de valence

Les angles de valence tendent à maintenir constant un angle $\theta$ entre trois atomes liés. En OPLS, les potentiels des angles de valence sont harmoniques : un *potentiel harmonique de valence* est un champ scalaire $\mathcal{U}$ qui définit l'énergie potentielle d'un angle $\theta$ par :

$$ \mathcal{U}(\theta) = k(\theta - \theta_0)^2 \qquad (2.2) $$
où $k$ est la raideur de l'angle et $\theta_0$ est l'angle d'équilibre (celui pour lequel aucune force n'est exercée sur les trois atomes par l'angle de valence).

La Fig.2.4 montre les courbes paraboliques des potentiels OPLS de valence entre trois carbones (CCC), entre un carbone et deux hydrogènes (HCH) et entre deux carbones et un hydrogène (CCH).

<figure><img src="image_5.png" /><figcaption>FIG. 2.4 – Potentiels OPLS des angles de valence entre trois atomes : CCC trois carbones, HCH un carbone et deux hydrogènes, CCH deux carbones et un hydrogène.</figcaption></figure>

## 2.3 Angles de torsion

Un angle de torsion (aussi appelé dièdre) tend à maintenir constant l'angle formé entre deux plans déterminés par quatre atomes liés. En OPLS, les potentiels des angles de torsion ont une forme en “triple cosinus” ce qui signifie que le potentiel $U$ d'un angle de torsion $\theta$ est donné par :

$$U(\theta) = 0.5[A_1(1 + \cos(\theta)) + A_2(1 - \cos(2\theta)) + A_3(1 + \cos(3\theta))] \quad (2.3)$$

La Fig.2.5 montre les courbes des potentiels OPLS de torsion entre quatre carbones (CCCC), entre deux carbone et deux hydrogènes (HCCH) et entre trois carbones et un hydrogène (CCCH). Il est à noter que pour tous les angles de torsion, les deux atomes centraux sont des carbones. Par ailleurs,
en OPLS, aucune force n'est exercée sur les quatre atomes lorsqu'ils sont dans un même plan (angle de torsion de $\pi$).

<figure><img src="image_6.png" /><figcaption>FIG. 2.5 – Potentiels OPLS (échelle AA) d'angles de torsion entre quatre atomes : CCCC quatre carbones, CCCH trois carbone et un hydrogène, HCCH deux carbones et deux hydrogènes.</figcaption></figure>

## 2.4 Forces inter-moléculaires

Les forces de van Der Waals qui s'exercent entre deux atomes sont extrêmement répulsives à faible distance et faiblement attractives à longue distance. En OPLS, les forces de van Der Waals sont décrites par des potentiels de Lennard-Jones 6-12 : un potentiel de Lennard-Jones 6-12 est défini par deux paramètres $\sigma$ et $\epsilon$; $\sigma$ est la distance à laquelle le potentiel est nul et $\epsilon$ est la profondeur du potentiel (maximum de l'énergie attractive). L'énergie potentielle $U(r)$ entre deux atomes à distance $r$ est définie par :

$$U(r) = 4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right] \quad (2.4)$$

La Fig.2.6 montre les courbes des potentiels OPLS de Lennard-Jones entre deux carbones (CC), entre un carbone et un hydrogène (CH) et entre deux hydrogènes (HH).
<figure><img src="image_7.png" /><figcaption>FIG. 2.6 – Potentiel Lennard-Jones OPLS entre deux atomes à l'échelle AA : CC deux carbones, CH un carbone et un hydrogène, HH deux hydrogènes.</figcaption></figure>
# Chapitre 3

## Application des forces

On décrit maintenant les forces à l'échelle $AA^1$. La définition des forces qui s'appliquent sur les atomes doit être très précise sous peine d'introduire ou de perdre de l'énergie lors des simulations de systèmes moléculaires clos. On utilise les notations d'algèbre vectorielle suivantes :

* On note $\overrightarrow{ab}$ le vecteur dont l'atome $a$ est l'origine et l'atome $b$ l'arrivée ; la distance entre les deux atomes est notée $|ab|$.

* Le vecteur nul est noté 0.

* La longueur du vecteur $\overrightarrow{u}$ est notée $|\overrightarrow{u}|$. On a donc $|\overrightarrow{ab}| = |ab|$.

* La multiplication du vecteur $\overrightarrow{u}$ par le scalaire $n$ est notée $n \cdot \overrightarrow{u}$, ou plus simplement $n \overrightarrow{u}$.

* Le produit vectoriel de $\overrightarrow{u}$ avec $\overrightarrow{v}$ est notée $\overrightarrow{u} \times \overrightarrow{v}$.

* Le produit scalaire de $\overrightarrow{u}$ avec $\overrightarrow{v}$ est noté $\overrightarrow{u} \cdot \overrightarrow{v}$.

* On écrit $\overrightarrow{u} \perp \overrightarrow{v}$ lorsque $\overrightarrow{u}$ et $\overrightarrow{v}$ sont orthogonaux ($\overrightarrow{u} \cdot \overrightarrow{v} = 0$).

* On note $norm(\overrightarrow{u})$ le vecteur normalisé de $\overrightarrow{u}$ (même direction, mais longueur égale à 1) défini par $norm(\overrightarrow{u}) = (1/|\overrightarrow{u}|) \cdot \overrightarrow{u}$.

* Lorsque $a$, $b$ et $c$ sont des atomes, on note $abc$ l'angle formé par $a$, $b$ et $c$.

### 3.1 Liaisons

Une liaison (bond) entre deux atomes engendre une force qui est la dérivée du potentiel de liaison défini entre les deux atomes. La figure 3.1 représente une force (attractive) engendrée entre deux atomes liés $a$ et $b$.

<footnote n="1">¹Ce chapitre est directement tiré de [17].</footnote>
<figure><img src="image_8.png" /><figcaption>FIG. 3.1 – Liaison attractive entre deux atomes liés</figcaption></figure>

Le potentiel $\mathcal{U}$ d'une liaison entre deux atomes est harmonique et vérifie l'équation Eq.2.1. La dérivée partielle de $\mathcal{U}$ par rapport à la distance est :

$$ \frac{\partial \mathcal{U}(r)}{\partial r} = 2k(r - r_0) \qquad (3.1) $$

Soient deux atomes $a$ et $b$, et $\vec{u} = \text{norm}(\vec{ba})$ la normalisation du vecteur $\vec{ba}$. La force engendrée sur l'atome $a$ est :

$$ \vec{f}_a = -2k(r - r_0) \cdot \vec{u} \qquad (3.2) $$

et celle sur $b$ est l'opposée :

$$ \vec{f}_b = -\vec{f}_a \qquad (3.3) $$

Ainsi, si $r > r_0$, la force sur $a$ est un vecteur de direction opposée à $\vec{u}$ et tend donc à rapprocher $a$ de $b$ (force attractive), tandis qu'elle tend à les éloigner (force répulsive) lorsque $r < r_0$. On note que la somme des deux forces est nulle :

$$ \vec{f}_a + \vec{f}_b = 0 \qquad (3.4) $$

## 3.2 Angles de valence

Un angle de valence tend à maintenir constant l'angle $\theta$ entre trois atomes $a$, $b$ et $c$ tels que $a$ soit lié à $b$ et $b$ à $c$, comme sur la figure 3.2. Les forces appliquées aux trois atomes appartiennent toujours au plan qu'ils définissent.

Le potentiel $\mathcal{U}$ d'un angle de valence est harmonique et vérifie l'équation Eq.2.2.

La dérivée partielle de $\mathcal{U}$ par rapport à l'angle $\theta$ est donc :

$$ \frac{\partial \mathcal{U}(\theta)}{\partial \theta} = 2k(\theta - \theta_0) \qquad (3.5) $$
<figure><img src="image_9.png" /><figcaption>FIG. 3.2 – Angle de valence</figcaption></figure>

La dérivée partielle de $\mathcal{U}$ par rapport à la position $r_a$ de $a$ est :

$$ \frac{\partial \mathcal{U}(\theta)}{\partial r_a} = \frac{\partial \mathcal{U}(\theta)}{\partial \theta} \cdot \frac{\partial \theta}{\partial r_a} \qquad (3.6) $$

c'est-à-dire :

$$ \frac{\partial \mathcal{U}(\theta)}{\partial r_a} = 2k(\theta - \theta_0) \cdot \frac{\partial \theta}{\partial r_a} \qquad (3.7) $$

Puisque $a$ décrit le cercle de rayon $|ab|$, centré en $b$, on $a^2$ :

$$ \frac{\partial \theta}{\partial r_a} = \frac{1}{|ab|} \qquad (3.8) $$

Soit $\vec{p}_a$ le vecteur normalisé du plan $abc$, orhogonal à $\vec{ba}$. La force appliquée en $a$ est :

$$ \vec{f}_a = -2k(\theta - \theta_0)/|ab| \cdot \vec{p}_a \qquad (3.9) $$

De même, la force appliquée en $c$ est :

$$ \vec{f}_c = -2k(\theta - \theta_0)/|bc| \cdot \vec{p}_c \qquad (3.10) $$

où $\vec{p}_c$ le vecteur normalisé du plan $abc$, orthogonal à $\vec{cb}$.

Enfin, la force appliquée à $b$ est :

$$ \vec{f}_b = -(\vec{f}_a + \vec{f}_c) \qquad (3.11) $$

On note que la somme des forces est nulle :

$$ \vec{f}_a + \vec{f}_b + \vec{f}_c = 0 \qquad (3.12) $$

<footnote n="2">²La longueur d'un arc de cercle est égale au produit du rayon par l'angle qui sous-tend l'arc.</footnote>
De plus, les deux moment exercés sur b par $\vec{f}_a$ et $\vec{f}_c$ sont inverses puisque :

$$ \vec{ab} \times \vec{f}_a = -\vec{cb} \times \vec{f}_d \qquad (3.13) $$

Aucune rotation autour de b n'est donc produite par $\vec{f}_a$ et $\vec{f}_c$.

## 3.3 Angles de torsion

Le potentiel d'un angle de torsion entre quatre atomes est en "triple-cosinus". Les forces sont représentées sur la Fig.3.3.

<figure><img src="image_10.png" /><figcaption>FIG. 3.3 – Angle de torsion $\theta$</figcaption></figure>

Le potentiel $\mathcal{U}$ d'un angle de torsion vérifie l'équation Eq.2.3. La dérivée partielle de $\mathcal{U}$ par rapport à la position $r_a$ de $a$ est :

$$ \frac{\partial \mathcal{U}(\theta)}{\partial r_a} = \frac{\partial \mathcal{U}(\theta)}{\partial \theta} \cdot \frac{\partial \theta}{\partial r_a} \qquad (3.14) $$

La dérivée partielle du potentiel par rapport à l'angle $\theta$ est :

$$ \frac{\partial \mathcal{U}(\theta)}{\partial \theta} = 0.5(-A_1 \sin(\theta) + 2A_2 \sin(2\theta) - 3A_3 \sin(3\theta)) \qquad (3.15) $$

$$ = -0.5(A_1 \sin(\theta) - 2A_2 \sin(2\theta) + 3A_3 \sin(3\theta)) \qquad (3.16) $$
Forces en a et d

Appelons $\theta_1$ l'angle $\widehat{abc}$. L'atome $a$ tourne autour de la direction $bc$, sur un cercle de rayon $|ab|\sin(\theta_1)$. La dérivée partielle de $\theta$ par rapport à la position de $a$ est :

$$ \frac{\partial \theta}{\partial r_a} = \frac{1}{|ab|\sin(\theta_1)} \qquad (3.17) $$

Et l'on a donc :

$$ \frac{\partial \mathcal{U}(\theta)}{\partial r_a} = \frac{-0.5}{|ab|\sin(\theta_1)} (A_1\sin(\theta) - 2A_2\sin(2\theta) + 3A_3\sin(3\theta)) \quad (3.18) $$

De même, pour l'atome $d$, en notant $\theta_2$ l'angle $\widehat{bcd}$ :

$$ \frac{\partial \mathcal{U}(\theta)}{\partial r_d} = \frac{-0.5}{|cd|\sin(\theta_2)} (A_1\sin(\theta) - 2A_2\sin(2\theta) + 3A_3\sin(3\theta)) \quad (3.19) $$

Soit $\vec{p}_1$ le vecteur unitaire orthogonal au plan $abc$ et $\vec{p}_2$ le vecteur unitaire orthogonal au plan $bcd$ (l'angle entre $\vec{p}_1$ et $\vec{p}_2$ est $\theta$).

La force appliquée en $a$ est :

$$ \vec{f}_a = \frac{0.5}{|ab|\sin(\theta_1)} (A_1\sin(\theta) - 2A_2\sin(2\theta) + 3A_3\sin(3\theta)).\vec{p}_1 \quad (3.20) $$

De même, la force appliquée en $d$ est :

$$ \vec{f}_d = \frac{0.5}{|cd|\sin(\theta_2)} (A_1\sin(\theta) - 2A_2\sin(2\theta) + 3A_3\sin(3\theta)).\vec{p}_2 \quad (3.21) $$

Forces en b et c

On cherche maintenant à déterminer les forces $\vec{f}_b$ et $\vec{f}_c$ à appliquer à $b$ et $c$. On a deux contraintes : (A) la somme des forces doit être nulle :

$$ \vec{f}_a + \vec{f}_b + \vec{f}_c + \vec{f}_d = 0 \qquad (3.22) $$

et (B) la somme des moments des couples exercés doit également être nulle<fnref n="3" />. En appelant $o$ le centre de la liaison $bc$, on doit avoir :

$$ \overrightarrow{oa} \times \vec{f}_a + \overrightarrow{oa} \times \vec{f}_d + \overrightarrow{o_b} \times \vec{f}_b + \overrightarrow{o_c} \times \vec{f}_c = 0 \quad (3.23) $$

<footnote n="3">³ On ne peut définir simplement $\vec{f}_b = -\vec{f}_a$ et $\vec{f}_c = -\vec{f}_d$ car on introduirait alors des couples qui entraîneraient un accroissement d'énergie potentielle.</footnote>
De (3.23), on tire :

$$ (\vec{ob} + \vec{ba}) \times \vec{f}_a + (\vec{oc} + \vec{cd}) \times \vec{f}_d + \vec{ob} \times \vec{f}_b + \vec{oc} \times \vec{f}_c = 0 \quad (3.24) $$

et :

$$ (-\vec{oc} + \vec{ba}) \times \vec{f}_a + (\vec{oc} + \vec{cd}) \times \vec{f}_d - \vec{oc} \times \vec{f}_b + \vec{oc} \times \vec{f}_c = 0 \quad (3.25) $$

d'où :

$$ \vec{oc} \times (-\vec{f}_a + \vec{f}_d - \vec{f}_b + \vec{f}_c) + \vec{ba} \times \vec{f}_a + \vec{cd} \times \vec{f}_d = 0 \quad (3.26) $$

De (3.22), on tire :

$$ - \vec{f}_a + \vec{f}_d - \vec{f}_b + \vec{f}_c = 2(\vec{f}_d + \vec{f}_c) \quad (3.27) $$

On doit donc avoir :

$$ \vec{oc} \times (2(\vec{f}_d + \vec{f}_c)) + \vec{ba} \times \vec{f}_a + \vec{cd} \times \vec{f}_d = 0 \quad (3.28) $$

donc :

$$ 2\vec{oc} \times \vec{f}_d + 2\vec{oc} \times \vec{f}_c + \vec{ba} \times \vec{f}_a + \vec{cd} \times \vec{f}_d = 0 \quad (3.29) $$

et :

$$ (2\vec{oc} + \vec{cd}) \times \vec{f}_d + 2\vec{oc} \times \vec{f}_c + \vec{ba} \times \vec{f}_a = 0 \quad (3.30) $$

d'où :

$$ 2\vec{oc} \times \vec{f}_c = -2\vec{oc} \times \vec{f}_d - \vec{cd} \times \vec{f}_d - \vec{ba} \times \vec{f}_a \quad (3.31) $$

En posant :

$$ \overrightarrow{cpl} = -(\vec{oc} \times \vec{f}_d + 0.5\vec{cd} \times \vec{f}_d + 0.5\vec{ba} \times \vec{f}_a) \quad (3.32) $$

on doit donc avoir :

$$ \vec{oc} \times \vec{f}_c = \overrightarrow{cpl} \quad (3.33) $$

L'équation $\vec{oc} \times x = \overrightarrow{cpl}$ a une infinité de solutions en $x$; on en choisit une
penpendiculaire à $\vec{oc}$ :

$$ \vec{f}_c \perp \vec{oc} \quad (3.34) $$

$$ \vec{f}_c = (1/|oc|^2) \overrightarrow{cpl} \times \vec{oc} \quad (3.35) $$

L'équation (3.33) est bien vérifiée car :
$$ \overrightarrow{\sigma c} \times \overrightarrow{f_c} = (1/|oc|^2) \overrightarrow{\sigma c} \times (\overrightarrow{cpl} \times \overrightarrow{\sigma c}) \quad (3.36) $$

donc<fnref n="4" /> :

$$ \overrightarrow{\sigma c} \times \overrightarrow{f_c} = (1/|oc|^2)|oc|^2 \overrightarrow{cpl} = \overrightarrow{cpl} \quad (3.37) $$

La valeur de $\overrightarrow{f_b}$ est finalement déduite de (3.22) :

$$ \overrightarrow{f_b} = -\overrightarrow{f_a} - \overrightarrow{f_c} - \overrightarrow{f_d} \quad (3.38) $$

## 3.4 Forces inter-moléculaires

Les potentiels de Lennard-Jones ont des formes données par l'équation (2.4). En posant $A = \sigma^{12}$ et $B = \sigma^6$, cette équation s'écrit :

$$ U(r) = 4\epsilon\left(\frac{A}{r^{12}} - \frac{B}{r^6}\right) \quad (3.39) $$

La dérivée partielle de $U$ par rapport à la distance est donc :

$$ \frac{\partial U(r)}{\partial r} = 4\epsilon(-12\frac{A}{r^{13}} + 6\frac{B}{r^7}) \quad (3.40) $$

$$ = 24\epsilon(-2\frac{A}{r^{13}} + \frac{B}{r^7}) \quad (3.41) $$

$$ = \frac{24\epsilon}{r}(-2\frac{A}{r^{12}} + \frac{B}{r^6}) \quad (3.42) $$

$$ = -\frac{24\epsilon}{r}(2(\frac{\sigma}{r})^{12} - (\frac{\sigma}{r})^{6}) \quad (3.43) $$

Soient deux atomes $a$ et $b$. La force engendrée sur l'atome $a$ est :

$$ \overrightarrow{f_a} = \frac{24\epsilon}{r} \left( 2\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} \right) \cdot \overrightarrow{u} \quad (3.44) $$

où $\overrightarrow{u}$ est la normalisation de $\overrightarrow{ba}$. La force sur $b$ est l'opposée de celle sur $a$ :

$$ \overrightarrow{f_b} = -\overrightarrow{f_a} \quad (3.45) $$

La somme des deux forces exercées sur $a$ et $b$ est donc nulle.

<footnote n="4">⁴ Si $u \perp v$, alors $u \times (v \times u) = |u|^2 v$.</footnote>
Résumé

Les résultats précédents sont concentrés dans le tableau suivant :

<table><tbody><tr><td>Liaison</td><td>3.2</td><td>$\vec{f}_a = -2k(r - r_0) \cdot \vec{u}$</td></tr><tr><td>ab</td><td>3.3</td><td>$\vec{f}_b = -\vec{f}_a$</td></tr><tr><td>Valence</td><td>3.9</td><td>$\vec{f}_a = -2k(\theta - \theta_0) / |ab| \cdot \vec{p}_a$</td></tr><tr><td>abc</td><td>3.11</td><td>$\vec{f}_b = -(\vec{f}_a + \vec{f}_c)$</td></tr><tr><td></td><td>3.10</td><td>$\vec{f}_c = -2k(\theta - \theta_0) / |bc| \cdot \vec{p}_c$</td></tr><tr><td>Torsion</td><td>3.20</td><td>$\vec{f}_a = \frac{0.5}{|ab|\sin(\theta_1)} (A_1\sin(\theta) - 2A_2\sin(2\theta) + 3A_3\sin(3\theta)) \cdot \vec{p}_1$</td></tr><tr><td>abcd</td><td>3.38</td><td>$\vec{f}_b = -\vec{f}_a - \vec{f}_c - \vec{f}_d$</td></tr><tr><td></td><td>3.35</td><td>$\vec{f}_c = (1/|oc|^2) \vec{cpl} \times \vec{o\dot{c}}$</td></tr><tr><td></td><td>3.21</td><td>$\vec{f}_d = \frac{0.5}{|cd|\sin(\theta_2)} (A_1\sin(\theta) - 2A_2\sin(2\theta) + 3A_3\sin(3\theta)) \cdot \vec{p}_2$</td></tr><tr><td>L-J</td><td>3.44</td><td>$\vec{f}_a = \frac{24\epsilon}{r} (2(\frac{\sigma}{r})^{12} - (\frac{\sigma}{r})^6) \cdot \vec{u}$</td></tr><tr><td>ab</td><td>3.45</td><td>$\vec{f}_b = -\vec{f}_a$</td></tr></tbody></table>

* Dans l'équation 3.2, $k$ est la raideur de la liaison, $r$ est la distance entre les atomes $a$ et $b$, et $r_0$ est la distance de consigne, pour laquelle l'énergie s'annule. Le vecteur $\vec{u}$ est défini par $\vec{u} = \text{norm}(\vec{ba})$.

* Dans les équations 3.9 et 3.10, $k$ est la raideur de l'angle, $\theta$ est l'angle $\widehat{abc}$, et $\theta_0$ est l'angle de consigne, pour lequel l'énergie s'annule. Dans 3.9, le vecteur $\vec{p}_a$ est défini par $\vec{p}_a = \text{norm}(\vec{ba} \times (\vec{ba} \times \vec{bc}))$. Dans 3.10, le vecteur $\vec{p}_c$ est défini par $\vec{p}_c = \text{norm}(\vec{cb} \times (\vec{ba} \times \vec{bc}))$.

* Dans les équations 3.20 et 3.21, $\theta$ est l'angle de torsion, $\theta_1$ est l'angle $\widehat{abc}$, $\theta_2$ est l'angle $\widehat{bcd}$ et $A_1, A_2$ et $A_3$ sont les paramètres qui définissent la forme en 3-cosinus du dièdre. Le vecteur $\vec{p}_1$ est défini par $\vec{p}_1 = \text{norm}(\vec{ba} \times \vec{bc})$ et $\vec{p}_2 = \text{norm}(\vec{cd} \times \vec{cb})$. Dans 3.35, $o$ est le milieu de $bc$ et $\vec{cpl}$ est défini par $\vec{cpl} = -(\vec{o\dot{c}} \times \vec{f}_d + 0.5\vec{cd} \times \vec{f}_d + 0.5\vec{ba} \times \vec{f}_a)$.

* Dans l'équation 3.44, $\sigma$ est la distance à laquelle le potentiel est nul et $\epsilon$ est la profondeur du potentiel (maximum de l'énergie attractive). Comme pour les liaisons, on a $\vec{u} = \text{norm}(\vec{ba})$.

Pour chaque cas (liaison, valence, torsion, Lennard-Jones), la somme des forces appliquées aux atomes est toujours nulle (Eq. (3.4), (3.12), (3.38), (3.45)).

De plus, aucun couple de forces ne résulte de l'application de ces forces : aucun couple ne peut résulter de l'application des forces dans le cas des liaisons ou des forces inter-moléculaires de Lennard-Jones car elles sont coli-
néaires; en ce qui concerne les angles de valence, nous avons vérifié l'absence
de couple dans la Sec. 3.2; pour les angles de torsion, nous avons choisi des
forces telles leur somme soit nulle et que la somme des couples induits le soit
également (3.38). En conséquence, on est sûr qu'aucune énergie supplémen-
taire ne risque d'être ajoutée à un système clos lors des simulations.

Il est à noter que ce sont les dièdres qui entraînent les changements de géométrie 3D des molécules, tous les autres composants produisant des forces restant systématiquement dans un même plan.
# Chapitre 4

## Implémentation

On considère dans ce chapitre la question de l'implémentation de la DM en décrivant un système implémenté dans le langage de programmation Java, avec la librairie de visualisation 3D, Java3D. Ce système a été présenté dans [6] duquel on reprend ici les principaux éléments. L'implémentation est un prototype qui ne prend pas en compte plusieurs fonctionnalités offertes généralement par les systèmes de DM, comme le contrôle de température (thermostat) ou de pression (barostat) ou encore la possibilité de définir des systèmes moléculaires par des conditions périodiques (cristaux). En fait, le système que l'on va considérer implémente uniquement le noyau de la DM, c'est-à-dire la mécanique newtonienne, et seuls les systèmes moléculaires composés de chaînes linéaires d'atomes de carbone et d'hydrogène (alkanes) introduits dans le Chap. 2 seront considérés.

L'objectif principal poursuivi est de fournir une implémentation de systèmes moléculaires multi-échelles (cf 1.2), permettant des changements d'échelle de description en cours de simulation (1.3). Sur le plan informatique, cela signifie que l'implémentation doit pouvoir traiter des systèmes moléculaires non définis une fois pour toute, non-figés dès le départ, mais au contraire doit permettre des définitions *modulaires* dans lesquelles des molécules peuvent être créées ou disparaître en cours d'exécution.

Un aspect central de la physique newtonienne sur laquelle la DM est fondée est le *déterminisme* ou, ce qui revient au même, la préservation de l'énergie d'un système clos, au cours du temps. Il est évidemment impératif que l'énergie totale d'un système moléculaire clos reste la même au cours de la simulation, y compris si des changements d'échelle ont lieu.

Le système de DM qui a été construit et qui est utilisé ici est fondé sur
un paradigme de programmation, appelé "programmation réactive"[5] qui permet de concilier modularité et déterminisme. L'outil central de modularité de la programmation réactive est le *parallélisme déterministe* que l'on utilise à travers une librairie Java appelée *SugarCubes*[7].

Dans la suite de ce chapitre, on décrit la programmation réactive, et la librairie *SugarCubes* qui l'implémente en Java. On présente ensuite la méthode de résolution utilisée pour obtenir des simulations stables au cours du temps (La stabilité de l'implémentation est discutée au Chap. 5).

## 4.1 Programmation Réactive

La programmation réactive (RP) propose un paradigme de programmation simple avec une sémantique claire et précise. La caractéristique centrale de RP est de fournir des primitives permettant d'exprimer le parallélisme directement au niveau de la programmation. Il s'agit en RP d'un *parallélisme logique* à bien distinguer du parallélisme d'exécution lié au système d'exploitation sur lequel les simulations sont exécutées. Le parallélisme logique est un outil extrêmement puissant de modularité, permettant de décomposer et de coder les systèmes complexes en sous-systèmes communiquants dont la structure peut évoluer dynamiquement (dynamicité). Le parallélisme logique de RP possède de plus une caractéristique fondamentale : il est déterministe. Concilier parallélisme, déterminisme et dynamicité peut sembler a-priori paradoxal, cependant RP donne un moyen de résoudre ce paradoxe dans un cadre informatique cohérent.

Dans l'approche réactive, les systèmes sont composés de composants parallèles partageant les mêmes instants qui définissent ainsi une *horloge logique* partagée par tous les composants. Les composants se synchronisent à la fin de chaque instant et s'exécutent ainsi au même rythme. Durant chaque instant, les composants peuvent communiquer entre eux par des signaux (appelés *événements*) diffusés instantanément. Ces signaux sont analogues à des messages radio, immédiatement reçus par tous les récepteurs à l'écoute, ceux-ci recevant tous, comme avec la radio, le même message. Dans l'approche réactive, la dynamicité n'est prise en compte qu'aux frontières d'instants.

Il existe plusieurs variantes de RP, étendant divers langages de programmation généralistes (par exemple, ReactiveC [5] qui étend C, et ReactiveML [14] qui étend le langage ML). RP est également fortement reliée au langage de programmation synchrone *Esterel*[4], la principale différence tenant
en ce que l'évolution dynamique des programmes n'est pas autorisée en Es-
terel alors qu'elle est permise en RP. L'une des variantes de RP qui étend le
langage Java s'appelle SugarCubes[7]. L'opérateur **merge** de parallélisme de
SugarCubes est totalement déterministe, ce qui signifie qu'à chaque instant
un programme SugarCubes a une sortie unique, fonction des entrées, et qu'en
plus la trace d'exécution est unique.

On va maintenant décrire les SugarCubes dans la section suivante<fnref n="1" />.

4.2 SugarCubes

Les deux classes principales des SugarCubes sont la classe Instruction
des instructions réactives définies en référence aux instants, et la classe Mach-
ine des machines réactives qui exécutent les instructions réactives.

Une machine réactive exécute les instructions réactives dont elle a la
charge d'une manière coordonnée, en leur permettant de communiquer à
l'aide d'événements diffusés (*broadcast events*) instantanément. Un événe-
ment généré au cours d'un instant est reçu au cours de ce même instant par
toutes les instructions en attente de cet événement. De plus, toutes les ins-
tructions recevant l'événement le "voient" de la même façon : il est présent si
il est généré au cours de l'instant par une des instructions réactives, et il est
absent sinon. Une conséquence est que la réaction à l'absence d'un événement
ne peut avoir lieu qu'à l'instant suivant. Ainsi, il ne peut y avoir de cas où
un événement serait vu comme présent au cours d'un instant par une des
instructions exécutée par la machine réactive et où il serait vu comme absent
par une autre instruction au cours du même instant.

Les événements sont automatiquement réinitialisés à absent par la ma-
chine réactive au début de chaque nouvel instant. Des valeurs peuvent être
associées aux événements et si c'est le cas, les instructions les recoivent toutes
et dans le même ordre.

Les machines réactives procèdent par instants : l'exécution a lieu durant
le premier instant, puis pendant le deuxième, puis pendant le suivant, et
ainsi de suite indéfiniment. L'ordre dans lequel les instructions réactives sont
exécutées par la machine est toujours le même. Toutes les instructions pré-
sentes dans la machine réactive sont exécutées à chaque instant et le passage
à l'instant suivant n'a lieu que lorsqu'il ne reste plus rien à exécuter durant

<footnote n="1">¹On présente les principales notions, une description complète étant disponible dans
[7].</footnote>
l'instant courant. En particulier, l'exécution dans l'instant est poursuivie s'il reste des instructions réactives bloquées sur des événements générés. Ainsi, on est sûr qu'à la fin de chaque instant toutes les instructions réactives ont fini leur exécution pour cet instant et ont donc bien réagi à la totalité des événements générés.

De nouvelles instructions réactives peuvent être introduite dans la machine réactive uniquement entre deux instants. De cette manière, le déterminisme de l'exécution durant un instant ne risque pas d'être remis en cause lors de l'introduction de nouvelles instructions en cours d'exécution.

Ainsi, les instructions réactives sont exécutées en parallèle et de manière déterministe par la machine réactive. Ce parallélisme déterministe et dynamique est également directement disponible sous forme d'une instruction réactive appelée **merge**.

Les principales instructions réactives des SugarCubes sont les suivantes :

**Instant suivant.** L'instruction *stop* suspend l'exécution pour l'instant courant de l'instruction réactive dans laquelle elle se trouve. L'exécution reprendra après le *stop* à l'instant suivant.

**Séquence.** L'instruction *seq* (i1,i2) se comporte comme i1 et l'exécution passe immédiatement (i.e. sans attendre l'instant suivant) à i2 dès que i1 termine.

**Parallélisme.** L'instruction *merge* (i1,i2) exécute un instant de l'instruction i1 et de l'instruction i2 et termine lorsque i1 et i2 sont elles-mêmes terminées. L'exécution commence toujours par i1 et passe à i2 lorsque i1 termine ou est suspendue.

**Boucle.** L'instruction *loop* (i) exécute *i* en boucle : l'exécution de *i* est immédiatement recommencée dès qu'elle termine. On suppose qu'il est impossible que *i* commence et termine au cours d'un même instant (autrement, on aurait une *boucle instantanée* qui pourrait cycler indéfiniment, en empêchant la machine réactive de détecter la fin de l'instant courant et donc de passer à l'instant suivant).

**Code Java.** L'instruction *action* (*act*) exécute la méthode *execute* de l'objet *Java act* (du type *JavaAction*) et termine immédiatement. L'exécution
tion de act peut avoir lieu plus d'une fois si l'instruction apparaît dans une boucle.

**Génération d'événement.** L'instruction `generate (event,value)` génère `event` avec `value` comme valeur associée, puis termine immédiatement.

**Attente d'événement.** L'instruction `await (event)` termine immédiatement si l'événement `event` est présent (i.e. a été généré précédemment durant l'instant courant). Sinon la même exécution recommence à l'instant suivant.

**Valeurs générées.** L'instruction `callback (event,call)` exécute la méthode `execute` de l'objet Java `call` (du type `JavaCallback`) pour chaque valeur de l'événement `event` générée durant l'instant courant. Pour ne pas risquer de perdre des valeurs, l'exécution se suspend après avoir traité chaque valeur et, dans tous les cas, ne termine qu'à l'instant suivant.

**Préemption.** L'instruction `until (event,i)` exécute *i* et termine soit parce que *i* termine, soit parce que l'événement `event` est présent.

Les instructions de séquence et de parallélisme sont naturellement éten-
dues à plus de deux branches : par exemple, l'instruction `seq (i1,i2,i3)`
met en séquence les trois instruction `i1,i2,i3`, avec pour signification `seq
(i1,seq (i2,i3))`.

Une machine réactive de la classe Machine exécute son programme qui est une instruction réactive. Les nouvelles instructions ajoutées dans la machine sont mises en parallèle (merge) avec le programme existant. Le rôle de la machine réactive est d'exécuter son programme, de détecter la fin de l'instant courant, c'est-à-dire quand toutes les branches parallèles de son programme sont soit terminées, soit suspendues, et lorsque c'est le cas, de passer à l'instant suivant. L'ajout de nouvelles instructions ne peut avoir lieu durant l'instant courant, mais uniquement lorsque celui-ci est terminé, avant le début de l'instant suivant.

L'exécution du programme durant un instant peut se décomposer en plusieurs phases successives, durant lesquelles de nouveaux événements sont générés ce qui a pour effet de faire réagir les instructions qui les attendent. Par exemple, considérons le code suivant, en supposant que l'événement e n'est pas généré :
```
merge (
    await e,
    generate (e,0)
)
```

L'exécution commence par l'instruction `await` (ligne 2) qui se suspend puisque `e` est absent. L'exécution passe alors à l'instruction `generate` (ligne 3) qui génère `e` (avec la valeur 0) et qui termine. La première phase d'exécution de l'instruction `merge` est terminée mais la machine réactive détecte qu'il existe une instruction suspendue en attente de `e` alors que celui-ci est présent. Une nouvelle phase débute alors et l'instruction `await` est à nouveau exécutée. Elle termine alors puisque `e` est maintenant présent. L'instruction `merge` termine également car ses deux branches sont terminées.

L'exécution du programme par une machine réactive est totalement déterministe : une seule trace d'exécution est possible. L'exécution des programmes `SugarCubes` est séquentielle et le parallélisme est purement logique. La définition précise (formelle) de la sémantique des `SugarCubes` et une comparaison avec l'approche standard par `threads` de la concurrence est disponible dans [8].

## 4.3 Méthode de résolution

La méthode de résolution que l'on utilise est connue sous le nom "velocity-Verlet"[23]. Pour chaque atome, elle fonctionne en deux étapes : la position de l'atome est déterminée à la fin de la première étape, et sa vitesse à la fin de la seconde. L'accélération est obtenue (loi de Newton) lors de la seconde étape en sommant toutes les forces appliquées à l'atome pendant cette étape. A la fin de la seconde étape, la vitesse et l'accélération sont mémorisées pour être utilisées à l'étape suivante.

On va décrire plus précisemment le fonctionnement de velocity-Verlet.

Soit $\mathbf{r}_t$ la position au temps $t$ d'un atome, $\mathbf{v}_t$ sa vitesse au temps $t$, et $\mathbf{a}_t$ son accélération au temps $t$. La méthode de résolution *Velocity-Verlet* est définie par les équations suivantes, où $\Delta t$ est un intervalle de temps :

$$ \mathbf{r}_{t+\Delta t} = \mathbf{r}_t + \mathbf{v}_t \Delta t + 0.5 \mathbf{a}_t \Delta t^2 \quad (4.1) $$

$$ \mathbf{v}_{t+\Delta t} = \mathbf{v}_t + 0.5(\mathbf{a}_t + \mathbf{a}_{t+\Delta t})\Delta t \quad (4.2) $$

L'accélération $\mathbf{a}_{t+\Delta t}$ d'un atome est obtenue en collectant les forces intra et inter-moléculaires qui s'exercent sur celui-ci lorsque sa position est $\mathbf{r}_{t+\Delta t}$.
Pour cela, deux étapes sont nécessaires : dans la première, l'atome est positionné et dans la seconde les forces qui s'exercent sur lui sont collectées; l'algorithme est le suivant :

1. Calcul de la vitesse à la moitié du pas de temps, à partir de la position et de l'accélération au pas-de-temps précédent :

$$v_{t+0.5\Delta t} = v_t + 0.5a_t\Delta t \quad (4.3)$$

Puis, utilisation du résultat pour calculer la position à la fin du pas-de-temps complet :

$$r_{t+\Delta t} = r_t + v_{t+0.5\Delta t}\Delta t \quad (4.4)$$

2. Calcul de l'accélération à la fin du pas-de-temps complet à partir des forces appliquées à l'atome, et calcul de la vitesse à la fin du pas-de-temps complet, en utilisant la vitesse précédemment calculée pour la moitié du pas-de-temps :

$$v_{t+\Delta t} = v_{t+0.5\Delta t} + 0.5a_{t+\Delta t}\Delta t \quad (4.5)$$

Les instants de la RP sont naturellement identifiés aux étapes de la résolution. Ainsi, deux instants réactifs consécutifs sont nécessaires pour implémenter chaque couple d'étape de résolution correspondant à un pas-de-temps $\Delta t$.

Il est à noter que le pas-de-temps peut varier suivant les objets simulés. C'est cette caractéristique qui permet de simuler, au sein d'un même système moléculaire, des molécules à des échelles de description différentes (et par exemple, de faire co-exister dans une même simulation des molécules UA et AA).

## Conclusion

En RP la programmation est à base de parallélisme logique, d'événements diffusés instantanément et de créations/destructions dynamiques d'instructions et d'événements. De plus, dans la version SugarCubes de RP, les programmes sont complètement déterministes "par construction". Le système de dynamique moléculaire réalisé en SugarCubes et utilisé dans ce texte montre que la programmation réactive est un outil intéressant pour implémenter des aspects essentiels de la physique classique. En particulier en ce qui concerne les points suivants :
* Les simulations numériques reposent sur une discrétisation du temps et sur une méthode de résolution implémentant un algorithme d'intégration. La découpe du temps en instants à la base de la programmation réactive discrétise naturellement le temps et permet de coder de façon très simple la méthode de résolution "velocity-Verlet" qui s'avère très efficace dans le contexte de la DM.

* Les forces de la physique classique sont transmises instantanément et s'implémentent naturellement comme des événements diffusés instantanément en RP.

* Les lois de Newton sont déterministes. Leur implémentation est donc grandement facilitée par l'utilisation des SugarCubes dont les programmes sont déterministes "par construction". Les lois de Newton sont également réversibles dans le temps. A l'aide d'une méthode de résolution réversible dans le temps (comme l'est velocity-Verlet), le déterminisme intrinsèque des SugarCubes est également un atout pour obtenir la propriété de réversibilité dans le temps de la physique classique.

* Des changements dans la structure chimique des systèmes simulés peuvent apparaître pendant les simulations. La dynamicité de RP, c'est-à-dire la possibilité de création/destruction d'objet en cours d'exécution, permet d'implémenter ces changements de manière naturelle, sémantiquement claire et préservant le déterminisme.

Le parallélisme logique et les événements diffusés de RP sont un puis-
sant outil de modularité permettant d'introduire de nouveaux composants
parallèles dans un système sans avoir à adapter les autres composants. Cette
modularité est effective aussi bien au niveau du codage qu’à celui de l’exécu-
tion. En particulier, des *observateurs* peuvent être introduits à divers stades
du codage, en parallèle avec le programme à tester, sans que celui-ci n’en soit
perturbé. Cette possibilité permet aux programmeurs d’adopter des straté-
gies de codage modulaire, en accord avec les approches modernes de la pro-
grammation.
# Chapitre 5

## Simulations en AA

On va maintenant simuler les molécules à l'échelle AA avec le potentiel OPLS.

Les formes des composantes du potentiel OPLS concernant les alkanes ont été décrites dans le Chap. 2. Les forces appliquées ont été définies dans le Chap. 4. Pour pouvoir simuler les alkanes à l'échelle AA il reste à donner les valeurs des paramètres des composantes OPLS.

Les unités du système sont les suivantes :

* temps : pico-seconde (ps, $10^{-12}$ seconde);

* distances : nano-mètre (nm, $10^{-9}$ mètre);

* masses : kg/mol;

* angles : radian (rad);

* énergies : kg/mol·nm²/ps²;

* vitesses : nm/ps (également km/s).

Les valeurs des paramètres de OPLS concernant les alkanes sont données, dans les unités système, dans la Fig. 5.1.

Dans la Fig. 5.1 :

* $k_l$ est la raideur de la liaison *l* et $r_l$ sa distance d'équilibre, pour $l = CC$ ou $l = CH$ (Eq. 2.1);

* $k_a$ est la raideur de l'angle de valence *a* et $r_a$ sa valeur d'équilibre, pour $a = CCC$ ou $a = CCH$ ou $a = HCH$ (Eq. 2.2);

* $A1_a$, $A2_a$ et $A3_a$ sont les paramètres des angles de torsion, pour $a = CCCC$ ou $a = CCCH$ ou $a = HCCH$ (Eq. 2.3);

* $\epsilon_p$ et $\sigma_p$ sont les paramètres des paires d'atomes, pour $p = CC, p = CH$ ou $p = HH$ (Eq. 2.4).
$$k_{CC} = 112.1312 \quad r_{CC} = 0.1529$$

$$k_{CH} = 142.256 \quad r_{CH} = 0.109$$

$$k_{CCC} = 0.2441364 \quad \theta_{CCC} = 1.9669860669976096$$

$$k_{CCH} = 0.1569 \quad \theta_{CCH} = 1.8920924724829928$$

$$k_{HCH} = 0.138072 \quad \theta_{HCH} = 1.9288037832146987$$

$$A1_{CCCC} = 0.00728016 \quad A2_{CCCC} = -6.56888.10^{-4}$$

$$A3_{CCCC} = 0.001167336$$

$$A1_{CCCH} = 0 \quad A2_{CCCH} = 0 \quad A3_{CCCH} = 0.001531344$$

$$A1_{HCCH} = 0 \quad A2_{HCCH} = 0 \quad A3_{HCCH} = 0.001330512$$

$$\epsilon_{CC} = 2.7614.10^{-4} \quad \sigma_{CC} = 0.35$$

$$\epsilon_{CH} = 1.86188.10^{-4} \quad \sigma_{CH} = 0.2958$$

$$\epsilon_{HH} = 1.2552.10^{-4} \quad \sigma_{HH} = 0.25$$

FIG. 5.1 – Valeurs des paramètres OPLS concernant les alkanes, exprimées en unités internes.

## 5.1 Stabilité

Le système de DM permet d'effectuer des simulation très longues, préser-
vant l'énergie. Pour illustrer cette stabilité, on considère la simulation d'une
molécule $C_8H_{16}$ à laquelle on a initialement donné de l'énergie en allongeant
légèrement ses liaisons CC. La simulation est effectuée avec un pas-de-temps
de $10^{-4}$ ps (c'est-à-dire 0.1 femto-seconde). La Fig.5.2 montre la molécule en
cours de simulation à gauche, et les énergies mesurées tous les $10^5$ instants.
En bleu, l'énergie cinétique, en vert, l'énergie interne, en rouge l'énergie to-
tale (l'énergie inter-moléculaire, en jaune, est toujours nulle). On constate la
grande stabilité de l'énergie totale.

Les $4.10^7$ instants correspondent à une durée de simulation de 2 ns.

Dans une nouvelle simulation, on considère deux molécules mises en face-
à-face à une distance de 0.9 nm. Les deux molécules sont initialement à
l'équilibre et l'énergie initiale du système est donc la seule énergie (attrac-
tive) inter-moléculaire. La Fig.5.3 montre les deux molécules en cours de
simulation, à gauche, et les énergies mesurées tous les $10^5$ instants, à droite.
<figure><img src="image_11.png" /><figcaption>FIG. 5.2 – À gauche : molécule en cours de simulation. À droite : énergies mesurées tous les $10^5$ instants.</figcaption></figure>

En rouge on a l'énergie totale qui est proche de 0, en jaune l'énergie inter-moléculaire (qui reste toujours négative).

<figure><img src="image_12.png" /><figcaption>FIG. 5.3 – À gauche : deux molécules initialement à l'équilibre en cours de simulation. À droite : énergies mesurées tous les $10^5$ instants.</figcaption></figure>

La Fig.5.4 montre le résultat obtenu en mettant en face-à-face deux molécules auxquelles on a donné une énergie potentielle initiale non nulle, en translatant les atomes vers le haut ($y > 0$).
<figure><img src="image_13.png" /><figcaption>FIG. 5.4 – En haut : deux molécules initialement non à l'équilibre en cours de simulation. En bas : énergies mesurées tous les $10^5$ instants.</figcaption></figure>

## 5.2 Chaos déterministe

Le contexte général des simulations moléculaires est celui du *chaos déterministe* (voir par exemple [13] pour une discussion éclairante de cette notion). Le chaos déterministe est souvent identifié à la sensibilité aux conditions initiales : dans les systèmes manifestant des comportements de chaos déterministe, la précision des conditions initiales n’est jamais suffisante pour empêcher l’apparition du chaos au bout d’un certain temps.

La physique sous-jacente à la Dynamique Moléculaire est la physique classique Newtonienne, qui est parfaitement déterministe et repose sur un ensemble de lois simples et réversibles au cours du temps. La réversibilité peut être simplement exprimée en prenant un pas-de-temps négatif; on obtient alors des simulations tout-à-fait semblables aux simulations avec un pas-de-temps positif. La réversibilité correspond en fait à la préservation de l’énergie totale des systèmes clos au cours du temps.
Le chaos déterministe peut être mis en évidence assez simplement, par un
couple de simulations présentées sur la Fig. 5.5. Dans les deux simulations, on
débranche la partie de l'implémentation correspondant aux forces de van Der
Waals. On se place dans un repère ortho-normé où l'axe des x est horizontal,
l'axe des y est vertical et l'axe des z est perpendiculaire aux deux autres.

On simule tout d'abord deux molécules identiques ($C_{10}H_{22}$), initialement superposées, à une température de 1000 K. La superposition des deux molécules est possible car les forces de van Der Waals ne sont pas prises en compte. Durant la simulation, on constate que la superposition des deux molécules se maintient au cours du temps (testée jusqu'à 5 nano-secondes ; partie centrale de la Fig. 5.5). La préservation de la superposition est le reflet du déterminisme de l'exécution de chacune des deux molécules.

Dans la seconde simulation, on change légèrement les conditions initiales : la coordonnée en y du premier carbone de l'une des deux molécules est très légèrement modifiée, d'une valeur de $10^{-17} nm$. La différence de placement est bien trop faible pour être observée sur la visualisation du système. Cependant, on observe rapidement (au bout de 2.7 pico-secondes ; partie droite de la Fig. 5.5) une nette divergence du positionnement des deux molécules. Cette divergence est la manifestation de la "sensibilité aux conditions initiales" spécifique du chaos déterministe.

<figure><img src="image_14.png" /><figcaption>FIG. 5.5 – À gauche, situation initiale : superposition des deux molécules. Noter que l'énergie initiale est due à la translation en y du premier carbone et de ses hydrogènes liés. Au milieu, situation finale après 5 ns : les deux molécules sont toujours superposées. À droite, après ajout à la translation initiale en y d'une seule molécule d'une quantité de $10^{-17} nm$ : divergence constatée après 2.7 ps.</figcaption></figure>

Un point important est à souligner : le déterminisme total de l'implémen-
tation de la DM est un atout essentiel pour la mise au point des programmes
de simulation. En particulier, un résultat surprenant ou mettant en évidence
une erreur peut toujours être reproduit, ce qui n’est pas le cas avec des im-
plémentations utilisant des *threads* d’exécution ou des communications par
réseau, par nature non-déterministe.
# Chapitre 6

## Approches multi-échelles

Dans ce chapitre on considère l'approche multi-échelles en définissant deux échelles construites à partir de AA. La première échelle est appelée UA et consiste à faire abstraction des atomes d'hydrogène : un grain UA est un atome de carbone avec les atomes d'hydrogène qui lui sont liés. Dans la seconde approche appelée CG on regroupe deux grains UA en un seul grain CG.

### 6.1 Echelle UA

Les molécules à l'échelle UA sont des chaînes de grains de deux sortes : un grain $G_2$ contient un atome de carbone et deux atomes d'hydrogène qui lui sont liés, tandis qu'un grain $G_3$ contient un atome de carbone et trois atomes d'hydrogène qui lui sont liés.

Comme en OPLS à l'échelle AA, les composants intra-moléculaires sont en UA des liaisons, des angles de valence ou des angles de torsion.

Comme on l'a fait pour l'échelle AA, on définit des fragments UA uniquement constitués de grains $G_2$. Dans la suite, pour simplifier, on ne considérera plus que ces fragments UA. Ainsi, en plus des grains $G_2$, on n'aura plus à considérer que des liaisons $G_2G_2$, des angles de valence $G_2G_2G_2$, des dièdres $G_2G_2G_2G_2$ et des forces inter-moléculaires $G_2G_2$.

Dans le système de DM les grains $G_2$ sont représentés par des boules de couleur cyan, comme sur la Fig.6.1 où le fragment considéré est formé de 8 grains $G_2$, de 7 liaisons $G_2G_2$, de 6 angles de valence $G_2G_2G_2$ et de 5 dièdres $G_2G_2G_2G_2$.
<figure><img src="image_15.png" /><figcaption>FIG. 6.1 – Fragment UA composée de 8 grains $G_2$.</figcaption></figure>

En accord avec l'interprétation des grains UA, la masse d'un grain UA est la somme des masses des atomes AA qu'il contient, donc 0.014 pour $G_2$.

Il existe un lien géométrique fort entre les molécules UA et les molécules AA : la molécule UA “équivalente” à une molécule AA est construite en “gommant” tous les hydrogènes, en remplaçant les carbones par des grains UA et les composants AA par les composants UA correspondants. Inversement, la molécule AA équivalente à une molécule UA est construite en remplacant les grains UA par des carbones, en rajoutant des hydrogènes à ces carbones et en remplaçant les composants AA par des composants UA.

La Fig.6.2 montre une molécule $C_8H_{16}$ à l'équilibre à gauche, et la molécule UA équivalente à droite.

Une question se pose : comment définir le potentiel UA à partir du potentiel AA ? La réponse sera donnée en profitant du lien géométrique existant entre les deux échelles au Chap. 11.

## 6.2 Echelle CG

A l'échelle CG, les grains CG sont formés de deux grains UA. On a donc trois sortes de grains CG qui diffèrent suivant le nombre d'hydrogènes qu'ils contiennent. Les grains CG sont désignés par $CG_n$ où n est le nombre d'atomes d'hydrogène ; n vaut donc 4, 5 ou 6. Le grain $CG_4$ est représenté
<figure><img src="image_16.png" /><figcaption>FIG. 6.2 – Molécule (fragment) AA $C_8H_{16}$, à gauche. Molécule UA équivalente, composée de 8 grains $G_2$, à droite.</figcaption></figure>

sur la Fig.6.3

<figure><img src="image_17.png" /><figcaption>FIG. 6.3 – Grain $CG_4$ composé de deux carbones et quatre hydrogènes. Par définition, le centre du grain CG est le milieu $m$ de la liaison entre les deux carbones.</figcaption></figure>

La masse d'un grain CG est la somme des masses de deux carbones, plus celles des hydrogènes qu'il contient. Ainsi, la masse du grain $CG_4$ est 0.028.

Dans la suite, de même qu'à l'échelle UA on ne considère que des grains $G_2$, à l'échelle CG on ne considérera pour simplifier que des grains $CG_4$ (donc formés de deux grains UA $G_2$).

Dans le système de DM, les grains $CG_4$ sont représentés par des boules jaunes. La Fig.6.4 montre une molécule CG formée de 4 grains $CG_4$.
<figure><img src="image_18.png" /><figcaption>FIG. 6.4 – Molécule CG composée de 4 grains $CG_4$.</figcaption></figure>

On établit un lien géométrique avec l'échelle AA en considérant que le grain CG est centré sur le milieu de la liaison AA de la molécule $C_2H_n$ équivalente. Ainsi, avec les notations de la Fig.6.3, la molécule AA équivalente à un grain $CG_4$ contient 5 liaisons ($h_0C_0, h_1C_0, C_0C_1, h_2C_1, h_3C_1$), 6 angles de valence ($h_0C_0h_1, h_0C_0C_1, h_1C_0C_1, h_2C_1C_0, h_3C_1C_0, h_2C_1h_3$) et 4 angles de torsion ($h_0C_0C_1h_2, h_0C_0C_1h_3, h_1C_0C_1h_2, h_1C_0C_1h_3$).

Une molécule UA à l'équilibre (donc, dont l'énergie intra-moléculaire est nulle) a pour équivalent CG une molécule dont les grains $CG_4$ sont tous alignés. La molécule CG de la Fig.6.4 est ainsi un équivalent CG des molécules AA et UA de la Fig.6.2.

A l'équilibre, en CG, les angles de valence sont plats (leur valeur est $\pi$).

Considérons maintenant le cas des angles de torsion à l'échelle CG. À un angle de torsion qui n'est pas à l'équilibre à l'échelle CG doit obligatoirement correspondre au moins un angle de valence CG non nul (sinon, si les deux angles de valence CG étaient nuls, l'angle de torsion CG serait à l'équilibre). Donc, en CG, il est possible de considérer que l'énergie d'un angle de torsion est répartie entre les deux angles de valence CG qu'il contient. Ce choix simplifie la détermination des potentiels CG : les angles de torsion CG ne nécessitent plus de traitement particulier.

On choisit donc de ne définir comme composants des molécules CG que les liaisons CG et les angles de valence CG, et donc de ne pas considérer des angles de torsion à cette échelle. Ainsi, la molécule CG de la Fig.6.4 a uniquement pour composants trois liaisons et deux angles de valence.

Il faut remarquer que l'“absorption” de l'angle de torsion par les angles
de valence qui lui sont associés est possible à l'échelle CG mais ne l'est pas aux échelles AA et UA. En effet, en AA et UA un angle de torsion peut avoir de l'énergie alors que les deux angles de valence associés sont à l'équilibre, ce qui est impossible à l'échelle CG.

De même que pour UA, la question qui se pose est de définir le potentiel CG à partir du potentiel AA, en profitant du lien géométrique existant entre les deux échelles. Cette question sera considérée au Chap. 12.

## 6.3 Complexité

Dans le passage d'une molécule AA à une molécule UA ou CG équivalente, le nombre d'atomes et de composants intra-moléculaires (liaisons, angles de valence, dièdres) diminue nettement. Ainsi, les simulations en UA ou CG ont une complexité plus faible, permettant de simuler des molécules plus longues, avec des ratios temps physique / temps réel plus importants.

La différence de complexité est encore plus manifeste dans le cas des forces inter-moléculaires qui sont en carré du nombre d'atomes.

La Fig.6.5 compare les complexités aux trois échelles. En abscisses, on a la taille *n* des molécules (*n* atomes de carbone pour AA, *n* grains $G_2$ pour UA, et *n*/2 grains $CG_4$ pour CG). Dans la figure du haut, en ordonnées, on a le nombre de composants intra-moléculaires (liaisons, angles de valence et dièdres) qui doivent être simulés à chaque instant. Dans la figure du bas, on a le nombre de paires d'atomes ou de grains qui doivent être analysées pour déterminer les forces inter-moléculaires entre deux molécules de taille *n*.

Pour illustrer le gain de temps réel obtenu avec les échelles UA et CG, on considère une situation où deux molécules sont mises face à face et lancées avec des vitesses opposées, l'une contre l'autre. Les deux molécules vont rebondir l'une contre l'autre et ensuite s'éloigner définitivement.

Les configurations initiales aux trois échelles sont montrées sur la Fig.6.6.

Les trois axes de coordonnées sont présents sur les trois images. À chaque échelle, la molécule de gauche est légèrement décalée en *y* vers le haut (de 0.3 *nm*). Le pas-de-temps est de $10^{-4}$ *ps* pour les trois simulations. La distance initiale entre chaque molécule et le centre des coordonnées est $\Delta = 2$ *nm* et la vitesse initiale de $1$ *nm*/*ps*. Dans chaque simulation, on mesure le temps réel (en milli-seconde) jusqu'à ce que les molécules, après rebond, soient éloignées du centre des coordonnées de plus de $\Delta + 0.5$ *nm*.
<figure><img src="image_19.png" /><figcaption>FIG. 6.5 – Complexités comparées des simulations à diverses échelles. Haut : forces intra-moléculaires ; Bas : forces inter-moléculaires.</figcaption></figure>

Les résultats sont les suivants : pour AA : 13.108 *ms* ; pour UA : 2.004 *ms* ; pour CG : 1.144 *ms*. On a donc un facteur 11 entre AA et CG.
<figure><img src="image_20.png" /><figcaption>FIG. 6.6 – Configurations initiales des simulations aux trois échelles.</figcaption></figure>
# Chapitre 7

## Méthode Boltzmann-inverse

Dans ce chapitre, on décrit une méthode, appelée Boltzmann-inverse, pour déterminer le potentiel à l'échelle UA à partir du potentiel AA. La détermination par la même approche du potentiel CG est considérée dans le Chap. 8.

La méthode Boltzmann-inverse repose sur un traitement statistique de données issues de simulations. Dans cette approche, on détermine l'énergie potentielle d'un oscillateur à partir des densités de probabilité de présence de l'oscillateur dans un état donné.

Pour l'implémentation, on découpe l'espace de variation parcouru en N classes de tailles égales. A chaque classe $C_i$ est associé un compteur $P_i$ qu'on incrémente à chaque étape de simulation lorsque la valeur de l'oscillateur (par exemple, dans le cas d'une liaison, sa longueur) appartient à celle associée à la classe $C_i$. Plus précisemment, la formule définissant l'énergie potentielle $U_i$ associée à la classe $C_i$ est :

$$U_i = -k_B \times T \times \ln(P_i/P_0) \qquad (7.1)$$

où $k_B$ est la constante de Boltzmann, $T$ est la température de simulation, $P_0$ est la valeur du compteur de la classe ayant le maximum d'éléments (la classe la plus "remplie") et $P_i$ la valeur du compteur de la classe $C_i$.

Dans le cas d'un oscillateur harmonique isolé, par exemple, on trouve que la densité de présence maximale correspond à l'énergie cinétique maximale. L'oscillateur étant un système clos, l'état d'énergie cinétique maximale est également l'état d'énergie potentielle minimale, car dans un système clos la somme des énergies cinétique et potentielle est constante.

Dans l'application de la méthode Boltzmann-inverse aux molécules on applique l'approche précédente aux divers composants des molécules, à partir
d'une même simulation. On produira donc trois ensembles de données : un pour les longueurs des liaisons, un pour les angles de valence et le dernier pour les angles de torsion.

Deux points sont à souligner immédiatement :

1. Les énergies obtenues dépendent de la température T qui devient ainsi un paramètre de simulation.

2. L'approche Boltzmann-inverse ne délivre que des énergies relatives ; il faut fixer la valeur de l'énergie potentielle associée à l'une des classes (par exemple, $C_0$) pour pouvoir en dériver les énergies des autres classes.

Deux questions de fond sont soulevées par une telle démarche. Premièrement, les oscillateurs élémentaires (par exemple, ceux associés aux liaisons) ne sont pas isolés : ils échangent de l'énergie avec les autres oscillateurs de la molécule. Dans quelle mesure ces échanges perturbent-ils la méthode ? Deuxièmement, la relation entre énergie potentielle et densité de probabilité de présence est-elle valide pour tous les cas d'oscillateurs ?

Dans la suite de ce chapitre, on applique la méthode Boltzmann-inverse à l'échelle UA.

## 7.1 Forces intra-moléculaires UA

Le traitement des liaisons, des angles de valence et des angles de torsion est effectué à partir d'une même simulation en dynamique moléculaire, à partir de la molécule $C_6H_{12}$ à la température de 218.152 K. La simulation commence par $2 \times 10^7$ étapes de stabilisation, sans remplissage des classes, puis les classes sont remplies à la fin de chacune des $10^8$ étapes suivantes. Les probabilités de présence dans chaque classe sont calculées à la fin de la simulation.

### Liaisons UA

Le traitement des liaisons UA est effectué en segmentant le domaine de variation des valeurs en 200 classes régulièrement réparties. La Fig. 7.1 montre la courbe obtenue.

Dans le cas des liaisons UA, la méthode Boltzmann-inverse donne un bon résultat pouvant être rendu plus précis en augmentant le nombre de classes et la durée de simulation.
<figure><img src="image_21.png" /><figcaption>FIG. 7.1 – Potentiel de liaison UA obtenu par Boltzmann-inverse à partir de la molécule $C_6H_{12}$ à la température de 218.152 K comparé au potentiel AA CC.</figcaption></figure>

Ainsi, on peut considérer que les échanges d'énergie de la liaison UA avec les autres composants (liaisons, angles de valence ou de torsion) ne perturbent pas l'approche qui associe à la liaison UA un potentiel harmonique très proche de celui du potentiel des liaisons CC à l'échelle AA.

### Angles de valence UA

Pour les angles de valence GGG, on segmente le domaine de variation en 180 classes. La Fig.7.2 montre la courbe obtenue.

Le résultat diffère légèrement du potentiel CCC. On expliquera cette divergence dans la suite par l'effet des angles de valence CCH, fortement corrélés à l'angle de valence considéré.

Les échanges d'énergie de l'angle de valence UA avec les autres composants de la molécule ne perturbent que peu l'approche qui associe à l'angle de valence UA un potentiel harmonique très proche de celui du potentiel des angles CCC à l'échelle AA.
<figure><img src="image_22.png" /><figcaption>FIG. 7.2 – Potentiel de valence UA obtenu par Boltzmann-inverse à partir de la molécule $C_6H_{12}$ à la température de 218.152 K, comparée avec la courbe CCC.</figcaption></figure>

## Angles de torsion UA

Pour les angles de torsion UA, on segmente le domaine de variation en 180 classes. Le potentiel de torsion UA est en bonne correspondance avec la somme des potentiels de torsion AA des angles de torsion partageant une même liaison centrale.

La Fig.7.3 montre une liaison CC (en rouge). Sur la Fig.7.3, on voit qu'il

<figure><img src="image_23.png" /><figcaption>FIG. 7.3 – Angles de torsion partageant la même liaison centrale (en rouge).</figcaption></figure>
y a neuf angles de torsion partageant la même liaison centrale $C_1C_2$ :

* l'angle CCCC $C_0C_1C_2C_3$;

* quatre angles HCCH : $h_2C_1C_2h_4$, $h_3C_1C_2h_4$, $h_2C_1C_2h_5$, $h_3C_1C_2h_5$;

* quatre angle HCCC : $h_2C_1C_2C_3$, $h_3C_1C_2C_3$, $C_0C_1C_2h_4$, $C_0C_1C_2h_5$.

La Fig.7.4 montre la courbe correspondant à l'angle de torsion UA obtenue par Boltzmann-inverse, et la somme des potentiels AA des neuf angles de torsion partageant une même liaison centrale. Les deux courbes sont en bonne correspondance.

<figure><img src="image_24.png" /><figcaption>FIG. 7.4 – Potentiel de torsion UA obtenu par Boltzmann-inverse à partir de la molécule $C_6H_{12}$ à la température de 218.152 K.</figcaption></figure>

Ainsi, les échanges d'énergie de l'angle de torsion UA avec les autres composants de la molécule ne perturbent que peu l'approche qui associe à l'angle de torsion UA la somme des potentiels des angles de torsion qui lui sont associés à l'échelle AA.

## 7.2 Forces inter-moléculaires UA

La réponse de la méthode Boltzmann-inverse aux données issues de la simulation de deux grains UA est donnée par la Fig.7.5.

On constate que la partie droite de la courbe de Boltzmann-inverse n'est pas asymptotique à une parallèle à l'axe $y = 0$ mais, au contraire, tend à
<figure><img src="image_25.png" /><figcaption>FIG. 7.5 – Potentiel inter-moléculaire UA obtenu par Boltzmann-inverse à partir de deux molécules $C_1H_2$.</figcaption></figure>

croître dans les valeurs positives des y ce qui n'est pas conforme à la dynamique. En effet, on s'attend à ce que les forces attractives deviennent de plus en plus faibles au fur et à mesure que les grains s'éloignent.

La courbe obtenue n'est manifestement pas celle d'un potentiel inter-moléculaire.

Pour ce potentiel, on a deux ensembles disjoints d'événements rares : les événements correspondant au domaine fortement répulsif, lorsque les grains sont très proches, d'une part. D'autre part, les événements correspondant au domaine faiblement attractif, lorsque les grains sont très éloignés. Or, dans le premier cas, le potentiel doit être très fort, et dans le second il doit être très faible, alors que la méthode Boltzmann-inverse ne prend en compte que leur rareté relative. La méthode tendra donc à identifier ces événements ce qui est en contradiction avec la dynamique.

La méthode Boltzmann-inverse pour déterminer le potentiel inter-moléculaire UA bute donc sur un obstacle majeur. En anticipant sur la suite, le raisonnement précédent s'applique également à la détermination du potentiel inter-moléculaire CG. La méthode Boltzmann-inverse semble donc en fait inadaptée à la détermination des potentiels inter-moléculaires en général.
# Chapitre 8

## Boltzmann-inverse CG

On applique maintenant la méthode Boltzmann-inverse à l'échelle CG.

La méthode Boltzmann-inverse étant inapte au traitement des forces inter-atomiques (cf. Chap. 7), on ne considérera que les composants intra-moléculaires. Le traitement des liaisons et des angles de valence est effectué à partir d'une même simulation en dynamique moléculaire, à partir de la molécule $C_6H_{12}$ utilisée dans le Chap. 7.

La simulation commence par $2 \times 10^7$ étapes de stabilisation, sans remplissage des classes, puis les classes sont remplies à la fin de chacune des $10^8$ étapes suivantes. Les probabilités de présence dans chaque classe sont calculées à la fin de la simulation.

### 8.1 Liaisons CG

Le traitement des liaisons CG est effectué en segmentant le domaine de variation des valeurs en 200 classes régulièrement réparties. La Fig. 8.1 montre la courbe obtenue pour une température de 218.152 K.

La courbe obtenue dépend de la température de simulation, comme on le voit sur la Fig. 8.2 où trois températures différentes sont considérées. A la température de 111 K, le domaine couvert ne contient pas le point d'inflexion : le dièdre central n'est jamais déclenché. A la température de 218 K, on observe l'inflexion due au déclenchement du dièdre. A la température de 322 K on observe que le potentiel est plus faible au point d'inflexion. A des températures supérieures l'inflexion devient de moins en moins marquée. On est face, ici, à une difficulté : quelle température faut-il choisir pour déter-
<figure><img src="image_26.png" /><figcaption>FIG. 8.1 – Potentiel de liaison CG obtenu par Boltzmann-inverse à partir de la molécule $C_6H_{12}$ à la température de 218.152 K.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>FIG. 8.2 – Potentiel de liaison CG déterminé à trois températures différentes.</figcaption></figure>

miner le potentiel de liaison CG ? La méthode Boltzmann-inverse n’apporte pas de réponse à cette question.
## 8.2 Angles de valence CG

Pour les angles de valence CG, on segmente le domaine de variation en 180 classes. La Fig.8.3 montre la courbe obtenue. La partie droite de la

<figure><img src="image_28.png" /><figcaption>FIG. 8.3 – Potentiel de valence CG obtenu par Boltzmann-inverse à partir de la molécule $C_6H_{12}$ à la température de 218.152 K.</figcaption></figure>

courbe (angle supérieur à 3 radians) est clairement erronée. En effet, l'énergie potentielle est minimale lorsque les grains CG sont alignés, ce qui correspond à un angle de valence CG égal à $\pi$ et n'est pas ce que montre la courbe. Cette anomalie se manifeste indépendemment de la température. De plus, la courbe obtenue dans la cas de l'angle de valence CG varie suivant la température : de même que lors de la détermination du potentiel de liaison CG, il n'est pas clair de savoir quelle température choisir pour déterminer le potentiel de valence CG.

## Conclusion

On peut tirer plusieurs conclusions concernant la méthode Boltzmann-inverse :

* Le traitement correct des événements rares nécessite des grands nombres d'événements, donc des simulations pouvant être très longues.
* Lorsque le nombre d'événements n'est pas suffisant pour traiter correctement un intervalle de valeurs, il se peut que l'énergie potentielle soit sous-estimée.

* Les courbes obtenues peuvent dépendre de la température ce qui pose alors la question du choix de la température à considérer.

* Par ailleurs, la température de simulation limite l'intervalle des valeurs analysées et donc limite le domaine de définition du potentiel obtenu.

* Enfin, la méthode donne de mauvaises réponses dans le cas des potentiels inter-moléculaires.

On va maintenant considérer un autre moyen pour déterminer les potentiels : la méthode par minimisation.
# Chapitre 9

# Méthode de minimisation

On propose dans ce chapitre une méthode pour déterminer les potentiels de composants moléculaires fondée sur une technique de minimisation sous-conditions. Cette méthode repose fondamentalement sur l'existence de liens géométriques avec l'échelle de base AA.

Pour déterminer le potentiel d'un composant moléculaire (par exemple, angle de valence UA), on choisit une molécule dans laquelle ce composant apparaît et on fixe la valeur *v* du composant. On calcule alors le potentiel global *p* de la molécule sous-jacente (aspect géometrique) à l'échelle AA, après minimisation préservant la valeur *v* du composant (aspect sous-conditions). Le potentiel du composant avec la valeur *v* est, par définition, *p*. En faisant varier *v*, on détermine ainsi le potentiel recherché.

Par exemple, pour l'angle de valence UA, on fixe l'angle entre trois grains UA c'est-à-dire entre trois atomes de carbone et on minimise tous les composants AA de la molécule autres que cet angle AA. Ainsi, les atomes de carbone des trois grains de l'angle UA restent immobiles tandis que les autres atomes, de carbone ou d'hydrogène, peuvent se déplacer pour minimiser l'énergie de la molécule globale. Les liaisons et angles impliquant les atomes d'hydrogène sont ainsi placés à des positions d'équilibre dont l'énergie est minimale. Il en est de même pour les atomes de carbone, à l'exception de ceux formant l'angle de valence UA. L'énergie finale de la molécule AA après minimisation complète est celle que l'on associe à l'angle de valence UA.

Dans la suite de ce chapitre, on applique la méthode de minimisation à l'échelle UA. On rappelle que les grains UA sont formés d'un atome de carbone et des atomes d'hydrogène qui lui sont liés. Le centre d'un grain UA coincide avec l'atome de carbone. La minimisation appliquée à CG est
décrite dans le Chap. 10.

## 9.1 Liaisons UA

Pour traiter le potentiel UA de liaison $G_2G_2$, on considère la molécule équivalente (en fait, un fragment) $C_2H_4$ et les minimisations sont effectuées en maintenant constante la distance entre les deux carbones de la molécule. La concordance avec le potentiel OPLS AA est parfaite et illustrée sur la Fig. 9.1. Ce résultat reflète l'indépendance qui existe entre la liaison CC,

<figure><img src="image_29.png" /><figcaption>FIG. 9.1 – Potentiel de la liaison $G_2G_2$ obtenu par minimisation, comparé au potentiel OPLS CC.</figcaption></figure>

d'une part, et les angles de valence HCH et HCC, d'autre part, dans la molécule $C_2H_4$.

## 9.2 Angles de valence UA

Pour traiter le potentiel de valence $G_2G_2G_2$, on considère le fragment $C_3H_6$ et les minimisations sont effectuées en maintenant constant l'angle de valence CCC. Les résultats sont exactement les mêmes à partir de la molécule $C_3H_8$.
Les résultats des minimisations sont montrés sur la Fig.9.2. On constate que la courbe produite par la minimisation ne coincide pas totalement avec le potentiel CCC, en particulier dans la partie correspondant à l'extension de l'angle de valence.

<figure><img src="image_30.png" /><figcaption>FIG. 9.2 – Potentiel de l'angle de valence $G_2G_2G_2$ obtenu par minimisation, comparé au potentiel OPLS CCC.</figcaption></figure>

Le potentiel UA de l'angle de valence $G_2G_2G_2$ diffère donc du potentiel AA de l'angle de valence CCC.

En enlevant l'énergie potentielle des huit angles de valence CCH, on retrouve le potentiel CCC, comme on le constate sur la Fig.9.3. La différence avec le potentiel de valence CCC tient à la corrélation partielle (en fait, uniquement en extension de l'angle de valence) existant entre l'angle CCC et les angles CCH.

Le potentiel $G_2G_2G_2$ est en fait la somme de deux fonctions, la première $f_{inf}$ définie pour les angles inférieurs à l'angle d'équilibre, la seconde $f_{sup}$ définie pour les angles supérieurs. La fonction $f_{inf}$ coincide sur son espace de définition avec la fonction harmonique CCC, tandis que $f_{sup}$ multiplie CCC par un facteur 1.142.

Les deux fonctions $f_{sup}$ et $f_{inf}$ sont toutes les deux approximées par une unique fonction harmonique qui est la fonction harmonique CCC multipliée par un facteur 1.1. C'est cette approximation que l'on utilisera par la suite;
<figure><img src="image_31.png" /><figcaption>FIG. 9.3 – Potentiel de l'angle de valence $G_2G_2G_2$ obtenu par minimisation, en enlevant les composantes dues aux angles de valence CCH, comparé au potentiel OPLS CCC.</figcaption></figure>

ainsi le potentiel UA de valence $G_2G_2G_2$ sera un potentiel harmonique égal au potentiel de valence AAA multiplié par le facteur 1.1.

## 9.3 Angles de torsion UA

Pour traiter le potentiel de torsion $G_2G_2G_2G_2$, on considère le fragment $C_4H_8$ et les minimisations sont effectuées en maintenant constant l'angle de torsion CCCC. Les résultats sont exactement les mêmes à partir de la molécule $C_4H_{10}$.

L'angle de torsion CCCC est en fait totalement corrélé aux angles de torsion HCCH et HCCC centrés sur le milieu de la liaison CC centrale. Quatre angles de torsion HCCH et quatre angles de torsion HCCC sont ainsi concernés. La Fig.9.4 montre le potentiel obtenu par minimisation comparé à la somme CCCC+4HCCH+4HCCC. On constate la correspondance des deux courbes.
<figure><img src="image_32.png" /><figcaption>FIG. 9.4 – Potentiel de l'angle de torsion $G_2G_2G_2G_2$ obtenu par minimisation, comparé à la somme des potentiels OPLS correspondant à $CCCC+4HCCH+4HCCC$.</figcaption></figure>

## 9.4 Forces inter-moléculaire UA

Pour déterminer le potentiel Lennard-Jones entre deux grains UA, on utilise deux molécules $C_1H_2$ placées à distance. La courbe obtenue en faisant varier cette distance est montrée sur la Fig.9.5. On obtient un potentiel Lennard-Jones 6-12 qui prend en compte la présence des hydrogènes, comme on le constate sur la Fig.9.6. Dans la Fig.9.6, les paramètres du potentiel Lennard-Jones 6-12 sont $\epsilon = 5.5 \times \epsilon_{CC}$ et $\sigma = 0.917 \times \sigma_{CC}$.

## Conclusion

La molécule $C_2H_4$ est utilisée pour le calcul de la liaison UA et on obtient exactement le même potentiel que celui de la liaison AA CC (Fig.9.1). En d'autres termes, les potentiels AA et UA des liaisons coincident et sont tous deux harmoniques.

La molécule $C_3H_6$ est utilisée pour la détermination du potentiel d'angle de valence UA (Fig.9.2). Le potentiel de valence UA est différent (et non-harmonique) mais cependant très proche du potentiel de valence AA.

A partir de la molécule $C_4H_8$ on obtient la courbe de la Fig.9.4 comme
<figure><img src="image_33.png" /><figcaption>FIG. 9.5 – Potentiel inter-moléculaire entre deux grains UA obtenu par minimisation, comparé au potentiel OPLS entre deux atomes de carbone.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>FIG. 9.6 – Potentiel inter-moléculaire entre deux grains UA obtenu par minimisation, comparé à une courbe Lennard-Jones 6-12.</figcaption></figure>

potentiel de l'angle de torsion UA. Le potentiel UA d'un angle de torsion $\theta$ est donc la somme des potentiels CCCC de l'angle $\theta$, plus quatre fois le potentiel HCCH de $\theta$, plus quatre fois le potentiel HCCC de $\theta$.
Enfin, la Fig.9.5 montre le potentiel inter-moléculaire UA entre deux
grains $G_2$. Ce potentiel a trois caractéristiques principales : il a la forme d'un
potentiel Lennard-Jones 6-12 ; il coincide avec le potentiel inter-moléculaire
AA lorsque les grains sont très proches ou suffisamment éloignés ; il est net-
tement plus “creusé” que le potentiel inter-moléculaire AA, dans sa partie
centrale.

La préservation de la forme Lennard-Jones 6-12 indique que le potentiel
inter-moléculaire UA se comporte globalement comme un potentiel inter-
moléculaire AA, c'est-à-dire que la présence des atomes d'hydrogène ne change
pas la forme des interactions entre grains UA.

L'aspect “creusé” est lié à la présence des atomes d'hydrogène qui à faible
distance contribuent de manière significative aux forces attractives entre
grains.
# Chapitre 10

## Minimisation CG

Comme on l'a vu en 6.2, la détermination du potentiel CG par minimisation se décline sous trois aspects uniquement : potentiel de liaison, potentiel de valence et potentiel inter-moléculaire; il n'y a pas, au niveau CG, de potentiel de dièdre.

### 10.1 Liaisons CG

Pour déterminer le potentiel de liaison CG, on utilise la molécule $C_4H_8$ de la Fig. 10.1.

<figure><img src="image_35.png" /><figcaption>FIG. 10.1 – Molécule $C_4H_8$. La liaison CG est représentée en rouge.</figcaption></figure>

On considère cette molécule comme formée de deux grains : le premier englobe les deux carbones $C_0$ et $C_1$ ainsi que les hydrogènes qui leurs sont
attachés, et le second les deux carbones $C_2$ et $C_3$ avec leurs hydrogènes attachés. Les deux grains sont centrés sur le milieu $m_0$ de la liaison AA entre $C_0$ et $C_1$, d'une part, et le milieu $m_1$ de la liaison AA entre $C_2$ et $C_3$, d'autre part.

Le potentiel de liaison CG est calculé en minimisant l'énergie de la molécule tout en maintenant constante la longueur $l$ de la liaison CG entre $m_0$ et $m_1$ (minimisation sous-condition).

On commence tout d'abord par définir l'énergie $eg_0$ du grain $G_0$ centré sur $m_0$ comme étant la somme des énergies des composants internes du grain :

* liaisons (5) : $h_0C_0, h_1C_0, C_0C_1, C_1h_2, C_1h_3$;

* angles de valence (6) : $h_0C_0C_1, h_0C_0h_1, h_1C_0C_1, h_2C_1C_0, h_2C_1h_3, h_3C_1C_0$;

* dièdres (4) : $h_0C_0C_1h_2, h_0C_0C_1h_3, h_1C_0C_1h_2, h_1C_0C_1h_3$.

On définit de la même manière l'énergie $eg_1$ du grain $G_1$ centré sur $m_1$.

L'énergie de liaison $el$ entre les deux grains est la somme des énergies des composants reliant les deux grains :

* liaisons (1) : $C_1C_2$;

* angles de valence (6) : $C_0C_1C_2, h_2C_1C_2, h_3C_1C_2, C_1C_2C_3, C_1C_2, h_4, C_1C_2h_5$;

* dièdres (13) : $h_0C_0C_1C_2, h_1C_0C_1C_2, C_0C_1C_2C_3, C_0C_1C_2h_4, C_0C_1C_2h_5, h_2C_1C_2h_4, h_2C_1C_2h_5, h_2C_1C_2C_3, h_3C_1C_2h_4, h_3C_1C_2h_5, h_3C_1C_2C_3, C_1-C_2C_3h_6, C_1C_2C_3h_7$.

L'énergie (potentielle) $ep$ de la molécule $C_4H_8$ est la somme des énergies de ses 50 composants AA (liaisons, angles de valence et dièdres). Cette énergie doit être égale à la somme de $eg_0$, $eg_1$ et $el$ :

$$ep = eg_0 + eg_1 + el \qquad (10.1)$$

## 10.2 Angles de valence CG

Pour déterminer le potentiel de valence CG, on utilise la molécule $C_6H_{12}$ composée de trois grains (Fig.10.2). Le premier grain englobe les carbones $C_0$ et $C_1$, le second les carbones $C_2$ et $C_3$, et le troisième les carbones $C_4$ et $C_5$. Les grains sont centrés sur le milieu $m_0$ de la liaison AA $C_0C_1$, le milieu $m_1$ de la liaison AA $C_2C_3$, et $m_2$ milieu de la liaison AA $C_4C_5$.

Le potentiel de valence CG est calculé en minimisant l'énergie de la molécule tout en maintenant constant l'angle formé par les trois centres $m_0, m_1$ et $m_2$ (il s'agit donc également d'une minimisation sous-condition).
<figure><img src="image_36.png" /><figcaption>FIG. 10.2 – Molécule $C_6H_{12}$. L'angle de valence CG est l'angle formé par les deux flèches rouges.</figcaption></figure>

L'énergie d'angle *ea* associée aux trois grains se résume à celle du dièdre central $C_1C_2C_3C_4$.

L'énergie (potentielle) *ep* de la molécule $C_6H_{12}$ est la somme des énergies de ses 86 composants AA (liaisons, angles de valence et dièdres). Cette énergie doit être égale à la somme de $eg_0$, $eg_1$, $eg_2$, plus les deux énergies $el_1$ et $el_2$ de liaison entre les grains, plus *ea* :

$$ep = eg_0 + eg_1 + eg_2 + el_1 + el_2 + ea \quad (10.2)$$

## 10.3 Pondération

L'énergie potentielle d'une molécule CG est la somme des énergies des grains, des énergies des liaisons CG et de celles des angles de valence CG la composant. Cette somme doit être égale à l'énergie potentielle de la molécule AA. Pour obtenir ce résultat, il est nécessaire de *pondérer* les énergies des grains, et éviter ainsi de compter plusieurs fois ces énergies. Le facteur de pondération *P* choisi est $1/5$ :

$$P = 1/5 \quad (10.3)$$

Pour justifier la valeur de *P*, considérons la molécule CG sur la Fig.10.3, plus précisemment le grain $G_2$. Celui-ci participe des deux liaisons $G_1G_2$ et $G_2G_3$ ainsi que des trois angles de valence $G_1G_2G_3$, $G_0G_1G_2$ et $G_2G_3G_4$. Ce sont
les cinq seuls composants CG dans lesquels $G_2$ apparaît. Il en va de même pour chacun des grains CG, à l'exception de ceux proches des deux extrémités de la molécule que l'on devra traiter différemment, c'est-à-dire avec des pondérations différentes. Pour simplifier, dans la suite, on ne considérera

<figure><img src="image_37.png" /><figcaption>FIG. 10.3 – Fragment de molécule CG.</figcaption></figure>

que des grains éloignés des extrémités (donc, des grains semblables à $G_2$, sur la Fig.10.3) et la seule pondération à utiliser sera $P = 1/5$<fnref n="1" />.

Dans l'implémentation, les minimisations lors des calculs des potentiels CG sont effectuées avec des perturbations allant de $10^{-3}$ nm à $10^{-9}$ nm.

## 10.4 Potentiel des liaisons CG

Le potentiel d'une liaison CG est la somme des énergies pondérées des deux grains liés et de l'énergie de liaison entre eux :

$$eb = P \times eg_1 + P \times eg_2 + el \quad (10.4)$$

On obtient pour le potentiel de liaison CG la courbe de la Fig.10.4.

La courbe obtenue ne dépend pas de la taille de la molécule. Elle montre un point d'inflexion aux environs de 0.24 nm. Cette inflexion résulte du déclenchement du dièdre central ($C_0C_1C_2C_3$).

Dans le calcul du potentiel de la liaison CG, la taille de la liaison CG est décrémentée jusqu'à 0.18nm, puis augmentée jusqu'à 0.28 nm, par pas de $10^{-4}$ nm.

<footnote n="1">¹ Pour un traitement complet des alkanes, il faudrait distinguer trois sortes de grains : les grains de type E aux extrémités ; les grains de type F connectés directement à ceux de type E ; et les autres, de type G. On aurait à définir les pondérations des potentiels de liaisons EF, FG, et celles des potentiels de valence EFG et FGG. La pondération 1/5 est celle des liaisons GG et des angles de valence GGG.</footnote>
<figure><img src="image_38.png" /><figcaption>FIG. 10.4 – Potentiel de liaison CG entre deux grains $C_2H_4$. Le potentiel est obtenu par minimisation à partir de la molécule $C_4H_8$.</figcaption></figure>

## 10.5 Potentiel des angles de valence CG

Le potentiel de valence CG est la somme des énergies pondérées des trois grains composant l'angle de valence, plus l'énergie de l'angle :

$$ev = P \times eg_1 + P \times eg_2 + P \times eg_3 + ea \quad (10.5)$$

Le potentiel de valence CG correspond à la courbe de la Fig.10.5.

Comme pour le potentiel de liaison, la courbe obtenue ne dépend pas de la taille de la molécule. La courbe obtenue montre trois points d'inflexion. Les deux points d'inflexion à gauche correspondent aux déclenchements des deux dièdres centraux des deux liaisons. Le point d'inflexion de droite correspond au déclenchement du dièdre central $C_1C_2C_3C_4$.

Dans l'implémentation, l'angle de valence CG est décrémenté de $\pi$ radians jusqu'à 0.9 radians, puis augmenté jusqu'à $\pi$, par pas de 0.017 radians (on obtient 262 mesures).

## 10.6 Forces inter-moléculaires CG

Pour déterminer la force s'exerçant entre deux grains CG, on considère le système formé de deux molécules $C_2H_4$ placées à une distance fixe $d$. Le po-
<figure><img src="image_39.png" /><figcaption>FIG. 10.5 – Potentiel de valence CG entre trois grains $C_2H_4$. Le potentiel est obtenu par minimisation à partir de la molécule $C_6H_{12}$. La valeur minimale de l'énergie est nulle et correspond à un angle de valence $\pi$.</figcaption></figure>

tentiel est calculé en minimisant l'énergie potentielle *ep* des forces de van Der Waals s'exerçant entre les deux grains, tout en maintenant constante la distance *d* (minimisation sous-condition). L'énergie *ep* est la valeur du potentiel entre les deux grains obtenue après minimisation, en maintenant la distance *d*. En faisant varier le paramètre *d*, on détermine la courbe de potentiel de van Der Walls à l'échelle CG montrée sur la Fig.10.6. L'énergie *ep* est la somme des trente-six énergies des potentiels de Lennard-Jones associés aux paires d'atomes entre les deux grains (les potentiels Lennard-Jones à l'échelle AA sont représentés sur la Fig.2.6).

Le potentiel obtenu par minimisation de la Fig. 10.6 n'a pas la forme d'un potentiel Lennard-Jones 6-12, défini par l'équation Eq.2.4. La Fig.10.7 compare le potentiel obtenu avec un potentiel Lennard-Jones 6-12 approchant ($\epsilon = 12.5 \times \epsilon_{CC}$, $\sigma = 0.9 \times \sigma_{CC}$) : on voit que le potentiel obtenu diffère sur un domaine significatif (entre 0.4 *nm* et 0.7 *nm*) d'un potentiel 6-12.

Le potentiel inter-moléculaire CG est linéaire dans le domaine où il ne coincide pas avec un potentiel Lennard-Jones 6-12. Dans ce domaine, la minimisation a pour effet de réorganiser les composants AA “inclus” dans les deux grains CG de façon à minimiser leur énergies. Au sein de chaque grain, les échanges d'énergie entre ces composants induisent une évolution linéaire
<figure><img src="image_40.png" /></figure>

FIG. 10.6 – Potentiel inter-moléculaire entre deux grains CG obtenu par minimisation.

<figure><img src="image_41.png" /></figure>

FIG. 10.7 – Potentiel de van Der Waals entre deux grains obtenu par minimisation à l'échelle CG, comparé avec un potentiel Lennard-Jones 6-12 approchant.

de l'énergie globale du grain. En dehors de la zone linéaire, le potentiel inter-moléculaire CG se comporte comme un potentiel Lennard-Jones 6-12 standard ce qui signifie l'absence de réorganisation des grains CG dans les zones
de forte proximité ou de fort éloignement des grains.

## 10.7 Comparaison avec l'approche Boltzmann-inverse

Pour comparer les approches par Boltzmann-inverse et par minimisation dans le calcul du potentiel de liaison CG, on superpose les courbes des Fig.10.4 et Fig.8.1 dans la Fig.10.8.

<figure><img src="image_42.png" /><figcaption>FIG. 10.8 – Potentiel de liaison CG obtenu par Boltzmann-inverse à partir de la molécule $C_6H_{12}$ à la température de 218.152 K comparé avec le potentiel obtenu par minimisation.</figcaption></figure>

On constate une bonne concordance entre les courbes, sauf en ce qui concerne le point d'inflexion, aux alentours de 0.24 *nm*. Ce point correspond au déclenchement du dièdre central, un événement rare auquel, donc, se trouve associée une énergie potentielle plus faible que celle qu'elle devrait être (ce qui peut être considéré comme une faiblesse intrinsèque de l'approche Boltzmann-inverse).

Pour comparer les approches par Boltzmann-inverse et par minimisation dans le calcul du potentiel de valence CG, on superpose les courbes des Fig.10.5 et Fig.8.3 dans la Fig.10.9.
10.7. COMPARAISON AVEC L'APPROCHE BOLTZMANN-INVERSE77

<figure><img src="image_43.png" /><figcaption>FIG. 10.9 – Potentiel de valence CG obtenu par Boltzmann-inverse à partir de la molécule $C_6H_{12}$ à la température de 218.152 K (en bleu) comparé avec le potentiel obtenu par minimisation (en rouge).</figcaption></figure>

La courbe Boltzmann-inverse est systématiquement plus haute que la courbe de minimisation pour les angles inférieurs à 1.9 radians. Cela est dû à la température trop faible pour que l'exploration angulaire au dessous de 1.9 soit réellement significative.

On peut considérer que la partie de la courbe Boltzmann-inverse sous la courbe de minimisation résulte d'un manque de précision dû à un ensemble d'événements trop réduit.

Pour augmenter le nombre d'événements pris en compte, on peut utiliser des molécules plus longues, par exemple la molécule $C_{30}H_{60}$ à la température de 206 K. Le résultat est montré sur la Fig.10.10 où l'on constate que la partie de la courbe Boltzmann-inverse sous la courbe de minimisation est alors fortement réduite.

Avec un nombre d'événements suffisants, on conjecture que la courbe Boltzmann-inverse sera toujours au dessus de la courbe obtenue par minimisation. On pourrait voir cela comme une validation (partielle...) de la minimalité de la courbe obtenue par minimisation, inférieure à toutes les celles obtenues à partir de simulations effectives.
<figure><img src="image_44.png" /><figcaption>FIG. 10.10 – Potentiel de valence CG obtenu par Boltzmann-inverse à partir de la molécule $C_6H_{12}$ à la température de 218.152 K (en bleu) comparé avec celui obtenu à partir de la molécule $C_{30}H_{60}$ à la température de 206.619 K (en vert).</figcaption></figure>
# Chapitre 11

## Simulations en UA

À l'aide de la méthode de minimisation, on a déterminé les potentiels UA des composants des molécules formées de grains $G_2$, ce qui permet de simuler les molécules à cette échelle. Ces potentiels UA obtenus à partir des potentiels AA sont résumés dans la Fig.11.1.

$$ 
\begin{aligned}
k_{G_2G_2} &= k_{CC} \qquad r_{G_2G_2} = r_{CC} \\
k_{G_2G_2G_2} &= 1.1 \times k_{CCC} \qquad \theta_{G_2G_2G_2} = \theta_{CCC} \\
A1_{G_2G_2G_2G_2} &= A1_{CCCC} \qquad A2_{G_2G_2G_2G_2} = A2_{CCCC} \\
A3_{G_2G_2G_2G_2} &= A3_{CCCC} + 4 \times A3_{CCCH} + 4 \times A3_{HCCH} \\
\epsilon_{G_2G_2} &= 5.5 \times \epsilon_{CC} \qquad \sigma_{G_2G_2} = 0.917 \times \sigma_{CC}
\end{aligned}
$$

FIG. 11.1 – Valeurs des paramètres UA associés aux grains $G_2$, exprimées en unités internes (les paramètres OPLS AA sont donnés dans la Fig.5.1).

Le potentiel de l'angle de valence $G_2G_2G_2$ est approximé comme un potentiel harmonique dont la raideur est celle de CCC multiplié par 1.1, comme on l'a vu en 9.2. Le potentiel de torsion a été déterminé en 9.3 et le potentiel intermoléculaire en 9.4.

Dans le système de DM, les grains $G_2$ (c'est-à-dire formés d'un atome de carbone et de deux atomes d'hydrogène) apparaissent comme les boules de couleur cyan. La Fig.11.2 montre une molécule UA composée de 8 grains $G_2$ à l'équilibre (son énergie est nulle).
<figure><img src="image_45.png" /><figcaption>FIG. 11.2 – Molécule UA composée de 8 grains $G_2$, à l'équilibre.</figcaption></figure>

Les simulations des molécules UA sont très stables, comme on le voit sur la Fig.11.3. Sur la partie gauche, une molécule à huit grains $G_2$ est montrée en cours de simulation. Initialement on donne à la molécule une énergie initiale en décalant les grains vers les y positifs (axe vertical). Sur la partie droite, on a l'évolution des énergies au cours du temps. La courbe rouge représente l'énergie totale de la molécule, dont on voit qu'elle reste stable.

<figure><img src="image_46.png" /><figcaption>FIG. 11.3 – Simulation d'une molécule UA composée de 8 grains $G_2$.</figcaption></figure>
La Fig.11.4 décrit la simulation de deux molécules UA à 8 grains, placées en vis-à-vis à une distance de 1.8 nm. Initialement, les deux molécules sont à l'équilibre. L'énergie des forces inter-moléculaires est en jaune, l'énergie cinétique en bleu, l'énergie intra-moléculaire est en vert et la somme des énergies est en rouge.

<figure><img src="image_47.png" /><figcaption>FIG. 11.4 – Simulation de deux molécules UA composées de 8 grains $G_2$.</figcaption></figure>

Les simulations à l'échelle UA sont plus efficaces qu'à l'échelle AA : le nombre de degrés de liberté à simuler est moindre et le pas-de-temps peut être augmenté. On va montrer cette efficacité accrue sur un exemple.

On simule tout d'abord une molécule AA $C_8H_{16}$ pendant $10^7$ instants avec un pas-de-temps de 0.1 femto-seconde. Le temps réel de simulation<fnref n="1" /> est de 1053 secondes. Cette simulation correspond à un temps physique de 0.5 nano-secondes. Les énergies sont montrées sur la partie gauche de la Fig.11.5.

On simule maintenant la molécule UA correspondante, composée de 8 grains $G_2$ avec un pas-de-temps de 0.5 femto-secondes. La simulation correspond donc à un temps physique de 2.5 nano-secondes. Le temps réel de simulation est de 175 secondes. Les énergies sont montrées sur la partie droite de la Fig.11.5.

Par rapport à la simulation précédente en AA, le temps-réel de simulation, est donc réduit d'un facteur 6 alors que le temps physique est, lui, augmenté d'un facteur 5. Le gain de UA par rapport à AA est donc bien réel.

<footnote n="1">¹sur un MacBook Pro, processeur i7 de 2.6GHz, avec 16GB de mémoire.</footnote>
<figure><img src="image_48.png" /><figcaption>FIG. 11.5 – Energies durant 10 millions d'instants. À gauche, molécule AA $C_8H_{16}$. À droite, molécule UA avec 8 grains.</figcaption></figure>
# Chapitre 12

## Simulations en CG

Contrairement aux potentiels AA et UA, le potentiel CG n'a pas une forme standard, exprimée en terme de fonctions harmoniques, de fonctions "3 cosinus" et de fonctions Lennard-Jones 6-12. On va donc commencer par déterminer les constituants du potentiel CG.

### 12.1 Détermination du potentiel CG

#### Liaisons

Le potentiel de liaison CG apparaît comme étant constitué de deux fonctions harmoniques, s'appliquant de part-et-d'autre d'une valeur ($split = 0.241 nm$). La Fig.12.1 compare le potentiel CG obtenu par minimisation (à partir de la molécule $C_4H_8$, en rouge) et celui construit "par cas" à l'aide de deux fonctions harmoniques (en bleu).

La fonction harmonique est définie par :

$$harm: k, r, d \rightarrow k \times (r - d)^2$$

où $k$ est la raideur, $r$ est la valeur d'équilibre et $d$ est un paramètre (une longueur en $nm$ dans le cas des liaisons et un angle en radians dans le cas des angles de valence).

Le potentiel CG de liaison est défini par la fonction cgBond suivante :

$$split = 0.241 \\ h_1: d \rightarrow harm(16.5, 0.225, d)$$
<figure><img src="image_49.png" /><figcaption>FIG. 12.1 – Le potentiel CG (en bleu) est défini “par cas” à l’aide de deux fonctions harmoniques.</figcaption></figure>

$$h_2 : d \rightarrow \text{harm}(38.5, 0.255, d)$$
$$\text{cgBond} : d \rightarrow \text{if } d < \text{split then } h_1(d) + 0.0049 \text{ else } h_2(d)$$

À gauche de 0.241 (c’est-à-dire pour les valeurs inférieures à 0.241 nm), s’applique la fonction harmonique $h_1$ définie par $k = 16.5$ et $r = 0.225$. Un décalage en $y$ constant de 0.0049 est ajouté au résultat.

Pour les valeurs supérieures à 0.241 nm, s’applique la fonction harmonique $h_2$ avec les paramètres $k = 38.5$ et $r = 0.255$.

## Angles de valence

La determination du potentiel de valence CG utilise quatre fonctions harmoniques $h_0, h_1, h_2, h_3$ de la manière suivante :

$$\begin{align*}
\text{split}_0 &= 2.56 & h_0 : a \rightarrow \text{harm}(0.0095, \pi, a) \\
\text{split}_1 &= 1.79 & h_1 : a \rightarrow \text{harm}(0.0055, \text{split}_0, a) + 0.004 \\
\text{split}_2 &= 1.6 & h_2 : a \rightarrow \text{harm}(0.041, \text{split}_1, a) + 0.0056 \\
& & h_3 : a \rightarrow \text{harm}(0.0059, 2.0, a) + 0.0047
\end{align*}$$
Le potentiel de valence CG est alors défini “par cas” par :

$$ \begin{array}{l} \text{cgValence} : a \rightarrow \text{if } a > \text{split}_0 \text{ then } h_0(a) \\ \text{else if } a > \text{split}_1 \text{ then } h_1(a) \\ \text{else if } a > \text{split}_2 \text{ then } h_2(a) \\ \text{else } h_3(a) \end{array} $$

Le potentiel de valence CG est représenté par la courbe bleue sur la Fig.12.2, la courbe rouge étant le résultat obtenu par minimisation.

<figure><img src="image_50.png" /><figcaption>FIG. 12.2 – Le potentiel de valence CG (en bleu) est défini “par cas” à l’aide de quatre fonctions harmoniques.</figcaption></figure>

## Forces de van Der Waals

Les forces inter-moléculaires sont généralement représentées par des fonctions Lennard-Jones de la forme :

$$ lj6\_12 : \epsilon, \sigma, r \rightarrow 4\epsilon \times ((\sigma/r)^{12} - (\sigma/r)^{6}) $$

Pour le potentiel CG, on commence par définir une fonction Lennard-Jones 6-12 $lj$ particulière, définie par :

$$ \begin{array}{l} \epsilon = 0.00336 \\ \sigma = 0.318 \\ \text{lj} : d \rightarrow \text{lj6\_12}(\epsilon, \sigma, d) \end{array} $$
La fonction $lj$ est représentée par la courbe verte de la Fig.12.3. Pour déterminer complètement le potentiel CG, on introduit deux points : le point $p_1 = (0.365, -\epsilon)$ dont la coordonnée en $y$ est minimale et le point $p_2 = (0.68, -000195)$. Entre ces deux points, le potentiel est globalement une droite, correspondant à la fonction affine passant par les deux points. La définition du potentiel inter-moléculaire CG est alors :

$$cgLj : d \rightarrow \textbf{if } d < 0.365 \textbf{ or } d > 0.68 \textbf{ then } lj(d) \textbf{ else } line(d)$$

La fonction *line* étant la fonction affine passant par les deux points $p_1$ et $p_2$. Le potentiel CG est représenté par la courbe bleue, sur la Fig.12.3, la courbe rouge représentant le potentiel CG obtenu par minimisation.

<figure><img src="image_51.png" /><figcaption>FIG. 12.3 – Le potentiel inter-moléculaire CG (en bleu) est défini “par cas” à partir d'une fonction Lennard-Jones 6-12 (en vert) et d'une fonction affine. En rouge, le potentiel obtenu par minimisation.</figcaption></figure>

## 12.2 Exécution CG

Pour illustrer l'exécution CG, on considère tout d'abord une molécule CG isolée composée de six grains $CG_4$. Les coordonnées en $x$ sont représentées horizontalement, celles en $y$ verticalement et celles en $z$ en hauteur.
Une énergie est fournie initialement à la molécule à l'équilibre en décalant son premier grain dans la direction des y croissants. Ainsi, la première liaison CG est étirée et c'est l'énergie due à son étirement qui est fournie initialement à la molécule. La partie droite de la Fig.12.4 montre l'évolution des diverses composantes de l'énergie au cours du temps. Les mesures sont prises tous les $10^5$ instants. L'énergie cinétique est en bleu, l'énergie intra-moléculaire en vert et l'énergie totale en rouge (l'énergie inter-moléculaire est toujours, ici, nulle). La molécule en cours d'exécution est montrée sur la partie gauche de la figure et on constate qu'elle reste linéaire. Cela est cohérent avec le fait que les forces, ici, sont dues uniquement aux liaisons.

<figure><img src="image_52.png" /><figcaption>FIG. 12.4 – Simulation d'une molécule CG composée de 6 grains $CG_4$. Initialement, seule la première liaison a de l'énergie. A gauche, la molécule après $3 \times 10^8$ instants.</figcaption></figure>

On change maintenant les conditions initiales en faisant tourner le premier grain autour du deuxième, tout en restant dans le même plan. L'énergie initiale est alors celle du premier angle de valence, formé des trois premiers grains. Les résultats d'exécution sont montrés sur la Fig.12.5 sur laquelle on constate que la molécule n'est plus linéaire mais reste toujours dans un même plan : en effet, toutes les forces qui s'exercent au cours de l'exécution sont co-planaires.

Finalement, on introduit une torsion dans la molécule, en décalant le premier grain dans la coordonnée x et le dernier grain dans la coordonnée z. Ces deux décalages font que la molécule n'est plus planaire. La simulation est montrée sur la Fig.12.6.
<figure><img src="image_53.png" /><figcaption>FIG. 12.5 – Simulation d’une molécule CG composée de 6 grains $CG_4$. Initialement, seul le premier angle de valence a de l’énergie. A gauche, la molécule après $3 \times 10^8$ instants.</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>FIG. 12.6 – Simulation d’une molécule CG composée de 6 grains $CG_4$. Initialement, le premier grain est décalé en $x$ et le dernier en $z$. La molécule n’est plus planaire. A gauche, la molécule après $3 \times 10^8$ instants.</figcaption></figure>

On va maintenant considérer le cas des forces inter-moléculaires. On met en face-à-face deux molécules identiques et à l’équilibre. Les deux molécules sont donc linéaires et co-planaires. Chaque molécule est composée de six grains.

La partie droite de la Fig.12.7 montre les diverses énergies : en bleu,
l'énergie cinétique; en vert, l'énergie intra-moléculaire; en jaune, l'énergie inter-moléculaire et en rouge l'énergie totale. La simulation dure $10^8$ instants, correspondant à 25 ns (le pas-de-temps est de $10^{-4}$ ps).

<figure><img src="image_55.png" /><figcaption>FIG. 12.7 – Simulation de deux molécules CG composées de 6 grains $CG_4$ et initialement co-planaires, à gauche. Variation des énergies à droite.</figcaption></figure>

Les deux molécules, initialement dans un même plan (co-planaires), le restent au cours de l'exécution, comme le montre la partie gauche de la figure. Les forces de van Der Waals restent dans le même plan que les molécules.

Ce n'est pas le cas lorsque les deux molécules ne sont pas co-planaires. Pour illustrer cela, on reprend la simulation précédente mais on fait pivoter d'un angle de un degré une des deux molécules autour de l'axe des $x$. Ainsi les deux molécules ne sont plus initialement co-planaires. Rapidement les molécules se mettent à décrire des trajectoires 3D.

Les molécules CG plus longues que celles représentées sur les images précédentes manifestent des repliements pouvant être très accentués. C'est le cas par exemple pour la molécule CG de 40 grains de la Fig. 12.9.
<figure><img src="image_56.png" /><figcaption>FIG. 12.8 – Simulation de deux molécules CG composées de 6 grains $CG_4$, initialement non co-planaires, après $20^6$ instants.</figcaption></figure>

<figure><img src="image_57.png" /><figcaption>FIG. 12.9 – Molécule CG composée de 40 grains $CG_4$, initialement non co-planaires, en cours de simulation.</figcaption></figure>
# Chapitre 13

## Reconstructions

Les reconstructions de molécules permettent de changer, au cours d'une même simulation, l'échelle à laquelle une molécule est simulée. Typiquement, lorsqu'une molécule est isolée, on préfèrera la simuler à l'échelle UA ou à l'échelle CG qui sont beaucoup plus efficaces que l'échelle AA. Par contre, lorsque la molécule arrive à proximité d'une autre, il devient nécessaire de repasser à l'échelle AA, pour prendre en compte les interactions entre elles, en particulier les interactions répulsives.

Pour gérer les reconstructions, plusieurs questions doivent être abordées. La première est que, lors d'une reconstruction, la molécule reconstruite doit remplacer "immédiatement" la molécule initiale. Ce qui signifie qu'au cours de la simulation, les deux molécules ne sont pas autorisées à coexister (leurs interactions ne correspondraient à rien). Dans l'approche réactive, ce remplacement, c'est-à-dire la destruction de la molécule initiale et l'introduction de la molécule reconstruite, est naturellement effectué entre deux instants.

Une seconde question concerne les possibilités mêmes de reconstructions. Il est difficile d'implémenter des mécanismes de reconstructions qui ne provoquent jamais aucun accroissement pathologique de l'énergie de la molécule reconstruite. Les mécanismes de reconstruction doivent être évalués en fonction de cet aspect.

Une troisième question est celle des interactions entre molécules à des échelles différentes. A ce propos, on peut introduire dans la simulation des mécanismes de *reconstruction automatique* déclenchés par la proximité des molécules, pour réaliser la mise à une même échelle des molécules interagissant.
## Reconstructions entre AA et UA

La reconstruction AA->UA consiste à construire une molécule à l'échelle UA à partir d'une molécule à l'échelle AA. Elle revient en fait à "gommer" les atomes d'hydrogène de l'échelle AA, en les "absorbant" dans les grains de l'échelle UA suivant l'équation :

$$C + H + H \rightarrow G_2 \qquad (13.1)$$

La masse est préservée, puisque la somme d'un carbone avec deux hydrogènes est 0.014 unités internes (la masse de C vaut 0.012 et celle de H vaut 0.01) qui est également la masse du grain $G_2$. On positionne le grain $G_2$ exactement sur l'atome C ce qui a pour effet de supprimer les énergies potentielles des deux liaisons CH et celle de l'angle HCH qui disparaissent dans la reconstruction. La somme des potentiels de la molécule AA est alors toujours plus importante que la somme des potentiels UA de la molécule reconstruite ce qui entraîne une perte systématique d'énergie potentielle lors de la reconstruction AA->UA.

La reconstruction inverse UA->AA consiste à construire une molécule AA à partir d'une molécule UA. Elle revient à remplacer les grains $G_2$ par des atomes de carbone avec deux atomes d'hydrogène, suivant l'équation :

$$G_2 \rightarrow C + H + H \qquad (13.2)$$

Comme dans l'équation (13.1), les masses sont cohérentes. On place les carbones exactement sur les grains $G_2$ et on cherche ensuite à placer les atomes d'hydrogène. Idéalement, on devrait placer ceux-ci de telle façon que l'énergie ajoutée (celle des liaisons CH, celle des angles HCH et celle des angles HCC) soit minimale. Pour cela, une phase de minimisation de la molécule AA préservant les positions des carbones devrait être lancée avant de reprendre la simulation. Une telle phase de minimisation globale est clairement très lourde et c'est pourquoi on n'a pas choisi de l'implémenter et de préférer une technique pragmatique mais non-sûre.

Dans cette technique, on greffe sur chaque carbone une paire de liaisons CH à l'équilibre avec deux hydrogènes formant avec le carbone un angle HCH également à l'équilibre (cf. Fig.13.1). Les deux liaisons CH et l'angle de valence HCH ainsi rajoutés n'introduisent pas d'énergie supplémentaire.

En ce qui concerne les angles de valence HCC et les angles de torsion HCCH et HCCC, on procède de même que lors de la création d'une molécule
<figure><img src="image_58.png" /><figcaption>FIG. 13.1 – Introduction dea hydrogènes, des liaisons HC et de l’angle de valence HCH.</figcaption></figure>

<figure><img src="image_59.png" /><figcaption>FIG. 13.2 – Introduction des angles de valence HCC et des angles de torsion HCCH et HCCC.</figcaption></figure>

à l'équilibre, en alternant les couples d'hydrogènes, dans des plans parallèles (cf. Fig.13.2).

Ainsi, contrairement à ce qui se passe avec la reconstruction AA->UA, de l'énergie supplémentaire peut apparaître, correspondant aux énergies des angles HCC, HCCH et HCCC introduits. Cette énergie supplémentaire peut même conduire à une explosion de la simulation en particulier lorsque des atomes d'hydrogène sont trop proches.

## Reconstructions entre UA et CG

La reconstruction UA->CG consiste à construire une molécule à l'échelle CG à partir d'une molécule à l'échelle UA. Elle revient en fait à regrouper deux grains UA en un seul grain CG suivant l'équation :

$$G_2 + G_2 \rightarrow CG_4 \qquad (13.3)$$

La masse d'un grain CG4 est de 0.028. Le grain CG4 est positionné au milieu de la liaison entre les deux grains G2 et une liaison CG est introduite entre deux grains CG4 consécutifs.

La reconstruction inverse CG->UA consiste à construire une molécule à l'échelle UA à partir d'une molécule CG. A chaque grain CG correspondent deux grains UA :

$$CG_4 \rightarrow G_2 + G_2 \qquad (13.4)$$
Les deux grains UA sont placés de part et d'autre du grain CG, de façon que celui-ci soit au milieu des deux grains UA. Une liaison UA est introduite entre les deux grains UA les plus proches. On adopte une technique similaire à celle utilisée pour la reconstruction UA->AA, dans le placement des grains UA, avec le même risque de mal positionner les grains UA en introduisant une énergie supplémentaire, pouvant faire "exploser" la simulation.

## Reconstructions entre AA et CG

Les reconstructions entre AA et CG peuvent être considérées comme des enchaînements de reconstructions entre AA et UA et entre UA et CG. La reconstruction AA->CG enchaîne deux reconstructions qui font systématiquement baisser l'énergie, et ne pose donc pas de problème d'explosion de l'énergie. Ce n'est par contre pas le cas de la reconstruction CG->AA inverse, pour laquelle une explosion peut survenir à la fois lors de l'introduction des grains UA et lors de l'introduction des atomes d'hydrogène durant la reconstruction UA->AA.

## Reconstructions manuelles

La Fig.13.3 montre une simulation dans laquelle l'utilisateur, en cliquant dans le panneau de contrôle, peut changer l'échelle de la molécule simulée (une molécule $C_{10}H_{20}$, à laquelle on a fourni initialement de l'énergie potentielle). Le passage de la molécule AA à la molécule UA a eu lieu en réponse au clic dans la partie ->ua du panneau. Le passage de la molécule UA à la molécule CG a eu lieu en réponse au clic dans ->cg. Les clics produisent des événements SugarCubes qui sont diffusés à tous les programmes s'exécutant.

## 13.1 Forces inter-moléculaires

Le traitement des forces inter-moléculaires est délicat lors des reconstructions pour deux raisons : la première est que l'introduction des hydrogènes lors des reconstructions UA->AA ou CG->AA peut conduire à des explosions, ceux-ci étant trop proches. Il en va de même, d'ailleurs, avec les grains UA lors des reconstructions CG->UA. La seconde raison est que les distances d'équilibre entre molécules ne sont pas les mêmes aux diverses échelles, ce qui
<figure><img src="image_60.png" /><figcaption>FIG. 13.3 – Haut/gauche : molécule AA initiale. Haut/droite : après exécution. Bas/gauche : après reconstruction de la molécule AA en UA. Bas/droit : après reconstruction de la molécule UA en CG.</figcaption></figure>

fait que les forces répulsives peuvent augmenter fortement lors des reconstructions de molécules trop proches. On va maintenant considérer ce cas.

La distance d'équilibre entre deux atomes (ou grains) régis par un potentiel Lennard-Jones 6-12 est :

$$rn = 2^{1/6}\sigma \qquad (13.5)$$

où $\sigma$ est la distance à laquelle le potentiel est nul. En effet, la force s'annule dans l'équation 3.44 lorsque $2(\sigma_r)^{12} - (\sigma_r)^6 = 0$ ce qui signifie, en posant
$X = (\sigma/r)^6$, que $2X^2 - X = 0$. D'où $X(2X - 1) = 0$ dont la seule solution valide est $X = 1/2$, c'est-à-dire $r = 2^{1/6}\sigma$.

À la distance d'équilibre, les forces répulsive et attractive entre les deux atomes s'équilibrent et le potentiel a la valeur minimale $-\epsilon$. Voici les valeurs de $rn$ pour les potentiels inter-atomiques associés aux alkanes :

$$rnCC = 0.39286171690828053$$

$$rnCH = 0.33202427388991257$$

$$rnHH = 0.28061551207734325$$

$$rnUA = 0.36025419440489326$$

$$rnCG = 0.35694293136238064$$

La valeur de $rnCG$ est obtenue à partir de celle de la courbe Lennard-Jones 6-12 de la Fig.12.3 dont le minimum coincide avec celui du potentiel inter-moléculaire CG.

Ces distances où les forces répulsives et attractives s'équilibrent ne sont pas les mêmes ce qui peut être problématique dans certains cas, par exemple lorsque deux molécules AA sont mises en vis-à-vis, comme sur la Fig.13.4 à une distance de 0.436 nm, distance à laquelle les forces répulsives et attractives globales s'annulent. La reconstruction AA->CG détruit l'équilibre des

<figure><img src="image_61.png" /><figcaption>FIG. 13.4 – Deux molécules $C_{10}H_{20}$ à distance de 0.436 nm. Les molécules sont à l'équilibre : la somme des forces attractives et la somme des forces répulsives sont égales.</figcaption></figure>

forces et résulte en une force répulsive qui éloigne les deux molécules CG à grande vitesse ce qui ne correspond certainement pas à la réalité.
## 13.2 Reconstructions automatiques

On introduit maintenant un mécanisme pour déclencher automatiquement les reconstructions en fonction de la proximité des molécules. Pour simplifier, on ne va considérer que les reconstructions CG->UA et UA->AA, les autres types de reconstructions pouvant être traitées de manière analogue.

Le mécanisme de déclenchement est une boucle exécutée à chaque instant qui mesure les distances aux autres molécules et qui génère un événement de reconfiguration lorsqu'une molécule est détectée à une distance inférieure à un seuil déterminé. Ce mécanisme est ajoutée aux molécules pouvant être reconstruites.

Les seuils de détection sont choisis de manière empiriques pour éviter autant que possible les explosions d'énergies :

$$detectUA = 0.39$$

$$detectCG = 0.71$$

Une molécule UA qui détecte que le plus proche atome de carbone d'une molécule AA est à une distance inférieure à detectUA se reconfigure automatiquement en AA. Une molécule CG qui détecte que le plus proche atome de carbone d'une molécule AA est à une distance inférieure à detectCG se reconfigure en AA. De même, lorsque la molécule CG détecte que le plus proche grain d'une molécule UA est à une distance inférieure à detectCG, elle se reconfigure en UA.

Il est à noter que les reconstructions ne peuvent concerner qu'un sous-ensemble des molécules simulées, ce qui a pour conséquence qu'il est nécessaire de définir les potentiels inter-moléculaires entre atomes et grains à des échelles différentes. Pour cela, on choisit l'échelle la plus haute (celles-ci étant ordonnées par AA < UA < CG) ce qui se traduit par le tableau suivant :

<table><thead><tr><th></th><th>AA</th><th>UA</th><th>CG</th></tr></thead><tbody><tr><th>AA</th><td>AA</td><td>UA</td><td>CG</td></tr><tr><th>UA</th><td>UA</td><td>UA</td><td>CG</td></tr><tr><th>CG</th><td>CG</td><td>CG</td><td>CG</td></tr></tbody></table>

Les Fig.13.5 et Fig.13.6 montrent une simulation de quatre molécules CG, composées chacune de deux grains. La situation initiale est celle de la Fig.13.5 où quatre molécules CG sont présentes. Les quatre molécules
<figure><img src="image_62.png" /><figcaption>FIG. 13.5 – Configuration initiale : quatre molécules CG situées aux quatre sommets d'un carré.</figcaption></figure>

se reconfigurent à l'échelle UA (image gauche/haut de la Fig.13.6). Assez rapidement, les molécules UA se détectent et se reconfigurent en AA (image droite/haut). Les deux dernières images (bas) montrent l'évolution des quatre molécules AA. La dernière image (bas/droit) indique que les molécules s'éloignent définitivement du centre de la figure, en se regroupant par paires.

## Conclusion

Concernant les reconstructions, on a adopté une approche pragmatique dans laquelle l'effectivité est privilégiée au détriment de la sûreté.

Les reconstructions présentées ne sont pas sûres dans le sens où elles peuvent échouer, en introduisant une trop forte énergie dans la molécule reconstruite. Il s'agit des cas où des grains d'une même molécule sont très proches, ce qui conduit à placer des grains ou bien des atomes dans des po-
<figure><img src="image_63.png" /><figcaption>FIG. 13.6 – Scénario des suites de reconfigurations des molécules initialement placées comme sur la Fig.13.5.</figcaption></figure>

sitions trop proches. Une reconstruction peut également échouer lorsque des grains de molécules différentes sont très proches et que les grains ou atomes introduits par la reconstruction produisent des forces inter-moléculaire trop importantes.

Par contre, les reconstructions présentées sont effectives dans le sens où, lorsqu'elles réussissent, l'énergie de la molécule initiale et celle de la molécule reconstruite diffèrent généralement peu. C'est particulièrement le cas des molécules avec peu d'énergie intra-moléculaire car dans ce cas, la reconstruction est proche de la création à l'équilibre de la molécule reconstruite.

Les reconstructions décrites ne sont donc pas adaptées aux molécules avec une grande énergie intra-moléculaire ou bien, en ce qui concerne les forces inter-moléculaires, aux molécules trop proches.

Il semble donc raisonnable de réserver les reconstructions aux cas où l'énergie des molécules n'est pas trop importante et lorsque les molécules sont suffisamment éloignées.

L'alternative à l'approche présentée, pour rendre les reconstructions sûres, serait d'introduire une phase de minimisation de la molécule reconstruite, un
processus qui peut prendre beaucoup de temps de simulation.

On pourrait envisager des reconstructions automatique AA->UA et UA->CG, lorsque les distances sont suffisamment grandes, ce qui serait un moyen supplémentaire pour accélérer les simulations. Ce type de reconstruction n'a pas été implémenté dans le système actuel mais il pourrait l'être sans grande difficulté.
# Chapitre 14

## Conclusion

L'approche réactive fournit un bon cadre de programmation des systèmes de DM. Elle possède une sémantique formelle claire et précise (Chap. 4). Elle permet de traiter simplement les aspects liés à la résolution temporelle, à la modularité de programmation, au déterminisme et à la dynamicité (struc-ture non figée des systèmes, pouvant évoluer dynamiquement en cours de simulation).

On s'est posé la question de la détermination des potentiels aux échelles UA et CG, pour les molécules alkanes. L'objectif était en fait de "déduire" les potentiels UA et CG à partir du potentiel AA.

La détermination des potentiels "coarse-grain" basée sur une méthode de physique statistique Boltzmann-inverse, traitant des données issues de simu-lations, bute sur plusieurs points, en particulier la dépendance des résultats à la température et le traitement incorrect des forces inter-moléculaires.

La détermination des potentiels par la méthode de minimisation se dif-férencie radicalement de la méthode Boltzmann-inverse. La minimisation ne repose pas sur des résultats de simulation, comme Boltzmann-inverse, mais sur des minimisations d'énergie à l'échelle AA (dans un sens, et bien que les minimisations soient effectuées également avec le même système de DM, il s'agit à proprement parler plus de *mécanique moléculaire* que de DM).

Un point fondamental est que les minimisations sont rendues possibles par l'établissement d'un lien géométrique fort avec l'échelle AA.

La stabilité du système de DM utilisé ici repose sur une explicitation com-plète des forces s'appliquant sur les atomes (Chap. 3) ainsi que sur l'emploi d'une résolution temporelle extrêmement stable (velocity-Verlet).
## Potentiel UA

Le potentiel UA déterminé par la méthode de minimisation possède plusieurs caractéristiques :

* Les potentiels UA déterminés ont la même forme que les potentiels AA correspondants : les potentiels UA de liaison et de valence sont harmoniques; le potentiel de torsion est un potentiel “triple-cosinus”; le potentiel inter-atomique est un potentiel Lennard-Jones.

* Le potentiel de liaison UA est identique au potentiel de liaison AA.

* Le potentiel de valence UA est très proche du potentiel de valence AA. On peut le voir comme prenant en compte l'angle CCC de valence des carbones, ainsi que les angles CCH associés.

* Le potentiel de torsion UA est la somme du potentiel CCCC et des potentiels HCCH partageant la même liaison CC centrale.

* Le potentiel inter-atomique UA est un potentiel Lennard-Jones "plus creusé" que le potentiel AA correspondant, prenant en compte les forces inter-moléculaires liées aux atomes d'hydrogène.

Ainsi, le potentiel UA est obtenu à partir du potentiel AA par un simple changement de paramètres.

## Potentiel CG

La détermination du potentiel CG est nettement plus difficile que celle du potentiel UA (bien qu'il n'y ait pas en CG de potentiel d'angle de torsion). La méthode de minimisation donne des potentiels qui n'ont pas la forme standard des potentiels AA ou UA et qui sont représentés sur les trois figures Fig.10.4, Fig.10.5, Fig.10.6. D'une certaine manière, le travail présenté dans ce texte peut être vu comme concentré dans ces trois courbes.

## Reconstructions

On a défini et implémenté des mécanismes de reconstruction, c'est-à-dire de changement d'échelle des molécules en cours de simulation. On a en particulier considéré les reconstructions automatiques en fonction de la proximité des molécules. Les reconstructions définies ne sont pas sûre dans le sens où elles préserveraient systématiquement l'énergie. Elles doivent donc être réservées aux molécules possédant une faible énergie intra-moléculaire ainsi
qu'aux molécules suffisamment éloignées pour ne pas induire des énergies inter-moléculaires trop importantes, lors des reconstructions.

## Perspectives

Plusieurs extension du travail présenté ici sont envisageables parmi lesquelle trois particulièrement.

Le premier axe de travail concernerait le “vrai parallélisme”, en utilisant une nouvelle version des SugarCubes (SugarCubesv5 [2]) dans laquelle l’utilisation de cartes graphiques (GPU) pour effectuer les calculs en parallélisme réel devient possible. L’utilisation de machines multi-processeurs serait bien sûr également d’un grand intérêt.

Le second axe de travail consisterait en l’implémentation d’une technique de reconstruction sûre, basée sur la minimisation des molécules reconstruites (cf. Chap.13).

Le troisième axe de travail concernerait l’extension des la classe des molécules prises en compte, en ne se limitant pas aux alkanes. Les “savons”, chaînes dont une extrémité contient des atomes d’oxygène semblent un bon candidat pour une telle extension. Les réactions chimiques provoquées par la présence d’atomes d’oxygène pourraient être naturellement prises en compte par des mécanismes très proches de ceux utilisés pour les reconstructions (de premières expérimentations ont été réalisées sur ce point).
Contexte

Ce texte a été rédigé par Frédéric Boussinot. Il décrit une partie du travail réalisé par Bernard Monasse et Frédéric Boussinot lorsqu'ils étaient chercheurs au CEMEF de l'École des mines de Paris. Les SugarCubes forment le travail de thèse de Jean-Ferdy Susini réalisé à l'Inria Sophia-Antipolis, dans l'équipe Mimosa. Merci à Pierre Montmitonnet pour ses nombreuses remarques, corrections et suggestions concernant ce texte.
Bibliographie

[1] DL_POLY. https://www.scd.stfc.ac.uk/Pages/DL_POLY.aspx.

[2] SugarCubesv5. http://cedric.cnam.fr/index.php/labo/membre/view?id=160.

[3] M. P. Allen and D. J. Tildesley. *Computer Simulation of Liquids*. Oxford, 1987.

[4] G. Berry and G. Gonthier. The Esterel Synchronous Programming Language : Design, Semantics, Implementation. *Science of Computer Programming*, 19(2) :87-152, 1992.

[5] F. Boussinot. Reactive C : An Extension of C to Program Reactive Systems. *Software Practice and Experience*, 21(4) :401-428, april 1991.

[6] F. Boussinot, B. Monasse, and J-F. Susini. Reactive programming of simulations in physics. *International Journal of Modern Physics C (IJMPC)*, 26(12) :1-16, 2015.

[7] F. Boussinot and J-F. Susini. The SugarCubes Tool Box - A Reactive Java Framework. *Software Practice and Experience*, 28(14) :1531-1550, december 1998.

[8] F. Boussinot and J-F. Susini. Java Threads and SugarCubes. *Software Practice and Experience*, 30(14) :545-566, 2000.

[9] W. Damm, A. Frontera, J. Tirado-Rives, and W. L. Jorgensen. OPLS All-Atom Force Field for Carbohydrates. *Journal of Computational Chemistry*, 18(16) :1955-1970, 1997.

[10] W. Damm, A. Frontera, J. Tirado-Rives, and W.L. Jorgensen. OPLS All-atom Force Field for Carbohydrates. *J. Comput. Chem.*, 16(18) :1955-1970, 1997.

[11] W.L. Jorgensen, J.D. Madura, and C.J. Swenson. Optimized Intermolecular Potential Functions for Liquid Hydrocarbons. *J. Am. Chem. Soc.*, (106) :6638-6646, 1984.

[12] A.D. MacKerell Jr., D. Bashford, M. Bellott, R. L. Dunbrack, Jr., J. D. Evansек, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy,
L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher III, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wirkiewicz-Kuczera, D. Yin, and M. Karplus. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. *J. Phys. Chem. B*, (102) :3586–3616, 1998.

[13] H. Krivine. *Comprendre sans prévoir, prévoir sans comprendre*. Cassini, 2018.

[14] L. Mandel and M. Pouzet. ReactiveML, A Reactive Extension to ML. In *ACM International conference on Principles and Practice of Declarative Programming (PPDP'05)*, Lisbon, Portugal, July 2005.

[15] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller. Equation of State Calculations by Fast Computing Machines. *Journal of Chemical Physics*, 21 :1087–1092, June 1953.

[16] R.A. Miron and K.A. Fichthorn. Multiple-time Scale Accelerated Molecular Dynamics: Addressing the Small-barrier Problem. *Phys. Rev. Lett.*, (93) :128301–128304, 2004.

[17] Bernard Monasse and Frédéric Boussinot. Détermination des forces à partir d'un potentiel en dynamique moléculaire (note). Technical report, November 2013. https://hal-mines-paristech.archives-ouvertes.fr/ hal-00880202.

[18] B. Mukherjee, L. Delle Site, K. Kremer, and Ch. Peter. Derivation of Coarse Grained Models for Multiscale Simulation of Liquid Crystalline Phase Transitions. *J. Phys. Chem. B*, (116) :8474–8484, 2012.

[19] University of Groningen. GROMACS.

[20] D. Reith, H. Meyer, and F. Muller-Plathe. CG-OPT: A Software Package for Automatic Force Field Design. *J. Comput. Phys. Commun.*, (148) :299–313, 2002.

[21] D. Reith, M. Putz, and F. Muller-Plathe. Deriving Effective Mesoscale Potentials from Atomistic Simulations. *J. Comput. Chem.*, (24) :1624–1636, 2003.

[22] K. Vanommeslaeghe, E. Hatcher, C. Acharya, S. Kundu, S. Zhong, J. Shim, E. Darian, O. Guvench, P. Lopes, I. Vorobyov, and A. D. MacKerell Jr. CHARMM General Force Field: A Force Field for Drug-like Molecules Compatible with the CHARMM All-atom Additive Biological Force Fields. *J. Comput. Chem.*, (31) :671–690, 2010.

[23] L. Verlet. Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules. *Phys. Rev.*, 159 :98–103, Jul 1967.

[24] A.F. Voter. A Method for Accelerating the Molecular Dynamics Simulation of Infrequent Events. *J. Chem. Phys.*, (106) :4665–4677, 1997.