【计算教程】在Pd(110)表面CO的吸附计算

该教程是CASTEP官方教程中有名的劝退教程。才读到开头就会发出这样的疑问：这是啥？这又是啥？为什么不能直接告诉我怎么算个吸附能就完事了？问题很简单就是涉及到的知识点不少。为了让大家能够更容易接受和吸收，我自己往里面加了一些东西。

教程目的：介绍使用CASTEP来计算气体在金属表面上的吸附能。

所用模块：MaterialsVisualizer, CASTEP阅读时间：★★★☆☆

阅读难度：★★★★☆

预备知识：吸附模型相关知识

在本教程中，将研究CO在Pd（110）上的吸附。Pd表面在各种催化反应中起着至关重要的作用。了解分子如何与此类表面相互作用是了解催化反应的第一步。在本例中，密度泛函理论(DFT)模拟能促进这一理解，可以解决以下问题：※ 分子吸附在何处？※ 多少分子会吸附在表面上？※ 什么是吸附能？※ 吸附结构如何？※ 吸附的机理是什么？

我们将关注一个吸附位点，即短桥位（the short bridge site），因为它是已知的能量优先位置。大家如果不是研究吸附的，可能会对这个概念比较懵，这里我给大家找了个简图。

在1 ML覆盖度下，CO分子彼此排斥，从而阻止了它们与表面完全垂直对齐。MonoLayer的缩写,ML,是覆盖度的单位。其表示的含义为吸附上的原子占基底原子的比例。例如基底是由100个原子组成的薄膜，它又吸附了100个其他原子在其表面形成一层新的薄膜,此时的覆盖度就是1.0(ML)还是懵，那我再来一张图（AppliedSurface Science 257 (2011) 6878–6883）

以上不是MS原版教程中的内容，是我自己为了方便大家理解加进去的。如果还是不懂，就找书来看看吧，一般讲Langmuir 吸附的部分多少会提到这个。我们将通过考虑一个（1×1）和（2×1）表面原胞来计算这种倾斜对化学吸附能的能量贡献。

Pd的bulk和Pd（110）表面的顶视图。（110）劈裂平面以蓝色突出显示。a0是体晶格常数，也称为晶格参数。以下开始本教程实战部分介绍：在本教程中，将使用CASTEP优化并计算几个不同系统的总能量。一旦确定了这些能量，就可以计算Pd（110）上CO的化学吸附能。 本教程包含以下内容1 Getting started（知道是什么意思，但依旧不知道该翻译成啥比较好）2 优化Pd晶格3 构建并优化CO4 构建Pd(110)表面5 弛豫Pd(110)表面6 添加CO分子到1x1 Pd(110)表面并优化结构7 建立并优化2 x 1 Pd(110)表面8 能量分析9 态密度分析

1 Getting started

打开MS并创建一个新的project，打开New Project对话框，并且输入CO_on_Pd作为project名字。然后点击OK按钮。本教程包括五个不同的计算。为了使项目容易管理，应该一开始就在项目中建立五个子文目录。在Project Explorer中根目录图标上右击，选择New | Folder。重复该操作四次。在NewFolder上右击，选择Rename，输入Pd bulk。对其他文件夹重复这一操作，并将它们分别命名为Pd(110), CO molecule, (1x1) CO on Pd(110) and (2x1)CO on Pd(110)。

2. 优化Pd晶格

Pd的晶体结构包含在Materials Studio提供的结构库中，直接导入即可。

在Project Explorer中，右击Pd bulk目录，并选择Import...，定位到Structures/metals/pure-metals，找到并导入Pd.msi。

在Pd 3D视图文件中右击，选择Display Style，在Atom选项卡上，选择Ball and stick，关闭对话框。

现在用CASTEP对Pd的几何形状进行优化。

学完上一节教程，这个应该已经会了吧（【计算教程002期】CASTEP—用第一性原理方法预测AlAs的晶格参数）

从菜单栏选择Modules | CASTEP | Calculation。

显示CASTEPCalculation对话框。

晶体的晶胞优化，比采用默认设置执行的优化需要更精确的计算。将Quality从Medium改为Fine。将Task从Energy改为Geometry Optimization。几何优化的默认值不包含对晶胞的优化。点击More…按钮，在CASTEP GeometryOptimization对话框里，选中Optimize Cell。点击Run按钮。显示一个转换到primitive cell的消息框，点击Yes按钮。此处新版和老版有差别，我出差了本本中是老版的MS官方教程，是这样选择，新版怎么选择请看上一期教程，切记看文字描述部分Cell optimization中选择full。提交并开始运行。此时可以转入下一节，构建CO分子。但是在计算完成时应该返回这里以显示晶格参数。这句意思是，不要以为做计算化学就是整好输入文件就可以开心玩一周。你应该在计算等待的时间为下一步做准备，或者看看文献看看书当优化结束时，必须把原胞转换回常规晶胞表现形式，以便在第4步时构建Pd(110)表面。在Project Explorer中，打开位于Pd CASTEP GeomOpt目录下的Pd.xsd文件。从菜单栏选择Build | Symmetry | Conventional Cell。从菜单栏选择File | Save Project，然后选择Window | Close All。保存是个好习惯，记不住的同学别担心，你会在某次宕机后记住的打开位于Pd CASTEP GeomOpt目录下的Pd.xsd文件。显示的是Pd的优化结构。在3D视图中右击，选择Lattice Parameters。a值大约是3.940 Å，相比较实验值是3.89 Å然后关闭 Lattice Parameters对话框和Pd.xsd.

3. 构建并优化CO

CASTEP只对周期系统起作用，所以，为了优化CO分子的几何结构，必须把它放到晶格中。在Project Explorer中，右击CO molecule目录，选择New | 3DAtomistic Document。右击 3DAtomistic.xsd，选择Rename，输入CO，按回车键。此时显示一个空白的3D视图。从菜单栏选择Build | Crystals | Build Crystal。选择LatticeParameters 选项卡，把晶胞长度值a、b 和c都设为8.00，点击Build 按钮。在3D视图中显示一个空的晶胞。从菜单栏选择Build | Add Atoms打开添加原子的对话框CO分子中C-O键长通过实验测量为1.1283Å，如果我们通过使用笛卡尔坐标添加原子，可以很精确地按照这个键长值建立CO分子。在Add Atoms对话框中，选择Options选项卡。确保Coordinate system被设为Cartesian。选择Atoms选项卡，按下Add 按钮。一个碳原子被添加到晶胞的原点。在Add Atoms对话框上，把Element改为O。保留x和y的值为0.000，把z 值改为1.1283。点击Add 按钮。然后关闭对话框。现在已经准备好优化CO分子了。从菜单栏选择Modules | CASTEP | Calculation保持以前的计算设置不变，但这次不必优化晶胞。在Setup选项卡上，点击More...按钮，取消选择Optimize Cell并关闭对话框。选择Electronic选项卡，把k-point设置由Medium改为Gamma。以后的教程会告诉大家k点是啥，以及这里为啥这么选取。现在就先照做完事了选择Properties选项卡，选择Density of states。把k-point 设为Gamma，勾选上Calculate PDOS。按下Run 按钮。这个算很快，先不着急分析数据，先整下一步。当然，保存操作还是要做的。

4. 构建Pd(110)表面

在本节中要使用前面优化的Pd结构。在Pd bulk/Pd CASTEP GeomOpt目录中打开Pd.xsd文件。创建一个表面需要两步。第一步是劈开表面，第二步是创建一个包含表面的真空层。从菜单栏选择Build | Surfaces | Cleave Surface。把Cleaveplane (h k l) 从-1 0 0 改为1 1 0，按TAB 键。增加Fractional Thickness值到1.5。点击Cleave按钮，关闭对话框。这步就是劈开表面此时打开一个新的3D 模型文件，它包含一个2D周期性表面。可是CASTEP需要一个3D周期性系统作为输入，它可以使用Vacuum Slab工具获得。想一下，我们为什么要2D结构，又是为什么要把2D表面做成3D模型？这个问题有点无厘头，细细想想你会感谢我这个傻问题从菜单栏选择Build | Crystals | Vacuum Slab。将Vacuumthickness值从10.00改为8.00，点击Build按钮。结构由2D周期改变为3D周期，并且有一个真空被加到原子的上方。在做接下来的事情之前，必须重新定位晶格。在3D视图中右击，从快捷菜单中选择Lattice Parameters，选择Advanced 选项卡，点击Re-orientto standard按钮，关闭对话框。此时改变晶格显示方式并旋转结构，使Z轴垂直于屏幕。在3D视图中右击，从快捷菜单中选择Display Style。选择Lattice选项卡。在Display style区域，将Style由Default改为Original。关闭对话框，按两次向上的方向键。3D视图就如下图所示。

具有最大Z坐标的Pd 原子被称为“最高Pd层”（the uppermost Pd layer）。在本教程的后面，需要知道Pd的层间距d0，可以使用原子坐标计算它。从菜单栏选择View | Explorers | Properties Explorer。选择FractionalXYZ为x = 0.5和 y = 0.5的Pd原子。注意该原子的z值是指来自于XYZ属性的值。Z值应该约为1.391 Å，这是层间距。这个z值是指XYZ(笛卡尔)属性的z坐标，而不是FractionalXYZ值大概就是下图那里（此图是未优化的Pd，所以数值和教程中不一样，我只是截个图）

注意：对fcc(110)系统，d0可以以下公式计算。

在弛豫表面之前，必须约束（constrain）Pd内部的原子，因为只需要弛豫表面。限制（restraint）、冻结（freeze）和约束的概念是不同的，请大家一定要注意按住SHIFT键，选中除了最顶层Pd原子之外的所有Pd原子。从菜单栏选择Modify |Constraints。选中Fix fractional position，关闭对话框。此时Pd内部的原子已经被束缚了，可以通过改变显示颜色查看被束缚的Pd原子。在3D视图中点击，以取消选择原子。右击选择Display Style。在Atom 选项卡中将Color by选项改为Constraint。此时3D视图如下图所示。

然后将Color by选项改回Element，关闭对话框。这个结构是做Pd(110)表面弛豫所必须的，同时也是(1x1)CO在Pd(110)表面优化的初始模型。从菜单栏选择File | Save As...，定位到Pd(110)目录，点击Save按钮。按住ctrl键然后拖拽这个文档进入(1x1) CO on Pd(110)文件夹中，并且重命名该文档为(1x1) CO on Pd(110)。然后选择File | SaveProject，然后then Window | Close All。 

5 弛豫Pd(110)表面

现在准备优化Pd (110)表面。从Project Explorer打开Pd(110)目录中的Pd(110).xsd。从菜单栏选择Modules | CASTEP | Calculation。将Task从Energy改为Geometry Optimization。按下More…按钮，确保Optimize Cell未被选中。关闭对话框。

为了保持即将执行的计算的一致性，需要在Electronic选项卡里做一些改变。选择Electronic选项卡，点击More…按钮。从CASTEP Electronic Options 对话框选择Basis选项卡，选中Use custom energy cutoff，并把其值改为300.0。

选择k-points选项卡并选中Custom grid parameters。在Grid parameters中，设置a值为3，b值为4，c值为1。关闭对话框。

还需要把系统的态密度计算了。在CASTEP Calculation对话框上选择Properties选项卡，选中Density of states。选中Calculate PDOS，并把k-point set改为Medium。现在已经准备好计算。点击Run按钮，然后关闭对话框。运行这个计算需要花费一些时间，我们将在最后对结果进行分析。现在应该继续构建下一组表面。从菜单栏选择Save Project，然后选择Window | Close All。

6 添加CO分子到1x1 Pd(110)表面并优化结构

现在我们将处理(1x1) CO on Pd(110)中的结构。这里将是劝退重灾区。在Project Explorer中，打开(1x1) CO on Pd(110)目录中的(1x1) CO on Pd(110).xsd文件。现在把CO分子添加到短桥位的上方，我们将在建模过程中使用实验方法测定了的键值。

上图为在yz平面中CO在Pd(110)面上的几何形状。第一步是添加碳原子。Pd-C键长(上图中用dPd-C表示)为1.93Å。当使用添加原子工具时，可以输入笛卡尔坐标，也可以是分数坐标，但在本例中，要使用分数坐标，xC、yC和zC。xC和yC的值是简单的分别为yC=0.5和xC=0。然而zC是比较困难的，可以从zPd-C和zPd-Pd这两个距离计算得到其值。这段看不明白，一定要上机边看边学。也不要略过上一期只看这一期。一定要边看边学！一定要边看边学！一定要边看边学！zPd-Pd可以简单地由晶格参数的a0除以得到(应当为2.77Å)。为方便大家理解，把这个图在这里再放一次

zPd-C可以由下式获得，其值应为1.35 Å。

zPd-C和zPd-Pd相加得到zC，其值约为4.12Å，可以使用Lattice parameters把这个距离转换成分数长度。在3D视图中右击，从快捷菜单中选择Lattice Parameters。注意c的值。为计算分数z坐标，可以简单地用晶格参数c除zC得到，其值约为0.382。从菜单栏选择Build | Add Atoms，选择Options选项卡。确保Coordinate system是Fractional。选择Atoms选项卡，把Element改为C。把a值改为0.0，b值改为0.5，c值改为0.382。点击Add按钮。如果想确认所建立的模型是否正确，可以使用Measure/Change工具。具体做法是：在工具栏上点击Measure/Change工具

旁的选项箭头，选择Distance。点击Pd-C键。下一步是添加O原子。在Add Atoms对话框上，把Element改为O。在实验中C-O键长已经被测定了，为1.15 Å。在分数坐标系中是0.107，把这个值与C的分数z坐标值(0.382)相加，得到O的分数z 坐标值0.489。把c的值改为0.489，然后点击Add按钮，关闭对话框。在Pd表面计算时使用的是是P1对称性，你可能经常听说make p1,把这个当成一句咒语。其实……（Crystals with P1 symmetry only possess translational periodicity.）但是加入CO分子后，系统有了更高的对称性。可以使用Find Symmetry工具找到并加上对称性，正确使用对称性特征可以加快计算速度。量子化学第一性原理计算都有类似的操作，但同时也要提醒，不该有对称性的时候，千万别施加类似约束，甚至在某些计算情形中要手动移动一些原子以破坏对称性。在工具栏上选择Find Symmetry工具

，点击Find Symmetry 按钮，然后点击Impose Symmetry按钮。对称性约束为PMM2。在3D视图中右击，从快捷菜单中选择Display Style。选择Lattice选项卡，将Style改为Default。在Atom选项卡上，选择Ball and stick显示样式，然后关闭对话框。现在的结构如下图所示。

在优化结构的几何形状之前，应该将结构保存到(2x1) CO on Pd(110)目录中。从菜单栏选择File | Save As...，定位到(2x1) CO on Pd(110)目录，把文件保存为(2x1) CO on Pd(110).xsd。现在已经准备好了优化结构。从菜单栏选择File | Save Project，然后选择Window | Close All。在Project Explorer中，打开(1x1)CO on Pd(110)文件夹内的(1x1)CO on Pd(110).xsd。从工具栏选择Modules | CASTEP | Calculation。前面为计算设置的参数在这里直接套用。点击Run按钮。在计算进行过程中，可以继续构建最后的结构。

7 建立并优化2 x 1 Pd(110)表面

第一步是在(2x1) CO on Pd(110)目录中打开3D Atomistic文件。在Project Explorer中，打开(2x1) CO on Pd(110)目录下的(2x1) CO on Pd(110).xsd 文件。当前是1 x 1的晶胞，所以需要使用Supercell工具把它改成2 x 1晶胞。从菜单栏选择Build | Symmetry | Supercell，把b值增加到2，点击Create Supercell按钮。关闭对话框。CO在Pd(110)面的(2 x 1)晶胞结构如下图所示。

现在使CO分子彼此翘起。为简化这一操作，把处于y=0.5的CO分子标记为A，y=0.0的CO分子标记为B。选择B分子中的碳原子。在Properties Explorer中，打开XYZ属性，把X的值减去0.6。对B分子中的O原子重复这一操作，但是把X值减去1.2。现在对分子A重复上述操作。选择A分子中的碳原子。在Properties Explorer中，打开XYZ属性，把X值加上0.6。对A分子中的O原子重复此操作，但是把X值加上1.2。翘起后沿z方向的结构如下图所示。

然而，你应该注意到了Pd-C和C-O键长已经从初始值发生了改变。选择A分子中的碳原子，使用Properties Explorer，把FractionalXYZ 属性中的Z改为0.369。对B分子重复这一操作。这将修正Pd-C键长，也可以使用Measure/Change工具来修正C-O键长。在工具栏上点击与Measure/Change工具 

旁的选项箭头，选择Distance，点击分子A的C-O键。在Properties Explorer内，把Filter改为Distance。把Distance属性改为1.15Å。对分子B重复这一操作。现在重新计算系统的对称性。在工具栏上选择Find Symmetry工具

，点击Find Symmetry按钮，然后点击Impose Symmetry按钮。对称性为PMA2。在原胞中Pd表面的CO分子由3个变为2个。现在已经可以开始优化几何形状了。从菜单栏选择Modules | CASTEP | Calculation。我一直没有提菜单上

也可以进行计算什么的，感觉教茴香豆有几种写法意义不大，而且有点无厘头。点击Run按钮。计算开始。当计算结束的时候，需要提取系统的总能量，就像在下一节中详细介绍的那样。可以转到下一节提取前面计算的能量。终于开始处理之前的能量数据了。

8 分析能量

在本节中将要计算化学吸附能（chemisorption energy）ΔEchem。它的定义为：

让CO分子相互倾斜，由此减少CO分子的自我排斥，将会使能量增加。排斥能（repulsion energy）可由下式计算：

为了计算这些属性，需要从每个计算的CASTEP文本输出文件提取总能量。在Project Explorer中，打开CO molecule/CO CASTEP GeomOpt目录中的CO.castep文件。按CTRL+F键搜索 Final Enthalpy，记录该值。重复这一操作，找出其它系统到总能量，不同系统的总能量列于下表中。

一旦有了这些能量，利用上面的公式就可以计算化学吸附能ΔEchem和排斥能ΔErep。它们的值分别约为1.96eV和 0.1eV。

9 分析T态密度

甭管有没有用，反正吸附了，态密度总要有点变化吧，变化的细节便是你的发现，你计算研究的意义。牛人可是能从能带图态密度等数据看出整个宇宙。接下来，您将检查状态密度（DOS）的变化。这将使您深入了解CO在Pd（110）上的键合机理。为了做到这一点，需要显示孤立CO分子的态密度和(2x1) CO on Pd(110)结构的态密度。在Project Explorer中，打开CO molecule/CO CASTEP GeomOpt目录下的CO.xsd 文件。如果下面点击分析的时候找不到Density of states，回头看看是不是没有打开CO.xsd从菜单栏选择Modules | CASTEP | Analysis。选择Density of states。选中Partial，取消选择f和sum，其他设置不变，点击View按钮。显示一个图表文件，示出了CO分子的PDOS。

对(2x1) CO on Pd(110).xsd文件重复这一操作。

很显然，孤立的CO分子在大约20、-5和-2.5eV处的电子态显著地低于吸附到表面的CO分子。这里有个note我也给翻译一下（挺担心翻错，大家自己有兴趣看一下原文）Pd的默认赝势（不同MS版本文件名不同）将4s和4p半芯态（semicore state）视为价层。这会在计算出的DOS中在-84和-49 eV处产生尖峰。上图不包括这些态；这可以通过使用Properties Explorer沿X和Y轴更改最大值和最小值来实现。作为一项独立的练习，您可以通过研究由C和O原子引起的吸附物络合物对PDOS的贡献来进一步分析PDOS。按住SHIFT键，选择(2x1) CO on Pd(110).xsd文件中的所有C和O原子，生成PDOS，它显示了与Pd态杂化的影响，能级加宽并且通常向低能量转化。还有一些与本次计算不相关的内容就不贴了。

科研很苦，坚持很苦

万事屋愿你顺利

--END--