Learn VASP The Hard Way

The 2nd edition

Ex43 p(1x1) Cu(111)的单点计算

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

上一节的内容主要简单介绍了一下Slab模型，具体的大家去看推荐的参考书！计算部分我们做了以下几个方面的操作：

Cu的Bulk计算，获取了稳定体相结构，晶胞参数；
将结果导入Material Studio，并切(111)面
将切出来的Cu(111) 面转换成POSCAR格式。

因此，VASP的四个输入文件，我们已经有了POSCAR和POTCAR，下面我们复习下POSCAR，INCAR和KPOINTS的准备工作，然后进行单点计算。

1 搭建slab模型

Cu(111) slab 模型（POSCAR）的制作流程（复习一下）。注意：前面我们用的是Conventional cell，下面用的是Primitive Cell。一般来说用Conventional cell， FCC的金属可以用Primitive Cell 但对于其他体系，通过Primitive Cell切出来表面模型有问题。

A）Bulk计算结束后，将CONTCAR用VESTA转换为cif文件，然后MS导入。
B）Build→Symmetry→primitive cell
C）Build→surfaces→cleave surface

D）Build→Crystals→Build Vaccum Slab 真空层选择15Å

E）导出为cif，通过VESTA转化为POSCAR格式，如下：

Cu\(1\1\1)
1.0
        2.5717000961         0.0000000000         0.0000000000
       -1.2858500481         2.2271576142         0.0000000000
        0.0000000000         0.0000000000        21.2994003296
   Cu
    4
Direct
     0.000000000         0.000000000         0.000000000
     0.333330005         0.666670024         0.098590001
     0.000000000         0.000000000         0.295760006
     0.666670024         0.333330005         0.197170004

前面的操作，我们生成 p(1x1) 的Cu(111)面，p是Primitive的简称，如果你不了解p(1x1)，去学习一下Wood Notation。当然，在此基础上，你也可以将p(1x1) 在晶胞a 和b方向上扩展，生成p(2x2)，p(3x3)， p(2x1)等不同大小的Cu(111)表面。本节，我们就用p(1x1)做为例子。

2 KPOINTS

POSCAR中，slab在平面方向上a 和 b 的大小为2.5717 Å，根据前面的经验规则，我们可以在两个方向上K点取值为：$13\times13$，z 方向上取值为1，最终KPOINTS如下图，脚本使用命令为：kpoints.sh 13 13 1 （不知道脚本的请往前翻）

qli@bigbro:~/test/cu$  kpoints.sh 13 13 1 
qli@bigbro:~/test/cu$  cat KPOINTS
K-POINTS
 0
Gamma
  13 13 1
  0 0 0

为什么在z方向上只用1个K点？查找参考书（DFT：A Practical Introduction， Page 88）的说明：
If the vacuum region is large enough， the electron density tails off to zero a short distance from the edge of the slab. It turns out that this means that accurate results are possible using just one k point in the b3 direction.

3 INCAR

在贴出来INCAR之前，大家首先回顾一下单点计算需要注意的几个方面：
1）System = Cu(111) 可有可无
2）Cu 是金属，可以使用ISMEAR = 1;SIGMA = 0.1
3）ALGO = FAST 或者使用默认值
4）纯净的Cu(111)体系中，Cu没有磁性，ISPIN 不用设置
注意：CuO中Cu具有磁性，如果你算的是Cu的氧化物，就需要考虑ISPIN了
5）单点计算：NSW = 0 也可以不设置，因为默认值就是0
6）截断能： ENCUT = 450 很多人喜欢测试这个参数，450直接用是一个很好的选择
7）EDIFF = 1E-5 电子步的收敛标准
8）偶极矫正：LDIPOL = .TRUE. ; IDIPOL = 3 Slab模型一般都需要加上这个参数
9）控制OUTCAR和其他文件的输出：NWRITE = 0 ; LWAVE = .FALSE.;LCHARG = .FALSE.
10）也可以加上NCORE这个参数来加快计算步骤。一般来说，不同的机子大家需要自己测试下NCORE参数对计算的影响，个人经验：如果你的机器一个节点有12个核，设置NCORE = 6；一个节点有24核，设置NCORE = 12。也就是节点中核数的一半。由于这个计算很快，本操作就不考虑了。最终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.

4 准备任务脚本，提交任务

计算完成后， grep ' without' OUTCAR （without前面有2个空格）或者grep E0 OSZICAR | tail-n 1，得到单点能为：-13.96918440 eV

5 扩展练习

5.1 计算的例子下载链接：链接: https://pan.baidu.com/s/1cHY-9g2Iswl5zV9bijBG7Q 提取码: ctxc
5.2 重复本节所有的练习
5.3 VASP官网查找本节中不了解的参数，然后学习
5.4 思考一下：我们得到的单点能量可以干吗用？

6 总结：

刚刚步入slab的计算，我们速度稍微慢一些，本节只讲单点计算，希望大家能够跟上步骤，亲自进行操作练习，如果出错了，请认真改正，直至得到和上传的结果一样为止。再次强调，大家平时一定要多浏览VASP官网，请把VASP手册下载放到桌面上，没事的时候就翻一翻。不要去网上找一些乱七八糟的教程，拿过来就直接照着练，官网才是最正宗的。

Ex42 slab模型

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

前面我们学习了Bulk计算的一些基本的计算，拟合以及优化晶格常数，DOS计算。在材料计算方面，还有能带，声子谱的计算等等，但由于本人方向不在这一块（主要是做表面反应机理研究的），就不能详细给大家解释了。如果有同学在这方面很擅长，可以把自己的计算心得写进本书，将不胜感激。说到表面反应的机理研究，我们首先要做的就是搭建模型。先理一下思路：

1）计算Bulk的结构，获取稳定的晶格常数
2）在此基础上，搭建Slab模型。

1 什么是Slab 模型？

这个问题，我们看一下参考书的第四章：Density Functional Theory: A Practical Introduction （注意：本书有中文版的，本人强烈不建议！请静下心来看原版的，也恳请不要在群里求中文的pdf版本！）

多相催化中，反应发生在催化剂表面上，也就是气（液）固两相的交界处。为了模拟固体的表面，我们需要一个在二维方向上具有周期性的结构来模拟固体的表面，第三维方向则不具有周期性，用来模拟气相或者液相。但是，更常用的则是在第三维方向上加上周期性，周期性的结构之间用一个什么都没有的空间（free space）来分割，也就是真空层。在这里我们需要注意的是要避免第三维方向上两个周期性结构的相互影响，也就是一个周期的电子密度沿着真空层方向要趋于零，这样才不会对相邻的结构产生影响。此时，我们需要做的就是：

1）增大真空层的厚度，厚度要适中，太小肯定不行，太大也不合适。想一想我们前面讲到的影响计算速度的一个因素，真空层太厚，也就意味着我们的模型尺寸变大，从而导致计算变慢。一般来说，对于表面反应的计算，15$\AA$的真空层厚度足够了。但是，对于功函数这一类对真空层敏感的计算来说，我们需要注意。

2）VASP中还需要加入：LDIPOL = .TRUE.和IDIPOL = 3(3指的是在z方向上)这两个参数来消除上下不对称的slab表面导致的偶极矩影响。https://cms.mpi.univie.ac.at/wiki/index.php/LDIPOL

对于表面结构，有以下几个需要注意的:
1）xy 方向上表面的大小;
A）这个影响表面吸附物种的覆盖度;
B）影响体系的尺寸大小和计算时间;
C）不同的大小需要选取对应的K点；回想一下我们前面提到的经验规则。

2）不同的晶面，(111)， (100)，(110);
A) 这取决于你研究的方向;
B) 不同晶面的表面能;
C) 不同晶面的表面结构，反应活性等。

3）表面结构的层数
A）层数多了，原子变多，体系在 z 方向的尺寸增加，也会影响计算速度；
B）计算中需要弛豫的层数；
C）不同层数对你要计算的性质会产生影响，比如表面能；
D）不同晶面需要的层数也不尽相同，一般开放的表面需要更多的层数；
E）根据自己感兴趣的性质，选择合适的层数，也就是需要你去测试一番。

4）Slab模型有两种，一种是上下表面对称的，一种是非对称的。对称性结构往往需要很多层，体系较大。非对称的结构体系较小，但存在偶极矩的影响，要注意加LDIPOL 和IDIPOL这两个参数来消除，后面我们会介绍到。

注意，为了更加顺利的进行下面的学习，参考书第4章的内容务必认真阅读。这会为你今后的表面研究打下坚实的基础

2 表面模型的搭建：

主要使用Material Studio进行切面操作，这一部分，本人不想做过多的解释，因为各大培训班，网上各种博客都有操作说明。本节大师兄通过Cu(111)面的例子给大家简单演示一下。

2.1 块体计算
1) 搭建Cu单胞模型，用VESTA导出VASP计算使用的POSCAR文件，这部分不清楚的筒子们可以复习前面bulk计算的相关内容。

2) INCAR文件的准备
A) Cu是金属，ISMEAR=1， SIGMA=0.1
B) ENCUT=700， Cu单胞含有的原子数很少，可以把ENCUT设置大一点
C) EDIFF和EDIFFG控制电子步收敛和离子弛豫的精度
D) PREC=High

SYSTEM = Cu-Bulk calculation RELAX
ISTART = 0
ICHARG = 2
ENCUT = 700
EDIFF = 1e-6
IBRION = 1
NSW = 200
EDIFFG = -0.01
ISIF = 3
POTIM = 0.1
ISMEAR = 1
SIGMA = 0.1
PREC = H
LREAL = Auto
ALGO = Fast
NCORE = 4

E）这里LREAL设置为Auto，是为了和后面表面能计算统一。计算时要根据原子数的多少，以及后续的计算选择LREAL的值。不能拿LREAL=.FALSE.和LREAL=ON/.TRUE.计算的结果进行能量比较。看下面

F) 注意：块体材料计算时，为了后续计算其他性质，常常要求整个计算中使用相同的ENCUT，ENAUG， PREC， LREAL， ROPT，这部分大家可以看看手册 8.3 What kind of “technical” errors do exist， overview

3) KPOINTS
a=b=c=3.6147 $\AA$，根据前面的经验，这里设置为$13\times13\times13$的KPOINTS(Gamma Centered).

4) 准备POTCAR文件进行计算：见前面讲的相关脚本—Ex30 进阶前的准备（三）

5) 最后计算完的Cu的晶格常数为 3.6370 $\AA$，实验值为3.615 $\AA$。大家可以去计算下。（计算例子可以QQ群中下载，也可通过该链接：http://pan.baidu.com/s/1dE84fhn 密码：nvwa）

2.2 结构优化的注意事项：

关于结构优化，大家在计算时要注意下面的部分，7.6.2 Accurate bulk relaxations with internal parameters (one)一定要详细的读。

2.3 小结
主要以金属Cu为例讲了块体材料的计算，以及块体材料计算时的一些注意事项：

体材料计算时首先要对要了解材料的性质，根据材料是金属、非金属，还是半导体，选择不同的ISMEAR和SIGMA；
计算时要根据原子数的多少，以及后续的计算选择LREAL的值，这方面不熟悉的朋友可以看看手册；
精确进行块体计算的过程以及注意事项；
回顾Ex36 Bulk的优化。

3 MS建立slab模型

3.1 具体步骤

VESTA将CONTCAR转成cif格式，
MS导入cif文件，
build→surface→cleave surface，如下图所示，在cleave plane(hkl)可以输入要切的晶面，如（111）面，
在Surface Mesh→surface vectors可以改变slab模型的vectors，之后选择Cleave进行切面，本例为（111）。这个时候也可以通过调节Surface Box→top和thickness，选择slab模型的表面终端和厚度。
确定好表面终端和slab厚度之后，加真空层。
Build→Crystals→Build vacuum slab，Vacuum orientation选择C指定真空层加在z方向上。可以通过改变Vacuum thickness设置真空层的厚度10 $\AA$，15 $\AA$或者更大。这可以参照之前的计算。最后选择Build完成slab模型的建立。
可以通过Lattice parameters→Advanced→Re-orient to standard把真空层设置沿Z轴。

4 总结：

从本节开始，大家开始步入表面相关的计算学习。请务必认真阅读参考书，重复计算示例。本节主要简单介绍了一下Slab模型以及基本的注意事项。通过Cu的计算复习了一下Bulk的计算以及展示了一下通过Material Studio进行Cu(111)表面的搭建过程。

Ex-41 DOS 计算(五)

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

前面我们讲解了手动求积分的原理，本节由QQ群中的小书童给大家演示一遍Origin求积分的操作，加深一下对积分原理的印象，也借此机会掌握另一项Origin的技能。Origin的 .opj文件和Excel数据已经打包上传到QQ群中，也可以百度网盘下载：http://pan.baidu.com/s/1bpvmXY3

1 求d带中心

1.1 复制粘贴能量、态密度、能量*态密度这三列数据到Origin中，数据见前一节练习：

下面我们取了d band 的一小部分区间作为示例进行讲解，大家学习完之后自己主动练习，获取最终的d band center.

1.2 选中数据并作图：

plot —> Line —> Line

1.3 求积分，先选中一条曲线，右键，Ungroup（将曲线从组中分离，便于选择）。

1.4 打开积分对话框：

1.5 再选择一条曲线进行积分，下图选择的是黑线。

点击OK，截图中的area就是我们要的积分，这一面积对应下图公式中的分母部分。

d band center 的积分区间为整个d带：即从负无穷到正无穷。本例我们用区间 [-6, 2]作为例子！有时候大家会对积分的区间产生疑问。一般来说积分区间是从负无穷到正无穷，但是：

1）如上图，[-9,-8]这个区间Density为0，积分从[-9,4.0] 和[-8,4.0]结果没什么区别。
2）如果我们关心的是[-6,4.0]这一部分的性质，倘若4.0以后的部分会产生很大的噪音，那么我们也可以只从-6积分到4.0。
所以，积分区间也不是很绝对，大家要根据自己的体系合理选择。

1.6 同理，选中另一条红线（能量和密度乘积），进行积分计算（上图中的分子部分）

结果如下:

最终，我们计算的d带区间[-6,2]的中心果为：-34.035/ 17.428 = -1.953

2 求电子数目

电子数目的积分是从负无穷到费米能级。本例中是从-10到0（费米能级）进行积分，选择黑线，重复前面的操作：

计算结果如下，电子数为13.00

3 扩展练习：

3.1 熟练Get本节Origin求积分的技巧；

3.2 阅读前面几节DOS计算的练习章节；熟练掌握DOS的提取，导出，作图和积分；

3.3 d band center 是对整个区间求积分，电子数目的积分区间是从负无穷到费米能级。

4 总结：

截止到本节，DOS的计算就暂时告一段落了，实际计算中，DOS的分析远远不止我们讲述的这么简单，而且网上DOS分析的脚本很多。大师兄了解到的有VTST的DOS处理脚本和vaspkit等，大家可以尝试练习使用下。如果你在科研工作中对DOS分析有自己的心得，DOS分析脚本的使用经验等，欢迎补充，分享给大家。

Ex-40 DOS 计算(四)

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

DOS积分的原理，积分公式相信大学的时候，大家已经掌握了，现在还有可能忘记了。不过没关系，虽然本节将积分，高等数学里面的暂时还用不到。对于DOS图来说，通过求积分可以得到相应的电子数。积分也就是求DOS曲线下的面积。但DOS图并不能通过一个函数的积分来获得，因为所对应的区域毫无规则可言。那么我们该怎么样积分呢？

1 积分的原理

积分的原理就是把区域分成N个小块，计算每一块的面积，然后求和。如下图：

大师兄…..

好吧，言归正传，大师兄亲自画了一幅图，讲解积分的原理，如下：

1)将每一个小区间近似看成一个小的梯形，

2) 梯形的上底和下底为：$y_2$ 和$y_1$，高为：$x_2-x_1$。

3, 4) 也可以看成一个高为：$（y_1+y_2）/2$，底为$x_2-x_1$ 的长方形。

5）面积为两者的乘积：$ΔA =（x_2-x_1）* （y_1+y_2）/2$。

6）将所有的小区间面积求和便得到了曲线所覆盖的面积，也就得到了我们所想要的电子数：$A =\sum\Delta A$

注意：

1）电子数的积分区间是从负无穷到费米能级；

2）p4vasp的图中的0点就是费米能级，已经扣除过了。

3）科学网的一个参考链接：<<http://blog.sciencenet.cn/blog-567091-736149.html/>>

注意：图中的DOS，如果你计算了，发现和大师兄的不一样，不用担心，大师兄本人把这个例子再运行了一遍，跟上面的图也不一样。（原因正在调查中）

2 d-band center 的公式

在这里：

ε为前面的我们提到的x， $n_d$(ε)为y。在DOS图中，ε为能量， $n_d$(ε)为电子的密度。

注意：积分区间是从 负无穷 到 正无穷

参考书籍： Fundamental Concepts in Heterogeneous Catalysis，Jens K. Nørskov, Felix Studt, Frank Abild-Pedersen, Thomas Bligaard， John Wiley & Sons, Inc.

QQ群文件的Books文件夹中下载

百度网盘下载链接：http://pan.baidu.com/s/1bWS2Gu

d 带穿过了费米能级，说明Ru是金属。（废话，等于白说！）

d带的积分从最左面到最右面。

关于d band 的介绍：知乎有个专栏，有兴趣的可以浏览一下：

https://zhuanlan.zhihu.com/p/26230183

3 Excel手动求积分：

对应的公式都写在第二行里面了，excel中计算完第一个数值，直接选中然后往下拽就可以得到后面的数据了。计算文件以及Excel表格已经上传到QQ群和百度网盘，大家可以下载仔细观看练习，掌握该技能。

计算下载链接：http://pan.baidu.com/s/1kVqWFQn

Step1: p4vasp 的 d band 数据！

Step2：A,B两列的乘积，也就是能量和密度的乘积：对应d band center公式中的分子：

Step3：

1）A4-A3: 小区间的能量间隔，也就是 dε，或者 xdx，$\Delta$x， $x_2-x_1$;

2）B4+B3: 小区间的两个起点和终点的高度：也就是对应的小区间两端的能量密度：$y_2+y_1$

3）C4+C3：小区间两端能量和密度乘积的加和： $x_1y_1 + x_2y_2 $

4）Step4：求小区间的面积：

$D4E40.5$：$(x_2-x_1) * (y_2+y_1)/2$

$D4F40.5$: $ (x_2-x_1) (x_1y_1 + x_2*y_2)/2$

5）Step5: 将所有小区间的面积求和：

电子数：$\sum( (x_2-x_1) * (y_2+y_1)/2 ) $，从G4到G137(G137对应的是费米能级的那个点！！！)求和。

再次强调：电子数目是从负无穷积分到费米能级！

d带中心公式中的分子：从H4到H303求和

d带中心公式中的分母：和电子数的积分一样，但积分区间不同，这里从G4到G303求加和。

再次强调：d带中心积分区间，按照公式，从负无穷积分到正无穷！

d带中心公式中的分子：从H4到H303求和

d带中心公式中的分母：和电子数的积分一样，但积分区间不同，这里从G4到G303求加和。

以上将公式分开写，是为了方便大家理解，实际操作中，直接使用导出来的数据，将所有的公式合并在一起就可以了。

4 扩展练习：

4.1 亲自上手操作，重现本节的结果；

4.2 尝试使用其他软件，方法求积分；

4.3 有好程序的，好脚本的分享给大家。

5 总结：

学习完本节，大家对DOS积分的数学知识应该可以掌握了，后面不管用什么软件或者脚本，都是将这个积分的原理写成机器语言来执行。知道了原理和通过excel自己操作一遍之后，基本上DOS积分的东西就懂了。下一节，QQ群里的小书童给大家演示一下Origin的积分操作。

Ex-39 DOS 计算(三)

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

本节主要介绍如何使用p4vasp获取DOS数据。

1 计算实例：

计算的体系是一个Ru的单胞，含有2个Ru原子。结构优化完毕后进行以下操作。

1.1 cp CONTCAR POSCAR

1.2 INCAR的设置（一步计算）：

IBRION = -1
NSW = 0
ISMEAR = -5
LORBIT = 11

1.3 KPOINTS 采用: 21 21 21 （gamma centered：Ru为六方堆积结构！！！）

1.4 提交任务等待结束。

1.5 该计算的所有文件可以通过百度网盘下载： http://pan.baidu.com/s/1kVebDWz

2 DOS可视化：

2.1 打开p4vasp

2.2 点击左侧open，导入vasprun.xml文件；

注意：前面的2步，Linux系统下面直接输入命令：

p4v &或者 p4v vasprun.xm即可。

2.3 查看DOS：Electronic –> DOS+bands

效果如下如：

2.4 获取d轨道的DOS:

这里我们取d band作为例子，大家根据自己的需要选择：

Electronic –> Local DOS+bands control

点击后如下图：

注意：

1）这里我们填的是Ru，表示选择所有的Ru原子，具体到大家自己的计算中，

换成自己感兴趣的原子即可，怎么选择，我们前面已经讲到过了。（通过左侧的show按钮，以及鼠标和空格选择）

2）这里我们没有考虑自旋，如果你的体系有自旋，可以选择up，down或者both。up和down的数据导出来的时候通过一个空行分割（后面会讲到）。

3）选择感兴趣的轨道：

如果计算中你设置的LORBIT=10，点击 dxy的时候，则其他的d轨道也会同时被选中，

如果你设置的LORBIT=11，点击dxy时，则仅仅选择dxy，其他的d轨道不会选中。如下图：

如果看不明白：

a）VASP官网查看 LORBIT = 10 和11 的含义；

b）亲自动手操作一下就明白了。

选择完毕后，点击 Add New line，效果如下：

图中增加了一条红色的曲线，为Ru的d轨道的DOS图，如果你选择了s，点击Add New line，就会增加另外一条曲线。不想要某个曲线了，可以点击remove line，大家自己随意操作，直至熟练为止！！

3 后处理：

3.1 导出数据：

p4vasp导出数据通过中间的一个按钮：（注意，这里不通过左上角的File按钮！），如下图：点击 Graph –> Export

弹出的窗口中，选择保存的目录和数据的文件名称

注意：

1）choose 选择保存的目录；

2）我们选择的D盘下的VASP目录；

3）填上保存数据的文件名：

a）文件名处 .dat 记得写上，要不然后面可能会导不出数据！

b）为了保证导出数据，Export 多点几遍！！！

c）去导出的目录里面查看一下，如有了数据的文件名，说明保存成功，没有的话，继续点Export.

图中说明：导出数据成功。

d）为了避免浪费时间，一定要保存成功后再关闭窗口，否则你还要重复一遍前面的步骤。

4）保存成功后，点击OK，关闭窗口。

3.2 查看数据：（Notepad++直接打开，无视.dat尾缀！）

问题1：为什么在301行后面有一个空行？（重点是301）

1) 这是因为DOS计算的点VASP默认是301个，该数目由NEDOS这个参数来确定。NEDOS决定了DOS曲线中间隔的点数，如果你想要自己的DOS曲线更加平滑，可以设置一个较大的NEDOS数值，比如：NEDOS = 3000

2) 如果你设置了NEDOS = 500，那么会在500行以后出现一个空行。

问题2：为什么会有一个空行？？？(一定要搞明白！！！)

我们看一下前面DOS的图：

图中有2个DOS曲线，黑色的为空行前面的部分，红色为后面的部分。

同样，如果DOS图中有N个曲线，则导出的数据中会有N-1个空行将它们分开。

N个曲线的横坐标坐标（也就是能量）是完全一样的。知道了这一点，你就知道了自己想要的数据在哪里了！！！

此外，如果有磁性的话，你要注意了，up和down也是通过空行分开的。

3.3 作图：

有了数据之后，就可以随心所欲地作图了。将数据导入excel或者origin作图，也可以使用其他小脚本，程序。这里我们使用Excel简单说下：

左侧为空行前面的数据，Ru总的DOS， DOS曲线1 （黑色的那个）

右侧为空行后面的数据，Ru d轨道的DOS， DOS曲线2（红色的那个）

作出来的图跟p4vasp的结果一样。图像就不美化了。关键是有了数据，并且知道数据里面空行前后代表的内容。

4 扩展练习：

4.1 下载本节的练习，可以在自己的服务器上提交，也可以直接使用结果；

4.2 通过p4vasp读取vasprun.xml文件，熟练进行本节的所有操作；

4.3 不要偷懒！！！

5 总结：

本节主要手把手交给大家如何用p4vasp这个软件查看DOS计算的结果。希望大家掌握这几点：

1 会熟练使用改软件进行本节的操作；
2 了解导出数据的格式（空行是怎么分割不同曲线的）；
3 知道LORBIT = 10 和 11 的区别，以及在可视化操作中的不同；
4 最关键的：要会分析DOS！！！

Ex-38 DOS 计算(二)

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

前面一节，我们留了一个问题：DOS为什么要算两步？

答：可以一步直接算！也就是前面一节的问题是瞎问的。那为什么大家总是有2步计算的疑问呢？本节，我们就解释这个。

为什么算两步？

官网说了，算DOS有2个方法，一个是直接进行self-consistent计算（大家常说的自洽或者静态计算。）然后处理DOSCAR和vasprun.xml文件即可，如图：

我们看一下这么做的前提：

1）由于高质量的DOS需要精细的K点，如果我们设置的K点很多，就会造成计算上的负担，前面我们讲过K点与计算时间的关系；而很久很久以前，计算能力并不如现在这样快，因此可以通过分步计算来解决这一问题。
2）另外一个原因就设计到能带的计算了，这里我们摘抄一下网上的解释：由于在能带计算时k点是一些在倒空间高对称线上的点，不能进行自洽计算。参考网址：https://http://blog.sciencenet.cn/blog-567091-675253.html/ 也就是计算能带的时候，自洽计算是必须的一步；
3）即使增加K点的数目，电荷密度和有效势能的收敛依然很快，也就是K点的变化对电荷密度的收敛影响不大。

分析：

能带计算我们暂不考虑，综合下1 和 3 ，在结构优化完成之后，我们可以这么算DOS：

第一步，用小的K点算个单点，生成CHGCAR文件；
第二步，读取上一步的CHGCAR文件（ICHARG=11）。

这样做就避免了直接用高K点网格所导致的计算负担。对于DOS计算的两个步骤，归根结底是节约时间的问题。因此，算2步并不是必须的！！！如果够土够豪，直接用高密度的K点，一步计算，没毛病！ 但是对于能带计算， 则必须算2步。

此外，VASP的说明书已经很古老了，以现在的计算能力，直接使用大K点一步计算，一般来说都可以承受的。所以，当你知道了为什么要算2步的时候，再去浏览网上的相关经验帖子，就很容易知道是怎么回事了。

LDOS 和 PDOS

LORBIT = 10 把态密度分解到每个原子以及原子的spd轨道上面，称为为局域态密度，Local DOS (LDOS)

LORBIT =11 在10的基础上，还进一步分解到px，py，pz等轨道上，称为投影态密度（Projected DOS）或者分波态密度(Partial DOS)，即PDOS。

所以LORBIT = 11可以提供我们更多的信息。

WAVECAR

那么，WAVECAR读不读呢？大师兄的观点是：有则读，无则不读。对于WAVECAR的读取，我们需要了解ISTART这个参数：

如果前面计算中保存了WAVECAR，且ISTART没有设置，VASP默认是读取的。
如果没有WAVECAR，即使你设置了ISTART=1或者2，VASP找不到可以读取的WAVECAR，也不会报错，而是继续算。

那么怎么控制WAVECAR的输出呢？

1)通过设置LWAVE这个参数

2)读取WAVECAR可以极大地减少自洽的时间，但是VASP的WAVECAR非常大，上百M或者几G都是很常见的。一不留神存储空间就被占满了。所以，在读取WAVECAR的时候，一定要确定自己的存储空间。
3)如果前面计算步骤中（优化的过程）保存了WAVCAR，那么后面DOS计算的时候（1步计算或者2步计算均可），都可以读取，这会加快计算速度。

扩展阅读

4.1 阅读DOS和能带计算的VASP官方手册；

4.2 了解DOSCAR的内容以及各行各列所代表的含义；

总结

结构优化完毕后：

一步计算DOS必须的参数：

1 ISMEAR = -5
2 LORBIT = 11
3 高密度的K点

两步计算DOS必须的参数:

第一步：

ISMEAR = -5
LCHARG = .TRUE.
稍微低密度的K点

第二步：

ISMEAR = -5
ICHARGE = 11
LORBIT = 11
高密度的K点

如果结构优化的时候，存了WAVECAR，计算DOS的时候可以读取WAVECAR，直接一步计算搞定。

Ex-37 DOS 计算(一)

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

使用VASP计算，很多时候都逃不掉DOS，能带计算的相关问题，尤其是对于计算材料的童鞋们，更是家常便饭一般。群里很多人，很多新手们都时常在讨论DOS的计算。这里我们通过VASP官网的说明，解释一下算DOS的具体步骤。前面我们学会了如何拟合或者优化稳定的晶胞结构。在此基础上，我们可以计算一下相关的DOS信息。

1 KPOINTS

1.1 K点数目

与结构优化相比，算DOS的时候，需要用到更多的K点数目，这是因为K点越多，画出来的DOS图质量越高。

引用官网的话：

1	A high quality DOS requires usually very fine k-meshes.

1.2 K点数目的选取

K点数目越多越好，我们该如何设置K点数目呢？

还记的前面我们讲到的K点选择的经验规则吗？那一个规则可以认为是我们平时计算时K点选择的标配。对于DOS计算，我们就需要把配置提高一个档次了。

一般来说，K * a = 45左右之间完全可以满足你的要求，大伙可以根据这个经验来选择K点。

2 NEDOS

NEDOS这个参数在DOS图的质量上面也有着很重要的作用。比如我们的DOS能量区间范围（DOS图的横坐标）为：[-10 eV，10eV]，VASP默认的将这个能量范围分成301点，然后作图。301也就是默认的NEDOS的取值。如果我们设置的NEDOS值够大，那么DOS区间就会被区分地越精确。NEDOS的取值一般来说：

NEDOS = 3000左右就足够好了。太大也没什么意义；
NEDOS越大，VASP输出的DOSCAR，vasprun.xml文件也就越大，占用存储空间。
经常有人抱怨说自己的DOS图有很多尖锐的峰，可以尝试着通过增加NEDOS这个办法来解决。
更多的信息，自己参考一下：https://cms.mpi.univie.ac.at/wiki/index.php/NEDOS

3 ISMEAR（一）

首先我们看官网的话： https://cms.mpi.univie.ac.at/wiki/index.php/ISMEAR

For the calculation of the total energy in bulk materials we recommend thetetrahedron method with Blöchl corrections (ISMEAR=-5). This method also givesa smooth nice electronic density of states (DOS).

也就是说 ISMEAR = -5 的时候（Blöchl修正的四面体方法），我们可以得到一个非常平滑的DOS图。

注意：

3.1 K点数目

设置ISMEAR = -5 的时候，如果K点数目K点的数目小于等于4 ，计算会出错，得到如下的错误结果：

VERY BAD NEWS! internal error in subroutineIBZKPT:

Tetrahedron method fails for NKPT<4. NKPT= 1

这也是很多人常见的错误。官网说的是K点数目小于三：the tetrahedron method is not applicable， if less than three k-points are used. （QQ群的恒驰一强发现官网的这个错误。）

K点不够，用ISMEAR = -5出错的解决办法：

既然Ｋ点不够，那么我就增加Ｋ点，然后再使用ISMEAR= -5　（简单粗暴，强烈推荐使用）
如果增加了K点，可能还是会出错。有时也会出现下面的错误（微信群的群友（Cu—Ni）提供）。我们先把解决方法列出来，错误部分大家慢慢看：
- 直接换一下ISMEAR的取值。
- 群友还发现：在保证K点数目大于4的时候，有时候减少K点数目或者增加K点数目都可以解决这个问题。如果你的服务器还算可以，建议增加K点数目，毕竟和K点数目越多，DOS的质量越高。这个办法大家可以参考一下。
1
2
WARNING: DENTET:can't reach specified precision
Number of Electronsis NELECT =

官方的解释：

http://cms.mpi.univie.ac.at/vasp-forum/viewtopic.php?t=416

http://www.error.wiki/The_old_and_the_new_charge_density_differ

出现此警告(DENTET)的原因是因为无法通过tetrahedron方法得到足够精确的费米能级。也就是将态密度积分到费米面的电子数和体系的价电子数目不一致。

3.2 适用体系：

ISMEAR = -5 适用于所有体系的DOS计算。非DOS计算的时候：
对于金属体系来说，结构优化的时候不能使用ISMEAR= -5（注意：是优化结构的时候不能用！），这是因为四面体方法不能很好地处理费米能级处的电子占据情况，导致算出来的力会有一定百分比的误差。所以，对金属体系结构优化的时候，ISMEAR > = 0 参考：https://cms.mpi.univie.ac.at/vasp/vasp/Partial_occupancies_different_methods.html
对于半导体和绝缘体的体系，ISMEAR>0 是不可以的。只能ISMEAR<=0。
四面体方法（ISMEAR = -5）不适合计算能带（对所有的体系来说的）。多谢wuli8老师帮忙完善！
使用ISMEAR= -5 的时候，SIGMA的取值没有影响，如果不放心，取SIGMA = 0.01。

3.3 小结：

算DOS，只要K点不少于3；用ISMEAR=-5是个不错的选择。

4 ISMEAR（二）

如果体系很大，只能适用gamma点来算，ISMEAR = -5的时候，肯定会出错，但服务器不给力，不能增加Ｋ点的时候，怎么办？

对于所有的体系（K点数目小于4也可以）：可以使用ISMEAR = 0；SIGMA = 0.01。对于大部分的体系都能得到理想的结果。原则上来说，使用GS方法的时候(ISMEAR=0)，SIGMA的数值要测试下，保证entropy T*S这一项平均到每个原子上小于0.001 eV也就是1meV。不想测试的话，直接用个很小的值，比如这里我们说的：SIGMA = 0.01。
对于金属体系来说，也可以使用ISMEAR = 1；SIGMA = 0.01。SIGMA取值太大，计算出来的能量可能不正确；SIGMA取值越小，计算越精确，需要的时间也就越多。值得注意的是：这里我们使用的0.01已经很小，没有必要设置的再小。因为对于金属体系，使用MP方法(ISMEAR=1..N)时，SIGMA= 0.10 差不多就足够了。官网给的参考值是0.20。

http://cms.mpi.univie.ac.at/vasp/guide/node124.html

在KPOINTS确定之后，使用多大的SIGMA值，大家最好测试一下。原则如下：SIGMA取值在保证OUTCAR中entropy T*S这项的能量平均到每个原子上小于 1 meV的前提下：尽可能地大。这样可以在保证准确度的同时，也加快收敛速度。

记住： VASP学习最快的途径就是不停地看官网，然后亲自上手去测试，测试，测试！并观察分析结果。

4 扩展练习:

4.1 阅读VASP官网关于ISMEAR和SIGMA的所有说明：

4.2 下载VASP的pdf说明书，搜索书中所有的ISMEAR和SIGMA关键词，阅读所有相关的内容；

4.3 思考SMEAR方法的意义？SIGMA的意义？

4.4 查看VASP说明书，查阅相关文献，了解MP和GS方法

4.5 测试MP和 GS方法中，SIGMA取值对计算时间，能量，收敛步数的影响。

4.6 分析下为什么算DOS的时候，要算两步： selfconsistent 和 none-selfconsistent calculations?

5 总结:

看完本节：你应该知道计算DOS的时候：

KPOINTS和NEDOS设置的一些内容。
ISMEAR要用-5。
KPOINTS因计算硬件限制不能设置的很大，数目小于4的时候：
- 对于所有体系均可以使用ISMEAR=0。
- 金属体系还可以用ISMEAR=1..N，官网建议SIGMA为0.20，太小的SIGMA值对收敛会产生影响。使用0.01-0.10的数值都是很安全的选择。
- SIGMA的数值需要测试一下，一般来说在0.01-0.05之间足够了。
非DOS计算的时候，对于金属来说ISMEAR不能等于 -5，优先使用ISMEAR= 1。非金属来说（半导体和绝缘体），不能 > 0 。对于所有的体系， ISMEAR= 0 则是一个很安全的选择，但SIGMA的数值要测试一下。说了这么多废话，还是官网简单明了：

For further considerations on the choice for the smearing method see sections 9.4， 10.6. To summarize， use the following guidelines:

For semiconductors or insulators use always tetrahedron method (ISMEAR=-5)， if the cell is too large (or if you use only 1 or two k-points) use ISMEAR=0.

For relaxations in metals always use ISMEAR=1 and an appropriated SIGMA value (the entropy term should less than 1 meV per atom). Mind: Avoid to use ISMEAR>0 for semiconductors and insulators， it might result in problems.

For metals a sensible value is usually SIGMA= 0.2 (that’s the value we use for most transition metal surfaces).

For the DOS and very accurate total energy calculations (no relaxation in metals) use the tetrahedron method (ISMEAR=-5).

Ex36 晶格参数的确定（直接优化晶格常数）

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

前面一节我们学习了如何通过提交一系列不同晶格常数的单点计算以及拟合BM方程获取稳定结构的方法，此外，我们还学习了BM方程拟合的一个python脚本以及分析了脚本的工作流程。

师兄，算个晶格常数这也太麻烦了吧？有没有简单点的，一步到位的方法？

有！！！前面我们提到说可以通过两个方法： BM方程拟合是一个，现在大家已经掌握了。今天我们学习第二个方法。

1 2	ISIF = 3 ENCUT = XXX

1 ISIF 参数

在学习ISIF = 3 前，我们先总结一下，之前不论是单点计算，结构优化，还是频率分析，我们都没有讲过这个参数。因此使用的是默认值。

查找VASP官网说明书，默认值是ISIF= 2 。是的，前面我们都是用的 ISIF = 2 。（废话，等于白说！）你还会发现有一堆其他的ISIF参数，从1 - 7，看的眼花缭乱。

坚持一下，你会发现它们有各自的特点，原子坐标是否变化，晶胞形状是否改变，以及晶胞体系是够变化。不同的参数大家根据自己研究体系的性质进行选取。

本文主要介绍ISIF = 3 时对晶胞的优化。其他的取值比如加压计算等，为避免对大家造成误解，在这里不做进一步地说明。因为本人计算的时候，只用到过2 和 3 的情况。

ISIF = 3 的时候，晶胞中原子的坐标，晶胞形状，以及体系都随着优化的过程发生变化。回到Fe单胞的计算，我们先看一下INCAR的输入，然后再详细解释其中的内容。

对比下，前面单点计算的INCAR，你会发现有几个不同：

1） ENCUT = 600 设置的比之前的450大。

2） EDIFFG = -0.015 结构优化的收敛标准

3） NSW = 100 结构优化，肯定不是一步！

4） ISIF = 3 （讲的就是ISF =3 ，肯定要加上）

5）此外， ISMEAR，SIGMA, ISPIN以及 MAGMOM大家自己主动回顾一下

PS：很多人在单点计算的时候，除了NSW = 0，还画蛇添足地设置了ISIF = 0。虽然对计算结果没有影响，但我还是想问一句：大哥，脑子去哪里了？单点计算就是原子不动弹，设置NSW = 0 或者 1，或者 IBRION = -1，其他的按照之前的即可。

POSCAR，KPOINTS，POTCAR直接将前面单点计算目录（1.00）下的复制过来。

2 ENCUT 取值

师兄，Fe的ENCUT 默认值是：ENMAX = 267.882 eV，前面单点计算设置的 450 eV。这里为什么设置成600 eV ？

如果你仔细看VASP官网的说明书，在底部有这么一段话：

红色的链接如下：

http://cms.mpi.univie.ac.at/vasp/vasp/Volume_vs_energy_volume_relaxations_Pulay_Stress.html

最后一句，如果计算时体积发生了变化，我们需要增加ENCUT的值，比如说：ENCUT = 1.3 * max(ENMAX) ， max(ENMAX) 的意思是，如果有N个元素，取最大元素的ENMAX值。

这么做的原因是为了尽可能地消除 Pulay stress 对计算的影响。那么什么是Pulay Stress 呢？

3 Pulay Stress

通过查看红色的链接内容，字面上的理解就是，当晶胞体积发生变化时，由于计算时的平面波基组不全（设置的ENCUT不够大）使得计算时的收敛不是绝对的，进而导致应力张量不对，最终导致体系的形状，体积以及能量不正确。（这个我不是很懂，字面上直译过来的。）维基百科的描述如下，链接：https://en.wikipedia.org/wiki/Pulay_Stress

A plane wave basis set is created for thehexagonal lattice (left), using the reciprocal lattice vectors inside the redcircle. Then the lattice relaxes into a cubic symmetry (right). Keeping the redcircle basis constant results in lattice vectors taken from an ellipsoidinstead of a spherical area (compare to the blue circle).

深层次的俺也不懂，但是PulayStess 有2个办法可以减小：

A) 体积不变，采用一系列的晶格常数计算，然后拟合得到准确的晶格常数；（你已经知道怎么去做了，这么做是为了避免计算时晶格体积发生变化！）；

注意：

1）在LVASPTHW的拟合流程中，拟合出来的是BM方程， BM方程的最高次数是3，也就是x的三次方。此时：x = (1/a$^2$) ，具体见Ex32-35的详细说明。

2) 如果直接用晶格常数去拟合，则是6次方的关系。

3) 如果你用体积去拟合BM方程的话，最高次数是2：

B) 既然是因为ENCUT不够导致的，那么我就使劲增加ENCUT来消除（上图中unless后面的那句）。借用维基百科的一句话：Pulay stress can be reduced by increasing the energy cutoff. 所以，当有人用ISIF= 3计算，而不去设置ENCUT的时候，尽情地去嘲笑他傻X吧！！！感觉不过瘾，还可以继续吓唬他说后面的计算也有问题。

4 ENCUT 取值

官网说设置到130% * max(ENMAX)，一般来说就可以了。但很多人还是不放心，那我们就继续疯狂一些，设置到600 或者700。本人一般习惯用600。

师兄，前面你不是讲到增加ENCUT会增加计算量么？而且要保持计算中使用同一个计算的参数？

是的，这些你说的都是正确的。但是，我们使用ISIF = 3 的目的是获取稳定的晶格常数，但为了避免Pulay stress导致的误差，就必须要增加ENCUT（No Choice！！！）。一旦我们计算完成晶格常数的计算后，可以在此结果基础上，统一使用其他的ENCUT值进行计算。也就是只有确定晶格常数的这一步是个例外，其他情况还是用原来的值。

举个不恰当的例子，这就好比当初你追女朋友的时候，为讨人家欢心，疯狂地买买买，什么要求都答应。等追到手了，结婚了，新鲜劲过去了，就原形毕露，不再这么大手大脚地去疯狂了。

讲到这里，使用ISIF = 3 时，你脑子第一反应就是疯狂地增加ENCUT，600不够，上700， 700感觉不够，上800，在计算承受的范围之内，越大越好！因此ISIF = 3 和 ENCUT两个参数必须同时出现。

5 计算结果

设置完后，提交任务：计算出来的Fe晶格参数为： 2.8318 Å，对比下前面拟合出来的结果： 2.8332 Å，两者相差：0.0014 / 2.8332 = 0.019 %，这么点误差，完全可以忽略掉了。如果不放心，将ISIF =3 的结果，使用前面的设置，算个单点，然后跟BM拟合结果对比下，两者能量相差甚微，完全可以忽略。

Methods	Lattice Constant	Energy
BM拟合	2.8332 $\AA$	-16.47179683 eV
直接优化（ISIF=3，ENCUT=600）	2.8318 $\AA $	-16.47181669 eV

6 总结：

学习完本节，怎么计算晶格常数的两个方法你应该都掌握了，希望大家能学到的东西有以下几个：

1 什么是BM方程？

2 BM方程拟合的方法（最小二乘法，脚本以及运行）；

3 BM方程获取晶格常数的原理；

4 ISIF = 3 和 ENCUT 必须同时设置；

5 Pulay Stress 是什么，如何避免；

6 （ISPIN，MAGMOM）：什么时候加，怎么设置；

7 （ISMEAR，SIMGA）：什么体系用什么值，怎么设置；

8 （IBRION，POTIM，NSW）：根据初始结构好坏进行设置；

9 （EDIFF，EDIFFG）：怎么决定收敛的；

10 KPOINTS的取值：经验方法（k * a），以及脚本；

11 POTCAR的生成：手动方法以及脚本。

12 以及怎么追女朋友，如何防止被追？

Ex35 晶格参数的确定（Birch-Murnaghan状态方程-3）

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

通过上一节（Ex34）的脚本，我们可以拟合出稳定晶体结构的晶格参数，但是，这个脚本是怎么运行的呢？大家在浏览之后，有没有什么疑问？本节我们重点介绍一下：这个脚本的运行原理。

1 数据

首先我们看下读取的data文件：

A) 共两列，两列之间用tab分开；

B) 第一列为缩放系数，第二列为对应结构的能量。

C) 在本文例子中，晶格参数等于= 缩放系数 * 2.8664 $\AA$ （记住这个公式）

2 脚本内容：浏览一遍即可，详解在下面：

详解：

2.1 1-15行：准备部分

A) 1-4 行为脚本的一些基本说明，为养成一个良好的写脚本习惯，以及避免后期不必要的麻烦，第一行必须要放进去。如果你用bash写的话，把python替换为bash即可，其他部分不用动，例如： #!/usr/bin/env bash。

B）2-3行注释，解释了一下脚本干什么用的，谁写的，怎么用。

C) 6-7 行调用python库中的math 和 numpy 模块，见下面两个链接：

https://docs.python.org/2/library/math.html

http://www.numpy.org/

D）9-15行：通过使用 print 命令输出一个警告，因为本文用的是缩放系数，很多人测试的时候保持缩放系数不变，直接采用不同的晶格参数。

2.2 读取data信息，准备拟合

A) 16-18行为脚本的注释说明，

B) 第20行中，我们使用 np.loadtxt读取data文件的内容，

注意：

1) data 用单引号括起来，表明我们要读以data命名的文件。

2) usecols=(0,1) 代表我们要读取第一列和第二列的内容。Python中的第一个用0表示，所以这里我们用的（0,1）来表示读取第1和第2列的内容；

3) delimiter='\t'，意思是第一列和第二列之间用tab（’\t’）分开；

4) unpack=True 意思是把读取的数据生成一个数列，并分配给a 和 E，具体见链接：

https://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html

通过第20行，我们可以获得两个数列，第一个是 a 包含了data文件中的缩放系数，第二个是E，包含了data文件中的能量信息

C) 第22-28行，把缩放系数转换为晶格常数(a2.8664)，然后平方取倒数（*(-2)），用来得到我们前面讲解的x。这是因为：当我们把BM方程写成： $y(x) = c_0 + c_1x + c_2x^2+c_3x^3 $ 的形式时，前提是令 (1/a)$^2$ = x。不懂的具体见Ex33中的内容。截止到现在，我们有了三个数列： a, E 和 x

2.3 第29-40行：拟合多项式方程：

A) 最小二乘法拟合 $y(x) = c_0 + c_1x + c_2x^2+c_3x^3 $这个方程，得到常数系列c，这里我们调用了numpy 的polyfit这个模块来进行拟合。具体见：

https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html

B) 第31行：本脚本的核心：p = np.polyfit(x, E, 3)

1) np.polyfit 读取x和 E的数值，拟合一个$y(x) = c_0 + c_1x + c_2x^2+c_3x^3 $，

注意：这里拟合的时候，用的数值是x和E，(1/a)$^2$ = x，所以，在进行拟合之前，务必确保你已经得到了正确的x数列。具体根据data文件的内容进行修改设置第24行的公式。

2) 括号中的3代表的是多项式方程中最大的那个次数。

扩展：如果你想拟合一个$ y = c_4x^4 +c_3x^4 …$的方程，把3改成4 即可。

3)拟合的结果输出到数列 p 中，数列p中共有4个常数（p[0], p[1]，p[2]，p[3]），按照次数从高到底的顺序排列，对应着我们Ex33中的（c3，c2, c1, 和 c0 ）

4) 33-36行, 我们将p 中的数值分别分配给c3，c2, c1, 和 c0。

2.4 38-40行，获取晶格参数：

到目前位置，多项式方程的常数我们已经获得了，下面就是求解问题了。

A) 第42行为下图中公式的代码描述，获取方程的解x，因为前面我们已经有了一个x变量了，这里方程的解用x1表示，

$x x = \frac{-2c_2 + \sqrt{(2c_2)^2 - 4\times3C_3\times C_1}}{2\times3C_3}$

B) 第43行获取晶格常数，公式为： (1/a)$^2$= x

C) 第45行，输出晶格常数的数值。

3 扩展练习：

3.1 根据解释，重新浏览脚本，知道每一行代表的含义。

3.2 本脚本使用的时候，需要自己根据data数列中的内容，设置更改一下获取 x 数列的公式（第24行）；

3.3 阅读Ex32-本节的所有内容，整个流程不要再存在什么疑问。

4 总结

本节我们详细讲解了BM方程拟合脚本，相信大家应该了解的差不多了，拟合的办法有很多，但BM方程只有一个。请务必掌握至少一个拟合的方法。后面的这句话本来想删掉，因为本人已经早就毕业了。暂且留下来吧，权当本书的一个时间证据。

此外，本人马上就要毕业答辩了，后面更新会有所放缓，努力争取保证每周一篇的速度。如果你有什么科研心得，欢迎分享。

Ex34 晶格参数的确定（Birch-Murnaghan状态方程-2）

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

前面一节，我们通过批量操作，获得了拟合BM方程的数据。这一节，我们通过这些数据，使用脚本进行拟合并获取Fe单胞的晶格参数。本节的脚本和读取的数据文件，可在QQ群文件和百度网盘中下载，见底部。

1 脚本使用方法：

在介绍脚本之前，先介绍一下脚本的使用方法：

首先：将前面一节得到的 data 文件和脚本放在一个目录下

其次：进入该目录，运行脚本的命令为（见下图）：

python bm0.py

注意： bm0.py 为脚本名，你也可以随意写成 birch, birch.py等。

通过脚本，我们可以得到Fe单胞的晶格参数为： 2.8332 Å

2 脚本的内容-1：

我们看一下bm0.py 这个脚本的内容：只有几行而已。大家先自己根据前面介绍的通过BM方程获得晶格参数的原理，阅读以下该脚本，看不懂不要紧，多看几遍。大体上了解是怎么回事。

3 脚本的内容-2：

前面的脚本没有任何的说明，阅读起来很是费力，大师兄把注释加上去了，并命名为: bm.py。大家再阅读一遍，看看有没有新的收获。

图中的链接为：(左下角查看原文)

https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html

该脚本运行效果如下图：

脚本的工作原理下一节介绍，本节中大家只看脚本，分析里面的内容，自己结合前面一节琢磨是怎么回事。脚本和读取的data文件可以在QQ群中下载，也可以通过百度网盘，链接： http://pan.baidu.com/s/1o8cCdBG

4 扩展练习

4.1 查找文献，获取Fe的晶格常数的实验数据，与脚本获取的结果进行对比；

4.2 Windows用户使用notepad++或者其他文本编辑软件，自己将bm0.py脚本重新写一遍，然后运行，如果出错了，请与图中对比并修改，直至完美匹配为止！！！

4.3 Linux用户，使用vim或者其他文本编辑软件，进行4.2中的操作；

4.4 结合Ex33的介绍，初步了解本脚本的运行原理。

Ex33 晶格参数的确定（Birch-Murnaghan状态方程-1）

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

O$_2$分子的计算中，我们知道了，单纯从数据库中获取的结构，只能作为一个合理的初始值，与计算所得到的理论结构还有一定的差距，因此我们需要对该结构进行优化才可以获取稳定的晶格参数信息。有两个方法可以实现：

1 Birch-Murnaghan状态方程拟合,
2 VASP计算中通过调节ISIF参数直接优化Bulk。

下面两节我们先讨论一下第一个方法：BM方程拟合。

1 什么是Birch-Murnaghan方程？

BM方程为块体材料的体积随着外压变化的状态方程。大师兄参考了维基百科的说明，大家下载维基百科导出的PDF文件参考。在学习下面的内容前，请务必阅读这个文件，了解什么是Birch-Murnaghan方程。下载：http://pan.baidu.com/s/1bUSYbs

2 如何通过BM方程获取晶格参数？

首先，我们将BM方程中的体积转化为晶格参数（a）的函数：

其中E(a)和E$_0$是晶格参数为a和a$_0$时bulk的能量，B$_0$，B$_0$’和V$_0$这些见BM方程的维基百科介绍。令$(1/a)^2 = x$，上图中的BM方程可以写成：

$y(x) = c_0 + c_1x + c_2x^2+c_3x^3 $ 的形式，

E$_0$, B$_0$,B$_0$’和V$_0$等写进c$_0$,c$_1$,c$_2$,和c$_3$这些常数里面，它们的具体形式此时并不重要，我们需要做的就是拟合这样的曲线去寻找使y(x)为最小值时的x取值。对y(x)这个方程求导数，dy(x)/dx = 0 的时候，便可以获取能量最低时的x值了，再由$(1/a)^2 = x$ 获得晶格参数。

dy(x)/dx = c$_1$+ 2c$_2$x + 3c$_3$x$^2$

令：dy(x)/dx = 0，即 c$_1$+ 2c$_2$x + 3c$_3$x$^2$ = 0

这是一个典型的二项式方程，求解很容易。

$x = \frac{-2c_2\pm\sqrt{(2c_2)^2 - 4\times3C_3\times C_1}}{2\times3C_3}$

x有两个值，其中负值被舍弃掉，最终：

$x = \frac{-2c_2 + \sqrt{(2c_2)^2 - 4\times3C_3\times C_1}}{2\times3C_3}$

3 获取拟合的数据：

经过前面的分析，我们需要的就是获得的c$_1$，c$_2$ 和c$_3$值。可以设置一系列不同的晶格常数进行计算，得到对应的能量，然后将这些数据带入下面的方程中：$y(x) = c_0 + c_1x + c_2x^2+c_3x^3 $进行拟合。本节我们主要讨论一下如何进行该计算。看到这里，你第一反应是我们之前测试K点或者SIGMA值时的那些Linux批量操作。是的，我们也需要进行一个批量操作的计算。

3.1 分数坐标

如果我们改变了晶格常数，那么Fe原子的坐标也要发生相应的变化。此时，用分数坐标更加方便，VESTA中在cif转VASP格式的时候，会提醒你选择哪个形式的坐标，我们选择分数坐标：

POSCAR如下图：

注意：

如果前面练习中，导出的为笛卡尔坐标形式的，你有两个选择：

A）重新按照前面的步骤，用VESTA导出为分数坐标形式；

B）将单点计算的结果直接复制过来（cp CONTCAR POSCAR），因为VASP的输出坐标为分数坐标形式，单点计算中，晶胞原子都没有动，CONTCAR可以直接使用。复制过来的POSCAR如下：

3.2 获得分数坐标的POSCAR后，我们可以批量对晶格常数a,b,c进行批量操作，也可以对缩放系数进行批量操作，而原子的坐标不用去管。

3.3 本练习对缩放系数进行批量操作

我们取10个点，将1.0 替换成0.95, 0.96, 0.97, 0.98, 0.99, 1.01, 1.02, 1.03, 1.04, 1.05。此时，晶格常数为：缩放系数*2.8664.（记住，下一节用得到）

a) 准备缩放系数为1.0 的计算文件夹，命名为1.00，里面包含我们前面一节已经讲过的Bulk单点计算相关文件；

1	for i in 0.95 0.96 0.97 0.98 0.99 1.01 1.02 1.03 1.04 1.05; do cp 1.00 $i ; sed -i "2s/1.0/$i/g" $i/POSCAR ; done

注意：这里大师兄在自己电脑的~/.bashrc文件里面设置了：alias cp='cp -r'

思考： sed 命令后面用的是双引号（为什么？）如果这里看不懂的话，建议从Ex00开始学习。

d.) 批量提交命令：

1	for i in *; do cd $i; qsub sub4; cd $OLDPWD; done

其中，qsub sub4 是大师兄提交单个任务的命令。

3.4 数据提取：

1	for i in *; do echo -e $i "\t" $(grep ' without' $i/OUTCAR \| tail -n 1\| awk '{print $7}'); done > data

3.5 以上计算文件以及相关的命令已经上传到QQ群文件和百度网盘了，大家对比下计算结果。 http://pan.baidu.com/s/1dEJr9rb

4 扩展练习

4.1 完成本练习中的相关操作，并学习VASP官网批量操作的例子：

https://cms.mpi.univie.ac.at/wiki/index.php/Fcc_Si

注意：FCC Si的例子用的原胞，并且该网址最后图中曲线不是二次方程！！！虽然用二次方程拟合出来的结果差不多。

4.2 探索学习曲线拟合的方法；

4.3 学习并了解BM方程以及如何拟合，如何通过BM方程获取稳定的晶格参数；

4.4 对晶格参数进行批量操作，并计算。

5 总结

本节从BM方程开始，简单介绍了什么是BM方程，如何通过BM方程获取稳定晶胞对应的晶格参数，最后，又复习了一下本书前面几节的批量操作，下一节我们讲解一下如何通过python脚本对BM方程拟合。

Ex32 Bulk的单点计算（Fe 单胞）

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

我们学习完了气相分子的相关计算后，下一步就是块体计算了，这也是VASP的强项所在，但块体的计算与气体分子的非常类似，如果前面掌握好了，后面的计算对你来说也就是轻而易举的事情了。

本着从简入繁的原则，我们先学习Bulk的单点计算（以Fe的单胞为例），姑且称之为本书中级篇的开始。大家需要学习并初步了解Material Studio和VESTA这两款软件。

Material Studio下载链接：http://pan.baidu.com/s/1i5or3ZR 对于这款软件怎么安装，大师兄就不指导了，网上全是相关的资料。本节默认大家已经安装好MS软件，并且可以打开界面。
VESTA下载官网：http://jp-minerals.org/vesta/en/download.html 只在官网下载，网上乱七八糟的版本不要去管！VESTA 下载后解压，直接打开就可以用了。注意：很多人问VESTA的使用说明，官网有，网上也有很多相关资料！大家耐着性子多练一天就能摸索个差不多了。

1 模型

1.1 课题的第一个难点

模型的选择，计算的准确性，以及结果分析的合理性，是一个课题是否成功的三个最主要的因素。Bulk的计算很简单，难点在于Bulk模型的获取。这也是一个课题最难的部分之一。很多人在计算的时候，模型不对，基本上这个课题就被一棍子打死了。这里大师兄提醒大家的是：模型，模型，模型！！！在进入计算的时候，一定要确保模型的合理性与正确性。这是计算中的第一个坎！

1.2 MS搭建Fe单胞的模型

A) 左上角file—> newproject，输入Fe，点OK;

B) file—> import—> Structures—>metals—>pure-metals—>找到Fe，选择打开即可;

C) file—>export—>选择cif格式（此时导出的是Fe的单胞，conventional cell）

1.3 MS导出Fe原胞的模型：

A)和B)步骤与前面一样，

C) Build —> Symmetry—> 选择primitive cell

D) file—>export—>选择cif格式（此时导出的是Fe的原胞，primitive cell）

1.4 VESTA 转换成VASP的格式

a) 打开VESTA软件，file—>open—>选择之前保存的cif文件

b) file—> export data —>保存类型选择VASP

c) 将保存的文件重新命名成POSCAR即可

1.5 需要掌握的知识：

a) 晶体学相关的基本知识，什么是primitive cell，什么是conventional cell?

b) MS和VESTA的基本操作

1.6 需要注意的部分（其他转换方法）：

a）可以使用openbabel转换：http://openbabel.org/wiki/Main_Page

b）可以使用MS生成的.cell文件转化，可以根据POSCAR的格式手动复制，也可以通过脚本。

c）或者通过其他脚本进行转化,例如：VTST中的cif2pos.pl。http://theory.cm.utexas.edu/vasp/scripts.html

d) 使用ASE: Atomic Simulation Environment的缩写。

2 VASP计算文件的准备

2.1 准备好INCAR文件：

A）Fe带有磁性，ISPIN和MAGMOM需要设置

B）Fe是金属， ISMEAR=1， SIGMA=0.1

C）ENCUT=450，统一起来，后面可能还需要计算其他元素，450是一个很安全的选择。

D）EDIFF控制电子步收敛的精度

2.2 准备KPOINTS文件：

以单胞为例： a=b=c=2.8664 $\AA$，根据前面的经验，我们可以设置$11\times11\times11$的KPOINTS（Gamma）

K-POINTS
 0
Gamma
  11 11 11
  0 0 0

2.3 准备好与POSCAR对应的POTCAR文件，提交任务的脚本

3 提交任务

3.1 提交任务之前，再次检查一遍我们的所有输入文件，确保无误；

3.2 提交任务进行计算；

3.3 大师兄的计算已经压缩放到百度网盘了，大家计算完毕后，与大师兄的进行对比。

链接：http://pan.baidu.com/s/1eRQHaX8

4 扩展练习

4.1 分析Fe体相的磁矩，并查找实验值；

4.2 复习之前O$_2$分子的相关计算，分析实验结果；

4.3 怎么在OUTCAR中查看磁矩等相关信息；（VASP 哪一个参数？）

4.4 复习晶体学的相关知识；

4.5 学会MS，VESTA的基本练习；

4.6 掌握2种以上通过MS获取POSCAR的方法。

5 总结

本节的计算我们需要掌握一定的晶体学基础知识，基本的MS和VESTA操作，计算过程和O$_2$分子的例子极为相似，大家在学习的时候，可以将这两部分相互结合。本节我们采用的Fe单胞的例子，计算本身不难。难点在于自己课题中模型的选取。

VASP计算遇到错误怎么办？

发表于 2017-12-22 更新于 2021-12-27 分类于 LVASPTHW附录阅读次数： Disqus：

搞科研，总会遇到各种各样的问题需要解决，同样VASP计算也不例外，本文大师兄列举了一些解决问题的步骤和方法，只提供大体的方向，没有具体事例，大家结合自己的计算研究亲自去实践体会，写本文的目的是告诉大家：自己亲自解决问题才是你快速成长的加速器。

第1阶段：

靠自己（最重要的！）

1.1 凭自己的经验解决；

1.2 按照VASP的输出错误，进行修改解决；

1.3 查找VASP官网的相关参数；

1.4 平时多阅读VASP手册；

1.5 不懂的参数不要往INCAR里面放。

第2阶段：

主动搜索（也是靠自己）

2.1 将错误信息复制到百度里面搜索，一般会得到一堆论坛或者博客的相关指导说明，按照里面的内容修改自己的输入文件，然后解决问题。

2.2 将错误信息复制到Google里面搜索，一般会得到VASP官网的论坛一节其他国外网址的相关错误示例，这时候，就要耐着性子把其中的内容读完，然后修改输入文件进行解决；一般来说，进行到这一步，至少60% 左右的问题可以解决!!!

第3阶段：

靠师兄师姐，师弟师妹，老板

如果时间过去2天，还没有解决。可以将自己的错误信息以及尝试的方法总结。

分享给组里的师兄弟，师姐妹们，寻求求助办法。如果还未解决，就要报告给自己的老师，求助指导！！！此时，解决的成功率应该能达到70-100%

大师兄的建议：自己主动去解决是关键，这需要培养主动解决问题的能力。老板是辅助的策略，除了指导课题的方向，在学生遇到问题时能出马解决，这也是他或者她为什么是老板的原因之一。

第4阶段（I）

论坛求助

如果老板不给力（很多做实验的老板要求学生去做计算，大师兄本人极度鄙视这样的老板，完全就是对学生的不负责任！），可以去尝试论坛求助：这里大师兄推荐：

将问题复制到百度里面，你就会发现这个问题会在很多论坛里面出现。不过论坛求助有一点需要注意：发帖的人水平参差不齐，要抱着怀疑的态度去接受所有人的答复（记住，是所有人！！！），别人的回复只能当作参考。

第4阶段（II）QQ群求助：

这里大师兄列举了三个管理员认真负责的QQ群：

遇见大师兄

第一性原理之家

物理化学材料互助群

进行到这一步，基本上95%的问题可以解决了。

第5阶段：VASP官方论坛发帖

以上都解决不了的话（极少的情况了），可以去VASP官网发帖求助，等待VASP的其他用户或者管理员解决。不过需要注册！

第 6 阶段

QQ群求助：

6.1 尊重别人的时间

A）QQ群求助的时候，提供的信息要全面：

1）算的什么体系

2）算的什么内容

3）INCAR长什么样子

4）错误信息

5）KPOINTS以及POSCAR 很多人把图片一截，或者拿手机拍张照片就直接在群里求助了，这是对大家时间的极度不尊重。而且被应助的概率很低！！！你不尊重大家的时间，大家也不会浪费时间在你身上。

B）可以将计算结果压缩，上传到QQ群文件进行求助（但不要把CHG，CHGCAR和WAVECAR等压缩到里面，因为它们太大，大家都懒得去下载！如果你上传了一个100M的求助文件，基本没人帮你！）

C）不要多群求助有很多人加了很多群，把遇到的问题在所有的群里面都发了一遍！这样的行为简直就是在浪费大家的宝贵时间，本人见到这种情况，一般都不会去帮助。

6.2 把你得小心眼收起来！！！

A) 除了涉密的一些计算，其他情况没有必要给大家隐瞒什么，大家也不会坏到去抢你的课题去算！心胸狭隘的人也不配得到帮助！

B) 要虚心接受别人的批评，也就是心眼要大。大师兄也遇到很多奇葩，在群里求助，没人回答，就直接退群了。被别人批评了几句，也退群了。只能说这些没心胸的人也不值得去帮助。

6.3 互相帮助

别人遇到的问题，如果你知道怎么解决，不要有所隐瞒。请认真帮助解决，用俺奶奶经常挂在嘴边的那句话：将心比心。这样大家的计算水平才可以得以普遍地提高。请认真按照这个流程来解决问题，不要一出现问题，脑子不假思索就去论坛或者QQ群求助!!! 请记住，你的问题也代表了你的老板的水平。一般来说，不负责任的老板带出来的学生总会有各种各样的低级问题。以上仅代表本人观点，不喜勿喷！祝大家踏踏实实做计算，多发实实在在的文章！

最后，新手的话，请认真学习Learn Vasp The Hard Way 这本书(www.bigbrosci.com) 按照大师兄写的内容，从头认真学习VASP（本书只作为参考，关键是去浏览官网！！！）；如果在学习中遇到问题，也可以把问题准备好，发送邮件给大师兄解决：lqcata@gmaill.com 最好不要直接QQ找我（拜托！）

自动生成KPOINTS的脚本

发表于 2017-12-21 更新于 2021-12-27 分类于 LVASPTHW附录阅读次数： Disqus：

前面学习了cd的一些用法，本节我们简单介绍下如何写自动生成KPOINTS文件的脚本。提到脚本，对于做计算的我们并不陌生，提交命令，分析数据，处理结果都会用到，脚本的存在使得一些繁重的工作极大地得到简化，节约了我们的时间和精力。但脚本怎么写出来的？怎么写脚本？很多人就望而却步了。

1 什么是脚本？

对于脚本的解释，大家可以浏览一下网上的解释，这里大师兄主要说一下自己的理解：

1）脚本是一个文本文件；我们可以用文本编辑器打开，查看，和修改。知道了这一点，当你从某处获得一个脚本时，就可以阅读里面的内容了，进而知道该脚本是怎么运行的，有哪些地方值得注意，为了实现另一个结果我们需要修改什么参数等；

2）脚本中的内容为程序语言；脚本中我们需要把我们期望的工作分解并转化为程序语言；所以你的选择有很多，bash， perl， python， java….， C++ 等等；

4）脚本可以执行；通过执行脚本，实现我们的目的，这里主要是给大家强调一下，在今后的学习中，仔细观察脚本的执行流程！

2 生成KPOINTS的脚本

这里说的KPOINTS文件指的是自动生成网格的KPOINTS文件，能带计算由于比较特殊，最好手动输入。编辑器用的是vim，还可以使用 notepad++或者其他文本编辑器进行练习。

2.1 kpoints.sh版本1

A）终端里面输入： vi kpoints.sh

B）在vi 界面里面输入下面几行内容：

输入完毕后保存。

C）脚本讲解：图中的除了第4行，所有命令都可以复制到 Terminal里面自己运行，然后查看结果。

C.1）前面我们学习了如何手动制作KPOINTS文件，并知道KPOINTS文件每一行所代表的内容。因此该脚本就是帮助我们自动生成每一行的内容。忘记的同学请复习Ex01的内容。

C.2）第一行： echo K-POINTS > KPOINTS （注意空格）

echo 在bash语言中，代表输出打印的意思，> 代表把前面的内容保存到后面的KPOINTS文件中；如果目录下不存在KPOINTS文件，那么会自动创建；

如果存在，KPOINTS文件之前的内容会被当前命令中 > 前面的内容替换。也就是KPOINTS文件以前的内容被清理掉，并换成了最新的内容。

因此，运行这一行命令，我们会创建一个KPOINTS文件，里面只有一行文字： K-POINTS

我们知道第一行的内容为说明，对计算不会产生影响，你也可以把 K-POINTS换成男朋友或者女朋友的名字….

C.3）第二行：

继续使用echo命令，将数字 0 保存到KPOINTS文件中。

这里 >> 两个箭头代表：将 >> 前面的内容保存到KPOINTS文件的最后一行。因此，使用 >> 不会将KPOINTS之前的内容替换掉。如果目录下没有KPOINTS文件，使用>>也会和>一样，创建一个。

这里的 0 代表自动生成网格：automatic generation scheme

C.4) 第三行：将 Gamma-centered 输出到KPOINTS的第三行里面，你需要知道G代表的是什么意思。如果你想用Monkhorst-Pack Grids，echo后面怎么写你要知道。建议一直用Gamma-centered。

C.5) 第4行：将$1 $2 $3输出到KPOINTS的第四行里面。

在这里，$1 $2 $3称为：arguments

代表的是我们运行脚本命令的时候，加入的三个方向上KPOINTS的数目。通过$1 $2 $3传递给 echo 命令。

如果不理解的话，跳过，等运行命令的时候就明白了。

C.6）第5行：Kpoints网格在三个方向的移动。一般来说0 0 0 即可。

2.2 脚本运行：

A) 调用bash运行脚本：

图中我们运行脚本的时候命令为：

1	bash kpoints.sh 3 3 1

bash 意思是我们调用bash来执行该脚本；

kpoints.sh 为脚本的名称, 当然啦，名字可以任意换， .sh也可以不用加；

3 3 1 为脚本后面的参数，第一个3 对应脚本里面的$1, 第二个3 为$2, 1 为脚本里面的$3

B）赋予脚本可执行权限：

赋予权限：chmod u+x kpoint.sh

取消权限：chmod u-x kpoint.sh

你会发现之前白色的脚本变成绿色的了，这说明本脚本可以执行了。不同电脑可能显示的颜色不容，有些甚至都不显示，大家不要纠结。但图中直接输入kpoints.sh 3 3 1 的时候失败了，我们必须用./ kpoints.sh 3 3 1 来执行。原因是系统不知道我们把脚本放在这里了。

怎么样才可以直接使用这个命令呢？将脚本转移到~/bin 文件夹中即可，如果你的系统中没有~/bin 这个文件夹，mkdir ~/bin 创建一下就完事。

OK，现在大功告成。我们可以根据之前的经验规则设置KPOINTS文件了。

2.3 脚本升级：版本2

前面的脚本实现了我们所期望的功能，这里我们所说的升级，无非就是在前面的基础上再增添点东西，使得脚本看起来更加高大上，具有可读性或者更加智能化。这里我们主要介绍一下使用# 注释自己的脚本，使其可以被广大吃瓜群众所理解和接受。

a）加入第一行，#!称为shebang （就是拼音，没错！）表示执行该脚本时会调用后面的程序；自行百度shebang的用途；

b） 3-5行中用 # 注释下，简单介绍本脚本的功能，作者，使用方法；

c） 7-11 行用 # 注释下每一行命令中的作用。

d）由于每个人写脚本的方式和风格不同，脚本的注释对于方便大家理解非常有用。有时候即使自己写的脚本，一年半载之后也会忘记里面的内容，加上注释，会让我们快速掌握脚本里面的内容，在后期的维护或者升级中也发挥着重大的作用。

2.4 其他升级：

除了注释外，可以将脚本写的更加智能化，通过调控 $1 $2 $3 等脚本后面的参数, 加入一些if语句等，来实现KPOINTS文件中所有行均可以通过命令进行修改。比如：kpoints.sh G 3 3 1 和 kpoints.sh M 3 3 1 分别生成Gamma centered 和正常的MP 网格。

Linux下面还可以使用其他的方式来生成脚本，总之选择的方法很多。

3 扩展练习：

3.1）重复本节所有操作，写人生中的第一个脚本；

3.2）网上搜索bash资料，认真练习。

3.3）bash在linux下处理日常工作非常有用，但数据分析的时候能力就有所欠缺了，建议大家学习python。Python有2.X和3.X的版本，2.X开始逐渐被抛弃，建议大家直接学习3.X的。

4 总结：

看完本节，你会初步了解脚本是怎么写出来的，它是怎么执行的以及对战胜了对写脚本的恐慌心理。当你获得一个脚本的时候，尝试着打开它，将其中的语言分解成若干命令执行，查看该脚本的工作流程。勤学苦练，就可以流利地写属于自己的脚本了。在这方面多花一分钟，以后节省的就不是一分钟的事情了，当然也有可能会多花几分钟。但多学习点东西总归是好的。

自动生成POTCAR的脚本

发表于 2017-12-21 更新于 2021-12-27 分类于 LVASPTHW附录阅读次数： Disqus：

前面我们介绍了KPOINTS的小脚本的写法，本节我们主要讲解一下生成POTCAR的脚本。学习完本节，脚本的神秘色彩基本烟消云散，大家都可以尝试写自己的脚本了。再强调说明三点：

1）Windows用户不要用记事本写脚本;

2）Windows用户写完脚本后，在服务器上先运行：dos2unix XXX (XXX 是你的脚本名)

3）Windows用户，如果计算出错，第一个要排除的是自己提交任务前有没有运行：dos2unix.

1 手动示例：

比如我们POSCAR中包含有Cu C H O 四种种元素，那么我们就需要按照顺序生成一个对应的POTCAR。首先回顾一下前面介绍的POTCAR的制备方法：

1.1）获取这四个元素各自的POTCAR：并命名成POTCAR_Cu,POTCAR_C,POTCAR_H,POTCAR_O

1.2）使用cat 命令将这四个POTCAR合并在一起：cat POTCAR_Cu POTCAR_C POTCAR_H POTCAR_O > POTCAR

前面我们学过了 > 的使用。分步练习：第一步将这些元素的POTCAR全部打印出来：

cat POTCAR_Cu POTCAR_C POTCAR_H POTCAR_O

第二步：将打印出来的内容通过 > 保存到最终的POTCAR里面。

脚本目的：

自动从赝势库中提取所需元素的POTCAR并将它们合并在一起。我们赝势库的目录为：/THFS/home/iciq-lq/bin/pot 这个pot文件夹（我自己随便命名的）就包含了VASP的所有元素的POTCAR，最好放在一个不经常变动的地方，这里我们放到了~/bin目录里面。

iciq-lq@ln3:/THFS/home/iciq-lq/bin/pot$ ls
Ac At_sv_GW  Bi_sv_GW  C_GW  C_s ....
iciq-lq@ln3:/THFS/home/iciq-lq/bin/pot$ ls B*
B:
POTCAR  PSCTR

Ba_sv:
POTCAR  PSCTR

Ba_sv_GW:
POTCAR  PSCTR

Be:
POTCAR  PSCTR

Be_GW:
POTCAR  PSCTR

所以我们可以使用下面的命令生成POTCAR：(不将这几个POTCAR先保存到当前目录下，直接调用赝势库中的POTCAR)

iciq-lq@ln3:/THFS/home/iciq-lq$ cat /THFS/home/iciq-lq/bin/pot/Cu/POTCAR  /THFS/home/iciq-lq/bin/pot/C/POTCAR /THFS/home/iciq-lq/bin/pot/H/POTCAR  /THFS/home/iciq-lq/bin/pot/O/POTCAR  > POTCAR
iciq-lq@ln3:/THFS/home/iciq-lq$ grep TIT POTCAR
   TITEL  = PAW_PBE Cu 22Jun2005
   TITEL  = PAW_PBE C 08Apr2002
   TITEL  = PAW_PBE H 15Jun2001
   TITEL  = PAW_PBE O 08Apr2002
iciq-lq@ln3:/THFS/home/iciq-lq$

但这样子也很麻烦。

2 脚本预览：

首先看下自动生成POTCAR的脚本，大家尽量先不看下面的解释，自己理解一下脚本的整个框架和运行的流程;

#!/usr/bin/env bash
# Create a GGA_PAW POTCAR file
# by BigBro  
# To Use it： potcar.sh Cu C H O

# Define local potpaw_GGA pseudopotentialrepository:
repo="/THFS/home/iciq-lq/bin/pot"

# Check if older version of POTCAR ispresent
if [ -f POTCAR ] ; then
 mv -f POTCAR old-POTCAR
 echo " ** Warning: old POTCAR file found and renamed to 'old-POTCAR'."
fi

# Main loop - concatenate the appropriatePOTCARs (or archives)
for i in $*
do
 if test -f $repo/$i/POTCAR ; then
  cat $repo/$i/POTCAR>> POTCAR
 else
 echo " ** Warning: No suitable POTCAR for element '$i' found!! Skipped thiselement."
 fi
done

讲解，这个脚本里面：

1）首先我们先设定赝势库的目录;repo="/THFS/home/iciq-lq/bin/pot" 等号前后不能有空格！！！

2）然后我们检查一下当前目录下是否存在POTCAR，如果存在的话，将其重命名为：POTCAR-old, 并输出警告信息。这里我们学到的是 if 句型：if [ ] ; then; fi 。[] 中间是if判断的条件，如果成立，则继续执行then后面的动作。

if [ -f POTCAR ]; then 意思是如果存在POTCAR这个文件，那么… 注意：

a）[ ] 和里面的内容 –f POTCAR 要通过空格分开。所以：

if [-f POTCAR ]（[ 和-f之间没有空格）
if [ -f POTCAR] (POTCAR 和] 之间没有空格)
if [-f POTCAR]（-f，POTCAR和 [] 之间都没有空格）

这三种写法都是错误的;if [ -f POTCAR ] 是正确的

b）then 后面如果另起一行的话，分号（;）可加可不加;

c）then 后面的内容如果和then 在同一行，必须加上分号（;），前后有无空格均可，建议加上，这样语言会很清晰;

d）then 后面的命令执行完毕后，要加上 fi 结束。同理，如果fi 另起一行，则前面可以不加;如果fi 和 then 后面执行的命令在一行，前面需要加上分号（;）。

3）生成新的POTCAR文件; 这里我们用到的是一个for循环：

a) for i in $*

这一行的意思是，对于命令后面的所有参数（arguments），前面一节我们已经学过了$1 $2 ... 的含义，这里用的$* 可以让我们在命令后面加任意数目的参数;因为另起了一行，后面加不加分号（;）不重要。

b） do 执行的意思

c）do后面有个if 语句，如果目录下面是POTCAR的话，那么使用cat 命令将所有$*对应的POTCAR输出到最新的POTCAR中，注意，这里我们用到的是>>; 因为for循环是对于后面的参数挨个执行的。

d) 有两个elif，判断目录下的另外两个POTCAR文件的格式，如果是POTCAR.Z 文件的话，则使用zcat 将其输出到最新的POTCAR;如果是POTCAR.gz 文件的话，则使用gunzip –c 命令。

e）else 指的是除了前面三种以外的情况，如果POTCAR，POTCAR.Z或者POTCAR.gz都不存在，那么使用echo命令，输出警告，提醒用户检查。

f）使用fi 结束if 语句

g）使用 done 结束 for 循环。

3 运行脚本：

3.1）赋予脚本可执行权限，然后将其转移到主目录下的bin文件夹中，前面几节已经讲过了。

1 2	chmod u+x potcar.sh mv potcar.sh ~/bin

3.2 运行实例：

iciq-lq@ln3:/THFS/home/iciq-lq$ ls
bin  LVASPTHW  old-POTCAR  POTCAR  potcar.sh  test_jobs
iciq-lq@ln3:/THFS/home/iciq-lq$ potcar.sh Cu C H O
 ** Warning: old POTCAR file found and renamed to 'old-POTCAR'.
iciq-lq@ln3:/THFS/home/iciq-lq$ ls
bin  LVASPTHW  old-POTCAR  POTCAR  potcar.sh  test_jobs
iciq-lq@ln3:/THFS/home/iciq-lq$ grep TIT POTCAR
   TITEL  = PAW_PBE Cu 22Jun2005
   TITEL  = PAW_PBE C 08Apr2002
   TITEL  = PAW_PBE H 15Jun2001
   TITEL  = PAW_PBE O 08Apr2002
iciq-lq@ln3:/THFS/home/iciq-lq$
iciq-lq@ln3:/THFS/home/iciq-lq$ potcar.sh BigBro
 ** Warning: old POTCAR file found and renamed to 'old-POTCAR'.
 ** Warning: No suitable POTCAR for element 'BigBro' found!! Skipped thiselement.
iciq-lq@ln3:/THFS/home/iciq-lq$

当前目录下存在POTCAR，脚本会给出警告，将其重命名为old-POTCAR，然后生成新的POTCAR文件，包含Cu C H O 这四种元素。如果有一个元素不存在（上面实例总的BigBro），potcar库里找不到，则输出错误提示。

注意！注意！注意！

上面例子中的错误，还可能是因为：

1）你的POSCAR的格式有问题（Windows用户常见的问题），如果发现这个错误，先运行下：dos2unix POSCAR 这个命令排除是Windows的格式问题。
2）如果还不能解决，去你的POTCAR库里面找找，是不是真的没有这个元素的POTCAR文件。

4 脚本分析：

4.1）生成一些特殊的POTCAR文件：

我们知道，每一个元素存在好几个不同的POTCAR，如图：

iciq-lq@ln3:/THFS/home/iciq-lq/bin/pot$ ls Cu*
Cu:
POTCAR  PSCTR

Cu_GW:
POTCAR  PSCTR

Cu_pv:
POTCAR  PSCTR

Cu_sv_GW:
POTCAR  PSCTR
iciq-lq@ln3:/THFS/home/iciq-lq/bin/pot$ ls H
H/        H1.33/    H1.66/    H.25/     H.42/     H.58/     H.75/     He/       He_GW/    Hf_pv/    Hf_sv_GW/ Hg_sv_GW/ H_h/      Ho/       H_s/
H1.25/    H1.5/     H1.75/    H.33/     H.5/      H.66/     H_AE/     He_AE/    Hf/       Hf_sv/    Hg/       H_GW/     H_h_GW/   Ho_3/

当我们需要使用 Cu_pv, C，H1.25，和O对应的POTCAR时，我们需要这样执行命令：

iciq-lq@ln3:/THFS/home/iciq-lq$ ls
bin  LVASPTHW  test_jobs
iciq-lq@ln3:/THFS/home/iciq-lq$ potcar.sh Cu_pv C H1.25 O
iciq-lq@ln3:/THFS/home/iciq-lq$ ls
bin  LVASPTHW  POTCAR  test_jobs
iciq-lq@ln3:/THFS/home/iciq-lq$ grep TIT POTCAR
   TITEL  = PAW_PBE Cu_pv 06Sep2000
   TITEL  = PAW_PBE C 08Apr2002
   TITEL  = PAW_PBE H1.25 07Sep2000
   TITEL  = PAW_PBE O 08Apr2002
iciq-lq@ln3:/THFS/home/iciq-lq$

因此，如果我们使用脚本的时候直接输入元素符号，默认读取的是与该元素符号相对应文件下的POTCAR。

5 读取POSCAR，自动生成POTCAR （防止出错的一个秘诀）

在我们的计算中，可以通过读取POSCAR中的元素信息，生成对应的POTCAR文件，这样可以避免POSCAR中元素和POTCAR中不一致的情况。那么我们需要怎么做呢？

5.1 首先我们要获取POSCAR中的元素信息：查看一下POSCAR文件,发现元素信息在第6行，因此加了一个tail命令提取出来：

iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/potcar$ ls
POSCAR
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/potcar$ head -n 6 POSCAR
Ru   H
   1.00000000000000
   8.1377999784000004   0.0000000000000000   0.0000000000000000
   4.0688999892000002   7.0475415119999996   0.0000000000000000
   0.0000000000000000   0.0000000000000000  21.5631999968999999
Ru   H
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/potcar$ head -n 6 POSCAR  | tail -n 1
Ru   H

5.2 使用我们刚刚写出来的脚本，通过使用$()调用前面命令的结果: （将获取的元素信息作为脚本命令的参数）

iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/potcar$  potcar.sh  $(head -n 6  POSCAR | tail -n 1)
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/potcar$ grep TIT POTCAR
   TITEL  = PAW_PBE Ru 04Feb2005
   TITEL  = PAW_PBE H 15Jun2001

大功告成！！！！

5.3 不想每次输入这么长的命令：可以将其写进 ~/.bashrc文件中：

1	alias pos2pot="potcar.sh $(head -n 6 POSCAR \| tail -n 1)"

iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/potcar$ ls
POSCAR  POTCAR
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/potcar$ vi ~/.bashrc
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/potcar$ . ~/.bashrc
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/potcar$ pos2pot
 ** Warning: old POTCAR file found and renamed to 'old-POTCAR'.

5.4 看完上面的部分，相信大家已经可以掌握很多内容了，额外福利，大师兄本人并没有使用head -n 6 POSCAR | tail -n 1 来获取POSCAR中的元素信息。而是使用的sed命令：

1 2	iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/potcar$ sed -n 6p POSCAR Ru H

6 扩展练习：

6.1）手动将这个脚本打一遍，理解里面的内容，如果出错了，对比我的脚本，进行改正，直至可以正确运行为止;

6.2）去官网查看POTCAR相关的内容，了解每一个元素不同的POTCAR之间的区别，例如：Cu 和 Cu_pv

6.3) 复习alias，学习本文第5节的方法，思考一下，如何通过操控脚本，避免计算中的一些失误，有了自己的想法后，通过脚本付诸实践。

7 总结：

写到这里，相信大家对于Linux系统下的一些基本操作已经熟练掌握了，脚本的怎么写，怎么去运行，心里也有了一个大致的了解，至少不会感觉写脚本是多么牛逼的工作了。这些内容都是我们保证计算准确，高效的很重要的一个因素。希望大家可以熟练掌握，最关键的是，要硬着头皮去练习，脚本看起来简单，等你写的时候就会出现各种各样的错误，只有不断练习，才能写出漂亮的脚本，提高自己的工作效率。脚本已经上传至QQ群文件，也可百度网盘下载：链接：https://pan.baidu.com/s/1eqeVVOo5nzGJsfPgClNJDg

Linux命令之cd使用心得

发表于 2017-12-20 更新于 2021-12-27 分类于 LVASPTHW附录阅读次数： Disqus：

科研狗必备技能

前面27节，我们已经通过气相分子的例子，带领大家顺利度过了菜鸟这一关。如果你掌握了其中的思想，勤加操练，看官网，认真阅读看大师兄推荐的参考书，相信你的能力足以甩出那些整天只知道搜集VASP学习资料，搜集DFT阅读书籍，而不去行动的人好几条街了。科研忌讳的就是眼高手低，你学到的知识和你下载的书成反比。

在菜鸟篇结束之后，大师兄需要把自己日常工作的一些小技巧展示给大家。一是作为大家摆脱菜鸟身份的奖励，二来为大家的进阶篇做好准备。今天我们主要讲解一下cd命令。

在Linux系统下，cd是必不可少的一个命令: cd 可以实现进入文件夹，返回上一级文件夹等操作等，日常操作中一些简单的做法可以极大的简化我们的工作方式，减轻手指的负担，有效避免键盘手的形成。大师兄结合自己多年来理论计算化学的操作实践，简单列举了一些组合使用 cd 命令的方法，希望大家可以熟练使用。

注意： 本节没有大师兄的截图操作，大家需要认真练习，思考这些命令是怎么运行的。希望大家阅读完本文后，按照内容练习一番，至少要达到高级篇的水平。此外，本节内容和前面的密切相关，如果前面没有掌握的话，在具体的应用部分中阅读起来会有些吃力。此时，建议从头学习。

1 菜鸟篇： 正常使用 cd 命令

进入文件夹： cd xxx/xxx/xxx

返回上一级： cd ..

返回上二级：cd ../..

2 初出茅庐： 懂得如何快速返回原来的文件夹目录

1)：cd -

2)：cd $OLDPWD

两者效果是一样的，区别是前者会在屏幕显示上一级的目录。

具体使用参见本书前面几节关于VASP批量操作的相关内容，尤其是在批量提交任务中的应用。这个命令一定一定一定要掌掌掌掌…握！！！

3 进阶篇： 知道如何结合 cd 与 alias

试想一下，如果你的文件系统下，一层套一层，好东西藏的很深（大师兄内心荡漾着邪恶清脆的笑声），问题是，当你进入这个文件夹，后退时需要不停地敲击：cd ../../../../../../../..

为避免这样的麻烦，可以这么做:

1
2
3

alias ..='cd ../../'
alias ...='cd ../../../'
alias ....='cd ../../../..'

每个点代表向上返1级目录，后面自己补上，大师兄手指头不想继续写了。

注意：

1）：不要设置 alias **.**='cd ../' 因为 **.** 在linux里面等于 source 命令。

2）：等号左右没有空格。

3）：等号右侧命令需要用单引号或者双引号括起来，最好是用单引号。

不会的请自己查找前面关于alias的介绍章节。

4 高级篇： cd 与 alias、& 、ls的结合

大师兄linux走火入魔，养成了一个进入文件夹就敲击 ls 命令的习惯，但是也舍不得折磨我那双精细的老手，于是，就组合了cd、alias、&、和ls 命令，如下alias ....='cd ../../../.. && ls'

这样你敲击4个点 .... 返回上面4级目录的时候，自动显示当前目录下的所有东西，避免再一次敲击ls 命令了。

在这里我们用到的是 && 这个命令，它表示如果前面的命令成功了就进行 && 的命令。

无聊的师弟师妹们可以这样玩一下：cd ../../../.. && cd -

但是如果 && 前面的命令失败了，直接退出，后面的便不会运行，这种情况下你就需要用到 双排管 这个命令了。

什么是双排管？　答：　||　（中间没有空格，大师兄给起的名字）

师兄，这玩意怎么使用？看下面的例子就明白了：

1	cd ../abc && ls

上一级目录没有abc文件夹，&& 前面命令的失败了，ls 不会运行。

cd ../abc || ls 上一级目录没有abc文件夹， || 前面命令的失败了，ls依然会运行。

此外，还有单管 | （pipe）这个命令，后面一般跟着前面一步结束的进一步操作。比如：

1） head –n 10 POSCAR | tail –n 1

2） grep ' without' OUTCAR | tail -n 1

3） grep ' without' OUTCAR | awk '{print $4}'

4）前面获取零点能的就是一个很好的例子

5） ……（大家自由发挥）

5 cd 私人定制篇

如果你的一个课题的计算都在一个文件夹里面，可以这么设置，下次敲击命令时直接进入这个文件夹

1	alias xxx='cd ~/bigbro/a/b/c/d/e/f/g/h/i/k'

藏得这么深（师弟们懂的！），肯定是好东西，不拿出来分享愧对大师兄谆谆教诲！学会了这个命令，就避免了用鼠标乱点进入目录了。还可以避免文件找不到….

6 cd 整蛊篇

如果你闲的蛋疼，还可以尝试着在办公室其他人的Linux终端下这样设置alias。二选一即可。

1 2	alias cd='cd ~' alias cd='echo "you are a dumbass" '

大家好好练习这些关于cd的基本命令，等手指头熟悉了cd的味道，大师兄下一节就会继续传授给大家其他的高级操作，进一步减轻手指头的负担。

Ex27 乙醇分子的振动频率计算(五)

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

通过前面读者的反馈，本节前面主要讲述一下零点能相关的注意事项以及频率计算的时间成本相关的知识。

1 参考书籍：

以下所提到的书籍都可以在QQ群文件和百度网盘里面下载：http://pan.baidu.com/s/1bKHMjG

1.1 Density Functional Theory：A Practical Introduction

大师兄极力推荐的菜鸟入门书！！！

1）不要网上随便下载（很多都是野鸡版，公式符号不全）

2）不要看中文版的；去看原汁原味的，困了就歇会，慢慢坚持看；

3）大师兄QQ群下载或者点击百度网盘链接；

4）第5章专门讨论的频率相关的问题！！！

1.2 Atkins的物理化学第9版

前面几节我们都是用的这一部分作为参考，忘了的童鞋们返回继续重修；这本书网上已经可以下载第十版的pdf文件了。

1.3 Jens. K. Norskov 的 Fundamental Concepts in Heterogeneous Catalysis

2 频率计算的细节问题

当你阅读完前面推荐书中的频率部分后，应该会明白很多。本节的第二部分我们主要讨论下频率计算的一些大家有顾虑的细节问题：

2.1 怎么进行零点能校正？

以乙醇的计算为例：

A）首先结构优化完毕后我们会得到分子的能量：$E_{0}$ （OSZICAR中的E0）

B）频率计算后我们会得到分子的零点能：$ZPE$

C）零点能校正之后的乙醇分子能量为： $E_{ZPE} = E_0 + ZPE$

D）A 和 B 得出的结果直接相加即可，不要想太多。

2.2 怎么计算过渡态和反应热的零点能校正？（菜鸟暂且跳过）

对一个反应：IS —> TS —> FS

IS: Initial State 反应物

TS: Transition State 过渡态

FS: Final State 产物

记住这几个定义，本书后面 IS, TS, FS 用的时候多着呢！

1）优化反应物 IS 和产物 FS 的结构，获得能量：$E(\textrm{IS})$, $E(\textrm{FS})$;

2）对反应物和产物进行频率计算，获得各自的零点能：$\textrm{ZPE(IS)}, \textrm{ZPE(FS)}$。

3）搜索过渡态，获得结构和能量 $E(\textrm{TS})$;

4）过渡态频率分析，获得零点能 $\textrm{ZPE(TS)}$

不考虑零点能的反应能垒 ($E_a$) 和反应热 ($\Delta E$)：

$E_a = E(\textrm{TS}) – E(\textrm{IS}) \\ \Delta E = E(\textrm{FS}) – E(\textrm{IS})$

考虑零点能校正：

$\begin{align} E{a(\textrm{ZPE})} &= E_{\textrm{ZPE}}(\textrm{TS}) – E_{\textrm{ZPE}}(\textrm{IS}) \\ &= E(\textrm{TS}) + \textrm{ZPE(TS}) – E(\textrm{IS}) – \textrm{ZPE(IS}) \\ &= E_a + \textrm{ZPE(TS)} – \textrm{ZPE(IS)} \end{align}$

同理:

$ \Delta E(\textrm{ZPE}) = \Delta E + \textrm{ZPE(FS)} – \textrm{ZPE(IS)} $

两个处理方式：

a. 先获取未校正的结果，然后把零点能各自相减；

b. 先将各个物种进行零点能校正，然后在计算反应能垒或者反应热

效果是一样的

2.2 频率计算的时候，是不是体系中所有的原子都放开？

不一定；这取决于你的体系，以及你想要关注的部分:

例子A:

乙醇在Cu(111)表面上的吸附，计算吸附热的零点能校正

$\ce{CH_3CH_2OH + Cu(111) -> Cu(111)-CH_3CH_2OH}$

此时，Cu(111)表面我们在计算频率的时候是要固定住的！只振动乙醇分子即可。

例子B:

计算CO在Cu(111)表面上的吸附：

同A，固定Cu(111) 表面，如果你只关心 CO 在垂直表面上的振动，那么CO的 $x, y$ 方向便可以固定住，在坐标后面为：F F T

例子C:

苯酚在Cu(111)表面上 O—H 键断裂活化能的零点能校正：

在这里我们拿苯酚作为例子，很多时候，计算的对象比较大，全部频率优化非常耗时，那我们就得选择性地固定住一部分，只关心关键的局域部分。这个例子中我们主要讨论零点能对 O—H 键断裂活化能的影响，因此我们可以把苯环的部分固定住，只放开 O 和 H 原子进行振动。

当然啦，这是一个简化的计算，肯定不如全部振动的结果好，但迫于计算时间的压力，这是一个折中的好办法。此外，由于在反应中苯环部分变化不大，所以即使我们全部放开，它对反应热和反应能垒的影响也很小，零点能部分相减的时候抵消了。在这里你要深刻体会到，意识到，把握到相对这两个字在理论计算中的巨大作用！！！

还有一点要注意的是：IS，TS和FS中，所固定和放开的原子必须一致！！！

不要在IS中固定苯环，而在TS或者FS中放开， 这样计算出来的结果只能用脑残来形容了。(ps. 本人犯的这种错误太多啦，已经不能用脑残来形容了！)

3 影响频率计算的因素测试:

前面提到过：ENCUT, KPOINTS, POTIM, EDIFF都会对频率计算产生影响。频率计算很多时候会有小的虚频出现，通俗点说就是数值噪音。我们可能会想到，凡是可以提高精度的办法是不是就可以降低噪音了呢？写到这里，大师兄突然想起初中的一个老师说过一句话：好记性不如烂笔头，意思是让我们多写，这样记得会比较快。还有一句名言：纸上学来终觉浅，绝知此事要躬行！学习 VASP，亲自上手实践及其重要，看遍官网，翻烂本书不去行动也是白搭，好听点说只能是纸上谈兵。

下面的测试主要从零点能和虚频，和时间三个方面来判断：测试的命令为：

获取虚频命令：

1	grep 'f/i' */OUTCAR \| awk '{print $1 "\t " $2 "\t" $8 "\t " $9 "\t" $10 "\t" $11}'

获取时间命令：

1	grep Elapsed */OUTCAR \| sort -n

获取零点能：

1	for i in * ; do echo $i $(cd $i ; fsum ; cd $OLDPWD);done \| sort –n

3.1 EDIFFG测试

3.1.1测试

$4, 5, 6, 7, 8$ 分别代表 EDIFF=1E-4 到 -8 的测试任务；

1	for i in * ; do echo $i $(cd $i ; fsum ; cd $OLDPWD); done

零点能的计算结果：

EDIFF	ZPE	Unit
4	2.14344	eV
5	2.12494	eV
6	2.11698	eV
7	2.11696	eV
8	2.11696	eV

a. 可以看出来在1E=-6以后，零点能变化基本为 0了；我们可以用这个设置；

b. 虽然你测出来一个很好的参数值，很不幸，这是默认的。

c. 从时间上考虑，本人会选择 EDIFF= 1E-5，这比 1E=-6 快了将近一倍；

d. $5$ 和 $6$ 之间的误差为：$ 0.008/2.117 = 0.38\% $；

e. 此外，单个的零点能貌似没有什么用处，一般都是两个结构的零点能相减。这个操作也会抵消一部分的计算误差。

3.1.2 虚频

（这个命令只展示一次，后面测试结果全部用下面的表格）

EDIFF	Nth	WN	Unit	Energy	Unit
4	25	10.183435	cm^-1	1.262585	meV
	26	43.891375	cm^-1	5.441839	meV
	27	61.646745	cm^-1	7.643225	meV
5	25	6.773094	cm^-1	0.839757	meV
	26	28.045092	cm^-1	3.477149	meV
	27	52.28955	cm^-1	6.48308	meV
6	25	6.132037	cm^-1	0.760276	meV
	26	18.62878	cm^-1	2.309675	meV
	27	65.657693	cm^-1	8.140519	meV
7	25	6.03848	cm^-1	0.748676	meV
	26	19.115009	cm^-1	2.36996	meV
	27	65.726299	cm^-1	8.149025	meV
8	25	6.036457	cm^-1	0.748425	meV
	26	19.156334	cm^-1	2.375084	meV
	27	65.711053	cm^-1	8.147135	meV

a. 改变收敛标准，并没有消除虚频，测试 $5$ 中，稍微降低了些，其他情况下，最大的仍然都再 1500 px-1 左右

b. $6$, $7$, $8$ 中的虚频基本一致，说明增加收敛标准，并没起到什么好的效果。

3.1.3 时间

EDIFF	Time	Unit
4	1702.463	s
5	2295.233	s
6	4298.571	s
7	5481.423	s
8	5998.518	s

复习下谁偷走了我的机时中EDIFF的因素；增加EDIFF，计算时间变长了。（离子步数一致，每个离子步的收敛步数增多了）一般来说EDIFF采用默认值或者1E-5都可以满足情况，任务比较繁重的时候议用1E-5。本节在后面的测试中，都采用 VASP 的默认EDIFF值进行（1E-6）

3.2 ENCUT测试

3.2.1 零点能

ENCU	ZPE	Unit
400	2.11698	eV
500	2.11448	eV
600	2.11880	eV
700	2.11971	eV
800	2.11928	eV

从表中可以得出结论：增加ENCUT,对零点能的影响很小；

在计算的过程中，ENCUT要保持一致，一但选定后就不要再轻易改动，这个表的数据告诉我们，即使改动了ENCUT，对零点能的影响也很小，所以ENCUT的影响可以忽略。

3.2.2 虚频：

ENCUT	Nth	WN	Unit	energy	unit
400	25	6.132038	cm^-1	0.760276	meV
	26	18.62878	cm^-1	2.309675	meV
	27	65.65769	cm^-1	8.140519	meV
500	24	0.489004	cm^-1	0.060629	meV
	25	16.7313	cm^-1	2.074418	meV
	26	23.491849	cm^-1	2.912619	meV
	27	66.062971	cm^-1	8.190767	meV
600	24	0.288462	cm^-1	0.035765	meV
	25	4.853849	cm^-1	0.601801	meV
	26	6.324527	cm^-1	0.784142	meV
	27	52.075614	cm^-1	6.456556	meV
700	25	0.747158	cm^-1	0.092636	meV
	26	6.175845	cm^-1	0.765708	meV
	27	32.540704	cm^-1	4.034535	meV
800	25	0.289332	cm^-1	0.035873	meV
	26	7.356367	cm^-1	0.912074	meV
	27	26.44812	cm^-1	3.27915	meV

增加ENCUT可以减小虚频；从 $400$ 增加到 $800$ 后，最大的虚频波数从 $65~cm^{-1}$ 减小到 $26~cm^{-1}$，但 A 中的结果是不是要求我们算频率的时候，必须要把 ENCUT增大呢？答：没必要！！

因为前面零点能变化甚微，且虚频从 $65~cm^{-1}$ 减小到 $26~cm^{-1}$, 只能说是在误差范围内变化，我们是可以承受的，此时ENCUT在减小虚频中的作用，大家记住这一点就可以了。此外，增加ENCUT也会使得计算时间变长（见下图）。但是如果虚频大于$100~cm^{-1}$，且不是反应路径的时候，你就得小心了。

3.2.3 时间

命令：

1	grep Elapsed */OUTCAR \| sort –n

ENCUT	Time	Unit
400	1832.405	s
500	2389.045	s
600	3019.433	s
700	4285.491	s
800	11943.02	s

3.3 PREC

这里我们设置了三个参数：Normal, Accurate和High

3.3.1 零点能变化甚微

PREC	ZPE	Unit
N	2.11698	eV
A	2.1147	eV
H	2.11363	eV

3.3.2 虚频

PREC	Nth	Wavenumber	Unit	Energy	Unit
N	25	6.132037	cm^-1	0.760276	meV
	26	18.62878	cm^-1	2.309675	meV
	27	65.657693	cm^-1	8.140519	meV
A	25	5.493108	cm^-1	0.681059	meV
	26	10.478293	cm^-1	1.299143	meV
	27	65.861524	cm^-1	8.165791	meV
H	24	0.552665	cm^-1	0.068522	meV
	25	16.954329	cm^-1	2.10207	meV
	26	29.95621	cm^-1	3.714098	meV
	27	68.287802	cm^-1	8.466611	meV

PREC同样对消除虚频不管用！！！！

3.3.3 时间

Prec	Time	Unit
N	4298.571	s
A	7596.994	s
H	7962.539	s

通过前面的虚频和零点能，结合此时的时间，我们可以得出结论：

PREC = Normal 完全够用了。

3.4 POTIM 参数

3.4.1 零点能：

POTIM	ZPE	eV
0.001	2.17125	eV
0.005	2.12616	eV
0.015	2.11769	eV
0.020	2.11698	eV
0.050	2.11671	eV
0.100	2.12481	eV

POTIM太小的时候，对零点能影响很大。（$0.001$ 和 $0.005$，$0.015$ 对比）

3.4.2虚频

POTIM	Nth	Wavenumber	Unit	Energy	Unit
0.001	25	21.746998	cm^-1	2.696285	meV
	26	39.656607	cm^-1	4.916794	meV
	27	56.54809	cm^-1	7.011072	meV
0.005	24	0.675194	cm^-1	0.083713	meV
	25	12.785399	cm^-1	1.585188	meV
	26	36.264179	cm^-1	4.496187	meV
	27	61.287932	cm^-1	7.598738	meV
0.015	25	5.626635	cm^-1	0.697614	meV
	26	23.117605	cm^-1	2.866219	meV
	27	59.01173	cm^-1	7.316525	meV
0.020	25	6.132037	cm^-1	0.760276	meV
	26	18.62878	cm^-1	2.309675	meV
	27	65.657693	cm^-1	8.140519	meV
0.050	24	0.327964	cm^-1	0.040662	meV
	25	23.575077	cm^-1	2.922938	meV
	26	24.337978	cm^-1	3.017526	meV
	27	95.239085	cm^-1	11.80815	meV
0.100	23	1.574379	cm^-1	0.195198	meV
	24	5.865037	cm^-1	0.727172	meV
	25	20.89784	cm^-1	2.591003	meV
	26	45.241296	cm^-1	5.609208	meV
	27	168.238904	cm^-1	20.85897	meV

POTIM太大的时候，会搞出来超大号的虚频（$0.050$ 和 $0.100$）！

3.4.3 时间

POTIM	Time	Unit
1	1530.591	s
5	2065.911	s
15	3845.977	s
20	4298.571	s
50	4829.897	s
100	5191.23	s

POTIM 越小，单个离子步收敛的越快，也就是需要更少的电子步！

综合前面的三个因素：POTIM太小或者太大都不好， POTIM = 0.015或者0.020 是很好的选择。 0.015是默认值。

3.5 POINTS

这里用 $1$, $2$, $3$ 来分别代表 $1 \times 1 \times 1$， $2 \times 2 \times 2$ 和 $3 \times 3 \times 3$ 的 K 点。

3.5.1 零点能：

KPOINTS	ZPE	Unit
1	2.11698	eV
2	2.11678	eV
3	2.11689	eV

3.5.2 虚频：

KPOINTS	Nth	Wavenumber	Unit	Energy	Unit
1	25	6.132037	cm^-1	0.76028	meV
	26	18.62878	cm^-1	2.30968	meV
	27	65.657693	cm^-1	8.14052	meV
2	25	5.726298	cm^-1	0.70997	meV
	26	19.193163	cm^-1	2.37965	meV
	27	64.645347	cm^-1	8.015	meV
3	25	4.826986	cm^-1	0.59847	meV
	26	19.138554	cm^-1	2.37288	meV
	27	64.83042	cm^-1	8.03795	meV

3.5.3 时间

KPOINTS	Time	Unit
1	4298.571	s
2	8997.709	s
3	13480.847	s

综合前面三点： Gamma点足矣！

4 扩展阅读：

4.1 下载本节的练习，按照本节的顺序操作分析；

4.2 阅读 VASP 官网中关于原子和分子的例子，去尝试回答里面的问题；

https://cms.mpi.univie.ac.at/wiki/index.php/Atoms_and_Molecules

5. 总结：

分析完前面的测试，我们总结一下高效频率计算的关键点：

5.1 IBRION = 5 （告诉 VASP 我们要算频率）

5.2 POTIM = 0.015

5.3 NFREE = 2

5.4 ENCUT和原来一样

5.5 PREC = Normal

5.6 EDIFF = 1E-5 或者 1E-6

5.7 KPOINTS Gamma点即可。

本节内容太多了，原因在于，这是菜鸟篇的最后一节！！！

大师兄本人亲自做了很多测试，原因只有一个，告诉大家如何把握计算结果与时间的关系。在最短的时间获取最多的有价值的结果，这是计算化学的一个精髓所在。体现了你对时间和生命的尊重，也体现了你对高效率生活的追求。在菜鸟的时候，可以尽情地去测试。

大师兄之所以用乙醇这个大分子去给大家展示，很多人直接用 $\rm{O_2}$，$\rm{CO}$ 或者水分子的计算作为例子。那样的简单例子不具有很强的代表性，当你把一个相对复杂的东西搞明白了，这些小分子的计算，就如同小儿科一般了。

Ex26 乙醇分子的振动频率计算(四)

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

前面我们学会了通过使用 Jmol 查看分子振动。但分子振动频率在OUTCAR是什么样子的呢？今天我们就分析一下OUTCAR文件中的频率信息，以及如何写脚本计算零点能矫正。

1 OUTCAR分析

1.1 回顾一下Jmol中的频率振动

Jmol 提取了OUTCAR中的振动信息，将每个振动模式的频率列了出来。

1.2 OUTCAR中的信息：

对比一下这两个振动频率和·中的前两个。

1）1 f = 代表第一个振动模式，细心的你仔细观察，会发现每一行有四个单位的数值：THz, 2PiTHz, cm-1,和 meV，这四个是完全等同的；

2）下面一行为坐标 X Y Z和每个原子在 $x, y, z$ 方向上的振动大小；

3）X Y Z下面的数字为结构的坐标信息（Cartesian坐标系），dx dy dz 为振动的具体数值；

4)后面的振动模式的频率和第一个的格式一样。

2 频率单位的换算

我们先讲一下这四个单位的换算公式：

Relevant Formulas:
$ E = hc / \lambda $
$ \nu = c / \lambda $
$\tilde{\nu} = 1 / \lambda $
$ T = 1 / \nu $

Definitions:
$E$ = energy ($eV$)
$\lambda$ = wavelength ($m$)
$\tilde{\nu}$ = wavenumber ($m^{-1}$)
$T$ = period ($s$)
$\nu$ = frequency ($s^{-1}$ or $\textrm{Hz}$)
$h$ = Planck’s constant = $4.135667516 \times 10^{-15}~eV \cdot s$
$c$ = speed of light = $ 299792458~m/s$

以第一个振动为例：

1f =
111.907	$\textrm{THz}$
703.134	2PiTHz
3732.82	$cm^{-1} $
462.811	$meV$

1） THz 和 2PiTHz 的换算： $2\pi$ 的关系：

$111.907230~(\rm{THz}) \times 2 \pi = 703.133866~(\rm{2PiTHz})$

2） THz 和 cm-1 的关系：$\nu = c / \lambda = c \tilde{\nu} $， $c$ 是光速，$\lambda$ 是波长，$\tilde{\nu}$ 是波数，单位是 $m^{-1}$。计算前先换算成标准单位： 1 THz = 1012 Hz, 1 $cm^{-1}$ = 100 $m^{-1}$，因此：

$111.907230 \times 1012~(\textrm{Hz}) = 299792458~(m/s) \times 3732.823294 \times 100 (m^{-1})$

3）THz 和能量 eV 的关系：

$ E = h\nu $， $h$ 为普朗克常数：$4.135667516 \times 10^{-15}~eV \cdot s$ (注意此时：Plank constant 的单位！！！)

$462.811270 (meV) = 0.4628 eV = 4.135667516 \times 10^{-15} (eV*s) * 111.907230 \times 1012 (Hz)$

4）波数（$cm^{-1}$）和能量（$meV$）的关系：
$E = hc / \lambda = h \tilde \nu $
$ 462.811270~meV = 0.4628~eV
= 4.135667516 \times 10^{-15}~(eV*s) \times 299792458~(m/s) \times 3732.823294 \times 100~(m^{-1}) $

大家根据上面的公式自己手动算一遍就明白了，还可以使用下面这个网址进行计算：http://halas.rice.edu/conversions 这里就不再多说了。

3 OUTCAR频率信息的提取：

我们可以使用grep命令提取，有以下两个方式。在进行之前，首先强调一点：下面的解释只适用于NWRITE= 0、1或者2 的时候。因为NWRITE = 3 的时候，会额外再输出一次频率的信息。

3.1 使用以下几个命令：

grep THz OUTCAR
grep 2PiTHz OUTCAR
grep cm-1 OUTCAR
grep meV OUTCAR；

这几个命令中，我们分别以振动的不同单位作为提取对象，便可以得到所有的振动信息（这里的所有指的是包含虚频）：以大师兄常用的 grep cm-1 OUTCAR 为例：

黄色标出来的第一列：$9 \times 3 = 27$ 个振动模式，第二列是以cm-1为单位的振动频率大小，最后三行 f/i= 指的是虚频。

前面我们提到过，虚频可以判断结构是否稳定。那这里，我们计算出的乙醇分子结构肯定不稳定喽？不一定。

因为频率计算和软件的数值积分有关（我也不清楚数值积分怎么进行的）；

计算过程中我们的设置对频率计算影响很大，KPOINTS, ENCUT, EDIFF, POTIM等都会影响计算的精度（下一节讨论）；综合这些因素，对于分子的振动频率来说（注意：声子谱不适用）一般低于 $100~cm^{-1}$ 的频率可以忽略。严格点可以降到 $50~cm^{-1}$，也就是说：如果你在计算中发现有个 $50~cm^{-1}$ 左右的虚频，完全可以不考虑。

3.2 `grep 'f =' OUTCAR`

注意：图中虚频部分没有显示出来！严格按照我用的这个命令

使用这个命令的时候，不提取虚频部分。查看虚频的时候，可以用之前的方式，也可以用这个命令：

在零点能的计算时，虚频是不能考虑在内的，因为它不是分子的真实的振动模式。在我们这个例子中，虚频的出现是软件的误差所导致。在过渡态中，虚频代表的是反应方向。从另一个角度去分析：乙醇分子的零点能（下面讲到）为：$2.117~eV$，图中三个虚频对应的能量为：$0.76 + 2.31 + 8.14 = 11.11~meV = 0.01111~ eV$，所占比例为：$0.0111/2.117 = 0.5\% $ 这个可以忽略不计。

4 零点能校正

4.1 明白什么是零点能：回顾频率计算第一节的内容：

4.2 获取振动能量数据：

分析下结构：

上图输出共有11列（列之间用空格分开）：我们要的零点能在第10列，使用下面的命令：

1	grep 'f =' OUTCAR \| awk '{print $10}'

如果想同时输出第 1 和 10 两列：

1	grep 'f =' OUTCAR \| awk'{print $1 " "$10}'

$1 和 $10 之间有 2个双引号：“ “，两个双引号里面有一个空格用来分开），否则两列会连在一起。

注意！注意！注意！

这里我们提取的能量为：$h\nu$ ！！！

而零点能为 $1/2 h \nu$！！！

4.3 将所有振动的能量求和：

1 2	qli@BigBroSci:~/Desktop/freq$ grep 'f =' OUTCAR \| awk '{print $10}' \| paste -sd+ \| bc qli@BigBroSci:~/Desktop/freq$ 4233.962325

输出的4233.96就是所有的 hv 之和。

1）不要忘记除以 $2$

2）此时单位是 $meV$，换算成 $eV$ 还需要除以 $1000$。

3）所以，我们的零点能是 $4233.962325/2/1000~ eV= 2.117~eV$

4.4 写脚本

将前面的命令写到一个文件里面就成为了脚本：怎么写呢？本人平时喜欢使用类似 $4233.962325$ 这个数字，除以 $2000$ 这一步在后面的工作中进行。因此只用到了

1	grep 'f =' OUTCAR \| awk '{print $10}' \| paste -sd+ \| bc

这个命令。写脚本的具体操作：

1）把这行命令写到一个文件中，文件名为 fsum （本人不喜欢 .sh 尾缀，加不加都一样）；

2）赋予可执行权限

1	chmod u+x ~/bin/fsum

3）移到 ~/bin 文件夹

1	mv fsum ~/bin

4）在任何一个频率计算的目录下运行：敲命令 fsum 即可；

5）不喜欢 fsum 这个命令，将文件名改成你自己喜欢的名字：

如果想在脚本里面直接完成除以 $2000$ 的任务，可以这么写：

1 2	hv_sum=$(grep"f =" OUTCAR \| awk '{print $10}'\| paste -sd+ \| bc) echo "scale =6; $hv_sum/2000" \| bc

此时结果的单位为：$eV$。

5 扩展练习：

5.1 OUTCAR中的频率输出要看明白是怎么回事；

5.2 频率的各个单位的换算要搞明白；

5.3 怎么提取信息，计算零点能要掌握；

5.4 怎么写脚本，从现在开始要练习了！

6 总结：

本节我们主要讨论了一下OUTCAR中的频率结果，能量换算，如何提取振动能量，以及如何计算零点能，最后简单介绍了一下脚本的写法，打消大家对脚本的崇拜心理，自己稍加琢磨也会写出实用的脚本！

Ex25 乙醇分子的振动频率计算(三)

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

本节推荐一款可视化程序：Jmol，可以用来看分子结构以及振动频率。Jmol 是一款Java语言编写的，开源，Linux，Windows均可的使用的可视化软件。还记的我们从ChemSpider获取乙醇分子的情景吗？

对了！RSC的网页版中，乙醇分子的3D结构就是通过 Jmol 展示给大家的（3D图上方），并且我们下载的结构也是jmol格式的文件！可以通过 Jmol 直接打开。现在我们通过在电脑上运行Jmol，并查看分子的结构和振动频率。

1. Jmol在windows和linux下的安装：

1.1 安装前准备（Java运行环境的安装）

Windows下：直接百度Java Runtime Environment（JRE），安装软件。

也可以官网下载：http://www.oracle.com/technetwork/java/javase/downloads/jre8-downloads-2133155.html

注意两点：

1：接受许可，

2：64位系统下载箭头所指的文件，然后安装即可。

Linux下Java运行环境安装移步下方链接：

怎样在Ubuntu 14.04中安装Java： https://linux.cn/article-3792-1.html

或者直接使用命令：

1	sudo apt-get install default-jre

1.2 Jmol软件的下载和安装

Windows：

https://sourceforge.net/projects/jmol/files/Jmol/

打开链接后，图中箭头指的地方下载最新版的（Linux其实和Windows下载的文件一模一样，见后面说明），也可以任选版本进行下载。

等待几秒后会弹出下载的窗口，其他浏览器也应该一样。

解压缩后就算安装好了，可以直接运行。

运行 Jmol 程序：

上图中 箭头1所指为 jmol.bat 文件，箭头2所指为 Java的可执行文件：jmol.jar。在Windows下任选一个双击即可打开 Jmol 程序。建议将解压缩之后文件夹中的jmol.bat或者jmol.jar 文件右键发送到桌面快捷方式。打开后如下图：

Linux下安装 Jmol 的方法:

在 https://sourceforge.net/projects/jmol/files/Jmol/下载 Jmol（和前面介绍的Windows下载的文件一模一样），解压后的文件复制到usr/bin中，打开文件时在命令行输入jmol.sh 文件名即可（注意，jmol.sh和文件名之间 有空格！）。将解压目录加入.bashrc路径中也可以实现，同时需要给jmol.sh加上可执行权限，即chmod u+x jmol.sh，打开文件的方法仍用jmol.sh 文件名。(QQ群友：连赞提供！)

大师兄在Linux系统下的操作如下：

a）下载和Windows版的过程一样，下载完毕后，解压，

b）终端里面进入解压后的文件夹：

注意一下几点：

b.1) jmol.sh就是我们在Linux系统下面的命令（Windows里面我们用jmol.bat 和 jmol.jar）；

b.2) 首先赋予它可执行的权限：

1	chmod u+x jmol.sh

相反地，取消权限：

1	chmod u-x jmol.sh

b.3) 这里大师兄电脑下，jmol.sh变成绿色的了（不同电脑显示会不同，不要纠结）

b.4) 尝试运行一下： ./jmol.sh OK

b.5) 下面设置环境变量：打开 ~/.bashrc 文件，并加入这一行：

1	export PATH=$PATH:~/Downloads/jmol-14.20.2:$PATH

注意：

等号=前后没有空格, 后面紧跟着$PATH；
再往后是一个冒号，冒号前后也不能有空格，填上Jmol的解压目录！！！

b.6) 保存.bashrc 文件并source一下: . ~/.bashrc （注意前面的 . ）或者使用命令：

1	source ~/.bashrc

c）进入其他目录，运行 Jmol 命令：下图中~/Destkop/freq 目录下有我们关于乙醇频率的计算结果。

大功告成！！！

此外，sudo apt-get install jmol这样安装的是旧版本，强烈不建议，原因如下：

可能图片有些模糊，但不重要！旧版本的不支持 OUTCAR…安了也是白搞！！！

2 使用Jmol可视化分子振动

2.1 载入振动文件到 Jmol

Windows：直接将频率计算得到的OUTCAR拖到jmol.bat就可以了

Linux：直接 jmol 文件名

两个系统下面均可使用左上角的文件à打开à选择OUTCAR导入。

得到如下界面。但是分子并没有开始振动，只是显示了其结构。

（该结构和ChemSpider上的一样！）

要查看振动模型，需要以下两步，选中振动模型和开启振动。这两步有多种方法可以实现，总结如下。

2.2 选择振动模型

A）可以在工具—>原子库选择器中选中要查看的振动模型

注意，先把右边的按钮拉倒底，然后双击Frequencies展开频率信息， 如下:

B) 也可以直接右键—>模型中选中

2.3 开启振动

A) 原子库选择器中最下方: 振动—>振动开 可实现，点击后，原子就开始振动了。

B) 在菜单栏—>工具—>震动 中开启（原子库选择器选中时采用这种方式，需要先开启再选中）

C) 也可以右键 —> 振动 —> 开启

2.4 查看不同振动

A) 原子集选择器中的频率列表，双击其中的一个就显示其振动方式了。

B) 要查看下一个振动还可以通过点击向右的这个箭头就可以了。

以上就是本节关于 Jmol 可视化的介绍了，更多信息，请查看 Jmol 的wiki：
http://wiki.jmol.org/index.php/Main_Page

3 扩展练习：

3.1 熟练掌握 Jmol 的软件安装，查看频率的基本操作。

3.2 通过在Linux系统下面安装 Jmol，尝试着安装其他软件；

3.3 分析各项频率值对应的分子或者原子的移动。

4 总结

本节主要介绍给大家一款除了 p4vasp 之外的另一款 VASP 可视化的软件，不仅仅局限在频率振动分析方面，查看POSCAR，CONTCAR也可以直接使用 Jmol 打开。该软件由连赞小朋友推荐，并编辑文中大部分内容，在此表示衷心的感谢。希望大家都可以熟练掌握这款软件。Linux下面可以直接用命令： jmol 文件名，打开查看结构，非常方便。

Ex23 乙醇分子的振动频率计算(一)

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

大师兄在本节开头放上这张图片的意思是，请大家再仔细体会Rubbish in， Rubbish out这句话，只有对我们要研究的体系有一定的理解，明白我们具体要计算什么内容，设置好输入文件的参数，才能得到我们需要的合理结果。接下来的几节会涉及到分子的振动频率计算，振动频率的可视化等内容，请大家跟着大师兄一起练习，掌握振动频率的基本计算和分析方法。

1 分子的振动

我们首先看回顾一下振动相关的基本知识，这里大师兄不具体解释，引用 $9$ 版 $Atkins$ 的《物理化学》书中的内容，书已上传至QQ群文件中，也可百度网盘下载: http://pan.baidu.com/s/1o8HlyOi 全是英文，大家耐心阅读下。

简谐振子、胡克定律、体系势能随着振动距离 $x$ 的关系

简谐振动的薛定谔方程描述:

薛定谔方程的解：振动的量子化，振动频率

零点能的数学和物理两个方面的解释:

2 频率计算的作用

频率计算有什么用？为什么要算频率? 大师兄稍微总结了一下频率计算的意义，大体有以下几个方面，没有提到的用处，烦请大家指出来，以便补充。

2.1 确定结构是否稳定;

2.2 看振动方式和大小，用来和实验对比，棋博士最新的文章就是一个非常好的例子;

2.3 反应热，反应能垒，吸附能等的零点能矫正;

2.4 确认过渡态(有一个振动的虚频)

2.5 热力学中计算entropy，用于计算化学势，微观动力学中的指前因子和反应能垒。

3 怎么用VASP计算频率?

3.1首先进行结构优化，获取稳定的构型，这个我们前面已经讲过了;

3.2 将原来的CONTCAR复制成POSCAR :

1	cp CONTCAR POSCAR

3.3 修改INCAR

修改后如下:

频率计算的INCAR

IBRION的值改成5
POTIM用一个更小的值，我们这里用的 0.02，默认值是 0.015
NSW 设置成1，这个可以直接不管，继续采用优化时的NSW值，因为你设置成 1， 2， 3， 4， 5， …， 1000 都不会影响计算；但不能不设置（因为默认值是0，这时算个单点后任务便停止了。）
NFREE=2 添加这一个参数，表明原子在某一方向上正反两个方向移动；
NCORE=4这一项要注释掉！大师兄这边的服务器，并行计算频率时 VASP 会罢工，只进行一步静态计算，注释掉就正常进行了；
此外，EDIFF也要设置一个严格的值（频率计算时，默认值为1E-6，足够了！下一节会讲到）

小结一下频率分析关键的参数:

1
2
3

IBRION=5
NFREE=2
POTIM=0.02

4 扩展练习

4.1 按照本节的流程新建一个文件夹 freq : 该文件夹中包含乙醇分子优化后的结构 (将CONTCAR复制成POSCAR)，以及优化时的POTCAR，INCAR， KPOINTS以及提交命令的脚本文件;

4.2 修改乙醇分子优化的INCAR为频率计算的INCAR；需要修改哪些参数心里要清楚；

4.3 运行乙醇分子频率计算，并查看频率分析的OUTCAR，OSZICAR等输出文件；

4.4 查看 VASP 官网对于IBRION=5 的解释，搜索网上相关频率计算的文章，帖子，初步了解NFREE，POTIM所代表的含义；

https://cms.mpi.univie.ac.at/wiki/index.php/IBRION

4.5 查找官网中频率计算的例子:

A) https://cms.mpi.univie.ac.at/wiki/index.php/CO_vibration

B) https://cms.mpi.univie.ac.at/wiki/index.php/H2O_vibration

5 总结:

5.1 熟悉频率计算初始文件的准备过程；

5.2 频率计算INCAR中的三个重要参数；

1
2
3

IBRION = 5
NFREE = 2
POTIM = 0.02

5.3 初步了解频率计算中各个参数的含义。