计算化学公社

标题: 一个辅助构建理论六方异质结晶格匹配超胞的Excel表格 [打印本页]
作者
Author: Uus/pMeC6H4-/キ    时间: 2025-12-19 22:24
标题: 一个辅助构建理论六方异质结晶格匹配超胞的Excel表格
本帖最后由 Uus/pMeC6H4-/キ 于 2025-12-21 15:22 编辑

一个辅助构建理论六方异质结晶格匹配超胞的Excel表格
An Excel spreadsheet for aid in constructing theoretical commensurate superlattice of hexagonal heterojunction
Uus_pMeC6H4- (\ Init post 2025-12-19 /) Last edit 2025-12-21


构建材料异质结是很多第一性原理计算研究的第一步，背后的数学原理涉及周期性平移的晶格矢量的线性变换。平面六方点阵的点阵点间距是可以枚举的，若再找到相应的变换矩阵，即可将有六方对称性的晶格扩胞使得晶胞参数大小合适，有助于构建六方异质结的晶格匹配超胞。（此处说的六方对称性的晶格相当于考虑二维周期性，包括六方、三方晶系的(001)面与面心立方的(111)面，如果连交界处的晶向都没确定就不是此处关心的内容了。）

本帖将分享一个处理该问题的Excel表格，虽然自动化程度有限，但没有复杂依赖且展示比较直观。这里的建模单纯考虑理论上的数字匹配，没有完整考虑如实际晶体材料中不同晶面的不同暴露终端与封端、表面重构与表面弛豫导致明显结构变化（比如经典的h-BN与Ru(111)晶格不匹配产生波纹状纳米网格）等等。
(, 下载次数 Times of downloads: 13)

该工作簿有两张工作表，其中WorkPlace作为填写参数与计算结果的区域，而HexGrid在距离原点64个单位的范围内枚举了格点及对应的模（距离平方，正比于向量所张成晶面的面积），并按模给格点着渐变色、标记六方点阵的坐标轴$\vec{a}$与$\vec{b}$的方向。由于范围比较大，HexGrid缩放40%来看，大致效果为
(, 下载次数 Times of downloads: 0)

在WorkPlace里指定的粉白色底格子内填写超胞晶格矢量上限、失配度上限及两种六方晶格的晶胞参数（要求二者之比不超过64），待运算完成给出最小超胞尺寸或最小失配度对应的模后，到HexGrid查找模对应的格子，即可写出变换矩阵。鼠标悬浮在格子右上角可以查看注释。文件中给出的数据示例分别对应COD-9008569的六方石墨（2.456埃）与COD-2016170的六方氮化硼（2.498埃），在失配度不超过0.15%且超胞尺寸不大于20埃的要求下可以找到两个解：将氮化硼做$\sqrt{27}$扩胞且将石墨做$\sqrt{28}$扩胞（指边长，下同）是超胞尺寸最小的解，将氮化硼做$\sqrt{61}$扩胞且将石墨做$\sqrt{63}$扩胞是超胞尺寸稍大但失配度更小的解。
(, 下载次数 Times of downloads: 1)

在HexGrid中放大视图使格子内数字可见，用ctrl+f查找61可看到三个结果。
(, 下载次数 Times of downloads: 1)

取最靠近右下方坐标轴的格子，可见其对应矢量$9 \vec{a} + 4 \vec{b}$。旁边的63对应矢量$9 \vec{a} + 3 \vec{b}$。则按行写出2x2变换矩阵分别为
  1. $\begin{bmatrix}
  2. 9 & 4 \\
  3. -4 & 5
  4. \end{bmatrix}$
复制代码
  1. $\begin{bmatrix}
  2. 9 & 3 \\
  3. -3 & 6
  4. \end{bmatrix}$
复制代码
（一般地，对于$p \vec{a} + q \vec{b}$，变换矩阵第一行即为p与q，第二行则为-q与p-q。）

有了变换矩阵，在如VESTA等可视化软件内建模也不难，之前我也在论坛分享过几个晶格矢量线性变换的例子，在这个回复帖的第二段有所总结。必须注意VESTA的晶格矢量线性变换矩阵的输入要求，得按列向量而不是行向量来书写。假如要建一层氮化硼与一层石墨烯的异质结，则操作流程为：
1. 打开六方氮化硼的结构2016170.cif，在Edit - Edit data窗口Unit cell栏点Remove symmtery去除对称性，切换到Structure parameters栏删除z=0.5处的原子而只保留一层，再切换回Unit cell栏点Transform按钮输入变换矩阵
  1. 9 -4 0
  2. 4 5 0
  3. 0 0 1
复制代码
最后File - Export导出所得超胞结构。
2. 打开六方石墨的结构9008569.cif，同样去除对称性、删除z=0处的原子而只保留一层，输入变换矩阵
  1. 9 -3 0
  2. 3 6 0
  3. 0 0 1
复制代码
同样导出超胞结构。
3. 最后用文本编辑器编辑结构文件，复制粘贴原子信息来合并结构。最终成果如下所示，盒子在a和b方向的尺寸约19.5埃、化学式为C126B61N61，异质结的晶格匹配不错。
(, 下载次数 Times of downloads: 1)


下面再举一个假设性的三层CeO2 (111)与双层g-C3N4的异质结作为例子，假定二者接触界面垂直于z方向。

面心立方的萤石型CeO2结构来自COD-9009008，石墨烯型氮化碳的结构来自Angew. Chem. Int. Ed., 2014, 53(29), 7450-7455的补充材料。前者对应的以埃为单位的六方晶格矢量长度应为立方晶胞边长的$\frac{1}{\sqrt{2}}$，即5.4110/1.4142 = 3.8262，后者则为5.0415。把两个数字填入上述Excel表格，发现放宽失配度上限为0.40%后可以找到一组较小的解，即把CeO2做$\sqrt{7}$扩胞、g-C3N4做$\sqrt{4}$扩胞。不过两个cif文件不能直接拿来合并，需要一点预处理。

在VESTA中打开CeO2的结构9009008.cif，这里先不去除面心立方的对称性，直接到Edit - Edit data - Unit cell输入变换矩阵
  1. 0.5 -0.5 2
  2. 0 0.5 2
  3. -0.5 0 2
复制代码
其中第三列是垂直于目标(111)晶面的[111]晶向，扩大为2倍方便选取层，而前两列是(111)晶面上相应于六方对称性的单位矢量，面心与顶点的等价性允许使用0.5而不是限定于整数。在Additional lattice point(s) are included in the new unit-cell窗口内选择Add new equivalent positions to a list of symmetry operations选项来确定添加原子的算法。点Remove symmetry去除对称性，到Structure parameters，单击z列的表头使原子按z方向分数坐标顺序排列，然后选中z方向分数坐标小于0.25或大于0.75的行点右边Delete删除，剩下3个Ce和6个O，符合化学计量比且是所需的三层结构。打开Utilities - Equivalent Positions窗口（忽略启动时的警告），在General equivalent positions里按键盘Del键删除多余的等价位置而只留下x y z一个，最后File - Export Data导出cif结构为CeO2-layers.cif。

在VESTA中打开g-C3N4的结构C3N4.cif（如果嫌晶胞外平移复制延展的原子太多难以观察，可以在Edit - Bonds选中唯一的成键项，把Boundary mode改成Do not search atoms beyond the boundary），用变换矩阵
  1. 1 0 0
  2. 0 1 0
  3. 0 0 2
复制代码
在c方向扩胞2倍而a、b方向不变，移除对称性、删除原子坐标和等价位置使得只有分数坐标小于等于0.25的保留，导出cif结构为C3N4-layers.cif。

预处理结束，接下来再构建晶格匹配的超胞。打开CeO2-layers.cif，施变换矩阵
  1. 3 -1 0
  2. 1 2 0
  3. 0 0 1
复制代码
获得a = b = 10.123，c = 18.744的超胞，导出记录笛卡尔坐标的POSCAR；打开C3N4-layers.cif，施变换矩阵
  1. 2 0 0
  2. 0 2 0
  3. 0 0 1
复制代码
获得a = b = 10.083, c = 13.153的超胞，导出记录笛卡尔坐标的POSCAR。

在用文本编辑器合并两个POSCAR时，应该采用什么晶格矢量？考虑到异质结建模本身只能提供初始结构、适应外压的变胞优化仍为必需，晶格矢量只要选得原子间距合理，不让周期镜像之间有异常重叠或缝隙即可。具体地说：
* 对a和b方向，两个超胞已经是按失配度较小的条件找到的，用二者之一或者二者的几何平均差别不大，不过我习惯用二者之中更大的那个；
* 对c方向，如果是导出cif文件或者记录分数坐标的POSCAR，当两个结构的c方向晶胞参数差别明显时相同的分数坐标差实际对应着不同的原子层间距，可能有拉长或压扁的问题，然而此处用了记录笛卡尔坐标的POSCAR而没有该问题，直接选一个足够大的即可，顺带着也可以扩充真空区。

由于构建CeO2-layers.cif与C3N4-layers.cif时保留的层位置适当，两种结构在c方向不太近也不太远。最终得到的结构有119个原子，化学式为C24N32Ce21O42，两个视角如下图所示。
(, 下载次数 Times of downloads: 0)


目前先写这么多，欢迎举出更多应用实例（若有讨论要求，必须提供来源清晰可追溯的cif文件下载），后面可以找时间再添加演示。
作者
Author: 玉米猫    时间: 2025-12-20 23:20
好像有个笔误 “旁边的28对应矢量”应该是“旁边的63对应矢量”对不对
作者
Author: Uus/pMeC6H4-/キ    时间: 2025-12-21 15:25
今日添加了一个g-C3N4与CeO2 (111)异质结的建模例子。

顺便再放几点相关讨论：
1. 表格里Results旁边有一块区域，原本是希望直接输出变换矩阵的，但在我这里返回HexGrids查找相应单元格位置的Excel函数有点故障，只能留空了。实际上对HexGrid用ctrl+f查找也不完全可靠，哪怕Excel里公式 - 计算 - 计算选项设置为自动来生成整个HexGrid，查找时也有可能返回0个结果。
2. 表格里用的失配率是按两种超胞的晶胞参数之差占晶胞参数的比例来计算的，设置一个如1%的失配率上限对于小晶胞还好，而对于比较大的晶胞则必须留意实际空间距离误差超过0.5埃（随便举的例子，不是严格的阈值）导致原子有异常接触的问题，建议降低失配率上限、增大超胞尺寸上限来找合适的解。
3. 所谓“根号扩胞”其实是英文名 Wood's notation 的记号的不太严谨的简化，而晶格矢量变换矩阵对应的记号就叫Matrix notation，这在讨论表面吸附的维基Overlayer页有简单介绍，用Google搜 Wood's notation 也能找到一些课件。
4. 对twistronics研究的双层二维转角体系如平行错位pi-pi堆积的转角石墨烯，建孤立团簇模型用量子化学软件做计算是比较常见的做法，我估计这可能与建周期性模型不好选取适应实际夹角对应的点阵矢量有关，比如$64 \vec{a} + \vec{b}$与$64 \vec{a} - \vec{b}$的夹角也有1.5627度，$63 \vec{a} + 31\vec{b}$与$63 \vec{a} + 32\vec{b}$的夹角为1.0501度，然而此时超胞边长已经是原胞的五六十倍了，原子数更是接近三四千倍，计算量极大。



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3