La présente invention concerne un procédé de simulation d'un ensemble d'éléments, selon lequel le comportement des éléments est déterminé, à des pas de simulation successifs, sur la base d'un Hamiltonien associé au système d'éléments (la somme de l'énergie cinétique et de l'énergie potentielle de l'ensemble)

H = -^ p ^T .M ^~ p + V , p étant un vecteur indiquant les moments des éléments, V étant l'énergie potentielle du système et M ^"1 une matrice diagonale fonction des masses des éléments (dans des cas, cette matrice peut être fonction des positions des éléments).

L'énergie potentielle V est, dans certains cas, une fonction des positions des éléments seulement. Dans d'autres cas, l'énergie potentielle V peut dépendre également des moments des éléments. Les forces agissant sur les éléments peuvent être dérivées de cette énergie potentielle.

La simulation d'un ensemble d'éléments permet d'étudier le comportement d'un tel ensemble et d'en analyser les propriétés : les déplacements en termes de positions et de moments successifs des éléments, les corrélations des déplacements entre éléments, les changements de structure, les hausses et baisses d'interactions entre éléments, les configurations adoptées en moyenne, les évolutions des énergies associées, etc. Les éléments peuvent représenter des corps mécaniques, par exemple célestes ou fluides, des particules telles que des atomes ou molécules, par exemple des protéines, des fluides etc.

Une façon usuelle de simuler un ensemble d'éléments est de considérer le Hamiltonien de l'ensemble, d'en dériver des équations de mouvement, et de déduire le mouvement des éléments selon ces équations.

WO 2009/007550 décrit par exemple une technique de simulation d'un ensemble d'éléments.

Les évolutions de l'ensemble d'éléments doivent parfois être simulées sur une longue période, en vue de pouvoir observer certains phénomènes ou de pouvoir calculer certaines statistiques. Les temps de calcul, et le coût en calcul, de ces simulations deviennent alors parfois très importants. De nombreuses méthodes ont été proposées pour accélérer les simulations d'un ensemble d'éléments et la collecte des statistiques.

La présente invention vise à proposer une nouvelle solution pour réduire ces problèmes. A cet effet, suivant un premier aspect, l'invention propose un procédé de simulation d'un ensemble d'éléments du type précité caractérisé en ce que ledit procédé comporte les étapes selon lesquelles :

- on représente le système d'éléments en un arbre k-aire comportant des nœuds feuilles, lesdits nœuds feuilles représentant chacun un élément respectif, et des nœuds internes dont un nœud racine

- lorsque des conditions déterminées sont vérifiées, on impose, pour le pas de simulation courant, un même mouvement en translation à au moins les éléments descendant d'un nœud donné de l'arbre en définissant la matrice

M ^"1 comme étant égale à Φ ₈ où R est le nœud racine de l'arbre, étant donné la formule récursive selon laquelle, pour tout nœud A de l'arbre comportant k nœuds fils A _Y , A ₂ , A _k ,

Φ _Αι 0 0

Φ _Α = ^ PA -rΕ, - +r (i1 -- ρ _AΑ ))< 0 ^' · . 0

0 0 Φ

la matrice E étant une matrice de taille dn _A ^* dn _A formée de n _A ^* n _A blocs de taille d ^* d égaux à la matrice identité de dimension d, n _A étant égal au nombre d'éléments descendant du nœud A, d étant la dimension de l'espace dans lequel évoluent les particules, m _A est la somme de la masse de ces n _A éléments descendant du nœud A, p _A étant une fonction de restriction du nœud A comprise entre 0 et 1 , dont la valeur est égale à 1 lorsque A est un nœud feuille ou lorsqu'un même mouvement a été imposé aux éléments descendants du nœud donné ; dans le cas d'un nœud feuille A _i , la matrice

Φ _Α est égale à la masse inverse de la particule représentée par le nœud A _i multipliée par la matrice identité de dimension d.

L'invention permet d'effectuer des simulations qui requièrent un plus petit volume de calcul et, par conséquent, moins de temps de calcul, pour déterminer le comportement des éléments selon ces simulations, par exemple l'énergie potentielle, les forces appliquées sur les éléments, les positions et/ou les moments des éléments.

Dans des modes de réalisation, le procédé de simulation d'un ensemble d'éléments suivant l'invention comporte en outre une ou plusieurs des caractéristiques suivantes : on impose, pour le pas de simulation courant, un même mouvement en translation à au moins les éléments descendants d'un nœud donné de l'arbre en fonction de la valeur d'une fonction des moments desdits éléments ; on impose, pour le pas de simulation courant, un même mouvement en translation à au moins les éléments descendants d'un nœud donné de l'arbre en fonction de la valeur prise par , ou

m _A est la somme de la masse des éléments du nœud A _; , et p _A ^s _i est le vecteur- somme des moments des éléments du nœud fils A _{ , \\.\\ est la norme du vecteur, C est une constante positive ;

- pour le pas de simulation courant pour un nœud interne A:

- si ε _Α est inférieur à un premier seuil, on fixe p _A égal à 1 pour imposer un même mouvement en translation aux éléments descendants du nœud A ;

- si ε _Α est supérieur à un second seuil supérieur au premier seuil, on fixe p _A égal à 0 ;

- pour le pas de simulation courant, si ε _Α est compris entre le premier seuil et le second seuil, on fixe p _A égal à la valeur prise par une fonction d'interpolation d'ordre 5 de la variable ε _Α ;

- le procédé comporte une étape de détermination des valeurs d'au moins une information, à des instants de simulation successifs sur la base dudit Hamiltonien, ladite étape tirant parti du fait que les valeurs de l'information relatives aux éléments auxquels un même mouvement en translation a été imposé pour l'instant de simulation courant, dépendent de la position relative desdits éléments et sont par conséquent inchangées ;

- l'information relative audit élément comprend l'énergie potentielle dudit élément et/ou la force d'interaction appliquée audit élément.

Suivant un deuxième aspect, la présente invention propose un programme d'ordinateur de simulation d'un système d'éléments, comprenant des instructions logicielles pour mettre en œuvre les étapes d'un procédé suivant le premier aspect de l'invention lors d'une exécution du programme par des moyens de calcul.

Ces caractéristiques et avantages de l'invention apparaîtront à la lecture de la description qui va suivre, donnée uniquement à titre d'exemple, et faite en référence aux dessins annexés, sur lesquels : - la figure 1 représente une simulation de trajectoires d'un système à 2 particules effectuée avec le Hamiltonien H standard ;

- la figure 2 représente des simulations de trajectoires d'un système à 2 particules avec le Hamiltonien adaptatif relatif H _A R selon l'invention ;

- la figure 3 représente un dispositif mettant en œuvre un mode de réalisation de l'invention ;

- la figure 4 est un organigramme des étapes d'un procédé dans un mode de réalisation de l'invention ;

- la figure 5 représente un organigramme des sous-étapes de l'étape 101 de la figure 4.

Considérons la simulation d'un ensemble E de N particules a,, i=1 à N.

Le Hamiltonien H(p, q) associé à l'ensemble E (cf. par exemple « Understanding molecular simulation : from algorithms to applications », Frenkel D., Smit B.) s'écrit souvent de la façon suivante :

p étant un vecteur indiquant les moments de toutes les particules, q un vecteur indiquant les positions de toutes les particules, M ^"1 une matrice diagonale fonction des masses de ces particules.

V(q) est le potentiel d'interaction entre les N particules ; il est fonction de leurs positions et on le considérera comme indépendant des moments.

Dans un espace en 3 dimensions par exemple, avec un référentiel de coordonnées (X, Y, Z), le moment p, de chaque particule a,, i=1 à N, s'écrit (p _{i x} , p _{i y} , p ) et la position q, de chaque particule a,, i=1 à N, s'écrit (q _{i x} , q _{i y} , q ).

P V

Les vecteurs p et q s'écrivent donc et q =

ΡΝ,Ζ Usuellement la matrice M ^"1 utilisée dans l'art antérieur est une matrice diagonale de dimension 3N ^* 3N, dont les termes M[3i-2, 3i-2] = M[3i-1 , 3i-1 ]= M[3i, 3i] = m _i , pour i=1 à N, où m _i est la masse de la particule a,.

Il s'agit là de la définition usuelle du Hamiltonien H, qu'on appellera ci-dessous, Hamiltonien standard.

Selon l'invention, on définit un Hamiltonien dit Hamiltonien adaptatif relatif H _AR , ainsi :

1 T

H _AR (p,q) = -p ¹ &(p,q) .p + V(q) , (formule 1 ) dans lequel (ρ,ο) , matrice bloc-diagonale 3N ^* 3N dite matrice inverse de masse adaptative relative, remplace M ^"1 et dépend du vecteur p, et éventuellement du vecteur q.

Plus précisément, les valeurs prises par Φ(ρ, ο) sont déterminées à chaque pas de simulation, la fonction Φ ayant été choisie de manière à pouvoir restreindre au moins un sous-ensemble des particules à suivre ensemble un même mouvement à un ou à certains pas de simulation. Les particules qui font partie du sous-ensemble sont déterminées par vérification d'une condition, par exemple visant à vérifier qu'avant le pas de simulation courant, elles étaient animées d'un mouvement proche.

Ces dispositions résultent en l'apparition d'un ou plusieurs sous-ensembles de particules, les particules d'un sous-ensemble étant animées du même déplacement en translation pendant le pas de simulation considéré (i.e. le vecteur joignant deux quelconques de ces particules reste constant pendant une plage de temps dans laquelle ces particules ont le même déplacement).

Les particules d'un tel sous-ensemble peuvent ainsi être assimilées, pour le pas de simulation considéré, comme faisant partie d'un corps rigide unique.

Les informations qui dépendent seulement des positions relatives des particules, par exemple, les forces d'interaction entre les particules, n'ont donc pas besoin d'être actualisées à l'issue du pas de simulation considéré au sein d'un même sous-ensemble rigide.

Lorsque le mouvement de particules d'un sous-ensemble n'est plus restreint entre elles par la matrice de masse inverse adaptative, chaque particule reprend son mouvement propre.

Dans un mode de réalisation de l'invention, une transition continue est appliquée entre les deux types de comportements (mouvement restreint/mouvement libre). De cet Hamiltonien adaptatif H _A R, on déduit des équations adaptatives de mouvement, définissant p et q qui sont les dérivées des vecteurs p et q par rapport au temps t.

Par exemple, dans un cas d'implémentation d'une simulation dans un ensemble NVE (ensemble E à nombre de particules, à volume et énergie constants) considérée ici à titre d'illustration, la valeur du Hamiltonien (adaptatif selon l'invention ou standard) est constante dans le temps, et les équations adaptatives du mouvement sont : dp ₌ dH _AR _ dV l _pT d<î>(p,q)

dt dq dq 2 dq

(formules 2)

Les exemples ci-après sont relatifs à un ensemble NVE, mais bien sûr l'invention peut être appliquée à d'autres ensembles, par exemple, l'ensemble NVT (ensemble E à nombre de particules, à volume et température constants).

Plusieurs solutions distinctes peuvent être utilisées pour définir les particules constituant un sous-ensemble rigide, sans avoir à tester toutes les possibilités de combinaison entre les particules de l'ensemble.

Une première solution, est d'imposer des contraintes sur des paires déterminées de corps rigides qui peuvent être considérés comme candidats à fusionner en un même corps rigide. Par exemple, en notant A _{1 ;} A _n , un ensemble de corps rigides indexés, on peut organiser les indices en un arbre binaire, de sorte que ne peut fusionner qu'avec A ₂ en un corps rigide (A _! +A ₂ ), quand une condition spécifique de fusion est satisfaite, A ₃ ne peut fusionner qu'avec A ₄ etc., et à un niveau supérieur dans l'arbre binaire, (A _! +A ₂ ) ne peut fusionner qu'avec (A ₃ +A ₄ ) etc. Un tel arbre comporte au rang hiérarchique maximum, un nœud racine, nommé nœud R, et au rang hiérarchique minimum, les nœuds feuilles correspondant aux particules considérées seules. Tous les nœuds de l'arbre, sauf les nœuds feuilles, s'appellent les nœuds internes.

Une autre solution est d'imposer des contraintes sur des sous-ensembles déterminés de corps rigides, comportant éventuellement plus de deux corps rigides. On pourrait ainsi organiser les indices en un arbre n-aire, de sorte que Ai ne peut fusionner qu'avec A ₂ , A ₃ ,..., A _n en même temps, etc.

Une autre solution est de considérer le graphe déduit d'une recherche des voisins, les nœuds du graphe étant des particules ou des corps rigides, et une arête entre deux nœuds indique que la distance entre les deux nœuds reliés par l'arête est inférieure à une distance-seuil.

Le mode de réalisation détaillé ci-dessous est basé sur une structure d'arbre binaire telle que mentionnée relativement à la première solution ; la structure d'arbre a l'avantage d'autoriser des mises à jour incrémentales et récursives des métriques de décision pour les nœuds internes dans la hiérarchie, et des mises à jour incrémentales, par exemple pour les forces et énergies partielles.

Dans un tel arbre, chaque nœud feuille représente une particule seule de l'ensemble, et chaque nœud interne représente un sous-ensemble de particules composé des particules représentées par les nœuds fils dudit nœud interne. On nommera corps C l'ensemble comprenant le nœud C et les descendants, directs ou indirects, du nœud C.

Selon l'invention, si un corps C, i.e. un nœud C interne de l'arbre, est composé de deux corps A et B (i.e. le nœud C comporte deux nœuds fils A et B), on définit la matrice de masse inverse adaptative Φ ₀ correspondant à ce nœud C avec la formule récursive suivante :

(formule 3)

où :

- la matrice de masse inverse adaptative Φ _Α correspondant au nœud A et la matrice de masse inverse adaptative Φ _Β correspondant au nœud B, sont elles- mêmes définies selon la formule 3 si les nœuds A et B ne sont pas des nœuds feuilles ; dans le cas d'un nœud feuille A, respectivement B, la matrice Φ _Α , respectivement Φ _Β , est égale à la masse inverse de la particule représentée par le nœud A, respectivement le nœud B, multipliée par la matrice identité de dimension 3 ;

- la matrice E est une matrice de taille 3n _c ^* 3n _c , formée de n _c ^* n _c blocs 3 ^* 3 égaux à la matrice identité de dimension 3, n _c étant égal au nombre de nœuds feuilles du nœud C (i.e. l'ensemble des nœuds feuilles dans les nœuds fils directs du nœud C ou dans la descendance de ces nœuds fils) ;

- m _c est la somme de la masse des particules représentées par ces n _c nœuds feuilles du nœud C ;

p _c est la fonction de restriction du corps C.

Ainsi on considère la fusion de deux corps rigides A et B en un corps rigide C (et inversement l'éclatement du corps rigide C composé des deux nœuds rigides A et B) sur la base de p _c , la fonction de restriction du corps C. Plus exactement, un corps C ou nœud interne C, est considéré comme rigide (i.e. toutes les particules feuilles comprises dans le nœud C ou les descendants du nœud C, sont animées du même mouvement en translation) lorsque la fonction p _c prend la valeur

1 , et sinon il est considéré comme non-rigide.

Tous les nœuds feuilles sont considérés comme rigides par définition, et ce tout le temps que dure la simulation, puisqu'ils sont constitués chacun d'une seule particule.

Lorsque la fonction p _c a pour valeur 0, le nœud interne C est considéré comme libre.

Dans le mode de réalisation de l'invention, la fonction p _c est définie récursivement, de façon à lisser le basculement entre les valeurs 0 et 1 , selon la formule suivante :

1 si C est un noeud feuille,

(formule 4)

[μ(ε )ρ _Α ρ _Β dans le cas contraire, où μ(ε)€ [θ,ΐ] est une fonction deux fois dérivable de ε

1 si 0≤e≤e ^r ,

μ(ε) 0 si ε≥ ε ^f , (formule 5) s(£ ) G [θ,ΐ] autrement, avec les seuils ε' et ε ^Γ définis de sorte que ε' > ε ^Γ et s étant une fonction deux fois dérivable de ε .

La nature deux fois dérivable de μ permet de préserver la stabilité de la simulation de l'ensemble E de particules.

Par exemple, une forme possible pour ε(ε) est une fonction d'interpolation d'ordre

5, par exemple égale à - 6η ⁵ + 15η ⁴ - 10η ³ + 1 , avec η = ou

δ = ε' - ε ^Γ représente la largeur de la région de transition entre les comportements rigide et libre du corps.

Dans le mode de réalisation de l'invention, ε ₀ est choisi dépendant d'une relation du moment des deux corps A et B, avec

s 2 s 2 s . s 2 2

1 PA 1 PB 1 PA + PB _ 1 PA ^SM B ^~ P ^M A

+- (formule 6)

2 m _A 2 m _B 2 m _c 2 m _A m _B m _c où A et B sont les nœuds fils du nœud C,

où le vecteur p _c ^s est la somme des moments de toutes les particules qui sont des nœuds feuilles dans le corps correspondant au nœud C : p _c ^s =∑Ρι , et est la somme ieC

des vecteurs p _A ^s et p _B ^s : p _c ^s = p _A ^s + p _B ^s -

Cette invention peut être généralisée pour un arbre k-aire. Dans ce cas, pour le nœud parent avec k nœuds fils A 1 < i≤ k , la matrice inverse de l'inertie est

Le coefficient ½ dans cette formule et la formule 6 peut être remplacé par une autre constante, par exemple 1 , ce qui ne changera pas la simulation si les seuils ε' et ε ^Γ sont modifiés de la même façon, i.e. multipliés par un facteur égal au double de l'autre constante).

Cette invention peut être généralisée aussi différemment : on peut choisir la fonction de restriction p _c différente pour des nœuds internes différents. Par exemple, jusqu'à certain niveau de l'arbre p _c est définie selon les formules 4 et 5, donc, les nœuds peuvent être réunis ensemble et séparés, et dans les niveaux de l'arbre supérieurs à ce niveau, p _c est toujours égale à zéro. Dans un mode de réalisation particulier, il est possible de définir un sous-ensemble du système E en permanence actif (les particules de ce sous-ensemble ne sont jamais fusionnées ni entre elles, ni avec le reste du système, donc elles suivent toujours le long de la simulation un mouvement non restreint, libre (un tel mode de réalisation de l'invention peut par exemple être appliqué à un polymère dans un solvant ) : toutes les particules de ce sous-ensemble devront être placées dans un sous-arbre séparé et les valeurs p _c devront être fixées à 0 pour tous les nœuds internes de ce sous-arbre et ne jamais être mises a jour. Dans tous les cas, les fonctions p _c doivent être définies pour tous les nœuds internes avant le début de la simulation et ne doivent pas changer pendant la simulation afin d'avoir une simulation stable. Selon l'invention, on choisit la matrice inverse de masse adaptative relative Φ utilisée dans la formule (1 ) du Hamiltonien adaptatif relatif H _A R, égale à Φ _¾ qui est la matrice définie conformément à la formule récursive 3 relativement au nœud racine R de l'arbre représentant l'ensemble E de particules : Φ(ρ,<?) = (ρ) = _Κ ■

1 τ

On a alors le Hamiltonien H _AR , tel que H _AR (p,q) =— p & _R (p).p + V(q) , qui est séparable car la matrice Φ _κ ne dépend que des moments (et pas des positions) des particules.

Les équations de mouvement définies par les formules 2 s'expriment alors ainsi :

^ _ dp _ dH _AR _ d V _^

dt dq dq '

(formules 7)

. dq dH _AR , . 1 _T 3Φ _κ (ρ)

at dp 2 dp

Il ressort en outre clairement de la formule 3 que si le nœud C est libre { p _c =0), alors la matrice Φ ₀ est composée de deux blocs sur la diagonale correspondant l'un au nœud fils A et l'autre au nœud fils B.

Dans un premier exemple de simulation selon l'invention, il est considéré un système à une dimension comprenant deux particules P1 , P2 de masse 1 connectées l'une à l'autre par un ressort de raideur 1 . L'arbre binaire correspondant au système est un nœud racine C dont les deux nœuds fils A et B correspondent à ces deux particules. Etant donné des moments initiaux, les deux particules oscillent dans l'espace avec le temps.

Les trajectoires d dans l'espace calculées pour ces particules selon une simulation de l'art antérieur sur la base du Hamiltonien H standard sont représentées en figure 1 en fonction du temps de simulation.

Les trajectoires d calculées pour ces particules selon une simulation sur la base du Hamiltonien adaptatif relatif H _A R selon l'invention sont représentées en figure 2 en fonction du temps de simulation, pour différentes valeurs prises par δ et ε ^f . Ainsi comme indiqué sur ces trajectoires, à certains moments de la simulation (parfois dès le début), les moments des deux particules deviennent proches selon la formule 6. Ceci survient lorsque le ressort est presque compressé et presque décompressé. Selon l'invention, le nœud racine C devient alors rigide. En conséquence, dans un cas général, les nœuds fils du nœud C rigide suivent le même mouvement de translation, dans le cas présent, du fait de la préservation des moments, les particules s'arrêtent. Les trajectoires des 2 particules deviennent donc parallèles. Les moments des particules continuent à évoluer cependant, ce qui a pour conséquence qu'à un pas donné de la simulation, la condition de rigidité du nœud C n'est plus satisfaite et que les deux particules reprennent leur mouvement propre. Cela se reproduit donc périodiquement au cours de la simulation.

Comme illustré par les trajectoires, pour une même valeur de ε ^f avec une valeur accrue pour δ , la région de transition de la trajectoire est plus large. Pour une même valeur de δ avec une valeur accrue pour ε ^f , la région de rigidification est plus longue et plus plate.

Ainsi selon l'invention, la matrice Φ , définie égale à Φ _¾ spécifie comment, et quand, des degrés de liberté en position relatives de sous-ensemble(s) de particules sont activés ou désactivés pendant la simulation.

Dans un mode de réalisation de l'invention, un dispositif informatique 1 représenté en figure 3 est utilisé pour mettre en œuvre une simulation d'un ensemble E de N particules, selon les principes de l'invention explicité ci-dessus.

Ce dispositif 1 comprend un ordinateur comprenant notamment une mémoire 2 adaptée pour stocker des programmes logiciels et des valeurs de paramètres calculées successivement décrites ci-dessous (forces d'interactions globales, partielles, potentiel d'interaction, positions, moments...), un microprocesseur 3 adapté pour exécuter les instructions de programmes logiciels et notamment du programme P décrit ci-dessous, et une interface homme/machine 4, comprenant par exemple un clavier et un écran, respectivement pour saisir des instructions d'un utilisateur et pour afficher des informations à destination de l'utilisateur, par exemple des courbes telles que celles illustrées en figure 2.

Dans le mode de réalisation de l'invention considéré, la mémoire 2 comprend le programme P simulant le comportement de l'ensemble E de particules a, i =1 à N sur la base de l'intégration dans le temps des équations correspondant aux formules 7.

Comme le Hamiltonien adaptatif relatif H _A R selon l'invention est séparable, de nombreuses techniques d'intégration sont utilisables.

En référence à la figure 4, le programme P comprend des instructions logicielles qui, lorsqu'elles sont exécutées sur le microprocesseur 3, sont adaptées pour mettre en œuvre l'étape préliminaire 100 et itérativement les étapes 101 à 104 lors d'une n+1 ' ^eme itération du programme P correspondant à l'instant de calcul h _n+1 = h ₀ + (n+1 )h, avec n entier≥0, h étant le pas de temps de simulation.

Dans une étape préalable 100 d'initialisation, un arbre binaire représentant les particules et les possibilités de rigidification de deux en deux, comme indiqué plus haut, est construit pour le système E, des valeurs initiales de moment et de position sont fixées pour les particules. Seuls les nœuds feuilles sont alors rigides.

Dans une étape 101 , les forces d'interaction s'exerçant entre les particules à la n+1 ' ^eme itération, sont déterminées.

Dans une étape 102, les moments des particules à la n+1 ' ^eme itération, sont déterminés.

Dans une étape 103, les métriques relatives aux nœuds internes de l'arbre à la n+1 ' ^eme itération, sont déterminées. Elles comprennent notamment p _c et e _c .

Dans une étape 104, les positions des particules à la n+1 ' ^eme itération, sont déterminées.

Toutes ces étapes d'actualisation des valeurs des forces, moments, métriques, positions, pour l'itération courante, sont mises en œuvre en fonction des valeurs calculées lors de l'itération précédente.

La mise en œuvre de ces étapes 101 à 104 est détaillée ci-après.

Dans l'étape 101 , l'actualisation de la valeur des forces d'interaction peut être accélérée, par rapport à l'art antérieur, en prenant en compte les rigidifications selon l'invention, si les forces d'interaction dépendent uniquement des positions relatives des particules.

Ainsi selon l'invention, il n'est pas nécessaire de recalculer les forces au sein de sous-ensembles de particules identifiés comme rigides (d'après la valeur prise par e _c conformément à la formule 6) à l'itération de simulation précédente (les forces sont déterminées avant e _c ). Pour efficacement mettre à jour les forces en tirant parti de cet avantage consécutif aux rigidifications selon l'invention, un algorithme incrémental de mise à jour des forces peut être utilisé.

Cet algorithme est basé, premièrement, sur une hiérarchie de volumes englobants

(un volume par nœud de l'arbre binaire), par exemple, sur la hiérarchie des AABB's (en anglais « Axis-Aligned Bounding Volumes », les volumes englobant le sous-système et ayant une forme de parallélépipède aux cotés parallèles aux axes des coordonnées, cf. par exemple Grudinin, S. and Redon, S., "Practical modeling of molecular Systems with symmetries", Journal of Computational Chemistry 31 , 9 (2010), pp. 1799—1814).

Cet algorithme de mise à jour des forces utilise deuxièmement des tables de forces partielles. Une table de forces partielles correspondant à chaque nœud interne C comprend des vecteurs 3D, dont le nombre est égal au nombre de nœuds feuilles descendant directement ou indirectement du nœud C. Chaque vecteur de cette table est la force d'interaction agissant sur la particule représentée par un de ces nœuds feuilles en provenance de toutes les autres particules représentées par les autres nœuds feuilles descendant de C. Chaque vecteur de la table de forces partielles du nœud racine R correspond à la force globale agissant sur une particule. Une table de forces partielles correspondant a un nœud feuille est un vecteur nul.

La mise à jour des forces à chaque itération comprend trois sous-étapes principales (l'énergie potentielle peut être mise à jour en même temps) :

Dans une étape 101 1 , la hiérarchie des volumes englobants est mise à jour, par exemple à partir du bas de l'arbre vers le haut (i.e. en partant des nœuds feuilles vers la racine). Pour chaque nœud feuille, la boîte AABB mise à jour coïncide avec la dernière position déterminée de la particule représentée par le nœud feuille. La boîte de chaque nœud interne est recalculée par exemple en englobant les boîtes de ses deux nœuds fils.

La structure de l'arbre est donc inchangée, les volumes englobants sont eux mis à jour selon les positions des particules qui ont changé entretemps.

Dans une étape 101 _2 , en utilisant cette hiérarchie mise à jour des AABB's, des listes d'interaction sont construites de façon récursive pour les nœuds internes de l'arbre.

Pour chaque nœud C, ayant deux nœuds fils A et B, la liste d'interactions ainsi construite contient toutes les paires de particules telles qu'une des particules de la paire appartient au corps A (i.e. est représentée par un nœud feuille descendant directement ou indirectement du nœud A), tandis que l'autre appartient au corps B. Cette liste d'interaction peut être construite pour chaque nœud interne comme décrit pour d'autres types de volumes englobants dans l'article R. Rossi, M. Isorce, S. Morin, J. Flocard, K.

Arumugam, S. Crouzy, M. Vivaudou, and S. Redon. "Adaptive torsion-angle quasi-statics: a gênerai simulation method with applications to protein structure analysis and design".

Bioinformatics 2007 23(13):i408-i417.

Lors de la première itération, les listes d'interactions sont construites pour tous les nœuds internes. Pour toutes les autres itérations, dans les nœuds rigides ces listes ne changent pas et donc peuvent ne pas être recalculées.

Dans une étape 101_3, une mise à jour incrémentale des tables de forces partielles est réalisée. Lorsqu'un nœud est identifié comme nœud rigide, les forces entre les particules du nœud n'ont pas changé. Par conséquent, seules les tables de forces partielles des nœuds non-rigides doivent être actualisées.

Elles peuvent être actualisées par exemple comme suit récursivement, en partant des nœuds feuilles vers le nœud racine.

Pour chaque nœud C non-rigide, dans la première moitié de la table de forces partielles relative au nœud C, les éléments de la table de forces partielles relative au nœud A sont copiés et dans la seconde moitié de la table de forces partielles relative au nœud C, les éléments de la table de forces partielles relative au nœud B sont copiés. Puis, pour chaque paire de la liste d'interaction établie à la sous-étape 101_2 pour le nœud C, les forces d'interaction entre les deux particules de chaque paire de cette liste sont calculées et ajoutées aux forces correspondantes de la table de forces partielles. Au final, les forces prenant en compte toutes les particules sont stockées dans la table de forces partielles du nœud racine R.

Les étapes 101_2 et 101_3 peuvent être combinées, par exemple dans une seule traversée de l'arbre, pour effectuer la mise à jour de la table de forces partielles correspondant à un nœud C dès lors que la liste d'interactions correspondant audit nœud C est disponible.

Les autres informations dépendant uniquement des positions relatives des particules peuvent être actualisées de façon similaire.

L'étape 102 de mise à jour des moments des particules est réalisée en utilisant une méthode classique.

Ainsi, pour chaque particule a,, en libellant p, _n+l la valeur du moment p, de la particule a,, lors de l'itération courante n+1 , cette valeur est obtenue par la formule suivante :

h étant le pas de temps de la simulation et f, _n+i étant la force d'interaction globale s'exerçant sur la particule a, _, i =1 à N, qui est due aux interactions exercées par toutes les autres particules du système E, figurant dans la table de forces partielles du nœud racine R telle qu'elle vient d'être déterminée à l'itération courante.

La complexité de la mise à jour « naïve », à chaque pas de temps, de la position des particules à partir de l'équation q = Φ _Κ (p)p +—p ^T ^- ^-p , est 0(N ³ ) en raison du

2 dp

ΘΦ ( p)

terme—-^— . Ce terme est une dérivée d'une matrice par rapport à un vecteur, donc dp

un objet à trois dimensions, où, pour i= 1 à N :

Dans le mode de réalisation considéré, la complexité peut être réduite à O(N.logN) par une mise à jour récursive depuis les nœuds feuilles de structures de données Q _c , Rc et S _c pour chaque nœud C interne avec deux nœuds fils A et B.

En remarquant que

ou p _A est le vecteur des moments de toutes les particules qui sont des nœuds feuilles du nœud A et p _B est le vecteur des moments de toutes les particules qui sont des nœuds feuilles du nœud B, le terme Q _c = <î> _c (p _c )p _c correspondant au premier terme de l'équation donnant q :

(formule 9)

où {x} est le vecteur de dimension égale à n _c et dont chacune des n _c composantes est égale à l'élément x.

De même, en considérant pour le nœud C, le terme R _C = dPc (Pc )

dPc

dPc (Pc )

R _C = avec p _c = μ(ε )ρ _Α ρ _Β , en différenciant,

p

(formule 10)

Le dernier terme de l'équation donnant q peut être mis à jour récursivement en utilisant la formule 9.

En effet :

g - r * c (Pc ) _

^J c ^_ Pc -, Pc ^~

àp _c

rmule 1 1 )

La fonction (·) représente le produit scalaire.

Ainsi l'étape 103 comporte une mise à jour des structures Q _c , R _c et S _c pour chaque nœud C comportant deux nœuds fils A et B.

On notera que l'étape 100 comporte en outre une étape d'initialisation de ces structures où pour chaque nœud interne C : p _c =0 et les valeurs de Q _c , Rc et S _c sont également nulles, et pour chaque nœud feuille F représentant une particule : p _F =1 ; Q _F = P _F l m _F ; S _F = 0 , R _F = 0 .

Dans l'étape 103, la mise à jour des métriques Q _c , Rc et S _c , ainsi que de e _c et p _c relativement au nœud C est effectuée récursivement, le long de la hiérarchie des nœuds depuis la racine:

- Si le nœud C est interne, avec deux fils A et B :

- faire la mise à jour des métriques du nœud A ;

- faire la mise à jour des métriques du nœud B ;

- mettre à jour e _c à l'aide de la formule 6 ;

- mettre à jour p _c : p _c = μ(ε _ε )ρ _Α ρ _Β ;

- mettre à jour Q _c à l'aide de la formule 9 ;

- mettre à jour R _c à l'aide de la formule 10 ;

- mettre à jour S _c à l'aide de la formule 1 1 ;

- Si le nœud C est un nœud feuille représentant une particule a _F : QF = PF NIF■ Pour optimiser les calculs, la mise à jour des vecteurs Q _c , Rc et S _c peut être effectuée seulement pour les nœuds pour lesquels p _c > 0 parce que ces métriques peuvent ne pas être utilisées pour les nœuds libres, il n'y a donc pas besoin de les mettre à jour.

Dans le cas général, ces métriques peuvent être mises à jour autrement qu'en traversant l'arbre, pour peu qu'elles soient mises à jour dans les sous-arbres avec le nœud rigide en haut.

Une fois les métriques déterminées, les positions des particules sont mises à jour dans l'étape 104.

La mise en œuvre de l'étape 104 peut être la suivante pour la mise à jour de la position d'un nœud C, chaque nœud interne C étant considéré comme comportant deux nœuds fils A et B :

- si p _c = 0 (le nœud C est un nœud libre) :

- mettre à jour la position du nœud A ;

- mettre à jour la position du nœud B ;

- sinon, pour chaque nœud feuille k descendant du nœud C, directement ou indirectement : q. _n+1 = q _n + (Q _c [k] + 0.5S _c [k]) h, q _{i n+i} étant la valeur de la position de la particule a,, pour l'itération n+1 considérée et h le pas de temps de simulation, Q _c [k] étant le k ^eme élément du vecteur Q _c et S _c [k] étant le k ^ème élément du vecteur S _c .

La position des particules est ainsi déterminée.

Dans cette réalisation, la correspondance biunivoque entre l'identifiant de la particule i et son numéro d'ordre parmi les nœuds fils du nœud C doit être établie, par exemple pendant l'initialisation de la simulation, pour toutes les particules et tous les nœuds internes de l'arbre.

Dans cette réalisation de mise à jour des positions, les métriques des nœuds pour lesquels p _c > 0 ne sont pas utilisés et donc peuvent ne pas être calculées. Si toutes les métriques sont quand même mises à jour pendant l'étape 103, les positions de toutes les particules a _t , l≤i≤N peuvent aussi être déterminées comme suit :

<? _;,«+i ⁼ 1i + (Q _R [k] + 0-5S _R [k])h , où l'index R correspond au nœud racine.

Pour chaque nœud interne C de l'arbre, 4 vecteurs de longueur égale au nombre de nœuds feuilles descendants directs ou indirects du nœud C, sont stockés : Qc, Rc et S _c et une table de forces partielles. La complexité en espace de mise à jour des positions est par conséquent O(N.logN). La complexité en temps est aussi O(N.logN) puisque ces structures doivent être mises à jour à chaque pas de temps, et que Q _c est toujours mise à jour pour tous les nœuds, y compris les nœuds feuilles.

Dans le cas général, les positions peuvent être mises à jour autrement qu'en traversant l'arbre.

Puis, si la durée maximum de la simulation n'est pas atteinte, on réalise une nouvelle itération du programme P.

A titre d'illustration, quatre simulations en 2 dimensions de l'évolution d'un ensemble NVE de N=5930 particules a, , i=1 à N, avec une masse de 1 g/mol chacune, en utilisant un potentiel Lennard-Jones ( E _m / k _B = 120 Kelvin, où E _m est le minimum d'énergie, distance d'équilibre S = 3.4 ângstroms, distance de coupure 8 ângstroms, le potentiel étant tronqué de manière lisse entre 7.5 et 8 ângstroms), ont été réalisées partant d'un choc déclenché par l'envoi à haute vitesse d'une particule sur le système initialement immobile : une simulation de référence, sur la base d'un Hamiltonien standard et trois simulations adaptatives, c'est-à-dire utilisant un Hamiltonien adaptatif relatif d'un procédé selon l'invention tel que décrit ci-dessus (pas de temps de taille 0,0488 femtosecondes (fs), 7000 pas de temps, temps de simulation total égal à 342 fs). Pour chacune des simulations utilisant un Hamiltonien adaptatif relatif, la racine carrée de la fluctuation par rapport à la simulation standard, nommée RMSD, est donnée, de même que l'erreur de déplacement de particule maximum Aq _max .

, où q, est le vecteur des coordonnées de la particule a, au dernier pas de la simulation adaptative, et qf est le vecteur des coordonnées de cette même particule au dernier pas de la simulation de référence.

Par exemple, pour la simulation adaptative où £ ^r =0,0l kcal/mol et f =1 ,01 kcal/mol ( <5 =1 kcal/mol), on obtient un facteur d'accélération du temps de calcul nécessaire pour mettre en œuvre cette simulation adaptative par rapport au temps de calcul relatif à la simulation de référence, d'une valeur égale à 1 ,5, RMSD=8,6 S et Aq _max =89,4 S, où S est la distance d'équilibre dans le potentiel de Lennard-Jones utilisé.

Pour la simulation adaptative où e ^r =\ kcal/mol et e ^f =2 kcal/mol ( <5 =1 kcal/mol), on obtient un facteur d'accélération du temps de calcul nécessaire pour mettre en œuvre cette simulation adaptative par rapport au temps de calcul relatif à la simulation de référence, d'une valeur égale à 1 ,55, et RMSD=8,7 S et Aq _max =89,4 S.

Pour la simulation adaptative où £ ^r =3 kcal/mol et e ^f =A kcal/mol ( <5 =1 kcal/mol), on obtient un facteur d'accélération du temps de calcul nécessaire pour mettre en œuvre cette simulation adaptative par rapport au temps de calcul relatif à la simulation de référence, d'une valeur égale à 1 ,6, et RMSD=12,3 S et Aq _max =89,4 S.

Ainsi un procédé selon l'invention permet d'accélérer les calculs, avec une altération potentiellement faible des comportements.

Les résultats dépendent de la représentation arborescente choisie pour le système de particules. Pour les simulations adaptatives ci-dessus, l'arbre binaire a été construit de bas en haut en divisant le système en moitiés.

Dans les modes de réalisation considérés de l'invention, on a considéré la fusion de sous-ensembles deux à deux en un ensemble rigide, dans d'autres modes de réalisation, le nombre de sous-ensembles fusionnés en un ensemble rigide est pris supérieur à deux.

La région de transition entre les deux états rigide et libre peut être prise de largeur plus ou moins grande.