CN101266627A - 一种检测计算机仿真结果中晶粒的取向和结晶程度的方法 - Google Patents

Abstract

本发明公开了一种检测计算机仿真结果中晶粒的取向和结晶程度的方法：将计算机仿真结果中待检测元素的所有原子的坐标读入计算机，按照设定的采样频率，沿某坐标轴的各采样点统计横轴原子密度分布函数，对其进行傅立叶变换，得到振幅－频率图，该图可反映大范围内的晶体取向、结晶程度；再通过分析归一化横轴原子密度分布研究局部的晶体取向和结晶程度。由于将三维的信息化简为一维，并且可以直接从“振幅－频率图”判断晶体的取向，使晶体取向信息得到充分表征，降低了判断晶体取向的难度。由于大部分步骤都可以编程达到自动化，极大提高了判断晶体取向的效率。表征的结果对选取的角度非常敏感，判定晶体取向的精确度很高。

Description

一种检测计算机仿真结果中晶粒的取向和结晶程度的方法

技术领域

本发明属于信息处理、计算机软件和计算材料科学的技术领域，具体涉及一种在大规模分子动力学仿真计算中检测晶体取向并判定结晶程度的方法。

背景技术

大规模分子动力学仿真是研究微纳尺度物质性质的有效方法。随着计算机技术的迅速提高，目前的模拟规模可以达到上亿个原子，相当于一个约100纳米的立方体。但是，如何表征所模拟的样品仍然是一个较为复杂的问题。一般而言，对于一个待模拟的单晶或多晶样品，经过驰豫(即不施加外力，令样品中的各原子自由运动)、拉伸等过程后，单晶样品可能转变为多晶甚至非晶态样品，原始的晶向也可能会发生变化。通常采用的键对分析(或称为“公共近邻分析”)可以检测出样品中的结晶区域和非结晶区域(J DanaHoneycutt，Hans C Andersen，Molecular Dynamics Study of Melting and Freezing of SmallLennard-Jones Clusters，Journal of Physical Chemistry，1987，91，4950-4963)不过目前尚无一个较好的、可以通过计算机自动实现的办法能够检测出经过模拟后的样品中各个晶粒(即被非结晶区域所分隔的各个结晶区域)的取向。如果分子动力学软件的显示功能较好，则有经验的研究者可以通过观察样品中的原子排布确定晶粒取向，但是这样费时费力，准确度有限，对于上百纳米的样品而言几乎没有可行性。

发明内容

本发明旨在针对大规模计算机仿真中的可能的周期性晶体结构，提出一种检测其中晶粒的晶体取向并判定其结晶程度的方法。该方法既可以检测局部的结晶情况，又可以反映大范围内的整体有序程度。

本发明所述的检测计算机仿真结果中晶粒的取向和结晶程度的方法，按以下步骤进行：

1)对于计算机仿真(如大规模分子动力学仿真)后得到的结果，将其中待检测元素的所有原子的坐标读入计算机中(如果仿真体系只有一种元素，就是读入所有原子的坐标)，并设定原子的直径和沿某一欲研究的坐标轴上的采样频率；

2)将整个样品沿要研究的坐标轴平移，直到样品边缘在此轴上的坐标为0；

3)只考虑所有原子在这一坐标轴上的投影(即只考虑所有原子在此坐标轴上的坐标)。按照设定的采样频率，逐个采样点统计能覆盖到该采样点的原子个数(可以称为“计数法”)、或者通过该采样点并以该坐标轴为法线的截面(亦即，该坐标轴与该截面垂直)所截到的各个原子面积之和(可以称为“计面积法”)，得到我们所定义的离散“横轴原子密度分布函数”；

4)以离散的横轴原子密度分布函数为原始输入信号，对其进行傅立叶变换(FourierTransform)，得到振幅-频率图，观察振幅-频率图中特定晶向的特征频率(该特征频率是垂直于该晶向的晶面间距的倒数，可以从晶体学中容易的得到)以及其二次或者多次谐波，可以得到初步判定各个沿该坐标轴的特定晶向的存在。

如果继续以下步骤：

5)以设定的截止频率对原始的横轴原子密度分布函数进行高通滤波和低通滤波，以高通滤波后的每一点的信号强度除以低通滤波后每一点的信号强度，得到归一化横轴原子密度分布函数；获得此图的目的是为了排除样品粗细不均匀、几何形状不规则所造成的影响；

6)将样品在空间进行旋转一定角度；可以将样品分别沿三维直角坐标系的x，y，z轴单独或者联合旋转(如先沿x轴旋转特定角度，再沿z轴旋转特定角度，等等)，再重复步骤2)-5)；

7)重复上述步骤6)，直到旋转到一个特定角度，此时在振幅-频率图的某一特征频率附近，其呈现出的振荡峰的振幅比旋转到其他角度都大，频率也比旋转到其他角度都高。这表明旋转到此角度某个晶粒的某个晶向与要研究的坐标轴平行，由此可确定原来某个晶粒的取向，并根据振荡频率与理论推算的该晶向特征频率振荡的差值，可以推断某个晶粒的晶格畸变程度。

8)直接观察旋转各角度后所得到归一化横轴原子密度分布函数图，如果可以直接找到振幅-频率图所反映出的周期性振荡信号的位置，则可以基本确定该取向晶粒的在所研究的坐标轴方向的位置以及长短，根据此处振荡振幅的大小，可以判定出该处结晶程度的高低(振幅越大，结晶越好)，如果难以在“横轴原子密度分布图”中直接确定各个周期性振荡信号的位置，则可能是在与所研究的坐标轴的垂直方向有多种取向的晶粒，它们的信号互相叠加导致难以分辨。

如果需要，可以再将研究的样品切分成数段，逐段研究其“横轴原子密度分布”及其傅立叶变换、“归一化横轴原子密度分布”。

如果需要，比如振幅非常小，或者“振幅-频率图”的杂峰太多，噪音太大，可以尝试改变所设定的原子直径(而不改变原子质心)，再重复以上的步骤。

如果需要，比如振幅非常小，或者“振幅-频率图”的杂峰太多，噪音太大，也可以尝试更改“计数法”为“计面积法”，或者更改“计面积法”为“计数法”，再重复以上的步骤。

本发明与现有技术相比，具有以下显著优点及突出性效果：既可以通过“振幅-频率图”反映大范围内的晶体取向、结晶程度、又可以通过分析“归一化横轴原子密度分布”研究局部的晶体取向和结晶程度。由于将三维的信息化简为一维，并且可以直接从“振幅-频率图”判断晶体的取向，使晶体取向信息得到充分表征，极大的降低了判断晶体取向的难度。由于大部分步骤都可以编程达到自动化，极大的提高了判断晶体取向的效率。表征的结果对选取的角度非常敏感，判定晶体取向的精确度很高。

附图说明

图1是横轴原子密度分布函数定义说明。

图2是[110]‖[111]的铜双晶示意图。

图3a是铜双晶在刚开始弛豫时横轴原子密度分布图。

图3b是铜双晶在刚开始弛豫时振幅-频率图。

图3c是铜双晶在刚开始弛豫时归一横轴原子密度分布图。

图4a是铜双晶在弛豫后拉伸一段时间时横轴原子密度分布图。

图4b是铜双晶在弛豫后拉伸一段时间时振幅-频率图。

图4c是铜双晶在弛豫后拉伸一段时间时归一横轴原子密度分布图。

图5a是铜双晶在拉伸断裂后横轴原子密度分布图。

图5b是铜双晶在拉伸断裂后振幅-频率图。

图5c是铜双晶在拉伸断裂后归一横轴原子密度分布图。

图6a是断裂后铜双晶的左段在未进行旋转时的横轴原子密度分布图。

图6b是断裂后铜双晶的左段在未进行旋转时的振幅-频率图。

图6c是断裂后铜双晶的左段在未进行旋转时的归一横轴原子密度分布图。

图7a是断裂后铜双晶的左段旋转到最佳位置处的横轴原子密度分布图。

图7b是断裂后铜双晶的左段旋转到最佳位置处的振幅-频率图。

图7c是断裂后铜双晶的左段旋转到最佳位置处的归一横轴原子密度分布图。

图8a是断裂后铜双晶的右段在未进行旋转时的横轴原子密度分布图。

图8b是断裂后铜双晶的右段在未进行旋转时的振幅-频率图。

图8c是断裂后铜双晶的右段在未进行旋转时的归一横轴原子密度分布图。

具体实施方式

下面结合具体实施例，进一步阐述以上步骤。

首先阐明本发明的原理：

众所周知，周期性的原子排布是晶体的特征。正是这种周期性的排布造成了晶体的各向异性。也就是，从不同的角度(晶向)去看，一个单晶的物理性质可能并不相同。“横轴原子密度分布”正是一个可以表征晶体各向异性的物理量，它可以很好的区分各个晶向。我们对它的定义如附图1所示：

图1的圆圈代表空间中一个晶粒中的原子(以仅包含同一种元素的单晶为例)。各坐标轴的含义都是三维实空间中的坐标。要注意y轴因为垂直于纸面而没有画出。我们所要研究的坐标轴是z轴，与z轴垂直的虚线正是每个采样点所在的截面。统计每个采样点所在的截面“穿”过的原子的个数(或者面积之和)，就得到z方向的“横轴原子密度分布函数”。可以明显的看出，这个函数必定是周期性振荡的。它的周期就是垂直于该晶向的晶面间距，它的振荡频率就是晶面间距的倒数(晶面间距是指相邻两个等价晶面之间的距离)。只要知道了z轴的晶向，就可以查到z轴的特征振荡频率。例如，对于简单立方晶胞的晶体，其晶面间距d和晶向指数的关系为：

$d_{\left(\mathrm{hkl}\right)}=a/\sqrt{h^2+k^2+l^2}$

其中a是晶格常数。对于面心立方晶体，要考虑到晶面增加的影响，比如对于[100]晶向， $d_{\left(\mathrm{hkl}\right)}=\frac{a/\sqrt{h^2+k^2+l^2}}{2}=\frac{a/\sqrt{1^2+0^2+0^2}}{2}0.5a,$ 其特征振荡频率是2/a。其他晶胞和晶向参见相关文献(例如，《结晶化学导论》(第二版)，钱逸泰编，中国科学技术大学出版社，1999)。

如果将晶体旋转，比如，将晶体旋转绕y轴旋转，使得图中的倾斜虚线于z轴垂直，再如上进行统计，就可以得到垂直于倾斜虚线方向的“横轴原子密度分布函数”。如果用一个3×n三维矩阵表示各原子的坐标，则旋转操作可如下进行：

绕z轴： $\left(\begin{array}{cccc}{x_1}^{\′}& {x_2}^{\′}& ...& {x_n}^{\′}\\ {y_1}^{\′}& {y_2}^{\′}& ...& {y_n}^{\′}\\ {z_1}^{\′}& {z_2}^{\′}& ...& {z_n}^{\′}\end{array}\right)=\left(\begin{array}{ccc}\cos \mathrm{\θ}& -\sin \mathrm{\θ}& 0\\ \sin \mathrm{\θ}& \cos \mathrm{\θ}& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{cccc}{x}_1& x_2& ...& x_n\\ y_1& y_2& ...& y_n\\ z_1& z_2& ...& z_n\end{array}\right)$

绕x轴： $\left(\begin{array}{cccc}{x_1}^{\′}& {x_2}^{\′}& ...& {x_n}^{\′}\\ {y_1}^{\′}& {y_2}^{\′}& ...& {y_n}^{\′}\\ {z_1}^{\′}& {z_2}^{\′}& ...& {z_n}^{\′}\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\θ}& -\sin \mathrm{\θ}\\ 0& \sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right)\left(\begin{array}{cccc}{x}_1& x_2& ...& x_n\\ y_1& y_2& ...& y_n\\ z_1& z_2& ...& z_n\end{array}\right)$

绕y轴： $\left(\begin{array}{cccc}{x_1}^{\′}& {x_2}^{\′}& ...& {x_n}^{\′}\\ {y_1}^{\′}& {y_2}^{\′}& ...& {y_n}^{\′}\\ {z_1}^{\′}& {z_2}^{\′}& ...& {z_n}^{\′}\end{array}\right)=\left(\begin{array}{ccc}\cos \mathrm{\θ}& 0& \sin \mathrm{\θ}\\ 0& 1& 0\\ -\sin \mathrm{\θ}& 0& \cos \mathrm{\θ}\end{array}\right)\left(\begin{array}{cccc}{x}_1& x_2& ...& x_n\\ y_1& y_2& ...& y_n\\ z_1& z_2& ...& z_n\end{array}\right)$

根据上面的定义，如果角度是正的，那么对于绕z、x、y，轴，分别是从x轴正方向往y轴正方向、从y轴正方向往z轴正方向、从z轴正方向往x轴正方向旋转。如果角度是负的，那么对于绕z、x、y，轴，分别是从y轴正方向往x轴正方向、从z轴正方向往y轴正方向、从x轴正方向往z轴正方向旋转。则而根据旋转晶体的角度，很容易反推原来晶向的角度。比如，做z轴的横轴密度分布图时，若发现绕x轴旋转+5°(即把晶体从y轴正方向往z轴正方向旋转了5°)并绕y轴旋转0°可以得到一个[100]晶向，那么这个晶向原来应该是从z轴正方向往y轴正方向偏了5°。

对于所得到的“横轴密度分布函数”，只要进行通常的傅立叶变换，就可以得到“振幅-频率”图。完全舍掉截止波长以下的成分，再进行反傅立叶变换，就得到“高通滤波后的横轴密度分布函数”。完全舍掉截止波长以上的成分，再进行反傅立叶变换，就得到“低通滤波后的横轴密度分布函数”。具体来说，由于经过采样后的“横轴密度分布函数”是一个离散的一维数列{x_n}，所以实际进行的是离散傅立叶变换(Discrete FourierTransformation，DFT)。变换的结果也仍然是一个一维数列。对于有N个采样点的横轴密度分布函数，其DFT可以如下进行：

$F_k=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_n{e}^{-j2k\frac{\mathrm{nk}}{N}},k=\mathrm{0,1,2},...N-1$

由于变换的结果可能是复数，所以做振幅-频率曲线时需要对上式求幅值：

A_k＝|F_k|，k＝1，2，3…N-1

{A_k}(或者{F_k})中的每一点所对应的频率为：

$f_n=\frac{{\mathrm{nF}}_s}{L},n=\mathrm{0,1,2},...N-1$

其中F_s是采样频率(即每nm设几个采样点)，L是样品在所测量轴上的投影长度(单位为nm)。用{A_k}对{f_n}作图，就得到“振幅-频率图”。

对DCT变换的结果{F_k}中对应于高于截止频率f_cutoff的每个点置零，然后再通过如下的逆变换(F_k和x_n带有“’”符号，表示进行过了滤波操作)，就得到低通滤波后的横轴密度分布函数。如果变换后的结果有复数出现，则只取复数的实部。

${x_n}^{\′}=\frac{1}{N}\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{F_k}^{\′}{e}^{j2\mathrm{\π}\frac{\mathrm{nk}}{N}},n=\mathrm{0,1,2},...N-1$

类似的，如果对{F_k}中对应于低于截止频率f_cutoff的每个点置零，再进行逆DCT变换，就得到了高通滤波后的横轴密度分布函数。

下面以[110]‖[111]的铜双晶为实施例进行分析：

如附图2所示，所谓双晶，由就是沿z方向取向不同的两个晶粒所组成一块多晶体。本文所仿真的双晶，左边为[110]晶向，右边为[111]晶向。铜是面心立方晶体，晶格参数为0.3614nm，由晶体学原理，知道[110]晶向上的晶面间距是0.1278nm，对应特征振荡频率7.83nm^-1，[111]晶向上的晶面间距是0.2087nm，对应特征振荡频率4.79nm^-1。此外，还有一个常见的[100]晶向，其晶向上的晶面间距是0.1807nm，对应特征振荡频率是5.53nm^-1。

这个双晶仿真的仿真过程是，弛豫100,000步，再缓慢拉伸直到断裂为止。我们首先用MATLAB程序表征刚开始弛豫的样品，如附图3所示：

图3a是“横轴原子密度分布图”，为低通滤波的结果。原子的直径选择为0.1nm，所用截止频率为3nm^-1。图3b是“振幅-频率”图，图3c是“归一横轴原子密度分布”。在“振幅-频率图”中，我们明显可以看到[110]晶向的特征频率7.83nm^-1及其高次谐波，[111]晶向的特征频率4.79nm^-1及其高次谐波。从“横轴原子密度分布”和“归一横轴原子密度分布”中，我们可以看到振幅很大而且均一，表明整个双晶结晶程度很高、很均匀。

再做出弛豫后拉伸一段时间的图像，此时一共仿真了240,000步。经过检测，仍然是绕x，y轴各旋转0度时各特征峰振幅最大。最终结果如附图4所示，三个子图含义同图3。从“振幅-频率”图中，可以看到[111]的特征振荡向低频移动，表明晶面间距增大，晶格被拉伸。此外，从两个“横轴原子密度分布图”都可以明显看到振幅普遍减小，表明晶体的晶化程度变差，或者晶体晶向有所变化。

最后表征拉伸断裂后的样品，此时仿真500,000步。经过检测，并非绕x，y轴都不旋转所得到的各特征峰强度最大、频率最高。不过还是可以先分析不旋转的表征结果，如附图5所示：

从两个“横轴原子密度图”可以看到，各处的振幅再次增大。进一步对“振幅-频率”图分析，发现导致振幅增大的是约频率为5.3-5.5nm^-1，这基本对应于[100]晶向的特征振荡，频率稍偏低可以解释为晶格被拉伸产生畸变所致。同时，[111]的特征振荡也很明显。表明同时存在[111]晶向和[100]晶向的两种晶粒。

由于双晶已经被拉伸断裂，下一步对双晶分段处理，先观察左端，即0到42nm。未进行旋转的表征如附图6所示，旋转到最佳位置(绕x轴旋转-1°，绕y轴旋转0°)的表征如图7所示。从这两幅图中，都没有发现[111]晶向在4.8附近的特征峰。可见，左侧基本上是一个较为完整的、z轴晶向为[100]晶向的晶粒。其中，8～16nm的部分晶格的晶向稍从z轴正向往y轴负向旋转了约1°。

再对样品的右侧进行分析，表征结果如图8所示。各频率振幅在不旋转时达到最大。而且，没有[100]晶面的特征峰，只有[111]晶面的特征峰及其高次谐波峰。这表明，经过拉伸后，右侧晶体仍然是一个较为完整的、沿z轴为[111]取向的晶粒。

总之，在本实施例中，采用技术方案中所描述的方法进行分析可以得到如下结论：双晶经过拉伸断裂后，左侧的z轴为[110]取向晶粒转变为了z轴为[100]取向的晶粒，右侧则仍为[111]取向的完整晶粒。

Claims (4)

1、一种检测计算机仿真结果中晶粒的取向和结晶程度的方法，其特征在于按以下步骤进行：

1)对于计算机仿真后得到的结果，将其中待检测元素的所有原子的坐标读入计算机中，并设定原子的直径和沿某一欲研究的坐标轴上的采样频率；

2)将整个样品沿要研究的坐标轴平移，直到样品边缘在此轴上的坐标为0；

3)只考虑所有原子在这一坐标轴上的投影，按照设定的采样频率，逐个采样点统计能覆盖到该采样点的原子个数；或者通过该采样点并以该坐标轴为法线的截面所截到的各个原子面积之和，得到离散的横轴原子密度分布函数；

4)以离散的横轴原子密度分布函数为原始输入信号，对其进行傅立叶变换，得到振幅-频率图，观察振幅-频率图中特定晶向的特征频率以及其二次或者多次谐波，可以初步判定沿该坐标轴的各个特定晶向的存在。

2、根据权利要求1所述的检测计算机仿真结果中晶粒的取向和结晶程度的方法，其特征在于还包括以下步骤：

5)以设定的截止频率对原始的横轴原子密度分布函数进行高通滤波和低通滤波，以高通滤波后的每一点的信号强度除以低通滤波后每一点的信号强度，得到归一化的横轴原子密度分布函数，可以排除样品粗细不均匀、几何形状不规则所造成的影响。

3、根据权利要求2所述的检测计算机仿真结果中晶粒的取向和结晶程度的方法，其特征在于还包括以下步骤：

6)将样品在空间旋转一定角度，重复步骤2)-5)；

7)重复上述步骤6)，直到旋转到一个特定角度，此时在振幅-频率图的某一特征频率附近，其呈现出的振荡峰的振幅比旋转到其他角度都大，频率也比旋转到其他角度都高，这表明旋转到此角度某个晶粒的某个晶向与要研究的坐标轴平行，由此可确定原来某个晶粒的取向，并根据振荡频率与理论推算的该晶向特征频率振荡的差值，可以推断某个晶粒的晶格畸变程度；

8)直接观察旋转各角度后所得到的归一化横轴原子密度分布函数图，如果可以直接找到振幅-频率图所反映出的周期性振荡信号的位置，则可以基本确定该取向晶粒的在所研究的坐标轴方向的位置以及长短，根据此处振荡振幅的大小，可以判定出该处结晶程度的高低即振幅越大，结晶越好，如果难以在横轴原子密度分布图中直接确定各个周期性振荡信号的位置，则可能是在与所研究的坐标轴的垂直方向有多种取向的晶粒，它们的信号互相叠加导致难以分辨。

4、根据权利要求1，2或3所述的检测计算机仿真结果中晶粒的取向和结晶程度的方法，其特征在于步骤4)中的特征频率是垂直于该晶向的晶面间距的倒数。

Images