分子动力学中晶体位错的三维可视化系统设计与实现

宁媛　张鹏　介军

摘要：在材料科学中，位错线是晶体或晶格内滑移面上已滑动区的边界，位错线的识别对于鉴定晶体的性能有着重要的影响。该文将三维可视化技术应用到分子动力学领域，改进了原有的只能通过计算结果分析数据、通过二维图像简单观察数据的方法，设计了晶体位错可视化系统，使得晶体中各原子之间形成的三维结构易于观察，便于研究晶体位错等重要的材质特性，提高位错线的识别效率。

关键词：位错线识别；三维可视化；骨架提取算法

中图分类号：TP311 文献标识码：A 文章编号：1009-3044（2014）16-3896-04

Abstract： In the material science， the dislocation lines are the boundary sliding of a crystal or a lattice slip plane，the identification of dislocation lines have a important impact for the property of crystals. The essay using the application of 3D visualization technology to the field of molecular dynamics and improving the conventional research ways， such as analyzing data by calculating results and observing the data by 2D images inconveniently， designing a visualization system for crystal dislocations. This way can make the 3D structure among atoms be observed easily， which is convenient to research the crystal dislocations and other important material properties， and improve the efficiency of the identification of dislocation lines.

Key words： the identification of dislocation lines； three-dimensional visualization； skeleton extraction algorithm

位错理论是金属力学性能的理论基础， 是材料科学研究中所不可缺少的。1934年，科学家泰勒· 奥罗万及波兰伊几乎同时提出了位错的假设， 即认为晶体中存在有一种线缺陷， 这一概念的提出成功地解释了以前所无法解决的现象[1]。位错线的形成和发展可用弗兰克-里德位错源（Frank-Read source）增值的原理解释[2]：当应力超过临界应力时，位错线扩张，形成内外两部分，外部位错逐渐扩大，内部位错线恢复原状，在外力作用下，不断产生新位错环，因而得到很大的滑移量。一个位错从材料内部运动到了材料表面，就相当于其位错线扫过的区域整体沿着该位错的伯格斯矢量方向[3]滑移了一个单位距离。在对材料进行“冷加工”时，其内部的位错密度会因为位错的萌生与增值机制的激活而升高。因此，从微观角度来看，位错线的识别对鉴定晶体的性能具有很大帮助。

目前，骨架提取算法的研究可以分为四大类：基于拓扑与几何分析的方法、基于拓扑细化的方法、基于距离场的方法以及广义势场方法[4]。综合以上四类算法的利弊，结合所研究数据对象的特点，该文所采用的改进的骨架提取算法有以下几个方面的优势：

1）能有效处理点云数据，适合离散模型，具有很高的执行效率；

2）能有效处理网格数据，有很好的连通性，对边界噪声不敏感；

3）能有效处理封闭的边界网格，易于实现，能较好的保证原晶体模型的拓扑结构。

本文所研究的对象为晶体的原子结构，通过IMD[5]计算出原子的三维坐标以及活动参数。例如原子的质量，速度，势能等。该算法能够直接利用原子数据模型，独立于特定的晶体结构，并且能够维持最大的信息数量，保证晶体拓扑结构的完整，从而很好的识别出位错线。

2 算法描述

本文所使用的算法是基于拓扑细化方法[4]和网格收缩法[6]的思想并结合晶体位错原子数据的特点采用了反距离插值法（IDW）[7]。通过读入晶体位错的原子数据显示出三维晶体，然后计算原子之间的距离与晶体中临近原子的关系构造出初始晶体的网格结构，再利用网格收缩的思想将已构造的网格结构收缩成折线，经过后期处理，从而得到位错线，其中线与线之间的交界处可以理解为位错的节点。

算法的核心步骤如下所述[8]：

1）构造初始网格，将原子替代为点，原子与邻近原子之间的欧几里得距离若小于该晶体完美最近距离的1.2倍时，在它们之间创建边（线）。

2）在不改变晶体网状拓扑结构的前提下，通过移动相邻节点使其相互接近收缩网格，当网格收缩为一个由关键节点组成的折线时，停止收缩。

3）如果两个节点之间的距离下降到低于某一临界阈值（例如最近邻原子距离的十分之一）时，合并两个节点以降低网格的复杂性。

4）分析处理由于收缩、合并所形成的折线，除去闭合的位错环路，最终达到骨架化。

算法的具体过程如图1所示：

2 系统设计与实现

2.1 系统结构设计

整个系统的设计是以上述的改进的骨架提取算法为主线再结合三维技术，调用三维图形库OpenGL 辅助实现的。在设计和实现上以低成本、扩展性好、易操作、易观察、专业性强等为特点，在应用功能上以较高的实用性，快速渲染，真实表达位错线特点为目的，为分子动力学领域中的位错晶体识别做了更好的诠释与演绎。endprint

晶体位错识别系统由六部分组成，如图2所示。首先通过读入由IMD计算好的数据文件显示出原晶体模型，再根据限制条件对读入的数据进行分析与处理，然后进一步建立收缩网格，当网格收缩成一个满足提取条件的位错线时，显示出最终的提取结果。

2.2 构造初始网格

显示出的晶体模型需要建立初始网格来为后续的骨架化过程做铺垫。由原子和近邻原子之间的关系创建初始网络（无向图）M =（V，E），为每一个原子创建一个节点，定义一组节点V，若每两个节点之间的距离小于一个给定的阈值时（通常为完美最近距离的1.2倍），在它们之间创建边（线），限定边的集合为E。

使用欧几里得距离公式计算每两个原子之间的距离：

[d=（x1-x2）2+（y1-y2）2+（z1-z2）2] （1）

其中x，y，z为原子的三维坐标。初始网格的建立过程如图3所示：

2.3 网格收缩

在不改变网格拓扑结构的前提下，通过移动中心原子的临近原子使原子相互接近实现网格收缩。这一过程不受约束，直到网格收缩到一个单一的点时停止。这时网格已局部转换成一个骨架。其中的中心原子（关键节点）通过如下所述的局部拓扑连接来确定：

设V'为从v节点出发且可以通过任意一条小于等于2的路径到达的节点的集合。当且仅当从v到v'有一条唯一确定的小于等于2的路径时，定义v为关键节点，即中心原子，其中v'[∈]V'。这一定义是为防止网格中任意四条或少于四条边形成环路，这些环路将会影响网格收缩和后来的节点合并，因此需要被删除。相反，若无进一步的收缩，任一节点和它们的邻边所形成的一个没有环路的连续路径即为位错线。

由于在网格收缩的过程中，我们需要中心原子周围的临近原子能够快速向其移动，而离中心原子较远的原子会由于距离的增大而作用力减小，因此采用反距离插值权重法（IDW）实现网格收缩，该方法主要依赖于反距离的幂值，幂参数可基于距输出点的距离来控制已知点对插值的影响，离散点距离插值点越远，其对插值点的影响力越低[7]。每个非关键节点v[∈]V从当前位置[cv]更新到新的位置[cv']，根据下式考虑所有非关键节点[Nv]：

[cv'=cv+qn∈Nv（cn-cv）/（||cn-cv||2）m∈Nv1/（||cm-cv||2）] （2）

其中因子[q]控制普遍收缩率，当[q]= 0.33时是收缩速度和精度之间的平衡值。

网格收缩的具体过程如图4所示：

2.4 合并节点

为了减少网格的复杂性，若两个节点之间的距离下降到低于某一临界阈值时（例如最近邻距离的十分之一），将两个节点合并在一起。为防止网格中拓扑结构的修改，当且仅当点[va]和[vb]中至少有一个点为非关键节点且它们共享第三个节点，并且它们所共享的节点能够通过小于等于2的路径到达时，合并这两个节点。这些要求能够防止当前网格的管道或环路塌陷，但不能保证在一系列的收缩和合并的过程中不改变拓扑结构。合并节点本身也是通过边收缩的方法来实现的。边（[va]，[vb]）从网格E中删除，[vb]从一组节点V中删除，任意节点vx与节点vb组成的边（vx，vb）由（vx，va）替代。在一个迭代中，若没有节点移动和合并则骨架化进程终止。节点合并的过程如图5所示：

2.5 后期处理

构造好的初始网格经过收缩和合并后，网格已经被矢量化成为以折线表示的骨架。此时，需要考虑两种情况进行修复处理：

1）位错开始和结束于位错结点、自由表面、或者形成一个闭和的环路。位错的起始和结束的交汇点可以认为是一个节点周围有不少于两个的邻近节点在自由表面上通过一个连结点连接。这两个节点之间的位错线以折线表示。在这种情况下，如果有任意不属于这些路径的边但属于其它闭合的回路时，也可以被提取。

2）所提取的位错线不一定在它们的位错中心处，因为一个节点的位置在骨架化中可能会有显著的改变。为达到拓扑结构稳定，每个偏离微小位置的节点归并到已经被合并的所有节点的坐标中去。

3 总结

读入的晶体模型数据首先通过调用OpenGL图形库显示出三维晶体模型，然后将晶体中的原子以点代替，原子间的关系以线代替构造出晶体网格结构。遍历网格中所有的点，对每一个点进行判断，若此点为非关键节点，则移动此点使网格收缩，当任意两个节点之间的距离达到临界阈值时，合并这两个节点。当所有节点不再移动时，对所形成的位错折线进行后期处理，使其变为光滑的曲线，即为最后的位错线。

本文与三维可视化技术相结合使用了一种可用于晶体位错线识别的骨架提取方法，并结合该方法设计并实现了晶体位错三维可视化系统，为材料科学中的位错线识别提供了更直观的帮助。在实际工作中，人们可以根据位错曲线方便的找出位错晶体的分布，从而更好的研究位错晶体的结构特征。

参考文献：

[1] 沈寅忠.位错理论的形成与发展[J].金属世界，1995（2）.

[2] 余勇，潘晓霞.Frank-Read 位错源的细观级模拟[J].金属学报，2009，45（11）：1309-1313.

[3] 石德珂.材料科学基础[M].北京：机械工业出版社，2003.

[4] 彭艺.三维模型骨架提取算法研究与实现[D].上海：华中科技大学，2008.

[5] J.ROTH.IMD-A Molecular Dynamics Program and Applications[J].World Scientific，1999，1-12.

[6] Au OKC， Tai CL， Chu HK， Or DC， Lee TY. Skeleton Extraction by Mesh Contraction[J]. 2008.

[7] 刘光孟，汪云甲，王允.反距离权重插值因子对插值误差影响分析[J].中国科技论文在线，2010，5（11）：879-884.

[8]Begau C，Hartmaier A，George E P，et al.Atomistic processes of dislocation generation and plastic deformation during nanoindentation[J].Acta Materialia，2011，59（3）：934-942.

[9] Ackland G J，Jones A P.Phys Rev B Applications of local crystal structure measures in experiment and simulation [J].2006，73（5）：54-104.endprint

晶体位错识别系统由六部分组成，如图2所示。首先通过读入由IMD计算好的数据文件显示出原晶体模型，再根据限制条件对读入的数据进行分析与处理，然后进一步建立收缩网格，当网格收缩成一个满足提取条件的位错线时，显示出最终的提取结果。

2.2 构造初始网格

显示出的晶体模型需要建立初始网格来为后续的骨架化过程做铺垫。由原子和近邻原子之间的关系创建初始网络（无向图）M =（V，E），为每一个原子创建一个节点，定义一组节点V，若每两个节点之间的距离小于一个给定的阈值时（通常为完美最近距离的1.2倍），在它们之间创建边（线），限定边的集合为E。

使用欧几里得距离公式计算每两个原子之间的距离：

[d=（x1-x2）2+（y1-y2）2+（z1-z2）2] （1）

其中x，y，z为原子的三维坐标。初始网格的建立过程如图3所示：

2.3 网格收缩

在不改变网格拓扑结构的前提下，通过移动中心原子的临近原子使原子相互接近实现网格收缩。这一过程不受约束，直到网格收缩到一个单一的点时停止。这时网格已局部转换成一个骨架。其中的中心原子（关键节点）通过如下所述的局部拓扑连接来确定：

设V'为从v节点出发且可以通过任意一条小于等于2的路径到达的节点的集合。当且仅当从v到v'有一条唯一确定的小于等于2的路径时，定义v为关键节点，即中心原子，其中v'[∈]V'。这一定义是为防止网格中任意四条或少于四条边形成环路，这些环路将会影响网格收缩和后来的节点合并，因此需要被删除。相反，若无进一步的收缩，任一节点和它们的邻边所形成的一个没有环路的连续路径即为位错线。

由于在网格收缩的过程中，我们需要中心原子周围的临近原子能够快速向其移动，而离中心原子较远的原子会由于距离的增大而作用力减小，因此采用反距离插值权重法（IDW）实现网格收缩，该方法主要依赖于反距离的幂值，幂参数可基于距输出点的距离来控制已知点对插值的影响，离散点距离插值点越远，其对插值点的影响力越低[7]。每个非关键节点v[∈]V从当前位置[cv]更新到新的位置[cv']，根据下式考虑所有非关键节点[Nv]：

[cv'=cv+qn∈Nv（cn-cv）/（||cn-cv||2）m∈Nv1/（||cm-cv||2）] （2）

其中因子[q]控制普遍收缩率，当[q]= 0.33时是收缩速度和精度之间的平衡值。

网格收缩的具体过程如图4所示：

2.4 合并节点

为了减少网格的复杂性，若两个节点之间的距离下降到低于某一临界阈值时（例如最近邻距离的十分之一），将两个节点合并在一起。为防止网格中拓扑结构的修改，当且仅当点[va]和[vb]中至少有一个点为非关键节点且它们共享第三个节点，并且它们所共享的节点能够通过小于等于2的路径到达时，合并这两个节点。这些要求能够防止当前网格的管道或环路塌陷，但不能保证在一系列的收缩和合并的过程中不改变拓扑结构。合并节点本身也是通过边收缩的方法来实现的。边（[va]，[vb]）从网格E中删除，[vb]从一组节点V中删除，任意节点vx与节点vb组成的边（vx，vb）由（vx，va）替代。在一个迭代中，若没有节点移动和合并则骨架化进程终止。节点合并的过程如图5所示：

2.5 后期处理

构造好的初始网格经过收缩和合并后，网格已经被矢量化成为以折线表示的骨架。此时，需要考虑两种情况进行修复处理：

1）位错开始和结束于位错结点、自由表面、或者形成一个闭和的环路。位错的起始和结束的交汇点可以认为是一个节点周围有不少于两个的邻近节点在自由表面上通过一个连结点连接。这两个节点之间的位错线以折线表示。在这种情况下，如果有任意不属于这些路径的边但属于其它闭合的回路时，也可以被提取。

2）所提取的位错线不一定在它们的位错中心处，因为一个节点的位置在骨架化中可能会有显著的改变。为达到拓扑结构稳定，每个偏离微小位置的节点归并到已经被合并的所有节点的坐标中去。

3 总结

读入的晶体模型数据首先通过调用OpenGL图形库显示出三维晶体模型，然后将晶体中的原子以点代替，原子间的关系以线代替构造出晶体网格结构。遍历网格中所有的点，对每一个点进行判断，若此点为非关键节点，则移动此点使网格收缩，当任意两个节点之间的距离达到临界阈值时，合并这两个节点。当所有节点不再移动时，对所形成的位错折线进行后期处理，使其变为光滑的曲线，即为最后的位错线。

本文与三维可视化技术相结合使用了一种可用于晶体位错线识别的骨架提取方法，并结合该方法设计并实现了晶体位错三维可视化系统，为材料科学中的位错线识别提供了更直观的帮助。在实际工作中，人们可以根据位错曲线方便的找出位错晶体的分布，从而更好的研究位错晶体的结构特征。

参考文献：

[1] 沈寅忠.位错理论的形成与发展[J].金属世界，1995（2）.

[2] 余勇，潘晓霞.Frank-Read 位错源的细观级模拟[J].金属学报，2009，45（11）：1309-1313.

[3] 石德珂.材料科学基础[M].北京：机械工业出版社，2003.

[4] 彭艺.三维模型骨架提取算法研究与实现[D].上海：华中科技大学，2008.

[5] J.ROTH.IMD-A Molecular Dynamics Program and Applications[J].World Scientific，1999，1-12.

[6] Au OKC， Tai CL， Chu HK， Or DC， Lee TY. Skeleton Extraction by Mesh Contraction[J]. 2008.

[7] 刘光孟，汪云甲，王允.反距离权重插值因子对插值误差影响分析[J].中国科技论文在线，2010，5（11）：879-884.

[8]Begau C，Hartmaier A，George E P，et al.Atomistic processes of dislocation generation and plastic deformation during nanoindentation[J].Acta Materialia，2011，59（3）：934-942.

[9] Ackland G J，Jones A P.Phys Rev B Applications of local crystal structure measures in experiment and simulation [J].2006，73（5）：54-104.endprint

晶体位错识别系统由六部分组成，如图2所示。首先通过读入由IMD计算好的数据文件显示出原晶体模型，再根据限制条件对读入的数据进行分析与处理，然后进一步建立收缩网格，当网格收缩成一个满足提取条件的位错线时，显示出最终的提取结果。

2.2 构造初始网格

显示出的晶体模型需要建立初始网格来为后续的骨架化过程做铺垫。由原子和近邻原子之间的关系创建初始网络（无向图）M =（V，E），为每一个原子创建一个节点，定义一组节点V，若每两个节点之间的距离小于一个给定的阈值时（通常为完美最近距离的1.2倍），在它们之间创建边（线），限定边的集合为E。

使用欧几里得距离公式计算每两个原子之间的距离：

[d=（x1-x2）2+（y1-y2）2+（z1-z2）2] （1）

其中x，y，z为原子的三维坐标。初始网格的建立过程如图3所示：

2.3 网格收缩

在不改变网格拓扑结构的前提下，通过移动中心原子的临近原子使原子相互接近实现网格收缩。这一过程不受约束，直到网格收缩到一个单一的点时停止。这时网格已局部转换成一个骨架。其中的中心原子（关键节点）通过如下所述的局部拓扑连接来确定：

设V'为从v节点出发且可以通过任意一条小于等于2的路径到达的节点的集合。当且仅当从v到v'有一条唯一确定的小于等于2的路径时，定义v为关键节点，即中心原子，其中v'[∈]V'。这一定义是为防止网格中任意四条或少于四条边形成环路，这些环路将会影响网格收缩和后来的节点合并，因此需要被删除。相反，若无进一步的收缩，任一节点和它们的邻边所形成的一个没有环路的连续路径即为位错线。

由于在网格收缩的过程中，我们需要中心原子周围的临近原子能够快速向其移动，而离中心原子较远的原子会由于距离的增大而作用力减小，因此采用反距离插值权重法（IDW）实现网格收缩，该方法主要依赖于反距离的幂值，幂参数可基于距输出点的距离来控制已知点对插值的影响，离散点距离插值点越远，其对插值点的影响力越低[7]。每个非关键节点v[∈]V从当前位置[cv]更新到新的位置[cv']，根据下式考虑所有非关键节点[Nv]：

[cv'=cv+qn∈Nv（cn-cv）/（||cn-cv||2）m∈Nv1/（||cm-cv||2）] （2）

其中因子[q]控制普遍收缩率，当[q]= 0.33时是收缩速度和精度之间的平衡值。

网格收缩的具体过程如图4所示：

2.4 合并节点

为了减少网格的复杂性，若两个节点之间的距离下降到低于某一临界阈值时（例如最近邻距离的十分之一），将两个节点合并在一起。为防止网格中拓扑结构的修改，当且仅当点[va]和[vb]中至少有一个点为非关键节点且它们共享第三个节点，并且它们所共享的节点能够通过小于等于2的路径到达时，合并这两个节点。这些要求能够防止当前网格的管道或环路塌陷，但不能保证在一系列的收缩和合并的过程中不改变拓扑结构。合并节点本身也是通过边收缩的方法来实现的。边（[va]，[vb]）从网格E中删除，[vb]从一组节点V中删除，任意节点vx与节点vb组成的边（vx，vb）由（vx，va）替代。在一个迭代中，若没有节点移动和合并则骨架化进程终止。节点合并的过程如图5所示：

2.5 后期处理

构造好的初始网格经过收缩和合并后，网格已经被矢量化成为以折线表示的骨架。此时，需要考虑两种情况进行修复处理：

1）位错开始和结束于位错结点、自由表面、或者形成一个闭和的环路。位错的起始和结束的交汇点可以认为是一个节点周围有不少于两个的邻近节点在自由表面上通过一个连结点连接。这两个节点之间的位错线以折线表示。在这种情况下，如果有任意不属于这些路径的边但属于其它闭合的回路时，也可以被提取。

2）所提取的位错线不一定在它们的位错中心处，因为一个节点的位置在骨架化中可能会有显著的改变。为达到拓扑结构稳定，每个偏离微小位置的节点归并到已经被合并的所有节点的坐标中去。

3 总结

读入的晶体模型数据首先通过调用OpenGL图形库显示出三维晶体模型，然后将晶体中的原子以点代替，原子间的关系以线代替构造出晶体网格结构。遍历网格中所有的点，对每一个点进行判断，若此点为非关键节点，则移动此点使网格收缩，当任意两个节点之间的距离达到临界阈值时，合并这两个节点。当所有节点不再移动时，对所形成的位错折线进行后期处理，使其变为光滑的曲线，即为最后的位错线。

本文与三维可视化技术相结合使用了一种可用于晶体位错线识别的骨架提取方法，并结合该方法设计并实现了晶体位错三维可视化系统，为材料科学中的位错线识别提供了更直观的帮助。在实际工作中，人们可以根据位错曲线方便的找出位错晶体的分布，从而更好的研究位错晶体的结构特征。

参考文献：

[1] 沈寅忠.位错理论的形成与发展[J].金属世界，1995（2）.

[2] 余勇，潘晓霞.Frank-Read 位错源的细观级模拟[J].金属学报，2009，45（11）：1309-1313.

[3] 石德珂.材料科学基础[M].北京：机械工业出版社，2003.

[4] 彭艺.三维模型骨架提取算法研究与实现[D].上海：华中科技大学，2008.

[5] J.ROTH.IMD-A Molecular Dynamics Program and Applications[J].World Scientific，1999，1-12.

[6] Au OKC， Tai CL， Chu HK， Or DC， Lee TY. Skeleton Extraction by Mesh Contraction[J]. 2008.

[7] 刘光孟，汪云甲，王允.反距离权重插值因子对插值误差影响分析[J].中国科技论文在线，2010，5（11）：879-884.

[8]Begau C，Hartmaier A，George E P，et al.Atomistic processes of dislocation generation and plastic deformation during nanoindentation[J].Acta Materialia，2011，59（3）：934-942.

[9] Ackland G J，Jones A P.Phys Rev B Applications of local crystal structure measures in experiment and simulation [J].2006，73（5）：54-104.endprint