CN104756115B - 用于模拟一组元件的方法 - Google Patents

Abstract

本发明涉及一种基于汉密尔顿算符(等式I)对由树表示的元件的系统进行模拟的方法，树包括叶节点和包括根结点R的内部结点，每个叶节点表示一个元件，p是力矩的矢量，q是元件的位置的矢量，以及V是系统的势能：－当预定的条件被验证时，通过把矩阵M^‑1定义为等于Φ_R，对树的给定结点的从属的至少一些元件施加同样的平移运动，对于包括k个子结点的树的任何结点A，有递归等式(等式II)。

Description

用于模拟一组元件的方法

技术领域

本发明涉及一种用于模拟一组元件的方法，根据该方法基于与元件系统(该组元件的动能和势能之和)相关联的汉密尔顿算符 来在相继的模拟步长确定元件的行为，其中，p是表示元件力矩的矢量，V是系统的势能，以及M^-1是这些元件的质量的函数的对角矩阵(在此情况下，该矩阵可以是元件位置的函数)。

背景技术

在此情况下势能V例如仅是元件位置的函数。在其它情况下，势能V还可以依赖于元件的力矩。作用于元件上的力可以从该势能推导。

一组元件的模拟允许研究这样的一组元件的行为并分析其属性：在元件的相继位置和力矩方面的位移，元件之间的位移的相关性，结构的变化，元件之间的相互作用的增加和减小，平均采用的配置，相关联能量的演变等等。这些元件可以表示机械体，例如天体或流体；颗粒，例如原子或分子，例如蛋白质，流体等。

模拟一组元件的常规方式是考虑该组元件的汉密尔顿算符并从中推导力矩方程，并根据这些方程推导元件的力矩。

WO2009/007550例如描述了模拟一组元件的技术。

该组元件的演变有时必须被长时段上模拟以便能够观察某些现象或能够计算某些统计量。这种模拟的计算时间、以及计算成本有时变得非常大。已经提出的很多方法来加速对一组元件的模拟并收集统计量。

发明内容

本发明旨在提出用于克服这些问题的新方案。

为此，根据第一方面，本发明提出一种对指定类型的元件的集合进行模拟的方法，其特征在于所述方法包括以下步骤：

－把元件的系统表示为包括叶节点和内部结点的k叉树，所述叶节点每个都表示相应的元件，内部结点中包括根结点，

－当所确定的条件被验证时，对于当前模拟的步长，通过把矩阵 M_- ¹定义为等于Φ_R，对树的给定结点的从属的至少一些元件施加同样的平移运动，其中R是树的根结点并且对于包括k个子结点A₁、A₂，…， A_k的树的任何结点A，有以下递归等式，根据该等式：

矩阵E是由各自尺寸为d*d的n_A*n_A个块形成的尺寸为dn_A*dn_A的矩阵，块等于尺寸为d的单位矩阵，n_A等于结点A的从属的元件的数量，d是其中颗粒演变的空间的尺度，m_A是从属于结点A的这n_A个元件的质量之和，ρ_A是介于0和1之间的结点A的约束函数，该函数的值当A是叶节点时或者当同样的运动被施加到从属于给定结点的元件时等于 1；在叶节点A_i的情况下，矩阵等于由结点A_i表示的颗粒乘以尺度为d的单位矩阵所得的逆矩阵。

本发明允许执行以下模拟：该模拟需要较小计算体积并因此较少计算时间来确定根据这些模拟的元件的性能，例如元件的势能、应用于元件的力、位置和/或力矩。

在实施方式中，根据本发明的对元件的集合进行模拟的方法还包括一个或更多个以下特征：

－对于当前模拟步长，对树的给定结点的从属的至少一些元件根据所述元件的力矩函数的值来施加同样的平移运动；

－对于当前模拟步长，根据由所取的值对树的给定结点的从属的至少一些元件施加同样的平移运动，其中是结点A_i的元件的质量之和，是子结点A_i的元件的力矩的矢量和，||·|| 是矢量的范数，C是正常数；

－对于针对内部结点A的当前模拟步长：

－如果ε_A小于第一阈值，则把ρ_A定为1以便对从属于结点A的元件施加同样的平移运动；

－如果ε_A大于比第一阈值大的第二阈值，把ρ_A定为0；

－对于当前模拟步长，如果ε_A介于第一阈值和第二阈值之间，则把ρ_A定为等于由变量ε_A的5阶插值函数所取的值；

－该方法包括在基于所述汉密尔顿算符的相继模拟的时刻确定至少一个信息的值的步骤，所述步骤利用了以下事实：对于当前模拟的时刻与已经施加了同样的平移运动的元件相关的信息的值取决于与所述元件相关的位置并且因此不变化；

－与所述元件相关的信息包括所述元件的势能和/或应用到所述元件的相互作用力。

根据第二方面，本发明涉及一种对元件的系统进行模拟的计算机程序，包括当由计算装置执行该程序时用于实施根据本发明第一方面所述的方法的步骤的软件指令。

附图说明

参照附图，通过对以下仅以示例方式提供的描述的阅读，本发明的其它特征和优点将会变得更加清楚，在附图中：

－图1示出用标准汉密尔顿算符H执行的具有2个颗粒的系统的轨迹的模拟；

－图2示出根据本发明的自适应相对汉密尔顿算符HAR进行的具有2个颗粒的系统的轨迹的模拟；

－图3示出实施本发明的实施方式中的设备；

－图4是根据本发明的实施方式中的方法的步骤的流程图；

－图5示出图4的步骤101的子步骤的流程图。

具体实施方式

考虑到N个颗粒a_i的一个组E的模拟：i＝1至N。

与该组E相关联的汉密尔顿算符H(p,q)常写作如下(参见例如“理解分子模拟：从算法到应用”，Frenkel D.,Smit B)：

p是表示所有颗粒的力矩的矢量，q是表示所有颗粒位置的矢量，M^-1是作为颗粒的质量的函数的对角矩阵。

V(q)是N个颗粒之间的相互作用势；它是颗粒的位置的函数并将被考虑为独立与力矩。

在3维空间中，例如，在坐标系(X,Y,Z)中，每个颗粒a_i的力矩被写为(p_i,x,p_i,y,p_i,z)并且每个颗粒a_i的位置被写为(q_i,x,q_i,y,q_i,z)， i＝1到N。

通常，在现有技术中使用的矩阵M^-1是3N*3N的对角矩阵，对于 i＝1到N，M[3i-2,3i-2]＝M[3i-1,3i-1]＝M[3i,3i]＝m_i，m_i是颗粒a_i的质量。

这就是汉密尔顿算符的通常定义，以下将被称为标准汉密尔顿算符。

根据本发明，称作自适应相对汉密尔顿H_AR的汉密尔顿算符是以下：

其中，Φ(p,q)是称作自适应相对逆质量矩阵的3N×3N块对角矩阵，替代M^-1并依赖于矢量p并可选地依赖于矢量q。

更准确地，Φ(p,q)的值是在每次模拟中确定的，函数Φ被选择为能够把颗粒的至少一个子集限制为在模拟的一个或某些步长全部跟随同一运动。构成子集的颗粒被通过条件验证来确定，例如旨在验证在当前模拟步长之前，颗粒被接近的运动所激励。

这些设置导致出现颗粒的一个或多个子集，在所考虑的模拟步长期间子集中的颗粒被同一平动位移激励(即，在其中颗粒具有相同位移的时间范围期间连接这些颗粒中的任两个的矢量保持恒定)。

针对所考虑的模拟步长，这样子集的颗粒可以因此被相似对待，如构成唯一刚性体。

仅取决于颗粒相对位置的信息，例如颗粒之间的相互作用力，因此在同一个刚性子集内在所考虑的模拟步长结束处不需要被更新。

当子集的颗粒的运动不再受自适应逆质量矩阵限制时，每个颗粒采取自己专有的运动。

在本发明的实施方式中，在两类行为(限制运动/自由运动)之间施加连续过渡。

根据该自适应汉密尔顿算符H_AR，推导出限定是的自适应运动方程，是是矢量p和q相对于时间t的导数。

例如，在此处举例考虑的NVE(具有恒定颗粒数、体积和能量的组E)组中实施模拟时，汉密尔顿算符的值(根据本发明或标准是自适应的)随时间是恒定的，并且自适应运动方程是：

下面的示例与NVE组相关，但是本发明当然可以应用于其它组，例如NVT组(具有恒定颗粒数、体积和温度的组E)。

可以使用多个不同方案来定义构成刚性子集的颗粒，而无需测试集合的颗粒之间的组合的可能性。

第一方案是对刚性体的预定对施加约束，这些刚性体可被认为是要合并成同一刚性体的候选。例如，注意A₁、…A_n，即一组编索引的刚性体，可以把索引组织成二叉树，以使得当具体合并条件满足时A₁仅与A₂合并成刚性体(A₁+A₂)，A₃仅与A₄合并等，并且在二叉树的上层级，(A₁+A₂)仅与(A₃+A₄)合并成等。这样的树在最高层级包括称作R 结点的根结点，以及在最低层级包括与各个考虑的颗粒对应的叶结点。除了叶结点外，树的所有结点称为内部结点。

另一方案是对刚性体的可能包括多于两个刚性体的预定子集施加约束。因此可以把索引组织为n叉树以使得在同一时间A₁仅与A₂、 A₃…、A_n合并等。

另一方案是考虑从邻域的研究得出的图，图的结点是颗粒或刚性体，并且两个结点之间的边表示由边连接的两个结点之间的距离小于距离阈值。

下面详述的实施方式是基于如关于第一方案描述的二叉树结构；二叉树结构的优点是准许针对层次结构中的内部结点的判决度量的增量且递归的更新，以及例如针对部分能量和力的增量更新。

在这样的树中，每个叶结点表示集合中的单个颗粒，每个中间结点表示由所述中间结点的叶结点表示的颗粒组成的颗粒的子集。称体 C为包括结点C和结点C的直接或间接附属物的集合。

根据本发明，如果体C，即树的内部结点C，由两个体A和B 组成(即，结点C包括两个子结点A和B)，则使用以下递归公式定义与该结点C对应的自适应逆质量矩阵Φ_C：

其中：

－如果结点A和B不是叶结点，则对应于结点A的自适应逆质量矩阵Φ_A和对应于结点B的自适应逆质量矩阵Φ_B本身是根据公式3 定义的；在叶结点A和B的情况下，矩阵Φ_A和Φ_B分别等于分别由结点A和B表示的颗粒乘以3维单位矩阵的逆质量；

－矩阵E是由n_c×n_c个等于3维单位矩阵的3×3块形成的尺寸为3n_c×3n_c的矩阵，n_c等于结点C的叶节点的数量(即，结点C的直接子结点中或这些子结点的附属结点中的叶节点的集合)；

－m_c是由结点C的这些叶结点n_c表示的颗粒的质量之和；

－ρ_C是体C的限制函数。

因此，可考虑把ρ_C，即体C的限制函数，把两个刚性体A和B 合并成一个刚性体C(并反之把由两个刚性结点A和B组成的刚性体 C分解)。

更确切地，当函数ρ_C取1的值时，体C或内部结点C被认为是刚性的(即结点C中包括的所有叶结点或结点C的附属结点被同一平移运动激励)；否则被认为是非刚性的。

所有叶结点根据定义被认为是刚性的，并且在模拟持续的所有时间上是这样，因为它们每个都由单个颗粒构成。

当函数ρ_C取0的值时，内部结点C被认为是自由的。

在本发明的实施方式中，函数ρ_C是根据以下公式递归地定义的，以便使在0和1之间的转换平滑：

公式(4)

其中μ(ε)∈[0，1]是ε的二次可导函数。

公式(5)

其中阈值ε^f和ε^r被定义为使得ε^f＞ε^r并且s是ε的二次可导函数。

μ的二次可导性质允许保持颗粒的集合E的模拟稳定性。

例如，针对S(ε)可能的形式是5阶插值函数，例如等于 -6η⁵+15η⁴-10η³+1，其中δ＝ε^f-ε^r表示在体的刚性行为和自由行为之间的过渡区域的宽度。

在本发明的实施方式中，ε_C被独立于两个体A和B的力矩的关系进行选择，其中：

其中：A和B是结点C的子结点；

矢量是作为与结点C对应的体中的叶结点的所有颗粒的力矩之和：并且是矢量和之和：

该发明可以被针对k叉树进行概括。在该情况下，对于具有k个子结点A_i的母结点，1≤i≤k，惯性逆矩阵是：

并且

在该公式和公式6中的系数1/2可以被别的系数代替，例如1，这在阈值ε^f和ε^r被同样修改(即，乘以等于其它常数的2倍的系数) 的情况下不改变模拟。

本发明还可以不同地概括：可以针对不同的内部结点选择不同的限制函数ρ_C。例如，直到树的某个级别，ρ_C被根据公式4和5定义，因此，结点可以集中为整体或分离，并且在大于该级别的树的级别中，ρC 总等于零。在具体实施方式中，可以定义持续激活的系统E的子集(该子集的颗粒既不彼此合并，也不与系统的其它部分合并，因此这些颗粒在整个模拟中总是遵循不受限制的自由的运动，本发明的这样的实施方式可以例如应用于溶剂中的聚合物)：该子集中的所有颗粒应该放置在单独的子树中并且针对该子树的所有内部结点值ρC应该固定为 0并不再更新。在所有情况下，函数ρ_C应该在模拟开始之前针对所有内部结点被定义并且不应在模拟期间改变以便具有稳定的模拟。

根据本发明，选择在相对自适应汉密尔顿算符H_AR的公式(1) 中使用的自适应相对逆质量矩阵Φ，等于Φ_R，其为相对于表示颗粒的集合E的树的根结点R而言根据递归公式3定义的矩阵：Φ(p，q)＝Φ(p)＝Φ_R。

因此具有汉密尔顿算符H_AR，如其是可分的，这是因为Φ_R仅取决于颗粒的力矩(而不是位置)。

由公式2定义的运动方程因此表示如下：

还从公式3清楚可见如果结点C是自由的(ρ_C＝0)，则矩阵Φ_C由对角线上的对应于子结点A和子结点B的两个块组成。

在根据本发明的模拟的第一示例中，考虑具有包括通过刚度为1 的弹簧彼此连接的质量为1的两个颗粒P1、P2的尺寸的系统。对应于系统的二叉树是根结点C，其两个子结点A和B对应于这两个颗粒。在给定初始力矩的情况下，两个颗粒在空间中随时间振荡。

根据现有技术基于汉密尔顿算符H的模拟针对这些颗粒计算出的在空间中的轨迹d在图1中被根据模拟时间表示。

根据本发明基于自适应相对汉密尔顿算符H_AR的模拟针对这些颗粒计算出的轨迹d在图2中被针对δ和ε^f的不同值根据模拟时间表示。因此如在这些轨迹上表示那样，在某些模拟时刻(有时候从头开始)，两个颗粒的力矩根据公式6变得接近。这随着弹簧几乎被压缩并释放而发生。根据本发明，根结点C因此变得刚性。结果，在一般情况下，刚性结点C的子结点遵循同一平移运动，在当前情况下，由于力矩的保持，所以颗粒停止。2个颗粒的轨迹因此变得平行。颗粒的力矩继续演变，然而结果是在给定的模拟步长，结点C的刚性条件不再满足并且两个颗粒采取自己的运动。这因此在模拟过程中周期性地产生。

如轨迹所示，针对δ值增大的ε^f的同一值，轨迹的过渡区域更大。针对ε^f值增大的δ的同一值，刚性区域更长并更平。

因此根据本发明，被定义为等于Φ_R的矩阵Φ指定如何并何时在模拟过程中激活或禁止颗粒子集的相关位置自由度。

在本发明的实施方式中，图3中表示的信息设备1被用于根据上述的本发明原理来实施N个颗粒的集合E的模拟。

该设备1包括计算机，计算机尤其包括：存储器2，适于存储下面相继描述的软件程序和计算出的参数值(全局相互作用力，局部相互作用力，相互作用势，位置，力矩等)；微处理器3，适于执行软件程序、尤其是下述程序P的指令；以及人机界面4，例如包括键盘和屏幕，分别用于输入用户的指令和向用户显示信息，例如图2所示的曲线。

在所考虑的本发明的实施方式中，存储器2包括基于与公式7对应的方程随时间的积分来模拟颗粒a_i的集合E的行为的程序P，i＝1 到N。

因此根据本发明的自适应相对汉密尔顿算符H_AR是可分的，所以可以使用若干个积分技术。

参照图4，程序P包括软件指令，软件指令当被微处理器3执行时适于实施预备步骤100并在程序P的与计算时刻h_n+1＝h₀+(n+1)h对应的第n+1次迭代时重复实施步骤101至104，整数n≥0，h是模拟时间步长。

在初始化预备步骤100中，针对系统E构建表示颗粒并且两两颗粒的刚性的可能性的二叉树，如上所述；针对颗粒而言力矩和位置的初始值是固定的。仅叶结点因此是刚性的。

在步骤101中，确定在第n+1次迭代时在颗粒之间施加的相互作用力。

在步骤102中，确定在第n+1次迭代时颗粒的力矩。

在步骤103中，确定在第n+1次迭代时与树的内部结点有关的度量。这尤其包括ρ_C和ε_C。

在步骤104中，确定在第n+1次迭代时颗粒的位置。

根据前次迭代时计算出的值来实施用于当前迭代的力、力矩、测度、位置的值的所有这些更新步骤。

下面详述这些步骤101至104的实施。

在步骤101中，如果相互作用力仅取决于与颗粒有关的位置，则相对于现有技术，通过考虑根据本发明的刚性，可以使相互作用力的值的更新加速。

因此根据本发明，不需要重新计算在前次模拟迭代时(在ε_C之前确定的力)被识别为刚性(按照公式6根据ε_C的值)的颗粒的子集内部的力。为了通过利用根据本发明的刚性得到的该优点来有效地更新力，可以使用力的递增更新算法。

该算法首先基于包围盒层次(每个二叉树结点一个盒)，例如基于AABB(轴对齐包围体)层次，盒包括子系统并且具有其中边平行于坐标轴的平行六面体形式，参照例如Grudinin,S.和Redon,S，“Practical modeling of molecular systems with symmetries(对称性分子系统的实际建模)”,Journal of Computational Chemistry 31,9 (2010),pp.1799-1814)。

力的这种更新算法第二要使用局部力的表。与每个内部结点C 对应的局部力的表包括3D矢量，其数量等于结点C的直接或间接从属的叶结点的数量。该表的每个矢量是来自由C的其它从属叶结点表示的所有其它颗粒作用于由这些叶结点之一表示的颗粒上的相互作用力。根结点R的局部力的表的每个矢量对应于作用于该颗粒上的全局力。与叶结点相对应的局部力的表是零矢量。

每次迭代中力的更新包括三个主要子步骤(势能可以同时更新)：

在步骤101_1中，例如基于树的向上(即从叶结点向根结点)的步长来更新包围盒层次。对于每个叶结点，更新的AABB盒与由叶结点表示的颗粒的所确定的最后位置相一致。每个内部结点的盒被例如在包括其两个子结点的盒的同时来重新计算。

因此树的结构是不变的，包围盒被根据此时已改变的颗粒的位置来更新。

在步骤101_2中，通过使用AABB的这种更新层次，针对树的内部结点来递归地构建相互作用列表。对于具有两个子结点A和B的每个结点C，如此构建的相互作用列表包括所有颗粒对，如对的颗粒之一属于体A(即由结点A的直接或间接从属叶结点表示)而另一个属于体B。该相互作用列表可以针对每个内部结点构建，如在文献R. Rossi,M.Isorce,S.Morin,J.Flocard,K.Arumugam,S.Crouzy,M. Vivaudou和S.Redon,"Adaptive torsion-anglequasi-statics:a general simulation method with applications to proteinstructure analysis and design",Bioinformatics 200723(13):i408-i417中针对其它类型的包围盒所描述那样。

当第一次迭代时，相互作用列表被针对所有内部结点构建。位于所有其它迭代，在刚性结点中，这些列表不改变并且因此可以不重新计算。

在步骤101_3中，实现局部力的表的递增更新。当结点被识别为刚性结点时，在结点的颗粒之间的力不改变。因此，仅非刚性结点的局部力应该被更新。

其可以被如后面那样从叶结点向根结点递归地更新。

对于每个非刚性结点C，在与结点C有关的局部力的表的第一半中，与结点A有关的局部力的表的元件被复制并且在与结点C有关的局部力的表的第二半中，与结点B有关的局部力的表的元件被复制。然后，针对用于结点C的在子步骤101_2中建立的相互作用表的每一对，该表中的每对的两个颗粒之间的相互作用力被计算并且添加到局部力的表的对应力。最后，在根结点R的局部力的表中存储考虑所有颗粒的力。

例如在树的单次遍历中，步骤101_2和101_3可以组合，以在与所述结点C对应的相互作用列表一旦可用时执行与结点C相对应的局部力的表的更新。

仅取决于颗粒的相对位置的其它信息可以类似地更新。

颗粒的力矩的更新的步骤102是通过使用常规方法实现的。

因此，对于每个颗粒a_i，通过为p_i,n+1填写颗粒a_i的力矩p_i的值，当当前迭代n+1时，该值通过以下等式获得：

p_i，n+1＝p_i，n+f_i，n+1h，

h是模拟的时间步长，f_i，n+1是作用于颗粒a_i上的全局相互作用力， i＝1到N，这是由于系统E中的所有其它颗粒施加的相互作用(出现在如当前迭代中刚确定的根结点的局部力的表中)引起的。

在每个时间步长，基于等式的颗粒位置的“自然”更新的复杂度为O(N³)，这是由于项该项是矩阵相对于矢量的倒数，因此是三维对象，对于1到N：

在所考虑的实施例中，对于具有两个子结点A和B的每个内部结点C，通过从数据结构Q_C、R_C和S_C的叶结点进行递归更新，复杂度可以减小为O(N.logN)。

注意：

其中p_A是作为结点A的叶节点的所有颗粒的力矩矢量，p_B是作为结点B的叶节点的所有颗粒的力矩矢量，项Q_C＝Φ_C(p_C)p_C对应于给出的等式的第一项：

其中{x}是维度等于n_C的矢量，其n_C个组成部分的每一个等于元素x。

类似地，对于结点C考虑项

给出的等式的最后一项可以通过使用等式9来递归更新。

事实上：

函数(·)表示标量函数。

因此，步骤103包括针对包括两个子结点A和B的每个结点C 的结构Q_C、R_C和S_C的更新。

应注意步骤100还包括这些结构的初始化步骤，其中对于每个内部结点C：ρ_C＝0并且Q_C、R_C和S_C还是0，并且对于表示颗粒的每个子结点F：ρ_F＝1；Q_F＝p_F/m_F；S_F＝0，R_F＝0。

在步骤103中，度量Q_C、R_C和S_C以及ε_C和ρ_C的更新相对于结点C的更新被从根结点沿着结点的架构递归地执行：

－如果结点C是内部的，具有两个子结点A和B：

－更新结点A的度量；

－更新结点B的度量；

－借助等式6来更新ε_C；

－更新ρ_C：ρ_C＝μ(ε_C)ρ_Aρ_B；

－借助等式9来更新Q_C；

－借助等式10来更新R_C；

－借助等式11来更新S_C；

－如果结点C是表示颗粒a_F的叶结点：Q_F＝p_F/m_F。

为了优化计算，矢量Q_C、R_C和S_C的更新可以仅针对ρ_C＞0的结点来进行，这是因为这些度量可以不用于自由结点，所以不需要更新它们。

在一般情况下，这些度量可以以不遍历树的方式来更新，只要在刚性结点在顶部的情况下在子树中更新即可。

一旦确定了这些度量，颗粒的位置就在步骤104中更新。

步骤104的实施可以在结点C的位置的更新之后进行，每个内部结点C被考虑作为包括两个子结点A和B：

－如果ρC＝0(结点C是自由结点)：

－更新结点A的位置；

－更新结点B的位置；

－否则，对于结点C的下属每个叶结点k，无论直接还是间接：

q_i，n+1＝q_i，n+(Q_C[k]+0.5S_C[k])h，

q_i，n+1是颗粒a_i的位置的值，针对被考虑的第n+1次迭代，并且h 是模拟时间，Q_C[k]是矢量Q_C的第k个元素，以及S_C[k]是矢量S_C的第k 个元素。

颗粒的位置因此就确定了。

在该实施方式中，对于所有颗粒以及树的所有内部结点，在颗粒 i的标识符及其在结点C的子结点中的序号之间的一一对应关系必需建立，例如在模拟初始化期间建立。

在位置的该实施方式中，ρ_C＞0的结点的度量不被使用并且因此可以不计算。如果在步骤103中更新所有度量，则还可以如下确定所有颗粒a_i(1≤a_i≤N)的位置：q_i，n+1＝q_i，n+(Q_R[k]+0.5S_R[k])h，其中索引 R对应于根结点。

对于树的每个内部结点C，存储与结点C的直接或间接从属的叶结点的数量相等的4个长度矢量：Q_C、R_C和S_C以及局部力的表。位置更新的空间复杂度因此是O(N.logN)。时间复杂度也是O(N.logN)，这是因为这些结构应在每个时间步长更新，并且Q_C总是对所有结点(其中包括叶结点)更新。

在一般情况下，位置可以不以遍历树的方式更新。

然后，如果未达到模拟的最大持续时间，则实施程序P的新的迭代。

作为示例，从把颗粒高速发送到初始不动的系统触发的震动开始使用兰纳-琼斯势(Lennard-Jones potential)(E_m/k_B＝120开尔文，其中E_m是最小能量，平衡距离S＝3.4埃，截止距离8埃，势在7.5和 8之间被平滑截短)来执行N＝5930个颗粒a_i(i＝1到N)的NVE集合的2个演化维度的4次模拟，每个颗粒具有1g/mol的质量：基于标准汉密尔顿算符的一个参考模拟和三个自适应模拟，也就是说使用根据如上所述的本发明的方法的自适应相对汉密尔顿算符(时间步长尺寸 0.0488飞秒(fs)，7000时间步长，总模拟时间342fs)。

对于使用自适应相对汉密尔顿算符的每次模拟，相对于标准模拟的波动的平方根，称为RMSD，被给出作为最大颗粒位移误差Δq_max。

其中，q_i是颗粒a_i在自适应模拟的上一步的坐标矢量，是该同一颗粒在参考模拟的上一步的坐标矢量。

例如，对于其中ε^r＝0，01kcal/mol和(δ＝1kcal/mol) 的自适应模拟，获得了用于执行该自适应模拟所需的计算时间相对于与参考模拟有关的计算时间的值等于1.5的加速因子，RMSD＝8.6S 并且Δq_max＝89.4S，其中S是所用的兰纳-琼斯势中的平衡距离。

对于其中和(δ＝1kcal/mol)的自适应模拟，获得了用于执行该自适应模拟所需的计算时间相对于与参考模拟有关的计算时间的值等于1.55的加速因子，RMSD＝8.7S并且Δq_max＝89.4S。

对于其中和(δ＝1kcal/mol)的自适应模拟，获得了用于执行该自适应模拟所需的计算时间相对于与参考模拟有关的计算时间的值等于1.6的加速因子，RMSD＝12.3S并且Δq_max＝89.4S。

因此，根据本发明的方法允许使计算加速，而潜在较小的改变性能。

结果取决于针对颗粒的系统所选择的树状表示。对于上述的自适应模拟，通过把系统一分为二来从上至下构建了二叉树。

在本发明的所考虑的实施方式中，已经考虑了把子集两两融合为刚性集合，在其它实施方式中，融合成一个刚性集合的子集的数量被取为大于2。

在刚性和自由的两个状态之间的过渡区域可以被取为更大或更小的宽度。

Claims (8)

1.一种对元件的系统进行模拟的方法，根据该方法所述元件的行为是基于元件的系统的汉密尔顿算符H在相继的模拟步长上确定的以使得其中p是表示元件的力矩的矢量，q是表示元件的位置的矢量，M^-1是随元件的质量而变化的矩阵，以及V是系统的势能，

其特征在于所述方法包括以下步骤：

-把元件的系统表示为包括叶节点和内部结点的k叉树，所述叶节点每个都表示相应的元件，内部结点中包括根结点，

-当所确定的条件被验证时，对于当前模拟的步长，通过把矩阵M^-1定义为等于Φ_R，对树的给定结点的从属的至少一些元件施加同样的平移运动，其中R是树的根结点并且对于包括k个子结点A₁、A₂，...，A_k的树的任何结点A，有以下递归等式，根据该等式：

矩阵E是由各自尺寸为d＊d的n_A＊n_A个块形成的尺寸为dn_A＊dn_A的矩阵，块等于尺寸为d的单位矩阵，nA等于结点A的从属的元件的数量，d是其中颗粒演变的空间的尺度，m_A是从属于结点A的这n_A个元件的质量之和，ρ_A是介于0和1之间的结点A的约束函数，该函数的值当A是叶节点时或者当同样的运动被施加到从属于给定结点的元件时等于1；在叶节点A_i的情况下，矩阵等于由结点A_i表示的颗粒乘以尺度为d的单位矩阵所得的逆矩阵。

2.按照权利要求1所述的对元件的系统进行模拟的方法，其中对于当前模拟步长，对树的给定结点的从属的至少一些元件根据所述元件的力矩函数的值来施加同样的平移运动。

3.按照权利要求1或2所述的对元件的系统进行模拟的方法，其中，对于当前模拟步长，根据由所取的值对树的给定结点的从属的至少一些元件施加同样的平移运动，其中是结点A_i的元件的质量之和，是子结点A_i的元件的力矩的矢量和，||·||是矢量的范数，C是正常数。

4.按照权利要求3所述的对元件的系统进行模拟的方法，其中，对于针对内部结点A的当前模拟步长：

-如果ε_A小于第一阈值，则把ρ_A定为1以便对从属于结点A的元件施加同样的平移运动；

-如果ε_A大于比第一阈值大的第二阈值，把ρ_A定为0。

5.按照权利要求4所述的对元件的系统进行模拟的方法，其中，对于当前模拟步长，如果ε_A介于第一阈值和第二阈值之间，则把ρ_A定为等于由变量ε_A的5阶插值函数所取的值。

6.按照权利要求1或2所述的对元件的系统进行模拟的方法，包括在基于所述汉密尔顿算符的相继模拟的时刻确定至少一个信息的值的步骤，所述确定至少一个信息的值的步骤利用了以下事实：对于当前模拟的时刻与已经施加了同样的平移运动的元件相关的信息的值取决于所述元件的相对位置并且因此不变化。

7.按照权利要求1或2所述的对元件的系统进行模拟的方法，其中，与所述元件相关的信息包括所述元件的势能和/或应用到所述元件的相互作用力。

8.一种计算机可读存储介质，上面存储有对元件的系统进行模拟的计算机程序，该计算机程序包括能由处理器存取并且当由处理器执行时用于实施根据权利要求1至7之一所述的方法的步骤的软件指令。