【入门基础教程】晶体与表面计算的一些基础知识

本教程将介绍：

• 三维晶体结构
  1. 晶体建模
  2. k点的选取与收敛性
  3. 晶格常数优化
• 二维体系
  1. 两种表面模型：Slab与QMMM
  2. 密勒指数
  3. 表面建模
  4. 表面钝化，极性、非极性表面
  5. 表面弛豫、重建
  6. 表面超胞及其矩阵表示
  7. 表面能
  8. 吸附位点和吸附能计算

本教程仅涵盖三维周期性晶体及其表面，不讨论准晶（非周期性晶体）。AMSinput可以导入（File → Import coordinates）、导出（File → Export coordinates）多种不同的晶体结构文件格式，包括.cif和.pdb等。导出*.xyz格式时，AMS还会将晶格矢量写入*.xyz文件中。

在AMS中，3D周期性晶体可以表示为某个单胞的周期性重复，这叫做周期边界条件（PBC）。单胞由三个晶格矢量定义，m每种晶体可以定义出无限种单胞供选择，其中最重要的两个是：

• 晶体学原胞（primitive cell），它是尽可能小的单胞（仅包含一个格点）
• 惯用单胞（conventional cell），遵循国际晶体学联合会概述的指导方针，

例如，铜是一种具有面心立方（fcc）密堆积（ccp）金属。晶体学原胞和惯用单胞是：

在AMSinput中，该晶体学原胞：

惯用单胞：

注意，伸向Cell边界的“棍”，表示该原子与相邻Cell内的原子存在“键”，但因为是周期性的，因此只显示一个Cell即可。不少晶体的惯用单胞和晶体学原胞是一样的。本例中，晶体学原胞只有1个原子，惯用单胞有4个原子，因此计算效率上晶体学原胞会快很多。

一个晶体，在AMSinput中，点击窗口底部的❄按钮 → Convert to Primate Cell即可得到晶体的晶体学原胞，❄按钮 → Convert to Conventional Cell可以得到惯用单胞。底部最右侧的四个圆圈按钮，可以切换显示单/多个周期，便于用户检查周期性是否正确，关于周期性如果理解不充分，强烈建议仔细阅读：周期边界条件

AMSinput中，View → Periodic → Show Lattice Vectors可以显示晶格矢量的方向，红绿蓝三条线分别表示\(\vec{a}\)、\(\vec{b}\)、\(\vec{c}\)的方向，View → Periodic → Show Lattice Vector Labels可以显示晶格矢量的这三个字符。View → Periodic → Show Unit Cell 可以显示惯用单胞/晶体学原胞的形状。View → Periodic → Unit Cell Style可以选择晶体学原胞/惯用单胞的显示风格，是否显示6个面，是否只显示12条边，边的粗细。

• 进行电子或声子带结构计算，建议尽量使用使用晶体学原胞（对于声子带结构的计算，应该使用晶体学原胞，但也需要指定声子计算内部使用的超胞的大小，这是为了考虑相邻晶胞原子的影响）。
• 根据密勒指数创建Slab时，应使用惯用单胞。
• 晶体存在掺杂、缺陷时，例如想引入空位缺陷时，应使用一些合适的超胞，超晶胞的大小决定了缺陷的浓度，超胞更大，可以得到的缺陷浓度就可以更低，但计算量也随之增大。如果原子偏离晶体位置，超胞也可以用来近似模拟非晶固体。

晶体学中，常用6个参数来描述惯用单胞/晶体学原胞：a、b、c、α、β、γ。其中

• a、b、c是\(\vec{a}\)、\(\vec{b}\)、\(\vec{c}\)的长度
• α是\(\vec{b}\)与\(\vec{c}\)的夹角
• β是\(\vec{c}\)与\(\vec{a}\)的夹角
• γ是\(\vec{a}\)与\(\vec{b}\)的夹角

在AMSinput中，Model → Lattice可以看到这些参数。调整晶格常数的时候，如果希望原子的分数坐标不变，整个晶体等比压缩、膨胀，可以勾选Adjust atoms when changing lattice vectors。

原子的位置可以用笛卡尔坐标，单位为Å（AMS中的默认值）或分数坐标（以晶格矢量为基础，在晶体学中很常见）来描述。AMSinput → Coordinates勾选Use fractional coordinates即可将笛卡尔坐标显示为分数坐标。

晶格矢量的三个分量是以Å为单位的笛卡尔（xyz）坐标。晶格矢量可以相对于笛卡尔轴以无限多的方式旋转。在AMS中，晶格矢量旋转到哪个方向通常无关紧要。不过一般有两种约定俗成二选一：

• \(\vec{c}\)平行于z轴，则\(\vec{b}\)一般处于yz平面
• \(\vec{a}\)平行于x轴，则\(\vec{b}\)一般处于xy平面

用户可以在AMSinpu中，按照上述两种规则，轻松调整晶格矢量：AMSinput → Model勾选Adjust atoms when changing lattice vectors，然后点击Align c to z-axis或Align a to x-axis按钮即可。

通过AMS使用DFTB、BAND、Quantum ESPRESSO和VASP，必须为周期性计算指定一些k点采样。ReaxFF、ForceField和ML Potential不使用k点。在DFTB和BAND中，可以使用值GammaOnly、Basic、Normal、Good、VeryGood和Excellent指定k空间质量，也可以手动指定沿每个倒空间矢量的k点数量，或者对于高度对称的体系，使用对称k空间格点。

在Quantum ESPRESSO和VASP via AMS中，必须手动指定Monkhorst Pack格点k点数量，所需的k点数量取决于晶胞材料和以及晶胞的大小：

• 金属需要的k点很多
• 绝缘体、半导体和分子晶体通常需的k点较少
• 小的晶胞需要更多k点
• 大型晶胞和超晶胞需要更少k点

为选取的方法和体系确定合适的k点设置，需要先进行k点收敛测试，逐渐提高k点数目，直到您感兴趣的量（例如能量或晶格参数）逐渐趋于一个收敛的结果（即，进一步提高k点数目，能量等不再显著变化）。参见示例：下文《纤锌矿ZnO与DFTB的晶格优化》一节。

一般不要去比较用不同的k点密度计算得到的能量。例如，如果从晶体学原胞创建超胞并引入间隙缺陷，则计算缺陷插入能，就需要与非缺陷超胞的能量比对，而不是直接把晶体学原胞能量的倍数拿来比较。具体参考：下文《表面能计算》一节。

在计算化学中，在构建超胞或表面之前，先优化晶体的晶格矢量是很常见的操作。在计算声子或弹性张量之前，必须首先优化晶格。在AMSinput中，Main窗口“Task”→ Geometry optimization 即进行几何优化，然后在Details → Geometry optimization面板勾选Optimize lattice。还可以设置最大迭代次数maximum number of iterations，即结构优化次数上限。

四个收敛标准：

• 最大梯度（原子力）分量
• 梯度（原子力）分量的均方根
• 从上一步到这一步原子的最大位移
• 每个原子的最大应力能

前三个定义为用于普通的几何优化。第四项被用作应力张量分量的收敛标准，它是“应力张量*晶胞体积/原子数”的张量矩阵元中的最大值。在AMSinput → Details → Geometry Optimization面板convergence criteria下面的值，可以修改收敛标准。

• 如果使用BAND、DFTB、Quantum ESPRESSO或VASP，请使用密集（Good）k点网格。k点网格的质量取决于晶格，在晶格优化过程中，晶格会发生变化，因此需确保质量高于正常水平。
• 对于晶格优化来说，BAND可能相当慢。考虑使用Quantum ESPRESSO优化晶格，然后切换回BAND进行其他计算。
• 若需要各向同性缩放优化晶格常数，Model → Geometry Constraints and PES Scan，并勾选Freeze Strain下面的非对角项（xy、xz和yz）的勾选框和Equal Strain下的对角项（xx、yy和zz）勾选框。
• 若要在外部压力下优化晶格，Details → Geometry Optimization，设置Pressure的值。

也可以使用ForceField（UFF）、BAND或Quantum ESPRESSO运行此示例。ZnO结晶为六方纤锌矿结构，其中a=b，α=β=90°，γ=120°。

1. AMSinput切换到DFTBP模块
2. File → Import Coordinates导入ZnO_wurtzite.xyz
3. 设置Model为SCC-DFTB
4. Parameter directory选择DFTB.org/znorg-0-1
5. K-space设为Good.
6. Task设为Geometry Optimization，并点击右侧>按钮进入结构优化详细设置窗口
7. 勾选Optimize lattice: Yes
8. File → Run保存作业名为zno_lattopt.ams并运行之
9. 作业运行完成后，AMSinput询问whether to update the coordinates，选择Yes, New job，从而根据优化好的结构重新创建了一个作业
10. Model → Lattice查看新的晶格常数
11. File → Export coordinates → .xyz导出当前窗口这个已经优化好的结构坐标文件，存储名为：ZnO_optimized_znorg.xyz

下表显示不同k点个数，优化晶格常数的情况：

基本上Good和Very Good之间已经没有什么变化，Good和Normal之间则尚有较大变化，因此选择Good K点是很好的。

• Band/DFTB/QE计算带隙：DFTB和BAND计算的*.logfile文件的末尾会显示均匀k点的带隙，从能带图中也可以得到沿高对称点的路径上的带隙，当k点足够密集时，二者一般会趋向一致。
• BAND/DFTB/QE计算电子态密度（DOS）：在Main面板上，勾选Calculate DOS复选框。这也计算了partial DOS、PDOS。
• BAND/DFTB/QE的能带结构：在Main面板上，勾选Calculate band structure复选框。能带结构通常是以晶体学原胞计算的。Properties → Band Structure可以通过布里渊区自定义路径。
• BAND计算费米表面：在AMSInput → Properties → Band Structure勾选Calculate Fermi surface，计算出的费米表面可以在SCM → Band Structure → View → Brillouin Zone看到。
• BAND/DFTB声子色散曲线和声子DOS计算：Properties → Phonons and Elastic tensor → Phonons勾选即可，如果计算声子，则需要在下方矩阵中填入超胞规模，例如2×2等，数字越大声子色散曲线越趋于一个收敛的形状。要注意需要首先优化晶格，另请参阅硅晶格优化和声子。
• BAND/DFTB体积模量、剪切模量、杨氏模量、弹性张量：Properties → Phonons and Elastic tensor → Elastic tensor勾选即可。注意应首先用更高的收敛准则对晶格进行优化。
• 应力-应变曲线，屈服点：请参见教程捕捉聚乙炔链。
• 热膨胀系数：在几个不同的温度下运行平衡良好的NPT分子动力学模拟：请参阅教程热固性聚合物的热膨胀系数。
• 原子化能和晶格能：参见教程键能计算。

在原子建模中，主要有两种表示表面的方法：

• Slab模型：在横向方向上是周期性的，也可以在垂直表面方向上也可以是周期性的，如此则需要用足够大的真空间隙将周期性的Slab隔开，从而得到一个周期重复的Slab模型）。AMS的BAND、DFTB和ReaxFF可以以二维周期性建模，并且不需要真空间隙，即Slab外为半无限大的真空。Quantum ESPRESSO和VASP则需要真空间隙，实际上使用的是三维周期性。
• embedded cluster (QM/MM)模型：表面使用高理论水平（例如DFT）模拟，其余部分用低理论水平（例如力场或点电荷）模拟。ADF、BAND和DFTB可用于这种类型的表面模型。通过AMS中的Hybrid模块，可以完成混合引擎的QM/MM模拟。

本教程的其余部分仅介绍Slab模型，这是最常见的表面模型类型。有关QM/MM的更多详细信息，请参阅教程有机框架材料中的无机链接器。

AMS中各个计算引擎对二维Slab、三维Slab、QMMM的支持情况：

表面通常是基于密勒指数定义的，表示为(hkl)。密勒指数是该平面与惯用单胞（Conventional Cell）晶格矢量的“分数截距”的倒数。如果指数为0，则表面与该方向平行。负指数用顶上的条形表示。

• 示例1：Cu(100)表面平行于Cu惯用单胞中的\(\vec{b}\)和\(\vec{c}\)矢量。
• 示例2：Cu(111)表面是对Cu惯用单胞的\(\vec{a}\)、\(\vec{b}\)、\(\vec{c}\)三个晶格矢量的顶点，三个点形成的平面。
• 示例3：Cu(511)表面是对Cu惯用单胞的\(\vec{b}\)和\(\vec{c}\)矢量的顶点，以及\(\vec{a}\)矢量五分之一处，共三个点形成的平面。

分别图示如下：

对于六方晶体（例如Zn或ZnO），经常使用四个米勒指数(hkil)。则第三指数i是冗余的：h+k+i＝0。对于具有菱形晶格的六方晶体（如α-Al2O3或方解石CaCO3），密勒指数是基于菱形晶体学原胞或六方单胞，你需要知道使用了哪种约定。在发布或与他人交流时，应明确说明使用的约定。

晶向指数通过[hkl]的格式描述，代表方向h\(\vec{a}\)+k\(\vec{b}\)+l\(\vec{c}\)。

该方向不一定垂直于平面(hkl)。例如：在α=γ=90°和β>90°的单斜晶胞中，[001]方向平行于\(\vec{c}\)，但不垂直于(001)平面。(001)平面平行于\(\vec{a}\)与\(\vec{b}\)。

对于六方晶格，晶向可以使用如上所述的三个指数来给出，也可以使用特殊的四指数[hkil]表示。材料科学或晶体学教科书中描述了如何在这两种类型的符号之间进行转换。注意，晶向指数一般都是基于惯用单胞（Conventional Cell）的。

AMS的图形窗口不显示晶向。

1. 创建Cu单晶：点击并复制坐标文件的内容，在AMSinput使用Ctrl V粘贴即完成。
2. 点击❄按钮 → Convert To Conventional Cell
3. 点击❄按钮 → Generate Slab
4. 输入Miller indices：5、1、1
5. Number of layers: 10
6. 点击Generate slab按钮即生成Slab Cu(511)，如上图中所示，Slab的法线方向平行于z轴。如果使用BAND、DFTB、ForceField、ReaxFF生成的是二维Slab（如果Periodicity改为Bulk则为三维Slab），QE、VASP则生成三维Slab，因为QE、VASP不支持二维Slab。生成的三维Slab有一定厚度的真空层，Model → Lattice可以看到晶格矢量\(\vec{c}\)的大小。

如果表面晶胞是六边形的（例如：Cu(111)和ZnO(0001)），则角度\(\vec{a}\)与\(\vec{b}\)的夹角γ为60°，不过更常用120°，可以在Model → Lattice修改，或者❄ → Convert To Conventional Cell。

密勒指数(hkl)虽然是基于惯用单胞（Conventional Cell）晶格矢量的截距定义的，但实际上任何与该平面平行的任何平面密勒指数都是(hkl)。生成Slab的时候，点击Generate Slab前，哪个原子是被选中的，Slab表面就从哪个原子那里截断，那个原子会保留到Slab的表面。生成的表面应该仔细检查，未必是稳定、合理的表面。

通过复制粘贴坐标，在AMSinpu中创建ZnO，并生成3层ZnO(0001)（对于六边形晶体，忽略第三个指数，因此键入0，0，1即可）。首先尝试在选择Zn(1)原子的时候生成Slab，然后尝试选择O(3)原子然后生成Slab，得到不同终断的表面。在第一种情况下，表面原子是3-配位，这是有利的。第二种情况下，表面原子为1-配位，就不那么稳定了。

类似地，您可以生成ZnO(\(10\bar{1}0\))的各种好、不好的表面。

下图所示的多种表面俯视图，其中a、c是良好的表面，b、d是不合理表面：

一般会创建比所需厚度更厚的Slab，然后手动删除原子，直到表面终断正确为止。

添加与表面原子配位的H原子，可以实现表面的钝化，这样的H原子可以钝化悬空键，悬空键是非常不稳定、活跃的。例如：金刚石(100)表面终端是2配位的sp$^3$ C原子，这是非常不稳定的。通过用H饱和，表面的悬空键将被钝化。软件操作：按住Shift键选中悬空的C原子，Ctrl e自动加氢。

未钝化的Slab的态密度，可以看到费米能级处，态密度非常大：

钝化的Slab的态密度，费米能级处，态密度正常（这不是导体，因此费米能级处DOS为0，费米能级处于带隙中）：

不过Ctrl e自动加氢，对主族元素一般没有问题，对富足元素往往需要手动添加H原子。加氢之后，H原子的位置较为随机，需要优化一下。

Slab的几何优化，表面截断位置的原子通常从原先偏移，例如可能移动到更靠近次表面原子的位置，次表层中的原子也会有一些变化。 要量化表面弛豫程度，三种方式如下：

1. SCM → Movie → 选中关心的原子 → View → Coordinates，然后查看第一帧、最后一帧，可以看到坐标的变化。
2. SCM → KFbrowser → File → Expert Mode → InputMolecule → Coords输入的坐标，SCM → KFbrowser → File → Expert Mode → Molecule → Coords可以找到最后一组坐标，即最终优化坐标。注意坐标单位是Bohr。1 Bohr=0.5291772083 Å
3. AMSJobs → Tools → Build spreadsheet可以看到坐标比对

表面重建是指与体材料相比，表面发生显著的结构变化。例如：Si(100)表面形成（2×1）重建，其中表面原子形成二聚行（在体相中不存在）。（2×2）和c（4×2）重建也是可能的。

未重构的Si(100)Slab，使用BAND/HSE06计算DOS，可以看到费米能级处，DOS非常高，体现出金属特征（这显然是不对的），因此计算的时候需要设置较高的电子温度（如前文所述），否则SCF不收敛：

重构之后的Si(100)Slab，存在带隙，PBE优化时SCF即可收敛，使用BAND/HSE06计算DOS，费米能级处于价带顶，体现出半导体特征（正确、符合实际）

极性Zn作为终端的ZnO(0001)表面，可以形成相当大的三角形凹坑和岛，化学计量发生了变化，这有助于稳定表面（原因详见下文“极性表面的稳定机制”）。

非极性表面是最常见的，通常比极性表面更稳定。不过对于离子化合物，表面可能就是极性的：如果垂直于Slab的Formal Dipole Moment随着平板厚度的增加而增加，则Slab是极性的。Formal Dipole Moment是根据离子的Formal Charges计算的，随着Slab厚度的增加，偶极矩不断增长，会导致能量发散。

不过极性表面是可以在现实中形成的，因为它们能以某种方式稳定下来。在计算中，极性Slab暴露出两个不同的极性表面，一个在“顶侧”，一个位于“底侧”。例如，前文中的ZnO(0001)面Slab顶部的Zn端和(\(000\bar{1}\))面底部的O端。

极性表面可以通过如下因素稳定下来：

• 化学计量的变化，例如空位或其他缺陷，如上面的ZnO
• 吸附
• 电荷从Slab的一侧转移到另一侧，可能导致表面金属化。这在电子结构计算中很常见，但在实验中并不常见。
• 对于薄Slab，整个结构都可能会发生变化

例如ZnO(0001)能够通过如下方式稳定：

• 1/4表面Zn空位
• 1/2 ML OH吸附

ZnO(\(000\bar{1}\))能够通过如下方式稳定：

• 1/4表面O空位
• 1/2 ML H吸附

案例：

1. 下载polar_stabilization.zip并解压，得到4个结构文件
2. AMSinput → File → Import Coordinates导入clean_good.xyz
3. 切换到DFTB模块，Task：Single Point，Model：SCC-DFTB
4. KSpace：Good
5. Parameter directory：DFTB.org/znorg-0-1
6. 保存并提交作业
7. SCM → Dos查看DOS
8. 选择DOS窗口左侧的1个/多个原子，右侧可以看到对应的partial DOS
9. 分别看看顶层（Zn）、中间、底层（O）的partial DOS
10. 对vacancies_good.xyz和adsorbates_good.xyz重复上述计算

下图a显示了DFTB计算的清洁ZnO(0001)Slab的DOS。从顶层到底层时，能带的位置会发生变化。Zn端侧的导带与O端侧上的价带重叠，导致电荷转移和表面金属化。如前文所述，此时SCC-DFTB（或DFT计算的SCF）可能难以收敛，需要使用电子温度。

b和c显示空位或吸附稳定后，DFTB计算的DOS。在这些情况下，没有表面金属化现象，ZnO Slab仍然是半导体。

理想（未稳定）的极性表面在DFT和DFTB计算中可能都难以收敛，可以尝试通过化学计量的变化或吸附来稳定表面。

就像晶体一样，表面也有一个晶体学原胞（只包含一个晶格点），不过在表面相关计算中，经常会使用表面超胞。通过创建超胞，可以对吸附质的覆盖范围进行微调。

AMSinput中生成的Slab一般是惯用单胞（Conventional Cell），可以通过❄ → Convert to Primitive Cell得到晶体学原胞，二者一般仅仅对中心正交晶格有区别。例：

1. Cu单晶（点击进入复制内容，到AMSinput粘贴完成建模）
2. ❄ → Convert to Conventional Cell
3. ❄ → Generate Slab，密勒指数设为1，1，1，厚度5，然后Generate Slab按钮。

关于Model → Lattice页面，勾选Adjust atoms when changing lattice，然后旋转曲面晶格矢量，然后窗口中部点击\(\vec{a}\)→x axis，将\(\vec{a}\)旋转到平行于x轴。

超胞通常使用Wood表示（AMSinput不支持）或矩阵表示（AMSinput支持）进行描述。

其中a) Cu(100)和b) Cu(111)表面的顶视图，其中仅顶层原子显示为球，基元原胞的超胞采用矩阵表示法给出的。晶体学原胞（该体系惯用单胞与晶体学原胞相同）显示为蓝色。

Wood表示法将超胞的晶格矢量长度原始晶格长度的倍数，然后搭配一定的转动度数。当创建超单元时，晶格矢量的夹角γ不发生改变时，才能使用Wood表示法。

Wood表示法的例子：

• Cu(100) (2×1)
• Cu(100) (√2×√2)-R45°
• Cu(111) (√3×√3)-R30°

AMSinput只支持矩阵表示法，矩阵表示法，可以表达任何倍数超胞，也能表达转动角度。上述Wood表示对应的矩阵表示为：

AMSinput中❄ → Generate SuperCell输入超胞矩阵即可。对于六边形表面，如Cu(111)或ZnO(0001)，γ=120°，使用γ=90°的超晶胞通常更方便，特别是用于分子动力学模拟的大Slab。

将六边形（γ=120°）表面晶胞转换为矩形（γ=90°）超胞：

六方体晶体（例如纤锌矿ZnO）可以类似地转化为正交超晶胞：

表面能E$^{surf}$给出了每单位面积大小的表面，形成所需的能量。如果是上下两个表面相同的Slab，能量E$^{slab}$可分为“体积”贡献和“表面”贡献：

\[E^{\text{slab}} = \frac{n^\text{slab}}{n^\text{bulk}}E^\text{bulk} + 2AE^\text{surf}\]

其中n$^{slab}$是Slab晶胞中配方单元的数量，n$^{bulk}$是体相晶胞中的配方单元数量（二者比值即原子个数比值），A是表面积，E$^{surf}$是表面能，因子2是因为Slab暴露出两个表面。E$^{slab}$和E$^{bulk}$都可以由BAND、QE、VASP等直接计算出。

随Slab增厚，绘制E$^{slab}$ v.s. n$^{slab}$曲线，拟合出一条直线，y轴的截距即E$^{surf}$（厚度为0，n$^{slab}$=0）。

这种方式中，E$^{bulk}$不需要明确（实际上对应于着直线的斜率），并且要求对每个厚度的Slab使用相同的k点采样设置（不是指k点个数一致，而是指k点密度一致，因此可以使用Good、Very Good等针对密度的定义），这意味着表面计算中的k点不需要进行收敛测试。

上面的Slab能量公式可以变形得到表面能E$^{surf}$计算公式：

\[E^{\text{surf}} = \frac{E^\text{slab} - \frac{n^\text{slab}}{n^\text{bulk}}E^\text{bulk}}{2A}\]

这就要求E$^{bulk}$和E$^{slab}$计算非常精确，需要对这两个能量的计算，进行k点收敛性测试，找到足够大的k点。且E$^{bulk }$需要非常仔细地进行结构优化，收敛标准应比一般略高。

注意：

• 从优化后的体相结构计算表面能是最容易的，这样在生成的Slab的横向上就不会产生应变。
• 如果是先模拟体相结构的表面，那么不要优化生成的Slab的晶格。只需要优化体相晶格，然后用它来构建Slab。
• 如果是模拟薄膜，就可以考虑优化Slab的晶格。

这里的ZnO(\(10\bar{1}0\))表面能计算使用DFTB来演示，结构是用DFTB优化好的。

1. AMSinput切换到DFTB模块，Model：SCC-DFTB
2. Parameter directory：DFTB.org/znorg-0-1
3. Kspace：Normal
4. Task：Geometry Optimization
5. Details → Geometry Optimization勾选Optimize lattice
6. ❄ → Generate Slab，Miller indices：1、0、0 (对六角晶系可以忽略第三个指数-1），Number of layers：4
7. 选中Zn(2)原子从而得到正确的表面终端
8. 点击Generate Slab按钮
9. Model → Lattice注意表面积area：\(|\vec{a} \times \vec{b}|\) = 17.68 Ų
10. 保存并提交作业，计算完毕SCM → Output，窗口底部搜索calculation results即可看到体系的能量
11. 对不同厚度的Slab重复上述过程

然后，将每个厚度的Slab能量E$^{slab}$制成表格：

根据方法一，这些点拟合出来的直线：

可以看到截距 = 2AE$^{surf}$，其中A = 17.68 Ų，对不同厚度Slab面积A均为该值。因此表面能:

\[E^\text{surf} = \frac{\text{intercept}}{2A} = \frac{0.106415\ \text{Ha}}{2 \times 17.68\ \text{Å}^2} = 0.00301\ \text{Ha/Å}^2 = 1.31\ \text{J/m}^2\]

如果使用方法2，k-space设置为Normal的最厚Slab，则E$^{surf}$变为1.04 J/m$^2$。当用Good k-Space重新优化最厚的Slab时，E$^{surf}$变为1.31 J/m$^2$，这与我们在上面使用方法1得到数值相同。

小贴士：

• 熟悉Python的用户，可以下载PLAMS环境下的python脚本以批处理设置、运行和后处理生成上述表面能图和数值。另请参阅：PLAMS编写Python脚本教程。
• 基于SCC-DFTB/znorg-0-1，使用方法1和方法2计算ZnO(\(11\bar{2}0\))的表面能。答案：E$^{surf}$=1.39 J/m$^2$。

解理能是从体相中产生Slab所需的能量，如果Slab暴露出两个相同的表面，则解理能是表面能的两倍。两个表面不一样，则只能计算二者之和，例如ZnO(0001)的Zn端面和ZnO(\(000\bar{1}\))的O端面一起的解理能，而不是单个表面能。

通过旋转和平移分子，可以将吸附分子以任意构象放置在表面（基质）上，有的表面存在高对称的吸附位点，如顶位、桥或穴位。

Cu(111)结构复制粘贴到AMSinput。为了便于可视化，我们将第二个Cu层渲染成绿色，第三个Cu层渲染成黄色：

1. AMSinput → View → View Direction → Along x-axis（转为侧视图）
2. Atoms → Details (Color, Radius, Mass, …)
3. 在第二层中选择一个原子。在右侧，单击该原子的棕色矩形框，并将其改为绿色
4. 选中所有第二层Cu原子，右键点击刚才的绿色矩形框，点击Use the same value for all selected atoms，则所有第二层Cu原子都成了绿色
5. 类似将第三层Cu原子设置为黄色
6. View → View Direction → Along z-axis（转为俯视图）

侧视图和俯视图如下：

1. Edit → Tune Geometry
2. 选中水分子，点击右侧窗口顶部的 ⬭ 按钮设置要调整的原子，此时水分子的三个原子成为了Tune Geometry要移动的对象
3. 选择顶层的任意3个Cu原子，点击Plan右侧的 ⬭ 按钮，从而设置表面所在平面
4. 设置Distance to plane为2.5 Å，后面的选项改为to closest atom，用于控制表面吸附的分子到表面的距离
5. 调整不同的吸附位点：
   • 顶位：选择顶层的一个Cu原子
   • 桥位：选择顶层两个相邻的Cu原子
   • fcc穴位：选择顶层三个相邻的Cu原子，俯视图中可以看到其中心围着一个黄色（第三层）Cu原子
   • hcp穴位：选择顶层三个相邻的Cu原子，俯视图中可以看到其中心围着一个绿色（第二层）Cu原子
6. 点击Above atoms右侧的 ⬭ 按钮.
7. 点击Move按钮，吸附质就按照距离和位点移动到该位置了。

如下图所示：

小贴士： 在Tune Geometry窗口中，您还可以：

• Translate atoms, direction：沿某个方向一定距离平移原子。使用 ⬭ 按钮来定义选定原子的方向，或者手动输入平移量xyz。单击◀ ▶ 按钮来移动原子。
• Translate atoms, target：向某个目标位置平移原子。使用 ⬭ 按钮，基于当前选定的原子定义目标位置。单击◀ ▶ 按钮来移动原子。
• Rotate atoms：以Center为中心围绕Axis轴旋转原子。可以通过一个或2个原子来定义Center、Axis，也可以手动输入。例如，将Angle设定为20°，将“轴”设定为0，0，1，将“中心”设定为0，0，0。那么原子将围绕z轴以20°为步长转动。
• 如果将鼠标光标放在输入字段上，您将在弹出的帮助信息中获得更多信息。

吸附能定义为吸附分子与基底的结合能：

\[ΔE^\text{ads} = E^\text{slab} + E^\text{mol} - E^\text{slab+mol}\]

右侧的所有能量都是几何优化过的能量，正吸附能意味着放热吸附。也可以使用不同的方式约定：

\[ΔE^\text{ads} = E^\text{slab+mol} - (E^\text{slab} + E^\text{mol})\]

如此，负吸附能意味着放热吸附。在报道吸附能时，请务必指定用于计算吸附能的公式。

上图a到c表示需要计算能量的Slab、分子、Slab+分子，d表示远距离未吸附状态，也可以用这个体系的能量表示Slab能量与分子能量之和。

在BAND中，如果基组不够大，吸附可能会因基组重叠误差（BSSE）而“显得”稳定。解决BSSE问题的最佳方法是增大基组，或者进行BSSE矫正。DFTB、Quantum ESPRESSO、VASP一般没有显著的BSSE。

示例：水分子在ZnO（101’0）上的吸附。

下图显示了包含/不包含BSSE校正的情况下计算的吸附能，每个几何结构都用相应的基组进行了优化，并对最终的Slab+分子体系进行了BSSE校正。

BAND的BBSE矫正参考教程：Ghost原子的设置与使用（用于周期性体系的BSSE修正）

在ReaxFF中，采用电负性均衡方法（EEM）确定原子的电荷，进而确定库仑能对总能量的贡献。在一个孤立的分子中，电荷只能在分子内各原子之间重新分配。类似地，对于孤立的Slab，电荷只能在Slab内的各原子之间重新分布。

在Slab+分子体系中，有更多的自由度来重新分配电荷，从而可能得到更低的能量，吸附显得更稳定，这个问题类似于BSSE。因此对于ReaxFF，可以选择使用In-Cell处理方法： \[E^\text{ads} = E^\text{slab+mol widely separated} - E^\text{slab+mol}\]

其中，第一项是在一个体系里面计算的，其中Slab和分子被放置在同一个Cell中（故而称之为In-Cell），类似上图d所示。从而两项对应的两个体系，电荷具有相同的重排自由度，从而消除这种误差。类似地，也可以使用不同的符号约定：

\[E^\text{ads} = E^\text{slab+mol} - E^\text{slab+mol widely separated}\]

例子：H2O在ZnO(\(10\bar{1}0\))上的吸附。（3×2）超胞，5个双层，晶格常数a=3.28Å，c=5.28Å。熟悉Python使用的用户，可以下载计算这些数值的PLAMS python脚本测试。

两种方法之间的吸附能差异为10 kcal/mol，这源于以下事实

\[E^\text{slab} + E^\text{mol} \ne E^\text{slab+mol widely separated}\]

对于更大的超胞和更厚的Slab，这种能量差变得更大，如下图所示：

训练吸附相关的力场时，训练集要注意使用In-Cell模型。

Slab暴露出两个表面，吸附质可以添加到两个表面上，也可以只添加到一个表面上。对于Quantum ESPRESSO和VASP，双面吸附可以防止在真空间隙中形成偶极矩，不过BAND、DFTB支持没有真实2D周期性，不存在这个问题。

一般通过固定Slab的中心层或底层，达到固定Slab的目的。当弛豫干净（没有吸附物）的Slab时，最好是弛豫所有原子（不固定任何原子），或者固定Slab中心层，这样Slab的顶部和底部都会弛豫。这种弛豫过的Slab可以用于：

• 单侧吸附（放松所有原子或固定中心层和底层）
• 双面吸附（放松所有原子或固定中心层）

即使在单侧吸附（在顶表面上）的情况下，使用DFTB、BAND、Quantum ESPRESSO或VASP时，具有松弛的底表面也是好的，这有助于SCF收敛，并提供“更清洁”的DOS，是对于共价或离子材料尤其如此。对于金属，这种影响较小。

固定原子的方法：在AMSinput中，在Task设定为Geometry Optimization后，选择要保持固定的原子。Model → Model → Geometry Constrained and PES Scan，单击selected atoms (fix positions)前面的+按钮，保存作业运行即可。

建模参考：【入门基础教程】创建固体表面，并均匀加入其它气体、液体分子