WO2013174923A1 - Procédé de simulation d'un ensemble d'éléments, programme d'ordinateur associé - Google Patents

Procédé de simulation d'un ensemble d'éléments, programme d'ordinateur associé Download PDF

Info

Publication number
WO2013174923A1
WO2013174923A1 PCT/EP2013/060622 EP2013060622W WO2013174923A1 WO 2013174923 A1 WO2013174923 A1 WO 2013174923A1 EP 2013060622 W EP2013060622 W EP 2013060622W WO 2013174923 A1 WO2013174923 A1 WO 2013174923A1
Authority
WO
WIPO (PCT)
Prior art keywords
elements
simulation
instant
value
current
Prior art date
Application number
PCT/EP2013/060622
Other languages
English (en)
Inventor
Svetlana ARTEMOVA
Stéphane REDON
Original Assignee
Inria Institut National De Recherche En Informatique Et En Automatique
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Inria Institut National De Recherche En Informatique Et En Automatique filed Critical Inria Institut National De Recherche En Informatique Et En Automatique
Priority to EP13725145.0A priority Critical patent/EP2856361A1/fr
Priority to RU2014146944A priority patent/RU2014146944A/ru
Priority to KR1020147036027A priority patent/KR102082777B1/ko
Priority to US14/402,116 priority patent/US20150134310A1/en
Priority to CN201380038556.3A priority patent/CN104508667B/zh
Publication of WO2013174923A1 publication Critical patent/WO2013174923A1/fr

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Definitions

  • the potential energy in cases is for example equal to the [or is a function of] the potential of interactions V (q) between the elements whose interaction forces can be derived, q being a vector indicating the positions of the elements (in a more general case the potential of interactions can be also dependent on the moments of the elements),
  • the simulation of a set of elements makes it possible to study the behavior of such a set and to analyze its properties: displacements in terms successive positions and moments of elements, correlations of displacements between elements, changes of structure, increases and decreases of interactions between elements, configurations adopted on average, evolutions of associated energies, etc.
  • the elements may represent mechanical bodies, for example celestial or fluid, particles such as atoms or molecules, for example proteins, fluids, etc.
  • a common way of simulating a set of elements is to consider the Hamiltonian of the set, and to derive motion equations.
  • WO 2009/007550 describes for example a simulation technique of a set of elements.
  • the present invention aims to propose a new solution to reduce these problems.
  • the invention proposes a method for simulating a set of elements of the aforementioned type, implemented by computer and characterized in that said method comprises a step according to which, when the moment vector p takes certain determined values relating to at least one element, assigning a zero value to at least one diagonal term of the matrix M "1 relative to said element.
  • the invention makes it possible to reduce the volume and, consequently, the calculation time required for the determination of the potential energy, the interaction potential, the interaction forces, the positions and / or the moments of the elements.
  • the method of simulating a set of elements according to the invention further comprises one or more of the following features:
  • said method comprises a step according to which, for at least one of said elements, if a parameter representative of the kinetic energy of said element has a value lower than a first strictly positive threshold, a zero value is assigned to at least one diagonal term of the matrix M "1 relating to said element;
  • it comprises a step of determining the values of at least one information, at successive simulation instants on the basis of said Hamiltonian, said step taking advantage of the fact that the values of the information relating to a k-tuple of elements , with k integer greater than or equal to 2, for which a zero value has been assigned to the diagonal terms of the matrix M "1 at a previous simulation instant, are therefore unchanged between at least said previous simulation instant and the instant simulation method, and calculating a value of the information relating to a given element, at a current simulation time by implementing the following steps when zero values have not been assigned to the diagonal terms of the matrix for each element of a k-tuple of elements to which said given element belongs:
  • a current list of pairs of elements separated by a distance less than a given threshold is set at a current simulation instant, and compared to a previous list of pairs of elements separated by a smaller distance.
  • the value of information relating to a given element, at a current simulation instant is calculated on the basis of the pairs comprising said given element by implementing the following steps :
  • the value of the information relating to said given element at a current simulation instant is determined by adding, to the work value, the values of the information relating to said given element and associated with the other element of the pairs considered if the considered pair is present only in the current list or if the vector connecting said element to the other element of the pair has varied between the previous simulation instant and the current simulation instant;
  • a current list of k-tuples of elements satisfying certain conditions, with k greater than or equal to two, is set at a current simulation instant, and compared to a previous list of k-tuples elements satisfying said conditions at a previous simulation time,
  • determining the value of the information relating to said element at the current simulation instant by adding, to the temporary value, the values of the information relating to said element and associated with said k-tuples at the current simulation instant, when said -uplets are present only in the current list or when the information associated with said k-tuplets has changed between the previous simulation instant and the current simulation instant (for example when relative positions of the k elements in the k-tuple have changed);
  • the space for locating the elements is partitioned into cells and each element, at each preceding simulation instant and a current simulation instant, is associated with a membership cell according to determined positional coordinates at said instant of simulation, and according to which, for the first elements such that the terms of the matrix M "1 relating to said first elements have not been assigned to a zero value at a current simulation time, the following steps are implemented:
  • the membership cell of the first elements is determined at the previous simulation instant
  • the second elements located at the preceding simulation instant at a distance less than a given threshold of said first element are determined. ; calculating a work value by subtracting from the value of information relating to said first element and determined at the previous simulation time, the values of said information relating to said first element and associated with said second elements;
  • the new membership cell of the first elements is determined at the current simulation instant
  • the third elements located at the current simulation instant at a distance less than a given threshold of said unit are determined; first element;
  • the information relating to said element comprises the potential energy of said element and / or the interaction force applied to said element;
  • step of determination of information I comprises, at certain instants of simulation, a step of determination of information I, said step advantageously taking advantage of the fact that this information I is unchanged and does not need to be determined again when it has been determined to a previous simulation instant and that a zero value has been assigned to a corresponding set of diagonal terms of the matrix M "1 between at least said previous simulation instant and the current simulation instant (" corresponding set diagonal terms "diagonal terms that influence the value of information I, ie information I does not change when these terms are nil);
  • the present invention proposes a computer program for simulating a system of elements, comprising software instructions for implementing the steps of a method according to one of claims 1 to 12 when a program execution by calculation means.
  • FIG. 1 represents a device embodying an embodiment of the invention
  • FIG. 3 is a flowchart of the steps of a method in one embodiment of the invention.
  • FIG. 4 illustrates an embodiment of step 103
  • FIG. 5 illustrates another embodiment of step 103
  • FIG. 6 represents trajectory simulations of a particle connected to a fixed point, in phase space (p, q), at constant Hamiltonian.
  • H (p, q) p ⁇ T .M - 1 .p + V (q), where p is a vector indicative of the time of the particles, q a vector indicating the position of the particles, M "1 a diagonal matrix depending on the masses of these particles.
  • V (q) is the interaction potential between the N particles; it is a function of their position and it will be considered as independent of moments.
  • Hamiltonian Hamiltonian adaptive H A a so-called Hamiltonian Hamiltonian adaptive H A is defined, thus:
  • H A (p, q) - p. ⁇ I> (p, q) .p + V (q), where ⁇ (/ ?, ⁇ ?), Diagonal matrix
  • 3N * 3N called inverse matrix of adaptive mass, replaces M "1 and depends on the vector p, and possibly the vector q.
  • adaptive motion equations defining p and q which are the derivatives of the vectors p and q with respect to time t.
  • p and q are the derivatives of the vectors p and q with respect to time t.
  • the value of the Hamiltonian (adaptive according to the invention or standard) is constant in time
  • the adaptive equations of motion are: __aH __ay _j_ ⁇ fà (p, g)
  • the invention therefore consists in "freezing” the particles, by assigning them an infinite “pseudo-mass", when their kinetic energy passes under a certain value, the amount of movement of these particles not being fixed.
  • the function p i is a function including as variable the moment (in the particular case considered as an example, it is not therefore dependent on the position
  • the function p i may be a function, depending on the moment of the particle a, (and possibly of its position) other than the kinetic energy of course.
  • a particle is frozen whose moment (or torque including the moment and the position) takes predetermined values (discrete values or ranges of values).
  • the matrix ⁇ > (p, q) specifies how, and when, degrees of freedom in position of one or more particles are activated or deactivated during the simulation.
  • the curves C1, C2, C3, C4 each correspond to a respective constant value of the adaptive Hamiltonian.
  • the curve C1 corresponding to a Hamiltonian equal to 1.
  • the circle D corresponding to a constant value equal to 1 of a standard Hamiltonian, ie non-adaptive.
  • the area of the phase space where the particle is frozen is between the dotted lines B2 and B3 (it corresponds to a moment value in [-1, 1]).
  • the area of the phase space between the lines B1 and B2 and between the dotted lines B3 and B4 corresponds to a transition zone between the free and frozen states of the particle.
  • a computing device 1 shown in FIG. 1 is used to implement a simulation of a set E of N particles.
  • This device 1 comprises a computer including in particular a memory 2 adapted to store software programs and calculated parameter values successively described below (values of the matrix coefficients ⁇ , global interaction forces, partial, interaction potential, positions, moments ...), a microprocessor 3 adapted to execute the software program instructions and in particular the program P described below, and a man / machine interface 4, comprising for example a keyboard and a screen, respectively to enter instructions of a user and to display information for the user, for example curves such as that illustrated in Figure 6.
  • a computer including in particular a memory 2 adapted to store software programs and calculated parameter values successively described below (values of the matrix coefficients ⁇ , global interaction forces, partial, interaction potential, positions, moments ...), a microprocessor 3 adapted to execute the software program instructions and in particular the program P described below, and a man / machine interface 4, comprising for example a keyboard and a screen, respectively to enter instructions of a user and to display information for the user, for example curves such as that illustrated in Figure 6.
  • the memory 2 comprises the program P simulating the behavior of the set E of particles of the NVE type.
  • the program P comprises software instructions which, when executed on the microprocessor 3, are adapted to perform the following steps, with reference to FIG.
  • the steps 101, 101 b, 102, 103 are intended for the determination of the updated values respectively of the moment, the position and the global interaction force relative to each of the particles a ,.
  • the present value of other characteristic parameters of the behavior of the particles, at the moment h n + 1 can also be calculated, for example the current value of the potential energy of the system E, the value of the autocorrelation between particle velocities.
  • a first technique comprises the following steps.
  • a current list of the pairs of particles is drawn up, such that the distance between the particles of each pair at initialization is less than a threshold dO (when the distance between two particles is greater than dO, the interaction between these two particles is neglected) and the interaction force f ij 0 of the particle a, on the particle a, of each pair present in the current list is further evaluated, according to the distance separating them and according to the simulated force field, and stored.
  • dO when the distance between two particles is greater than dO, the interaction between these two particles is neglected
  • an element e ij Q also stored in memory 2, having the identifier of each of the two particles a ,, a, of the pair, the coordinates of the vector r 0 ! joining the two particles and starting from the particle a ,, and the interaction force f ij 0 exerted by the particle a, on the particle a, (which is equal to - f ji 0 , f jifi being the force of interaction exerted by the particle a, on the particle a).
  • step 103 During each iteration of step 103, the steps below are then implemented, with reference to FIG. 4.
  • a 103_a1 step is assigned as a starting value for the strength of overall interaction f i n + l exerted on each particle, the value of the force f in overall interaction calculated in the previous iteration.
  • a current list L a n + 1 of the pairs of interacting particles is drawn up, ie these are the pairs of particles such that the distance between the particles of each pair considered at computation time h n + i is below the threshold dO.
  • a step 103_c1 the current list of pairs L n + 1 is compared with the previous year list L pairs, ie made in the previous iteration (the n th iteration).
  • n is associated an element e ij n having the identifier of each of the two particles a ,, a, of the pair, the coordinates of the vector r n l ⁇ joining the two particles and starting from the particle a ,, calculated according to positions determined during the previous iteration for the particles a ,, a, and the value of the interaction force f ij n exerted by the particle a, on the particle a, (which is equal to - f ji n , f ji n being the interaction force exerted by the particle a, on the particle a,).
  • the interaction force f ij n + 1 exerted by the particle a, on the particle a is calculated, according to their respective position in particular; it is saved in memory in the element e ij n + l .
  • the invention by freezing particle positions, generates an increased number of pairs for which the vector between two particles, and therefore the interaction force between these two particles, remain unchanged.
  • step 103 make it possible not to recalculate all the components of the global interaction forces by taking advantage of the characteristics of a method according to the invention.
  • a second embodiment of step 103 makes it possible to exploit the advantages conferred by the invention without using a comparison of the lists of pairs of particles in interaction with the current iteration with respect to the previous iteration, but using a three-dimensional grid (if particle motion is considered in three-dimensional space, if the particles move in a plane, a two-dimensional grid may be sufficient).
  • an initial grid is created, considering a parallelepiped comprising all the particles and subdividing it into cells, for example cubic cells whose size of one side is greater than or equal to dO .
  • each particle a ,, i 1 to N, is assigned to the cell to which it belongs, depending on the position of the particle at the initialization step.
  • This force is equal to - f i , where Q is the interaction force exerted by the particle a, on the particle a ,.
  • the neighboring cells considered are the immediately adjacent cells, ie those which have at least one side in common with the given cell; in other embodiments, the neighboring cells considered are those located at r cells of distances of a cell immediately adjacent to the given cell.
  • a 103_a2 step is assigned as a starting value for the strength of overall interaction f i n + l exerted on each particle a ,, the value of the global interaction force calculated in the previous iteration in f .
  • a step 103_b2 for all the particles a, for which p i> n + 1 ⁇ 1 (ie the particles not considered as frozen), the particles a are determined, verifying the following conditions:
  • these particles a were located at the preceding iteration (n) in the cell in which the particle a, was positioned at the preceding iteration n, or the cells which are close to it (26 cells considered at maximum);
  • these particles a were at the previous iteration n at a distance from the particle a, less than dO;
  • composition of the grid so far considered is therefore that corresponding to the positions updated at the previous iteration (iteration n).
  • a step 103_c2 the composition of the grid is updated, by determining the current membership cells of all the particles a, for which p i> n + 1 ⁇ 1 (ie the particles not considered as fixed), according to the positions q i n + i of these particles corresponding to the iteration n + 1.
  • a step 103_d2 for all the particles a, for which p i> n + 1 ⁇ 1 (ie the particles not considered as frozen), the particles a are determined, satisfying the following conditions:
  • these particles a are situated at the current iteration (n + 1) in the cell in which the particle a, is positioned at the current iteration, or the cells which are close to it (26 cells considered at maximum);
  • these particles a are at the current iteration (n + 1) at a distance from the particle a, less than dO;
  • composition of the grid considered here is therefore that corresponding to the positions updated at the current iteration (iteration n + 1).
  • this second technique exploits the fact of having fixed some of the particles by not recalculating the interaction forces between frozen particles. It performs the subtraction of the forces corresponding to the old positions and the addition of those corresponding to the new positions. It does not involve the long process of listing and comparing pairs of each list. On the other hand, the volume of the interaction forces between two particles to be calculated is greater than that to be achieved in the first technique.
  • the invention also makes it possible to reduce the corresponding calculation load when the calculation of the potential involves interaction forces between k particles, k being strictly greater than 2.
  • the current interaction potential is calculated from the interaction potential determined in the previous simulation step, advantageously taking advantage of the fact that the interaction force between k particles is unchanged between the step of current simulation and the previous step (and therefore not to recalculate) when the k particles are frozen particles. Then the total forces exerted on the particles are subtracted from the forces calculated in the previous step which are related to kuplets of particles comprising particles that have moved between the preceding simulation step and the current step. To the total forces exerted on the particles thus obtained, the current forces relating to the kuplets of particles having particles which have moved, according to their new positions, are calculated and added to the total forces exerted on the particles thus obtained.
  • Similar operations can be implemented to update the potential energy of the system, considering the potential energy as the sum of the potential energies between at most k particles. Similar operations can also be implemented to update values or data structures that depend on the positions of at most k particles, for k any integer greater than or equal to 1.
  • the information to be calculated comprises the center of gravity of the particles in question, which changes over time, but that it is desired to determine only every 10 steps. simulation time. If the terms of the inverse adaptive mass matrix corresponding to these first 5 particles have been set to zero between the time when the center of gravity was last determined and the current time, then the particles do not have moved, and it is not necessary to update the center of gravity.
  • the invention provides a method and a device for accelerating the computation of object set simulations.
  • the use of the adaptive Hamiltonian allows, during the simulation, to activate or deactivate degrees of freedom, in position, objects verifying certain criteria. The volume of calculations necessary to update the forces or potential energy relative to these objects can thus be reduced.
  • pH A (q t , p t ) is the gradient of the adaptive Hamiltonian with respect to the variable P;
  • V q H A (q t , p t ) is the gradient of the adaptive Hamiltonian with respect to the variable q.
  • the calculation of a time step can be performed as follows: a half-step of time for the Langevin part of the equations, a time step for the part Hamiltonian equations and again a half-step of time for the Langevin part of the equations.
  • ⁇ G k ⁇ is a sequence of independent Gaussian random vectors identically distributed with a zero mean and a covariance equal to the Identity matrix.
  • a program similar to the program P described above is adapted to put steps similar to the steps 101, 102, 103, replacing, in these steps, the taking into account of the equations (3) specific to the NVE case, by those of the equations (5) specific to the NVT case, for updating the values of p n and q n on the basis of the adaptive Hamiltonian H A according to the invention.
  • Mean values can furthermore be calculated during simulation in an NVT set, using adaptive simulation (ie using an adaptive Hamiltonian) according to the invention, so as to determine the values that would have been calculated by performing the simulation with a Classical Hamiltonian.
  • the average values obtained using an adaptive Hamiltonian are equal to those obtained with the standard Hamiltonian, which is advantageous.
  • the movement of a particle was frozen in all the dimensions of the displacement space considered.
  • the motion of a particle is fixed on only 1 or some of the displacement axes, which may be useful for studying certain types of motion.
  • the position of the particle is frozen when its kinetic energy is below a threshold.
  • a particle is frozen for at least one simulation time step, when its moment p takes certain determined values (discrete values, or one or more ranges of values), or even when a pair comprising the moment p and the position q takes certain fixed values.
  • the position of groups of particles is frozen.
  • q is the vector of the coordinates of the particle a. at the last step of the adaptive simulation
  • q ⁇ is the coordinate vector of the same particle at the last step of the reference simulation.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

Procédé de simulation d'un système d'éléments, selon lequel le comportement desdits éléments est déterminé sur la base d'un Hamiltonien H du système d'éléments, tel que (formule I) avec p étant un vecteur indiquant les moments des éléments, q un vecteur indiquant les positions des éléments, M-1 étant une matrice diagonale fonction des masses des éléments et V étant l'énergie potentielle du système, ce procédé comportant 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.

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é 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 = -^pT .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) = ^ pT .M -1.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 (pi x, pi y, pi z) et la position q, de chaque particule a,, i=1 à N, s'écrit (qi x, qi y, qi z).
Les vecteurs p et q s'écrivent donc : p
Figure imgf000008_0001
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] = mi , pour i=1 à N, où mi 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 HA, 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 HA, 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 dHA 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](pi3C|i) = — (l-pt(qt,pt)), où mi est la masse de la particule a,.
On nomme ί(ρί,ςί) la valeur— (1- >;(<?;,/?;)) .
mi
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 à mi
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
Figure imgf000009_0001
]0,l[, la particule a, effectue une transition lisse entre ces deux comportements.
La nature deux fois dérivable de pi 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 < £,· <
Figure imgf000010_0001
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 ei figurant en abscisse, avec ε =1 et e{ =1 .1 .
Par exemple, une forme possible pour st (ei ) est - 6η5 + 15η4 - 10η3 + 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 et=
2mi
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 pi 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
Figure imgf000010_0002
Dans d'autres modes de réalisation, la fonction pi 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 à mi
l'Hamiltonien H standard du système E.
Lorsque A (<?; , .?; ) = 1 (ie position de la particule constante figée), Φ; {pi ,qi ) = 0 , par conséquent, q a une valeur nulle (en effet, pi é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 pi qi , pi ) 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 :
Figure imgf000012_0001
Ainsi, selon cette méthode, les formules (2) peuvent s'écrire ainsi :
Figure imgf000012_0002
Formules (3) qn+l = qn + (M -l (l - p(Pn+l ))Pn+l ~ Ρη Τ +ιΜ ~ι ^ ^ pn+l )h .
2 opn+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,.
P2
Dans le cas présent, les fonctions définies ci-dessus yoi (- -^) 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 pi 0 , la position qi 0 et la force d'interaction fi 0 de chaque particule a,, i= 1 à N, correspondant à un instant initial h0 .
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 hn+1= h0 + (n+1 )h, avec n entier≥0, h étant le pas de temps de simulation.
On nommera ci-dessous :
fij n+l , la force d'interaction exercée par la particule a, sur la particule a, (qui est égale à - fji n+l ), à l'étape de calcul hn+i ;
fi 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 hn+1 ; elle est donc égale à pi n+l , la valeur du moment de la particule a,, à Γ instant de calcul hn+1 ; qi n+l , la valeur de la position de la particule a,, à Γ instant de calcul hn+i ; pi n+l , la valeur prise par la fonction pi à Γ instant de calcul hn+i (comme vu ci- p2
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 pi 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 fi 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 pi n
Figure imgf000014_0001
Dans une étape 101 b, la valeur pi n+l est recalculée à partir de la nouvelle valeur du moment p n+1 : A +i =
Figure imgf000014_0002
Dans une étape 102, la valeur courante qi n+l de la position de chaque particule a,, i=1 à N, est à présent déterminée, d'après les formules (3) :
Figure imgf000014_0003
■
Dans une étape 103, la valeur courante fi 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 qi et des moments pi de chaque particule ai , avec i=1 à N, déterminée pour l'instant de calcul hn+1.
La valeur actualisée d'autres paramètres caractéristiques du comportement des particules a, à l'instant hn+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 < nmax, 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 fij 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 eij 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 r0 !' joignant les deux particules et partant de la particule a,, et la force d'interaction fij 0 exercée par la particule a, sur la particule a, (qui est égale à - fji 0 , fjifi étant la force d'interaction exercée par la particule a, sur la particule a, ).
La force d'interaction globale fi 0 s'exerçant sur la particule a, est égale à la somme des forces d'interaction fij 0 exercée sur la particule a, par les particules aj; 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 fi n+l s'exerçant sur chaque particule a, la valeur de la force d'interaction globale fi n calculée lors de l'itération précédente.
Dans une étape 103_b1 , une liste courante La 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 hn+i est inférieure au seuil dO.
A chaque paire dans la liste La n+1 est associé un élément eij 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 hn+1 , et la force d'interaction fij n+l exercée par la particule a, sur la particule a, (qui est égale à - i„+1 , ffii„+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 La n+1 est comparée à la liste précédente La n des paires, i.e. établie lors de l'itération précédente (soit la neme itération). A chaque paire dans la liste La,n est associé un élément eij n comportant l'identifiant de chacune des deux particules a,, a, de la paire, les coordonnées du vecteur rn 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 fij n exercée par la particule a, sur la particule a, (qui est égale à - fji n , fji 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 La 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 fi n+l s'exerçant sur la particule a,, respectivement de fj n+l s'exerçant sur la particule aj; la force d'interaction fij n calculée à l'étape 100 lors de l'itération précédente, respectivement la force d'interaction fji n = -fij n .
Si une paire est présente seulement dans la liste courante La 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 fij 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 eij n+l .
On ajoute la force d'interaction fij n+l à la force d'interaction globale sur la particule a en cours de détermination fi n+l , et on ajoute la force d'interaction fji n+l = -fij n+l à la force d'interaction globale en cours de détermination fj n+l s'exerçant sur la particule a,.
Si une paire est comprise à la fois dans la liste courante La n+i et dans la liste précédente La n, alors on compare les vecteurs rn l } et joignant les deux particules a,, a,. S'ils sont différents, on calcule la force d'interaction fij 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 fij n+l dans l'élément eij n+l . On ajoute à la force d'interaction globale en cours de détermination fi n+l s'exerçant sur la particule a,, la force d'interaction fij,n+i■ On ajoute à fj n+l s'exerçant sur la particule a,, la force d'interaction fji,n+i = ~fij,n+i■ En outre, on retire de la force d'interaction globale en cours de détermination fi n+l s'exerçant sur la particule a,, respectivement de fj n+i s'exerçant sur la particule a,, la force d'interaction fij n calculée à l'étape 100 lors de l'itération précédente, respectivement la force d'interaction fji n = -fij 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 fij 0 de la particule a, sur la particule a,.
Cette force est égale à - fji Q , fji 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 fi 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 fi n .
Dans une étape 103_b2, pour toutes les particules a, pour lesquelles pi>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 pj>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 fij n exercée par la particule a, sur la particule a, (et par là-même la force d'interaction fji 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 fij 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 fji n exercée par la particule a, sur la particule a, de la force d'interaction globale exercée sur a : ainsi on calcule fi n+y = fi n+y - fij n et fj >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 pi>n+1 < 1 (ie les particules non considérées comme figées), selon les positions qi n+i de ces particules correspondant à l'itération n+1 . Dans une étape 103_d2, pour toutes les particules a, pour lesquelles pi>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 pj>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 fij n+l exercée par la particule a, sur la particule a, (et par là-même la force d'interaction fji 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 fij 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 fji n+l exercée par la particule a, sur la particule a, de la force d'interaction globale exercée sur ¾ : ainsi on calcule fi n+l = fi n+l +fij n+l et fj n+l = fj 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, = VpHA(qt, pt)dt,
Formules (4)
dp, = -^q HA (Q, . Pt )dt - f7pHA (qt, Pt )dt + adWt
où :
- t -> dWt 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 = \lkBT , fcB étant la constante de Boltzmann et T la température ;
pHA (qt, pt) est le gradient du Hamiltonien adaptatif par rapport à la variable P ;
VqHA (qt, pt) 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 (dHA(qn , pn+l/2) dHA (qn , pn+l/2) h ïh_
àg„ àpn+1/2 2 2
Figure imgf000022_0001
où {Gk} 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 pn+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 pn et qn sur la base du Hamiltonien adaptatif HA 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
HA =— p .Φ(ρ,ς) .p + V(q) s'écrit aussi : H A = ^pT .M~l.p + V(q) + ^pT .(<î>(p,q).p -M-l ).p soit
H A = H + VA (p,q) , où VA (p,q) = -^pT .(<î>(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 HA, on peut retrouver la valeur de (A H à partir de la valeur de (A FL par l'égalité suivante : H où kBT , est le produit de la
Figure imgf000023_0001
constante de Boltzmann kB 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 eY à el0 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 ( Em / kB = 120 Kelvin, où Em 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 Aqmax.
Figure imgf000024_0001
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 Aqmax=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 Aqmax=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 Aqmax=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.

Claims

REVENDICATIONS
1 . Procédé de simulation d'un système d'éléments, selon lequel le comportement desdits éléments est déterminé sur la base d'un Hamiltonien H du système d'éléments, tel que H(p,q) = -^pT .M~l.p + V avec p étant un vecteur indiquant les moments des éléments, q un vecteur indiquant les positions des éléments, M"1 étant une matrice diagonale fonction des masses des éléments et V étant l'énergie potentielle du système, ledit procédé étant mis en œuvre par ordinateur et étant 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.
2. Procédé de simulation d'un système d'éléments selon la revendication 1 , selon lequel 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.
3. Procédé de simulation d'un système d'éléments selon la revendication 1 ou 2, selon lequel 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.
4. Procédé de simulation d'un système d'éléments selon l'une quelconque des revendications précédentes, selon lequel 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.
5. 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 à 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.
6. Procédé de simulation d'un système d'éléments selon l'une quelconque des revendications précédentes, selon lequel, à 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 selon lequel 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.
7. Procédé de simulation d'un système d'éléments selon l'une quelconque des revendications précédentes, selon lequel, à 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 selon lequel 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;
8. Procédé de simulation d'un système d'éléments selon l'une quelconque des revendications précédentes, selon lequel 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
9. Procédé de simulation d'un système d'éléments selon l'une quelconque des revendications précédentes, 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.
10. Procédé de simulation d'un système d'éléments selon l'une quelconque des revendications précédentes, comprenant, à 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.
1 1 . Procédé de simulation d'un système d'éléments selon l'une quelconque des revendications précédentes, comprenant, à 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.
12. Procédé de simulation d'un système d'éléments selon l'une quelconque des revendications précédentes, comprenant, à 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.
13. 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.
PCT/EP2013/060622 2012-05-25 2013-05-23 Procédé de simulation d'un ensemble d'éléments, programme d'ordinateur associé WO2013174923A1 (fr)

Priority Applications (5)

Application Number Priority Date Filing Date Title
EP13725145.0A EP2856361A1 (fr) 2012-05-25 2013-05-23 Procédé de simulation d'un ensemble d'éléments, programme d'ordinateur associé
RU2014146944A RU2014146944A (ru) 2012-05-25 2013-05-23 Способ моделирования совокупности элементов и компьютерная программа, применяемая при его осуществлении
KR1020147036027A KR102082777B1 (ko) 2012-05-25 2013-05-23 엘리먼트들의 세트를 시뮬레이팅하기 위한 방법 및 연관된 컴퓨터 프로그램
US14/402,116 US20150134310A1 (en) 2012-05-25 2013-05-23 Method for Simulating a Set of Elements, and Associated Computer Program
CN201380038556.3A CN104508667B (zh) 2012-05-25 2013-05-23 用于模拟一组元件的方法

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
FR1254838 2012-05-25
FR1254838A FR2991081A1 (fr) 2012-05-25 2012-05-25 Procede de simulation d'un ensemble d'elements, programme d'ordinateur associe

Publications (1)

Publication Number Publication Date
WO2013174923A1 true WO2013174923A1 (fr) 2013-11-28

Family

ID=46963807

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/EP2013/060622 WO2013174923A1 (fr) 2012-05-25 2013-05-23 Procédé de simulation d'un ensemble d'éléments, programme d'ordinateur associé

Country Status (7)

Country Link
US (1) US20150134310A1 (fr)
EP (1) EP2856361A1 (fr)
KR (1) KR102082777B1 (fr)
CN (1) CN104508667B (fr)
FR (1) FR2991081A1 (fr)
RU (1) RU2014146944A (fr)
WO (1) WO2013174923A1 (fr)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112052516B (zh) * 2020-08-13 2021-11-23 中国人民解放军军事科学院国防科技创新研究院 基于序列摆放的组件布局随机采样方法
CN113371231B (zh) * 2021-06-25 2022-03-08 四川大学 一种带约束的航天器姿态控制方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2009007550A2 (fr) 2007-06-20 2009-01-15 Inria Institut National De Recherche En Informatique Et En Automatique Dispositif informatique pour la simulation d'un ensemble d'objets en interaction et procédé correspondant

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH04192008A (ja) * 1990-11-27 1992-07-10 Pentel Kk 多関節ロボット旋回制御方式
JP4192008B2 (ja) * 2003-02-18 2008-12-03 株式会社渡辺商行 気化器及び気化器の洗浄方法並びに気化器を用いた装置

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2009007550A2 (fr) 2007-06-20 2009-01-15 Inria Institut National De Recherche En Informatique Et En Automatique Dispositif informatique pour la simulation d'un ensemble d'objets en interaction et procédé correspondant

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
BENNETT ET AL: "Mass tensor molecular dynamics", JOURNAL OF COMPUTATIONAL PHYSICS, LONDON, GB, vol. 19, no. 3, 1 November 1975 (1975-11-01), pages 267 - 279, XP024751538, ISSN: 0021-9991, [retrieved on 19751101], DOI: 10.1016/0021-9991(75)90077-7 *
HAIRER E.; LUBICH C.; WANNER G: "Geometric numerical integration : structure preserving algorithms for ordinary differential equations", vol. 31, 2006, SPRINGLER VERLAG
PLECHAC P ET AL: "Implicit mass-matrix penalization of Hamiltonian dynamics with application to exact sampling of stiff systems", MULTISCALE MODELING & SIMULATION, SOCIETY FOR INDUSTRIAL AND APPLIED MATHEMATICS, US, vol. 8, no. 2, 1 January 2009 (2009-01-01), pages 498 - 539, XP009167383, ISSN: 1540-3459 *

Also Published As

Publication number Publication date
FR2991081A1 (fr) 2013-11-29
RU2014146944A (ru) 2016-06-10
US20150134310A1 (en) 2015-05-14
CN104508667B (zh) 2018-09-14
KR102082777B1 (ko) 2020-02-28
EP2856361A1 (fr) 2015-04-08
CN104508667A (zh) 2015-04-08
KR20150013880A (ko) 2015-02-05

Similar Documents

Publication Publication Date Title
FR3052891A1 (fr) Procede d&#39;estimation du facteur d&#39;intensite des contraintes et procede de calcul de duree de vie associe
FR2994451A1 (fr) Procede hybride de non-appariement local pour la simulation d&#39;ecoulements multiphasiques dans les milieux fractures heterogenes
WO2011007058A1 (fr) Simulation d&#39;un agrégat évolutif du monde réel, notamment pour gestion de risque
EP3912099A1 (fr) Mise à l&#39;échelle de modèle composite pour réseaux neuronaux
Arnaiz-González et al. MR-DIS: democratic instance selection for big data by MapReduce
Ahmad Suhaimi et al. Integrated species distribution models: A comparison of approaches under different data quality scenarios
FR2911976A1 (fr) Procede de conception d&#39;un systeme comprenant des composants materiels et logiciels
WO2013107819A1 (fr) Procédé d&#39;optimisation de traitement parallèle de données sur une plateforme matérielle.
EP2856361A1 (fr) Procédé de simulation d&#39;un ensemble d&#39;éléments, programme d&#39;ordinateur associé
CN111158901B (zh) 计算图的优化方法、装置、计算机设备和存储介质
CN106648883B (zh) 基于fpga的动态可重构硬件加速方法及系统
Badri et al. Investigating the accuracy of test code size prediction using use case metrics and machine learning algorithms: An empirical study
Cheptsov HPC in big data age: An evaluation report for java-based data-intensive applications implemented with Hadoop and OpenMPI
CN109408722A (zh) 社区划分方法、装置、计算设备及存储介质
WO2014037236A1 (fr) Procédé de simulation d&#39;un ensemble d&#39;éléments
CN114461390A (zh) 结合多维度分析和关键路径法的评估方法及相关装置
Slavova Parallel triangular solution in the out-of-core multifrontal approach for solving large sparse linear systems
Javadi et al. Bandwidth modeling in large distributed systems for big data applications
WO2014053606A2 (fr) Procédé de traitement de données définissant un élément dans un espace e de dimensions d, programme d&#39;ordinateur associé
Lambert On the Effect of Replication of Input Files on the Efficiency and the Robustness of a Set of Computations
WO2011131248A1 (fr) Procédé et appareil de compression/décompression de données sans perte
FR3039677A1 (fr) Procede de conception de pieces mecaniques, notamment d&#39;aubes de turbomachine
WO2019129958A1 (fr) Procede de stockage de donnees et procede d&#39;execution d&#39;application avec reduction du temps d&#39;acces aux donnees stockees
CN116938769B (zh) 流量异常检测方法、电子设备和计算机可读存储介质
EP3224656B1 (fr) Procede de determination d&#39;un cube de proportion

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 13725145

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 2013725145

Country of ref document: EP

WWE Wipo information: entry into national phase

Ref document number: 14402116

Country of ref document: US

ENP Entry into the national phase

Ref document number: 2014146944

Country of ref document: RU

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

ENP Entry into the national phase

Ref document number: 20147036027

Country of ref document: KR

Kind code of ref document: A