Ex58 吸附能的计算（三）
前面一节我们将Cu的slab表面扩展成p(2x2), 然后放了一个O原子在上面吸附。最后发现与p(1x1)表面上吸附相比，O原子的吸附能从+1.2 降低到0.2 eV. 说明覆盖度高的时候，O原子的吸附很差。吸附能受覆盖度的影响很大。
那什么是覆盖度？ Surface coverage（θ）
从字面上理解，就是表面被吸附物(原子或分子等)覆盖的程度。一个固定大小的表面，如果上面某一吸附物越多，那么说它的覆盖度也就越大。也可以理解为：该吸附物在表面上的浓度也就越大。在多相催化的过程中，吸附是第一步。而在催化剂表面上，其活性位点数量是有限的（N个），如果N1个活性位点被吸附物种占据了，那么此时的覆盖度：θ= N1/N。如果所有的活性位点都被占据了，那么覆盖度就是1。 这个相信学过物理化学，以及相关多相催化基础知识的都能理解。 如果不知道的，请自己查找相关参考书学习，主要是在Langmuir 吸附那一部分。
在表面计算中，覆盖度我们通常使用原子数目的比来表示。首先，在p(2x2)的Cu(111)表面上，有4个Cu原子，我们可以认为活性位点是4个。当我们放一个O上去后，占据了一个Cu，那么此时的覆盖度是1/4。前面我们在p(1x1)的表面上放了一个O，可以认为是Monolayer的O吸附，也就是覆盖度为1。
1）之所以用原子数来表示，而不用表面积。是因为吸附的物种往往不像O原子这么简单，比如甲醇，乙醇，等等。这些表面上吸附物种的面积很难算，将它们的面积除以表面的面积，可行性不大。
2）师兄，如果我的吸附分子很大，比如苯环。 1个Cu原子上面放不下，那该怎么表示?
首先你需要扩胞，如上一节所讲的，将表面扩大到能容得下苯环为止。
1）比如我们将表面扩展成p(3x3)后才能放上去一个苯环，那么此时的覆盖度就是1/9。虽然数值很小，但对于苯环来说，Monolayer的吸附，按照这个定义，其覆盖度就是1/9。
2） 如果我们知道一个苯环可以占据9个表面的原子，那么我们将这9个原子看做一个活性位点，那么覆盖度就是 1/1。第一个1指的是1个苯环，第二个1指的是1个活性位点（9个原子）。
大家在计算的时候，一方面参考专业书，文献中的表示，另一方面根据自己的体系需要，只要将问题描述清楚就可以了。
继续回到上一节留的另一个问题：扩胞时有什么需要注意的事项。
1）扩胞的操作已经展示给大家了，相信这不再是什么难事了；
扩胞后的结果见下图：
2）扩胞后计算，我们需要注意的主要是KPOINTS这个文件。前面p(1x1)-Cu(111)，我们用的K点是：13x13x1。那我们扩完胞后，晶格变大了，还能继续用$13\times13\times1$吗？
答：可以用来计算。但结果不能用来和其他K点的进行比较。
什么意思呢？
这里大师兄想要表达的是：如果你想要对比两个计算结果，首先要保证它们在一个计算的等级上， 即两个计算K点的网格密度是一样的，或者差不多的。在Reciprocal space中，晶胞的长度为$2\pi/a$。然后我们将这个长度分割成k段。每一段长度为：$2\pi/ka$
我们扩胞后，晶胞长度变为2a，则在Reciprocal Space中，晶胞的长度为$2\pi/2a$ =$\pi/a$。
如果我们用之前的KPOINTS，那么就会将π/a分成k段，每段长：$\pi/a$， 由于分的更细，这会导致扩胞后K点的密度更大。
为了保持和扩胞前的密度一致，我们需要将$\pi/a$分割成 k/2 段，这样每段长为$2\pi/ka$。
因此当我们扩胞后，KPOINTS如果保持不变，会导致前后两个计算的K点密度不一致时，从而计算的精度也不一样，所以两个结果没有可比性。
师兄，你说了这么多废话。那我们该怎么办呢？
很简单，将K点数目减小就可以了。 前面p(1x1)我们使用的是$13\times13\times1$，那么扩胞成p(2x2)后，K点就需要相应的减小到$7\times7\times1$。更简单一点，回想我们前面说到的那个经验选择。K * a 的取值范围。
https://wiki.fysik.dtu.dk/gpaw/exercises/surface/surface.html
我们只要保证所有计算中： K * a 的数值差不多就OK了。所以，当你看到这里，你就可以随便扩胞了，什么p(2x2), p(3x3), p(4x4)啦，扩胞后都知道该怎么选择K点并修改KPOINTS了。
那对于其他输入文件呢？
INCAR：
此时，对于Cu体系来说，可以保持不变。但有些与原子数目相关的参数，你还是需要在INCAR里面修改，保持与新的POSCAR对应。比如：
MAGMOM, 但如果你的体系是具有磁性，需要重新设置MAGMOM。
更复杂点的体系，如果你的体系是不仅有磁性，还是反铁磁，当你扩胞后，你还要给不同的原子设置磁矩的方向。
POTCAR： 保持不变
提交任务的脚本： 扩胞后，由于体系变大， 计算会变慢，为了加快速度，我们可能会需要增加节点数目。
计算步骤：
老实人的做法：
扩胞后，很多人问计算吸附能的步骤，不清楚该怎么样进行。
第一步：优化扩胞后的Slab 模型： 因为扩抱前已经relax过了，所以这一步很快就会结束；得到新的slab能量。
第二步：在第一步优化的基础上搭建吸附模型； 然后优化，得到slab+吸附物的能量。
第三步：套公式，计算吸附能。
老司机的做法
当然了，如果你经验丰富，计算的够多了，可以尝试这样做：
第一步：优化扩胞后的Slab 模型： 因为扩抱前已经relax过了，所以这一步很快就会结束；得到新的slab能量：
同时：在扩胞后的slab模型上（没有优化）直接搭建吸附模型，然后优化，得到slab+吸附物的能量。
第二步：套公式，计算吸附能。
这两种做法，最后的结果应该是一样的，如果有差别，也是微乎其微，可以忽略不计。注意：在计算纯的slab能量时，表面应该是放开的。
另一个很常见的问题是，slab的表面我应该放开还是固定呢？
严格来说，slab的表层原子应该是弛豫的，也就是坐标后对应的是 T T T，这样更符合物理化学意义，因为催化剂表面的原子不可能会不动。厉害的时候，吸附物还会导致表面重构。所以你在阅读文献，会发现理论部分都会写：Top two layers 放开来模拟表面，底部的原子固定来模拟体相等等这样类似的话。
但是，将表面原子放开后，会导致结构优化步数的增加，从而引起计算量的增加。如果你的模型很大，吸附结构很多，可以先从小模型，低K点出发，先固定住表面，算一遍吸附，大体有个感觉后，再一步一步增加计算的精度。
如果体系太大了，可以尝试着固定表面来算吸附能，步骤如下：
- 1）先优化slab模型，这一步表面原子是要放开的；
- 2）固定表面原子，
- 3）放吸附物种优化。
目前计算能力基本上可以完全满足放开表面原子的计算。所以，大家不要纠结这个问题。表面的原子能放开，尽量放开。如果实在是不放心，那么就找一篇和你工作很相近，体系类似的计算文章，看看别人怎么设置的，照猫画虎就可以了。
参考阅读：
Density Functional Theory：A Practical Introduction, Chapter 3：Nuts and bolts of DFT calculations
扩展练习：
另一个无关痛痒，但对后面工作效率会造成影响的问题就是：
当扩展完之后，你会看到，Cu原子的坐标，是按照每个单胞中的坐标来排列的。在p(2x2)中有4个p(1x1)的单胞，先列出来第一个(1x1)中4个原子的坐标，然后再列出来第二个(1x1)的，以此类推。但是，比如我们想仅仅放开表面的第一层原子（也就是纵坐标为6.26的所有原子），就得把纵坐标为4.189的原子挨个由T T T 改成 F F F （图中箭头所指的部分）。这样会很麻烦，还浪费很多时间。
怎么解决这个问题，快速将表层，次表层原子固定或者放开呢 ? 看下一节的内容。