FR2995109A1 - Procede de simulation d'un ensemble d'elements, programme d'ordinateur associe - Google Patents

Abstract

Procédé de simulation d'un système d'éléments représenté par un arbre comportant des noeuds feuilles représentant chacun un élément, et des noeuds internes dont un noeud racine R, sur la base d'un Hamiltonien avec p vecteur des moments, q vecteur des positions des éléments, et V étant l'énergie potentielle du système: lorsque des conditions déterminées sont vérifiées, on impose un même mouvement en translation à au moins les éléments descendant d'un noeud donné de l'arbre en définissant la matrice M égale à Phi , étant donné la formule récursive selon laquelle, pour tout noeud A de l'arbre comportant k noeuds fils

Description

Procédé de simulation d'un ensemble d'éléments, programme d'ordinateur associé 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 = -1pT.M-1 p +V , p étant un vecteur indiquant les moments des éléments, V étant 2 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. VVO 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 noauds feuilles, lesdits noauds feuilles représentant chacun un élément respectif, et des noauds internes dont un noaud 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 noaud donné de l'arbre en définissant la matrice AJ-1 comme étant égale à (DR où R est le noaud racine de l'arbre, étant donné la formule récursive selon laquelle, pour tout noaud A de l'arbre comportant noauds fils , 142 , (DA = PA E+(1- pA) mA (I)A 0 0 0 0 (1)4_ la matrice E étant une matrice de taille dnA*dnA formée de nA*nA blocs de taille d*d égaux à la matrice identité de dimension d, nA étant égal au nombre d'éléments descendant du noaud A, d étant la dimension de l'espace dans lequel évoluent les particules, mA est la somme de la masse de ces nA éléments descendant du noaud A, pA étant une fonction de restriction du noaud A comprise entre 0 et 1, dont la valeur est égale à 1 lorsque A est un noaud feuille ou lorsqu'un même mouvement a été imposé aux éléments descendants du noaud donné ; dans le cas d'un noaud feuille A,, la matrice (DA, est égale à la masse inverse de la particule représentée par le noaud A, 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 noaud 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 noaud donné de 2 .\ , où P AS 112 1-1 MA, l'arbre en fonction de la valeur prise par E A = C ( est la somme de la masse des éléments du noaud A, , et pAs, est le vecteur-somme des moments des éléments du noaud fils A,, 11.11 est la norme du vecteur, C est une constante positive ; - pour le pas de simulation courant pour un noaud interne A: - si EA est inférieur à un premier seuil, on fixe pA égal à 1 pour imposer un même mouvement en translation aux éléments descendants du noaud A ; - si EA est supérieur à un second seuil supérieur au premier seuil, on fixe pA égal à 0 ; - pour le pas de simulation courant, si EA est compris entre le premier seuil et le second seuil, on fixe pA égal à la valeur prise par une fonction d'interpolation d'ordre 5 de la variable EA ; - 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 oauvre 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 : système à 2 - 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 particules avec le Hamiltonien adaptatif relatif HAR selon l'invention ; - la figure 3 représente un dispositif mettant en oeuvre 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 : 1 II(P,q)=-2PT M-119 +V(q) 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 lo p 1,z, et la position q, de chaque particule a' i=1 à N, s'écrit (a a a 1 , ,i,z, - Pl,x Pl,y Pl,z Les vecteurs p et q s'écrivent donc : p= et q = PN,x PN,y _PN,z _ 25 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, , pour i=1 à N, où m, est la masse de la particule ai. 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 HAR, ainsi : H 1 AR(P,q)=-2PT .(1)(P, q) .P ± V (q) (formule 1) dans lequel (I)(p,q) , 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 (I)(p, q) sont déterminées à chaque pas de simulation, la fonction (I) 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 HAR, on déduit des équations adaptatives de mouvement, définissant p et 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 : Ô11 ÔT7 1 7' ae'(P q) p aq aq 2 p P, (formules 2) 1 T ad)(1 (1)4 P)P ± 2 P P at ap 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 A1, ..., An, un ensemble de corps rigides indexés, on peut organiser les indices en un arbre binaire, de sorte que A1 ne peut fusionner qu'avec A2 en un corps rigide (A1FA2), quand une condition spécifique de fusion est satisfaite, A3 ne peut fusionner qu'avec A4 etc., et à un niveau supérieur dans l'arbre binaire, (A1FA2) ne peut fusionner qu'avec (A3+A4) etc. Un tel arbre comporte au rang hiérarchique maximum, un noeud racine, nommé noeud R, et au rang hiérarchique minimum, les noeuds feuilles correspondant aux particules considérées seules. Tous les noeuds de l'arbre, sauf les noeuds feuilles, s'appellent les noeuds 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 A1 ne peut fusionner qu'avec A2, , Ar, en même temps, etc. Une autre solution est de considérer le graphe déduit d'une recherche des voisins, les noeuds du graphe étant des particules ou des corps rigides, et une arête entre deux noeuds indique que la distance entre les deux noeuds 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 noeuds internes dans la hiérarchie, et des mises à jour incrémentales, par exemple pour les forces et énergies partielles. Dans un tel arbre, chaque noeud feuille représente une particule seule de l'ensemble, et chaque noeud interne représente un sous-ensemble de particules composé des particules représentées par les noeuds fils dudit noeud interne. On nommera corps C l'ensemble comprenant le noeud C et les descendants, directs ou indirects, du noeud C. Selon l'invention, si un corps C, i.e. un noeud C interne de l'arbre, est composé de deux corps A et B (i.e. le noeud C comporte deux noeuds fils A et B), on définit la matrice de masse inverse adaptative Oc correspondant à ce noeud C avec la formule récursive suivante : (DA 0 0 (DB (formule 3) où : la matrice de masse inverse adaptative (PA correspondant au noeud A et la matrice de masse inverse adaptative OB correspondant au noeud B, sont elles- mêmes définies selon la formule 3 si les noeuds A et B ne sont pas des noeuds feuilles ; dans le cas d'un noeud feuille A, respectivement B, la matrice (I)A, respectivement OB, est égale à la masse inverse de la particule représentée par le noeud A, respectivement le noeud B, multipliée par la matrice identité de dimension 3; - la matrice E est une matrice de taille 3nc*3nc, formée de rIc*nc blocs 3*3 égaux à la matrice identité de dimension 3, nc étant égal au nombre de noeuds feuilles du noeud C (i.e. l'ensemble des noeuds feuilles dans les noeuds fils directs du noeud C ou dans la descendance de ces noeuds fils) ; - mc est la somme de la masse des particules représentées par ces nc noeuds feuilles du noeud C; - pc 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 noeuds rigides A et B) sur la base de pc , la fonction de restriction du corps C.

Plus exactement, un corps C ou noaud interne C, est considéré comme rigide (i.e. toutes les particules feuilles comprises dans le noaud C ou les descendants du noaud C, sont animées du même mouvement en translation) lorsque la fonction pc prend la valeur 1, et sinon il est considéré comme non-rigide.

Tous les noauds 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 pc a pour valeur 0, le noaud interne C est considéré comme libre. Dans le mode de réalisation de l'invention, la fonction pc est définie récursivement, de façon à lisser le basculement entre les valeurs 0 et 1, selon la formule suivante : Pc =I1-1(ec)PAPB dans le cas contraire, 1 si C est un noeud feuille, (formule 4) où ,u(s) G [0,11 est une fonction deux fois dérivable de E : (formule 5) Si 0 E , 1-1(e) 0 Si E > E f , 4E) G [0,11 autrement, avec les seuils Ef et Er définis de sorte que Ef > Er et s étant une fonction deux fois dérivable de E .

La nature deux fois dérivable de ji permet de préserver la stabilité de la simulation de l'ensemble E de particules. Par exemple, une forme possible pour s(E) est une fonction d'interpolation d'ordre s - sr 5, par exemple égale à - 6/75 +15/74 -10/73 +1 , avec = où = - Er 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, Ec est choisi dépendant d'une relation du moment des deux corps A et B, avec ec (Pc ) = 1 11P,S411 1 VpBS 112 11pAS pBS 112 11pASniB pBSrnA 112 (formule 6) 2 2 mA 2 mB 2 mc 2 mAmBnic où A et B sont les noauds fils du noaud C, où le vecteur pcs est la somme des moments de toutes les particules qui sont des noauds feuilles dans le corps correspondant au noaud C : pcs , et est la somme i eC s sss des vecteurs pA et pB : pc=pA+pB. Cette invention peut être généralisée pour un arbre k-aire. Dans ce cas, pour le noaud parent avec k noauds fils AI, la matrice inverse de l'inertie est (I) A = PA E+(1- pA) A ( ) ° ° 0 - . 0 0 0 D Ai, 2 s 1 k 11PAs 112 1 2 et E A (pA)= ' 1 k P 112 1011 2 1=1 niA, 2 2 mA, 2 mA Le coefficient 1/2 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 Er 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 pc différente pour des noauds internes différents. Par exemple, jusqu'à certain niveau de l'arbre pc est définie selon les formules 4 et 5, donc, les noauds peuvent être réunis ensemble et séparés, et dans les niveaux de l'arbre supérieurs à ce niveau, pc 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 pc devront être fixées à 0 pour tous les noauds internes de ce sous-arbre et ne jamais être mises a jour. Dans tous les cas, les fonctions pc doivent être définies pour tous les noauds 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 cl) utilisée dans la formule (1) du Hamiltonien adaptatif relatif HAR, égale à cl)R, qui est la matrice définie conformément à la formule récursive 3 relativement au noeud racine R de l'arbre représentant l'ensemble E de particules : (1)(p, q)=c1)(p)=(DR On a alors le Hamiltonien HAR, tel queI - 1(p,q)=-1pT.e^R(p).p+V(q), qui est 2 séparable car la matrice (1)R 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 : .p T3HAR V &:1 (formules 7) ae^R(p) p. * HAR =a)R(p)p1 T + 2 P at ap Il ressort en outre clairement de la formule 3 que si le noeud C est libre (pc =0), alors la matrice cl) c est composée de deux blocs sur la diagonale correspondant l'un au noeud fils A et l'autre au noeud fils B.

Dans un premier exemple de simulation selon l'invention, il est considéré un système à une dimension comprenant deux particules Pl, 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 noeud racine C dont les deux noeuds 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 HAR selon l'invention sont représentées en figure 2 en fonction du temps de simulation, pour différentes valeurs prises par ê et sf . 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 noeud racine C devient alors rigide. En conséquence, dans un cas général, les noeuds fils du noeud 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'a un pas donné de la simulation, la condition de rigidité du noeud 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 sf 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 sf , la région de rigidification est plus longue et plus plate. Ainsi selon l'invention, la matrice t), définie égale à (I)R, 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 oeuvre 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 HAR 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 oeuvre l'étape préliminaire 100 et itérativement les étapes 101 à 104 lors d'une n+1Ième itération du programme P correspondant à l'instant de calcul hi,+1= ho + (n+1)h, avec n entier kO, 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 noauds feuilles sont alors rigides. Dans une étape 101, les forces d'interaction s'exerçant entre les particules à la n+1Ième itération, sont déterminées.

Dans une étape 102, les moments des particules à la n+1Ième itération, sont déterminés. Dans une étape 103, les métriques relatives aux noauds internes de l'arbre à la n+1Ième itération, sont déterminées. Elles comprennent notamment pc et Ec Dans une étape 104, les positions des particules à la n+1Ième 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 oauvre en fonction des valeurs calculées lors de l'itération précédente. La mise en oauvre 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 Ec conformément à la formule 6) à l'itération de simulation précédente (les forces sont déterminées avant Ec). 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 noaud 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 noaud interne C comprend des vecteurs 3D, dont le nombre est égal au nombre de noauds feuilles descendant directement ou indirectement du noaud C. Chaque vecteur de cette table est la force d'interaction agissant sur la particule représentée par un de ces noauds feuilles en provenance de toutes les autres particules représentées par les autres noeuds feuilles descendant de C. Chaque vecteur de la table de forces partielles du noeud racine R correspond à la force globale agissant sur une particule. Une table de forces partielles correspondant a un noeud 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 noeuds feuilles vers la racine). Pour chaque noeud 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 noeud feuille. La boîte de chaque noeud interne est recalculée par exemple en englobant les boîtes de ses deux noeuds 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 noeuds internes de l'arbre. Pour chaque noeud C, ayant deux noeuds 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 noeud feuille descendant directement ou indirectement du noeud A), tandis que l'autre appartient au corps B. Cette liste d'interaction peut être construite pour chaque noeud interne comme décrit pour d'autres types de volumes englobants dans l'article R. Rossi, M. lsorce, S. Morin, J. Flocard, K. Arumugam, S. Crouzy, M. Vivaudou, and S. Redon. "Adaptive torsion-angle quasi-statics: a general 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 noeuds internes. Pour toutes les autres itérations, dans les noeuds 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 noeud est identifié comme noeud rigide, les forces entre les particules du noeud n'ont pas changé. Par conséquent, seules les tables de forces partielles des noeuds non-rigides doivent être actualisées. Elles peuvent être actualisées par exemple comme suit récursivement, en partant des noeuds feuilles vers le noeud racine. Pour chaque noeud C non-rigide, dans la première moitié de la table de forces partielles relative au noeud C, les éléments de la table de forces partielles relative au noeud A sont copiés et dans la seconde moitié de la table de forces partielles relative au noeud C, les éléments de la table de forces partielles relative au noeud B sont copiés. Puis, pour chaque paire de la liste d'interaction établie à la sous-étape 101_2 pour le noeud 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 noeud 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 noeud C dès lors que la liste d'interactions correspondant audit noeud 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 pi,'+1 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 : 19,,, fi,'111, h étant le pas de temps de la simulation et f,'i étant la force d'interaction globale s'exerçant sur la particule 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 noeud 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 (Di a R(p) des particules à partir de l'équation = e^R(p)p+-pT p, est 0(N3) en raison du 2 (p) terme R . Ce terme est une dérivée d'une matrice par rapport à un vecteur, donc un objet à trois dimensions, où, pour i= 1 à N : ad)," a.:1)1N al) al) ad.'7, api Dans le mode de réalisation considéré, la complexité peut être réduite à 0(N.logN) par une mise à jour récursive depuis les noeuds feuilles de structures de données Qc, Rc et Sc pour chaque noeud C interne avec deux noeuds fils A et B.

En remarquant que Pc = PA _PB_ ou IDA est le vecteur des moments de toutes les particules qui sont des noeuds feuilles du noeud A et pB est le vecteur des moments de toutes les particules qui sont des noeuds feuilles du noeud B, le terme Qc=e,c(pc)pc correspondant au premier terme de l'équation donnant : PcU9c Qc = E±(1-Pc(Pcs)) Pc(P) y bec c P,}±(1-Pc(Pcs» c ( A(PA) (I) B(PB)_ (I)A(PA)PA _( )B(PB)PB _ Pc = Pc(P) sl QA _QB _ (formule 9) 1Pc f+ -Pc(Pcs)) mc où {x} est le vecteur de dimension égale à n et dont chacune des n composantes est égale à l'élément x. Pc(P) . De même, en considérant pour le noeud C, le terme Rc, = aPc Rc avec Pc =,u(sc)pApB, en différenciant, ( pB(pBs) du(ec)apA(pAs) ' s )a,u(Ec)aEc 1.-ApA aEc aPA, ( B, .1, c, .1 P A(P As ) I-1(e B(P Bs + pB (ns )al-1(e c) E (--PB Rc= ( ( s s s all(EC) PA Pc pB(pBs) ,u(Ec)RA+PA(PA) aec 111A 111c )aPA, ( ( s s ap(Ec) PB Pc p,(p,$) ,u(Ec)RB+ pB(pB) aE m m c B C )_ (formule 10) = PB(PBs) I-1(ec)RA± PA(PAs)ajl(ec) pA(pAS) du(ec)RB ± ,BI ,S %all(EC 1 PB PC 1 \ c ) _ ( \ C ) P "B) aE niB IllC aE mA MC PA Pc S S Le dernier terme de l'équation donnant peut être mis à jour récursivement en utilisant la formule 9.

En effet : T a(1)C(PC) SC Pc Pc= 0 (DB_ (DA 0 Pc±(1-Pc(Pcs))PcT +0-pc(pcs)) a (DA 0 Pc= 0 (DB _)) 0 (1)B_ Pc= ))) (DA 0 0 (DB _) apc ( a 1 1 1 T a E -Pc Pc(Pcs) Pc m c 1 ( T a Pc Pc(Pcs) Pc ( (DA 11- 1c 0 = aPc(p'â( T( u a(DA(pA) Pc - pAT PA aPc 111c T a(DB(pB) PB PB (DA 0 0 (DB_ Pc + -pc(pcs)) ) ) ( QA _QB_ SA aPc(Pcs)(Pc 1 _ B_ 11)1- .111C ieC +0-pc(pcs)) ) ) ( Ils Pc- Pc c =Rc, QA +0-pc(pcs)) SA (formule 11) _ B_ _QB_)) La fonction 0 représente le produit scalaire.

Ainsi l'étape 103 comporte une mise à jour des structures Qc, Rc et Sc pour chaque noaud C comportant deux noauds fils A et B. On notera que l'étape 100 comporte en outre une étape d'initialisation de ces structures où pour chaque noaud interne C : pc=0 et les valeurs de Qc, Rc et Sc sont également nulles, et pour chaque noeud feuille F représentant une particule : pF=1; QF PF I IIIF S F - 0 RF - 0 . Dans l'étape 103, la mise à jour des métriques Qc, Rc et Sc, ainsi que de sc et pc relativement au noeud C est effectuée récursivement, le long de la hiérarchie des noeuds depuis la racine: - Si le noeud C est interne, avec deux fils A et B: - faire la mise à jour des métriques du noeud A; - faire la mise à jour des métriques du noeud B; - mettre à jour sc à l'aide de la formule 6; - mettre à jour pc : pc ,u(sc)pApB ; - mettre à jour Qc à l'aide de la formule 9; - mettre à jour Rc à l'aide de la formule 10; - mettre à jour Sc à l'aide de la formule 11; - Si le noeud C est un noeud feuille représentant une particule aF QF = pF/mF Pour optimiser les calculs, la mise à jour des vecteurs Qc, Rc et Sc peut être effectuée seulement pour les noeuds pour lesquels pc> 0 parce que ces métriques peuvent ne pas être utilisées pour les noeuds 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 noeud 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 oeuvre de l'étape 104 peut être la suivante pour la mise à jour de la position d'un noeud C, chaque noeud interne C étant considéré comme comportant deux noeuds fils A et B : - si pc = 0 (le noeud C est un noeud libre) : - mettre à jour la position du noeud A; - mettre à jour la position du noeud B; - sinon, pour chaque noeud feuille k descendant du noeud C, directement ou indirectement : q1,1 = ql,n+ (Qc[k] + 0.55c[k]) h, q1,'+1 étant la valeur de la position de la particule ab pour l'itération n+1 considérée et h le pas de temps de simulation, Qc[k] étant le kème élément du vecteur Qc et Sc[k] étant le kè" élément du vecteur Sc . 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 noeuds fils du noeud C doit être établie, par exemple pendant l'initialisation de la simulation, pour toutes les particules et tous les noeuds internes de l'arbre. Dans cette réalisation de mise à jour des positions, les métriques des noeuds pour lesquels pc> 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" 1 N peuvent aussi être déterminées comme suit : qi, (QR[k]+ 0.5S R[k])h où l'index R correspond au noeud racine. Pour chaque noeud interne C de l'arbre, 4 vecteurs de longueur égale au nombre de noeuds feuilles descendants directs ou indirects du noeud C, sont stockés : Qc, Rc et Sc et une table de forces partielles. La complexité en espace de mise à jour des positions est par conséquent 0(N.logN). La complexité en temps est aussi 0(N.logN) puisque ces structures doivent être mises à jour à chaque pas de temps, et que Qc est toujours mise à jour pour tous les noeuds, y compris les noeuds 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 1g/mol chacune, en utilisant un potentiel Lennard-Jones (E. /kB = 120 Kelvin, où En, est le minimum d'énergie, distance d'équilibre S = 3.4 àngstrôms, distance de coupure 8 àngstrôms, le potentiel étant tronqué de manière lisse entre 7.5 et 8 àngstrôms), 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 A.qmax. RMSD 112 , où q, est le vecteur des coordonnées de la particule a, au dernier pas de la simulation adaptative, et q{ 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ù s'=001 kcal/mol et El =1,01 kcal/mol (ê =1 kcal/mol), on obtient un facteur d'accélération du temps de calcul nécessaire pour mettre en oeuvre 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 Aimax=89,4 S, où S est la distance d'équilibre dans le potentiel de Lennard-Jones utilisé. Pour la simulation adaptative où Er =1 kcal/mol et s=2 kcal/mol (ê =1 kcal/mol), on obtient un facteur d'accélération du temps de calcul nécessaire pour mettre en oeuvre 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 Aimax=89,4 S. Pour la simulation adaptative où Er =3 kcal/mol et s=4 kcal/mol (ê =1 kcal/mol), on obtient un facteur d'accélération du temps de calcul nécessaire pour mettre en oeuvre 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 Aimax=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.

Claims (8)

1. REVENDICATIONS1. Procédé de simulation d'un système d'éléments, selon lequel le comportement desdits éléments est déterminé, à des pas de simulation successifs, sur la base d'un Hamiltonien H du système d'éléments, tel que H (p,q) -1p7. +y avec p étant un vecteur 2 indiquant les moments des éléments, q un vecteur indiquant les positions des éléments, M-1 étant une matrice fonction des masses des éléments et V étant l'énergie potentielle du système, caractérisé en ce que ledit procédé comporte les étapes suivantes : on représente le système d'éléments en un arbre k-aire comportant des noeuds feuilles, lesdits noeuds feuilles représentant chacun un élément respectif, et des noeuds internes dont un noeud 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 noeud donné de l'arbre en définissant la matrice Ar comme étant égale à OR où R est le noeud racine de l'arbre, étant donné la formule récursive selon laquelle, pour tout noeud A de l'arbre comportant k noeuds fils A1, A2, Ak, F- ) 0 (13 A = »221-E (1-PA) 0 -"- mA 0 0 CD A la matrice E étant une matrice de taille dnA*dnA formée de nA*nA blocs de taille d*d égaux à la matrice identité de dimension d, nA étant égal au nombre d'éléments descendant du noeud A, d étant la dimension de l'espace dans lequel évoluent les particules, mA est la somme de la masse de ces nA éléments descendant du noeud A, PA étant une fonction de restriction du noeud A comprise entre 0 et 1, dont la valeur est égale à 1 lorsque A est un noeud feuille ou lorsqu'un *même mouvement a été imposé aux éléments descendants du noeud donné ; dans le cas d'un noeud feuille A,, la matrice C A. est égale à la masse inverse de la particule représentée par le noeud A, multipliée par la matrice identité de dimension d.30
2. 2. Procédé de simulation d'un système d'éléments selon la revendication '1, selon lequel on impose, pour le pas de simulation courant, un même mouvement en translation à au moins les éléments descendants d'un nceud donné de l'arbre en fonction de la valeur d'une fonction des moments desdits éléments.
3. 3. Procédé de simulation d'un système d'éléments selon la revendication 1 ou 2, selon lequel on impose, pour le pas de simulation courant, un même mouvement en translation à au moins les éléments descendants d'un noeud donné de l'arbre en fonction de la valeur 2 prise par eA=C P A, , où m4 est la somme de la masse des vdc Zri=lin éléments du nceud A ' et 19;4, est le vecteur-somme des moments des éléments du nceud fils A, , In est la norme du vecteur, C est une constante positive.
4. 4. Procédé de simulation d'un système d'éléments selon la revendication 3, selon lequel, pour le pas de simulation courant pour un noeud interne A: - si C A est inférieur à un premier seuil, on fixe pA égal à 1 pour imposer un même mouvement en translation aux éléments descendants du noeud A - si eA est supérieur à un second seuil supérieur au premier seuil, on fixe pA égal à O.
5. 5. Procédé de simulation d'un système d'éléments selon la revendication 4, selon lequel, pour le pas de simulation courant, si eA est compris entre le premier seuil et le second seuil, on fixe pA égal à la valeur prise par une fonction d'interpolation d'ordre 5 de la variable eA
6. 6. Procédé de simulation d'un système d'éléments selon l'une quelconque des revendications précédentes, comprenant 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.
7. 7. Procédé de simulation d'un système d'éléments selon la revendication 6, selon lequel l'information relative audit élément comprend l'énergie potentielle dudit élément et/ou la force d'interaction appliquée audit élément.
8. 8. Programme d'ordinateur (P) de simulation d'un système d'éléments, comprenant des instructions logicielles pour mettre en oeuvre les étapes d'un procédé selon l'une des revendications 1 à 7 lors d'une exécution du programme par des moyens de calcul.