EP2893473A1

EP2893473A1 - Method for simulating an assembly of elements

Info

Publication number: EP2893473A1
Application number: EP13756385.4A
Authority: EP
Inventors: Svetlana ARTEMOVA; Stéphane REDON
Original assignee: Institut National de Recherche en Informatique et en Automatique INRIA
Current assignee: Institut National de Recherche en Informatique et en Automatique INRIA
Priority date: 2012-09-06
Filing date: 2013-08-26
Publication date: 2015-07-15
Also published as: FR2995109A1; CN104756115B; CN104756115A; US20150254378A1; WO2014037236A1

Abstract

The invention relates to a method for simulating a system of elements represented by a tree comprising leaf nodes each representing an element, and inner nodes including a root node R, on the basis of a Hamiltonian (formula I), p being the vector of momentum, q being the vector of the positions of the elements, and V being the potential energy of the system: when predetermined conditions are verified, the same translational movement is imparted on at least the descending elements of a given node of the tree, by defining the matrix M^-1 equal to Φ_R, given the recursive formula according to which, for any node A of the tree comprising k child nodes (formula II).

Description

Simulation method of a set of elements

The present invention relates to a method of simulating a set of elements, according to which the behavior of the elements is determined, at successive simulation steps, on the basis of a Hamiltonian associated with the element system (the sum of the kinetic energy and potential energy of the whole)

H = - ^ p ^T ^~ .M + V p, where p is a vector indicating elements the moments, V is the potential energy of the system and M ^{"1 of} a diagonal matrix depending on the masses of the elements (in the case, this matrix may be a function of the positions of the elements).

The potential energy V is, in some cases, a function of the positions of the elements only. In other cases, the potential energy V may also depend on the moments of the elements. The forces acting on the elements can be derived from this potential energy.

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 of positions and successive moments of the elements, correlations of displacements between elements, changes of structure, the increases and decreases of interactions between elements, the configurations adopted on average, the evolutions of the 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, to derive motion equations, and to deduce the motion of the elements according to these equations.

WO 2009/007550 describes for example a simulation technique of a set of elements.

The evolutions of the set of elements must sometimes be simulated over a long period, in order to be able to observe certain phenomena or to be able to calculate certain statistics. The calculation times, and the cost in calculation, of these simulations then become sometimes very important. Many methods have been proposed to accelerate the simulations of a set of elements and the collection of statistics.

The present invention aims to propose a new solution to reduce these problems. For this purpose, according to a first aspect, the invention proposes a method for simulating a set of elements of the aforementioned type, characterized in that said method comprises the steps according to which:

the system of elements is represented in a k-ary tree comprising leaf nodes, said leaf nodes each representing a respective element, and internal nodes including a root node

when determined conditions are verified, for the current simulation step, the same translational movement is imposed on at least the descending elements of a given node of the tree by defining the matrix

M ^"1 as being equal to Φ ₈ where R is the root node of the tree, given the recursive formula according to which, for any node A of the tree having k child nodes A _Y , A ₂ , A _k ,

Φ _Αι 0 0

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

0 0 Φ

the matrix E being a matrix of size dn _A ^* dn _A formed of n _A ^* n _A blocks of size d ^* d equal to the identity matrix of dimension d, n _A being equal to the number of elements descending from the node A, d being the dimension of the space in which the particles evolve, m _A is the sum of the mass of these n _A falling elements of the node A, p _A being a function of restriction of the node A between 0 and 1, whose value is equal to 1 when A is a leaf node or when the same movement has been imposed on the descending elements of the given node; in the case of a leaf node A _i , the matrix

Φ _Α is equal to the inverse mass of the particle represented by the node A _i multiplied by the identity matrix of dimension d.

The invention makes it possible to perform simulations that require a smaller calculation volume and, consequently, less computing time, to determine the behavior of the elements according to these simulations, for example the potential energy, the forces applied to the elements, positions and / or moments of the elements.

In embodiments, the method of simulating a set of elements according to the invention further comprises one or more of the following features: for the current simulation step, the same translational movement is imposed on at least the descending elements of a given node of the tree as a function of the value of a function of the moments of said elements; for the current simulation step, the same translational movement is imposed on at least the descending elements of a given node of the tree as a function of the value taken by, or

m _A is the sum of the mass of the elements of node A _; , and p _A ^s _i is the sum vector of the moments of the elements of the child node A _{ , \\. \\ is the norm of the vector, C is a positive constant;

for the current simulation step for an internal node A:

if ε _Α is smaller than a first threshold, p _A equal to 1 is fixed to impose the same translational movement on the descending elements of the node A;

if ε _Α is greater than a second threshold greater than the first threshold, p _A equal to 0 is fixed;

for the current simulation step, if ε _Α is between the first threshold and the second threshold, p _A equal to the value taken by an interpolation function of order 5 of the variable ε _Α is fixed;

the method comprises a step of determining the values of at least one piece of information, at successive simulation instants on the basis of said Hamiltonian, said stage taking advantage of the fact that the values of the information relating to the elements to which the same movement in translation has been imposed for the moment of current simulation, depend on the relative position of said elements and are therefore unchanged;

the information relating to said element comprises the potential energy of said element and / or the interaction force applied to said element.

According to a second aspect, the present invention provides a computer program for simulating a system of elements, comprising software instructions for implementing the steps of a method according to the first aspect of the invention in a execution of the program by calculation means.

These features and advantages of the invention will appear on reading the description which follows, given solely by way of example, and with reference to the appended drawings, in which: FIG. 1 represents a simulation of trajectories of a 2-particle system performed with the standard Hamiltonian H;

FIG. 2 represents trajectory simulations of a 2-particle system with the relative adaptive Hamiltonian H _A R according to the invention;

FIG. 3 represents a device embodying an embodiment of the invention;

FIG. 4 is a flowchart of the steps of a method in one embodiment of the invention;

FIG. 5 represents a flowchart of the substeps of step 101 of FIG. 4.

Let us consider the simulation of a set E of N particles a ,, i = 1 to N.

The Hamiltonian H (p, q) associated with set E (see for example "Understanding molecular simulation: from algorithms to applications", Frenkel D., Smit B.) is often written in the following way:

p being a vector indicating the moments of all the particles, q a vector indicating the positions of all 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 positions and it will be considered as independent of moments.

In a 3-dimensional space, for example, with a coordinate reference system (X, Y, Z), the moment p, of each particle a ,, i = 1 to N, is written (p _ix , p _iy , p) and the position q, of each particle a ,, i = 1 to N, is written (q _ix , q _iy , q).

P V

The vectors p and q are therefore written and q =

ΡΝ, Ζ Usually the matrix M ^"1 used in the prior art is a diagonal matrix of dimension 3N ^* 3N, whose terms M [3i-2, 3i-2] = M [3i-1, 3i-1] = M [3i , 3i] = m _i , for i = 1 to N, where m _i is the mass of the particle a ,.

This is the usual definition of the Hamiltonian H, which will be referred to below as the standard Hamiltonian.

According to the invention, a Hamiltonian, said Hamiltonian, is defined that is adaptive relative to H _AR , thus:

1 T

H _AR (p, q) = -p ¹ & (p, q) .p + V (q), (formula 1) in which (ρ, ο), 3N ^* 3N diagonal block matrix referred to as inverse adaptive mass matrix relative, replaces M ^"1 and depends on the vector p, and optionally on the vector q.

More precisely, the values taken by Φ (ρ, ο) are determined at each simulation step, the function Φ having been chosen so as to be able to restrict at least a subset of the particles to follow together the same movement with one or some no simulation. The particles that are part of the subset are determined by verification of a condition, for example to verify that before the current simulation step, they were driven by a close movement.

These arrangements result in the appearance of one or more subsets of particles, the particles of a subset being animated with the same displacement in translation during the simulation step considered (ie the vector joining any two of these particles remains constant during a time period in which these particles have the same displacement).

The particles of such a subset can thus be assimilated, for the simulation step considered, as part of a single rigid body.

The information that depends only on the relative positions of the particles, for example, the interaction forces between the particles, does not need to be updated at the end of the simulation step considered within the same sub-group. rigid set.

When the movement of particles of a subset is no longer restricted by the adaptive inverse mass matrix, each particle resumes its own motion.

In one embodiment of the invention, a continuous transition is applied between the two types of behaviors (restricted movement / free movement). From this adaptive Hamiltonian H _A R, we deduce adaptive motion equations, defining p and q which are the derivatives of the vectors p and q with respect to time t.

For example, in a case of implementation of a simulation in a set NVE (set E with number of particles, constant volume and energy) considered here by way of illustration, the value of the Hamiltonian (adaptive according to the invention or standard) is constant over time, and the adaptive equations of motion are: dp ₌ dH _AR _ dV l _pT d <i> (p, q)

dt dq dq 2 dq

(formulas 2)

The following examples relate to a set NVE, but of course the invention can be applied to other sets, for example, the whole NVT (set E particle number, constant volume and temperature).

Several distinct solutions can be used to define the particles constituting a rigid subset, without having to test all the possibilities of combination between the particles of the set.

A first solution is to impose constraints on determined pairs of rigid bodies that can be considered as candidates to merge into a single rigid body. For example, noting A _1; At _n , a set of rigid indexed bodies, the indices can be organized into a binary tree, so that it can only merge with A ₂ into a rigid body (A _! + A ₂ ), when a specific fusion condition is satisfied, A ₃ can only merge with A ₄ etc., and at a higher level in the binary tree, (A _! + A ₂ ) can only merge with (A ₃ + A ₄ ) etc. Such a tree comprises at the maximum hierarchical rank, a root node, named node R, and at the minimum hierarchical rank, the leaf nodes corresponding to the particles considered alone. All the nodes of the tree, except the leaf nodes, are called the internal nodes.

Another solution is to impose constraints on determined subsets of rigid bodies, possibly comprising more than two rigid bodies. We could thus organize indices into an n-ary tree, so that Ai can merge only with A ₂ , A ₃ , ..., A _n at the same time, and so on.

Another solution is to consider the graph derived from a search of the neighbors, the nodes of the graph being particles or rigid bodies, and an edge between two nodes indicates that the distance between the two nodes connected by the edge is less than a threshold distance.

The embodiment detailed below is based on a binary tree structure as mentioned with respect to the first solution; the tree structure has the advantage of allowing incremental and recursive updates of the decision metrics for the internal nodes in the hierarchy, and incremental updates, for example for the forces and partial energies.

In such a tree, each leaf node represents a single particle of the set, and each inner node represents a subset of particles composed of particles represented by the child nodes of said internal node. C will be called the set comprising node C and the descendants, direct or indirect, of node C.

According to the invention, if a body C, ie an internal node C of the tree, is composed of two bodies A and B (ie the node C comprises two child nodes A and B), the inverse adaptive mass matrix is defined. Φ ₀ corresponding to this node C with the following recursive formula:

(formula 3)

or :

the adaptive inverse mass matrix Φ _Α corresponding to the node A and the adaptive inverse mass matrix Φ _Β corresponding to the node B, are themselves defined according to the formula 3 if the nodes A and B are not leaf nodes; in the case of a leaf node A, respectively B, the matrix Φ _Α , respectively Φ _Β , is equal to the inverse mass of the particle represented by the node A, respectively the node B, multiplied by the identity matrix of dimension 3 ;

the matrix E is a matrix of size 3n _c ^* 3n _c , formed of n _c ^* n _c blocks 3 ^* 3 equal to the identity matrix of dimension 3, n _c being equal to the number of leaf nodes of the node C (ie set of leaf nodes in the direct child nodes of node C or in the descendants of these child nodes);

- m _c is the sum of the mass of the particles represented by these n _c leaf nodes of the node C;

p _c is the restriction function of body C.

Thus we consider the fusion of two rigid bodies A and B into a rigid body C (and conversely the bursting of the rigid body C composed of the two rigid nodes A and B) on the basis of p _c , the restriction function of the body C . More exactly, a body C or internal node C, is considered as rigid (ie all the leaf particles included in the node C or the descendants of the node C, are animated of the same movement in translation) when the function p _c takes the value

1, and otherwise it is considered non-rigid.

All the leaf nodes are considered rigid by definition, and this all the time that the simulation lasts, since they consist each of a single particle.

When the function p _c is 0, the internal node C is considered as free.

In the embodiment of the invention, the function p _c is defined recursively, so as to smooth the changeover between the values 0 and 1, according to the following formula:

1 if C is a leaf node,

(formula 4)

[μ (ε) ρ _Α ρ _Β in the opposite case, where μ (ε) € [θ, ΐ] is a function twice differentiable from ε

1 if 0≤e≤e ^r ,

μ (ε) 0 if ε≥ ε ^f , (formula 5) s (£) G [θ, ΐ] otherwise, with the thresholds ε 'and ε ^Γ defined so that ε'> ε ^Γ and s being a function two times differentiable from ε.

The twice-differentiable nature of μ makes it possible to preserve the stability of the simulation of the set E of particles.

For example, a possible form for ε (ε) is an order interpolation function

5, for example equal to - 6η ⁵ + 15η ⁴ - 10η ³ + 1, with η = or

δ = ε '- ε ^Γ represents the width of the transition region between the rigid and free behaviors of the body.

In the embodiment of the invention, ε ₀ is chosen depending on a relation of the moment of the two bodies A and B, with

s 2 s 2 s. s 2 2

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

+ - (formula 6)

2m _A 2m _B 2m _c 2m _A m _B m _c where A and B are the child nodes of node C,

where the vector p _c ^s is the sum of the moments of all the particles which are leaf nodes in the body corresponding to the node C: p _c ^s = ΣΡι, and is the sum ieC

vectors p _A ^s and p _B ^s : p _c ^s = p _A ^s + p _B ^s -

This invention can be generalized for a k-ary tree. In this case, for the parent node with k child nodes A 1 <i k k, the inverse matrix of the inertia is

The coefficient ½ in this formula and the formula 6 can be replaced by another constant, for example 1, which will not change the simulation if the thresholds ε 'and ε ^Γ are modified in the same way, ie multiplied by an equal factor double the other constant).

This invention can be generalized as differently: one can choose the restriction function p _c different for different internal nodes. For example, up to a certain level of the tree p _c is defined according to formulas 4 and 5, therefore, the nodes can be joined together and separated, and in the levels of the tree above this level, p _c is always equal to zero. In a particular embodiment, it is possible to define a subset of the system E which is permanently active (the particles of this subset are never merged with each other or with the rest of the system, so they always follow the During the simulation, an unrestricted, free movement (such an embodiment of the invention may for example be applied to a polymer in a solvent): all the particles of this subset must be placed in a separate subtree. and the values p _c must be set to 0 for all the internal nodes of this subtree and never be updated In all cases, the functions p _c must be defined for all the internal nodes before the beginning of the simulation and should not change during the simulation in order to have a stable simulation. According to the invention, the relative adaptive mass inverse matrix Φ used in the formula (1) of the relative adaptive Hamiltonian H _A R is chosen, equal to Φ _¾ which is the matrix defined according to the recursive formula 3 relative to the root node R of the tree representing the set E of particles: Φ (ρ, <?) = (ρ) = _Κ ■

1 τ

We then have the Hamiltonian H _AR , such that H _AR (p, q) = - p & _R (p) .p + V (q), which is separable because the matrix Φ _κ depends only on the moments (and not the positions) of the particles.

The equations of motion defined by formulas 2 then express themselves as follows:

^ _ dp _ dH _AR _ d V _^

dt dq dq '

(formulas 7)

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

at dp 2 dp

It is also clear from the formula 3 that if the node C is free {p _c = 0), then the matrix Φ ₀ is composed of two blocks on the diagonal corresponding one to the child node A and the other to the node son B.

In a first example of simulation according to the invention, it is considered a one-dimensional system comprising two particles P1, P2 of mass 1 connected to one another by a spring of stiffness 1. The binary tree corresponding to the system is a root node C whose two child nodes A and B correspond to these two particles. Given initial moments, the two particles oscillate in space with time.

The space-time trajectories calculated for these particles according to a simulation of the prior art based on the standard Hamiltonian H are shown in FIG. 1 as a function of the simulation time.

The trajectories d calculated for these particles according to a simulation on the basis of the relative adaptive Hamiltonian H _A R according to the invention are represented in FIG. 2 as a function of the simulation time, for different values taken by δ and ε ^f . Thus as indicated on these trajectories, at certain moments of the simulation (sometimes from the beginning), the moments of the two particles become close according to the formula 6. This occurs when the spring is almost compressed and almost decompressed. According to the invention, the root node C then becomes rigid. Consequently, in a general case, the son nodes of the rigid node C follow the same translational movement, in this case, because of the preservation of the moments, the particles stop. The trajectories 2 particles therefore become parallel. The moments of the particles continue to evolve however, which has the consequence that at a given pitch of the simulation, the rigidity condition of the node C is no longer satisfied and that the two particles resume their own motion. This is repeated periodically during the simulation.

As illustrated by the trajectories, for the same value of ε ^f with an increased value for δ, the transition region of the trajectory is wider. For the same value of δ with an increased value for ε ^f , the stiffening region is longer and flatter.

Thus according to the invention, the matrix Φ, defined as Φ _¾, specifies how, and when, degrees of freedom in the relative position of the subset (s) of particles are activated or deactivated during the simulation.

In one embodiment of the invention, a computer device 1 shown in FIG. 3 is used to implement a simulation of a set E of N particles, according to the principles of the invention explained above.

This device 1 comprises a computer including in particular a memory 2 adapted to store software programs and calculated parameter values successively described below (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 from a user and to display information intended for the user, for example curves such as those illustrated in FIG. 2.

In the embodiment of the invention considered, the memory 2 comprises the program P simulating the behavior of the set E of particles a, i = 1 to N on the basis of the integration over time of the equations corresponding to the formulas 7.

Since the relative adaptive Hamiltonian H _A R according to the invention is separable, many integration techniques can be used.

With reference to FIG. 4, the program P comprises software instructions which, when they are executed on the microprocessor 3, are adapted to implement the preliminary step 100 and iteratively the steps 101 to 104 during a n + 1 ^'th iteration of the program P corresponding to the moment of calculation _{n + 1} h = h ₀ + (n + 1) h where n entier≥0, where h is the pitch of simulation time.

In a previous initialization step 100, a binary tree representing the particles and the possibilities of stiffening two in two, as indicated above, is built for system E, initial moment and position values are set for the particles. Only the leaf nodes are then rigid.

In a step 101, the interaction forces exerted between the particles at the n + 1 ' ^th iteration, are determined.

In a step 102, the moments of the particles at the n + 1 ' ^th iteration, are determined.

In a step 103, the metrics relating to the internal nodes of the tree at the n + 1 ' ^th iteration, are determined. They include in particular p _c and e _c .

In a step 104, the positions of the particles at the n + 1 ' ^th iteration, are determined.

All these steps of updating the values of the forces, moments, metrics, positions, for the current iteration, are implemented according to the values calculated during the previous iteration.

The implementation of these steps 101 to 104 is detailed below.

In step 101, the updating of the value of the interaction forces can be accelerated, compared to the prior art, taking into account the stiffenings according to the invention, if the interaction forces depend solely on the positions related particles.

Thus according to the invention, it is not necessary to recalculate the forces within subsets of particles identified as rigid (according to the value taken by e _c according to formula 6) at the previous simulation iteration (the forces are determined before e _c ). In order to effectively update the forces by taking advantage of this advantage as a result of the stiffenings according to the invention, an incremental force update algorithm can be used.

This algorithm is based, first, on a hierarchy of bounding volumes

(a volume per node of the binary tree), for example, on the hierarchy of AABB's (in English "Axis-Aligned Bounding Volumes", the volumes including the subsystem and having a parallelepiped shape with the sides parallel to the axes of the coordinates, see for example Grudinin, S. and Redon, S., "Practical Modeling of Molecular Systems with Symmetries", Journal of Computational Chemistry 31, 9 (2010), pp. 1799-1814).

This force update algorithm uses secondly partial force tables. A partial force table corresponding to each internal node C comprises 3D vectors, the number of which is equal to the number of leaf nodes descending directly or indirectly from the node C. Each vector of this table is the interaction force acting on the particle represented by one of these leaf nodes in from all the other particles represented by the other leaf nodes descending from C. Each vector of the partial force table of the root node R corresponds to the global force acting on a particle. A partial force table corresponding to a leaf node is a null vector.

The force update at each iteration includes three main substeps (the potential energy can be updated at the same time):

In a step 101 1, the bounding volume hierarchy is updated, for example from the bottom of the tree upward (i.e. from the leaf nodes to the root). For each leaf node, the updated AABB box coincides with the last determined position of the particle represented by the leaf node. The box of each internal node is recalculated for example by including the boxes of its two child nodes.

The structure of the tree is unchanged, the surrounding volumes are updated according to the positions of the particles that have changed meanwhile.

In a step 101 _2, using this updated AABB's hierarchy, interaction lists are built recursively for the inner nodes of the tree.

For each node C, having two child nodes A and B, the interaction list thus constructed contains all the pairs of particles such that one of the particles of the pair belongs to the body A (ie is represented by a falling leaf node directly or indirectly from node A), while the other belongs to body B. This interaction list can be constructed for each inner node as described for other types of bounding volumes in R. Rossi, M. Isorce, S. Morin, J. Flocard, K.

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

Bioinformatics 2007 23 (13): 1408-1417.

During the first iteration, the interaction lists are built for all internal nodes. For all the other iterations, in the rigid nodes these lists do not change and therefore can not be recalculated.

In a step 101_3, an incremental update of the partial force tables is performed. When a node is identified as a rigid node, the forces between the particles of the node have not changed. Therefore, only the partial force tables of the non-rigid nodes must be updated.

They can be updated for example as follows recursively, starting from the leaf nodes to the root node.

For each non-rigid node C, in the first half of the partial force table relating to node C, the elements of the partial force table relating to node A are copied and in the second half of the partial force table relative to node C, the elements of the partial force table relating to the node B are copied. Then, for each pair of the interaction list set in sub-step 101_2 for node C, the interaction forces between the two particles of each pair in this list are calculated and added to the corresponding forces of the force table. partial. Finally, the forces taking into account all the particles are stored in the partial force table of the root node R.

The steps 101_2 and 101_3 can be combined, for example in a single traversal of the tree, to update the partial force table corresponding to a node C since the list of interactions corresponding to said node C is available.

The other information depending only on the relative positions of the particles can be updated in a similar way.

The step 102 of updating the moments of the particles is carried out using a conventional method.

Thus, for each particle a ,, by denoting p, _{n + 1} the value of the moment p, of the particle a ,, during the current iteration n + 1, this value is obtained by the following formula:

h is the time step of the simulation and f, _{n + i} is the overall interaction force acting on the particle, _i = 1 to N, which is due to the interactions performed by all the other particles of the system E, appearing in the partial force table of the root node R as it has just been determined at the current iteration.

The complexity of the "naive" update, at each time step, of the particle position from the equation q = Φ _Κ (p) p + -p ^T ^ - ^ -p, is 0 (N ³ ) because of the

2 dp

ΘΦ (p)

term - ^ -. This term is a derivative of a matrix with respect to a vector, so dp

a three-dimensional object, where, for i = 1 to N:

In the embodiment considered, the complexity can be reduced to O (N.logN) by a recursive update from leaf nodes of data structures Q _c , Rc and S _c for each internal node C with two child nodes A and B.

Noticing that

where p _A is the vector of the moments of all the particles which are leaf nodes of the node A and p _B is the vector of the moments of all the particles which are leaf nodes of the node B, the term Q _c = <i> _c (p _c ) p _c corresponding to the first term of the equation giving q:

(formula 9)

where {x} is the dimension vector equal to n _c and each of the n _c components is equal to the element x.

Similarly, considering for the node C, the term R _C = dPc (Pc)

dPc

dPc (Pc)

R _C = with p _c = μ (ε) ρ _Α ρ _Β , differentiating,

p

(formula 10)

The last term of the equation giving q can be updated recursively using the formula 9.

Indeed :

g - r * c (Pc) _

^Jc ^_ Pc -, Pc ^~

atp _c

rmule 1 1)

The function (·) represents the dot product.

Thus, step 103 includes updating the structures Q _c , R _c and S _c for each node C comprising two child nodes A and B.

It will be noted that step 100 further comprises a step of initializing these structures where for each internal node C: p _c = 0 and the values of Q _c , Rc and S _c are also zero, and for each leaf node F representing a particle: p _F = 1; Q _F = P _F lm _F ; S _F = 0, R _F = 0.

In step 103, the update metric Q _c, R c and S _c and e _c and p _c with respect to the node C is performed recursively along the nodes of the hierarchy from the root:

- If node C is internal, with two wires A and B:

- update the metrics of node A;

- update the metrics of node B;

- update e _c using the formula 6;

update p _c : p _c = μ (ε _ε ) ρ _Α ρ _Β ;

- update Q _c using the formula 9;

update R _c using the formula 10;

- update S _c using the formula 1 1;

- If the node C is a leaf node representing a particle a _F : QF = PF NIF ■ To optimize the computations, the update of the vectors Q _c , Rc and S _c can be carried out only for the nodes for which p _c > 0 because these metrics may not be used for free nodes, so there is no need to update them.

In the general case, these metrics can be updated other than by traversing the tree, as long as they are updated in the subtrees with the rigid node at the top.

Once the metrics are determined, the positions of the particles are updated in step 104.

The implementation of step 104 may be the following for updating the position of a node C, each internal node C being considered as having two child nodes A and B:

if p _c = 0 (the node C is a free node):

- update the position of node A;

- update the position of node B;

otherwise, for each leaf node k descending from node C, directly or indirectly: q. _{n + 1} = _n + q (Q _c [k] + 0.5S _c [k]) h q _{i n + i} being the value of the position of the particle a ,, for the iteration n + 1 considered and h the step of simulation time, Q _c [k] being the k ^th element of the vector Q _c and S _c [k] being the k ^th element of the vector S _c .

The position of the particles is thus determined.

In this embodiment, the one-to-one correspondence between the identifier of the particle i and its serial number among the child nodes of the node C must be established, for example during the initialization of the simulation, for all the particles and all the nodes. internal of the tree.

In this embodiment of updating the positions, the metrics of the nodes for which p _c > 0 are not used and therefore may not be calculated. If all the metrics are still updated during step 103, the positions of all the particles _t , l≤i≤N can also be determined as follows:

<? ,, _{"+ I} ⁼ 1i + (Q _R [k] + 0-5S _R [k]) h, where the index R corresponds to the root node.

For each internal node C of the tree, four vectors of length equal to the number of direct or indirect descending leaf nodes of the node C, are stored: Qc, Rc and S _c and a partial force table. The complexity in position update space is therefore O (N.logN). Time complexity is also O (N.logN) since these structures must be updated at each time step, and Q _c is always updated for all nodes, including leaf nodes.

In the general case, the positions can be updated other than by traversing the tree.

Then, if the maximum duration of the simulation is not reached, a new iteration of the program P is carried out.

By way of illustration, four 2-dimensional simulations of the evolution of an NVE set of N = 5930 particles a,, i = 1 to N, with a mass of 1 g / mol each, using a Lennard-potential. Jones (E _m / k _B = 120 Kelvin, where E _m is the minimum energy, equilibrium distance S = 3.4 angstroms, cutoff distance 8 angstroms, the potential being truncated smoothly between 7.5 and 8 angstroms), were carried out starting from a shock triggered by the high speed sending of a particle on the initially immobile system: a reference simulation, based on a standard Hamiltonian and three adaptive simulations, that is to say ie using a relative adaptive Hamiltonian of a method according to the invention as described above (no time of size 0.0488 femtoseconds (fs), 7000 time steps, total simulation time equal to 342 fs). For each of the simulations using a relative adaptive Hamiltonian, the square root of the fluctuation from the standard simulation, named RMSD, is given along with the maximum particle displacement error Aq _max .

, where q, is the coordinate vector of the particle a, at the last step of the adaptive simulation, and qf is the coordinate vector of the same particle at the last step of the reference simulation.

For example, for adaptive simulation where ^{r r} = 0.01 kcal / mol and f = 1.01 kcal / mol (<5 = 1 kcal / mol), we obtain an acceleration factor of the computation time necessary to put implement this adaptive simulation with respect to the calculation time relative to the reference simulation, of a value equal to 1, 5, RMSD = 8.6 S and Aq _max = 89.4 S, where S is the distance of balance in the potential of Lennard-Jones used.

For adaptive simulation where e ^r = \ kcal / mol and e ^f = 2 kcal / mol (<5 = 1 kcal / mol), we obtain an acceleration factor of the computation time necessary to implement this adaptive simulation by relative to the calculation time relative to the reference simulation, a value equal to 1, 55, and RMSD = 8.7 S and Aq _max = 89.4 S.

For adaptive simulation where £ ^r = 3 kcal / mol and e ^f = A kcal / mol (<5 = 1 kcal / mol), we obtain an acceleration factor of the computation time necessary to implement this adaptive simulation by compared to the calculation time relative to the reference simulation, of a value equal to 1, 6, and RMSD = 12.3 S and Aq _max = 89.4 S.

Thus, a method according to the invention makes it possible to accelerate the calculations, with a potentially low alteration of the behaviors.

The results depend on the tree representation chosen for the particle system. For adaptive simulations above, the binary tree was built from the bottom up by dividing the system into halves.

In the considered embodiments of the invention, the fusion of two subsets into a rigid set has been considered, in other embodiments, the number of subsets merged into a rigid set is higher than together.

The transition region between the two rigid and free states can be taken of greater or smaller width.

Claims

1. A method of simulating a system of elements, according to which the behavior of said elements is determined, at successive simulation steps, on the basis of a Hamiltonian

H the system of elements, such that H (p, q) = - ^ p ^{~ l} ^T .M .p + V p is a vector indicating elements the moments, q indicating a vector element positions, M ^{" 1} being a matrix function of the masses of the elements and V being the potential energy of the system,

characterized in that said method comprises the following steps:

Φ _Αι 0 0

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

0 0 Φ

2. A method of simulation of a system of elements according to claim 1, wherein it imposes, for the current simulation step, the same movement in translation to at least the descendant elements of a given node of the tree in function of the value of a function of the moments of said elements.

3. A method of simulating a system of elements according to claim 1 or 2, wherein it imposes, for the current simulation step, the same movement in translation to at least the descendant elements of a given node of the tree as a function of the value taken by ε _Α =, where m, is the sum of the mass of

elements of node A _; , and p ^s is the sum vector of the moments of the elements of the child node A _i , ||. || is the norm of the vector, C is a positive constant.

4. A method of simulating a system of elements according to claim 3, wherein for the current simulation step for an internal node A:

if ε _Α is greater than a second threshold greater than the first threshold, p _A equal to 0 is fixed.

5. Simulation method of a system of elements according to claim 4, wherein, for the current simulation step, if ε _Α is between the first threshold and the second threshold, p _A is fixed equal to the value taken by an interpolation function of order 5 of the variable ε _Α .

6. A method of simulating a system of elements according to any one of the preceding claims, comprising a step of determining the values of at least one information, at successive simulation instants on the basis of said Hamiltonian, said step pulling This is because the information values relating to elements to which the same translational movement has been imposed for the current simulation instant depend on the relative position of said elements and are therefore unchanged.

7. A method of simulating a system of elements according to any one of the preceding claims, wherein the information relating to said element comprises the potential energy of said element and / or the interaction force applied to said element.

8. Computer program (P) for simulating a system of elements, comprising software instructions for implementing the steps of a method according to one of claims 1 to 7 during execution of the program by calculation means.