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é 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 ^~l .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 dans des cas est par exemple égale au [ou est fonction du] potentiel d'interactions V(q) entre les éléments dont les forces d'interactions peuvent être dérivés, q étant un vecteur indiquant les positions des éléments (dans un cas plus général le potentiel d'interactions peut être dépendant également des moments des éléments), La simulation d'un ensemble d'éléments permet d'étudier le comportement d'un tel ensemble et d'analyser ses propriétés : les déplacements en termes de positions successives et des moments 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, et d'en dériver des équations de mouvement.

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. Des méthodes ont été proposées pour accélérer les simulations d'un ensemble d'éléments.

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é, mis en œuvre par ordinateur et caractérisé en ce que ledit procédé comporte une étape selon laquelle, lorsque le vecteur moment p prend certaines valeurs déterminées relatives à au moins un élément, on affecte une valeur nulle à au moins un terme diagonal de la matrice M ^"1 relatif audit élément.

L'invention permet de réduire le volume et, par conséquent, le temps de calcul, requis pour la détermination de l'énergie potentielle, du potentiel d'interaction, des forces d'interaction, des positions et/ou des 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 :

- ledit procédé comporte une étape selon laquelle, pour au moins un desdits éléments, si un paramètre représentatif de l'énergie cinétique dudit élément a une valeur inférieure à un premier seuil strictement positif, on affecte une valeur nulle à au moins un terme diagonal de la matrice M ^"1 relatif audit élément ;

- les termes diagonaux de la matrice M ^"1 qui sont fonction de la masse d'un élément sont affectés à une valeur maximale lorsque l'énergie cinétique dudit élément est supérieure à un second seuil strictement positif ;

- on affecte une valeur nulle à au moins un terme diagonal de la matrice M ^"1 relatif audit élément si le couple comprenant le moment de l'élément et la position de l'élément prend certaines valeurs déterminées ;

- il comprend 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 à un k-uplet d'éléments, avec k entier supérieur ou égal à 2, pour lesquels une valeur nulle a été affectée aux termes diagonaux de la matrice M ^"1 à un instant de simulation précédent, sont par conséquent inchangées entre au moins ledit instant de simulation précédent et l'instant de simulation courant, et calculant une valeur de l'information relative à un élément donné, à un instant de simulation courant en mettant en œuvre les étapes suivantes lorsque des valeurs nulles n'ont pas été affectées aux termes diagonaux de la matrice concernant chaque élément d'un k-uplet d'éléments dont fait partie ledit élément donné :

- calculer une valeur de travail de ladite information relative audit élément donné en retranchant à la valeur de l'information relative audit élément donnée et déterminée à l'instant de simulation précédent, au moins les valeurs de l'information relative audit élément donnée et associée auxdits k-uplets d'éléments à l'instant de simulation précédent, et/ou - ajouter à ladite valeur de travail au moins les valeurs de l'information relatives audit élément donné et associées aux k-uplets d'éléments, déterminées à l'instant de simulation courant ;

- à un instant de calcul courant, une liste courante des paires d'éléments séparés par une distance inférieure à un seuil donné est dressée à un instant de simulation courant, et comparée à une liste précédente des paires d'éléments séparés par une distance inférieure à un seuil donné dressée à un instant de simulation précédent, et la valeur d'une information relative à un élément donné, à un instant de simulation courant, est calculée sur la base des paires comportant ledit élément donné en mettant en œuvre les étapes suivantes:

- calculer une valeur de travail en retranchant, à la valeur d'information relative audit élément donné et déterminée à l'instant de simulation précédent, les valeurs d'information relative audit élément donné associées à l'autre élément des paires considérées si la paire considérée est présente seulement dans la liste précédente ou si le vecteur reliant ledit élément donné à l'autre élément de la paire a varié entre l'instant de simulation précédent et l'instant de simulation courant ;

- la valeur de l'information relative audit élément donné à un instant de simulation courant est déterminée en ajoutant, à la valeur de travail, les valeurs de l'information relatives audit élément donné et associées à l'autre élément des paires considérées si la paire considérée est présente seulement dans la liste courante ou si le vecteur reliant ledit élément à l'autre élément de la paire a varié entre l'instant de simulation précédent et l'instant de simulation courant ;

- à un instant de calcul courant, une liste courante de k-uplets d'éléments satisfaisant certaines conditions, avec k entier supérieur ou égal à deux, est dressée à un instant de simulation courant, et comparée à une liste précédente de k-uplets d'éléments satisfaisant lesdites conditions à un instant de simulation précédent,

et la valeur d'une information relative à un élément à un instant de simulation courant, est calculée sur la base des k-uplets comportant ledit élément en mettant en œuvre les étapes suivantes:

- calculer une valeur temporaire en retranchant à la valeur de l'information relative audit élément et déterminée à l'instant de simulation précédent, les valeurs des informations relatives audit élément et associées auxdits k-uplets à l'instant de simulation précédent, lorsque lesdits k-uplets sont présents seulement dans la liste précédente ou lorsque les valeurs des informations relatives audit élément et associées auxdits k-uplets ont changé entre l'instant de simulation précédent et l'instant de simulation courant;

- déterminer la valeur de l'information relative audit élément à l'instant de simulation courant en ajoutant, à la valeur temporaire, les valeurs des informations relative audit élément et associées auxdits k-uplets à l'instant de simulation courant, lorsque lesdits k-uplets sont présents seulement dans la liste courante ou lorsque les informations associées auxdits k-uplets ont changé entre l'instant de simulation précédent et l'instant de simulation courant (par exemple lorsque des positions relatives des k éléments dans le k-uplet ont changé) ;

- l'espace de localisation des éléments est partitionné en cellules et chaque élément, à chacun d'instant de simulation précédent et un instant de simulation courant, est associé à une cellule d'appartenance selon des coordonnées de position déterminées audit instant de simulation, et selon lequel, pour les premiers éléments tels que les termes de la matrice M ^"1 relatifs audits premiers éléments n'ont pas été affectés à une valeur nulle à un instant de simulation courant, les étapes suivantes sont mises en œuvre :

- on détermine la cellule d'appartenance des premiers éléments à l'instant de simulation précédent ;

- pour chaque premier élément, on détermine dans ladite cellule d'appartenance ou ses cellules dans un voisinage donné de la cellule d'appartenance, les seconds éléments situés à l'instant de simulation précédent à une distance inférieure à un seuil donné dudit premier élément ; on calcule une valeur de travail en retranchant de la valeur d'une information relative audit premier élément et déterminée à l'instant de simulation précédent, les valeurs de ladite information relatives audit premier élément et associées auxdits seconds éléments ;

- on détermine la nouvelle cellule d'appartenance des premiers éléments à l'instant de simulation courant ;

- pour chaque premier élément, on détermine dans la nouvelle cellule d'appartenance ou les cellules dans le voisinage donné de la nouvelle cellule d'appartenance, les troisièmes éléments situés à l'instant de simulation courant à une distance inférieure à un seuil donné dudit premier élément ;

on détermine la valeur de l'information relative audit premier élément, à l'instant de simulation courant, en ajoutant à la valeur de travail, les valeurs de l'information relatives au premier élément et associées aux troisièmes éléments ; - l'information relative audit élément comprend l'énergie potentielle dudit élément et/ou la force d'interaction appliquée audit élément ;

- il comprend, à certains instants de simulation, une étape de détermination d'une information I, ladite étape tirant avantageusement parti du fait qu'une valeur nulle a été affectée à certains termes diagonaux de la matrice M ^"1 à certains instants de simulation ;

- il comprend, à certains instants de simulation, une étape de détermination d'une information I, ladite étape tirant avantageusement parti du fait que cette information I est inchangée et ne nécessite pas d'être déterminée de nouveau lorsqu'elle a été déterminée à un instant de simulation antérieur et qu'une valeur nulle a été affectée à un ensemble correspondant de termes diagonaux de la matrice M ^"1 entre au moins le dit instant de simulation antérieur et l'instant de simulation courant (on entend par « ensemble correspondant des termes diagonaux » les termes diagonaux qui ont une influence sur la valeur de l'information I, ie l'information I ne change pas quand ces termes sont nuls) ;

- il comprend, à certains instants de simulation, une étape de détermination de l'énergie potentielle, respectivement des forces d'interaction, ladite étape tirant avantageusement parti du fait qu'une valeur nulle a été affectée à au moins un élément diagonal de la matrice M ^"1 à certains instants de simulation.

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é selon l'une des revendications 1 à 12 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 un dispositif mettant en œuvre un mode de réalisation de l'invention ;

- la figure 2 représente en ordonnée l'évolution de la fonction /λ (£ _; ) en fonction de ε _i figurant en abscisse ;

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

- la figure 4 illustre un mode de réalisation de l'étape 103 ;

- la figure 5 illustre un autre mode de réalisation de l'étape 103 ;

- la figure 6 représente des simulations de trajectoire d'une particule relié à un point fixe, dans l'espace de phase (p,q), à Hamiltonien constant.

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 :

H(p, q) = ^ p ^T .M - ¹ .p + V(q) , p étant un vecteur indiquant le moment des particules, q un vecteur indiquant la position des 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 leur position 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 _{i z} ) et la position q, de chaque particule a,, i=1 à N, s'écrit (q _{i x} , q _{i y} , q _{i z} ).

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

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 H _A , ainsi :

1 T

H _A (p, q) =— p .<î>(p, q) .p + V(q) , dans lequel Φ(/?, <?) , matrice diagonale

3N*3N dite matrice inverse de masse adaptative, remplace M ^"1 et dépend du vecteur p, et éventuellement du vecteur q.

De cet Hamiltonien adaptatif H _A , 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 d'abord 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 : __aH __ay _j_ _τ fà(p,g)

dt dq dq 2 dq

formules (1)

. dq dH _A 1 _τ ΒΦ(ρ,ρ)

at dp 2 dp

Selon l'invention, les termes de la matrice Φ relatifs à la particule a,, pour i=1 à N, sont les suivants :

Φ[3ί-2, 3i-2](pi,qi) = Φ[3ί-1, 3i-1](Pi,qi) =Φ[3ί, 3i](p _i3 C|i) = — (l-p _t (q _t ,p _t )), où m _i est la masse de la particule a,.

On nomme _ί (ρ _ί ,ς _ί ) la valeur— (1- > _; (<? _; ,/? _; )) .

m _i

Selon l'invention, /λ(<? _; ,/? _; )€ [θ,ΐ] et est une fonction deux fois dérivable, qui est utilisée pour que la position de la particule soit constante pendant un certain temps.

Lorsque A (<? _;> /> _; ) = 0 (i.e. il n'y a pas de fixation de la position de la particule a,),

Φ _; (Pi'Qi ) =— - ^et ' ^a particule a, suit les lois dynamiques usuelles correspondant à m _i

l'Hamiltonien H standard du système E.

Lorsque p^q^p^ = 1 (i.e. la position de la particule est figée), Φ _; ) = 0 et la particule a, ne bouge pas, quelles que soient les forces qui s'appliquent sur elle (sa masse est considérée comme infinie).

Lorsque ]0,l[, la particule a, effectue une transition lisse entre ces deux comportements.

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

Dans un mode de réalisation de l'invention, on définit que : 1 « 0 < £,· <

i ₍ (£ ₍ ) e [0,l] « e ₍ e

où s(£ _; ) est une fonction de p, et éventuellement de q,, et est deux fois dérivable.

En figure 2, les valeurs prises par /λ (£ _; ) dans un mode de réalisation sont représentées en ordonnée, en fonction des valeurs prises par e _i figurant en abscisse, avec ε =1 et e{ =1 .1 .

Par exemple, une forme possible pour s _t (e _i ) est - 6η ⁵ + 15η ⁴ - 10η ³ + 1 , avec η = ^L et δ = ef - e. ^r .

Dans le mode de réalisation considéré plus loin, la fonction e relative à la particule a,, est choisie égale à l'énergie cinétique de la particule a,, soit e _t =

2m _i

L'invention consiste donc à « figer » les particules, en leur affectant une « pseudo »-masse infinie, lorsque leur énergie cinétique passe sous une certaine valeur, la quantité de mouvement de ces particules n'étant, elle, pas figée.

La fonction p _i est une fonction comportant comme variable le moment (dans le cas articulier considéré en exemple, elle n'est pas donc dépendante de la position

Dans d'autres modes de réalisation, la fonction p _i peut être une fonction, dépendant du moment de la particule a, (et éventuellement de sa position) autre que l'énergie cinétique bien sûr.

Dans des modes de réalisation, on fige une particule dont le moment (ou le couple comprenant le moment et la position) prend des valeurs prédéterminées (valeurs discrètes ou plages de valeurs).

Les équations adaptatives de mouvement (1 ) deviennent ainsi : dq dq

formules (2) dp 2 dp

où p(p) est une matrice diagonale 3N*3N indiquant les p^q^ p^ , pour i= 1 à N : p(p) [3i-2,3i-2]= p(p) [3i- l,3i- l]= P(p) [3i,3i]= , pour i= 1 à N. Comme indiqué auparavant, lorsque Pi iq^ Pi ) = 0 (ie pas de fixation de position),

Φ; (/?; _> <?; ) ^{= ~ ~ et} ' ^a particule a, suit les lois dynamiques usuelles correspondant à m _i

l'Hamiltonien H standard du système E.

Lorsque A (<?; , .?; ) = 1 (ie position de la particule constante figée), Φ _; {p _i ,q _i ) = 0 , par conséquent, q a une valeur nulle (en effet, p _i étant alors constant égal à 1 , la valeur du terme dA ( - ) _{est nu} || _e ) Dans une interprétation, la masse est considérée comme dp,

infinie, la particule a, est considérée comme figée.

Lorsque p _i q _i , p _i ) e ]0,l[ , la particule a, effectue une transition lisse entre ces deux comportements.

Ainsi selon l'invention, la matrice <ï> (p,q) spécifie comment, et quand, des degrés de liberté en position d'une ou plusieurs particules sont activés ou désactivés pendant la simulation.

Pour donner un exemple qui explique le comportement du système dans le cas de l'invention, en figure 6 sont représentées des simulations de trajectoire à une dimension dans l'espace de phase (p,q), d'un système comprenant une particule de masse 1 , attachée à un ressort de raideur 1 à un point fixe. Dans ce cas, N=1 . Les isolignes du

Hamiltonien adaptatif sont représentées, pour ε[ =0,5 et ε[ =0,8.

Notamment les courbes C1 , C2, C3, C4 correspondent chacune à une valeur constante respective du Hamiltonien adaptatif. Par exemple, la courbe C1 correspondant à un Hamiltonien égal à 1 . Le cercle D correspondant à une valeur constante égale à 1 d'un Hamiltonien standard, ie non-adaptatif.

La zone de l'espace des phases où la particule est figée se trouve entre les lignes en pointillés B2 et B3 (elle correspond à une valeur de moment comprises dans [-1 , 1 ]).

La zone de l'espace des phases où la particule est libre et suit une trajectoire conforme à un Hamiltonien standard se trouve au-dessus de la ligne B1 et en-dessous de la ligne B4.

La zone de l'espace des phases entre les lignes B1 et B2 et entre les lignes en pointillés B3 et B4 correspond à une zone de transition entre les états libre et figé de la particule.

Sur chaque isoligne C1 , C2, C3, C4, dans la zone où la particule est figée, la position q ne change pas, mais le moment p change. Pour mettre en œuvre une simulation du système E, une intégration numérique dans le temps de ces équations de mouvement indiquées en formules (2) est à réaliser. Dans le cas d'implémentation d'une simulation dans un ensemble NVE considérée ici à titre d'illustration, une méthode d'Euler partitionnée est par exemple utilisée (cf. par exemple « Géométrie numerical intégration : structure preserving algorithms for ordinary differential équations », Hairer E., Lubich C, Wanner G ; Volume 31 ; Springler Verlag 2006) pour cette intégration numérique.

Selon cette méthode d'Euler, les équations de la forme :

ù = a(u, v),

v = b(u,v),

résultent, pour l'intégration numérique, en l'ensemble suivant d'équations où h est un pas de temps :

Ainsi, selon cette méthode, les formules (2) peuvent s'écrire ainsi :

Formules (3) q _n+l = q _n + (M - ^l (l - p( _Pn+l )) _Pn+l ^~ _Ρη ^Τ _+ι Μ ^~ι ^ ^ p _n+l )h .

² op _n+l

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

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 (valeurs des coefficients de la matrice Φ, 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 celle illustrée en figure 6.

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 de type NVE. Le programme P comprend des instructions logicielles qui, lorsqu'elles sont exécutées sur le microprocesseur 3, sont adaptées pour opérer les étapes suivantes, en référence à la figure 3.

Dans une étape préalable 100, les fonctions /λ (/? _; , <? _; ) de la matrice Φ sont préalablement définies pour chaque particule a,.

P ²

Dans le cas présent, les fonctions définies ci-dessus _y o _i (- -^) ont ete choisies et

2m _;

définissent ainsi Φ(/?) βη fonction du vecteur des moments p, et des valeurs fixées pour ε\ et ε{ .

Des valeurs initiales sont également déterminées pour le moment p _{i 0} , la position q _{i 0} et la force d'interaction f _{i 0} de chaque particule a,, i= 1 à N, correspondant à un instant initial h ₀ .

Les étapes suivantes sont alors mises en œuvre itérativement 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.

On nommera ci-dessous :

f _{ij n+l} , la force d'interaction exercée par la particule a, sur la particule a, (qui est égale à - f _{ji n+l} ), à l'étape de calcul h _n+ i ;

f _{i n+i} , 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, considérée à l'instant de calcul h _n+1 ; elle est donc égale à p _{i n+l} , la valeur du moment de la particule a,, à Γ instant de calcul h _n+1 ; q _{i n+l} , la valeur de la position de la particule a,, à Γ instant de calcul h _n+ i ; p _{i n+l} , la valeur prise par la fonction p _i à Γ instant de calcul h _n+ i (comme vu ci- p ²

dessus ; elle est fonction de la valeur de l'énergie cinétique ) ;

Les étapes 101 , 101 b, 102, 103 sont destinées à la détermination des valeurs actualisées respectivement du moment, de la position et de la force d'interaction globale relative à chacune des particules a,. Dans une étape 101 , la valeur courante p _{i n+i} du moment de chaque particule a,, i=1 à N, est déterminée, d'après les formules (3), la valeur de la force f _{i n} d'interaction globale exercée sur la particule a, déterminée à l'instant de calcul précédent ayant été stockée en mémoire 2 comme celle du moment p _{i n}

Dans une étape 101 b, la valeur p _{i n+l} est recalculée à partir de la nouvelle valeur du moment p _n+1 : A ₊ i =

Dans une étape 102, la valeur courante q _{i n+l} de la position de chaque particule a,, i=1 à N, est à présent déterminée, d'après les formules (3) :

■

Dans une étape 103, la valeur courante f _{i n+l} de la force d'interaction globale, s'exerçant sur la particule a, et qui est due à toutes les autres particules au moins, i=1 à N, est déterminée, selon par exemple une des deux méthodes indiquées plus bas.

Le résultat de cette intégration numérique est la valeur actualisée des positions q _i et des moments p _i de chaque particule a _i , avec i=1 à N, déterminée pour l'instant de calcul h _n+1 .

La valeur actualisée d'autres paramètres caractéristiques du comportement des particules a, à l'instant h _n+1 , peut en outre être calculée, par exemple la valeur courante de l'énergie potentielle du système E, la valeur de l'autocorrélation entre les vitesses de particules.

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

Différentes techniques de calcul des valeurs actualisées des forces d'interaction dans l'étape 103, permettant de tirer avantageusement parti de la définition de la matrice Φ , peuvent être utilisées.

Une première technique comprend les étapes suivantes.

Lors de l'étape d'initialisation 100, une liste courante des paires de particules est dressée, telle que la distance entre les particules de chaque paire à l'initialisation est inférieure à un seuil dO (lorsque la distance entre deux particules est supérieure à dO, l'interaction entre ces deux particules est négligée) et la force d'interaction f _{ij 0} de la particule a, sur la particule a, de chaque paire présente dans la liste courante est en outre évaluée, en fonction de la distance les séparant et selon le champ de force simulé, et mémorisées.

A chaque paire dans la liste est associé un élément e _{ij Q} également mémorisé en mémoire 2, comportant l'identifiant de chacune des deux particules a,, a, de la paire, les coordonnées du vecteur r ₀ ^! ' joignant les deux particules et partant de la particule a,, et la force d'interaction f _{ij 0} exercée par la particule a, sur la particule a, (qui est égale à - f _{ji 0} , f _jifi étant la force d'interaction exercée par la particule a, sur la particule a, ).

La force d'interaction globale f _{i 0} s'exerçant sur la particule a, est égale à la somme des forces d'interaction f _{ij 0} exercée sur la particule a, par les particules a _j; j=1 à N.

Lors de chaque itération de l'étape 103, les étapes ci-dessous sont ensuite mises en œuvre, en référence à la figure 4.

Considérons que l'itération en cours soit la (n+1 ) ^eme .

Dans une étape 103_a1 , on affecte comme valeur de départ à la force d'interaction globale f _{i n+l} s'exerçant sur chaque particule a, la valeur de la force d'interaction globale f _{i n} calculée lors de l'itération précédente.

Dans une étape 103_b1 , une liste courante L _{a n+1} des paires de particules en interaction est dressée, i.e. ce sont les paires de particules telles que la distance entre les particules de chaque paire considérée à l'instant de calcul h _n+ i est inférieure au seuil dO.

A chaque paire dans la liste L _{a n+1} est associé un élément e _{ij n+l} comportant l'identifiant de chacune des deux particules a,, a, de la paire, les coordonnées du vecteur r ^ joignant les deux particules depuis la particule a,, conformément à leur localisation à l'instant h _n+1 , et la force d'interaction f _{ij n+l} exercée par la particule a, sur la particule a, (qui est égale à - _i „ ₊₁ , f _fii „ _+l étant la force d'interaction exercée par la particule a, sur la particule a,), dont la valeur à ce stade n'est pas encore renseignée.

Puis, dans une étape 103_c1 , la liste courante des paires L _{a n+1} est comparée à la liste précédente L _{a n} des paires, i.e. établie lors de l'itération précédente (soit la n ^eme itération). A chaque paire dans la liste L _a,n est associé un élément e _{ij n} comportant l'identifiant de chacune des deux particules a,, a, de la paire, les coordonnées du vecteur r _n ^{l }} joignant les deux particules et partant de la particule a,, calculé en fonction des positions déterminées lors de l'itération précédente pour les particules a,, a, et la valeur de la force d'interaction f _{ij n} exercée par la particule a, sur la particule a, (qui est égale à - f _{ji n} , f _{ji n} étant la force d'interaction exercée par la particule a, sur la particule a,).

Si une paire est présente seulement dans la liste précédente L _{a n} , cela signifie que l'interaction entre les deux particules de la paire a disparu entre l'itération n et l'itération n+1 .

Pour chaque paire a,, a, présente seulement dans la liste précédente, on retire de la force d'interaction globale en cours de détermination f _{i n+l} s'exerçant sur la particule a,, respectivement de f _{j n+l} s'exerçant sur la particule a _j; la force d'interaction f _{ij n} calculée à l'étape 100 lors de l'itération précédente, respectivement la force d'interaction f _{ji n} = -f _{ij n} .

Si une paire est présente seulement dans la liste courante L _{a n+} i , cela signifie que l'interaction entre les deux particules de la paire est apparue entre l'itération n et l'itération n+1 .

Pour chaque paire a,, a, présente seulement dans la liste courante, on calcule alors la force d'interaction f _{ij n+l} exercée par la particule a, sur la particule a,, en fonction de leur position respective notamment ; on la sauvegarde en mémoire dans l'élément e _{ij n+l} .

On ajoute la force d'interaction f _{ij n+l} à la force d'interaction globale sur la particule a en cours de détermination f _{i n+l} , et on ajoute la force d'interaction f _{ji n+l} = -f _{ij n+l} à la force d'interaction globale en cours de détermination f _{j n+l} s'exerçant sur la particule a,.

Si une paire est comprise à la fois dans la liste courante L _{a n+} i et dans la liste précédente L _{a n} , alors on compare les vecteurs r _n ^{l }} et joignant les deux particules a,, a,. S'ils sont différents, on calcule la force d'interaction f _{ij n+l} exercée par la particule a, sur la particule a,, en fonction de leur position respective notamment et on sauvegarde cette force d'interaction f _{ij n+l} dans l'élément e _{ij n+l} . On ajoute à la force d'interaction globale en cours de détermination f _{i n+l} s'exerçant sur la particule a,, la force d'interaction fi _j,n+ i■ On ajoute à f _{j n+l} s'exerçant sur la particule a,, la force d'interaction f _j i _,n+ i = ^~ fi _j,n+ i■ En outre, on retire de la force d'interaction globale en cours de détermination f _{i n+l} s'exerçant sur la particule a,, respectivement de f _{j n+i} s'exerçant sur la particule a,, la force d'interaction f _{ij n} calculée à l'étape 100 lors de l'itération précédente, respectivement la force d'interaction f _{ji n} = -f _{ij n} .

L'invention, en figeant des positions de particules, génère un nombre accru de paires pour lesquelles le vecteur entre deux particules, et donc la force d'interaction entre ces deux particules, demeurent inchangés.

Dans un tel cas, les méthodes proposées pour l'étape 103 permettent de ne pas recalculer toutes les composantes des forces d'interaction globale en tirant parti des caractéristiques d'un procédé selon l'invention.

Cette technique de détermination des valeurs des forces d'interaction globale est optimale en termes de volume de calcul. Toutefois, la construction des listes requiert du temps.

Une deuxième technique de réalisation de l'étape 103 permet d'exploiter les avantages conférés par l'invention sans utiliser de comparaison des listes de paires de particules en interaction à l'itération courante par rapport à l'itération précédente, mais en se servant d'une grille à trois dimensions (si le mouvement des particules est considéré dans en espace à trois dimensions ; si les particules se déplacent dans un plan, une grille à deux dimensions peut être suffisante).

Lors de l'étape d'initialisation 100, en outre, une grille initiale est créée, en considérant un parallélépipède comprenant toutes les particules et en le subdivisant en cellules, par exemple cubiques dont la taille d'un côté est supérieure ou égale à dO.

Puis chaque particule a,, i= 1 à N, est affectée à la cellule à laquelle elle appartient, selon la position de la particule à l'étape d'initialisation.

Puis pour chaque particule a, qui se trouve dans une cellule donnée, on considère la cellule donnée et/ou les cellules qui lui sont voisines (26 cellules considérées au maximum) dans lesquelles des particules a, sont à une distance de la particule a, inférieure à dO. Pour ces particules a, situées à une distance de la particule a, inférieure à dO et telle que j>i, on calcule la force d'interaction f _{ij 0} de la particule a, sur la particule a,.

Cette force est égale à - f _{ji Q} , f _{ji Q} étant la force d'interaction exercée par la particule a, sur la particule a,.

On notera que dans le mode de réalisation décrit, les cellules voisines considérées sont les cellules immédiatement voisines, i.e. celles qui ont au moins un côté commun avec la cellule donnée ; dans d'autres modes de réalisation, les cellules voisines considérées sont celles situées à r cellules de distances d'une cellule immédiatement voisine de la cellule donnée.

Les étapes suivantes sont ensuite mises en œuvre lors de chaque itération n+1 du programme P, avec n≥0, en référence à la figure 5.

Dans une étape 103_a2, on affecte comme valeur de départ à la force d'interaction globale f _{i n+l} s'exerçant sur chaque particule a,, la valeur de la force d'interaction globale calculée lors de l'itération précédente f _{i n} .

Dans une étape 103_b2, pour toutes les particules a, pour lesquelles p _i>n+1 < 1 (ie les particules non considérées comme figées), on détermine les particules a, vérifiant les conditions suivantes :

- ces particules a, étaient situées à l'itération précédente (n) dans la cellule dans laquelle la particule a, était positionnée à l'itération précédente n, ou les cellules qui lui sont voisines (26 cellules considérées au maximum) ; et

- ces particules a, étaient à l'itération précédente n à une distance de la particule a, inférieure à dO ; et

- ces particules a, vérifient que leur indice j est strictement supérieur à i ou vérifient que p _j>n+1 = 1.

La composition de la grille jusqu'ici considérée est donc celle correspondant aux positions actualisées à l'itération précédente (itération n).

Puis pour chacune de ces particules a, ainsi déterminées, on calcule la force d'interaction f _{ij n} exercée par la particule a, sur la particule a, (et par là-même la force d'interaction f _{ji n} exercée par la particule a, sur la particule a,), en fonction de la distance les séparant à l'étape précédente n.

Et on retranche cette force d'interaction f _{ij n} exercée par la particule a, sur la particule a, de la force d'interaction globale exercée sur a, ; de même on retranche cette force d'interaction f _{ji n} exercée par la particule a, sur la particule a, de la force d'interaction globale exercée sur a : ainsi on calcule f _{i n+y} = f _{i n+y} - f _{ij n} et f _{j >B+1} = f _{j n+y} - f f f

Dans une étape 103_c2, la composition de la grille est mise à jour, en déterminant les cellules courantes d'appartenance de toutes les particules a, pour lesquelles p _i>n+1 < 1 (ie les particules non considérées comme figées), selon les positions q _{i n+} i de ces particules correspondant à l'itération n+1 . Dans une étape 103_d2, pour toutes les particules a, pour lesquelles p _i>n+1 < 1 (ie les particules non considérées comme figées), on détermine les particules a, vérifiant les conditions suivantes :

- ces particules a, sont situées à l'itération courante (n+1 ) dans la cellule dans laquelle la particule a, est positionnée à l'itération courante, ou les cellules qui lui sont voisines (26 cellules considérées au maximum) ; et

- ces particules a, sont à l'itération courante (n+1 ) à une distance de la particule a, inférieure à dO ; et

- ces particules a, vérifient que leur indice j est strictement supérieur à i ou vérifient que p _j>n+1 = 1.

La composition de la grille considérée ici est donc celle correspondant aux positions actualisées à l'itération courante (itération n+1 ).

Puis, pour chacune de ces particules a, ainsi déterminées, on calcule la force d'interaction f _{ij n+l} exercée par la particule a, sur la particule a, (et par là-même la force d'interaction f _{ji n+l} exercée par la particule a, sur la particule a,), en fonction de la distance les séparant à l'étape courante n+1 .

Et on ajoute cette force d'interaction f _{ij n+l} exercée par la particule a, sur la particule a, de la force d'interaction globale exercée sur a, ; de même on ajoute cette force d'interaction f _{ji n+l} exercée par la particule a, sur la particule a, de la force d'interaction globale exercée sur ¾ : ainsi on calcule f _{i n+l} = f _{i n+l} +f _{ij n+l} et f _{j n+l} = f _{j n+l} + /,,· _Β+1 = f j,n+l f ij,n+l '

Comme la première technique, cette deuxième technique exploite le fait d'avoir figé certaines des particules en ne recalculant pas les forces d'interactions entre particules figées. Elle effectue la soustraction des forces correspondant aux anciennes positions et l'ajout de celles correspondant aux nouvelles positions. Elle ne comporte pas la longue opération de dressage des listes et de comparaisons des paires de chaque liste. En revanche, le volume des forces d'interaction entre deux particules à calculer est plus important que celui à réaliser dans la première technique.

On notera par ailleurs que dans les exemples ci-dessus, des forces d'interaction ont été considérées entre les particules considérées deux à deux pour calculer les potentiels d'interaction et la mise à jour de ces potentiels. Néanmoins l'invention permet aussi de réduire la charge de calcul correspondante lorsque le calcul du potentiel fait intervenir des forces d'interaction entre k particules, k étant strictement supérieur à 2. Dans ce cas, on calcule le potentiel d'interaction courant à partir du potentiel d'interaction déterminé à l'étape de simulation précédente, en tirant avantageusement parti du fait que la force d'interaction entre k particules est inchangée entre l'étape de simulation courante et l'étape précédente (et par conséquent n'est pas à recalculer) lorsque les k particules sont des particules figées. On retranche alors aux forces totales exercées sur les particules les forces calculées à l'étape précédente qui sont relatives aux k uplets de particules comportant des particules qui ont bougé entre l'étape de simulation précédente et l'étape courante. On calcule et on ajoute aux forces totales exercées sur les particules ainsi obtenues les forces courantes relatives aux k uplets de particules comportant des particules qui ont bougé, en fonction de leurs nouvelles positions.

Pour k=2, ou pour k étant différent de 2, des opérations semblables peuvent être mises en œuvre pour mettre à jour l'énergie potentielle du système, en considérant l'énergie potentielle comme la somme des énergies potentielles entre au plus k particules. Des opérations semblables peuvent également être mises en œuvre pour mettre à jour des valeurs ou des structures de données qui dépendent des positions d'au plus k particules, pour k entier quelconque supérieur ou égal à 1 .

Par exemple, dans une simulation relative à un ensemble de 5 particules ou plus, les informations à calculer comportent le centre de gravité des 5 particules considérées, qui évolue au cours du temps, mais que l'on veut ne déterminer que tous les 10 pas de temps de simulation. Si les termes de la matrice inverse de masse adaptative correspondant à ces 5 premières particules ont été mis à zéro entre l'instant où on a déterminé le centre de gravité pour la dernière fois et l'instant courant, alors les particules n'ont pas bougé, et il n'est pas nécessaire de mettre à jour le centre de gravité.

Dans un autre exemple, si les termes de la matrice inverse de masse adaptative correspondant aux quatre premières particules considérées ont été mis à zéro entre l'instant où on a déterminé le centre de gravité pour la dernière fois et l'instant courant, - mais - que les termes correspondant à la cinquième particule ne l'ont pas été à un certain instant (et donc que la 5eme particule a bougé depuis le dernier calcul du centre de gravité), on met à jour le centre de gravité g de manière incrémentale :

- multiplier g par 5 ;

- retrancher de g l'ancienne position de la 5eme particule ;

- ajouter à g la nouvelle position de la 5eme particule ;

- divise g par 5.

On tire ainsi parti du fait que les termes de la matrice inverse de masse adaptative correspondant aux 4 premières particules étaient à zéro. Dans un autre mode de réalisation, l'invention propose un procédé et un dispositif permettant d'accélérer le calcul des simulations d'ensemble d'objets. L'utilisation de l'Hamiltonien adaptatif permet, pendant la simulation, d'activer ou de désactiver des degrés de liberté, en position, d'objets vérifiant certains critères. Le volume de calculs nécessaire pour mettre à jour les forces ou l'énergie potentielle relative(s) à ces objets peut ainsi être diminué.

Dans le cas d'une implémentation d'une simulation du comportement d'un ensemble E dans un ensemble statistique NVT (ensemble E à nombre de particules, volume et température constants) considérée maintenant ici à titre d'illustration, on utilise par exemple, sur la base du Hamiltonien adaptatif décrit ci-dessus, une simulation suivant la dynamique de Langevin (cf. par exemple « Free energy computations : a mathematical perspective », T. Lelievre et al., Impérial Collège Pr, 2010).

Les équations de la dynamique de Langevin sont :

dq, = V _p H _A (q _t , p _t )dt,

Formules (4)

dp, = -^ _q ^H A (Q, . P _t )dt - f7 _p H _A (q _t , _Pt )dt + adW _t

où :

- t -> dW _t est une fonction de mouvement brownien standard à 3N dimensions, et - σ et γ sont des matrices réelles 3N ^* 3N, satisfaisant la relation de fluctuation-dissipation suivante : σσ ^τ = 2γΙ β , avec fi = \lk _B T , fc _B étant la constante de Boltzmann et T la température ;

pH _A (q _t , p _t ) est le gradient du Hamiltonien adaptatif par rapport à la variable P ;

V _q H _A (q _t , p _t ) est le gradient du Hamiltonien adaptatif par rapport à la variable q.

Dans ce cas d'une simulation NVT, pour l'algorithme d'intégration, le calcul d'un pas de temps peut être effectué ainsi : un demi-pas de temps pour la partie Langevin des équations, un pas de temps pour la partie Hamiltonienne des équations et à nouveau un demi-pas de temps pour la partie Langevin des équations.

On obtient alors les formules (5) suivantes :

_ _{n (} dH _A (q _n , p _n+l/2 ) dH _A (q _n , p _n+l/2 ) h ïh_

àg„ àp _n+1/2 2 2

où {G _k } est une séquence de vecteurs aléatoires Gaussien indépendants et identiquement distribués avec une moyenne nulle et une covariance égale à la matrice Identité.

Une équation des formules (5) comporte le terme p _n+ i dans les termes de droite et de gauche. Pour résoudre cette équation implicite, un algorithme à point fixe est utilisé, par exemple.

Un programme similaire au programme P décrit ci-dessus est adapté pour mettre des étapes similaires aux étapes 101 , 102, 103, en remplaçant, dans ces étapes, la prise en compte des équations (3) propres au cas NVE, par celles des équations (5) propres au cas NVT, pour l'actualisation des valeurs de p _n et q _n sur la base du Hamiltonien adaptatif H _A selon l'invention.

Des valeurs moyennes peuvent en outre être calculées lors de la simulation dans un ensemble NVT, en utilisant la simulation adaptative (i.e. utilisant un Hamiltonien adaptatif) selon l'invention, de manière à déterminer les valeurs qui auraient été calculées en réalisant la simulation avec un Hamiltonien classique.

En effet,

1 T

H _A =— p .Φ(ρ,ς) .p + V(q) s'écrit aussi : H _A = ^p ^T .M ^~l .p + V(q) + ^p ^T .(<î>(p,q).p -M- ^l ).p soit

H _A = H + V _A (p,q) , où V _A (p,q) = -^p ^T .(<î>(p,q).p -M ^~l ).p et H est le Hamiltonien standard (ie non adaptatif).

Si on nomme (A la valeur moyenne de la variable A obtenue à l'aide des paramètres p et q obtenus avec une simulation utilisant un Hamiltonien standard H, et (A la valeur moyenne de la variable A obtenue à l'aide des paramètres p et q obtenus avec une simulation utilisant un Hamiltonien adaptatif H _A , on peut retrouver la valeur de (A _H à partir de la valeur de (A _FL par l'égalité suivante : H où k _B T , est le produit de la

constante de Boltzmann k _B et de la température T ,

On peut démontrer, à partir de cette égalité, que lorsque la variable A dépend uniquement des positions et que l'Hamiltonien adaptatif est séparable (ie Φ ne dépend que des moments et pas des positions), alors (A) _H = (A)

Donc les valeurs moyennes obtenues à l'aide d'un Hamiltonien adaptatif sont égales à celles obtenues avec le Hamiltonien standard, ce qui est avantageux.

On notera que les intégrateurs d'Euler ou Langevin ont été considérés ci-dessus pour la mise en œuvre de l'invention. Celle-ci peut néanmoins être utilisée avec tout type d'algorithmes d'intégration.

Dans le cas considéré ci-dessus, le mouvement d'une particule était figé dans toutes les dimensions de l'espace de déplacement considéré. Dans un autre mode de réalisation, le mouvement d'une particule est figé sur seulement 1 ou certains des axes de déplacement, ce qui peut être utile pour étudier certains types de mouvements.

On réduit alors le nombre des calculs nécessaires pour les actualisations des composants des forces parallèles à cet ou ces axes.

Dans le mode de réalisation décrit ci-dessus, la position de la particule est figée lorsque son énergie cinétique est inférieure à un seuil.

Dans un autre mode de réalisation, combinable avec les précédents ou non, on fige une particule pendant au moins un pas de temps de simulation, lorsque son moment p prend certaines valeurs déterminées (valeurs discrètes, ou une ou différentes plages de valeurs), voire lorsqu'un couple comportant le moment p et la position q prend certaines valeurs fixées.

Dans un mode de réalisation, on fige la position de groupes de particules. Par exemple, on n'affecte la valeur nulle aux termes diagonaux de la matrice inverse de masse adaptative relatifs à une particule a,, d'énergie cinétique e i =1 à 10, que si toutes les valeurs d'énergie cinétique e _Y à e _l0 sont inférieures à un certain seuil. Cette disposition peut accélérer les calculs, ou permettre des calculs plus précis, dans le calcul de certains potentiels.

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 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, 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 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ù e =0 et ef =0.01 kcal/mol (i=1 à

N), 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 à 2,5, RMSD=0,01 14 S et Aq _max =0,18 S, où S est la distance d'équilibre dans le potentiel de Lennard-Jones utilisé.

Pour la simulation adaptative où e\ = 0,05 et ef =0,1 (i=1 à N), 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 à 5, et RMSD=0,0612 S et Aq _max =0,3 S.

Pour la simulation adaptative où e = 0,625 et ef =0,7 (i=1 à N), 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 à 10, et RMSD=0,359 S et Aq _max =13,94 S.

Des gains appréciables sont également constatés en mettant en œuvre l'invention relativement à des ensembles de particules NVT. Ainsi un procédé selon l'invention permet d'accélérer les calculs de manière importante, avec une altération faible des comportements.