分子动力学模拟中通用索引方法及其应用

分子动力学模拟中的通用索引方法及应1引分子动力学模拟中的通用索引方法及应1引不同分法不同分法2通用索直观想法与直观算法的直接对应，改变了传统的直观想法对应不直观算法的格局，有利于序实现。考虑完整性，本节首先介绍目前广直观想法与直观算法的直接对应，改变了传统的直观想法对应不直观算法的格局，有利于序实现。考虑完整性，本节首先介绍目前广泛采用的背景网格索引，然后介绍空间多级树据结构以及基于空间多级树的快速搜索索引方法方法如有限元计算中常用的连接关表方法不适合整体动态计算，在这里不再介绍背景网格索对于点物体或尺寸均匀的小物体果其空间分布不是特别适合用背景网行索引。这种方法广泛用于分子动力学、随机模拟方法（MC、光滑粒子动力学（SPH、质点(MPM)、离散元方法中用来搜索紧邻分子、粒子、物体。对于每个物体，可以给定一点作为这个物体的代表位置，一般选这个物体的中心或特征点。这样，每个物体都可用空的点来定位背景网格索引[6]的思路是：根据空间所有物体的位置设定一个箱子能够存放这些物将这个箱子按照长、宽及高分成规则的小箱子，将物体根据其空间位置放入相应的小箱中。这样就对空间物体建立了一种索引。在搜索物体时，根据搜索条件选择合适的几个小网格，搜索圆区域中的物体时，不用全部搜索，仅需计算与圆区域重叠的有哪些网格（粗框绿色矩形进入这些网格搜索这些网格内的物体即可图 物体的背景网格索引示意建立背景网格索引时，由于对进入网格的物体可以用数组保存，因此需要设计数组小，当物体在空间分布不均匀时会浪费内存。更方便的方法是采用链表,如图2所示。建背景网格索引的算法如下(以二维空间为例图 背景网格上的物体链表示小，当物体在空间分布不均匀时会浪费内存。更方便的方法是采用链表,如图2所示。建背景网格索引的算法如下(以二维空间为例图 背景网格上的物体链表示意⑴根据搜索要求设定网格尺寸⑵检查空间物体确定系统尺寸(x0,x1(y0,y1x0y1y01计算网格数组大小(I,J)，其中I1，⑶动态开设网格数组⑷对每个物体A，根据其空间位置(xA,yA)计算其所属的网格单元(i,j)，其xAx0yAy0ij将这个物体从链表尾部插入到此网格单元的链表中，如图3所示，其中红色圆点表物体所在网格单元格点g(i,j)，蓝色圆点代表物体Ag(i,j)obj代表格单元格点指向的物体。具体插入步骤分两步，如下Anextg(i,j)objg(i,j)objA图 将物体插入网格链表的示意图，（a）插入链表前，（b）插入链表搜索给定区域内物体的算法如下⑴计算与给定搜索区域重叠的图 将物体插入网格链表的示意图，（a）插入链表前，（b）插入链表搜索给定区域内物体的算法如下⑴计算与给定搜索区域重叠的网格单元，对这些网格单元循环⑵对每个网格单元内的物体循环，查找符合要求的物单元包含的物体数。背景网格索引的搜索计算量与关联网格的数量(Ngrid)和网格平均内含N体数(Nobj)有关，一般为(Ngrid)，其 1。由于网格单元越大，则关联网格元越少，网格平均内含物体数越多；网格单元越小，则关联网格单元越多，网格平均内含体数越少，同时随着网格单元变小，内存量会增加，因此需要选择合理的网格大小从而最程度减小计算量。一般选择网格平均内含物数目在3~5之间的背景网格。可以通过预先1/x)(y算物体的平均密度给出网格尺寸，例如，二维选择 N1/x)(yz N，其中为网格平均内含物体数背景网格索引的优点在于数据结构简单，在物体尺寸均匀和物体空间分布均匀时计（（2）由于物体的尺寸不进入索引，当物体尺寸不均匀时，难以支持；（3）维空间或体空间分布极不均匀时，需要的内存量和搜索计算量很大，代价很高空间多级树及通用快速搜空间多级树及通用快速搜2.2.1空间多级离散点离散点物体的空间多级树管理区域示意（a）二维离散点，（b）三维离散点图图5空间多级树管理物体的示意图。每行的矩形方块表示枝（左边框内，体，黑色水平虚线箭头表示属于同一枝的子枝间的链⑴建立树就是从一个物体出发寻找一个能装下物体的最小⑴建立树就是从一个物体出发寻找一个能装下物体的最小枝，其尺寸是最小分辨图6建立三角形树示意图（a）三角形物体，（b）成一个分辨率为σ的枝，（c）将的变成翻倍到能容纳三角形物体，（d）建立的三角形物体树箭头代表枝上悬挂的三角“/”代表树根，灰⑵向树上添加一个物体的思路是先扩展树以包含此物体，然后将物体挂到合适的枝上心所在象限产生当前枝的干s。将r连接为其子枝，将s设置为树根。(b)回到上一步(a)。7A和B为三角形物体，O1识，a是三级子枝的标识，B0B7-c（左边是枝的标识，右边是枝容纳的物体图树上添加物体示意图，(a)调整树根已容纳物体图树上添加物体示意图，(a)调整树根已容纳物体B，(b)将B物体挂到相应子枝，(c)应的空间多级树结构从树上删除一个物体8图从树上删除一个物体示意图，(a)找到需删除物体所在枝a，(b)删除a枝上物体，(c)除空枝a，(d)删除空枝I体所在象限，产生当前枝的子枝。通过调用这四个函数，对任意空间维数中任何形状物体2.2.2SHT的快速搜2.2.2SHT的快速搜条件搜结束。否则，出栈一枝b;③检查b的每个子枝，将子枝中的候选枝入栈；④检查b挂的物IVAB在圆内。这样,在各级树枝层次上避免搜索很多树枝中的物体,速减小.这样搜索计算量仅与多级树的总级别成正比平面三角形分布和空间划9速减小.这样搜索计算量仅与多级树的总级别成正比平面三角形分布和空间划9对应的空间多级快速搜索圆形区域中物体的过程极值搜极值搜距离A最近的点。快速搜索不是依次计算比较中每一点到A距离，而LBBBBo循环:计算oS，如果比候选值S，oLSS，oLSBZZZZ的上限比L某个值域枝C的下限小，删除C后面的所有同级值域枝；如果L中某个值域枝的上限比Z的下限小，放弃Z，否则将Z插入到链表L中。在插入过程中，链表L中的值域枝的排序,12平面点分布和空间划图 与图12对应的空间多级图 快速搜索距离A点最近点的过图 快速搜索距离A点最近点的过程计算一个物体的函数值，用3应空间多级树和基于空间多级树的快速搜索不仅提供了空间物体集合高效3应空间多级树和基于空间多级树的快速搜索不仅提供了空间物体集合高效通用的空间引，而且提供了数据有效管理手段。基于这种技术，大部分应用问题都可以用统一计算框处理。文中将带有物理属性的物体称为物元，比如带有能量、质量、动量、角动量属性的子、分子团，带有密度、相成分的颗粒，带有晶向的晶粒，带有密度、能量属性的有限大部分分析问题可以归结为从已有的物元集合构造另一种具有应用要求属性的其他物元合，例如离散动力学问题可以归结为：已有物元集合按照物理规则演变（物理体属性改物理体重新再构造）的过程。下面将分别介绍通用索引在最短路径问题、空间剖分、团簇构分析、原子类型鉴别以及界面构造等方面的应用最短路径问最短路径问题[9]是图论研究中的一个经典问题，旨在寻找图（由结点和路径组成）中结点之间的最短路最短路径不仅仅指一般地理意义上的最短距离引申到其它度量，如时间、费用、线路容量等。最短路径算法的选择与实现是通道路线设计的基础是计算机科学与地址信息科学等领域的研究热点，很多网络相关问题都可以纳入最短路问题的范畴之中。经典的图论与不断发展完善的计算机数据结构及算法的有效结合使得的最短路径算法不断涌现邮递员最短路径问题是给定N个城市，寻找一条最短的环路，可以遍历这些城市。这代表性的NP问题，因为最优的路径需要比较所有的环路。最为简单应用实例，此处给出种近似解称为迷宫。近似解的算法可以简单描述为：依次将点加入回路中，将其中一条边改造为过加入点的两条边线，改进后回路不自相交，且回路最短，如图15所示BPA图采用空间多级树实现算法的预备工作是设计快速搜索器。这里提供了包含三个极值索和一个条件搜索的设计方案，其中，每个极值搜索需要引入一个求值函数和图采用空间多级树实现算法的预备工作是设计快速搜索器。这里提供了包含三个极值索和一个条件搜索的设计方案，其中，每个极值搜索需要引入一个求值函数和两个评估数，条件搜索需要引入一个条件函数和一个判别函数⑴极小值快速搜索器MS1：在二维点树上搜索离给定固定点r0最近的点；求值函数value(p)rp，评估函数M(b)m(b)cbmax(cb2db别是枝的值域上、下⑵极小值快速搜索器MS2：在二维点树上搜索离给定线段AB两端相连的折线最短点；其求值函数value(P)rPrP，2db，那么评估函数形式如下A,BAC，BBC，A22drM(b)b 2drb rrBA，B ， ABA，Bm(b)cr ⑶极小值快速搜索器MS3：给定一点P，在二维线段树上搜索线段AB，使得它与线端点连接的折线减去原线段得到的值最小。求值函数为value(AB)rPrPrA;同上，按照枝的外接C评估，圆心为cb,半径为2db。评估函数形式如下M(b)2(2dbrPcb，Pm(b)max(c32d P Pm(b)max(c32d P b⑷条件快速搜索器CS1：给定线段AB，在二维线段树上搜索线段ab，使得它与给定ABxrA1x)rByra1y)rb，得到(xy0x0ycondition(ab)判别函数如下ABbmaycontain(b)至此，预备工作完毕。迷宫构造算法如下⑴初始用平面上给定离散点构造点树tp；从树上任意摘下一点P0；用MS1类搜索器在tpP3（P0P1）；构造线段P0-P1，用它建立线段树tl；构造线段P1-P2，将其挂在tl上构造线段P2-P0，将其挂在树tl上⑵生成连tp上摘下一个点P,s；②如果P是空的,结束程序；③用MS3tlLL1L1LtlLP线并tl上，回到步骤①；如果非空L堆栈s中,回到步骤③16是连40000个随机点构造出的迷宫，其中(a)是整体图，(b)是局域放大图图连接40000个随机点构造出的迷宫。(a)是整体图，(b)是局域放图连接40000个随机点构造出的迷宫。(a)是整体图，(b)是局域放大图Delaunay剖分的快速构造图 二维Delaunay剖分中加入一个新的三角形；(b)将这些三角形图 二维Delaunay剖分中加入一个新的三角形；(b)将这些三角形的内部线掉，保留外界线段；(c)将每一个外界线与相连，形成新的三角形⑴条件快速搜索CS1：给定点r0，在球接四面体树上搜索包含此点的球接四面体CTcondition(CT),半径为maycontain(b)cb ⑵条件快速CS2：给定三个P1P2P3，在球接四面体树上搜索以此三点为顶点P1,bCTP2,cCTcondition(CT)maycontain(b)图二维正方形[0,4]×[0,4]区域中随机分布的离散点构造的Delaunay剖分，(b)是中黑色小矩形区域的局部放大图三维球型区域中随机分布的离散点构造的Delaunay剖nn+1{r1,r2rn1构造外接球。球心坐标ccAB，球半径为(c图三维球型区域中随机分布的离散点构造的Delaunay剖nn+1{r1,r2rn1构造外接球。球心坐标ccAB，球半径为(cr2。1 2r11112r22，，A Bi210r210rnn0nn10团簇构造算法及其分析方应的空间多级树为t；②删除t上所有存在大于分辨率边长的四面体；构造单个团簇：③C->s上。⑤继续过程④，直到iC->s上。⑤继续过程④，直到i树为空。构造所有团簇：⑥将构造好的团簇C添加到团图二维空间[0,1]×[0,1]区域中1000个随机离散点计算得到的团簇结构，(a)是团簇（b）是相应的团簇边图三维空间中5000个随机离散点计算得到的团簇结图三维空间中5000个随机离散点计算得到的团簇结缺陷原子分3.4.1能量法及对称参数3.4.1能量法及对称参数s距的2倍，因此球半径可以选为最紧邻间距的13.4.2键对分析技结构分析的键对分析技术就可以更精密地鉴别缺陷类型。键对分析技术的原理是：根据两端原子的间接连接情况标志结构分析的键对分析技术就可以更精密地鉴别缺陷类型。键对分析技术的原理是：根据两端原子的间接连接情况标志键类型；根据一个原子的所有键的类型标志原子类型原子间距小于一个给定尺R（成键距离）的两个原子之间的连接定义为键。这里键的概念为了方便叙述借用的化学中一个名词，并没有电子波函数交叠等量子化学的涵。成键间距一般选为完整晶格最近间距的 倍。在预先不知道晶格常数的情况下，择径向分布函数第一峰对应的间距的 倍。这样选择使拓扑分析可以具有一定抗随机扰能力键类型标键类型取决于键两端原子的间接成键的特性。可以将与这两个原子都成键的原子选出来，用它们间的连接来标识键类型具体来讲可以用三个数字编码标识一个键的类型。设键两端原子为A,B；与A,B都键的原子构成一个集合c。第一c的元素数目；第二位表示c中原子中成键的目；第三位表示c中原子的最大配位数与A和B原子都成键的原子共有4个（H、I、C、D）4表示这些原子中成键的数DCEABFGIHJKLM图22键类型标识的示意完整晶体中键的标识很容易计算出来：对于fcc原子中的421；hcp晶体中的键661,200,100,311,211,411432542,300,400,原子类型标,661,200,100,311,211,411432542,300,400,原子类型标,键对分析算rorbrcmaycontain(b)其中，dbb立方体边长的一半&rprcrc&rprcrcdb&&rbrbrprc其中，dbb立方体边长的一半分布函数，根据径向分布函数的第一极大值对应的原子间距，给出成键距离。a，VCS1tpa（rc的原子）ba-b,CS2tpabrc的原子对mna-b100l+10m+nB的键类型比较，如果是新类型，存入键类型数组B，将键类型分布矢量V的相应的分量上aVAV子类型数组，记录a原子类型为V。效果展是在图采用键对分析技术计算得到的层错结构和位错结图采用键对分析技术计算得到的层错结构和位错结图采用键对分析技术计算得到的空洞表面和位错结界面构造算 由空间无序点构造物体表面的包裹-雕刻算由空间无序点构造物体表面是计算几何中的重要问题，目前的算法可以分为四类、变形法包裹算⑴极小值搜索器包裹算⑴极小值搜索器图搜索半面转动时最先遇到value(p)angle(x,其中xxˆrPrAyyˆrPrA为点的局部坐标Pxy(rArBrBˆrBPxy(rArBEˆˆ为投影算符，E为恒等算符ASorBM(b)，Eˆˆ为投影算符，E为恒等算符ASorBM(b)，rA,B/Sbb 0ABm(b)max(0,vaule(c)3d/P))A,Brb bcb是b，db⑵条件快速搜索P1,DP2BP1,P2①初始化：P1 图26包裹算法的初始化过程示意图。(a)绕P1P2旋转碰到Q1点绕P1P1旋转碰到Q2，(c)绕Q1Q2旋转碰到Q3，(d)初始三角形界②形成P1 图26包裹算法的初始化过程示意图。(a)绕P1P2旋转碰到Q1点绕P1P1旋转碰到Q2，(c)绕Q1Q2旋转碰到Q3，(d)初始三角形界②形成包从三角形边树tl摘下一个三角形边L，L为空，结束程序，否则找到对应的三形。用MS1搜索器在点树tp上搜索三角形绕L旋转最先碰到的点P,L的顶La、LbP构造三角形LbLaP，。LbLaP挂到三角形树ttCS1搜索器在三角形tb索与PLa共边的三角形边，如果找到（见图27-a）tb上摘下，如果找不到（见图b）LaP生成三角形边，挂到三角形边树上。对LaP采用同样操作PPPP图27两种情况。(a)LaP不在边树上，(b)LaP在边树雕刻算雕刻算法涉及到PPPP图27两种情况。(a)LaP不在边树上，(b)LaP在边树雕刻算雕刻算法涉及到雕刻方法和雕刻标准。雕刻过程中要保证没有点在雕刻时被雕刻掉而且雕刻后表面尽量光滑。光滑的表面曲率很小，对应于雕刻掉的四面体的外接球半径大。对于待雕刻的三角形ABC来说，雕刻过程意味着寻找当前面包围区域内部点P，使角形顶点ABC和P成的球冠具有最小高度28所示对一个三角形面来说，雕刻的停止条件是外接球的半径小于给定值。对于表面均匀样的物体这个判据符合要求，而对于非均匀采样（曲率大，采样点密）的表面构造，应保证曲率与局部采样点间距相关。采用球冠高与三角形外接圆的比值，即雕刻ch/判据比较合适，这个判据能同时适应两种情形。雕刻算法是图28雕刻算法的原(rrrr (rB图28雕刻算法的原(rrrr (rBrA)(rCrA(r，n r 给出ABC外接圆半径r，圆心o，和法线方向n；通过求解方rcro(rr)22 r value(p) 2rb)MMM 2r0m(b)mm，rn2r2 其中，2 obbmP点为Q2点时，球冠具有最大球冠高M(rBrA)(rCrA图29b枝的外接球与球冠ABCP相切的两种情况。图29b枝的外接球与球冠ABCP相切的两种情况。其中，黑色圆代表b的外接球，蓝色直线代表三角形ABC的外接圆，绿色和红色圆代表球冠ABCP，Q1和Q2为交点雕刻算法如下tt。设置雕刻结束临界值ccri2.0，设置界面树tt0为空。否则，用MS搜索器查找最浅的点p，计算雕刻度c，判断是否满足cccri，如果满足，不雕刻，将三角形T挂在界面树tt0上。如果不满足，将三角形T的每条边与点p分别③回到步骤②界面tt0上所有的表面构成物体的表面算法展30fcc晶体中两个纳米空洞和一个点缺陷边界上的均匀点集由包裹算法给出的壳及雕刻出的表面。图31是由多连通单个物体表面随机采样点构造出的包裹凸壳和雕刻的表面。图32是由多个物体表面随机采样点构造出的表面和包裹面图30包裹-雕刻算法过程，（a）匀离散点（b）包裹过程包图30包裹-雕刻算法过程，（a）匀离散点（b）包裹过程包裹凸壳(d)雕刻过程(e)表（a）均匀离散点（b）包裹凸壳(e)表图31由两个物体随机采样9000点重构的表（a）均匀离散点（b）包裹凸(e)表图32由一个复杂物体随机采样点重构的表 滚球法由空间无序点寻找物理区域界滚球法算⑴构造极值搜索器给定定向三角形ABC和 滚球法由空间无序点寻找物理区域界滚球法算⑴构造极值搜索器给定定向三角形ABC和指定的AB，在点树上搜索三角形上方的半径为r的球旋转时遇到的第一个点。求值函数构造过程如下：计算初始旋转球的球心位置ro及局rrroAP(rrroA,ˆ ,ˆr)2，oBr 2roC计算旋转后的球心位置rn及在局部坐标中的坐标x,y,z，具体公式Pxy(rorAr2(r r2(rnrB)2，xxˆ(rnrAr2(r r2(rnrB)2，xxˆ(rnrA),yyˆr2(r value(p)angle(x,y图33三角形ABC的旋转示评估函数构造过程如下：按照枝的外接球与旋转球（球心在）相切计算切点位置22,rrrrrr3d 1B b此方程有两个根rM1,rm1,对应空间两个点M1和M(b)三角形外接圆在外接球，其他情m(b)三角形外接圆在外接球其他情绿色和红色为滚球，M1m1M1旋转后对应最小旋转角的滚旋转前的滚旋转旋转后对应最大旋转角的滚图滚球与b的外旋转后对应最小旋转角的滚旋转前的滚旋转旋转后对应最大旋转角的滚图滚球与b的外接球相切的两种情况示意图，图中黑色圆点代表旋转轴⑵构造极值搜索MS2：给定一个点，在点树上搜索离此点最近的点。求值函数为value(p)，其中rpp的坐标；评估函数为M(b)cb3db,m(b)max(cb 3db,0)其中cbb的中心坐标dbb方体半宽⑶构造极值搜索MS3：给定一个点和转轴，在点树上搜索大小固定的球绕指定转旋转时遇到的第一个点。算法与MS1基本相同，此处不再重复⑷构造条件搜索器CS1:给定两个p1,p2，在边树上搜索对应的边BD，p1和p2的端点。条件函数为Bp1,Dp2按照枝的外接球S进行评估，判别函数如下 p1,p2滚球法算法如下①初始设定滚球的半径，球心初始P0。用给定的离散点构造点树tp;用搜索器MS2搜②滚球法构造表面算法效果展图35由离散空间点构造空洞表面效果(a)给定离散点，(b)用滚球法给出图36由离散空间点构造空图36由离散空间点构造空洞表面的空间多级树通用索引C++类结构及实数据结构利c++面向对象的优势，根据通用性要求使类的继承结构从抽象到具体分呈现，树的动态构造在抽象类中实现，快速搜索在形式类中实现。这种处理使得应用类的A1空间多级树图A-1A1空间多级树图A-1树枝的继承结抽象枝类Branchclass{intBranch*father; Branch*son; Branch*brother;Object空间枝类2{Branch2D(intlevel,doublesz,doublexxc,doubleyyc);doublexc,double抽象物体类class{voidbuild(Branch抽象物体类class{voidbuild(Branch*&tree,double //由当前物体以dd的分辨率构造树booladd_to(Branch*&Object*next_obj;//virtualBranch*create_branch_include_me(doublemsize)=0//产生一个尺寸为msize*2m的枝,//以枝b为参考位置，在当前物体所在象限位置产生b的干//以枝b在象限位置产生bvirtualboolin_branch(Branch*b)=0；//判断当前物体是否包含于枝b},形式物体类形式物体类的数据结构含有空间参数（位置，尺寸{virtualBranch*create_branch_uplevel(Branch*);virtualBranch*create_branch_uplevel(Branch*);virtualboolin_branch(Branch*);classBody3D:public{Body3D(doublexx,doubleyy,doublezz,doublexs,doubleys,doublezs)；virtualBranch*create_branch_uplevel(Branch*);virtualboolin_branch(Branch*);double形式物体类在实现四个虚函数create_branch_include_me{A2快速搜索class{voidscreen(Objectclass{voidscreen(Object voidObject*search_obj()=0;BranchintObject条件搜索类{virtualObject*search_obj();virtualboolcondition(Object*o)=0;intBranchBranchObject应用条件搜索类{virtualboolcondition(Object*o);virtualboolmaycontain{virtualboolcondition(Object*o);virtualboolmaycontain(Branch*b);InCirTet3D_Searcher搜索类中condition函数判别给定球接四面体是否包含点{virtualObject*search_obj();virtualvoidmin_max_value(Branch*b,double&min,double&max)=0;virtualdoublevalue(Object*o)=0;Object*desired; doublemvalue; …应用极值搜索类{virtualvoidmin_max_value(Branch*b,double&min,double&max);virtualdoublevalue(Object*o); A3空间多级树及快速搜索的使用方A3空间多级树及快速搜索的使用方classAtom3D:publicPoint3Ddoublevx,vy,vz, Brancha-…a-//将物体挂在树treettreeBranch*ttree=0;…//将物体挂在树tree使用快速搜索器搜索物体的方法：通过继承形式搜索器（Condition_Searcher对于conditionmay_be_contained_inMinimum_Searchervaluemax_min函数下面程序段通过继承Condition_Searcher定义了一个应用条件搜索器 bool{下面程序段通过继承Condition_Searcher定义了一个应用条件搜索器 bool{Bool A_Condition_Searcher*cs=new 下面程序段通过继承Minimum_Searcher定义了一个应用极值搜索器doubleA_Minimum_SearchervoidA_Minimum_Searcher::max_min(Branch*b,double&max,A_Minimum_Searcher*ms=new//initializethesearch参考文DierkRaabe,ComputationalMaterialsScience:TheSimulationofMaterialsMicrostructuresandProperties,Wiley-VCH1998.K.S.Kumar,H.VanSwygenhoven,S.Suresh,Mechanicalbehaviorofnanocrystallinemetalsandalloys-ActaMaterialia51(2003)5743–5774.Kelchner,Plimpton,Hamilton,PhysRevB,58,11085Ackland,Jones,PhysRevB,73,054104Faken,Jonsson,ComputMaterSci,2,279G.Michael，K.Stephan，Z.Gerhard，“NumericalSimulationinMolecularDynamics”，Springer-VerlagBerlinHeidelberg2007.MarkdeMarkdeBerg,MarcvanKreveld,MarkOvermars,andOtfriedSchwarzkopf,ComputationalGeometry(2ndrevisededitioned.),Springer-Verlag,2000.DonaldKnuth,SortingandSearching,TheArtofComputerProgramming,Volume3(Thirded.),Addison-Wesley,1997.Moore,E.F.(1959)."Theshortestpaththroughamaze".ProceedingsofanInternationalSymposiumontheTheoryofSwitching(Cambridge,Massachusetts,2–5April1957).Cambridge:HarvardUniversityPress.pp.deBerg,Mark;OtfriedCheong,MarcvanKreveld,MarkOvermars(2008).ComputationalGeometry:AlgorithmsandApplications.Springer-Verlag.ISBN978-3-540-77973-5.AmentP.EffectsofionizationgradientsonInertial-Confinement-Fusioncapsulehydrodynamicsstability.PhysRevLett,2008,101:115004deVries1PC,HuaMD,McDonaldDC,etal.ScalingofrotationandmomentumconfinementinJETplasmas.NuclFusion,2008,48:065006BowlerBP,WallerWH,MegeathST,etal.Aninfraredcensusofstarformationinthehorseheadnebula.AstronJ,2009,137:3685—3699HernquistL.HierarchicalN-bodymethods.ComputPhysCommun,1988,48:107—MakinoJ.Vectorizationofatreecode.JComputPhys,1990,87:148—HidakaY,ChoiEM,MastovskyI,etal.Observationoflargearraysofplasmafilamentsinairbreakdownby1.5-MW110-GHzgyrotronpulses.PhysRevLett,2008,100:035003NogaretT,RodneyD,FivelM,etal.Clearbandformationsimulatedbydislocationdynamics:Roleofhelicalturnsandpile-ups.JNuclMater,2008,380:22—29KotsiantisSB,PintelasPE.Recentadvancesinclustering:Abrief