LAMMPS 晶界能与取向分析

本篇继续镍双晶建模，把实际松弛末态转成可检查的能量、取向与空间分布。下载完整输入、原始数据与分析包。

1. 先生成表格#

自己的 GPU 作业完成后：

bash analyze.sh work

如果只想先学习配套真实数据，在新解压目录直接运行：

bash analyze.sh

默认读取 reference/，不重新提交 LAMMPS 松弛。两种方式均会生成能量 CSV、取向 CSV 和四张 PNG，原生坐标、逐原子能量与取向值保留。

head -n 9 tables/bicrystal-energies.csv
head -n 6 tables/orientation-summary.csv

2. 为什么不能直接比较两个总能量#

完美晶体有 1280 原子，双晶删除过近原子后为 1248 原子，原子数不同。直接比较总能量会混入少了 32 个原子的影响。

本例先用相同晶胞的完美 A 晶体得到每原子参考能 e_bulk≈−4.450000003577 eV/atom；完美 B 的结果与其一致。再按：

excess_energy = E_bicrystal - N_bicrystal * e_bulk
gamma = excess_energy / (2 * Lx * Lz)

分母的 2 来自两个周期界面，单个面积为 313.45508096 Ų。eV/Ų 乘 16.02176634 转换成 J/m²。

P100、初始平移 0 的实际例子：

另一个初始平移得到约 1.341657082 J/m²，在本例精度下几乎相同。这里只比较两种起始平移，不能把这个数称为经过全局结构搜索与尺寸收敛的 Ni 晶界能；也没有振动、磁性或有限温度熵，因此不是晶界自由能。

右图的高能区域集中在中间界面和周期接缝附近。它是分箱后的每原子势能均值，不是每层局部晶界能；不能把曲线高度直接替换进面积公式。

3. OVITO 怎样检查晶粒取向#

脚本调用现有 OVITO Python 3.16.1 / PolyhedralTemplateMatchingModifier，输出 Structure Type 与 Orientation。GUI 中可打开：

work/p100-GB-shift0.0-prune1.4/final.dump

添加 Polyhedral Template Matching，开启取向输出，保持 xyz 周期边界。观察 FCC 与 Other 的空间分布，再比较晶粒内部的取向。Other 可以来自晶界局部畸变，不自动表示一个新体相。

本例的取向检查步骤是：

1. 按原子 ID 与初始晶粒标签对应。
2. 选择距离两个 y 界面至少 8 Å 的晶粒内部原子。
3. 对 FCC 的原生四元数先转换为双精度、明确单位化，再转成旋转矩阵；原始单精度值另存。
4. 考虑立方晶体的 24 个正旋转对称操作，再求最小错取向，不直接平均四元数。

两个晶粒内部各取到 176 个 FCC 原子。松弛后相对各自输入晶向的中位偏差约 0.09°；理想两晶粒的最小立方错取向为 36.86989765°。结构标签从初始 992 FCC/256 Other 变为末态 1056 FCC/192 Other，总原子数保持 1248。

左图把晶粒内部取向与输入 A 比较，因此 A、B 分别集中在接近 0° 和 36.87°；它与“相对各自输入晶向的约 0.09° 偏差”是两个不同的角度定义。

4. 能量差分怎样核对力#

同一份 P100 松弛末态交给 V100 16G 做静态单点，比较全部原子的能量与力，避免把两个独立松弛路径的差异直接当作 GPU 数值误差。

随后，对预先选定的三个原子、每个 x/y/z 分量施加 ±10⁻⁵ Å 微小位移：

F_component = -[E(x+h)-E(x-h)] / (2*h)

四个独立模型共检查 36 个力分量，本次能量导数与原生力最大差约 9.1×10⁻⁸ eV/Å。这说明本例在相同几何上的能量/力提取对应，不替代晶界宽度、盒长或势函数的研究验证。

5. 常用输出#

bash plot.sh reference

这会用提供的原始实算数据重画参考图；要保留自己新计算的图片，请先复制 figures。正式研究时，继续扫描晶粒平移、边界原子密度、晶胞厚度及边界条件，再考虑温度与所用势的适用范围。

参考：Atomsk 双晶建模、OVITO PTM。