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

Abstract

一种对元件的系统进行模拟的方法，根据该方法所述元件的行为是基于元件的系统的哈密顿算符H确定的以使得(公式I)，其中p是表示元件的力矩的矢量，q是表示元件的位置的矢量，M^‑1是随元件的质量而变化的对角矩阵，以及V是系统的势能，该方法包括以下步骤：根据该步骤，当力矩矢量p采取与至少一个元件有关的某个预定值时，零值被分配给矩阵M^‑1的与该至少一个元件有关的至少一个对角项。

Description

用于模拟一组元件的方法

技术领域

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

背景技术

在此情况下势能例如等于元件之间的可从其推导相互作用力的相互作用势或是相互作用势V(q)的函数，q是表示元件的位置的矢量 (在此情况下，更一般而言相互作用势还可以依赖于元件的力矩)。

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

模拟一组元件的常规方式是考虑该组元件的哈密顿算符并从中推导力矩方程。

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

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

发明内容

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

为此，根据第一方面，本发明提出一种对上述类型的一组元件进行模拟的方法，所述方法是由计算机执行的并且其特征在于所述方法包括以下步骤：根据该步骤，当力矩矢量p采取与至少一个元件有关的某个预定值时，零值被分配给矩阵M^-1的与该至少一个元件有关的至少一个对角项。

本发明使得可以减少计算量，因此减少确定元件的势能、相互作用势、相互作用力、位置和/或所需力矩所需的计算时间。

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

-其中所述方法包括以下步骤：对于所述元件中的至少一个元件，如果表示所述元件的动能的参数具有低于第一严格为正的阈值，则把零值分配给矩阵M^-1的与该至少一个元件有关的至少一个对角项；

-当所述至少一个元件的动能大于第二严格为正的阈值时随元件的质量而变化的矩阵M^-1的对角项被分配最大值；

-如果包括所述至少一个元件的力矩和所述至少一个元件的位置的数偶采用某些预定值则把零值分配给与该至少一个元件有关的矩阵M^-1的至少一个对角项；

-该方法包括基于所述哈密顿算符在相继的模拟时刻确定至少一个信息的值的步骤，所述步骤利用与在前一模拟时刻为矩阵M^-1的对角项分配零值所针对的元件k元组有关的信息的值因此在至少所述前一模拟时刻和当前模拟时刻之间未改变的事实，k是大于等于2的整数，并且在零值未被分配给矩阵的与元件k元组中的每个元件有关的对角项时通过执行以下步骤来计算在当前模拟时刻与构成k元组一部分的给定元件有关的信息的值：

-通过从在前一模拟时刻确定并且与所述给定元件相关的信息的值中减去在前一模拟时刻至少与所述元件k元组相关联并且与所述给定元件有关的所述信息的值来计算与所述给定元件相关的所述信息的工作值；以及/或者

-向所述工作值添加在当前模拟时刻确定的至少与所述元件k元组相关联并且与所述给定元件有关的信息的值。

-在当前计算时刻，由低于给定阈值的距离分隔的元件对的当前列表被在当前模拟时刻准备并且与在前一模拟时刻准备的由低于给定阈值的距离分隔的元件对的前一列表相比较，

并且在当前模拟时刻与给定元件有关的信息的值被基于包括所述给定元件的各元件对通过执行以下步骤来计算：

-如果所考虑的对仅存在于前一列表中或者如果把所述给定元件连接到对中的另一元件的矢量已经在前一模拟时刻和当前模拟时刻之间变化，则通过从在前一模拟时刻确定并且与所述给定元件相关的信息的值中减去与所述给定元件相关的信息的与所考虑的对中的另一元件相关联的值，来计算工作值；

-如果所考虑的对仅存在于当前列表中或者如果把所述给定元件连接到对中的另一元件的矢量已经在前一模拟时刻和当前模拟时刻之间变化，则通过向所述工作值添加与所考虑的对中的另一元件相关联并且与所述给定元件相关的信息的值来确定在当前模拟时刻与所述给定元件相关的信息的值。

-在当前计算时刻，满足某些条件的元件k元组的当前列表被在当前模拟时刻准备并且与在前一模拟时刻满足所述条件的元件k元组的前一列表相比较，k是大于等于2的整数，

并且在当前模拟时刻与一个元件有关的信息的值被基于包括所述元件的k元组对通过执行以下步骤来计算：

-当所述k元组仅存在于前一列表中时或者当与所述k元组相关联并且与所述元件相关的信息的值已经在前一模拟时刻和当前模拟时刻之间变化时，通过从在前一模拟时刻确定并且与所述元件相关的信息的值中减去在前一模拟时刻与所述k元组相关联并且与所述元件相关的信息的值，来计算临时值；

-当所述k元组仅存在于当前列表中时或者当与所述k元组相关联的信息已经在前一模拟时刻和当前模拟时刻之间变化时(例如当k元组中的k个元件的相对位置已变化时)，则通过向所述临时值添加在当前模拟时刻与所述k元组相关联并且与所述元件相关的信息的值来确定在当前模拟时刻与所述元件相关的信息的值。

-元件的位置空间被划分为单元并且在前一模拟时刻和当前模拟时刻中的每一个，每个元件根据在模拟时刻确定的位置坐标与归属单元相关联，以及对于第一元件执行以下步骤，使得在当前模拟时刻与所述第一元件相关的矩阵M^-1的项还未被分配零值：

-确定在前一模拟时刻第一元件的归属单元；

-对于每个第一元件，在所述归属单元中或所述归属单元的给定邻域中的单元中，确定在前一模拟时刻与所述第一元件相距小于给定阈值的第二元件；通过从在前一模拟时刻确定并且与所述第一元件相关的信息的值中减去与所述第二元件相关联并且与所述第一元件相关的所述信息的值来计算工作值；

-确定在当前模拟时刻的第一元件的新归属单元；

-对于每个第一元件，在所述新归属单元中或所述新归属单元的给定邻域中的单元中，确定在当前模拟时刻与所述第一元件相距小于给定阈值的第三元件；

-通过向工作值添加与所述第三元件相关联并且与第一元件相关的信息的值来确定在当前模拟时刻与所述第一元件相关的信息的值。

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

-该方法包括在某些模拟时刻确定信息I的步骤，所述步骤优选地利用在某些模拟时刻零值被分配给矩阵M^-1的某些对角项的事实。

-该方法包括在某些模拟时刻确定信息I的步骤，所述步骤优选地利用如下事实：该信息I不改变并且当信息I在前一模拟时刻已确定并且在至少所述前一模拟时刻和当前模拟时刻之间零值被分配给矩阵M^-1的对角项的对应集合时不必再次确定(“对角项的对应集合”是指具有对信息I的值的影响的对角项，即当这些项为零时信息I不变化)；

-该方法包括在某些模拟时刻确定势能或相互作用力的步骤，所述步骤优选地利用在某些模拟时刻零值已被分配给矩阵M^-1的至少一个对角元素的事实；

根据第二方面，本发明提出一种用于模拟元件的系统的计算机程序，包括用于在该程序被计算装置执行时执行根据权利要求1至12 中任一项的方法的步骤的软件指令。

附图说明

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

-图1示出实施本发明的实施例的设备；

-图2示出作为x轴上的ε_i的函数的函数ρ_i(ε_i)在y轴上的演变；

-图3是本发明的实施例中的方法步骤的流程图；

-图4示出步骤103的一个实施例；

-图5示出步骤103的另一实施例；

-图6示出用常数哈密顿算符的与相空间(p，q)中的固定点相关联的颗粒的轨迹的模拟。

具体实施方式

考虑到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。

矢量p和q因此写为：和

通常，在现有技术中使用的矩阵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_A的哈密顿是以下：

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

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

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

根据本发明，矩阵的与颗粒a_i相关的项(i＝1至N)如下：

其中m_i是颗粒a_i的质量。

的值被称为Φ_i(p_i，q_i)。

根据本发明，ρ_i(q_i，p_i)∈[0，1]并且是二次可导函数，这被使用以使得颗粒的位置在某段时间上恒定。

当ρ_i(q_i，p_i)∈[0，1](也就是说颗粒a_i的位置不固定)时，并且颗粒a_i遵循与系统E的标准哈密顿算符H对应的通用动态法则。

当ρ_i(q_i，p_i)＝0(也就是说颗粒a_i的位置不固定)时，并且颗粒a_i遵循与系统E的标准哈密顿算符H对应的通用动态法则。

当ρ_i(q_i，p_i)＝1(也就是说颗粒a_i的位置固定)时，并且颗粒a_i不移动，而无论作用于其上的力如何(其质量被考虑无限大)。

当ρ_i(q_i，p_i)∈]0，1[时，颗粒a_i在这两个行为之间进行平滑过渡。

ρ_i的二次可导性质允许保持颗粒的组E的稳定性。

在本发明的一个实施例中，定义：

其中，s(ε_i)是的p_i函数并且可选地是q_i的函数，并且是二次可导的。

在图2中，在一个实施例中ρ_i(ε_i)采取的值被示于y轴上，作为由示于x轴上的ε_i采用的值的函数，其中并且

例如，用于s(ε_i)的可能形式是-6η⁵+15η⁴-10η³+1，其中并且

在以下考虑的实施例中，涉及颗粒a_i的函数ε_i被选择为等于颗粒 a_i的动能，即

本发明因此包括通过向颗粒分配几乎无限大的质量来固定颗粒，此时颗粒的动能低于某个值，颗粒的运动量不固定。

函数ρ_i是包括力矩作为变量的函数(在特别情况下，举例来说，其因此不依赖于位置q_i：ρ_i(q_i，p_i)＝ρ_i(p_i))。

在其它实施例中，函数ρ_i可以是依赖于颗粒a_i的力矩(可选地依赖于位置)而当然不是依赖于动能的函数。

在一些实施例中，固定其力矩(包括力矩和位置的数偶)采用预定值(离散值或值的范围)的颗粒。

力矩(1)的自适应方程因此变为：

其中，ρ(p)是表示ρ_i(q_i，p_i)的3N＊3N对角矩阵，对于i＝1到N：

ρ(p)[3i-2，3i-2]＝ρ(p)[3i-1，3i-1]＝ρ(p)[3i，3i]＝ρ_i(q_i，p_i)，对于i＝1到N。

如上所述，当ρ_i(q_i，p_i)＝0(也就是说位置不固定)时，并且颗粒a_i遵循与系统E的标准哈密顿算符H对应的通用动态法则。

当ρ_i(q_i，p_i)＝1(也就是说恒定颗粒的位置固定)时，Φ_i(p_i，q_i)＝0，因此，具有0值(事实上，ρ_i是等于1的常数，项的值为0)。在一个解释中，质量被考虑无限大，颗粒a_i被考虑为固定。

当ρ_i(q_i，p_i)∈]0，1[时，颗粒a_i在这两个行为之间进行平滑过渡。

因此，根据本发明，矩阵Φ(p，q)指定在模拟期间如何、并且何时在一个或更多个颗粒的位置方面的自由度被激活或禁止。

为了给出解释在本发明情况下的系统行为的示例，图6示出包括质量为1、附加到在固定点处刚度为1的弹簧的颗粒的系统在相空间 (p，q)中的一个维度的轨迹模拟。在此情况下，N＝1。自适应哈密顿算符的等值线被针对和示出。

特别地曲线C1、C2、C3、C4每个都对应于自适应哈密顿算符的各对应值。例如，对应于哈密顿算符的曲线C1等于1。对应于常数的圆D等于标准哈密顿算符的1，即非自适应。

颗粒被固定的相空间的区域位于虚线B2和B3之间(其对应于包括于[-1，1]之中的力矩值)。

颗粒自由并且遵循根据标准哈密顿算符的轨迹的相空间的区域被发现位于线B1之上并且在线B4之下。

在线B1和B2之间并且在虚线B3和B4之间的相空间的区域对应于颗粒的自由和固定状态之间的过渡区域。

在颗粒固定的区域中在每个等值线C1、C2、C3、C4上，位置q 不变化，但是力矩p变化。

为了执行系统E的模拟，要执行公式(2)中指示的力矩方程随时间的数值积分。在此处举例考虑的NVE组中实施模拟的情况下，针对该数值积分例如使用分区的欧拉方法(参见例如″Geometric numerical integration：structure preserving algorithms forordinary differential equations″，Hairer E.，Lubich C.，Wanner G.；Volume 31；Springer Verlag 2006)。

根据该欧拉方法，如下形式的方程：

针对数值积分导致以下组方程，其中h是时间步长：

u_n+1＝u_n+a(u_n+1，v_n)h，

v_n+1＝v_n+b(u_n+1，v_n)h.

因此，根据该方法，公式(2)可以写作：

在本发明的一个实施例中，图1所示的计算设备1用于执行N 的颗粒的组E的模拟。

设备1包括计算机，计算机尤其具有：适于存储软件程序以及后文描述的相继计算出的参数值(矩阵Φ的系数值，全部、部分相互作用力，相互作用势，位置，力矩等)的存储器2、适于执行软件程序的指令和特别地后文描述的程序P的微处理器3，以及包括例如用于输入用户指令和用于显示用于用户的信息(例如图6中示出的那样的曲线)的键盘和屏幕的人机/界面4。

在所考虑的本发明的实施例中，存储器2包括模拟NVE类型的颗粒的组E的行为的程序P。

程序P包括软件指令，当程序在微处理器3上执行时软件指令适于执行参照图3的以后步骤。

在前一步骤100中，矩阵Φ的函数ρ_i(p_i，q_i)被预先针对每个颗粒 a_i定义。

在本案例中，以上定义的函数被选择并因此根据力矩的矢量p以及针对和固定的值来定义Φ(p)。

还对于每个颗粒a_i的力矩p_i，0、位置q_i，0和相互作用力f_i，0来确定与初始时刻h₀对应的初始值，i＝1到N。

然后在与计算时间时刻h_n+1＝h₀+(n+1)h相对应的程序P的n+1 次迭代中重复执行以后步骤，其中n是≥0的整数，h是模拟时间步长。

下文中：

f_ij，n+1表示在计算步骤h_n+1中由颗粒a_i对颗粒a_j施加的相互作用力(等于-f_ji，n+1)；

f_ij，n+1表示在计算时刻h_n+1对颗粒a_i施加的全部相互作用力，这是由于由系统E中的其它颗粒施加的相互作用引起的，因此等于

p_i，n+1将表示颗粒a_i在计算时刻h_n+1处的力矩值；

q_i，n+1将表示颗粒a_i在计算时刻h_n+1处的位置值；

ρ_i，n+1将表示由函数ρ_i在计算时刻h_n+1处采取的值(如上可见，它是动能的值的函数)。

步骤101、101b、102、103旨在确定与每个颗粒a_i相关的力矩、位置和总相互作用力各自的更新值。

在步骤101中，每个颗粒a_i的力矩的当前值p_i，n+1是根据公式(3)、已经如力矩p_i，n一样存储于存储器2中的在前一计算时刻确定的施加于颗粒a_i上的全部相互作用力f_i，n确定的

p_i，n+1＝p_i，n+f_i，nh。

在步骤101b中，根据力矩p_i，n+1的新值来重新计算值ρ_i，n+1：

在步骤102中，每个颗粒a_i的位置的当前值q_i，n+1现在根据公式 (3)来确定，i＝1到N：

在步骤103中，根据例如以下指出的两种方法之一来确定由所有其它颗粒引起的作用于颗粒a_i上的全部相互作用力的当前值f_i，n+1，至少i＝1到N。

该数值积分的结果是针对计算时刻h_n+1确定的每个颗粒a_i的位置 q_i和力矩p_i的更新值，其中i＝1到N。

颗粒a_i在时刻h_n+1的行为的其它特性参数的更新值可以被进一步计算出来，例如系统E的势能的当前值，颗粒速度之间的自相关性的值。

然后，如果模拟的最大持续时间还未达到，也就是说，如果n+1≤ n_max，则可以执行程序P的进一步重复。

可以使用在步骤103中用于计算相互作用力的更新值的各种技术，其有利地允许利用矩阵Φ的定义。

第一技术包括以下步骤：

在初始化步骤100中，准备颗粒对的当前列表，使得每对颗粒之间在初始化时的距离低于d0(当两个颗粒之间的距离大于d0时，这两个颗粒之间的相互作用被忽略)，并且当前列表中存在的每对的颗粒a_j对颗粒a_i的相互作用力f_ij，0进一步被根据它们之间相隔的距离以及根据被模拟的力场来评价，并且存储。

在列表中每一对还与也存储在存储器2中的一个元素e_ij，0有关联，该元素包括：该对中的两个颗粒a_i、a_j中每一个的标识符，从颗粒a_i连接两个颗粒的矢量的坐标，以及由颗粒a_j对颗粒a_i施加的相互作用力f_ij，0(等于-f_ji，0，f_ji，0等于由颗粒a_i对颗粒a_j施加的相互作用力)。

作用于颗粒a_i上的总相互作用力f_i，0等于由颗粒a_j对颗粒a_i施加的相互作用力f_ij，0之和，j＝1到N。

在步骤103的每次迭代中，然后执行参照图4的以下步骤。

令当前迭代为第n+1次。

在步骤103_a1中，把前次迭代中计算出的总相互作用力f_i，n的值分配给作用于每个颗粒a_i上的总相互作用力f_i，n+1作为初始值。

在步骤103_b1中，准备相互作用的颗粒对的当前列表L_a，n+1，也就是说这些是使得在计算时刻h_n+1所考虑的每对颗粒之间的距离低于阈值d0的颗粒对。

在列表L_a，n+1中每一对还与一个元素e_ij，n+1有关联，该元素包括：该对中的两个颗粒a_i、a_j中每一个的标识符，根据在时刻h_n+1的位置的从颗粒ai连接两个颗粒的矢量的坐标，以及由颗粒a_j对颗粒a_i施加的相互作用力f_ij，n+1(等于-f_ji，n+1，f_ji，n+1等于由颗粒a_i对颗粒a_j施加的相互作用力)，其值在该阶段未知。

然后，在步骤103_c1，当前对列表L_a，n+1与前一对列表L_a，n比较， L_a，n也就是说，前次迭代(即第n次迭代)中建立的列表。在列表L_a，n中每一对还与一个元素e_ij，n有关联，该元素包括：该对中的两个颗粒 a_i、a_j中每一个的标识符，根据在用于颗粒a_i、a_j的前次迭代中确定的位置从颗粒a_i连接两个颗粒的矢量的坐标，以及由颗粒a_j对颗粒 a_i施加的相互作用力f_ij，n(等于-f_ji，n，f_ji，n等于由颗粒a_i对颗粒a_j施加的相互作用力)。

如果一对仅存在于前一列表L_a，n中，这意味着该对的两个颗粒之间的迭代在第n和第n+1次迭代之间消失了。

对于仅在前一列表中存在的每对a_i、a_j，从当前确定的总作用力中拿去了作用于颗粒a_i的f_i，n+1或者作用于颗粒a_j上的f_j，n+1，在前次迭代中在步骤100中计算出的相互作用力f_ij，n或相互作用力f_ji，n＝-f_ji，n。

如果一对仅存在于当前列表L_a，n+1中，这意味着该对的两个颗粒之间的迭代在第n和第n+1次迭代之间出现了。

对于仅在当前列表中存在的每对a_i、a_j，因此尤其根据各自相应的位置计算由颗粒a_j作用于颗粒a_i的相互作用力f_ij，n+1；其被存储在存储器的元素e_ij，n+1中。相互作用力f_ij，n+1被添加到当前被确定的对颗粒a_i的总相互作用力f_i，n+1，并且相互作用力f_ji，n+1＝-f_ij，n+1被添加到当前被确定的对颗粒a_j的总相互作用力f_j，n+1。

如果一对被包括在当前列表L_a，n+1和前一列表L_a，n中，则比较连接两个颗粒a_i、a_j的矢量和如果不同，则尤其根据各自相应的位置来计算由颗粒a_j作用于颗粒a_i的相互作用力f_ij，n+1，并且相互作用力被存储在元素e_ij，n+1中。相互作用力f_ij，n+1被添加到作用于颗粒a_i的当前被确定的总相互作用力f_i，n+1。相互作用力f_ji，n+1＝-f_ij，n+1被添加到作用于颗粒a_j的f_j，n+1。此外，前次迭代中在步骤100中计算出的相互作用力f_ji，n或者相互作用力f_ji，n＝-f_ij，n被从作用于颗粒a_i上的当前被确定的总相互作用力f_i，n+1或者作用于颗粒a_j上f_i，n+1的中拿去了。

通过固定颗粒的位置，本发明生成增大数量的对，针对这些对而言，两个颗粒之间的矢量、并因此这两个颗粒之间的相互作用力保持不变。

在这样的情况下，步骤103提出的方法意味着，通过利用根据本发明的方法的特征，不必重新计算总相互作用力的所有分量。

总相互作用力的这种确定技术在计算量方面是最优的。然而，列表的构造需要时间。

用于执行步骤103的第二技术允许利用本发明提供的优点而无需使用当前迭代中的相互作用的颗粒对列表与前次迭代的比较，但是使用三维网格(如果颗粒的运动在三维空间中被考虑；如果颗粒在平面中移动的话，则二维网络可以足够)。

在初始化步骤100中，此外，通过考虑包含所有颗粒的平行六面体并把它细分为单元，例如其一边的尺寸大于或等于d0的立方单元，来创建初始网格。

每个颗粒a_i(i＝1到N)然后根据在初始化步骤的颗粒位置被分配给其所属于的单元。

然后，对于位于给定单元中的每个颗粒a_i，考虑该给定单元或与之相邻的单元(考虑最大26个单元)，在所述相邻的单元中，颗粒a_j与颗粒a_i的距离小于d0。对于位于与颗粒a_i的距离低于d0并使得 j＞i的那些颗粒a_j，计算颗粒a_j对颗粒a_i的相互作用力。该相互作用力等于-f_ji，0，f_ji，0等于由颗粒a_i对颗粒a_j施加的相互作用力。

应当注意在所描述的实施例中，所考虑的相邻单元是紧邻的单元，也就是说，与给定单元具有至少一条公共边的单元；在其它实施例中，所考虑的相邻单元是与给定单元紧邻的单元相距r个单元的那些单元。

以下步骤然后在程序P的每个n+1迭代中执行，其中n≥0，参照图5。

在步骤105_a2中，在前次迭代中计算的总相互作用力f_i，n的值被分配给作用于每个颗粒a_i上的总相互作用力f_i，n+1，作为初始值。

在步骤105b2中，对于ρ_i，n+1＜1的所有颗粒a_i(也就是说不被认为是固定的颗粒)，确定满足以下条件的颗粒a_j：

-在前次迭代(n)中，这些颗粒a_j位于在前次迭代n中颗粒a_i所位于的单元中，或者与该单元相邻(最多考虑26个单元)的单元中；和

-在前次迭代n中，这些颗粒a_j与颗粒a_i相距小于d0的距离；和

-这些颗粒a_j满足其索引j严格大于i或满足ρ_j，n+1＝1。

所考虑的网格的组成因此是与前次迭代(迭代n)中更新的位置对应的组成。

然后，对于这样确定的这些颗粒a_j中每一个，由颗粒a_j对颗粒 ai施加的相互作用力f_ij，n(并因此由颗粒a_i对颗粒a_j施加的相互作用力f_ji，n)根据在前一步骤n中分隔它们的距离来计算。

由颗粒a_j对颗粒a_i施加的相互作用力f_ij，n被从施加于a_i的总相互作用力中减去；类似地，由颗粒a_i对颗粒a_j施加的相互作用力f_ji，n被从施加于a_j的总相互作用力中减去：因此计算了f_i，n+1＝f_i，n+1-f_ij，n和f_j，n+1＝ f_j，n+1-f_ji，n＝f_j，n+1+f_ij，n。

在步骤103_c2中，通过确定符合ρ_i，n+1＜1的所有颗粒(即不被认为固定的颗粒)的当前所属单元、根据这些颗粒的与迭代n+1对应的位置q_i，n+1来更新网格的组成。

在步骤103_d2中，对于符合ρ_i，n+1＜1的所有颗粒(也就是说不被认为是固定的颗粒)，确定满足以下条件的颗粒a_j：

-在当前迭代(n+1)中，这些颗粒a_j位于在当前迭代中颗粒a_i所位于的单元中，或者与该单元相邻(最多考虑26个单元)的单元中；和

-在当前迭代(n+1)中，这些颗粒a_j与颗粒a_i相距小于d0的距离；和

-这些颗粒a_j满足其索引j严格大于i或满足ρ_j，n+1＝1。

此处所考虑的网格的组成因此是与当前迭代(迭代n+1)中更新的位置对应的组成。

然后，对于这样确定的这些颗粒a_j中每一个，由颗粒a_j对颗粒 a_i施加的相互作用力f_ij，n+1(并因此由颗粒a_i对颗粒a_j施加的相互作用力f_ji，n+1)根据在当前步骤n中分隔它们的距离来计算。

由颗粒a_j对颗粒a_i施加的该相互作用力f_ij，n+1被添加到施加于a_i的总相互作用力；类似地，由颗粒a_i对颗粒a_j施加的相互作用力f_ji，n+1被添加到施加于a_j的总相互作用力中：因此计算了f_i，n+1＝f_i，n+1-f_ij，n+1和 f_j，n+1＝f_j，n+1-f_ji，n+1＝f_j，n+1+f_ij，n+1。

如同第一技术，第二技术利用一些颗粒不固定不用重新计算固定颗粒之间的相互作用力的事实。其进行与以前位置对应的力的减法和对于新位置对应的力的加法。其不包括准备列表的冗长操作和执行每个列表的各对之间的比较。作为对比，要计算的两个颗粒之间的相互作用力的量大于第一技术中要进行的量。

还应注意，在上述示例中，在两两考虑的颗粒之间考虑相互作用力以便计算相互作用势并更新那些势。不过，本发明还允许在计算涉及k个颗粒之间相互作用力的势时减小相应的计算负荷，k严格大于 2。在此情况下，根据前一模拟步骤确定的相互作用势通过有利地利用当k个颗粒是固定颗粒时这k个颗粒之间的相互作用力在当前模拟步骤和前一步骤之间不变化(因此不必重新计算)的事实来确定当前的相互作用势。然后从施加于颗粒的总力中减去与包括在前一模拟步骤和当前步骤之间已移动的颗粒的颗粒的k元组(k_uplets)有关的前一步骤中计算出的力。因此根据与已移动的颗粒的新位置来计算并向如此获得的作用于颗粒的总力添加包括与已移动的颗粒的颗粒的k元组有关的当前力。

对于k＝2，或者k不同于2，可以执行类似的操作以便更新系统的势能，该势能被认为是不超过k个颗粒之间的势能之和。还可以执行类似的操作以便更新依赖于不大于k个颗粒的位置的数据结构或值，其中k是大于或等于1的任何整数。

例如，在涉及5个或更多个颗粒的集合的模拟中，要计算的信息包括所考虑的5个颗粒的重心，所述重心随时间变化但仅要在每10 个模拟时间步长中确定。如果在重心上次确定的时刻和当前时刻之间与前5个颗粒对应的自适应逆质量矩阵的项被设置为零，则颗粒并未移动并且不必更新重心。

在另一示例中，如果在重心上次确定的时刻和当前时刻之间与所考虑的前4个颗粒对应的自适应逆质量矩阵的项被设置为零，但是与第5个颗粒对应的项未被在某个时刻设置为零(因此第5个颗粒已经从上次计算重心时移动)，则重心如下逐步更新：

-把g乘以5；

-从g中减去第5个颗粒的前一位置；

-向g添加第5颗粒的新位置；

-把g除以5。

因此利用了与前4个颗粒对应的自适应逆质量矩阵的项为零的事实。

在另一实施例中，本发明提出允许计算和模拟要被加速的一组对象的方法和设备。使用自适应哈密顿算符允许在模拟期间激活或禁止满足某些标准的对象在位置方面的自由度。因此可以减少更新与这些对象有关的力或势所需的计算量。

在实施对举例考虑的NVT统计集中的一组E(具有恒定数量的颗粒，容积和温度)的行为的模拟中，例如基于上述自适应哈密顿算符使用根据朗格文动力学的模拟(例如参见“Free energy computations：a mathematical perspective”，T.Lelievre等人，Imperial College Pr，2010)。

朗格文动力学方程是：

其中：

-t-＞dW_t是标准3维布朗运动函数，-σ和γ是满足以下涨落耗散关系(fluctuation-dissipation relation)的3N＊3N矩阵：σσ^T＝2γ/β，其中β＝1/k_BT，k_B是玻尔兹曼常数并且T是温度；

-是自适应哈密顿算符相对于变量p的梯度；

-是自适应哈密顿算符相对于变量q的梯度。

在NVT模拟的情况下，对于积分算法，时间步长的计算可以如下进行：用于方程的朗格文部分的半个时间步长，用于方程的哈密顿部分的一个时间步长，以及再一次用于方程的朗格文部分的半个时间步长。

然后获得以下公式(5)：

其中，{G_k}是独立同分布(具有零均值和等于单位矩阵的协方差)的随机高斯矢量的序列。

公式(5)的一个方程包括右手和左手项中的项p_n+1。为了求解该隐式方程，使用例如固定点算法。

通过用适于NVT情况的公式(5)替代适于NVE情况的公式(3) 的考虑进行的，与上述程序P相似的程序适于执行与步骤101、102、 103类似的步骤，用于基于根据发明的自适应哈密顿算符H_A来更新 p_n和q_n的值。

使用根据本发明的自适应模拟(即使用自适应哈密顿算符)来还在NVT中的模拟期间计算平均值，以便确定否则将会使用常规哈密顿算符进行模拟来计算的值。

事实上，

还写作：

或

H_A＝H+V_A(p，q)，

其中且H是标准哈密顿算符。

如果借助使用标准哈密顿算符H进行的模拟获得的参数p和q 获得的变量A的平均值被表示为<A>_H并且借助使用自适应哈密顿算符 H_A进行的模拟获得的参数p和q获得的变量A的平均值被表示为，则<A>_H的值可以从的值通过以下等式找到：

其中k_BT是玻尔兹曼常数k_B和温度T的乘积。

可以从该等式示出当变量A仅依赖于位置并且自适应哈密顿算符是可分的(即Φ仅依赖于力矩而不是位置)时，

因此，借助自适应哈密顿算符获得的平均值等于借助标准哈密顿算符获得的平均值，这是有利的。

应注意，以上考虑的欧拉或朗格文积分器来执行本发明。然后可以使用任何类型的积分算法。

在以上考虑的情况下，颗粒的运动在所考虑的位移空间的所有维度上都是固定的。在另一实施例中，颗粒的运动仅在位移的一个或一些轴线上是固定的，这对于某些类型的运动可以是有用的。

用于更新与该轴线或那些轴线平行的力的分量所需的计算数量因此减少了。

在上述实施例中，颗粒的位置当其动能低于阈值时是固定的。

在可以或可以不与前面实施例结合的另一实施例中，当一个颗粒的力矩p采用某些预定值(离散值，或者一个或不同的范围值)，或者当包括力矩p和位置q的数偶采用某些固定值时，在至少一个模拟时间步长期间该颗粒是固定的。

在一个实施例中，颗粒组的位置是固定的。例如，仅当所有动能值e₁到e₁₀低于某个阈值的情况下，把动能e_i(i＝1到10)的零值分配给与颗粒a_i有关的自适应逆质量矩阵的对角项。该设置可以在某些势的计算中使计算加速或者允许更精确的计算。

举例来说，从把颗粒高速发送到初始不动的系统触发的震动开始使用兰纳-琼斯势(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在自适应模拟的上一步骤的坐标矢量，是该同一颗粒在参考模拟的上一步骤的坐标矢量。

例如，对于其中和(摩尔千卡，i＝1到N) 的自适应模拟，获得了用于执行该自适应模拟所需的计算时间相对于与参考模拟有关的计算时间的值等于2.5的加速因子，RMSD＝ 0.0114S并且Δq_max＝0.18S，其中S是所用的兰纳-琼斯势中的平衡距离。

对于其中和(i＝1到N)的自适应模拟，获得了用于执行该自适应模拟所需的计算时间相对于与参考模拟有关的计算时间的值等于5的加速因子，RMSD＝0.0612S并且Δq_max＝0.3S。

对于其中和(i＝1到N)的自适应模拟，获得了用于执行该自适应模拟所需的计算时间相对于与参考模拟有关的计算时间的值等于10的加速因子，RMSD＝0.359S并且Δq_max＝13.94S。

通过执行与NVT颗粒集合有关的本发明还发现可观的增益。因此，根据本发明的方法允许在对行为较小改变的情况下大大加速计算。

Claims (12)

1.一种对元件的系统进行模拟的方法，根据该方法所述元件的行为是基于元件的系统的哈密顿算符H确定的以使得其中p是表示元件的力矩的矢量，q是表示元件的位置的矢量，M^-1是随元件的质量而变化的对角矩阵，以及V是系统的势能，所述方法是由计算机执行的并且其特征在于所述方法包括以下步骤：当力矩矢量p采取与至少一个元件相关的某些预定值时，零值被分配给矩阵M^-1的与该至少一个元件相关的至少一个对角项。

2.按照权利要求1所述的对元件的系统进行模拟的方法，其中所述方法包括以下步骤：对于所述元件中的至少一个元件，如果表示所述元件的动能的参数具有低于第一严格为正的阈值，则把零值分配给矩阵M^-1的与该至少一个元件相关的至少一个对角项。

3.按照权利要求1或2所述的对元件的系统进行模拟的方法，其中，当所述至少一个元件的动能大于第二严格为正的阈值时随元件的质量而变化的矩阵M^-1的对角项被分配最大值。

4.按照权利要求1或2所述的对元件的系统进行模拟的方法，其中，当包括所述至少一个元件的力矩和所述至少一个元件的位置的数偶采用某些预定值时把零值分配给矩阵M^-1的与该至少一个元件相关的至少一个对角项。

5.按照权利要求1或2所述的对元件的系统进行模拟的方法，包括基于所述哈密顿算符在相继的模拟时刻确定至少一个信息的值的步骤，所述确定至少一个信息的值的步骤利用与在前一模拟时刻为矩阵M^-1的对角项分配零值所针对的元件k元组相关的信息的值因此在至少所述前一模拟时刻和当前模拟时刻之间未改变的事实，k是大于等于2的整数，并且在零值未被分配给矩阵的与元件k元组中的每个元件相关的对角项时通过执行以下步骤来计算在当前模拟时刻与构成k元组一部分的给定元件相关的信息的值：

-通过从在前一模拟时刻确定并且与所述给定元件相关的信息的值中减去在前一模拟时刻至少与所述元件k元组相关联并且与所述给定元件相关的所述信息的值来计算与所述给定元件相关的所述信息的工作值；以及

-向所述工作值添加在当前模拟时刻确定的至少与所述元件k元组相关联并且与所述给定元件相关的信息的值。

6.按照权利要求1或2所述的对元件的系统进行模拟的方法，其中，在当前计算时刻，由低于给定阈值的距离分隔的元件对的当前列表被在当前模拟时刻准备并且与在前一模拟时刻准备的由低于给定阈值的距离分隔的元件对的前一列表相比较，

并且在当前模拟时刻与给定元件相关的信息的值被基于包括所述给定元件的各元件对通过执行以下步骤来计算：

-如果所考虑的对仅存在于前一列表中或者如果把所述给定元件连接到对中的另一元件的矢量已经在前一模拟时刻和当前模拟时刻之间变化，则通过从在前一模拟时刻确定并且与所述给定元件相关的信息的值中减去与所述给定元件相关的信息的与所考虑的对中的另一元件相关联的值，来计算工作值；

-如果所考虑的对仅存在于当前列表中或者如果把所述给定元件连接到对中的另一元件的矢量已经在前一模拟时刻和当前模拟时刻之间变化，则通过向所述工作值添加与所考虑的对中的另一元件相关联并且与所述给定元件相关的信息的值来确定在当前模拟时刻与所述给定元件相关的信息的值。

7.按照权利要求1或2所述的对元件的系统进行模拟的方法，其中，在当前计算时刻，满足某些条件的元件k元组的当前列表被在当前模拟时刻准备并且与在前一模拟时刻满足所述条件的元件k元组的前一列表相比较，k是大于等于2的整数，

并且在当前模拟时刻与一个元件相关的信息的值被基于包括所述元件的k元组通过执行以下步骤来计算：

-当所述k元组仅存在于前一列表中时或者当与所述k元组相关联并且与所述元件相关的信息的值已经在前一模拟时刻和当前模拟时刻之间变化时，则通过从在前一模拟时刻确定并且与所述元件相关的信息的值中减去在前一模拟时刻与所述k元组相关联并且与所述元件相关的信息的值，来计算临时值；

-当所述k元组仅存在于当前列表中时或者当与所述k元组相关联的信息已经在前一模拟时刻和当前模拟时刻之间变化时，则通过向所述临时值添加在当前模拟时刻与所述k元组相关联并且与所述元件相关的信息的值来确定在当前模拟时刻与所述元件相关的信息的值。

8.按照权利要求1或2所述的对元件的系统进行模拟的方法，其中，元件的位置空间被划分为单元并且在前一模拟时刻和当前模拟时刻中的每一个模拟时刻，每个元件根据在所述模拟时刻确定的位置坐标与归属单元相关联，以及对于第一元件执行以下步骤，所述第一元件为使得在当前模拟时刻矩阵M^-1的与所述第一元件相关的项还未被分配零值的元件：

-确定在前一模拟时刻第一元件的归属单元；

-对于每个第一元件，在所述归属单元中或所述归属单元的给定邻域中的单元中，确定在前一模拟时刻与所述第一元件相距小于给定阈值的第二元件；通过从在前一模拟时刻确定并且与所述第一元件相关的信息的值中减去与所述第二元件相关联并且与所述第一元件相关的所述信息的值来计算工作值；

-确定在当前模拟时刻的第一元件的新归属单元；

-对于每个第一元件，在所述新归属单元中或所述新归属单元的给定邻域中的单元中，确定在当前模拟时刻与所述第一元件相距小于给定阈值的第三元件；

-通过向工作值添加与所述第三元件相关联并且与第一元件相关的信息的值来确定在当前模拟时刻与所述第一元件相关的信息的值。

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

10.按照权利要求1或2所述的对元件的系统进行模拟的方法，包括：在某些模拟时刻确定信息I的步骤。

11.按照权利要求1或2所述的对元件的系统进行模拟的方法，包括在某些模拟时刻确定信息I的步骤，所述确定信息I的步骤利用如下事实：该信息I不改变并且当信息I在前一模拟时刻已确定时并且在至少所述前一模拟时刻和当前模拟时刻之间零值被分配给矩阵M^-1的对角项的对应集合时不必再次确定。

12.按照权利要求1或2所述的对元件的系统进行模拟的方法，包括在某些模拟时刻确定势能或相互作用力的步骤。