CN104756115A - 用于模拟一组元件的方法 - Google Patents
用于模拟一组元件的方法 Download PDFInfo
- Publication number
- CN104756115A CN104756115A CN201380056652.0A CN201380056652A CN104756115A CN 104756115 A CN104756115 A CN 104756115A CN 201380056652 A CN201380056652 A CN 201380056652A CN 104756115 A CN104756115 A CN 104756115A
- Authority
- CN
- China
- Prior art keywords
- node
- particle
- tree
- matrix
- partiald
- Prior art date
- Legal status (The legal status 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 status listed.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
- G16C10/00—Computational theoretical chemistry, i.e. ICT specially adapted for theoretical aspects of quantum chemistry, molecular mechanics, molecular dynamics or the like
-
- 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
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Computing Systems (AREA)
- Health & Medical Sciences (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- General Health & Medical Sciences (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明涉及一种基于汉密尔顿算符(等式I)对由树表示的元件的系统进行模拟的方法,树包括叶节点和包括根结点R的内部结点,每个叶节点表示一个元件,p是力矩的矢量,q是元件的位置的矢量,以及V是系统的势能:当预定的条件被验证时,通过把矩阵M-1定义为等于ΦR,对树的给定结点的从属的至少一些元件施加同样的平移运动,对于包括k个子结点的树的任何结点A,有递归等式(等式II)。
Description
技术领域
本发明涉及一种用于模拟一组元件的方法,根据该方法基于与元件系统(该组元件的动能和势能之和)相关联的汉密尔顿算符 来在相继的模拟步长确定元件的行为,其中,p是表示元件力矩的矢量,V是系统的势能,以及M-1是这些元件的质量的函数的对角矩阵(在此情况下,该矩阵可以是元件位置的函数)。
背景技术
在此情况下势能V例如仅是元件位置的函数。在其它情况下,势能V还可以依赖于元件的力矩。作用于元件上的力可以从该势能推导。
一组元件的模拟允许研究这样的一组元件的行为并分析其属性:在元件的相继位置和力矩方面的位移,元件之间的位移的相关性,结构的变化,元件之间的相互作用的增加和减小,平均采用的配置,相关联能量的演变等等。这些元件可以表示机械体,例如天体或流体;颗粒,例如原子或分子,例如蛋白质,流体等。
模拟一组元件的常规方式是考虑该组元件的汉密尔顿算符并从中推导力矩方程,并根据这些方程推导元件的力矩。
WO2009/007550例如描述了模拟一组元件的技术。
该组元件的演变有时必须被长时段上模拟以便能够观察某些现象或能够计算某些统计量。这种模拟的计算时间、以及计算成本有时变得非常大。已经提出的很多方法来加速对一组元件的模拟并收集统计量。
发明内容
本发明旨在提出用于克服这些问题的新方案。
为此,根据第一方面,本发明提出一种对指定类型的元件的集合进行模拟的方法,其特征在于所述方法包括以下步骤:
-把元件的系统表示为包括叶节点和内部结点的k叉树,所述叶节点每个都表示相应的元件,内部结点中包括根结点,
-当所确定的条件被验证时,对于当前模拟的步长,通过把矩阵M- 1定义为等于ΦR,对树的给定结点的从属的至少一些元件施加同样的平移运动,其中R是树的根结点并且对于包括k个子结点A1、A2,…,Ak的树的任何结点A,有以下递归等式,根据该等式:
矩阵E是由各自尺寸为d*d的nA*nA个块形成的尺寸为dnA*dnA的矩阵,块等于尺寸为d的单位矩阵,nA等于结点A的从属的元件的数量,d是其中颗粒演变的空间的尺度,mA是从属于结点A的这nA个元件的质量之和,ρA是介于0和1之间的结点A的约束函数,该函数的值当A是叶节点时或者当同样的运动被施加到从属于给定结点的元件时等于1;在叶节点Ai的情况下,矩阵等于由结点Ai表示的颗粒乘以尺度为d的单位矩阵所得的逆矩阵。
本发明允许执行以下模拟:该模拟需要较小计算体积并因此较少计算时间来确定根据这些模拟的元件的性能,例如元件的势能、应用于元件的力、位置和/或力矩。
在实施方式中,根据本发明的对元件的集合进行模拟的方法还包括一个或更多个以下特征:
-对于当前模拟步长,对树的给定结点的从属的至少一些元件根据所述元件的力矩函数的值来施加同样的平移运动;
-对于当前模拟步长,根据由所取的值对树的给定结点的从属的至少一些元件施加同样的平移运动,其中是结点Ai的元件的质量之和,是子结点Ai的元件的力矩的矢量和,||·||是矢量的范数,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个颗粒ai的一个组E的模拟:i=1至N。
与该组E相关联的汉密尔顿算符H(p,q)常写作如下(参见例如“理解分子模拟:从算法到应用”,Frenkel D.,Smit B):
p是表示所有颗粒的力矩的矢量,q是表示所有颗粒位置的矢量,M-1是作为颗粒的质量的函数的对角矩阵。
V(q)是N个颗粒之间的相互作用势;它是颗粒的位置的函数并将被考虑为独立与力矩。
在3维空间中,例如,在坐标系(X,Y,Z)中,每个颗粒ai的力矩被写为(pi,x,pi,y,pi,z)并且每个颗粒ai的位置被写为(qi,x,qi,y,qi,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]=mi,mi是颗粒ai的质量。
这就是汉密尔顿算符的通常定义,以下将被称为标准汉密尔顿算符。
根据本发明,称作自适应相对汉密尔顿HAR的汉密尔顿算符是以下:
其中,Φ(p,q)是称作自适应相对逆质量矩阵的3N×3N块对角矩阵,替代M-1并依赖于矢量p并可选地依赖于矢量q。
更准确地,Φ(p,q)的值是在每次模拟中确定的,函数Φ被选择为能够把颗粒的至少一个子集限制为在模拟的一个或某些步长全部跟随同一运动。构成子集的颗粒被通过条件验证来确定,例如旨在验证在当前模拟步长之前,颗粒被接近的运动所激励。
这些设置导致出现颗粒的一个或多个子集,在所考虑的模拟步长期间子集中的颗粒被同一平动位移激励(即,在其中颗粒具有相同位移的时间范围期间连接这些颗粒中的任两个的矢量保持恒定)。
针对所考虑的模拟步长,这样子集的颗粒可以因此被相似对待,如构成唯一刚性体。
仅取决于颗粒相对位置的信息,例如颗粒之间的相互作用力,因此在同一个刚性子集内在所考虑的模拟步长结束处不需要被更新。
当子集的颗粒的运动不再受自适应逆质量矩阵限制时,每个颗粒采取自己专有的运动。
在本发明的实施方式中,在两类行为(限制运动/自由运动)之间施加连续过渡。
根据该自适应汉密尔顿算符HAR,推导出限定是的自适应运动方程,是是矢量p和q相对于时间t的导数。
例如,在此处举例考虑的NVE(具有恒定颗粒数、体积和能量的组E)组中实施模拟时,汉密尔顿算符的值(根据本发明或标准是自适应的)随时间是恒定的,并且自适应运动方程是:
。
下面的示例与NVE组相关,但是本发明当然可以应用于其它组,例如NVT组(具有恒定颗粒数、体积和温度的组E)。
可以使用多个不同方案来定义构成刚性子集的颗粒,而无需测试集合的颗粒之间的组合的可能性。
第一方案是对刚性体的预定对施加约束,这些刚性体可被认为是要合并成同一刚性体的候选。例如,注意A1、…An,即一组编索引的刚性体,可以把索引组织成二叉树,以使得当具体合并条件满足时A1仅与A2合并成刚性体(A1+A2),A3仅与A4合并等,并且在二叉树的上层级,(A1+A2)仅与(A3+A4)合并成等。这样的树在最高层级包括称作R结点的根结点,以及在最低层级包括与各个考虑的颗粒对应的叶结点。除了叶结点外,树的所有结点称为内部结点。
另一方案是对刚性体的可能包括多于两个刚性体的预定子集施加约束。因此可以把索引组织为n叉树以使得在同一时间A1仅与A2、A3…、An合并等。
另一方案是考虑从邻域的研究得出的图,图的结点是颗粒或刚性体,并且两个结点之间的边表示由边连接的两个结点之间的距离小于距离阈值。
下面详述的实施方式是基于如关于第一方案描述的二叉树结构;二叉树结构的优点是准许针对层次结构中的内部结点的判决度量的增量且递归的更新,以及例如针对部分能量和力的增量更新。
在这样的树中,每个叶结点表示集合中的单个颗粒,每个中间结点表示由所述中间结点的叶结点表示的颗粒组成的颗粒的子集。称体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是由nc×nc个等于3维单位矩阵的3×3块形成的尺寸为3nc×3nc的矩阵,nc等于结点C的叶节点的数量(即,结点C的直接子结点中或这些子结点的附属结点中的叶节点的集合);
-mc是由结点C的这些叶结点nc表示的颗粒的质量之和;
-ρ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η5+15η4-10η3+1,其中δ=εf-εr表示在体的刚性行为和自由行为之间的过渡区域的宽度。
在本发明的实施方式中,εC被独立于两个体A和B的力矩的关系进行选择,其中:
。
其中:A和B是结点C的子结点;
矢量是作为与结点C对应的体中的叶结点的所有颗粒的力矩之和:并且是矢量和之和:
该发明可以被针对k叉树进行概括。在该情况下,对于具有k个子结点Ai的母结点,1≤i≤k,惯性逆矩阵是:
并且
在该公式和公式6中的系数1/2可以被别的系数代替,例如1,这在阈值εf和εr被同样修改(即,乘以等于其它常数的2倍的系数)的情况下不改变模拟。
本发明还可以不同地概括:可以针对不同的内部结点选择不同的限制函数ρC。例如,直到树的某个级别,ρC被根据公式4和5定义,因此,结点可以集中为整体或分离,并且在大于该级别的树的级别中,ρC总等于零。在具体实施方式中,可以定义持续激活的系统E的子集(该子集的颗粒既不彼此合并,也不与系统的其它部分合并,因此这些颗粒在整个模拟中总是遵循不受限制的自由的运动,本发明的这样的实施方式可以例如应用于溶剂中的聚合物):该子集中的所有颗粒应该放置在单独的子树中并且针对该子树的所有内部结点值ρC应该固定为0并不再更新。在所有情况下,函数ρC应该在模拟开始之前针对所有内部结点被定义并且不应在模拟期间改变以便具有稳定的模拟。
根据本发明,选择在相对自适应汉密尔顿算符HAR的公式(1)中使用的自适应相对逆质量矩阵Φ,等于ΦR,其为相对于表示颗粒的集合E的树的根结点R而言根据递归公式3定义的矩阵:Φ(p,q)=Φ(p)=ΦR。
因此具有汉密尔顿算符HAR,如其是可分的,这是因为ΦR仅取决于颗粒的力矩(而不是位置)。
由公式2定义的运动方程因此表示如下:
还从公式3清楚可见如果结点C是自由的(ρC=0),则矩阵ΦC由对角线上的对应于子结点A和子结点B的两个块组成。
在根据本发明的模拟的第一示例中,考虑具有包括通过刚度为1的弹簧彼此连接的质量为1的两个颗粒P1、P2的尺寸的系统。对应于系统的二叉树是根结点C,其两个子结点A和B对应于这两个颗粒。在给定初始力矩的情况下,两个颗粒在空间中随时间振荡。
根据现有技术基于汉密尔顿算符H的模拟针对这些颗粒计算出的在空间中的轨迹d在图1中被根据模拟时间表示。
根据本发明基于自适应相对汉密尔顿算符HAR的模拟针对这些颗粒计算出的轨迹d在图2中被针对δ和εf的不同值根据模拟时间表示。因此如在这些轨迹上表示那样,在某些模拟时刻(有时候从头开始),两个颗粒的力矩根据公式6变得接近。这随着弹簧几乎被压缩并释放而发生。根据本发明,根结点C因此变得刚性。结果,在一般情况下,刚性结点C的子结点遵循同一平移运动,在当前情况下,由于力矩的保持,所以颗粒停止。2个颗粒的轨迹因此变得平行。颗粒的力矩继续演变,然而结果是在给定的模拟步长,结点C的刚性条件不再满足并且两个颗粒采取自己的运动。这因此在模拟过程中周期性地产生。
如轨迹所示,针对δ值增大的εf的同一值,轨迹的过渡区域更大。针对εf值增大的δ的同一值,刚性区域更长并更平。
因此根据本发明,被定义为等于ΦR的矩阵Φ指定如何并何时在模拟过程中激活或禁止颗粒子集的相关位置自由度。
在本发明的实施方式中,图3中表示的信息设备1被用于根据上述的本发明原理来实施N个颗粒的集合E的模拟。
该设备1包括计算机,计算机尤其包括:存储器2,适于存储下面相继描述的软件程序和计算出的参数值(全局相互作用力,局部相互作用力,相互作用势,位置,力矩等);微处理器3,适于执行软件程序、尤其是下述程序P的指令;以及人机界面4,例如包括键盘和屏幕,分别用于输入用户的指令和向用户显示信息,例如图2所示的曲线。
在所考虑的本发明的实施方式中,存储器2包括基于与公式7对应的方程随时间的积分来模拟颗粒ai的集合E的行为的程序P,i=1到N。
因此根据本发明的自适应相对汉密尔顿算符HAR是可分的,所以可以使用若干个积分技术。
参照图4,程序P包括软件指令,软件指令当被微处理器3执行时适于实施预备步骤100并在程序P的与计算时刻hn+1=h0+(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-angle quasi-statics:ageneral simulation method with applications to protein structureanalysis 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是通过使用常规方法实现的。
因此,对于每个颗粒ai,通过为pi,n+1填写颗粒ai的力矩pi的值,当当前迭代n+1时,该值通过以下等式获得:
pi,n+1=pi,n+fi,n+1h,
h是模拟的时间步长,fi,n+1是作用于颗粒ai上的全局相互作用力,i=1到N,这是由于系统E中的所有其它颗粒施加的相互作用(出现在如当前迭代中刚确定的根结点的局部力的表中)引起的。
在每个时间步长,基于等式的颗粒位置的“自然”更新的复杂度为O(N3),这是由于项该项是矩阵相对于矢量的倒数,因此是三维对象,对于1到N:
在所考虑的实施例中,对于具有两个子结点A和B的每个内部结点C,通过从数据结构QC、RC和SC的叶结点进行递归更新,复杂度可以减小为O(N.logN)。
注意:
其中pA是作为结点A的叶节点的所有颗粒的力矩矢量,pB是作为结点B的叶节点的所有颗粒的力矩矢量,项QC=ΦC(pC)pC对应于给出的等式的第一项:
其中{x}是维度等于nC的矢量,其nC个组成部分的每一个等于元素x。
类似地,对于结点C考虑项
给出的等式的最后一项可以通过使用等式9来递归更新。
事实上:
函数(·)表示标量函数。
因此,步骤103包括针对包括两个子结点A和B的每个结点C的结构QC、RC和SC的更新。
应注意步骤100还包括这些结构的初始化步骤,其中对于每个内部结点C:ρC=0并且QC、RC和SC还是0,并且对于表示颗粒的每个子结点F:ρF=1;QF=pF/mF;SF=0,RF=0。
在步骤103中,度量QC、RC和SC以及εC和ρC的更新相对于结点C的更新被从根结点沿着结点的架构递归地执行:
-如果结点C是内部的,具有两个子结点A和B:
-更新结点A的度量;
-更新结点B的度量;
-借助等式6来更新εC;
-更新ρC:ρC=μ(εC)ρAρB;
-借助等式9来更新QC;
-借助等式10来更新RC;
-借助等式11来更新SC;
-如果结点C是表示颗粒aF的叶结点:QF=pF/mF。
为了优化计算,矢量QC、RC和SC的更新可以仅针对ρC>0的结点来进行,这是因为这些度量可以不用于自由结点,所以不需要更新它们。
在一般情况下,这些度量可以以不遍历树的方式来更新,只要在刚性结点在顶部的情况下在子树中更新即可。
一旦确定了这些度量,颗粒的位置就在步骤104中更新。
步骤104的实施可以在结点C的位置的更新之后进行,每个内部结点C被考虑作为包括两个子结点A和B:
-如果ρC=0(结点C是自由结点):
-更新结点A的位置;
-更新结点B的位置;
-否则,对于结点C的下属每个叶结点k,无论直接还是间接:
qi,n+1=qi,n+(QC[k]+0.5SC[k])h,
qi,n+1是颗粒ai的位置的值,针对被考虑的第n+1次迭代,并且h是模拟时间,QC[k]是矢量QC的第k个元素,以及SC[k]是矢量SC的第k个元素。
颗粒的位置因此就确定了。
在该实施方式中,对于所有颗粒以及树的所有内部结点,在颗粒i的标识符及其在结点C的子结点中的序号之间的一一对应关系必需建立,例如在模拟初始化期间建立。
在位置的该实施方式中,ρC>0的结点的度量不被使用并且因此可以不计算。如果在步骤103中更新所有度量,则还可以如下确定所有颗粒ai(1≤ai≤N)的位置:qi,n+1=qi,n+(QR[k]+0.5SR[k])h,其中索引R对应于根结点。
对于树的每个内部结点C,存储与结点C的直接或间接从属的叶结点的数量相等的4个长度矢量:QC、RC和SC以及局部力的表。位置更新的空间复杂度因此是O(N.logN)。时间复杂度也是O(N.logN),这是因为这些结构应在每个时间步长更新,并且QC总是对所有结点(其中包括叶结点)更新。
在一般情况下,位置可以不以遍历树的方式更新。
然后,如果未达到模拟的最大持续时间,则实施程序P的新的迭代。
作为示例,从把颗粒高速发送到初始不动的系统触发的震动开始使用兰纳-琼斯势(Lennard-Jones potential)(Em/kB=120开尔文,其中Em是最小能量,平衡距离S=3.4埃,截止距离8埃,势在7.5和8之间被平滑截短)来执行N=5930个颗粒ai(i=1到N)的NVE集合的2个演化维度的4次模拟,每个颗粒具有1g/mol的质量:基于标准汉密尔顿算符的一个参考模拟和三个自适应模拟,也就是说使用根据如上所述的本发明的方法的自适应相对汉密尔顿算符(时间步长尺寸0.0488飞秒(fs),7000时间步长,总模拟时间342fs)。
对于使用自适应相对汉密尔顿算符的每次模拟,相对于标准模拟的波动的平方根,称为RMSD,被给出作为最大颗粒位移误差Δqmax。
其中,qi是颗粒ai在自适应模拟的上一步的坐标矢量,是该同一颗粒在参考模拟的上一步的坐标矢量。
例如,对于其中εr=0,01kcal/mol和(δ=1kcal/mol)的自适应模拟,获得了用于执行该自适应模拟所需的计算时间相对于与参考模拟有关的计算时间的值等于1.5的加速因子,RMSD=8.6S并且Δqmax=89.4S,其中S是所用的兰纳-琼斯势中的平衡距离。
对于其中和(δ=1kcal/mol)的自适应模拟,获得了用于执行该自适应模拟所需的计算时间相对于与参考模拟有关的计算时间的值等于1.55的加速因子,RMSD=8.7S并且Δqmax=89.4S。
对于其中和(δ=1kcal/mol)的自适应模拟,获得了用于执行该自适应模拟所需的计算时间相对于与参考模拟有关的计算时间的值等于1.6的加速因子,RMSD=12.3S并且Δqmax=89.4S。
因此,根据本发明的方法允许使计算加速,而潜在较小的改变性能。
结果取决于针对颗粒的系统所选择的树状表示。对于上述的自适应模拟,通过把系统一分为二来从上至下构建了二叉树。
在本发明的所考虑的实施方式中,已经考虑了把子集两两融合为刚性集合,在其它实施方式中,融合成一个刚性集合的子集的数量被取为大于2。
在刚性和自由的两个状态之间的过渡区域可以被取为更大或更小的宽度。
Claims (8)
1.一种对元件的系统进行模拟的方法,根据该方法所述元件的行为是基于元件的系统的汉密尔顿算符H在相继的模拟步长上确定的以使得其中p是表示元件的力矩的矢量,q是表示元件的位置的矢量,M-1是随元件的质量而变化的矩阵,以及V是系统的势能,
其特征在于所述方法包括以下步骤:
-把元件的系统表示为包括叶节点和内部结点的k叉树,所述叶节点每个都表示相应的元件,内部结点中包括根结点,
-当所确定的条件被验证时,对于当前模拟的步长,通过把矩阵M- 1定义为等于ΦR,对树的给定结点的从属的至少一些元件施加同样的平移运动,其中R是树的根结点并且对于包括k个子结点A1、A2,…,Ak的树的任何结点A,有以下递归等式,根据该等式:
矩阵E是由各自尺寸为d*d的nA*nA个块形成的尺寸为dnA*dnA的矩阵,块等于尺寸为d的单位矩阵,nA等于结点A的从属的元件的数量,d是其中颗粒演变的空间的尺度,mA是从属于结点A的这nA个元件的质量之和,ρA是介于0和1之间的结点A的约束函数,该函数的值当A是叶节点时或者当同样的运动被施加到从属于给定结点的元件时等于1;在叶节点Ai的情况下,矩阵等于由结点Ai表示的颗粒乘以尺度为d的单位矩阵所得的逆矩阵。
2.按照权利要求1所述的对元件的系统进行模拟的方法,其中对于当前模拟步长,对树的给定结点的从属的至少一些元件根据所述元件的力矩函数的值来施加同样的平移运动。
3.按照权利要求1或2所述的对元件的系统进行模拟的方法,其中,对于当前模拟步长,根据由所取的值对树的给定结点的从属的至少一些元件施加同样的平移运动,其中是结点Ai的元件的质量之和,是子结点Ai的元件的力矩的矢量和,||·||是矢量的范数,C是正常数。
4.按照权利要求3所述的对元件的系统进行模拟的方法,其中,对于针对内部结点A的当前模拟步长:
-如果εA小于第一阈值,则把ρA定为1以便对从属于结点A的元件施加同样的平移运动;
-如果εA大于比第一阈值大的第二阈值,把ρA定为0。
5.按照权利要求4所述的对元件的系统进行模拟的方法,其中,对于当前模拟步长,如果εA介于第一阈值和第二阈值之间,则把ρA定为等于由变量εA的5阶插值函数所取的值。
6.按照前述权利要求中任一项所述的对元件的系统进行模拟的方法,包括在基于所述汉密尔顿算符的相继模拟的时刻确定至少一个信息的值的步骤,所述步骤利用了以下事实:对于当前模拟的时刻与已经施加了同样的平移运动的元件相关的信息的值取决于与所述元件相关的位置并且因此不变化。
7.按照前述权利要求中任一项所述的对元件的系统进行模拟的方法,其中,与所述元件相关的信息包括所述元件的势能和/或应用到所述元件的相互作用力。
8.一种对元件的系统进行模拟的计算机程序(P),包括当由计算装置执行该程序时用于实施根据权利要求1至7之一所述的方法的步骤的软件指令。
Applications Claiming Priority (3)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
FR1258330A FR2995109A1 (fr) | 2012-09-06 | 2012-09-06 | Procede de simulation d'un ensemble d'elements, programme d'ordinateur associe |
FR1258330 | 2012-09-06 | ||
PCT/EP2013/067630 WO2014037236A1 (fr) | 2012-09-06 | 2013-08-26 | Procédé de simulation d'un ensemble d'éléments |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104756115A true CN104756115A (zh) | 2015-07-01 |
CN104756115B CN104756115B (zh) | 2018-09-04 |
Family
ID=47624208
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201380056652.0A Active CN104756115B (zh) | 2012-09-06 | 2013-08-26 | 用于模拟一组元件的方法 |
Country Status (5)
Country | Link |
---|---|
US (1) | US20150254378A1 (zh) |
EP (1) | EP2893473A1 (zh) |
CN (1) | CN104756115B (zh) |
FR (1) | FR2995109A1 (zh) |
WO (1) | WO2014037236A1 (zh) |
Families Citing this family (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP6844930B2 (ja) * | 2017-09-05 | 2021-03-17 | 住友重機械工業株式会社 | シミュレーション方法及びシミュレーション装置 |
CN112179236B (zh) * | 2020-09-24 | 2022-06-24 | 青岛科技大学 | 一种面向装配的基于最小势能的平面装配性能评价方法 |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2004102159A3 (en) * | 2003-05-13 | 2005-01-27 | Penn State Res Found | Quantum mechanics based method for scoring protein-ligand interactions |
TW200601154A (en) * | 2004-04-20 | 2006-01-01 | Univ Mcgill | A nano molecular modeling method |
CN101019122A (zh) * | 2004-07-12 | 2007-08-15 | 阿托米斯蒂克斯公司 | 在非平衡条件下用于分子的量子化学模拟的方法和计算机系统 |
US20070233440A1 (en) * | 2006-03-29 | 2007-10-04 | International Business Machines Corporation | Reduced message count for interaction decomposition of N-body simulations |
US20090076780A1 (en) * | 2007-08-10 | 2009-03-19 | Fujitsu Limited | Method, apparatus and computer program for molecular simulation |
US20100030534A1 (en) * | 2008-07-28 | 2010-02-04 | Fujitsu Limited | Method, apparatus and computer program for simulating behaviour of thermodynamic systems |
Family Cites Families (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
KR101210339B1 (ko) * | 2006-10-10 | 2012-12-18 | 삼성전자주식회사 | 트리 구조에서의 노드 식별자 생성 방법 |
US20080275685A1 (en) * | 2007-05-01 | 2008-11-06 | Thomas Michael Gooding | Miss-accumulation in a binary space partitioning tree |
FR2917866B1 (fr) | 2007-06-20 | 2009-09-04 | Inst Nat Rech Inf Automat | Dispositif informatique pour la simulation d'un ensemble d'objets en interaction et procede correspondant |
-
2012
- 2012-09-06 FR FR1258330A patent/FR2995109A1/fr active Pending
-
2013
- 2013-08-26 EP EP13756385.4A patent/EP2893473A1/fr not_active Withdrawn
- 2013-08-26 US US14/426,413 patent/US20150254378A1/en not_active Abandoned
- 2013-08-26 CN CN201380056652.0A patent/CN104756115B/zh active Active
- 2013-08-26 WO PCT/EP2013/067630 patent/WO2014037236A1/fr active Application Filing
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2004102159A3 (en) * | 2003-05-13 | 2005-01-27 | Penn State Res Found | Quantum mechanics based method for scoring protein-ligand interactions |
TW200601154A (en) * | 2004-04-20 | 2006-01-01 | Univ Mcgill | A nano molecular modeling method |
CN101019122A (zh) * | 2004-07-12 | 2007-08-15 | 阿托米斯蒂克斯公司 | 在非平衡条件下用于分子的量子化学模拟的方法和计算机系统 |
US20070233440A1 (en) * | 2006-03-29 | 2007-10-04 | International Business Machines Corporation | Reduced message count for interaction decomposition of N-body simulations |
US20090076780A1 (en) * | 2007-08-10 | 2009-03-19 | Fujitsu Limited | Method, apparatus and computer program for molecular simulation |
US20100030534A1 (en) * | 2008-07-28 | 2010-02-04 | Fujitsu Limited | Method, apparatus and computer program for simulating behaviour of thermodynamic systems |
Non-Patent Citations (2)
Title |
---|
ALLEN M P: "《Introduction to molecular dynamics simulation》", 《COMPUTATIONAL SOFT MATTER: FROM SYNTHETIC POLYMERS TO PROTEINS》 * |
SCHNEIDER R ET AL: "《Introduction to Molecular Dynamics》", 《COMPUTATIONAL MANY-PARTICLE PHYSICS》 * |
Also Published As
Publication number | Publication date |
---|---|
FR2995109A1 (fr) | 2014-03-07 |
CN104756115B (zh) | 2018-09-04 |
US20150254378A1 (en) | 2015-09-10 |
WO2014037236A1 (fr) | 2014-03-13 |
EP2893473A1 (fr) | 2015-07-15 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP7133894B2 (ja) | データに基づくインタラクティブ3dエクスペリエンス | |
Vanaerschot et al. | Stochastic multi-scale modelling of textile composites based on internal geometry variability | |
Tromme et al. | System-based approaches for structural optimization of flexible mechanisms | |
CN111195906B (zh) | 用于预测机器人的运动轨迹的方法和系统 | |
CN103049592A (zh) | 拟真式空间变化 | |
Bujny et al. | Topology optimization of 3D-printed joints under crash loads using evolutionary algorithms | |
Valero‐Lara | Reducing memory requirements for large size LBM simulations on GPUs | |
Nguyen et al. | Topology optimization using the p-version of the finite element method | |
JP6137175B2 (ja) | 階層型確率モデル生成システム、階層型確率モデル生成方法、およびプログラム | |
Yang et al. | Simulearn: Fast and accurate simulator to support morphing materials design and workflows | |
Zheng et al. | An improved local remeshing algorithm for moving boundary problems | |
Xie et al. | Adaptive topology optimization under suitably graded THB‐spline refinement and coarsening | |
JP5405641B2 (ja) | 挙動解析システム、挙動解析方法及び挙動解析プログラム | |
CN104756115A (zh) | 用于模拟一组元件的方法 | |
EP3374941A2 (en) | Dynamic design of complex system-of-systems for planning and adaptation to unplanned scenarios | |
Fox et al. | Learning Everywhere: Pervasive machine learning for effective High-Performance computation: Application background | |
Bonaventura et al. | A self adjusting multirate algorithm based on the TR-BDF2 method | |
CN104820729A (zh) | 钢轨全轧程仿真分析的温度继承方法 | |
Akçay Perdahcıoğlu et al. | Updating the Craig–Bampton reduction basis for efficient structural reanalysis | |
Li et al. | A hamiltonian formulation on manifolds for dynamic modeling of constrained mechanisms and energy-preserving integration | |
Zeng et al. | Efficient parallelization strategy for real-time FE simulations | |
CN104508667A (zh) | 用于模拟一组元件的方法及相关的计算机程序 | |
Alexandrov | Opportunities for breakthroughs in large-scale computational simulation and design | |
Li | Topology optimization of compliant mechanisms based on the BESO method | |
Grim-McNally et al. | Preconditioning parametrized linear systems |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
EXSB | Decision made by sipo to initiate substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |