Learn VASP The Hard Way

The 2nd edition

Ex63 吸附能的计算（八）

发表于 2018-10-11 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

在进行本文的内容之前，大家先玩个人尽皆知的小游戏：

找不同

（一）找出下面两幅图中不同的部分。

相信这个大家很快都能找出来，找不出来的话，不建议继续浏览下面的部分。

（二）下面两个截图中，找出不同的部分

（三）下面两个截图中，找出不同的部分

聪明美丽能干活泼善良可爱宇宙超级无敌的你会发现：

在找不同2中，两个图中内容的区别在于： POSCAR 的第6行
在找不同3中，两个图中内容的区别在于： POTCAR 的第二个元素

当你发现了这些不同之处，就应该差不多明白了，这是2个计算任务文件，一个是C原子在Cu(111), 另一个是O原子在Cu(111)上的计算。

仔细分析不同的地方，你就应该想到: 既然我们有了O吸附相关的计算文件，那么就可以在此模板基础上，快速编辑C的吸附模型，然后进行计算。

在进行下面的快速编辑模型的内容，先给大家介绍一个Linux下面常用的比较不同的命令：diff。简单操作如下，更高级的大家自由去百度学习，并加以发挥拓展。

批量搭结构

快速准备C原子在Cu(111)上吸附的计算文件：

1）INCAR 保持一样
2）KPOINTS 保持一样
3）POSCAR: 将 O 直接替换成C就OK了。
i）直接编辑，无数种编辑器任你选。不过 Windows下一定要记得编辑完后运行下: dos2unix
ii）使用sed命令： sed -i '6s/O/C/g' POSCAR

4）POTCAR:

i) 使用前面我们介绍的脚本：potcar.sh，脚本运行例子如下图：（去本书的附录章节里面找脚本）

https://www.bigbrosci.com/2017/12/21/A05/

直接输入脚本名，后面跟你想要的元素名。

看下本人在天河II号服务器上的设置。

1） potcar.sh 在本人账号的： ~/bin 目录下
2） VASP所有的POTCAR文件在： ~/bin/pot 目录下
3） pot目录下，每个元素的文件夹中有2两个文件：POTCAR 和 PSCTR （PSCTR这个我也不知道有什么用。）
4）新手，没有经验，不会用超算的，不会Linux的，可以效仿下：

脚本下载后，将POTCAR所在的目录替换成你的： path=”你的目录”。
怎么知道目录呢？使用上图中的pwd 这个命令。

ii) 我很多时候都懒得去看POSCAR中的元素以及顺序，便打算根据POSCAR中的内容生成对应的POTCAR，于是便写了一个脚本：pospot.sh。使用方法如下：

上图运行了这个脚本两次：
1）第一次运行的时候，目录下没有POTCAR，于是脚本直接生成新的
2）第二次运行的时候，目录下有了一个POTCAR，于是脚本将POTCAR中的元素顺序读出来，然后删掉旧的POTCAR，重复1）的操作，生成新的POTCAR。（大师兄本人的扯淡逻辑，大家可以忽略，脚本文章末尾下载直接使用即可。）
3）注意：这个脚本读取的仅仅是以元素命名文件夹中的POTCAR。什么意思呢？举个例子：Cu的POTCAR有：Cu，Cu_GW, Cu_pv，Cu_sv_GW这四种。脚本读取的仅仅是Cu中的POTCAR。如果你想用Cu_pv，那么请使用： potcar.sh Cu_pv C H O这个命令。
4）脚本内容如下：大家可以自己手动敲一遍，找找感觉。

#!/usr/bin/env bash  
# show old POTCAR 
if [ -e POTCAR ]; then
  echo 'old POTCAR containes....'
  grep TIT POTCAR| awk '{print $4}' | xargs
# Remove old POTCAR  
  echo '************************'
  echo "Removing OLD POTCAR....."
  echo '************************'
  rm POTCAR -f
fi

# generate new POTCAR 
echo "Generating NEW POTCAR..."
echo '************************'
potcar.sh $(sed -n 6p POSCAR)
echo "Done"

# compare elements in new POTCAR and POSCAR
echo '************************'
echo 'NEW POTCAR containes....'
grep TIT POTCAR| awk '{print $4}' | xargs
echo '************************'
echo "Elements in POSCAR" 
echo '************************'
sed -n 6p POSCAR

5）提交任务的脚本保持一样

前面所说的都是一些具体的做法，用化学的思维来说，应该叫基元反应。当你掌握了这些具体的做法后，可以通过批处理的方式将它们集成起来。本人的做法如下：

1）批量编辑POSCAR：

1	sed -i '6s/O/C/g' */POSCAR

i）将O原子的计算直接复制成C原子的
ii）使用sed 批量编辑POSCAR中的元素行。

2）完成前面的一步，直接批量提交任务：提交任务的时候，脚本自动根据POSCAR生成对应的POTCAR.

我是怎么做的呢？

在~/.bashrc 文件中将pospot.sh 写到 qsuball里面了。大家可以看下我的~/bashrc文件：

alias q='yhqueue'
alias qsuball='for i in * ; do cd $i ; pospot.sh && qsub ; cd $OLDPWD; done'
alias qdel='yhcancel'
alias fenergy='grep "  without" OUTCAR'
alias ..='cd ..'
alias ...='cd ../..'
alias ....='cd ../../..'
alias gl='grep LOOP OUTCAR'

批量删除任务

前面批量提交了任务，但这些任务已经算起来了，但这些结果对我没有用出，纯属浪费机时，所以我要杀死它们。这里介绍给大家一个简单的批量删除任务的命令： {..} 操作如下：

qdel 就是天河II号中的yhcancel命令，本人根据自己的习惯，修改了一下alias中的内容（见前面~/bashrc文件的内容）

如果你的任务不是连续的，也可以使用这个命令，只要把第一个和最后一个任务的ID号写在花括号里两个点..的两端即可。非常简单，暴力和直接。

中间如果有些任务不是你的，即使你想杀死也没有权限，系统会报错，直接忽视即可。

找不同

前面我们弄完了C原子的计算任务，下面我们再玩一次找不同的游戏

找不同4：下面2个截图中找出不同的地方：

师兄，这个题很简单，你要算H原子的吸附。跟前面的操作是一样的。

但是，请看最后一行。

仔细的你会发现，这个题跟前面的有些不一样，原因在于最底部的坐标也稍微变化了一下，从7.935变成 7.635了。这是为什么呢？
原因在于和C，O原子比起来，H原子很小，它在表面上吸附的时候，会更贴近表面。因此我们把z方向的坐标稍微修改了一下：减小了 0.3 A。这样做的话可以保证我们的初始结构更具有合理的物理化学意义，从而加快计算的收敛，节约时间。大家在搭建模型的时候，结构合理是必须要考虑的，也是脑子里要时刻思考的事情。

通过sed命令批量操作在H原子的吸附上，只能完成一半的任务，也就是把元素行中的C或者O改成H。而另一半的任务：坐标的修改，则需要大家自己手动操作一下（记得在Cartesian坐标下修改，如果你脑子转的很快，可以迅速将Direct转化为Cartesian的话，也可以直接编辑Direct坐标）。

扩展练习：

1）完成本节的所有操作；
2）计算C,H，O在Cu(111)表面上，不同位点上的吸附；
3）思考脚本的运行原理；
4）思考如何通过：（找不同 + 写脚本）来节省自己的体力；
5）思考自己平时常见的人工错误，并尝试用脚本来解决或者避免。

总结：

本文的题目是闭着眼算，不是让你真正闭着眼去提交任务，而是不通过可视化的界面进行结构的批量搭建，输入文件的批量处理以及提交任务。希望大家在搭结构的时候，多一些思考，少一些操作，在搭建合理结构的路上突飞猛进。

Ex62 吸附能的计算（七）

发表于 2018-10-10 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

前面我们得到了不同吸附位点上吸附能的顺序，但结构是什么样子的呢？我们优化完的结果对不对？这还是一个问题。因此我们需要查看一下优化完的结果。

1 获取（下载）CONTCAR

没有结构，我们看个屁啊？所以第一步就是把超算中心的计算结果下载到自己的电脑里面。这里我们说获取或者下载CONTCAR，而不是OUTCAR等其他VASP的输入文件，原因在于本人这边网速传输太慢了。所以我的策略是能量等信息在服务器里面直接获取，结构的话只下载CONTCAR。如果网速允许的话，可以把所有的计算结果下载到自己电脑里面，这样查看更加方便。

上图中，我们先挂载超算中心到本地电脑上，然后将计算目录下的CONTCAR复制到本地桌面上。（ccall 这个命令）

1	for i in * ; do if [ -e $i/CONTCAR ]; then mkdir ~/Desktop/$i; cp $i/CONTCAR ~/Desktop/$i ; fi; done

备注：
由于本人这边传输很慢，即使挂载了超算中心到本地电脑上，访问内容的时候，后台依然有数据传输。所以先下载再查看。

2 使用ASE查看结构：

ASE 是Atomic Simulation Environment的简称，下载安装见：https://wiki.fysik.dtu.dk/ase/ 本人只会在Linux下面使用，Windows用户自行解决。解决不了，我也没有办法。如果Linux用户解决不了，那么使用后面的第二种方法：p4vasp查看结构。

如果你的网速很给力，可以直接通过自己电脑进入超算中心的目录，进行下面的操作。

这个软件的优点就是：我们可以一次性打开当前目录下，所有计算的CONTCAR，从而避免了使用软件挨个导入结构查看。无形中会减少我们很多的工作量。

1	ase-gui */CONTCAR

3 使用p4vasp 查看结构

Ubuntu下面唯一推荐的软件：下面图片拍的不是很好，大家凑活着看吧。

这里，p4vasp和前面说的ASE一样，也可以一个命令打开所有的计算结果。

作为p4vasp的忠实粉丝，这也是本人唯一推荐的Ubuntu系统下查看，搭建结构的软件。

i）使用p4vasp可以非常容易地进行原子替换，平移，旋转等基本操作。
ii）可以查看VASP的结算结果，DOS，能带，优化过程等等。
iii）Windows 系统下p4vasp的功能有些弱，除了不能批量打开文件外（可能是本人不会用），其他的和Ubuntu差不多。
iv）这个软件也有很多其他细微不尽人意的地方，但不影响我们的正常使用。
如果你刚开始接触这个软件，认真用鼠标各个地方点点操作一下，查看各个功能按钮的作用。
此外，VASP官网的ppt教程中也有一些零星的p4vasp操作教程，大家可以参考一下。

4 其他软件：

当然了，不论在linux还是Windows下面，都有很多查看结构的软件，比如：Jmol，Xcrysden，Molden， VESTA， Material Studio等等。这里就不再详细介绍了，主要原因是本人不太会使用这些软件操作。目前大家需要做的就是根据自己的喜好，掌握一个软件：学会查看结构，键长，键角等信息即可。切记不可贪多，等一个软件掌握好了之后，有余力的话再去学习另一个的操作。

5 扩展练习：

1）自己优化O在Cu(111)表面上不同位点的吸附，计算吸附能
2）选择一款自己喜欢的软件，查看不同的吸附结构。
3）思考其他单原子在其他金属表面上的吸附，该如何计算？
4）思考原子在表面上，为什么不同吸附位点的吸附能不一样？

6 总结：

本节没有什么技术难度，全靠自己亲自手动操作，使用一个软件并不是一蹴而就的过程，大家先把基本的简单操作掌握了，后面再逐渐提高自己的其他技能。此外，本节学习完之后，单原子在表面上的吸附对大家来说应该不是什么困难的事情了。

Ex61 吸附能的计算（六）

发表于 2018-10-09 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

前面一节我们搭建好了p(3x3)-Cu(111)表面上不同吸附位点的O原子吸附模型。（fcc, hcp, top以及bri）也就是我们有了4个POSCAR。下面我们计算O原子在不同位点上的吸附能。根据吸附能的公式，我们还要有Cu(111) slab的能量。总共5个计算。O$_2$的前面已经算了，这里就不重复了。

强调一下吸附能的定义：就是物种吸附到表面上所放出的能量，放出能量意味着吸附能是负值。很多人把这个过程颠倒了，然后来定义吸附能，这是错误的！虽然数值一样，但表示的物理意义完全不一样。颠倒过来是脱附，而不是吸附！！！而基本概念搞不清楚，最直接的后果就是给审稿人留下一个坏印象，增加拒稿的几率。当然，基本概念正确也会被拒稿，概念搞错了文章也有可能会接收。

1 准备工作:

INCAR：和前面p(1x1)的计算一样
KPOINTS：根据前面我们经常提到的经验规则：k*a=30

链接：https://wiki.fysik.dtu.dk/gpaw/exercises/surface/surface.html
POSCAR：已经准备好
POTCAR: 根据POSCAR生成对应的POTCAR，脚本见本书的附录章节。
提交任务的脚本或者命令。

A rule of thumb for choosing the initial k-point sampling is, that the product, ka, between the number of k-points, k, in any direction, and the length of the basis vector in this direction, a, should be:

    ka ~ 30 Å, for d band metals
    ka ~ 25 Å, for simple metals
    ka ~ 20 Å, for semiconductors
    ka ~ 15 Å, for insulators

Remember that convergence in this parameter should always be checked.

2 天河II号提交任务

前面我们已经讲解了如何在天河2号上提交任务，Learn VASP 系列的前面几节也介绍了一些相关的批量操作知识。下面简单通过实例操作展示下批量提交任务的一个流程。

其中qsuball是本人在~/.bashrc文件中，给批量提交任务的命令行随便起的一个名字，详见前面关于alias的使用介绍。

1	alias qsuball='for i in *; do cd $i ; qsub ; cd $OLDPWD; done'

qsub 是我们在超算中心提交单个任务的脚本，内容如下：

#!/usr/bin/env bash
rm job_sub 
echo '#!/bin/bash' >> job_sub
echo 'export LD_LIBRARY_PATH=/THFS/opt/intel/composer_xe_2013_sp1.3.174/mkl/lib/intel64:$LD_LIBRARY_PATH' >> job_sub 
echo 'yhrun -p gsc -n 24 /THFS/opt/vasp/5.4.4/vasp.5.4.4/bin/vasp_std ' >> job_sub 
yhbatch -p gsc -N 1 -J test job_sub

将上面内容复制到一个文件里面，将文件命名为qsub，然后走一遍下面的2个步骤。
1）chmod u+x qsub (赋予可执行权限)
2) mv qsub ~/bin （放到 ~/bin 文件夹下）

3 查看结果

等待任务排队，计算，结束。一个任务结束后，我们需要做的肯定就是检查结果了。那么改怎么样检查结果呢？这里列举出本人常用的几个方法和步骤。

1）查看OSZICAR

i) 在查看OSZICAR之前，我们回顾下前面的INCAR。

这里我们设置了NSW = 500,也就是允许本计算最大的离子步数。

ii）再看下OSZICAR的末尾部分，命令： tail OSZICAR

从这里可以看出来，经过25个离子步，我们的计算就已经停止了，说明收敛了。

2）查看OUTCAR

对于VASP的所有计算来说，只要你看到这样的结果，说明你的任务已经算完了。这里你要记住：算完和收敛不是一回事。比如你设置的NSW = 500， VASP优化了500步没有收敛的话，会自动停止，也会出现上图的结果。

怎么查看收敛呢？

i) 可以对比计算的离子步和我们设置的离子步。如果实际计算得了离子步小于INCAR中的，任务就收敛了。这也是前面查看OSZICAR的目的。
ii）可以通过OUTCAR里面的关键词： reached

如果你看到最后面的： Stopping structural energy minimisation 这就表明你的优化任务已经完成。
而事实上，的确有很多人把这句话误认为是计算出错了，然后在群里面求助，本人实在是想不通。

4 删除VASP乱七八糟的输出文件

计算完成之后，VASP会出现很多平时我们用不到的文件，CHG CHGCAR WAVECAR 这三个占用的空间很大，容易把应胖撑爆。

由于计算目录下的存储空间有限，本人计算的时候，一般都默认VASP不写入这些文件。

此外，在表面结构优化的时候，EIGENCAL, IBZKPT, REPORT,DOSCAR``以及超算中心与任务相关的输出文件：slurm-XXX 文件，一般来说正常计算结束后，我们也不需要。

因此，本人写了一个小脚本，用来删除这些乱七八糟的文件。脚本运行效果如上，脚本内容如下：

#!/usr/bin/env bash
for  i  in $(find . -name INCAR); do
    cd $(dirname $i)
    rm CHG  CHGCAR  job_sub  PCDAT REPORT  slurm-*.out  WAVECAR  XDATCAR -f
    cd $OLDPWD;
done

大家可以针对自己的计算任务，修改rm XXX 哪一行中的文件名。

上图中还运行了一个脚本（ta.sh），就是提取当前目录下，所有计算结果的能量。也是本人临时写的。有兴趣的可以自己照着写一遍，然后运行下。

#!/usr/bin/env bash 
for i in *; do
if [ -e $i/OUTCAR ]; then
echo -e  $i "\t" $(grep '  without' $i/OUTCAR |tail -n 1 | awk '{print $7}') 
fi
done

5 获取能量，计算吸附能

有了上面的数据，我们就可以进行吸附能的计算，如下图：

一般来说，使用eV作为单位，精确到小数点后面2位即可。整理下表格：

从上图我们可以看出来，这四个位点的吸附能强弱顺序为： fcc ≈ hcp > bridge > top

6 扩展练习

1）学会如何判断计算正常结束，并且优化达到收敛的标准了？
2）总结并思考不同位点的吸附能的计算。
3）学会判断稳定的吸附位点。
4）思考下计算结束之后，除了我们讲的这些，还有那些需要检查的？
5）本节的三个脚本：
i) qsub （天河II号提交任务）
ii） rmvasp.sh (删除乱七八糟文件)
iii）ta.sh (获取当前目录下，所有计算任务的能量)

Ex60 吸附能的计算（五）

发表于 2018-10-08 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

前面我们发现将p(1x1)-Cu(111)的表面扩展成p(2x2)后，由于O的覆盖度降低了，O原子的吸附能能+1.2 eV降低到 0.2 eV作用，说明O原子更加容易吸附在表面上了。当你仔细观察p(2x2)的表面，你会发现，表面上不仅仅有Cu原子的上方可以放O原子，还有其他的位点。如下图：

如果我们想分析下不同位点对吸附能的影响，首先我们需要搭建模型，然后重复O原子的吸附优化步骤。本节主要简单介绍一下如何快速搭建不同的吸附位点模型。

1 不同表面位点：

观察不同位点的结构特征，在面心立方金属的(111)表面上有4种不同的位点:

1) Top位，前面已经花了很多时间介绍了，本节就不再啰嗦
2) Bridge位：从名字大家就可以推断这个位点的吸附为两个Top原子的中间。
3) Fcc位和Hcp位：这两个吸附位点都是在三个原子的中心，为hollow位。

那么这两个位点怎么区分呢？

如果表面的Hollow位正下方（即第二层，图中黄色的）有原子，为Hcp位，在第三层的话（图中红色的），则这个hollow为Fcc位。
由于Fcc位点的原子在第三层，离表面最远（far），大家把英文单词的far 和 Fcc位的第一个字母关联起来即可。远的那个就是Fcc，近的是Hcp。
注意： Fcc 和 Hcp 分别是 face centered cubic 和 hexagonal close-packed 的缩写。

2 不同位点的吸附模型

怎么搭建Bridge和Fcc，以及Hcp位的吸附模型？

1) Bridge位:

这个位点在两个原子的中心。设两个原子的坐标分别为（x1，y1， z1）和（x2， y2， z2），它们中心得坐标为：（x1+x2）/ 2，（y1+y2）/ 2，（z1+z2）/ 2。在这里，z1和z2 是相等的。所以我们可以通过两个原子的坐标，计算一下中心的坐标，然后按照之前O原子搭建的步骤，添加O原子的Bridge位的坐标即可。
2) Fcc位和Hcp

这两个位点就更容易了，前面我们知道它们正下方分别为第三层和第二层的原子。那么我们在这两层随便取个原子，x和y 方向的坐标就确定了。
关于z方向的坐标，大家可以想象一下，如果O原子从Top位移到Bridge，Fcc和Hcp，由于和O成键的原子多了，也就是O和表面结合地更强了，那么它在z方向的坐标肯定比在Top位上要小。所以在这三个构型的z坐标，我们可以在Top位O的基础上稍微调小一点。那么小多少呢？一般来说0.1-0.2Å就可以了。

3 使用p4vasp搭结构

前面介绍的方法是直接修改POSCAR来搭建模型，对于O原子的吸附这种简单的模型，大家都还可以胜任，但如果体系复杂了，比如一个苯环，乙醇分子等等，这样我们的计算坐标的工作量就有点大了。此时，p4vasp的优势就开始慢慢显示出来了。

3.1 从Top到Bri

1）首先，我们要有一个O原子在Top位上的吸附结构
2）然后选择Top位的原子和相邻的那个原子，（p4vasp选中的结构看起来不是很清楚，见黄圈标出来的那两个）

3）选则 Edit –> Move Atoms 后会弹出这样的界面：

按照图中的1 2 3依次点击，你便会获得一个Vector。也就是从5到7号原子的一个向量。

但我们如果按照这个向量平移O原子的话，只会把它从5号原子移动到7号原子上，而不是Bridge的位点。

前面我们说了，Bridge的坐标就是这两个原子连线的中心。所以，我们把Vector减小一半就可以了。

修改完Vector之后，我们便可以移动O原子了。

这里注意

要把前面选中的5 和 7 号原子换成 O，要不然你移动的是 5和7 号原子而不是O。最后点move，效果如下：

这样我们便获得了Bridge的吸附结构。

前面提到，我们在搭建Bridge结构的时候，O原子的坐标可以比Top位的稍微低一些，这个操作怎么实现呢？
1）首先选中氧原子，
2）在Vector那一行修改移动的位置即可。如下图：

点击move ，保存POSCAR即可。

3.2 从Bridge 到 Hcp：

1）选中O原子和Hcp位下面的Cu原子，并获取Vector

2）修改Vector在z方向的大小。

点击 Move，效果如下：

保存成Hcp对应的POSCAR即可。

3.3 从Hcp到Fcc

1）选中O原子和Fcc的一个原子，获取Vector：

2）修改Vector在z方向的大小，由于Fcc和Hcp的吸附很接近，直接平移即可，也就是z方向大小为0

点击move，最终的Fcc结构，如下图，保存POSCAR。

注意：

保存结构的时候，可以自己建对应的文件夹，在文件夹里面将结构保存成POSCAR，也可以保存成类似POSCAR-Fcc， POSCAR-Hcp， POSCAR-Top这一种的名字，计算的时候，将它们重新命名为POSCAR就行了。

总结：

本节，大师兄主要给大家介绍了如何在表面上，使用p4vasp平移原子来搭建不同的吸附位点模型。使用p4vasp搭结构非常方便简单，也会加深你对模型结构的理解，另一方面，使用这一种操作的话准确性更高，搭建出来的模型比自己手动放原子要合理很多。所以强烈推荐大家使用p4vasp。

Ex59 吸附能的计算（四）

发表于 2018-10-07 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

丑话说在前头，本节内容有点多，大家慢慢消化。掌握本节的核心思想，对于后面提高你的工作效率，节省时间有着很大的帮助。前面一节，我们在扩展练习中提到了扩胞后，原子坐标的后面T T T 和 F F F 修改的问题。因为扩展后原子坐标是按照单胞为顺序排列的，如果想要将表面的原子快速固定或者放开，那么我们就需要知道表面原子所在的行数。

方法一：找规律+sed命令

上图中，p(2x2)的POSCAR命名为POSCAR-2x2，最表层原子所在的行数为：13，17， 21， 25，次表层的原子所在的行数为： 12， 16， 20， 24。仔细分析下，你会发现，表层原子的行数为： 9 + 4n ；次表层的行数为：8 + 4 n，其中n = 1-4，4表示在扩完后的slab模型里面有4个(1x1)的单元。有了这个关系，我们就可以通过sed命令任意固定这两层的原子了, 命令如下：

1	for i in {1..4}; do sed -i "$((8+4$i)), $((9+4$i))s/T T T/F F F/g" POSCAR

或者：

1	for i in {1..4}; do sed -i "$((8+4$i)), $((9+4$i))s/T/F/g" POSCAR

效果如下图：通过sed命令快速将表层原子批量固定住。

我们还可以只放开最表层的原子，使用sed命令快速将最表层的原子批量放开！

在前面结果的基础上可以使用命令：

1	for i in {1..4}; do sed -i "$((9+4*$i))s/ F F F/ T T T/g" POSCAR

或者：

1	for i in {1..4}; do sed -i "$((9+4*$i))s/F/T/g" POSCAR

Soooo Easy！！！再也不用挨行修改T T T 或者F F F 了。

方法二：分层处理+sed命令

什么，原来还有第二种方法？这种方法是本人常用的。前面我们提到，在p(2x2)的slab模型中，原子坐标是按照p(1x1)的单元排列的。既然我们想根据层数固定原子，那么我们可以先将原子按照层数来排列，这样就可以方便的选择同一层的原子了。

那么该怎么将它们按照层数排列呢？这时候你应该想到，不同层之间的原子，它们的坐标在z方向是不同的。有了这个依据，我们就可以根据z方向坐标的大小来排列原子了。

首先我们先用一个笨方法实现所预想的小目标。

第一步：将坐标复制到excel表格里面。本人用的是libre office(Linux下面的Microsoft)。不过功能大同小异，Windows用户按照做就可以了。

第二步：选中所复制的数据，然后根据图中的D列来排序。

排序完的结果, 如下图：

从图中，我们可以看出，原子根据在z方向上的坐标分成了4部分。

注意：这里我们所有原子都是Cu原子，因此可以全部选中然后直接排序。如果你的体系中有好几种不同的原子，你需要按照元素，逐一进行排序。比如体系中有C H O，我们先将C原子按照Z方向坐标排序，然后再排H原子的，最后再排O原子的。要保持和POSCAR前面元素顺序一致。

排序完成之后，将POSCAR中的坐标替换为Excel中的数据。

这样的话，我们的模型就搞定了。剩下的就是使用sed命令批量替换T T T和 F F F了

但是，师兄等等，上图看起来坐标很乱啊，结构能用吗？放心，绝对可以用。

如果不放心的话，可以用p4vasp打开一下，然后重新保存成POSCAR即可。最终我们的POSCAR如下图：

从图中可以看出来，次表层的原子在18-21行，最表层的原子在22-25行。所以可以这样使用sed命令来快速实现表层原子的固定与弛豫。

方法三：使用脚本分层

师兄，这个问题竟然还可以使用脚本来解决？

当然了，使用鼠标拖拖点点的事情，一般都可以用脚本来实现。脚本是用python写的，本人给它起了一个高大上的响亮名字：sortcar.py。（适用于python2.6以及以上的版本，低于2.6或者python3可能会出错。）

由于p(2x2)的slab坐标前面已经用过了，下面我们用p(3x3)的slab来演示一下脚本的使用。

新鲜出炉的，等待被sort的p(3x3)slab的POSCAR。被sort过的POSCAR如下图：

使用sortcar.py 将坐标按照z方向大小排列。输出文件为XXX_sorted，其中 XXX为我们想要排列的POSCAR或者CONTCAR。如果你的VASP坐标结构文件名字为：BigBro，那么被sort之后，输出文件就是: BigBro_sorted。

脚本这么神奇，内容如下：

#!/usr/bin/env python 
from collections import defaultdict
import numpy as np
import sys 

in_file = sys.argv[1]

### READ input file ###
def read_inputcar(in_file):
    f = open(in_file, 'r')
    lines = f.readlines()
    f.close()
    ele_name  = lines[5].strip().split()
    ele_num = [int(i) for i in lines[6].strip().split()]
    dict_contcar =  {ele_name[i]:ele_num[i] for i in range(0, len(ele_name))} 
    dict_contcar2 = defaultdict(list)
    for element in ele_name: 
        indice  = ele_name.index(element)
        n_start = sum(ele_num[i] for i in range(0, indice+1)) - dict_contcar.get(element) +1
        n_end = sum(ele_num[i] for i in range(0, indice+1)) +1
        dict_contcar2[element].append(range(n_start, n_end)) 
    return lines, ele_name, ele_num, dict_contcar2, dict_contcar

def get_elements(ele):
    lines, ele_name, ele_num, dict_contcar2, dict_contar = read_inputcar(in_file)
    coord_total = []
    my_list = []
    my_dict = {}
    for j in dict_contcar2.get(ele)[0]:
        coord_list = lines[j+8].strip().split()[0:3]
        tf_list = lines[j+8].strip().split()[3:]
        my_list.append(coord_list)
        dict_key = '-'.join(coord_list)
        my_dict[dict_key] = tf_list
       
    data = np.array(my_list)
    data=data[np.argsort(data[:,2])]

    for k in data:
         coord = '  '.join(k)
         tf = '  '.join(my_dict.get('-'.join(k)))
         coord_total.append(coord + '  ' + tf )
    return coord_total

## Generate the New POSCAR file

def Get_and_Save_lines(file_name, start_line, end_line):
    f = open(file_name)
    lines =  f.readlines()
    for line in lines[int(start_line):int(end_line)]:
        file_out.write(line.rstrip()+'\n')

out_name = in_file + '_sorted'
file_out = open(out_name, 'w')
Get_and_Save_lines(in_file, 0, 9)

ele_name = read_inputcar(in_file)[1]
dict_contcar = read_inputcar(in_file)[-1]

for i in ele_name:
    if dict_contcar.get(i) > 1 :
        file_out.write('\n'.join(get_elements(i)))
    else: 
        file_out.write('\n %s \n' %('  '.join(get_elements(i))))

注意事项：

1）这里我们用的是Cu(111)的例子，比较简答，脚本的威力展示不出来。
2）如果你的体系有很多不同的原子，脚本首先会将同一种元素原子的坐标排列，然后依次类推。
3）有些地方可能有些冗余。但本人一直在用，效果还不错。
4）本人建议，只针对表面使用这个脚本。也就是吸附前的slab模型。如果你放了吸附物种，不建议使用，因为吸附物种间不同原子之间的连接顺序可能会被打乱。
5）在放吸附物种之前，最好将表面先排序完，怎么做，前面三种方法任选。
6）脚本的话已经上传，大家可以通过链接下载：https://pan.baidu.com/s/1X5xLRsvRmFfNE8IlKwvbFQ）

7）如果使用脚本出错的话，可能是因为你的python版本太低。这个脚本只适用于python2.6及以上的。Python3可能也不适用。出现错误，直接放弃，掌握了这个方面的思想就足够了。

总结：

本节讨论的内容与科研关系说大不大，说小也不小，主要是在解决工作效率和时间上。想象一下，如果你的slab中有很多原子（金属氧化物，硫化物等），而且还不是同一种的，你用前两种方法的时候，就得十分小心了，此时脚本的作用就显示出来了。

Ex58 吸附能的计算（三）

发表于 2018-10-05 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

前面一节我们将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 （图中箭头所指的部分）。这样会很麻烦，还浪费很多时间。

怎么解决这个问题，快速将表层，次表层原子固定或者放开呢 ? 看下一节的内容。

Ex57 吸附能的计算（二）

发表于 2018-10-04 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

$\require{mediawiki-texvc}$

本节我们简单讨论一下覆盖度对吸附能的影响。首先回顾一下，前面我们讲到的内容。

吸附能为初态和末态的能量差：
- 初态为纯净的slab 和 O$_2$分子或者O原子
- 末态为O吸附的slab
O$_2$分子优化计算的时候，要注意:
- 从数据库获取O$_2$分子的键长，作为初始值；
- ISPIN = 2
- ISMEAR = 0; SIMGA = 0.05
- IBRION= 2; POTIM = 0.1
- 用gamma点。
O 原子优化计算的时候，注意的计算细节与前面相同。额外注意的是：
- O 原子再怎么优化还是一个O原子，所以上一行表述是不恰当的，应该是O原子的单点计算；
- Box的尺寸为：$13\times14\times15~\AA^3$ 或者$13.1\times13.2\times13.3 \AA^3$ 但绝对不能为$13\times13\times13$这种正方体的
怎么优化Bulk的结构，扫晶格拟合和ISIF = 3直接优化（ENCUT要大，不妨取个700~eV）
slab怎么优化，Selective 和坐标后面 T F 的关系
O吸附模型怎么搭建，如何得到合理的初始模型（有实验值优先使用，没有的话用粗精度的初算一下）。
此外，还要学会怎么判断计算结束，不同类型任务结束的特征是什么？

打个比方

回到上一节我们留的问题，为什么O在p-(1x1)-Cu(111)表面上的吸附能是正的呢？在解释这个问题先，我们先看下面一组照片：

第一张里面，地铁车厢里面是空的，这时候你上去，座位随便挑。

第二张里面，有些人了已经占了，但还是有座位的，你的选择并没有那么多了。

第三张里面，人太多，能不能出来都是个问题，有没有座位想都不敢想了。

那么回到我们的这个问题，我们把O的吸附结构在xy平面上重复一下：

你会发现，表面上密密麻麻全是O原子（绿色的为Cu原子）！空间这么小，如果你是O原子，你愿意在这上面挤么？你是选择空旷的车厢还是人满为患的呢？如果我是氧原子，我肯定不愿意。所以，如果非得在这个表面上挤的话。我们就需要一些外力的作用，比如下图中的这俩年轻力壮的小伙子协警。

这俩人使的劲就是O原子的吸附能。正值的意思表明我们需要额外的力才能将其吸附在表面上。表面越挤，正值越大。

如果你还是不理解的话，我们看另外一个例子：

早上，开学了。老师有1个苹果，需要将这个苹果分给教室里面的学生。而学生还没有来齐，当你第一个到达教室时，发现整个苹果都是你的。开心程度 100%。正要吃的时候，小明进来了，老师说你俩要分一下，你只能得到一半了，开心程度 50%; 你俩正要吃的时候，小红和小白也进来了，老师说，打住，你们四个分，于是每人只能拿1/4，你的开心程度掉到 25%了。

表面就是这个苹果，O原子就是学生，总共就这一个苹果，学生越多，平均分到的就越少。学生之间可能还会因为分的不均匀而打架（相互排斥），导致开心程度越来越低。当开心程度成负的时候，老师就需要施展神通费心思（额外的功）逗你们开心。

看完这两个例子，你也许就明白了，为什么这个时候O原子的吸附能这么大了。太他妈挤了或者人太多，分不到苹果吃了，心里肯定不爽！下面，我们把车厢清空，学生数目控制一下，看看心情能不能好点。

搭建p-(2x2)-Cu(111)表面

1 用p4vasp打开p-(1x1)-Cu(111)的CONTCAR

2 选择Edit —> Multiply cell

3 弹出的窗口，如图填入扩展的倍数：

这代表我们在x和y方向上分别扩展到原来的2倍。点击Multiply，效果如下图：

4 关闭Multiply这个窗口，p4vasp主界面左上角选择：File —> Save system as

5 在弹出的窗口，选择目录，保存成POSCAR，点击OK即可。

6 这样我们就有了一个p-(2x2)的slab模型了。

搭建O在p-(2x2)-Cu(111)上的吸附结构

1 我们可以用之前计算的结果搭建这个初始结构,O吸附优化的CONTCAR：

O原子在$z$方向的坐标为8.0928$\AA$。先记住这个数字。

p-(2x2)的Slab坐标（在p-(1x1)基础上扩展的）:

在这个坐标上面修改，搭建结构：

1）添加O原子元素

2）添加O原子数目

3）添加O原子坐标

4）弄完后，保存成POSCAR

到现在，p-(2x2)的slab以及O原子的吸附模型都有了，准备INCAR,KPOINTS,POTCAR以及任务脚本，提交这两个结构的优化计算。等待结束。

计算吸附能：

公式1：

$E_{ads1}(O) = E_{slab+O} – E_{slab} – E_{O_2}/2$

公式2：

$E_{ads2}(O) = E_{slab+O} – E_{slab} – E_{O}$

带入数据：

$E_{ads1}(O) = -60.61682891~eV – (-55.96049190~eV )- (-9.85498627~eV)/2 = 0.271~eV$ $E_{ads2}(O) = -60.61682891~eV – (-55.96049190~eV )- (-1.89237882~eV)= -2.764~eV$

发现吸附能$E_{ads1}(O)$从正的1.216降低到0.271 eV了。看来空间对O原子的吸附影响很大。

地铁里面空间稍微大了些，现在不需要年轻小伙子暴力推搡了，大妈指挥下就可以了。

分的苹果稍微多了些，虽然还是有些不情愿，稍微诱导一下，就开心了。

前面的解释只是用现实的例子瞎扯一顿。但发表文章的时候不能这样说啊，我们就需要一些科学和专业的解释。而在实际表面上，覆盖度影响吸附能的因素有很多。体相中电子的分布转移，与吸附物种的成键，吸附物种之间的排斥和吸引等等。如果你想了解的更多，更专业！阅读下面这篇文献：

谷歌上可以下载免费的，有权限的也可以直接下载。类似这样已经注明的文献，请不要花时间找我要，自己去下载。

（每设一道坎，就会把很多懒家伙挡住！跨过一道又一道坎，你就把别人远远甩在后面了。所以，请自己主动起来！）

扩展训练：

1 重复p-(2x2)-Cu(111)的练习，心里面默念一遍INCAR, KPOINTS, POTCAR, POSCAR。看看都有些什么需要注意的。

2 阅读推荐的文献，学习表面物种之间的相互作用，覆盖度是怎么影响吸附能的。有哪几方面？

3 本文提到了一个专业名词：覆盖度，查找相关的文献，了解这个名词是怎么回事，文献报道的覆盖度都是怎么算出来的。

总结：

本节，大师兄带你做了趟地铁，当了回幼儿园的小朋友。来体验下人多的感觉。覆盖度高的时候吸附能很小或者是正值，说明物种在表面不稳定，需要额外的力来促使它们老老实实待在表面上。表面就那么大点的地方，空间有限，p-(1x1)上吸附的时候，表面太挤了，O原子不愿意待在上面。而表面扩成p-(2x2)后，空间明显大了很多，O原子虽然不情愿，但明显不像在(1x1)上那么不爽了。从吸附能就可以看出来，降低了1~eV左右。从另一个角度来看，这明显就是一场O原子争夺表面资源的战争。

Ex56 吸附能的计算（一）

发表于 2018-10-03 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

我们已经优化完成了O原子在Cu(111)表面上的吸附。总结并数一下，到目前为止我们从前面的计算中可以获取哪些能量。

Cu(111) slab 的能量
slab上有O吸附时，体系优化完的能量

再往前想一下，我们也计算了

O$_2$分子的能量
以及O原子在气相中的能量。

有了这几个能量，我们就可以计算O在p(1x1)的Cu(111)表面上的吸附能了。这也是本节的主要内容。

1 什么是吸附能？

从字面上不难理解，就是分子或者原子从气相中吸附到表面上所释放的能量。计算吸附能的时候，需要注意的有两点：起始和终态。

1） 终态：这里说的终态就是O原子在表面吸附，并且优化完的构型。我们标记为 slab+O，它的能量为 E(slab+O)。

2） 初态：我们计算O的吸附能，起始状态是O$_2$ 和纯净的slab表面。因为在实际反应中，O$_2$在气相中解离成2个O原子，然后再吸附，这种可能性微乎其微。前面我们讲了O$_2$分子的优化，这里我们直接把O$_2$分子的能量拿来用。起始的两个结构能量分别标记为：E(O$_2$) 和 E(slab)。

由于在终态的吸附结构上面，我们只有1个氧原子，所以初始态我们要用O$_2$能量的一半: E(O$_2$)/2。

3）注意：文献报道里面，也有很多人用O原子在气相中的能量来计算O的吸附能。这样做的话，其物理意义为：O$_2$需要先解离成O原子，然后O原子再吸附。此时，初始状态为分解的O原子，因此O$_2$分子解离的能量没有考虑在内。由于O$_2$解离是吸热反应，忽略掉解离能，会导致O的吸附能很强。单个O原子的能量标记为：E(O)。

4）所以，大家在计算的时候，一定要把自己的计算公式标出来，此外，大家看到直接用O原子能量计算出来的吸附能，不要用O$_2$作为参考能量的结果进行对比。

5）这两个计算方法，哪个对，哪个错呢？答：都是对的。因为O$_2$的解离能是一个常数，加上去（用O$_2$能量作为起始）或者忽略掉（用O原子能量作为起始）得到的结果之间的区别无非也是这个常数。看下面的计算。

2 计算吸附能（单个O原子）

公式1:

$E_{ads1}(O) = E_{slab+O} – E_{slab} – E_{O_2}/2$

公式2:

$E_{ads2}(O) = E_{slab+O} – E_{slab} – E_{O}$

$E_{slab+O}$ ：O 在Slab上，优化完结构的能量；

$E_{slab}$: 优化的slab的能量；

$E_{O_2}$: O$_2$ 分子在气相中的能量；

$E_{O}$: O原子在气相中的能量；

带入数据后：

$E_{ads1}(O) = -17.68238420~eV – (-13.97082922~eV) – (-9.85498627~eV) = 1.216~eV$ $E_{ads2}(O) = -17.68238420~eV – (-13.97082922~eV) – (-1.89237882~eV) = -1.819~eV$

前面我们说了，如果直接用O原子在气相中的能量，则忽略了O$_2$分子解离的能量（或者是O原子结合能）。

我们知道O$_2$ 分子解离是一个吸热反应，那么，解离能可以通过以下两个方式来获取：

1) 两个吸附能量相减来得到：

$E_{ads1}(O) - E_{ads2}(O) = 1.216~eV – (-1.819~eV) = 3.035~eV$

2) 直接通过O$_2$分子和O原子的能量得到。单个O原子的结合能为：

$E_{binding} = [ E_{O_2} – E_O * 2] /2，$

代入数据后得到：

$E_{binding} =[ -9.85498627~eV – (-1.89237882~eV *2) ] / 2 = -3.035~eV$

3 思考问题：

1）前面我们计算的$E_{ads1}(O)$ = 1.216 eV 是一个正值，为什么会这样？算错了吗？还是有其他物理意义？

2）如何计算H，CO，CH$_4$, CH$_3$OH 等的吸附能，把计算的公式列出来。每个能量计算的细节，注意部分有哪些？

3）如何计算脱附能？

4 扩展练习：

1）推荐阅读下面这两篇文章（本人的，当然喜欢大家多多引用了。）讲的是关于甲醇，乙醇，乙二醇，以及甘油在Cu, Ru, Pd以及Pt的(111)上面的解离。Ex系列的后续文章，表面反应的优化，过渡态的计算等，基本都以这两篇文章的内容为主。希望大家能够好好学习，并多多引用。此外，计算的很多细节都在支持信息里面，大家不要错过。一时半会看不完没有关系，下载后结合本书慢慢看。

链接：https://pubs.acs.org/doi/abs/10.1021/cs501698w

链接：https://www.nature.com/articles/s41467-018-02884-y

2）回顾并重复本节所涉及到的计算，以及注意细节部分。

A. 单个O原子的能量；

B. O$_2$分子优化；

C. Cu bulk的优化；

D. p(1x1)-Cu(111)表面的搭建及优化；

E. O原子在p(1x1)-Cu(111)的top位吸附的模型搭建及优化；

G. 不满半年的新手，建议从Ex0 从头学习。

5 总结：

本节我们主要介绍了O原子在p(1x1)-Cu(111)表面top位上的吸附能的计算公式。顺便带大家回顾一下整个计算流程。如果你从头跟着学过来的话，每个流程中的细节，需要注意的部分肯定都已经掌握了，后面举一反三，就慢慢渐入佳境了。也就是所谓的开始入門了。

Ex55 简单粗暴地获取初始构型（二）

发表于 2018-10-02 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

前面一节，我们通过一个简单的方法获取了一个Cu-O键的键长。因为我们固定了Cu的坐标，当计算完成之后，O原子初始构型的坐标可以直接将计算后的结果复制过来。此时需要注意Direct和Cartesian的坐标问题。

1 搭建合理的O吸附模型

1）取Cu(111)表面模型的POSCAR (Cartesian坐标)

2) 复制快速获取初始构型计算(Cu-O dimer)的结果

上图中的最后一行。

3) 构建O的吸附模型：

当然了,我们也可以根据键长或者O在z方向的坐标直接在Cu(111)表面基础上搭建。

这个简单粗暴的方法到此就讲解完了。有以下几个需要注意的地方：

A）这里我们用的是Cu(111)的 p(1x1) slab模型,表面只有一个原子。一般来说，大家计算吸附的时候，模型都比这个大，我们可以取一层原子。记住，要固定住这层原子。

B）这个方法基于的是气相中的计算，因此偶极矫正不要加！因为加上去之后，收敛会变得很困难。这个大家可以用自己测试测试。

下面我们讲一个基于slab模型的方法，虽然比这个计算量大，但也异常的快，相对于后面的计算，花费的时间也可以忽略不计。

2 快速获取初始构型方法（二）

1）在这个方法中，第一步我们先取slab模型的结构，这里就不讲了。

2）搭建一个初始的吸附模型。此时我们要根据原子半径，大体确认一下键长的合理范围，前面也讲过了，就不再细说了。由于前面一节我们已经有了一个O原子的坐标数据，可以直接拿过来。如下：

注意：实际操作的时候，这个方法可以使用前面一节的方法的出来的结果；也可以直接设置初始值，这两个方法之间联系不太大。

3）固定slab模型中的所有原子！所有原子！所有原子！

重要的事情说三遍，是把slab模型中的原子全部固定!只放开吸附的分子。这里大师兄在vim里面操作，使用了下面的命令：:12,13s/T/F/g 意思是把12和13行中的T全部替换成F。当然也可以使用sed进行操作，相信大家经过这么长时间的练习，已经掌握基本的Vim和sed操作，就不浪费时间了。效果如下：

4）设置INCAR，这里是slab模型，把之前Cu(111)表面优化计算的INCAR直接拿过来了。

EDIFF、EDIFFG 可以适当放宽，下图大师兄懒得改了。

5）生成对应的POTCAR，K点只用gamma点！K点只用gamma点！K点只用gamma点！

这里需要注意的有两点：

A）计算前一定要检查POTCAR和POSCAR是否对应，养成这个好习惯

B）一定要用gamma点。因为我们知道，K点越多，计算量越大。

6）提交任务等待结束。

O原子的坐标，从最初的7.950优化到了8.009。

到此，另外一种快速获取稳定构型的方法，也介绍完了。希望大家可以举一反三，运用到自己实际的课题计算中。

3 正式优化O的吸附

前面我们已经获取的初始的一个合理的构型，下面就开始正常计算了。想一想，在上一步的基础上，正式计算的时候，我们有哪几个需要修改的部分。

1）结构

这个是必须要想到的。怎么做呢？直接将CONTCAR 复制成 POSCAR即可。

PS：怎么完成复制这一步还不会的话，请从头开始学。

2）复制完了并不代表完事，前面一步我们还固定了表层的原子，此时我们需要放开。使用sed命令：

1	sed -i '12,13s/F/T/g' POSCAR

3）KPOINTS：使用前面Cu slab优化时用的 13x13x1

4） INCAR中，如果前面一步的计算精度较低，这时，我们就要提高一下了。

5）POTCAR 不变

6）提交任务的脚本：由于前面的计算量很小，一般可以使用小队列进行计算，如果计算量增大了，这时我们也需要稍微修改下脚本。

7) 提交任务，等待结束，查看结果。

O原子的坐标从8.009 优化到了8.093。键长从1.75 Å 优化到 1.80 Å，变化很小，这说明我们的初始构型已经很接近优化结果了。

4 扩展练习：

1) 本节的实例中，使用了 && 这个连接前后的命令，大家学习下是怎么回事

2) 重复本节练习，掌握其中的关键点。

5 小结

一个合理的初始构型可以极大地降低我们的工作量，这一点是大家务必要记在心里的。而获取合理的初始构型，一方面需要我们用化学的知识去搭建结构，另一方面，也需要配合一些简单的低精度的计算去完成。希望大家可以掌握这两节的精髓，运用到自己的课题中，提高自己的计算效率。

Ex54 简单粗暴地获取初始构型（一）

发表于 2018-10-01 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

$\require{mediawiki-texvc}$

Ex54-Ex56主要介绍一下如何计算：单个氧原子在p(1x1)-Cu(111)表面top位上的吸附。在实际的计算过程中，一个好的初始结构会极大地加快并节约你的计算时间，这不仅仅体现在服务器运行的时间上，也会避免很多错误的结果，因为这些错误结果的纠正耗时费力，通常折磨得新手们心力交瘁。可以回顾下前面我们O$_2$ 分子的优化过程：不合理的初始结构导致的错误结果以及计算时间的增加。本节主要介绍一个快速获取优化初始构型的方法。该方法简单粗暴，可以让你在极短的时间内快速获得一个化学的feeling。

1 Top site 的构型

在Cu(111)表面上，top位指的是直接吸附在Cu原子上方。其他的位点是什么，先不要着急，后面会慢慢学。如果一个O原子吸附在Cu原子上方，在空间坐标上，大家很容易想到，O和Cu在xy方向的坐标应该是差不多的，z方向O的坐标为：Cu的z坐标 + Cu-O键的键长。如下图：

2 确定O原子在z方向上的坐标

知道了前面的基本原理，O原子的坐标我们就可以通过Cu-O的键长来初步获得了。那么Cu-O的键长怎么获得呢？这里列举了一下几个常用的方法。

1）查数据库：

2）查文献：自己去查，不要在QQ群里让别人给你发文献。

这两条主要考验的是大家查询资料的能力，这里暂且不详细介绍了。有兴趣的可以加入大师兄文献互助超级群，跟众多文献检索大牛学习。(群号:157099073)

3）自己估算

估算的话，可以根据Cu和O原子的半径：Cu和O成键，键长肯定要小于两者的半径之和。搜一下维基百科，Cu的半径为1.35 Å，O的为0.60 Å。(维基百科，搜索Atomic Radius就可以得到下面这个表)

从上面的数据我们可以得出：Cu-O键要小于1.95 Å。在吸附构型搭建的时候设置比1.95 Å小点的值作为初始，进行优化。一般来说，原子半径之和减去0.1-0.3 Å都是可以的。但是不能太小，否则第一步优化的时候排斥力太大，会导致计算出错。

4）自己初算：

初算就是初步采用一个小的模型，简单优化一下，得到一个合理的键长数值。小的模型主要指2个方面：

结构简单；
计算参数简单。

下面我们主要在估算的基础上，介绍一个初算的方法：直接优化一下气相中Cu-O双原子分子的结构。

这个结构不难搭建，将前面练习中O$_2$分子中的一个O原子换成Cu进行优化即可。

为了加深大家对搭建模型的印象，我们从p(1x1)-Cu(111)表面的结构出发，然后一步一步搭建CuO的气相结构模型，并计算。

3 搭建初算的模型(Cu-O双原子结构)

1）p(1x1)-Cu(111)表面的结构

将前面几节计算的一个真空层为15 $\AA$的例子直接拿过来。表层原子在z方向的坐标为6.259 (第十二行)。

2）修改格子大小(3-5行)

这里修改的很随意，第三行中直接把2.571改成了12.571；

第四行中把2.227改成了12.227。

第五行中将21.299改成了12.299

效果如下:

注意：

A）当然也可以用其他大小的格子；(例如：8x8x8 $\AA^3$)

B）格子大小直接影响计算量和时间。(回顾下前面所学)

3）保留表层的Cu原子，删除其他的Cu原子。

注意：

A）我们这里将表层以下的三个原子删掉了，只保留z坐标为6.259的这个原子；

B）第7行中原子数目也要相应的改变。从4 改成1。

C）效果如下图：

4）加入O原子

注意：

A）第6行中，在Cu的后面加上了一个O，不是数字0；

B）Cu 和 O 之间用1个或者几个空格隔开，不要用Tab；

C）第7行中，记得加上O原子的数目

D）最后一行添加O原子的坐标，这里我们直接把Cu原子的复制过来了。

由于两个原子坐标一样，Cu和O堆在一起了。

5）修改O原子坐标

注意:

A）这里我们把Cu固定住了;为了方便下一步的计算。

B）Cu-O之间的键长，设置的为1.8 Å；

C）O的z坐标为：6.259 + 1.8 = 8.059。

6）到这里，Cu-O双原子分子的气相结构就搭建完毕了，保存成POSCAR即可。

上面的效果图是本人每一步打开展示给大家的。实际操作中，完全没有必要。我们需要学习的是：

A）怎么将脑子中的模型转化为VASP的POSCAR文件；

B）格子大小怎么修改；

C）怎么添加原子，添加或者删除原子后，原子数目怎么弄；

D）怎么添加原子坐标。

E）怎么通过改坐标修改原子位置。

4 INCAR检查

注意:

1）算的是气相中的：ISMEAR = 0；SIGMA = 0.05 (本书前面就讲过了)

2）IBRION = 2；POTIM = 0.1；NSW =100 是优化的参数

3）EDIFF和EDIFFG是电子步和离子步的收敛标准。

师兄，磁性呢？对称性呢？

因为本例子是一个初算，几步算完，看下键长。这个任务的使命就完成了。

很多细节的东西可以暂时不用考虑。这里EDIFFG用的也有点小，本人忘了修改了。大家可以设置为 -0.05或者直接使用能量作为标准，这些都可以的。

虽然本节我们很多地方都没有考虑，这是由于任务的性质所决定的。我们对其定位就是瞎算下，得到一个初始的构型。这也是一个课题中为数不多的，可以为所欲为的计算了。但你正儿八经开始算的时候，各种细节问题都要考虑进去。而这一步的计算，也可以作为一个缓冲期，让你有充足的时间去思考一下正式计算时其他需要的注意事项。

5 POTCAR、KPOINTS检查

A) 根据POSCAR获取对应的POTCAR：见附录自动生成POTCAR的脚本。

B) 使用gamma点计算。cat POTCAR

注意：

在提交一个任务前，一定要将INCAR、KPOINTS、POSCAR、POTCAR以及脚本在心里默念一遍。然后对应的检查一下，是否有些遗漏的地方。否则等算完了，发现错了，又得浪费很多时间重新计算了。

6 检查结果

1）上个厕所的功夫，任务就结束了，共算了6步。(其实是切换电脑系统的功夫)

2）Cu—O 的键长：7.95-6.26 = 1.69 Å

注意：

此时的键长只作为下一步的初始猜测，如果和文献去比有点差异，也不要太较真。

7 总结:

本节，你会学到：

1）如何通过原子半径，估算一个初始的键长值。

2）如果搭建一个简单的模型，初算一个键长值。

3）复习气相中的优化。

4）复习下计算前准备工作：INCAR、POSCAR、POTCAR、KPONTS的检查。

如果本节的内容，你有不理解的地方；想获取脚本等信息；请查阅前面讲的以及附录中的内容。搭建几何模型，无非就是按照软件的格式，修改原子坐标罢了。这里我们强调的是几何模型。但，几何模型所具有或者要表达的物理、化学意义，这才是最关键的。

Ex53 表面吸附以及模型的搭建思想

发表于 2018-09-30 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

吸附这个词，对于做表面化学的人来说，是一个普通的不能再普通的名词了。从今天这一节开始，我们将逐步进入表面化学相关计算领域。

1) 我们先看一下维基百科的解释:

a）吸附是指某种气体，液体或者被溶解的固体的原子，离子或者分子附着在某表面上。这一过程使得表面上产生由吸附物构成的膜。吸附不同于吸收，吸收是指作为吸附物的液体浸入或者溶解于另一液体或固体中的过程。吸附仅限于固体表面，而吸收同时作用于表面和内部。

b) 吸附也属于一种传质过程，物质内部的分子和周围分子有互相吸引的引力，但物质表面的分子，其中相对物质外部的作用力没有充分发挥，所以液体或固体物质的表面可以吸附其他的液体或气体，尤其是表面面积很大的情况下，这种吸附力能产生很大的作用，所以工业上经常利用大面积的物质进行吸附，如活性炭、水膜等。吸附过程有两种情况:

一种是物理吸附，在吸附过程中物质不改变原来的性质，因此吸附能小，被吸附的物质很容易再脱离，如用活性炭吸附气体，只要升高温度，就可以使被吸附的气体逐出活性炭表面。

另一种是化学吸附，在吸附过程中不仅有引力，还运用化学键的力，因此吸附能较大，要逐出被吸附的物质需要较高的温度，而且被吸附的物质即使被逐出，也已经产生了化学变化，不再是原来的物质了，一般催化剂都是以这种吸附方式起作用。

2) IPUAC GoldBook中的定义:

Adsorption: An increase in the concentration of a dissolved substance at the interface of a condensed and a liquid phase due to the operation of surface forces. Adsorption can also occur at the interface of a condensed and a gaseous phase.

https://goldbook.iupac.org/html/A/A00155.html

3) 参考书

吸附在物理化学中也是一个非常重要的概念，只要学过物理化学，肯定不会错过这一部分的内容。主要列举有关吸附的几本参考书，具体内容本节则不再详细介绍，这也超出了本书所包含的范围。

a) Atkins’ Physical Chemistry

(链接: https://pan.baidu.com/s/1hsJ95sg 密码:nfxw)

b) Concepts of Modern Catalysis and Kinetic (By I. Chorkendorff）

链接: https://pan.baidu.com/s/1hrA5LHi 密码:aqvq)

c) 南大物理化学(傅献彩主编)

吸附在催化研究中的地位

催化剂发挥催化作用的前提是:反应物需要和催化剂有接触（也就是反应物吸附到催化剂表面上并且得到活化）。写这句话的意思就是:研究一个催化反应，第一步我们需要搞懂的就是反应物在表面的吸附情况。个人计算的心得就是：反应的初始态，过渡态以及末态，它们都是表面的吸附物种。这句话的意思是：在我们的计算中，吸附贯穿着我们对催化现象的整个研究过程中。吸附不仅仅是把反应物放到表面上就完事了，更与后面的反应息息相关。切记不要在思考的时候把后面我们要学习的过渡态和吸附隔离开。

吸附模型的搭建

前面我们已经讲过slab模型了，搭建吸附模型的本质就是在slab模型基础上，再添加吸附分子的xyz坐标。搭建的模型无非就是一个xyz的坐标文件。需要注意的有3个方面:

1) 熟悉VASP的POSCAR格式: 每一行的内容所代表的意思以及它们之间的相互联系:

元素行，原子数目以及原子的坐标之间的关系；
Selective Dynamic 行与 T T T， F F F 等原子固定的关系；
Cartesian 坐标和 Direct 坐标的区别等等；
POSCAR与KPOINTS的关系。

熟悉POSCAR或者CONTCAR的格式可以使你在搭建模型过程中游刃有余。

2) Slab模型在Z方向上:

slab的厚度
真空层的厚度，

这两点一方面决定了我们计算量的大小。另一方面，对于不同的体系，我们需要不同厚度的slab模型来保证计算的准确性。例如：对于金属体系来说，越开放的表面往往需要更多的层数。所以，在准确性和计算量上，我们要合理地权衡和取舍。可以通过一些基本的测试工作，以及参考他人发表的计算参数来确定。

3) Slab模型在XY方向上:

表面的大小：这个主要影响覆盖度，计算的工作量；
表面的吸附位点：比如面心立方金属的111表面，有top(t)，bridge(b)，fcc(f) 和hcp(h)这几个位置，在搭建模型的时候要考虑吸附物种分别在这几个位置上的情况；
吸附物种与表面的结合情况：不同的分子构型?用什么原子?哪一部位和表面接触?初始猜测的键长多少? 等等

这些对于你搭建合理的模型至关重要。一个合理的模型可以极大地降低计算工作量，提高你的计算成功率。当然，这需要一个扎实和良好的化学基本功来作为保障。

小结

本节简单介绍了一下吸附的概念，吸附模型搭建的一些基本的思路，大家可以自己主动思考一下:本节提到的相关内容之间的相互关系以及具体的注意事项。在后面具体的实例操作计算中，我们会逐步慢慢展开详细介绍。

Ex52 影响功函数计算的因素（二）

发表于 2018-09-06 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

前面功函数的计算流程( VASP 计算细节以及结果分析)，我们已经掌握了。这一节介绍一下影响功函数的主要因素：真空层的厚度。首先个人经验：计算功函数的时候，真空层厚度一定要足够。

这里有两点：

多厚才算足够?具体要根据自己的体系测试下，一般来说 20 Å 左右足以。
是不是越厚越好? 越厚的话，计算出来的结果差别不大，但会增加计算量。前面我们讲过影响计算量的因素。

我们以 Cu(111) 表面的例子，简单介绍下真空层厚度对功函数的影响。

1 不同真空层厚度的测试(上一节内容)

先准备一个真空层为 10 Å 的计算文件夹： POSCAR 如上图；
将前面两节的 INCAR,KPOINTS, POTCAR, 以及提交任务的脚本复制过来；
批量生成测试文件，每个文件夹以真空层的厚度命名。命令如下：
1
for i in $(seq 12 2 36); do cp 10 $i; sed -i "5s/16/$((6+$i))/g"$i/POSCAR ; done
批量提交任务。

注意：

这样的批量操作本书刚刚开始的时候就已经讲到了，这里不要再问我。
如果不会用这个办法批量制作输入文件，想其他的办法批量操作，条条大路通罗马，不要死磕在我的这个命令上。

2 结果处理

前面的测试计算中，共有 14 个任务。在 Ex50 这一节中vtotav-v5.2f, vtotav.py 或者 p4vasp 这三种可视化的方法，我们查找真空能级的数值，都必须读取数据文件，比如 VLINE, LOCPOT_Z 或者 p4vasp 导出来的 *dat 文件。所以，我们要分别在这 14 个任务的文件夹中重复 Ex50 的操作14次。手动挨个弄可以完全实现，但本人不想这么做，因为这样弄除了让你的操作更熟练些外，并不会有其他的收获。那就写个脚本批量处理吧。

首先我们分析下上图：我们获取的真空能级在图中直线部分中间的一个点或者一个小区间的平均值。大约是真空层一半的高度。比如，我们的表层原子 z 方向坐标为 6，真空层厚度为 16，那么我们可以取 6 + (16/2)= 14 时对于的点来得到真空能级。也可以通过 [14-1, 14+1] 这个区间中所有的真空能级求平均数来获得。下面我们看一下脚本的主要内容：

注意部分：

使用这个脚本的时候，需要注意的内容。本人写了一堆废话在脚本里面，主要是提醒大家使用脚本的时候，需要注意的一些事项。

#!/usr/bin/env python

# Welcome to visit our website: www.bigbrosci.com to get more useful information.
 
'''
Written By Qiang 
This script has two functions:
 1) workfunction Visualization from LOCPOT_Z file
 2) Roughly estimate the Vaccum Energy

Please Contact lqcata@gmail.com if you have questions.

Note1: Read POSCAR with Cartesian coordinations! 
  This means 
  1) you have to convet the direct to cartesian firstly 
  2) otherwise you will get an error.

Note2: Read LOCPOT_Z to plot the figure and calculate the vaccum energy
  This means 
  1)1you have to run command:  work.py LOCPOT z firstly to get the LOCPOT_Z file before use this script

Note3: the idea to calculate the vaccum energy is following:
  1) get the middle z value in the plateau: 
     middle z value = ( coordination of the highest atom + length of the slab in z direction) / 2 
  2) selecte the area to calculate the vaccum energy: 
     from  (middle z  value - 1) to (middle z value + 1)  unit is in angstrom 
  3) do the average of all points y direction.

Note 4: this is only a rough estimation but useful.

Note 5: Check the figure firstly and then use the numbers calculated from this script. 
  1) If the midlle z value is far from the plateau in the figure, you have to calculate the energy by hand.

'''

脚本是谁写的，主要功能是什么；
脚本读取 POSCAR, POSCAR 要为 Cartesian 坐标；
脚本默认读取 LOCPOT_Z 文件，如果没有该文件，你需要做的一些事情(生成以及如何将 VLINE 文件转换为 LOCPOT_Z 文件)
该脚本是如何获取真空能级的；
一定要查看功函数可视化的结果，用来判断我们的这个方法是否试用。

写这么多废话，原因主要有 2 个：

提醒大家在使用别人脚本的时候，一些需要注意的事项你要清楚。如果不清楚，可以咨询下脚本的作者；
如果你写了一个脚本放到网上供大家使用，请花些时间写清楚脚本的运行原理，以及尽可能详细的注意事项。

脚本正文：

import matplotlib.pyplot as plt
import numpy as np 
x = []
y = []
dic = {}
with open("LOCPOT_Z", mode='r') as f:
    first_line = f.readline()
    name_x = first_line.split()[1]
    name_y = first_line.split()[2]
    for line in f:
        xy=line.rstrip().split()
        x.append(float(xy[0]))
        y.append(float(xy[1]))
        dic.update({xy[0]:xy[1]})
plt.plot(x,y) 
plt.xlabel(name_x) 
plt.ylabel(name_y) 
plt.savefig('workfunction' + '.pdf', dpi=400)
plt.show()

#%%%%%%%%%%%%%%%%%%%%%%#
# Get the Vaccum Energy
#%%%%%%%%%%%%%%%%%%%%%%#

# Get the total line numbers of POSCAR 
num_lines = sum(1 for line in open('POSCAR'))

# Read POSCAR 
pos = open('POSCAR', mode = 'r')
line = pos.readlines()

# Get  the  slab length in z direction 
vaccum = float(line[4].split()[2])

# Get all atoms' coordination in z direction and store them in the list
z_list = []
for i in range(9, num_lines): 
    z_list.append(float(line[i].split()[2]))

#  max(z_list) is highest atoms' coordination in z direction
# Get the vaccum lenth: 
l_vaccum = vaccum - max(z_list)
print 'The Vaccum in this calculation is:\t\t %s'  %(l_vaccum)

# Choose the middle z value in the workfuntion.pdf 
num_middle = (max(z_list) + vaccum) / 2 
#print num_middle

middle_y = []

for i in dic.keys():
    i = float(i)
# Select the date area within 1 angstrom from the middle point:
    if i > num_middle -1 and i <  num_middle + 1:
        middle_y.append(float(dic.get(str(i))))
# Get the average value in the selected area
print 'The Vaccum Energy in this calculation is:\t %s'  %(np.mean(middle_y))

pos.close()

plt.show() 之前为我们 Ex50 这一节可视化的脚本： 44 和 53 行中，额外将 LOCPOT_Z 中两列的数据保存到一个字典里面，便于后面根据 z 方向的坐标查找真空能级。

获取 POSCAR 文件总的行数

1	num_lines = sum(1 for line in open('POSCAR'))

读取 POSCAR 文件

1 2	pos = open('POSCAR', mode = 'r') line = pos.readlines()

获取 slab 在 z 方向的数值
1
vaccum = float(line[4].split()[2])

将所有原子的 z 坐标保存到 z_list 数列中

1
2
3

z_list = []
for i in range(9, num_lines): 
  z_list.append(float(line[i].split()[2]))

获得真空层的厚度
1
l_vaccum = vaccum - max(z_list)
计算距离表层原子一半真空层厚度的z方向的数值： num_middle
1
num_middle = (max(z_list) + vaccum) / 2

将 [num_middle-1, num_middle+1] 区间范围内所有的真空能级提取出来，并保存到middle_y 这个数列中

middle_y = []

for i in dic.keys():
    i = float(i)
# Select the date area within 1 angstrom from the middle point:
    if i > num_middle -1 and i <  num_middle + 1:
        middle_y.append(float(dic.get(str(i)))

求平均值，并打印出来

1	print 'The Vaccum Energy in this calculation is:\t %s' %(np.mean(middle_y))

运行脚本：

图中第一个命令中，输出比较复杂，并且会和前面示例一样，展示功函数可视化的结果。

但这不便于批量提取数据，可注销脚本中一下三行，直接获取真空能级
1
2
3
#plt.show()
#print 'The Vaccum in this calculation is:\t\t %s' %(l_vaccum)
#print 'The Vaccum Energy in this calculation is:\t %s' %(np.mean(middle_y))
使用脚本批量获取真空能级

先将脚本复制到bin文件夹中，赋予可执行权限。
1. 批量生成 LOCPOT_Z 文件：
2. 使用脚本批量读取
批量获得费米能级：

我们需要的为第 4 列的数据。(图中的冒号也算作一列)

3 数据分析：

真空能级，费米能级以及功函数随真空层厚度的变化

图中的测试曲线看起来很漂亮，但纵坐标范围太大，我们将功函数随真空层的变化作图如下：

从图中，可以看出： 20-30 Å 范围内，功函数的波动变化较小，为 4.81 eV。其他的波动值在 4.85-4.87 eV 之间。Cu(111) 表面的功函数为：4.94 eV。参考的是 CRC handbook 中的数值。

100	5.10	FE
110	4.48	PE
111	4.94	PE
112	4.53	PE

实验与理论偏差为： -0.13 / 4.94 = 2.63 %

师兄，实验值结果是 4.94 eV，如果我们采用 30 Å 时的结果 (4.87 eV)，跟实验的差别更小，是不是更好?

可以这么说。但是

从我们测试的结果可以看出来，功函数随着真空层的变化，存在一个计算的误差；
实验测量值也会有一定的误差范围。

所以，如果你硬要把理论和实验结果完美地吻合，这是很难做到的，这当然也是大家所梦寐以求的终极目标。一般来说，只要在合理的误差范围之内 (5% 的样子)，跟实验值结果一致就可以了，此外，如果你发表文章的时候，为了保证计算的可重复性。计算参数设置，模型尺寸等都要尽可能地详细。

师兄，不是说真空层对功函数的影响很大吗，为什么在 10 Å 的时候，得到的功函数 (4.80 eV) 和 20 Å (4.81 eV) 的时候一样?

这是因为，在当前的体系下，4 层的 p(1x1)-Cu(111) 表面，10 Å 已经可以了。但 10-20 Å 时，误差波动较大。12 Å 的数据也可以用，如同前面所提到的，他们都在一个合理的误差范围之内，但必须注明你的模型尺寸。对于其他的体系，还需要大家自己动手去测试一下。

此外，给大家推荐一篇文章：Fermi level, work function and vacuum level

本文刚刚开始的图片就是从这篇文章中偷过来的。

4 扩展练习：

复习功函数的相关内容；
计算，测试自己研究体系的功函数；
弄清获取真空能级的脚本原理(切记，搞不明白之前，不可乱用！！！)

5 总结：

到本节为止，功函数的计算就先简单介绍到这里了，计算步骤，可视化以及测试都大体讲解了一下。但我们用的模型例子比较简单，如果你的体系复杂，计算仍有疑问，请在 QQ 群中交流或者发邮件给我 (lqcata@gmail.com) 。

Ex51 影响功函数计算的因素（一）

发表于 2018-09-05 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

前面我们通过 Cu(111) 表面作为例子，学会了功函数计算的基本步骤和可视化过程。那么功函数计算的时候，需要注意的事项有哪些呢？

1 INCAR中的参数

计算功函数的参数：LVHAR =.TRUE.

加入这一参数时，VASP 只将静电势能写入 LOCPOT 文件中。

在早期的 VASP 版本中，静电势的写入是通过设置 LVTOT 这个参数的。

在 5.2.12 版本之后，如果你设置 LVTOT= .TRUE.，那么静电势，交换相关势都会写入到 LVTOT 中。由于我们计算功函数的时候，只需要静电势这一部分。所以，如果你用的是 5.2.12 版本之后的 VASP，设置 LVHAR= .TRUE. 即可。

如果有疑惑的话，不妨做个测试，分别设置 LVTOT 和 LVHAR = .TRUE. 然后做个单点计算对比下结果。

设置 LVTOT = .TRUE. 的结果如下：

设置 LVHAR = .TRUE. 的结果如下：

我们对比下 13-15 $\AA$范围内纵坐标的大小和平均值

很显然，两个参数对功函数的影响是不可忽略的。使用 LVTOT 这个参数，由于加入了交换相关势，曲线变得不再那么光滑，并且与 LVHAR 的结果有一定的偏差。所以，在计算功函数的时候，LVHAR 这个参数一定要注意。

这一点也体现在VASP官网最新的ppt中，如下如：(自己主动根据下图中的关键词找这个 ppt，别问我要，也不要在大师兄群里求助浪费别人的时间。)

注意：

在 Hand-on-session (老版本的官方教程)中，使用的是 LVTOT 这个参数。老版本就是过时的意思。希望大家的以新版本的计算为准。

老版本中功函数的计算例子。

2 真空层的厚度及修改

真空层的厚度：指的是 slab 在 z (或者 c )方向上的长度减去表面原子在 z 方向的坐标。

Slab 方向的长度，指的是 POSCAR 或者 CONTCAR 中第 5 行中的数值，上图箭头所指的地方。

那么我们怎么修改真空层的厚度呢？

由于 slab 模型中的原子部分就在那边乖乖地待着，我们只需改变 slab 中晶格常数在 z 或者 c 方向的长度即可。

例子1：上图中真空层的厚度为 15 $\AA$，我们需要一个 20 $\AA$ 的 slab 模型，也就是在 21.2994 的基础上再加 5 个 $\AA$，等于 26.2994。但这样做，对不对呢？修改之后的 POSCAR：

结构如下图：

我们发现 slab 的 Cu 原子部分之间好像也被拉长了。测量了一下两个 Cu 原子之间的距离为： 2.956 Å。

修改之前为： 2.547 Å。

所以：我们在前面的操作中，直接修改的 z 方向的数值，方法是错误的。

原因在于：前面的结构中坐标为分数坐标： Direct

我们修改完成之后，Cu 原子在 c 方向的距离也会发生相应的改变。

所以，如果直接修改 POSCAR 或者 CONTCAR 改变真空层厚度的话，我们一定一定要先将它们转化为 Cartesian 坐标。怎么转化呢？

方法1：用软件操作，比如 p4vasp。

我们可以切换坐标通过鼠标点一下即可，然后保存成 Cartesian 的 POSCAR。

当然啦，也可以使用其他软件，比如 VESTA 等等，更好的选择，也欢迎留言补充。

方法2：使用脚本转换：

VASP 官网在 POSCAR 的解释部分，提到了怎么进行坐标切换的公式。

链接如下：https://cms.mpi.univie.ac.at/vasp/vasp/POSCAR_file.html

在此基础上，本人写了一个 python 的小脚本，可以实现 Direct 到 Cartesian 的转换。

运行如下：

图中流程的解释：

将一个计算中的 CONTCAR 复制过来；
使用 head -n 10 看一下这个 CONTCAR 的文件结构。(10 指的是前面 10 行，如果你想看前面 5 行，使用 head -n 5 )；
dire2cart.py CONTCAR 使用脚本进行转换，转换的对象为 CONTCAR；
转换完成后，Cartesian 的保存为 CONTCAR_C 文件；
cat CONTCAR_C 这个命令查看转化后的内容。(当然也可以继续使用前面的 head 命令)

本脚本下载链接： https://pan.baidu.com/s/1eRMJ7m6 密码：btsl

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#Convert direc coordiation to cartesian Writen By Qiang 
import sys
 
script, file_to_be_converted = sys.argv 
 
print """
###################################
#                                 #
#for VASP 5.2 or higher versions  #
#                                 #
###################################
"""
 
file_read = open(file_to_be_converted, 'r')
 
line = file_read.readlines()
a1 = float(line[2].split()[0])
a2 = float(line[3].split()[0])
a3 = float(line[4].split()[0])
b1 = float(line[2].split()[1])
b2 = float(line[3].split()[1])
b3 = float(line[4].split()[1])
z1 = float(line[2].split()[2])
z2 = float(line[3].split()[2])
z3 = float(line[4].split()[2])
 
num_atoms = sum([int(x) for x in line[6].split()])
 
x_cartesian = []
y_cartesian = []
z_cartesian = []
tf = []
 
start_num = 9 # Default: With Selected T T T, coordination starts from line 9
 
def convert():
     
    for i in range(start_num, num_atoms + start_num):
        x_cartesian.append(float(line[i].split()[0]) * a1 + float(line[i].split()[1]) * a2 + float(line[i].split()[2]) * a3)
        y_cartesian.append(float(line[i].split()[0]) * b1 + float(line[i].split()[1]) * b2 + float(line[i].split()[2]) * b3)
        z_cartesian.append(float(line[i].split()[0]) * z1 + float(line[i].split()[1]) * z2 + float(line[i].split()[2]) * z3)
        if len(line[i].split()) > 3:   # if  T T T exist, there are more than 3 elements in the list line[i].split()
            tf.append((line[i].split()[3]))
        else:
            tf.append(' ')   # if there is no T T T, use space instead. 
     
    file_out = open(file_to_be_converted+'_C', 'w')
     
    for i in range(0,7):
        file_out.write(line[i].rstrip() + '\n')  # first 7 lines are kept the same 
     
    if 'S' in line[7]:
        file_out.write(line[7].rstrip()+ '\n')  # if  T T T exists, write the Selective line 
    file_out.write('Cartesian' + '\n')          # Coordination system is Cartesian now. 
     
    for i in range(0,len(x_cartesian)):
        file_out.write("%+-3.10f   %+-3.10f   %+-3.10f   %s %s %s\n" 
        %(x_cartesian[i], y_cartesian[i], z_cartesian[i], tf[i], tf[i], tf[i]))
     
    file_out.close()
    print '-----------------------------------------------------\n'
    print 'POSCAR with Cartesian Coordiations is named as %s_C\n' %(file_to_be_converted)
    print '-----------------------------------------------------\n'
 
if line[7][0]  == 'S' or line[7][0]  == 's':  # # With Selected T T T, coordination starts from line 9
    start_num = 9
 
    if  line[8][0]  == 'D' or line[8][0]  == 'd':
        print """
This POSCAR has Direct Coordinations, Conversion is starting....
 
              """
        convert()
     
    elif  line[8][0]  == 'C' or line[8][0]  == 'c':
        print """
This POSCAR has Cartesian Coordinations! Process is aborted!
               
              """
 
 
else : 
    print """
----------------------------------------------------
Pay Attetion! There is no TTT in coordinations part!
----------------------------------------------------
"""
     
    start_num = 8 # without Selected, No  T T T , coordination starts from line 8 
     
    if  line[7][0]  == 'D' or line[7][0]  == 'd':
        print """
This POSCAR has Direct Coordinations, Contersion starts....
 
"""
        convert()
     
    elif  line[7][0]  == 'C' or line[7][0]  == 'c':
        print """
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
This POSCAR has Cartesian Coordinations already! 
 
Process is aborted!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
"""
 
file_read.close()

那么怎么把 Cartesian 转化为 Direct 呢？

使用 VASP，因为 VASP 的默认输出就是 Direct 坐标，算个单点就可以啦（笑话，别当真！）；
使用 p4vasp 等其他软件；
写个脚本倒着再做一遍，不过本人经常使用的是 Cartesian 坐标，懒得再写了，有兴趣的可以自己试试。

3 批量处理POSCAR

当我们完成转化后，就可以批量处理 Cartesian 坐标的 POSCAR 了。

先准备一个文件夹，名字为 10，这个文件夹中有一个真空层为：10 Å 的 POSCAR 以及 INCAR, KPOINTS, POTCAR，任务脚本；

运行命令：

1	for i in $(seq 12 2 36); do cp 10 $i ; sed -i "5s/16/$((6+$i))/g" $i/POSCAR; done

示例演示：

运行完这个命令后,会生成从 12 到 36 的 N 个文件夹，每个文件夹之间间隔为 2。我们通过下面这个命令查看所有文件夹中 POSCAR 中 z 方向的大小：
1
for i in *; do head -n 5 $i/POSCAR | tail -n 1 ; done

解释：

head -n 5 $i/POSCAR 获取POSCAAR前5行中的内容，
后面跟着一个 |tail -n 1，这个命令的意思是，显示前面5行中的最后一行。
head 和 tail 这两个命令之间用 | (pipe) 连在一起，表示将前面 head 命令的结果传递给后面的 tail 命令。
任务准备好之后，批量提交就可以了。(不会的话，请前面自己翻批量操作的介绍。)

4 扩展练习：

使用 LVTOT 和 LVHAR 测试一下，加深下印象；
使用不同的测试不同真空层对功函数的影响以遍下一节的学习。

5 总结：

下一节我们讨论一下真空层对计算功函数的影响，以及如何批量获取真空能级。

Ex50 功函数（work function）的计算（二）

发表于 2018-09-04 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

前面一节我们介绍了计算功函数的具体流程，以及用p4vasp获取真空能级的方法。

本节我们介绍另外两个方法，从本质上来说，这些方法是一样的，都是基于对 LOCPOT 中数据的处理。

1 使用 vtotav-v5.2f 脚本

操作流程如下：

注意：

gfortran -o vtotav.x vtotav-5.2f

这个命令将 vtotav-5.2f 这个脚本编译成可执行文件 vtotav.x。

gfortran 不管用的话可以试试ifort -o vtotav.xvtotav-5.2f这个命令。
./vtotav.x 在当前目录下运行这个可执行文件;会让你选择计算的方向，这里我们 slab 的真空层是在 z 方向，所以我们输入 3 (回想一下，上一节，我们在 p4vasp 中选择z方向的鼠标操作，它们是一样的)，都是讲 LOCPOT 文件处理，然后生成z方向的相关数据。
运行完毕后，当前目录下会多出来一个 VLINE 文件。这个文件就是 p4vasp 中选择 z 方向之后所显示的红线数据，(不过横坐标有些不同)。
我们看一下 VLINE 文件的数据内容，第一行中 224 代表有 224 个点，3 代表 z 方向。
VLINE 文件中横坐标为当前点的序号，从 1-224，而 p4vasp 中导出来的文件，横坐标为 z 方向的坐标，不过这个不影响我们对真空能级的计算。

上图中最后的两行命令，

从上一级目录下拷贝了一个 python 脚本，名为：wplot-f.py；
运行这个 python 脚本，效果如下：

我们又得到了前面一节的功函数的图像。真空能级，大家从 150 前后取个值即可，也可以在一段直线的区间范围内求个平均数。

对比下前面一节和本节中的数据：

左图为本节的，右图为上节的。

Python 脚本这么神奇，是怎么写出来的呢？脚本内容如下：

#!/usr/bin/env python
# Written By Qiang for workfunction Visualization from VLINE file

import matplotlib.pyplot as plt

x = []
y = []
with open("VLINE", mode='r') as f:
    next(f)
    for line in f:
        xy=line.rstrip().split()
        x.append(float(xy[0]))
        y.append(float(xy[1]))

plt.plot(x,y) 
plt.savefig('workfunction' + '.pdf', dpi=400)
plt.show()

思路如下：

读取 VLINE 文件；
第一行画图的时候用不着，需要跳过，next(f)；
将剩下的那些行中，第一个数存到 x 列表里面，第二个数存到 y 列表里面；
使用 matplotlib 中的 pyplot 读取 x， y 列表中的数据，进行画图： plt.plot(x,y)；
将生成的图保存成 workfunction.pdf 文件，分辨率 dpi 为 400。

注意：

如果你想保存成 eps, png, jpg,只需将 .pdf 改成 .eps, .png, .jpg 即可；
如果想要分辨率更高，可以修改400这个数值；
如果只想保存成 pdf 文件，不想弹出图片的查看窗口，把 plt.show() 这一行注释掉就可以了， # plt.show()；
如果只想看一下图，不想保存 pdf 文件，将第 16 行注释掉就可以了。

2 使用 vtotav.py 脚本

第一个脚本就讲到这里了，这个脚本不错，唯一的缺点就是横坐标不是 z 方向的数值，那么我们看一下第二个脚本的操作。

注意：

这里用到的脚本名字为：vtotav.py；
vtotav.py LOCPOT z 这个命令意思是，用 vtotav.py 读取 LOCPOT 文件，并处理生成 z 方向的数据；
脚本运行完毕后，目录中多了一个 LOCPOT_Z 文件，这个文件内容如下：
图中最后两行，
1. 大师兄又从前面的目录下拷贝了一个脚本： wplot.py；
2. 然后运行这个 python 脚本，效果如下：

这个脚本中，横坐标为 z 方向的大小，比前面的脚本有所改进，和 p4vasp 的数据是一样的。首先我们先看一下可视化的脚本，然后再讲解 vtotav.py 这个脚本的获取和使用。

#!/usr/bin/env python
# Written By Qiang for workfunction Visualization from LOCPOT_Z file

import matplotlib.pyplot as plt

x = []
y = []
with open("LOCPOT_Z", mode='r') as f:
    first_line = f.readline()
    name_x = first_line.split()[1]
    name_y = first_line.split()[2]
    for line in f:
        xy=line.rstrip().split()
        x.append(float(xy[0]))
        y.append(float(xy[1]))

plt.plot(x,y) 
plt.xlabel(name_x) 
plt.ylabel(name_y) 
plt.savefig('workfunction' + '.pdf', dpi=400)
plt.show()

这个跟前面的那个脚本很像，不同的地方在于这个脚本中读取了第一行的内容，并将它们作为横坐标和纵坐标的名字。(Distance 和 Potential)

3 vtotav.py 脚本的获取和使用

这个脚本可以在github网站免费下载，网址为：

https://github.com/compphys/ase_tools/blob/master/scripts/vtotav.py

本人的做法是：

把图中的代码复制下来，然后保存到一个文本里面，将文本命名为：vtotav.py；
赋予脚本可执行权限： chmod u+x vtotav.py；
将脚本移至 ~/bin 目录下： mv vtotav.py ~/bin；
然后直接用就可以了，跟前面图中的操作一样。

不过本人刚刚发现这个脚本可以按照下面的方法下载，下载完成后，重复前面的 2-4 步即可：

点击 scripts
鼠标右击这个脚本
另存为，然后选择保存的目录
这里本人保存到电脑桌面上了

注意：

如果你运行脚本的时候，没有得到下面图中的这个结果：

原因可能在与：

这个脚本需要调用 ase 这个程序中的模块，那么什么是 ase 呢？ase 怎么安装呢？

对于 ase，本书的 ex0 中就已经提到了，可能大家也忘记了。ASE 是 Atomic Simulation Environment 的缩写，这是一个非常强大的工具。官方网址：

https://wiki.fysik.dtu.dk/ase/

是大大师兄(本人的师兄)做博后的课题组开发出来。Linux 系统下面，安装 ASE 很简单：

一个命令就可以搞定了，Soooo eeeeeasy! 至于 Windows 下面嘛，本人不会。不想用 Linux，还想用 ASE 的就只能自己捣鼓捣鼓了。捣鼓好了也不要给我发教程，因为本人不推荐用 Windows 做计算。

4 扩展练习：

下载本节的例子，以及所有的脚本，按照讲解的操作一遍流程，链接：https://pan.baidu.com/s/1dEHvHr7 密码：ponm；
Linux 用户自行安装 ASE，使用 vtotav.py 这个脚本走一遍过程；
熟悉下用 python 画图的基本方式。

5 总结：

讲解到这里，功函数的可视化也就差不多了。当然了，只要你有了数据，用什么画图都可以， excel, origin, matlab, gplot… 不过貌似作图也没什么用，我们大体上看下，然后找个区间取静电势的数值才是正事。

Ex49 功函数（work function）的计算（一）

发表于 2018-09-03 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

前面我们讲完了表面的计算，后面两节我们主要讨论一下功函数的计算。

1 功函数的定义

首先，我们先看一下维基百科的解释：

功函数(又称功函、逸出功)是指要使一粒电子立即从固体表面中逸出，所必须提供的最小能量(通常以电子伏特为单位)。这里“立即”一词表示最终电子位置从原子尺度上远离表面但从宏观尺度上依然靠近固体。功函数不是材料体相的本征性质，更准确的说法应为材料表面的性质(比如表面暴露晶面情况和受污染程度)功函数是金属的重要属性。功函数的大小通常大概是金属自由原子电离能的二分之一。

The work function W for a given surface is defined by the difference

W = - eϕ - E_F

where −e is the charge of an electron, ϕ is the electrostatic potential in the vacuum nearby the surface, and E_F is the Fermi level (electrochemical potential of electrons) inside the material. The term −eϕ is the energy of an electron at rest in the vacuum nearby the surface. In words, the work function is thus defined as the thermodynamic work required to remove an electron from the material to a state at rest in the vacuum nearby the surface.

我们再看下 IUPAC 官网的解释：

The minimum work needed to extract electrons from the Fermi level of a metal M across a surface carrying no net charge. It is equal to the sum of the potential energy and the kinetic Fermi energy taken with the reverse sign:

$\phi^{M} = -(V_{e}+\varepsilon_e^ F)$

where $V_{e}$ is the potential energy for electrons in metals and $\varepsilon_e^ F$ is the kinetic energy of electrons at the Fermi level.

2 VASP 计算功函数的过程

从前面的定义中可以看出，计算功函数我们只需要得到体系的费米能级和电子所处的静电势能，然后求差即可。

2.1 费米能级的计算：

VASP 计算结束之后，通过命令：

grep E-fermi OUTCAR 即可提取出来。

2.2 静电势能的计算：

通过在 INCAR 中添加： LVHAR =.TRUE. 这个参数。

加入这个参数，计算结束后，VASP 会输出一个文件： LOCPOT 文件。我们可以通过脚本，或者程序对这个文件后处理来获取静电势能，怎么处理后面会讲的很清楚，大家先不要着急心慌。

2.3 VASP的计算流程：

优化获取稳定的结构
将 CONTCAR 复制成 POSCAR
在 INCAR 中设置单点计算：

方法1：设置 NSW=1 或者 0

方法2：设置 IBRION = -1

方法3：直接删除 INCAR 中 NSW 这一行，采用默认值
提交任务

注意：
1. 如果第一步计算中有保存 WAVECAR，则静态计算的时候可以设置：ISTART = 2 读取一下，以便节省计算时间;
2. 如果没有保存，直接算即可。不要在读不读 WAVECAR 这个问题上瞎操心；
3. 电子步数的默认值是 40，如果你的体系不容易收敛，计算的时候请设置一个较大的 NELM 值，比如 NELM = 500。

3 Cu(111) 表面的功函数计算

3.1 准备计算文件：

新建一个 workfunction 的文件夹，将前面优化计算的： INCAR KPOINTS POTCAR CONTCAR 和提交任务的脚本复制到这个文件夹中。

3.2 修改 INCAR 改为静态（单点）计算

3.3 修改 INCAR，添加 LVHAR = .TRUE. 这一行

System = Cull
ISMEAR = 1 
SIMGA = 0.1
ALGO = FAST
ENCUT = 450
EDIFF = 1E-5

IDIPOL =.TRUE.
IDIPOL = 3

NWRITE = 0
LWAVE = .FALSE.
LCHARG = .FALSE.
IBRION = 2
POTIM = 0
EDIFFG = -0.01
NSW = 1 

LVHAR = .TRUE.

3.4 将 CONTCAR 重命名为 POSCAR

3.5 提交任务，等待结束。

4 后处理方法-1：p4vasp

p4vasp 一直是本人强烈推荐的VASP计算相关的建模，可视化，查看结构，结果的软件。本节我们讲解一下如何通过 p4vasp 对功函数的计算结果可视化。

4.1 导入数据

Linux 系统下面，进入计算的目录直接敲命令： p4v & 即可 (我们稍微等下 Windows 用户的操作，后面可视化过程是一样的)
Windows 下面，首先要将计算生成的 LOCPOT 文件保存到本地电脑上。

注意：有时候这个文件可能会很大，下载的时候比较蛋疼，不想下载的话等待下节的讲解。
打开 p4vasp 并导入 LOCPOT 文件：file —> loadsystem，然后找到 LOCPOT 所在的目录，点解导入，效果如下

4.2 开始可视化 (Linux, Windows 都一样了)

Electronic –> Local potential

效果如下：

WOC，这是什么东东？是不是算错啦？

别急，看下面：

点击下图中的 Show，会出现三个方向让你选择，一般都是沿着 z 方向的，按照红色箭头点下即可。

效果如下：

Wow, Wonderful!!! 是不是很有成就感？

But，怎么才能获取静电势能的数值呢？图中的纵坐标范围太大了，目测取值的话，误差应该在 0.5 eV 左右，太大了，不行。那么我们只能导出数据，慢慢分析了。

4.3 导出数据

还记得前面我们讲的 DOS 的计算吗？我们求 d band center，电子数目积分的时候将数据导出来了。同样，图中组成曲线的点的坐标我们也可以导出来。
点击： Graph—>Export

选择保存的目录，并给保存的文件一个名字。然后点击，Export，图中的 cu.dat 就会保存你指定的目录下了。

注意：

名字后面一定要有尾缀.dat或者 .data，否则会导出失败；
如果按照我说的步骤做的话，还不会弄，请自行查看前面关于 DOS 的计算，已经讲的很清楚了，不要再问我这个问题了。

4.4 分析数据

有了数据，我们就底气十足，什么都不怕了。用 Notepad ++ 打开 cu.dat 文件。

注意：这里大师兄说的是用 Notepad++，别问我从哪里下载这个神奇的软件。百度会告诉你的。

打开之后，我们会发现，和 DOS 计算导出来的结果有几分神似。数据共有三部分，每两部分之间由空行分开。

这三部分的顺序，如下图，

我们从图中可以看出，基本上在 10-18 的这个区间范围，三条线的数据差不多是一样的。静电势能的话，可以通过红线( Average 那条线)，也就是第二部分来获取。可以查看下 14-15 $\AA$ 时红色曲线的数值，如下图：

4.5 获取费米能级

套公式，功函数为： 4.84-0.03 = 4.81 eV

对比下维基百科中的数值：

注意：虽然维基百科给了功函数的一些数值，在发表文章的时候，你不能直接放 wiki 的链接到论文里面，这显得有些不专业。建议大家去参考 CRC handbook 里面的数值，然后引用。别问我 CRC handbook 是啥玩意，前面我也讲过了。

5 扩展练习：

重复本节的计算；
重复 VASP 官网 Ni(111) 表面的计算。

6 总结：

本节我们主要

介绍了下功函数的定义(从别的地方拷过来的)；
介绍了下 VASP 怎么计算功函数；
怎么用 p4vasp 可视化，导出数据，并获取静电势能；
怎么获取 fermi level；
举例介绍了p(1x1)Cu(111) slab 的功函数计算流程，计算结果下载：

链接：https://pan.baidu.com/s/1c2m52UG 密码：kfse

Ex48 表面能的计算（二）

发表于 2018-09-02 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

前面一节，我们讲解了一下 VASP 官网 Hand-On-Session-III 中表面能的计算。

1 首先我们要指出一点 VASP 官网的错误

最后的表面能应该是： 0.86 - 0.015 = 0.845 eV 而不是 0.71 eV。

2 为什么要除以 2？

很多人困惑为什么我们在计算的时候，上图中要除以 2。这是因为我们的单点计算中包含有 2 个面。前一节解释说是因为把 Bulk 切成了两半，以至于把大家给弄糊涂了。为什么我们的计算中会有两个面呢？首先我们看下 Slab 模型：

问题：上图的结构中有几个表面呢？

答案是 2 个，因为我们的 slab 模型在 z 方向上也具有周期性，虽然我们把底部的两层原子固定来模拟体相，但实际上它们还是表面。如下图：

那么，为什么 E^rel 没有除以 2 呢？

这是因为我们只优化了一个表面啊。

如果我们两个表面都要优化的时候，是不是 E^rel 也除以 2 了呢？

答：是的。

如果你喜欢上下两个面都优化的话，那么这将引申出 slab 模型的两种情况，对称性和非对称性模型。前面我们讲的都是非对称的 slab 模型。那么什么是对称的 slab 模型呢？如下图：

在这个模型中，中间的三层原子被固定用来模拟体相，上下三层的原子放开用来模拟表面。此时，如果你计算表面能的时候，E^rel 的数值就要除以 2 了。计算公式可以简化为：

$\begin{align*} \sigma &=\frac{1}{2A}(E_{surf}-N_{atoms}\times E_{bulk}+E^{rel})\\ &=\frac{1}{2A}(E_{surf}^{rel} -N_{atoms}\times E_{bulk})\\ \end{align*}$

其中： $E_{surf}^{rel}$ 直接就是我们优化完毕之后的对称slab模型的能量。

单个表面的面积 A， slab 模型中的原子数 $N{atoms}$ 和 $E{bulk}$ 这些能量的获取和前面一节讲的内容一致。

注意：将这个公式和参考书 96 页底部的公式对比一下。然后仔细阅读参考书 97 页的公式各项的说明。

3 对称和非对称slab模型的区别：

对称 slab 明显具有更多的原子数，在 z 方向更长，需要更多的真空层(上下两层)，这就不可避免地出现了计算量大的缺点，尤其是对于机时捉襟见肘的筒子们来说，这个模型就不用考虑啦；
非对称模型由于其非对称性，在计算过程中会产生偶极矩，注意这个偶极矩是由于模型引入的，并非我们常说的分子中的偶极矩。尤其是当我们要计算一些表面吸附或者反应的时候，两个 slab 之间偶极矩的相互作用会对我们的结果产生影响。但这不是什么大问题，很多计算软件都可以通过控制计算参数来消除或者将这一影响减至最小。
VASP 中可以通过设置：
1. LDIPOL=.TRUE.(打开偶极校正)
2. IDIPOL = 1-4来解决。 1，2，3分表代表在x，y和z 方向上进行校正。4代表在所有方向上。

注意：1 和 2 必须同时加在 INCAR 里面。详见VASP官网的解释： LDIPOL 和 IDIPOL

4 表面能计算的参数影响：

从公式中可以看出 Slab 的能量和 Bulk 的能量是主要因素：

Slab 能量：
- Slab的层数：如参考书中的结果。
- slab 表面的大小，一般来说 p(1x1) 的即可。当然啦，你可也比较下 p(1x1) 和 p(2x2) 的区别：这里要注意的是：改变表面的大小，KPOINTS 也要发生相应的变化，只有这样，两个计算的结果才具有可比性。
- 真空层的厚度这些大家都可以测试一下。
Bulk 的能量前面一节我们已经讲过 Bulk 的能量计算了。前面有同学问到： Bulk 计算的时候，模型中应该包含多少个原子？
个人认为：单胞或者原胞的计算均可，当然你也可以扩一下晶胞。但需要注意的是如果你要比较这些不同大小模型中单个原子的能量时，一定要注意 K 点的选择。如果 K 点密度不一致，导致单个原子的能量之间存在细小的差别，这会令你不知所措。有兴趣的可以测试一下。
此外，大师兄还要推荐一本进阶版的参考书：Theoretical Surface Science: A Microscopic Perspective，作者是 Axel Groβ。课题组主页：https://www.uni-ulm.de/en/nawi/institute-of-theoretical-chemistry/

本书关于表面能的介绍在第68页，具体内容大家认真阅读，必定受益匪浅。

5 扩展练习：

请务必认真阅读参考书中表面能， slab 模型的相关介绍；
认真测试，重复 Cu(111), Cu(110) 和 Cu(100) 面的计算；
仔细阅读 Theoretical Surface Science: A Microscopic Perspective 这本书中关于表面能的内容；
结合计算结果和参考书中的解释，自己参悟一下。

6 总结：

本节主要回复了一下大家提问的比较多的问题，简单介绍了一下对称和非对称 slab 模型以及一些影响表面能计算的因素。学习 VASP 没有捷径，大家多练习，多出错，多总结分析，时间长了就会慢慢理解了。

Ex47 表面能的计算（一）

发表于 2018-09-01 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

继续前面的计算，这一节我们主要讲一下表面能的计算过程。

1 表面能定义

首先，我们先介绍一下表面能的定义。为了偷懒，本人直接摘抄的维基百科的解释，如下：

表面能是创造物质表面时，破坏分子间化学键所需消耗的能量。在固体物理理论中，表面原子比物质内部的原子具有更多的能量，因此，根据能量最低原理，原子会自发的趋于物质内部而不是表面。表面能的另一种定义是，材料表面相对于材料内部所多出的能量。把一个固体材料分解成小块需要破坏它内部的化学键，所以需要消耗能量。如果这个分解的过程是可逆的，那么把材料分解成小块所需要的能量就和小块材料表面所增加的能量相等。但事实上，只有在真空中刚刚形成的表面才符合上述能量守恒。因为新形成的表面是非常不稳定的，它们通过表面原子重组和相互间的反应，或者对周围其他分子或原子的吸附，从而使表面能量降低。

2 VASP计算

具体到VASP的计算中，表面能是怎么算出来的呢？下面我们看一下VASP官网的表面能的计算过程。

公式先大体浏览一遍，我们分析下里面的各项所表示的内容：

E^rel： Relaxation energy，即弛豫能量。这个前面一节我们已经讲过了，从刚刚切好的表面优化到稳定的表面所释放的热量。图中的第一个点就是刚刚切好表面的能量（E^surf），最后一个点就是优化完成之后的能量。这个相信大家现在都能理解了。
σ：在handonsession-III里面，这个用来表示表面能。
E^surf：刚刚切好的 slab 能量；可以直接算个单点，也可以取优化过程中第一个离子步的能量。
N_atoms： Slab 中的原子数目
E_bulk： Bulk结构中单个原子的能量。

注意1：

前面说的最后一点。E_bulk，这个不是指的体相的能量，而是体相中单个原子能量。即体相的能量除以原子数目。我们顺便回顾一下前面单胞的优化：

获取稳定的晶格常数，两个方法：
- ISIF + Large ENCUT
- 通过计算不同的晶格常数的晶胞能量，用BM方程拟合
得到稳定的晶格常数之后，如果用的 ISIF= 3 这个方法，则将 CONTCAR 复制成 POSCAR ，然后在 INCAR 中设置： ISIF = 2 + 正常的 ENCUT ，做个单算，算下计算晶胞的能量。
用晶胞的能量除以晶胞中的原子数目，就得到了这里的 E_bulk。比如前面 Cu 的单胞中有4个原子，我们则用计算后体系的能量除以4。

注意2：

Handonsession-III 里面的表面能，严格来说，这个定义其实是不对的。因为表面能，从字面上理解，它的单位应该是：能量/面积。而图中只有能量，没有面积。所以，图中的表面能指的是纯能量。

σ^unrel：将块体直接切开所用的功，因为是切成了2份，公式中除以2。

注意，在这个计算里面，暂且不考虑表面的弛豫，单纯从体相切成2个表面的过程。
能量是刚刚切出来的slab的能量： -25.560
从体相到切开表面直至达到稳定状态，总的能量变化为：
$\sigma = \sigma^{unrel + E^{rel}}$

注意3：

E^rel 这里没有除以2。为什么呢？这是因为我们只放开了上面的两层原子。也就是只优化了一个表面。
至此，从体相切成表面，然后直至表面稳定，这一个过程的能量就算完了。
既然前面说了，表面能里面有面积，那么面积指的是什么面积？怎么求面积呢？

面积就是你计算的slab的那一层表面积，很多人问我这个问题，我也是无语了！
怎么求呢? 我们先看一下Slab模型(Notepad++这个软件打开的！没有的话赶紧去安装一个)

前面几行放大一下，看看晶格尺寸。

如果面积用A来表示，则：

A = 2.57170 2.22716 – 0.000 (-1.28585)

注意4：

这里的单位是 Å²
原理如下图：(自画自拍，不喜勿喷!)

菱形的面积等于： x1 * y2
如果坐标系不是那么正儿八经对在一起的时候，如下图：

则面积等于：x1 y2 – y1 x2
Handonsession-III 中 Ni(100) 表面的面积为：

A = (0.50.5– 0.5 (-0.5)) = 0.5 Å²
这里我们还忘了乘以前面的缩放系数：所以，
A = 0.5 3.53 3.53 = 6.23045 Å²
最终表面能为： 0.71 eV / 6.23045 Å² = 0.114 eV/Å²
一般来说表面能的单位为： J/m²
换算系数为： 1 eV/Å² = 16.02 J/m²
换算完毕后，表面能为： σ = 0.114 * 16.02 J/m² = 1.826 J/m² 。

注意5：我们的参考书中(第97页中)给弄反了。

这是因为：
1 eV = 1.60217648710^-19 J
1 Å² = (110^-10) (110^-10) m²
大家自己算算就知道了。讲到这里，相信金属表面能，大家都可以照着去算一遍了。

注意6：

获取E^surf 的时候，一定要检查：1）单点计算是不是收敛了？或者2）第一个离子步是不是收敛了？

如果没有收敛的话，表面能的计算就是错的。σ^unrel和E^rel都跟这个量密切相关。

3 下面我们看一下这篇参考文献中的计算介绍：

Shapecontrol in concave metal nanoparticles by etching
Nanoscale, 2017,9, 13089-13094

这是本人博士期间的一个工作，主要讨论了不同盐酸浓度下Pt纳米颗粒的成长，腐蚀过程以及形貌的预测。有兴趣的可以看看，当然啦，本人也希望大家多多引用。对于表面能的计算，大家可以参考支持信息中的公式1，如下图：

相信到这里，大家可以完全明白这段话的意思，以及计算的整个流程了。

4 扩展练习：

重复一下 Hand on Session -III 中 Ni(100) 表面能的计算;
计算下我们前面的 Cu(111) 表面能的计算;
计算下 Cu(100), Cu(110) 的表面能。
对比 Cu(111), Cu(110) 和 Cu(100) 的表面能，分析下规律。
思考1：表面能的计算有什么需要注意的，请公众号留言回答！
思考2：表面能能用来干什么？请公众号留言回答！

5 总结：

本节，我们主要讲了表面能的计算，参考的资料有：

VASP 官网的 Hand on session (如果你自己下载不到，请认真反思自己！)
我们的参考书： Density Functional Theory: A Practical Inroduction
以及本人发表的文章一篇。
希望大家可以认真学习，不要急躁！！！

本书中只要是用到的，提到的资料，本人都已经明确说明出处，获取方式以及下载链接了，请认真浏览，不要再浪费时间找我索要！

Ex46 表面弛豫的计算结果分析以及Direct坐标转换脚本

发表于 2018-08-30 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

前面一节我们回顾了一下表面弛豫的基本概念，对于大多数的金属体系来说，表面弛豫一般发生在垂直方向上。将金属沿表面方向切成两份后，由于表面原子的配位数发生了变化，我们可以想象:该表层原子与下面一层的结合能力会更强，从而降低体系的能量。具体体现在层间距和能量两个方面。我们先回顾一下Ex44中Cu(111)表面的练习。

1 结构优化前后变化

首先观察一下结构的变化， POSCAR的结构如下：(刚刚切好的，热乎的表面)

Cu\(1\1\1)
   1.00000000000000
     2.5717000960999998    0.0000000000000000    0.0000000000000000
    -1.2858500480999999    2.2271576141999998    0.0000000000000000
     0.0000000000000000    0.0000000000000000   21.2994003296000010
   Cu
   4
Selective
Cartesian
+0.0000000000   +0.0000000000   +0.0000000000   F F F
-0.0000128765   +1.4847792201   +2.0999078998   F F F
+1.2858629618   +0.7423784587   +4.1996028482   T T T
+0.0000000000   +0.0000000000   +6.2995107693   T T T

优化后的结构：冷却后的，转化为Cartesian坐标的CONTCAR)

Cu\(1\1\1)
   1.00000000000000
     2.5717000960999998    0.0000000000000000    0.0000000000000000
    -1.2858500480999999    2.2271576141999998    0.0000000000000000
     0.0000000000000000    0.0000000000000000   21.2994003296000010
   Cu
     4
Selective dynamics
Cartesian
+0.0000000000   +0.0000000000   +0.0000000000   F F F
-0.0000128765   +1.4847792201   +2.0999078998   F F F
+1.2858542192   +0.7423835062   +4.1890255872   T T T
+0.0000182442   -0.0000105333   +6.2592676588   T T T

表层原子的坐标从6.2995 $\AA$减小到6.2593 $\AA$。说明该层原子向体相收缩了。
第一、二层的原子间距（单位： $\AA$）为:
6.2995-4.1996 = 2.0999 (POSCAR)
6.2593-4.1890 = 2.0703 (CONTCAR)，
前后变化了(2.0999-2.0703) / 2.0999 * 100% = 1.4 %
注意：这里没有区分正负号，一般来说，用负值表示向体相收缩！！！(详见参考书)

2 弛豫能量

那么怎么知道弛豫前后的能量变化呢? 这个问题很容易回答:

第一步：我们对刚刚切好的表面算个单点计算；获得一个能量。
第二步：用优化完的能量减去前面单点能结果即可。

实际计算过程中，第一步可以免去，因为在优化的时候，VASP会对初始结构计算一下，对应的为第一个离子步的能量。

qli@bigbro:~/test/cu/cu111/opt1$ grep '  without' OUTCAR  | awk '{print $7}'
-13.96892338
-13.96932426
-13.96974724
-13.97030280
-13.97016827
-13.97041579
-13.97047584
-13.97063103
-13.97082922
-13.97086603

所以，弛豫前后的能量（eV）变化为:
-13.97086603 - (-13.96892338) = -0.00194265

能量为负值，说明弛豫这个过程是放热的，也就是说刚刚切开的表面不稳定，表层原子向体相收缩后，体系能量降低，变得更稳定。
注意：这么做的时候，一定要先检查第一个离子步中的电子步是否收敛。比如：VASP默认的单个离子步中的电子步数为60步（NELM = 60），https://cms.mpi.univie.ac.at/wiki/index.php/NELM
如果你设置的的NELM为默认值，且第一个离子步中，到了60步还没有收敛，则该这一步的能量是不可以用的。
此时，你需要做的就是：

尝试增大 NELM = 100 或者其他更大的数值，对未优化的结构重新做个单点计算。
收敛很困难的话，可以调节下ALGO这个参数。

3 dire2cart.py

本节讲的太少，留个转化Direct为Cartesian的脚本(dire2cart.py）作为补充。下载链接：

https://pan.baidu.com/s/1ScsLWhLAPpSul0SzYNPvpg。

上面CONTCAR的转化就是用的这个脚本。使用如下图:

qli@bigbro:~/test/cu/cu111/opt1$ ls
CONTCAR  INCAR  KPOINTS  OSZICAR  OUTCAR  POSCAR  POTCAR  sub4  vasprun.xml
qli@bigbro:~/test/cu/cu111/opt1$ dire2cart.py POSCAR
###################################

# #

#for VASP 5.2 or higher versions  #

# #

###################################

## This POSCAR has Direct Coordinations， Conversion is starting....

## POSCAR with Cartesian Coordiations is named as POSCAR_C

注意：

如果你的POSCAR或者CONTCAR已经是Cartesian坐标了，这个脚本可以自动识别并终止转换。
转化后的文件为：AAA_C, AAA是你要转化的文件（POSCAR或者CONTCAR）。

4 扩展练习

下载VASP官网的关于表面相关计算的pdf文档：VASP Tutorial: A bit of surface science
根据自己所学，尝试着练习里面的相关内容。
其他相关的材料：http://www.nersc.gov/users/training/events/3-day-vasp-workshop/

Ex45 表面弛豫

发表于 2018-08-28 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

答辩心得

顺利完成博士答辩，现在我们开始继续前面的章节。
首先，想跟大家分享下自己答辩的体会：基础知识太不扎实，前瞻性又太不够。导致回答问题的环节，被大牛们问的就想找个缝钻进去了。这对于我的启示就是，以后写这本书的时候，每个练习尽量分成4个部分:
1) 第一部分讲基础知识，
2) 第二部分讲计算过程的细节部分，
3) 第三部分讨论计算结果和基础知识;
4) 相关发表的文章中的计算。
这对大家的要求就是:
1) 有什么任何有疑问的地方，自己要记下来，主动查阅资料，然后反馈给我;
2) 有什么欠妥的地方，恳请批评指正;
3) 如果你有相关发表的文章，可以放在前面提到的第四部分供大家研究，欢迎分享;
4) 这本书的作者不是我自己一个人，更是我们大家的心血，希望我们一起努力;
5) 本人只接受邮件咨询或者批评: lqcata@gmail.com 因为群里讨论的很多宝贵的经验总是没法及时完整的保留下来，以至于造成了很多损失。

表面弛豫

我们先想想一下，有一块金属，我们从中间切开后，新生成的两个表面上的原子和体相中的有什么区别?
首先，表面上的原子配位数比体相中的少了，即表层原子具有悬键。因为切开的过程就等同于把两个表面的原子键切开了。其次，配位数的减少使得表面层与内部的相互作用会更强，与体相中有所不同。这就体现在层间距(键长)上面。

1 维基百科的解释

理想情况下的晶体向各方向无限延伸，其中任一原子的平衡位置由晶体中其他所有原子对其作用力的总和决定。因此，每个原子所在的位置在理想晶体中是等价的，生成的晶体结构是周期性的。现实情况下的晶体大小有限，这就导致了靠近晶体表层的原子受到的作用力不同于晶体内部的原子，从而造成晶体表面原子排列方式不同于晶体内部的结果。晶体表层原子的这种行为可被分为弛豫(relaxation)和重构(reconstruction)两种情况。
弛豫指表面的原子层整体相对于内部的本体(bulk)的位置变化。较为常见的情况是垂直方向上的上下位移，即法向弛豫(normalrelaxation)。大部分金属表面上的弛豫都是这种类型。某些材料的表面也可能在发生法向弛豫的同时有切向的弛豫。

2 Surface Science: An Introduction 这本书的解释

3 参考书的解释:
Density Functional Theory: A Practical Introduction，第94页，第四章，第4.5节 SURFACE RELAXATION 这一部分(请务必认真看完这一部分)

4 扩展练习:

1) 了解表面弛豫的相关内容;
2) 查找相关书籍;
3) 认真阅读参考书4.5节的内容，并重复4.5 节的练习
4)下载本节的两本参考书：QQ群共享或者百度网盘链接：http://pan.baidu.com/s/1eSB2Nx4

5 总结

表面弛豫相信大家都很了解了，借此机会供大家复习一下。关于表面化学的经典书籍，如果你有特别喜欢的，也请推荐给大家。

Ex44 p(1x1)Cu(111) 的优化计算

发表于 2018-08-27 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

前面我们学习了Slab模型的单点计算，在此基础上，本节我们主要学习下怎么优化Slab模型。

1 修改Slab模型

前面单点计算的POSCAR不可以直接拿过来用吗？为什么要修改POSCAR文件？

这是因为在Slab模型中，要固定一部分原子来模拟体相，然后放开一部分原子来模拟表面。一般来说，表面的两层原子允许弛豫，而下面的那些原子则直接固定住。

这样做的物理意义在于：在真实的环境下，催化剂体相看做是不变的，只有表面的原子参与催化反应。

VASP中固定原子需要在POSCAR中进行操作。有两个关键点：

1）在POSCAR的第7行后添加一行，改行内容为Selective Dynamics，VASP只认第一个字母，你可以直接在这一行只加 S或者s。也可以换成其他的S开头的单词，比如SB，Sexy BigBro ， Super BigBro等等；
2）在原子坐标后面加上 F 或者 T 表示固定或者放开，因为坐标有xyz三个数值，因此我们需要三个F或者T。我们可以通过设置允许原子在某一方向上移动，而其他方向上固定。

F F F表示xyz全部固定；
T T T 表示xyz全部放开，
F F T 表示 xy方向固定，只允许原子在z方向上移动。

这一点我们前面也已经讲过了，详见Ex24 频率计算的输出与POSCAR原子的固定。

注意：

如果前面的POSCAR直接拿过来直接用，那么在优化过程中所有的原子都会被放开。回想一下我们计算Fe, Cu单胞的例子，计算中我们没有设置F或者T，晶胞尺寸和原子的坐标在优化过程前后都发生了变化。因此你不设置Selective 以及 F或者T，VASP默认是所有的原子都放开的。
Selective 和 F 或者 T同时出现，不要只设置Selective忘了FT，也不要只设置FT忘了Selective.

2 如何修改？

根据前面说，我们需要做的无非有两步，
1）在第7行后加入Selective Dynamics
2）在坐标后面加入T T T或者 F F F
我们在本节中，将底部的两层原子固定，上面两层放开。

实现本目的，有以下几种方式：

2.1）这个用文本编辑器就可以做到，很多编辑器都有列块模式，允许你同时添加N列到每行的最后。你可以使用下Notepad++，Atom这些免费开源的软件来操作

2.2）使用p4vasp，详见前面的练习Ex24

2.3）本节我们重温一下sed的用法，直接使用命令进行操作。

A) 在第 7行后添加一行：

1	sed -i '7aSelective Dynamics' POSCAR

7a 代表在第7行后面添加一行，a 是append的缩写。注意，a 和 Selective之间可以没有空格，也可以有空格，这里我们没有加空格。

B）添加了Selective这一行之后，原子坐标的行数也发生了变化，为10-13行。

C）第10，11行为底部的两层原子坐标，我们在最后面加上F F F

1	sed -i '10,11s/$/ F F F/g' POSCAR

注意： 第一个F前面有个空格。

d）在12-13行最后面加上 T T T

1	sed -i '12,13s/$/ T T T/g' POSCAR

这样，我们就顺利完成POSCAR的改造了，可以用来进行slab的优化计算了。

3 KPOINTS 和 POTCAR

保持原来的样子不变

4 INCAR

优化的时候，我门回想一下要设置的东西：
4.1IBRION 和 POTIM
IBRION可以使用 1 和2，我们这里用2。（1的话适合初始结构比较好的情况，这里我们也可以用。）
POTIM =0.1 官网默认值是0.5，个人感觉比较大， 0.1是一个很不错的选择。
4.2 EDIFFG = -0.01优化时以力作为收敛标准。注意 -0.01 已经很严格了，除非特殊需要，不要设置-0.001，这样你就会在QQ群里求助大家，为什么我优化的时候一直不收敛啊。如果设置-0.01，很多步过后仍然不收敛，可以尝试一下-0.02，也是一个很安全的选择。
4.3 NSW = 500 设置最大的优化步数。
4.4 ISIF = 2 slab 模型优化的时候，我们用ISIF =2 ，也就是VASP的默认值，不用设置，最终INCAR如下：

System = Cu111

ISMEAR = 1
SIGMA = 0.1
ALGO = FAST
ENCUT = 450
EDIFF = 1E-5

LDIPOL = .TRUE.
IDIPOL = 3

NWRITE = 0
LWAVE = .FALSE.
LCHARG = .FALSE.

IBRION = 2
POTIM = 0.1
EDIFFG = -0.01
NSW = 500

注意：
上一节我们提到Cu没有磁性，表述是不准确的，应该是纯净的Cu(111) slab模型中，Cu不需要考虑自旋。如果你算的是CuO(111)表面，则需要考虑自旋了。

5 提交任务

计算结束后，查看OUTCAR，OSZICAR的最后几行检查下是不是正常结束。

6 扩展练习

6.1 掌握本节提到的固定原子的三种做法：文本编辑器，p4vasp, sed，当然啦，更可以直接使用Vim 这个强大的工具来执行
6.2 查看结构优化相关的这几个参数，复习下它们的VASP官网说明。
6.3 运行本文的例子，然后将EDIFFG设置成 -0.02或者-0.001进行测试，查看与-0.01的区别。
6.4 VASP官网查找表面能计算的案例，结合本节以及前面的学习，尝试算一下Cu(111)面的表面能。

7 总结

截止到本节为止，Slab模型的优化，你已经可以初步掌握了。请务必跟着节奏认真练习。后面我们会继续学习Slab模型的优化，现在感觉陌生没上手的也不用担心。认真总结一下，我们如何从Bulk模型，然后一步一步操作到Slab模型优化的。每一步中都有哪些该注意的细节。这些都掌握了，就可以大大提高你计算的成功率，进而节省时间和精力。

找不同

批量搭结构

批量删除任务

找不同

扩展练习：

总结：

1 获取（下载）CONTCAR

2 使用ASE查看结构：

3 使用p4vasp 查看结构

4 其他软件：

5 扩展练习：

6 总结：

1 准备工作:

2 天河II号提交任务

3 查看结果

1）查看OSZICAR

2） 查看OUTCAR

4 删除VASP乱七八糟的输出文件

5 获取能量，计算吸附能

6 扩展练习

1 不同表面位点：

2 不同位点的吸附模型

3 使用p4vasp搭结构

3.1 从Top到Bri

3.2 从Bridge 到 Hcp：

3.3 从Hcp到Fcc

总结：

方法一：找规律+sed命令

方法二：分层处理+sed命令

方法三： 使用脚本分层

总结：

参考阅读：

扩展练习：

打个比方

搭建p-(2x2)-Cu(111)表面

搭建O在p-(2x2)-Cu(111)上的吸附结构

扩展训练：

总结：

1 什么是吸附能？

2 计算吸附能（单个O原子）

3 思考问题：

4 扩展练习：

5 总结：

1 搭建合理的O吸附模型

2 快速获取初始构型方法（二）

3 正式优化O的吸附

4 扩展练习：

5 小结

1 Top site 的构型

2 确定O原子在z方向上的坐标

3 搭建初算的模型(Cu-O双原子结构)

4 INCAR检查

5 POTCAR、KPOINTS检查

6 检查结果

7 总结:

1) 我们先看一下维基百科的解释:

2) IPUAC GoldBook中的定义:

3) 参考书

吸附在催化研究中的地位

吸附模型的搭建

小结

1 不同真空层厚度的测试(上一节内容)

2 结果处理

3 数据分析：

4 扩展练习：

5 总结：

1 INCAR中的参数

2 真空层的厚度及修改

3 批量处理POSCAR

4 扩展练习：

5 总结：

1 使用 vtotav-v5.2f 脚本

2 使用 vtotav.py 脚本

3 vtotav.py 脚本的获取和使用

4 扩展练习：

5 总结：

1 功函数的定义

2 VASP 计算功函数的过程

2.1 费米能级的计算：

2.2 静电势能的计算：

2）查看OUTCAR

方法三：使用脚本分层