LVASPTHW

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

大师兄在本节开头放上这张图片的意思是，请大家再仔细体会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 初步了解频率计算中各个参数的含义。

LVASPTHW

Ex22 乙醇气相分子的优化

前面我们终于讲完了 $\rm{O_2}$ 分子优化的例子。相信大家对 VASP 计算已经有了一个初步的理解。这一节我们继续学习气相分子的优化。为了让大家进一步了解计算的过程，我们选取一个稍微复杂的分子作为例子：乙醇（$\rm{CH_3CH_2OH}$）。

问：把大象装进冰箱需要几步？

大师兄，这个问题我早就知道答案了。为啥还问这样低智商的问题？

大师兄要求的不是让你回答开冰箱门，装进去，关门的这三步。而是让你尝试回想一下：当你第一次接触这个问题的时候，你的反应是什么？

大师兄比较笨，我的第一反应是，这怎么可能？ 冰箱那么点，大象那么大。

当朋友告诉我答案的时候，才恍然大悟，这跟大小没关系。

同样的，怎么用 VASP 计算乙醇分子？很多童鞋就如同第一次被问到大象这个问题一样不知所措。答，也是三步！

1） 打开冰箱：准备 VASP 文件

2） 把大象塞进去：准备乙醇分子模型

3）关上冰箱：运行 VASP

1. 打开冰箱：

我们可以直接用O$_2$分子计算的输入文件，

1.1 复习一下前面学到的INCAR和KPOINTS的内容：

1）乙醇分子是闭壳层的分子，没有磁性，不需要ISPIN=2

2）气相分子计算，我们要用ISMEAR=0，SIGMA取值要小，SIGMA=0.02；

# 表示注释，这个符号后面的内容，VASP在运行的时候不考虑。

示例 1 ：

1
2
ENCUT = 400
# ENCUT = 500

和
1
ENCUT = 400

效果是一样的。

示例 2：# 的用法：

VASP 计算中一些常见的错误，以及注意事项，你可以通过 # 写在INCAR里面，方便计算的时候进行设置。可以 # 开头，单起一行，也可以在参数的最后面加上注释。新手们刚刚开始，可以结合 VASP 官网参考书，把用到的INCAR参数注释下来，时间长了慢慢就掌握了。

3）气相分子计算，K点使用Gamma点就够了；

1
2
3
4
5
K-POINTS
0
Gamma
1 1 1 
0 0 0

2 把大象塞到冰箱里面 :POSCAR的准备（本节重点）

前面我们学习到，一个好的初始结构会加快计算，获取准确的计算结果。因此我们需要去找一个合理的模型。对于气相分子的结构，一是手动搭建，大师兄推荐用 GaussView，另一个办法就是找数据库，大师兄推荐英国皇家学会的 ChemScpider。网址：http://www.chemspider.com

下面我们把手动搭建乙醇分子的模型具体解释一下：

2.1 打开网站

搜索框中输入： 分子名称Ethanol 或者 分子式C2H6O…. 点击搜索，等待结果：

上面箭头指的两个地方随便点，效果是一样的，如下：

分子默认显示 2D 的结构，点击箭头指的 3D，切换。第二个箭头所指的有对应结构的维基百科链接，翻墙的筒子们可以查阅下相关结构的知识。选择3D后，如下：

点击箭头指的地方保存，浏览器会下载对应的 .mol 文件，文件名为该结构在数据库中的编号。

.mol文件也是纯文本，使用Notepad++ 打开如下：看到这么多内容，不要害怕！

注： 不仅仅是 .mol 文件，很多结构文件都是文本格式，直接打开就是。

从里面找到乙醇结构的 $x, y, z$ 坐标信息

删除上图中除红色框之外的所有行，1-4行，以及 14 至最后一行，只保留 $x, y, z$ 坐标信息的那几行。

手动写POSCAR，首先你要熟记的格式，知道POSCAR从头开始，往下每一行代表的内容。

我们的模型是把乙醇分子放到一个 $ 20 \times 20 \times 20~A^3 $ 的格子里面。输入完之后，坐标后面的元素符号以及那些 $0$ $0$ $0$ 可以删掉，也可以不管。

关键点1： .mol文件中，坐标那几行中的第 4 列写到POSCAR中的第 6, 7行！

关键点2： 注意坐标为 Cartesian 或者是 Direct. Direct 是分数坐标，其 $x, y, z$ 值都小于 1.

如果你想删掉图中坐标第3列后面的内容：大师兄推荐notepad++里面的列块模式：如下：

点击之后，可能会弹出对话框，告诉你如何使用，列块模式

方法1：

摁住 Alt 键，然后用鼠标选择文本，不同电脑可能不一样，大师兄这边同时摁住 Ctrl 和 Alt 两个键，然后用鼠标选择的。

方法2：

同时摁住 Alt和Shift键，通过键盘上前后左右的箭头选择文本

大家可以尝试下，选中效果如图：

然后点键盘上的Delete键删除.

然后另存为POSCAR即可。

我们可以使用 p4vasp 来查看一下模型的结构：如下图

我靠，师兄，结构怎么跑到格子外面啦？前面 $\rm{O_2}$ 分子的学习中，你已经知道了这是因为周期性导致的显示问题。对计算不会产生影响。这个结构可以拿来直接用。

如果感觉不爽，想把结构放到中间，可以这么做：

数学上，把 $x, y, z$ 坐标统统加上 $10$ 即可；

软件使用上，我们讲一下 p4vasp 的操作方法：

选择 edit –> Move atoms

Move group 是你要移动的原子，这里大师兄直接输入了 C H O 三个元素符号（中间有空格）表示选择所有元素的原子。然后在 Vector 中选择 $x, y, z$ 三个方向上移动的大小。你也可以写 $1$ $1$ $1$，然后点击 Move 按钮 $10$ 次。

如果你想通过选择原子来实现移动的话(不直接在 Move group 里面输入 C H O)，需要按照大师兄说的步骤走：

1) 空格键结合鼠标选中所有的原子；鼠标指到原子上就点一下空格键

2) 选中原子后，主界面会显示一些数字，这些数字和 POSCAR 中元素的顺序是一致的；

3) 所有原子选中后，左下角的框中点击 Get group，会显示选择的那些原子；

4) 在 Vector 中选择 $x, y, z$ 三个方向上移动的大小，然后点击 Move 按钮。

上面一堆废话就此打住，结果就是这样子的：

然后点击 file -> Save system as -> 选择目录 -> 保存成POSCAR

POTCAR

POSCAR 讲完了，我们就要按照里面的元素顺序制备POTCAR了。

首先: 我们要准备 O C H 三个元素的POTCAR，去POTCAR的数据库中去找：

然后复制到当前目录下，三个元素的POTCAR分别命名为：POTCAR-O, POTCAR-C和 POTCAR-H， 把三个元素的POTCAR合并在一起，命令就是

1
cat POTCAR-O POTCAR-C POTCAR-H >> POTCAR

3 提交任务

提交任务之前，需要再次检查自己的输入文件一遍，没有问题，提交直至结束。

4 扩展练习：

4.1 从头到尾，认真重复本节中大师兄的操作；

4.2 记住本节讲解的内容，自己重复一遍操作，直至自己通过文本编辑器会搭建结构模型；

4.3 运行乙醇计算的例子；

4.4 尝试使用 GaussView，Materials Studio (MS)，以及其他可视化界面搭建模型；

4.5 学会使用 VESTA 导出POTCAR格式的结构数据。

5 总结：

本节，我们主要讨论了一下分子模型的数据库搜索和搭建工作。希望大家能够完全掌握本节的所有内容和细节。在计算中，你会遇到各种各样的结构文件，其实都是$xyz$坐标的衍生物而已。不要害怕，直接打开它们，学会提取里面有价值的信息。

LVASPTHW

Ex19 谁偷走的我的机时？（四）

模型对计算时间的影响

上一节介绍的KPOINTS对计算的影响，相信大家已经认真阅读参考书的第三章部分了。本季我们讨论一下模型的大小对计算的影响。主要体现在晶胞的尺寸，对称性以及对K点的影响上。

1 测试工作：

为了方便处理，我们把O$_2$计算的格子设置为长宽高均为8.0 $\AA$。


重复之前KPOINTS的批处理操作，我们可以获得一系列不同大小格子的文件夹。如下图：

命令： for i in $(seq 10 2 20); do cp 8 $i; sed –i "3,5s/8.0/$i/g" $i/POSCAR ;done

2 测试结果分析

2.1 模型大小对计算时间的影响

注意：在后面加入 sort –n 后输出的变化。

从图中可以看出来，计算时间随着格子的大小，需要的计算时间增加的很快。

注意：在测试中，KPOINTS一直保持不变（因为只有一个Gamma点）。而在我们实际的计算操作中，使用1x1x1 KPOINTS的机会并不多。如果格子在某个方向增加了2倍，那么对应的改方向的K点就需要除以2。重复一下上节的经验指导。也就是在计算过程中，保持ka保持不变。当然，ka是我们提前测试好的。

举例：

一个10x10x10 $\AA^3$的体相材料，我们计算的时候K点设置为：6x6x6。

当我们将材料在x方向增加1倍，变为20x10x10 $\AA^3$。为保持一致的精确度，那么我们的K点需要设置为：3x6x6。

这是因为倒易晶格矢量和实际的晶格矢量之间存在着倒数的关系：

注： 类似的图，不加说明，均出自我们的参考书。也就是说，我们选取的晶格越大，倒易晶格矢量越小。用同等数目的K点分布到倒易晶格中，网格的密度也会越大，从而造成计算量的增加。

2.2 体系的对称性对计算速度的影响：

2.2.1 K点保持不变：

这一点前面关于氧原子的计算就已经介绍到了，降低体系的对称性会增加额外的计算时间。如图：将12x12x12 $\AA^3$ （计算需要186.86 s）的格子修改如下：

计算结束后，查看时间，为194.9 s, 计算时间增加了8秒。

2.2.2 对称性对K点的影响：

体系的对称性不仅仅提现在前面的计算中，更可以在计算中极大地减少K点的数目，从而加快计算，节省时间。这一点我们引用参考书中的一段话：

2.3.3 模型对称性与K点对称性的关系

在这里，体系的对称性与K点对称性的匹配问题，尤其是对于hexagonal的结构来说，必须要使用 gamma centered points. 也就是第三行的第一个字母必须为G或者g。我们看一下官网的原话：

We strongly recommend to use only Gamma centered grids forhexagonal lattices. Many tests we have performed indicate that the energy convergessignificantly faster with centered grids than with standard Monkhorst Pack grids. Grids generated with the “M” setting in the third line, in fact do not have full hexagonal symmerty.

如果你不确定自己的体系，直接用G就可以了。

For reasons of safety it might be a good choice to use only meshes with theirorigin at (switch “G” or “g” on third line or odd divisions) if the tetrahedron method is used.

3 扩展练习：

1 认真阅读： Density Functional Theory: A Practical Introduction: 第三章的前两节；

2 VASP官网查找K点相关的说明。

4 总结：

学习完本节，大家应该掌握的内容有：

4.1 晶格大小对计算时间的影响；

4.2 体系的对称性对计算时间的影响；

4.3 掌握K点和晶格大小的经验规则；

4.4 晶格对称性和K点对称性的一致性。

LVASPTHW

Ex20 谁偷走的我的机时？（五）

并行

前面几节讨论的都是一些涉及到模型以及计算细节对时间的影响。本节我们讨论一下服务器节点设置（并行）的影响。从字面上不难理解，并行就是多个节点同时计算同一个任务。好比之前用一匹马拉车，现在改用两匹，三匹或者更多的马拉同一辆车。一匹马拉车的时候，对马的要求是足够强壮，能拉得动还要跑得快。但你想跑的更快，就需要驾驭两匹或者多匹马，但这个时候对于驾驶马车的你就需要提出技术要求了：如何控制马儿之间的节奏。控制好了，并驾齐驱，得儿驾得儿驾，爽歪歪。控制不好结构，你跑你的，我跑我的，整体下来，马儿也累，马车行驶的反而更慢了，搞不好还会栽跟头。

同样的道理，如果想加快计算速度，或者在最短的时间内获取最多的计算数据，我们就需要知道并行在计算中的作用。首先我们先测试一下，不同节点同时运行O$_2$分子计算所需要的时间。注意：本练习讨论的是多个节点下的并行！

1 测试活动：

为了观测更加明显，我们使用8x8x8 $\AA^3$的格子，K点使用3x3x3， ENCUT = 400 。

测试的服务器每个节点有4个核，每个核有两个进程。

文件夹用核数命名，4,8,12,16 分别代表使用了1,2,3,4个节点进行计算。

在这里，每个节点可以看成1匹马，4个核可以看作马的四条腿。每个核的两个进程可以看做每条腿的上下两部分。

1.1 设置计算的节点数目

使用多少个节点进行计算，在提交任务的脚本里面设置。每个课题组可能不太一样，也有通过命令设置的，不过都大同小异。大师兄晒一晒自己组里提交命令的脚本，如果和你们的不一样，不要纠结：

这个任务的脚本里面，我们使用32个核（32条腿），也就是8个节点（8匹马）同时计算。提交所有的测试任务后，等待结束。

2 测试结果分析：

2.1 提取结果

for i in *; do echo -e $i "\t" $(grep Elapsed $i/OUTCAR | awk '{print $4}'); done

这里复习上一节： sort –n 命令的用法，将数列从小到大排列。

注意：这里我们没有用之前的User time提取时间, 而是用的Elapsed time。这一项是计算真正花费的时间。所有的时间里面：Total CPU time = User Time + System time，一般来说， Elapsed time 总会比 Total CPU time 多上那么几秒。

User time 和 System time 是干嘛的？怎么的出来的？大师兄也不知道具体的含义。希望知道的童鞋们给大师兄上上课。之前几节中，用User time的分析结果可靠吗？ 没问题的，大师兄检查了一下，前面的System time均为 2s左右，Elapsed Time 中绝大部分由User time来贡献。

2.2 作图分析

核数的单位为 个 ， 对应图中的拼音 ge（个数）

1 从4核到8核，计算时间并没有减少，反而增加了15s

2 从4核到12核，计算时间也仅仅减少了20s

3 从12核以后，增加核数，计算时间反而增加了。

从12核开始，增加的核数白白占了你的任务，缺丝毫没有提高计算速度。因此，增加核数反而成了计算的负担。这是因为并行计算的时候，不同核之间的数据传输浪费了大量的计算时间。这好比是你多安排了几匹马来拉车，但马儿们却彼此交流，各跑各的，不愿意拉车了。

2.3 既然是这样，VASP并行还有个卵用？

答：并行是有用的，只不过我们需要调教这些不听话的马儿们，让他们服从我们的命令。这里我们就会需要一些其他的参数，NCORE和NPAR。

NCORE：控制多少个核同时计算；
NPAR：如何把计算任务分配到计算资源上面计算。

它们之间的关系是：NCORE= 计算使用的核数 / NPAR

注意：这两个参数只能选取一个来使用:

测试数据结果（给喜欢看表格的人准备的）

图中红色部分为前面计算的结果，

  • 蓝色为NCORE=4的结果（测试的服务器每个节点有4个核，所有的计算军用NCORE=4）；
  • 绿色为NPAR=4的结果（所有的计算中均用的NPAR=4）；
  • 黄色为NPAR=Core/4的结果，

根据NPAR和NCORE的关系，该设置与NCORE=4等同，但是实际测试中蓝色和黄色为什么没有完全一样，我也说不清楚。

1）使用NCORE以后，单节点运行也加快了；（单匹马儿被你调教的更听话了）

2）加入NCORE和NPAR参数后，计算时间明显提高了。在20核以后的计算中尤为明显；

3）如果你想用多个节点计算，NCORE或者NPAR，不要忘了加（二者选其一）。

４）可以肯定的是：NPAR和NCORE的乘积就是我们计算所用的核数，这一点大家要记在心里。

5) 在本测试中（这是前提！！），增加节点个数并未实现 1+1 = 2 的效果，1+N = 2 的效果也没有实现。但在多个节点的计算中，加入NCORE或者NPAR，节点数越多，效果越明显。

手册上解释说NPAR，或者NCORE的取值可以为总核数的开方值。很多人对此会有疑虑，大师兄的建议是：

  • 测试，测试，测试!….测试一下自己的体系这是主要的，不要完全相信本测试的结果，体系不同，结果可能会差很多！！；

  • NPAR的取值可以设置为节点的数目，（默认值为计算的核数！）

  • NPAR实在不懂的话，直接设置NCORE=单节点的核数，单节点的核数/2，单节点的核数/4…….

  • 个人使用经验是：NCORE = 单个节点核数 / 2 的时候，运行最省时间，设置也最方便。

3 扩展练习：

1 阅读VASP官网关于NCORE，NPAR的内容；

2 对于自己的计算体系，测试这两个参数；

3 总结前面关于计算时间的影响因素

4 总结：

关于VASP计算时间的影响，我们暂且告一段落了，希望大家能在前面几节的学习中熟练掌握影响计算时间的这些因素，在实际的计算过程中合理把握，提高计算效率，节约机时。浪费的机时不仅仅属于你自己，也属于组里其他人的。

LVASPTHW

Ex18 谁偷走的我的机时？（三）

KPOINTS对计算时间的影响

继续前面的学习，本节讨论KPOINTS文件中K点的设置对计算时间的影响。本节图中的Linux命令不再详细介绍，大师兄默认大家已经基本掌握了其中的原理和窍门。VASP官网上还有其他的批量测试的脚本，大家现在也可以差不多能看懂了。比如Si计算的例子：https://cms.mpi.univie.ac.at/wiki/index.php/Fcc_Si

1 KPOINTS测试

1.1 准备测试模板

新建文件夹: 3，将之前的文件夹0复制过来后重命名为 1。

1.2 批量制备测试文件：

命令： for i in {2..6}; do cp 1 $i ; sed -i"s/1 1 1/$i $i $i/g" $i/KPOINTS ; done

文件夹1 代表KPOINTS为1 1 1 ,6 代表KPOINTS为 6 6 6 其他的类推。

1.3 批量提交任务：

备注：这里的qsuball 命令是把前面的批量命令放到 .bashrc 文件中了。不懂的请看Ex17的批量操作命令和Ex13中.bashrc文件中alias的使用方法。

2 测试结果分析

2.1 查看OUTCAR 中的K点信息

图中，我们找出可以通过grep 查询的字符：irreducible

2.2 批量查看所有测试的K点信息：

命令 grep irreducible 3/OUTCAR

从图中可以看出：

1）K点2 2 2 和3 3 3 的计算中，生成的K点数目是一样的；类似地，4 4 4 和5 5 5 具有同样地K点数目；

2）不难理解，相同的K点数目，其计算时间也是一样的；

3）计算时间随K点数目的增加也增加了。

师兄，为什么K点数目会存在奇数和偶数相同的关系？

原因在于K点生成的方法。当KPOINTS为偶数的时候，K点都在布里渊区的内部，而为奇数的时候，部分K点处在布里渊区的边界上。引用参考书中的一个表格和一段话：第三章第56-57页：

2.3 不同K点对能量的影响：


先不画图了，直接看能量吧：

当K点为 1 1 1 和 6 6 6 的时候，O$_2$分子的能量差别为： 0.0015 eV。可以忽略不计。在这里，你可以清晰地知道：为什么算气相分子的时候 gamma点（1x1x1）足够了。

对于其他slab或者体相材料的计算，K点怎么选择呢？我们看下面这一段话：

出处： https://wiki.fysik.dtu.dk/gpaw/exercises/surface/surface.html

再参考一下Quantumwise 网站的说明：

出处：http://quantumwise.com/forum/index.php?topic=2628.0

再次强调一下：浏览网站说明的时候，要养成这样的一种习惯，凡是看到Note这个单次，就要跟打了鸡血一般！前面的东西看不懂不要紧，Note后面跟的都是重点易出错的地方。

以上只是经验参数的说明，给我们提供一个大体的指导。这个参数在使用中，要注意我们前面提到的奇数和偶数的情况。而具体到我们的计算中，需要用什么数值，我们需要认真地测试检查一下，而不能直接就用图中的经验参数。

1）通过测试不同K点对体系能量的变化；（参考书中的例子）

2）查找参考文献的取值；

3）此外，不同K点之间的数据不能混用。比如计算CO在一个 (3x3) Cu(111)表面上的吸附能：

等号后面的前两项，必须要用同一个K点下计算出来的能量，如果$E{CO+slab}$ 用5x5x1 的K点， E(slab) 采用 3x3x1的K点能量，得出的结果必然是错的。

2.4 来自VASP官方的提醒：

出处： https://www.vasp.at/vasp-workshop/slides/accuracy.pdf

1）常见错误（一）： 体系中ENCUT的取值不统一；

2）常见错误（二）：采用不同KPOINTS计算出来的结果。

3）在关于Accuracy的这个pdf文件中，最后一行大家要谨记：TEST，TEST，TEST ….

3 扩展练习：

1 下载大师兄分享的压缩文件： 本节中所有的pdf文件和链接；大师兄QQ群文件下载，或者百度网盘：http://pan.baidu.com/s/1eSCGWeA

2 阅读参考书中第三章的内容，掌握K点的基本概念和一些选取的注意事项；

3 浏览本节中所有的网址。

4 总结：

1 K点数目越多，计算越准确，需要的时间也会相应地越多，大家要把握好准确度和时间的关系；

2 K点的确定，需要经验和测试相结合，经验为辅，测试为主。

5 补充：

VASP的新版本中，可以直接在INCAR中设置K点。也就是说，没有KPOINTS文件也可以正常计算，但这可能只适用于简单的K点情况，能带结算等需要制定K点路径的计算，我们还是需要KPOINTS文件的。

LVASPTHW

Ex17 谁偷走的我的机时？（二）

这一节和后面几节，我们会复习一下前面学到的Linux批量操作知识，然后对一些影响计算时间的参数进行测试。 本节主要考虑ENCUT的影响。前面我们一直在说EDIFF和EDIFFG对计算时间的影响。为了给大家一个感性的认识，现在我们用O$_2$分子的计算作为测试例子，将结果展示出来。这两个参数主要是在计算精度上影响计算时间，很容易想到，精度越高，收敛的越慢，需要的时间对应的也会更长。

1 调节EDIFF和EDIFFG

左侧为之前计算的INCAR，右侧为提高精度后的INCAR。O$_2$的初始距离设置的为1.207$\AA$。提交计算，等待任务结束。前面我们知道VASP计算完成后，OUTCAR最后输出的是计算时间，内存等信息。我们现在查看一下：文件夹 0 对应的是之前的O$_2$计算，文件夹1 中是提高精度后的计算。

注意：grep User OUTCAR 后得出的结果被空格分成了4部分，时间信息在第4部分里面。User（1） time（2） （sec：）（3） 44.20 （4）

讲解：

1.1）通过tail OUTCAR 这个命令，可以找到用grep命令查看时间的关键词：User 或者 Elapsed，这里我们采用User 后面的时间作为参考；

1.2）右下方黄色图框中的时间表面：提高了精度后，计算时间从44秒增加到78秒。

1.3）我们看一下，改变精度后体系的能量变化：

从-9.8609降低到-9.8611，变化大小为： -0.0002 eV。 这么小的能量变化，我们可以认为忽略不计。

在这里，我们要认真思考收敛标准对于我们计算体系能量的影响，选取一个合适的标准而又不会浪费太多的机时。一般来说，结构优化的时候，EDIFF=1E-5, EDIFFG =-0.01至 -0.03 都是被认可的。

2 ENCUT 测试：

2.1 制备ENCUT测试模板

解释：

2.1.1）新建测试目录 2 ，进入后，将前面的文件夹0复制过来，文件夹名为 400

2.1.2）使用sed命令，在INCAR中最后一行加入ENCUT参数，值设置为400

2.2 快速制作测试任务

图中的命令行为：

1
for i in {1..8}; do cp 400 $((400+$i*50)); sed -i "s/400/$((400+$i*50))/g"  $((400+$i*50))/INCAR ; done

复习前面学到的linux操作：

2.2.1）明白{1..8} 是怎么回事；

2.2.2） $i 变量的调用；

2.2.3）新学：$((加减乘除))，注意数学运算用2个括号 括起来；

2.2.4）sed 命令进行文本中某一项的替换；

2.2.5）运行完毕后，会获得一系列ENCUT值的文件夹，且每个里面INCAR已经对应地设置完毕。

2.3 批量提交任务


命令：
for i in *; do cd $i ; qsub sub4; cd $OLDPWD; done

讲解：

2.3.1）do后面执行的是，进入for循环中的文件夹，然后提交任务，（大师兄提交任务的命令是 qsub sub4， qsub 是命令，sub4是脚本名），任务提交后，返回原来的目录下（cd $OLDPWD），然后再进入下一个for循环中的文件夹，重复之前的操作，直至遍历所有for循环的变量文件夹；

2.3.2）图中的Single是任务的说明，这个在提交任务的脚本里面自己随意设置：

2.4 批量查看结果的命令：

此处需要注意 * 的用法

我们只需要图中方框标出来的信息，其他的结果可以直接扔掉。为简化输出结果，这里大师兄用了另外一个强大的命令： awk

for i in *0; do echo -e $i "\t" $(grep User $i/OUTCAR | awk '{print $4}'); done

2.4.1） 学会echo 命令：

2.4.2） echo –e 后面加上”\t”后（双引号），可以直接输出 tab，方便导入excel；

2.4.3） 本例中，awk 后面用{}将打印的内容括起来；

2.4.4） print $4 意思是输出前面结果的第4项；

2.4.5） awk命令极其强大，强烈建议大家尝试着去网上查找资料，主动去学习；

2.4.6）将结果复制到excel里面作图。

2.5 计算时间随着ENCUT增加的变化曲线

从图中可以看出，计算时间随着ENCUT的增加也相应地增加了。因此，在保证计算准确度和ENCUT的值的选取，你要学会合理取值，加快自己的计算速度。此外，图表的比较形象直观，大家在今后的学习中，多多思考将自己的数据转化为图表的形式！可以很好的表达自己的计算结果。

设想一下，同学A和B把同样的数据结果给老师看，A单纯把数据列到表格里面，B做成了上图的形式，老师会喜欢谁？此外，我们也可以通过写一个简单的python脚本来实现作图的功能。

2.6 作图的 Python 脚本：

2.6.1 获取数据，并保存成文件

注意echo $i 后面的逗号，输出的文件 data.dat 中，逗号用于将两列数据分开。

2.6.1.1）图中圈出来的部分中， > 代表将前面命令的输出保存到 data.dat 文件中；

2.6.1.2） data.dat 后缀可以随便写，data.txt， data.out, data.export, 也可以不写：data 。 因为输出后的都是文本格式，直接可以编辑打开。

2.6.1.3） 查看一下保存的数据结果

2.6.2 作图 脚本， 名为 plt.py

vim打开后，内容如下图：

python程序学习参考书： learn python the hard way（免费），网上还有很多学习资料，这里就不介绍了。python读取data.dat 文件，根据 delimiter 后面的参数（此处为逗号）将数据分成若干列。

2.6.3 运行脚本：（python plt.py）

2.7 体系的能量随ENCUT的变化：

获取图上结果的命令行：

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

作图过程如下：

2.7.1 生成数据：

前面的命令中， “\t” 改为 逗号， done 后面 加上 > data.dat， 重新运行一遍。

2.7.2 修改脚本如下：

增加了 一列（z）

2.7.3 选择合适的画图区间：

图中蓝色的O$_2$分子的能量随ENCUT值增加的变化，为一直线！

大师兄问：结果对不对？

答：结果是对的，但不具有任何的合理性。

因为它不能体现出能量变化的趋势：因为 y 轴的取值范围太大了，而我们能量的变化又太小。这也是很多人在做模拟的时候容易犯的错误，看到一个直线，或者平滑的曲线，就如同哥伦布发现了新大陆一番。其实则不然，很多时候，真相往往被我们的粗糙的观测范围给掩盖了。

这个称之为作图骗人的一个Trick！很多人在发文章的时候将这一技巧运用的如火纯情。所以，当你看文献的时候，如果发现别人的结果也是一条直线或者平滑曲线，第一直觉是去看坐标轴的范围，而不是感觉别人的工作是多么地牛逼。

重新作图，这次只用ENCUT和O$_2$分子的能量，如下图：

能量一直在降低， 结果对不对呢？对的，但变化却被过度放大了！

很多新手看到这个图以为能量降低的很快，感觉自己的计算不收敛，因此还要继续往下做或者测试。其实则不然，在这里，我们还是需要主要 y 轴的取值范围。图中整个y轴的变化为0.016 eV。每两个点之间的变化仅仅为0.002 eV。这个变化我们也是可以接受的，虽然图中感觉下降的很厉害的样子。

上面两个图中，走了两个极端，一个极大，一个极小。所以，在我们处理数据时，一定要仔细观察坐标轴的取值范围。

2.8 小结一下：

但从数据上可以看出来，O$_2$分子的能量随着ENCUT值的增加还是会有些少许的波动。比如，从400到800 eV，能量降低了0.01 eV。大师兄，那我们该取哪个ENCUT值呢？

首先，我们要明确几个要点：

2.8.1）ENCUT值越大，计算的越精确，花费的时间也就越多；

2.8.2）ENCUT的取值仅仅通过一个例子的测试来确定，这是不对的；

  • A) 在VASP的计算中，单个结构的能量所具有的意义不大，也就是说相对能量最重要！
  • B) 为什么A）中说相对能量最重要：这是因为：没有对比，就没有伤害！设想一下，你把一个O$_2$的能量算的再精确或者能量再低，而不去使用它，它也就是个数字而已，不具有任何的物理化学意义。
  • C) B)中要表达的意思是，算出来的数值必须要应用到我们的物理化学概念中，也就是物理化学的概念体现在这些能量的使用过程中，也就是相对能量里面：比如，O$_2$的结合能， O$_2$的吸附能；某一化学反应的能量，反应能垒，表面能，功函数等等，无一不是多个能量的数学运算所的出来的。
  • D) O$_2$的结合能示例：

公式： EB = E(O$_2$) –2E(O)

O$_2$的结合能随着ENCUT值增加的变化情况：

第一、二、三、四列分别为ENCUT，O原子能量，O$_2$分子能量和O$_2$分子中O的结合能。从图中，我们可以分析出来，ENCUT在450 eV时，O$_2$的结合能和ENCUT为800 eV的时候差别很小，因此我们可以选取450 eV进行计算。

注： O原子不同ENCUT的单点能，自己根据前面的过程，补充计算。

2.8.3）ENCUT的取值与体系中所有的元素有关；

  • A）体系中含有不同元素的时候，查看这些元素POTCAR中的ENMAX值，找出最大的那个；
  • B）ENCUT的最小值为所有元素中ENMAX的最大值！
  • C）也就是说，找到最大的ENMAX（max）， ENCUT值大于等于ENMAX（max）。

示例图中，
potcar.sh 为生成POTCAR的脚本，可在本书的附录中获取。ENCUT 的取值至少为 400 eV。

2.8.4） ENCUT的值，也要查阅相关的参考文献进行确定。你的计算体系大家都用400 eV，你也可以设置该值或者稍微高一些。ENCUT很多人在刚开始做计算的时候，都会测试一个数值。如果不想测试，查找参考文献其实是一个很好的办法。

3 扩展练习：

3.1 复习前面学到的批量操作方式；

3.2 熟练运用前面的操作，并理解命令的工作原理；

3.3 改变KPOINTS的大小，查看计算时间，能量的变化；

4 总结：

1 EDIFF和EDIFFG对计算的影响通过实例强调了一次；

2 ENCUT 测试的操作流程；

3 ENCUT对计算时间和能量的影响结果分析；

4 通过excel 和python作图的两个方法；

5 能量和计算时间随ENCUT的变化；

6 看图时坐标轴的区间范围要注意；

7 总结ENCUT取值的一些注意事项。

LVASPTHW

Ex16 谁偷走的我的机时？（一）

前面O$_2$初始结构为0.9 $\AA$时，如果使用BRION=2，相对于初始结构为1.07$\AA$的时候，为避免过度矫正，我们需要设置一个更小的POTIM，大师兄尝试过了，POTIM = 0.01 的时候得到了正确的计算结果。到此，O$_2$分子的计算我们暂且告一段落。做计算，机时很重要，尤其是缺钱，捉襟见肘的时候，我们就更加需要珍惜，保证并提高自己的计算成功率，避免重复计算。这一节我们总结下前面讲到的影响计算时间的一些细节。大家在计算的时候知道如何把握时间，在有限的机时内获取更多的有效的计算结果。

1 体系的磁性 (EX8)

考虑自旋后（ISPIN=2），VASP计算时会将电子分为两部分处理，一部分是$\alpha$电子，另一部分是$\beta$电子，我们在氧原子的自旋极化计算中提到过了，忘了的可以查看Ex8，Ex11中的内容。

补充一下：如果原子的MAGMOM为负值：应该这么写：

1
MAGMOM = 10*-2 # 有10个原子，每个原子的初始磁矩为 -2

注意： -2 不用括号（）括起来

2 对称性 (EX8)

体系的对称性降低，会增加相应的计算量。查看O原子计算时，改变晶胞大小，取消对称性前后的计算时间。

3 SIGMA的取值 (EX1)

我们看一下VASP练习手册：handonsession-I 里面的一句话说明：

因此，设置一个较小的SIGMA值会使收敛变慢。大家对于ISMEAR的选择，一定要多看官网的说明！

3.1 半导体和绝缘体：

1） K点小于4 的时候，用ISMEAR=0, SIGMA取值小一些；比如SIGMA=0.05; 此时用ISMEAR = -5 会出错；

2）K点数目大于4的时候，可以使用ISMEAR= -5；

3） 注意： 我们算的气体分子，是绝缘体，且只用了gamma点，所以我们参数的设置如下： 

1
2
ISMEAR=0 
SIGMA=0.01

3.2 金属体系：

ISMEAR一般用ISMEAR=0 或者整数1，2即可。

SIGMA =0.1 足够。

4 合理的初始结构 （Ex13-15）

前面我们刚讲到，一个合理的初始结构，可以避免很多意外的错误以及快速得到正确的结果。

如果你的初始结构不合理导致的计算出错，首先应该去调整结构，最后才是去调节参数，比如上一节中IBRION=2时的POTIM值。

5 EDIFF （Ex9）

EDIFF 控制了电子迭代的收敛标准，如果你设置的标准比较严格，则每一个离子步需要更多的电子步数，需要的时间也会随着离子步的增加成线性关系增长。

6 EDIFFG （Ex9）

EDIFFG控制了结构收敛的标准，同样严格的标准需要更多的结构优化步骤来实现。

扩展练习：

1 本文所提到的内容，必须去官网查看，并认真阅读，如有不懂的，请在群里自由提问；

2 思考一下，还有其他因素影响我们的机时吗？

LVASPTHW

Ex15 论合理结构的重要性（三）

继续前面一节，我们分析一下结构优化过程的细节问题，以及合理结构的重要性。

哪里出错了？

查看每一离子步后的能量结果，注意命令中的单引号里面有两个空格！如果你看到这样的结果时，意味着大事不妙，结构优化失败了，SCF 也失败了！！！也就是结构没有优化好，离子步中电子步收敛也同样失败。

查看一下OSZICAR:

1
2
3
4
5
6
7
8
9
10
DAV:  52     0.142641482694E+03   -0.87447E-01   -0.50432E-03   120   0.175E-01    0.188E+00
DAV:  53     0.142675598279E+03    0.34116E-01   -0.20261E-03   120   0.869E-02    0.185E+00
DAV:  54     0.142703773182E+03    0.28175E-01   -0.26393E-03   144   0.121E-01    0.190E+00
DAV:  55     0.142724228087E+03    0.20455E-01   -0.11280E-03   120   0.795E-02    0.193E+00
DAV:  56     0.142809666479E+03    0.85438E-01   -0.29739E-02   120   0.433E-01    0.211E+00
DAV:  57     0.142939449864E+03    0.12978E+00   -0.29353E-02   192   0.367E-01    0.196E+00
DAV:  58     0.142943989012E+03    0.45391E-02   -0.90264E-04    96   0.810E-02    0.193E+00
DAV:  59     0.142966025527E+03    0.22037E-01   -0.43114E-04   120   0.704E-02    0.196E+00
DAV:  60     0.142973535062E+03    0.75095E-02   -0.53590E-05   144   0.234E-02
   3 F= 0.14297354E+03 E0= 0.14297354E+03  d E =0.142711E+03  mag=     0.0441

这里说明，第三步中电子收敛的步数达到了默认值，被强制停止了。每一个离子步中电子步数的最大值是由NELM这个参数控制的，VASP中默认值为60，也就是最大的电子收敛为60步，过了60步还不收敛就必须停止。自己VASP官网查看下相关的NELM参数。

下面两行如果你能看懂，说明优化的基本情况已经摸清了。

电子步（SCF）: EDIFF <====> NELM

离子步（结构优化）：EDIFFG <====> NSW

电子步不收敛怎么办？

遇到图中这种情况，大师兄分享一下自己的经验：

1） 首先检查自己的初始结构是不是合理的。这是关键点之一。如果合理，那么再进行下面的步骤：

2）如果第一个离子步中：SCF（也就是电子步）的计算不收敛，尝试下增加NELM的值；

  • 对于一般普通的体系可以设置NELM = 100
  • 对于一些电子结构比较难收敛的体系，可以设置更大一些：NELM = 200
  • 增加NELM后依然不收敛，尝试下改变AMIX，BMIX，官网推荐的参数如下：不过个人的感觉，调这些参数好像没什么用，效果甚微。
1
2
3
4
AMIX = 0.2
BMIX = 0.00001
AMIX_MAG = 0.8
BMIX_MAG = 0.00001

3) 第一个离子步中的电子步收敛了，后面的不收敛，能量变的极大（本例），首先应该想到的是去检查结构，一般在结构不合理的时候会出现类似的情况；调整结构再提交任务。

4）如果前面几个离子步中电子步都不收敛，且能量变化正常，可以尝试着让任务再多算几步，后面跑着跑着可能就收敛了。

5）如果跑了很长时间，每一个离子步中的电子步都不收敛，可以尝试着换一个更加稳定的电子步迭代算法（ALGO参数）。这种情况：ALGO = ALL 结合 NELM = 200 可以解决大部分的问题。

6）但需要注意的是：ALGO= ALL 这个算法虽然稳定，但比较耗时。自己要权衡一下。可以尝试着两步走的战略：step1）先用ALGO =ALL的办法算一个单点，保存WAVECAR。step2）然后将ALGO改回原来的，读WAVECAR继续优化。

大师兄尝试了很多种办法，意图把初始结构从0.9 $\AA$一步直接计算正确，但是大部分时候都失败了！而且出现了各种各样的问题和错误的结果，因此，出现这种情况，第一直觉是去看结构而不是想着调节参数去怎么解决这个错误！！！从这里可以看到，如果你的初始结构不合理，会出现各种各样的问题，这也是大师兄们解决不了的。

为什么会这样子

这是因为：两个氧原子之间距离在初始结构中很小，导致第一步估算之间的作用力过大，以至于后面没有办法再矫正过来。这里我们看一下IBRION=2时的计算步骤。

第一步，从初始结构出发，计算体系中离子间的作用力，

第二步，VASP尝试着把离子沿着前面估算的方向移动，尝试移动的大小由POTIM这一项决定，

第三步，计算尝试移动后能量和力的大小，据此加入一个矫正项来控制真实移动的大小；

第四部， 移动后，重新计算能量和力，重复前三步直至能量或者力收敛到我们设置的EDIFFG值。

IBRION = 2 时，对POTIM的依赖性很强，因此我们计算的时候要设置一个合理值。在我们的计算中，由于初始的原子间距离很小，第一步计算时，得到的原子间的初始排斥力很强，第二步中，VASP默认的POTIM值是0.50，前面两步导致了尝试步中离子的移动过大，以至于后面没有办法矫正回来，最后导致O$_2$分子计算出错。

POTIM 显神威

大师兄，用0.9 $\AA$的时候可以调节POTIM来获得正确的结果吗？答案是肯定的。

如果想要正确计算的话，可以设置POTIM一个更小的值。POTIM=0.1，虽然从初始值算出来的力很大，我们通过POTIM强制VASP一点一点调节，来保证计算的准确。

1
2
3
4
5
6
7
8
SYSTEM = O atom
ISMEAR = 0
SIGMA = 0.01
ISPIN = 2
MAGMOM = 2*2
IBRION = 2
POTIM = 0.1
NSW = 10

提交任务，查看结果：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex15$ tail OSZICAR
DAV:   3    -0.984191032147E+01   -0.12318E-01   -0.38005E-02    96   0.891E-01    0.516E-01
DAV:   4    -0.984475378528E+01   -0.28435E-02   -0.74281E-03   144   0.367E-01    0.198E-01
DAV:   5    -0.984632399820E+01   -0.15702E-02   -0.51626E-04    96   0.101E-01    0.914E-02
DAV:   6    -0.984676440618E+01   -0.44041E-03   -0.15294E-04   144   0.511E-02    0.169E-02
DAV:   7    -0.984690803808E+01   -0.14363E-03   -0.99478E-06   144   0.122E-02    0.829E-03
DAV:   8    -0.984706391665E+01   -0.15588E-03   -0.95528E-06   120   0.898E-03    0.489E-03
DAV:   9    -0.984722304729E+01   -0.15913E-03   -0.99988E-06   120   0.862E-03    0.292E-03
DAV:  10    -0.984733039102E+01   -0.10734E-03   -0.63466E-06    96   0.690E-03    0.215E-03
DAV:  11    -0.984738716468E+01   -0.56774E-04   -0.27665E-06    96   0.463E-03
   7 F= -.98473872E+01 E0= -.98473872E+01  d E =-.101099E+02  mag=    -2.0000
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex15$ cat CONTCAR
O
   1.00000000000000
     7.5000000000000000    0.0000000000000000    0.0000000000000000
     0.0000000000000000    8.0000000000000000    0.0000000000000000
     0.0000000000000000    0.0000000000000000    8.9000000000000004
   O
     2
Direct
  0.0000000000000000  0.0000000000000000 -0.0196261760084087
  0.0000000000000000  0.0000000000000000  0.1207497715140281

  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex15$ python
Python 2.6.6 (r266:84292， Sep  4 2013， 07:46:00)
[GCC 4.4.7 20120313 (Red Hat 4.4.7-3)] on linux2
Type "help"， "copyright"， "credits" or "license" for more information.
>>>(0.1207497715140281--0.0196261760084087)*8.9
1.2493459329496877
>>>

算出来的能量和磁矩等信息与之前正确计算的结果一致，键长也对上了。 说明我们修改POTIM生效了。这里POTIM的作用相信大家有了一个大体的了解。在IBRION=2 时（这是前提！！），如果初始结构很差，设置较小的POTIM可以有效的避免过度矫正。初始结构越好，POTIM的选择也就越随意。

Python计算器

在上面的演示中，大师兄还教给你了一个终端里面的计算器，python！加减乘除非常好用，注意，使用python计算时，>>> 和数字之间不能有空格，否则python会提示出错！

1
2
3
4
5
6
7
8
9
10
11
12
13
Python 2.6.6 (r266:84292， Sep  4 2013， 07:46:00)
[GCC 4.4.7 20120313 (Red Hat 4.4.7-3)] on linux2
Type "help"， "copyright"， "credits" or "license" for more information.
>>>  5 + 2
  File "<stdin>"， line 1
    5 + 2
    ^
IndentationError: unexpected indent
>>>  (0.1207497715140281--0.0196261760084087)*8.9
  File "<stdin>"， line 1
    (0.1207497715140281--0.0196261760084087)*8.9
    ^
IndentationError: unexpected indent

扩展练习

1 阅读IBRION 参数说明：https://cms.mpi.univie.ac.at/vasp/vasp/IBRION_2.html 知道优化过程的具体细节；

2 尝试不同的初始结构配合POTIM值，完成O$_2$的正确计算；

3 使用0.9 $\AA$作为初始结构，改变POTIM的参数值，直至计算正确；

4 尝试用0.9 $\AA$作为初始结构，POTIM采用默认值，调节其他自己认为会影响计算的参数，查看相关输出，如有错误，复制到google里面查找相关的原因。

总结：

体系的初始结构越合理，不仅仅是节约我们的时间，还会我们的计算过程越省心。从前面出现的问题可以看出来，结构不合理，对应出错的地方越多，我们也要绞尽脑汁去纠正。O$_2$分子的例子比较简单，我们简单通过POTIM实现了正确的计算。但是大家计算的体系比O$_2$ 分子复杂的多，处理起来也会更加棘手。所以需要在初始结构上下功夫。更严肃的说：初始结构就是我们的计算模型，如果模型不合理，后续的计算都会出错，甚至失败。

如果初始结构不合理，用IBRION=2的时候，POTIM可以很好的控制收敛，还是建议初始结构搭建的合理些，省时省力。认真学习化学基础知识，搭建好合理的初始模型才是王道，正应了那句古话：磨刀不误砍柴工！

LVASPTHW

Ex14 论合理结构的重要性（二）

在ex13中，当初始值为1.5 $\AA$ 的时候，计算共进行了9步。对比下之前我们采用实验值（Ex11）作为初始结构计算时用了4步。从这里我们可以看出来，如果你一个合理的初始结构，可以加快优化的速度，减少机时，节约你的时间。当然，计算的具体时间可以通过OUTCAR尾部的信息查看。

不同初始结构对结果的影响。

那么不同的初始结构，除了在时间上，对计算的结构还有影响吗？ 首先我们对比下Ex11和Ex13的结果。 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW$ tail ex13/nsw10/CONTCAR
     0.0000000000000000    8.0000000000000000    0.0000000000000000
     0.0000000000000000    0.0000000000000000    8.9000000000000004
   O
     2
Direct
  0.0000000000000000  0.0000000000000000  0.0149145061380336
  0.0000000000000000  0.0000000000000000  0.1536248197046605

  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW$ tail ex11/opt/CONTCAR
     0.0000000000000000    8.0000000000000000    0.0000000000000000
     0.0000000000000000    0.0000000000000000    8.9000000000000004
   O
     2
Direct
  0.0000000000000000  0.0000000000000000 -0.0019537557431563
  0.0000000000000000  0.0000000000000000  0.1367852164173130

  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW$

在Ex11和Ex13中，两个计算的能量和CONTCAR中的键长值也几乎相等，看下面的结果。

1
2
3
4
5
6
7
8
9
10
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW$ python
Python 2.6.6 (r266:84292, Sep  4 2013, 07:46:00)
[GCC 4.4.7 20120313 (Red Hat 4.4.7-3)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> 0.1536248197046605-0.0149145061380336
0.13871031356662689
>>> 0.1367852164173130--0.0019537557431563
0.13873897216046929
>>>

另一个bad结构的计算

上次我们用了一个大于实验值的初始结构1.5$\AA$。下面我们看下小于实验值的初始情况：0.9 $\AA$ 。

1
2
3
4
5
6
7
8
9
10
11
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex14$ cat POSCAR
O
1.0
7.5 0.0 0.0
0.0 8.0 0.0
0.0 0.0 8.9
O
2
Cartesian
0.0 0.0 0.0
0.0 0.0 0.9

INCAR，KPOINTS，POTCAR等不变，提交任务，等待计算结束，查看结果。

1
2
3
4
5
6
7
8
9
10
11
12
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex14$ tail OSZICAR
DAV:  18     0.141396114006E+03    0.20879E-02   -0.88760E-04   144   0.130E-01    0.579E-01
DAV:  19     0.141397297948E+03    0.11839E-02   -0.52751E-04    96   0.893E-02    0.455E-01
DAV:  20     0.141388687914E+03   -0.86100E-02   -0.28657E-03   144   0.231E-01    0.206E+00
DAV:  21     0.141395321236E+03    0.66333E-02   -0.62073E-03    96   0.361E-01    0.112E+00
DAV:  22     0.141389669824E+03   -0.56514E-02   -0.95109E-03    96   0.429E-01    0.161E-01
DAV:  23     0.141382943017E+03   -0.67268E-02   -0.22808E-03    96   0.216E-01    0.118E+00
DAV:  24     0.141390086243E+03    0.71432E-02   -0.39635E-03    96   0.307E-01    0.825E-01
DAV:  25     0.141385767811E+03   -0.43184E-02   -0.44984E-03    96   0.305E-01    0.480E-01
DAV:  26     0.141385726226E+03   -0.41585E-04   -0.10329E-04    96   0.429E-02
   5 F= 0.14138573E+03 E0= 0.14138680E+03  d E =0.141123E+03  mag=    -2.0000

计算共进行了5步，且最后的磁矩看起来是正确的。思考： 从这里得出的信息，能确定我们的计算结果是正确的吗？

答： 不知道。因为我们还要要去查看一下结构。判断结构是否合理。如果结构不合理，则收敛的计算也是失败了。

用p4vasp打开CONTCAR后，如下图：

师兄，这是神马情况，两个原子怎么跑这么远？

不用担心，这是因为周期性的原因。p4vasp中可以进行如下的操作：

首先，点击左侧的Build按钮，然后再点击右侧的 To unit cell。这样你会发现结构调整到下图的样子：

两个原子之间的距离还是很长（7.821 $\AA$），但实际键长不是这么长的。

而是8.9-7.821 = 1.079 $\AA$。

师兄你为什么这么算?

因为我们的体系是周期性的，也就是图中的格子在三维方向上可以无限重复，如果我们向左重复一个单元，那么在新的单元中右侧的氧原子与原来左侧的氧原子距离很短。已知格子在z方向的长度为8.9 $\AA$，减去7.821就是剩下的两个氧原子之间的键长了。

如果，你还不明白，进行下图的操作：

点击左侧的Control 选项，然后在下面红色框中，将格子在三维方向上重复，效果如下：

注意，该操作只是展示三维方向的结构，如果此时你保存结构，不管你在三维方向上重复了多少次，保存的结构则还是原来的尺寸大小。

键长为1.0785 $\AA$。

周期性的显示问题

暂且抛开对错不说，由于周期性导致的原子不在一个格子里面的情况，在今后的计算中你会经常碰到。如果你遇到这种情况，不要立即在群里问：师兄，为什么优化之后，体系中的原子不见了？为什么之前左面原子不见了，右侧本来没有原子，优化完多了？

归根结底都是周期性导致的显示问题。你需要做的就是把结构在三维方向上重复一下，查看结构是对还是错。

结果对还是错？

现在我们分析下对错，已知O$_2$分子的键长为1.2075 $\AA$，因此该计算与实验值偏差为：（1.0785-1.2075）/1.2075 =10.68 %，这么大的偏差，是不可以忍受的。

检查一下能量： 为 -8.54642426 eV。 之前正确的能量为： -9.85498627eV。

在第二版的改进中，大师兄又计算了这个任务一次，得到了另外一个结果：2个O原子距离很短。

查看能量：

1
2
3
4
5
6
7
8
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex14$ grep '  without' OUTCAR
  energy  without entropy=        0.26254047  energy(sigma->0) =        0.26254047
  energy  without entropy=       11.76187113  energy(sigma->0) =       11.76187113
  energy  without entropy=      142.97353506  energy(sigma->0) =      142.97353506
  energy  without entropy=  1359169.21650280  energy(sigma->0) =  1359169.21650280
  energy  without entropy=      141.38788142  energy(sigma->0) =      141.38680382
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex14$

思考： 能量为什么会差这么多呢（1.11 eV 或者最新结果140多eV）？

答： 我们需要知道体系的能量随键长的变化关系：如下图：

图中X 处是O$_2$的稳定结构，两个原子间距离小于X处的键长时，它们之间的排斥力导致了体系的能量快速升高。

由于我们已经知道了正确的计算结果，通过分析后，这次的计算失败！但对于不知道结果的时候，怎么判断计算是否成功失败呢？首先根据VASP计算的收敛情况，也就是计算至少应该正常结束，其次，这是远远不够的，我们还要查看输出结构的几何构型，判断是否具有物理或者化学的意义，还要看每一步收敛的能量信息。这些就需要我们化学基础知识了。

扩展练习及思考

1 计算为什么会失败?

2 分析该计算中每一步收敛的情况，以及能量的变化。

总结：

1 不合理的结构会增加计算时间；

2 不合理的结构会导致计算结果没意义；

3 知道怎么处理周期性结构中，原子不在一个晶格里面的情况；

4 学会判断计算结果的物理或者化学意义。

LVASPTHW

Ex13 论合理结构的重要性（一）

前面我们学会了O$_2$分子的优化，分析了其电子构型，并且知道了理论结果和实验结果间存在的偏差。这一节我们依然研究O$_2$分子的结构优化，虽然模型简单，但熟练掌握其中的技巧，对今后的计算工作，意义重大。

测试Bad 结构

首先，我们将O$_2$的初始键长设置为：1.5 Å。其他文件与Ex11中的保持不变。 

1
2
3
4
5
6
7
8
9
10
O
1.0
7.5 0.0 0.0
0.0 8.0 0.0
0.0 0.0 8.9
O
2
Cartesian
0.0 0.0 0.0
0.0 0.0 1.5

然后提交命令进行计算，等待结束后。查看结果，并对比前面Ex11练习的结果。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13$ tail OSZICAR
DAV:   1    -0.985412177027E+01   -0.57432E-02   -0.12528E+00   144   0.457E+00    0.433E-01
DAV:   2    -0.985315825637E+01    0.96351E-03   -0.76431E-03   144   0.367E-01    0.166E-01
DAV:   3    -0.985344344928E+01   -0.28519E-03   -0.78039E-04   120   0.129E-01    0.707E-02
DAV:   4    -0.985368960614E+01   -0.24616E-03   -0.14892E-04   144   0.510E-02    0.268E-02
DAV:   5    -0.985378634076E+01   -0.96735E-04   -0.12727E-05    96   0.167E-02
   8 F= -.98537863E+01 E0= -.98537863E+01  d E =-.142913E+01  mag=     2.0000
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.985528191871E+01   -0.15923E-02   -0.11782E-01   168   0.140E+00    0.130E-01
DAV:   2    -0.985521232023E+01    0.69598E-04   -0.82384E-04   120   0.126E-01
   9 F= -.98552123E+01 E0= -.98552123E+01  d E =-.143055E+01  mag=     2.0000
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13$ tail ../ex11/opt/OSZICAR
DAV:   8    -0.985382407600E+01   -0.31391E-03   -0.20867E-05   120   0.151E-02    0.891E-03
DAV:   9    -0.985418412704E+01   -0.36005E-03   -0.33607E-05   120   0.162E-02    0.520E-03
DAV:  10    -0.985440719048E+01   -0.22306E-03   -0.20675E-05   120   0.115E-02    0.367E-03
DAV:  11    -0.985450935523E+01   -0.10216E-03   -0.85511E-06   120   0.750E-03    0.174E-03
DAV:  12    -0.985454505708E+01   -0.35702E-04   -0.14739E-06   120   0.339E-03
   3 F= -.98545451E+01 E0= -.98545451E+01  d E =-.451963E-01  mag=     2.0000
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.985517058251E+01   -0.66123E-03   -0.43217E-02    96   0.851E-01    0.792E-02
DAV:   2    -0.985515527518E+01    0.15307E-04   -0.49070E-04   168   0.969E-02
   4 F= -.98551553E+01 E0= -.98551553E+01  d E =-.458065E-01  mag=     2.0000

在本练习中，我们设置了键长为：1.5 $\AA$，距离稳定结构较大，用了9个离子步才收敛。在Ex11练习中，我们使用的是数据库中O$_2$的键长，只优化了4步就收敛了。由此可见，一个好的初始结构可以加快我们的计算。那么有多快呢？ 我们看一下OUTCAR中的时间信息，对比看下一目了然。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13$ tail OUTCAR
                            User time (sec):       20.196
                          System time (sec):        7.877
                         Elapsed time (sec):      144.523

                   Maximum memory used (kb):      116704.
                   Average memory used (kb):           0.

                          Minor page faults:        72067
                          Major page faults:            0
                 Voluntary context switches:         3191
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13$ tail ../ex11/opt/OUTCAR
                            User time (sec):        8.387
                          System time (sec):        7.789
                         Elapsed time (sec):       17.433

                   Maximum memory used (kb):      114276.
                   Average memory used (kb):           0.

                          Minor page faults:        26287
                          Major page faults:            3
                 Voluntary context switches:         3202
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13$

判断计算是否结束以及结构是否收敛：上面大师兄用了2个命令：

1） tail OSZICAR

左下角显示了结构优化进行了多少步，这里是9步！注意，我们的INCAR 里面设置的NSW = 10。如果优化的步数小于我们设置的，说明结构已经收敛到了我们所期望的标准。

2） tail OUTCAR

会显示：VASP结束后计算的内存，时间等信息，如果你看到与上图中类似的信息，说明计算结束了。

优化的小思考

思考1：

上面的计算中，我们在INCAR中设置的NSW =10，实际计算用了9步，说明计算刚刚收敛了。但是：如果结构优化的步数等于NSW设置的步数，说明了什么，是什么原因？又该怎么办？

出现这种情况，说明可能还没有达到我们的收敛标准，可能是以下几个原因造成的：

1) NSW设置的偏小； 如果我们在本节例子中设置了NSW=8;

2) 初始结构不合理，计算需要更多的离子驰豫过程；

3) 设置的收敛标准太严格， 比如：-0.01 或者 -0.001；

4) 结构很复杂，每一离子步中的电子步骤收敛很困难。

除此之外，还有一种可能，即刚刚进行到NSW设置的步数时，计算恰好收敛了。比如本例中，我们设置了NSW = 9。但这种可能几率很低，但如果你计算的足够多，还是有机会碰到的。（碰到这种情况时，一定不要慌张，要立马出门去买彩票，中了奖后要记得和大师兄平分！）

思考2：

怎么判断上述这种特殊情况呢？

首先，我们要知道计算收敛结束后VASP所输出的内容，也要知道未收敛时VASP结束后输出的内容。目前，现在我们知道进行了9步的时候，该计算正常结束。那我们就可以设置NSW为一个较小的值（比如 NSW=5），然后查看下未收敛时候的结果，并进行对比分析！

结构收敛OUTCAR （NSW=10）：

结构未收敛OUTCAR （NSW=2）：

通过对比，我们可以发现，结构收敛的结果里面多了这一行：　 

1
reached required accuracy - stopping structural energy minimisation

因此我们在今后的计算中，可以通过这一行来判断计算是否收敛结束。怎么判断呢？前面我们讲到了通过分析OUTCAR并结合grep命令来提取信息的方法，这里我们稍微复习一遍。下图中大师兄尝试了几个grep命令：

看下输出结果:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13/nsw10$ grep accuracy   OUTCAR
 reached required accuracy - stopping structural energy minimisation
 iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13/nsw10$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13/nsw10$ grep stopping OUTCAR
   EDIFF  = 0.1E-03   stopping-criterion for ELM
   EDIFFG = 0.1E-02   stopping-criterion for IOM
 reached required accuracy - stopping structural energy minimisation
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13/nsw10$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13/nsw10$ grep structural OUTCAR
Analysis of structural, dynamic, and magnetic symmetry:
 reached required accuracy - stopping structural energy minimisation
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13/nsw10$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13/nsw10$ grep 'reached required accuracy' OUTCAR
 reached required accuracy - stopping structural energy minimisation
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13/nsw10$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13/nsw10$ grep 'reached required accuracy - stopping structural energy minimisation' OUTCAR
 reached required accuracy - stopping structural energy minimisation
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13/nsw10$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13/nsw10$ grep reached OUTCAR
------------------------ aborting loop because EDIFF is reached ----------------------------------------
------------------------ aborting loop because EDIFF is reached ----------------------------------------
------------------------ aborting loop because EDIFF is reached ----------------------------------------
------------------------ aborting loop because EDIFF is reached ----------------------------------------
------------------------ aborting loop because EDIFF is reached ----------------------------------------
------------------------ aborting loop because EDIFF is reached ----------------------------------------
------------------------ aborting loop because EDIFF is reached ----------------------------------------
------------------------ aborting loop because EDIFF is reached ----------------------------------------
------------------------ aborting loop because EDIFF is reached ----------------------------------------
 reached required accuracy - stopping structural energy minimisation

从图中可以看出来，选择不同的grep 参数，会得到不同的输出结果。一般来说，grep后面的内容越详细，得到的输出结果也就会越匹配。我们提取信息的首要原则是：精确匹配并提取最有价值的信息！

其中：grep ‘reached required accuracy’ OUTCAR 是大师兄常用的一个命令。

（要提取的内容通过单引号扩起来了。）

~/.bashrc

那么你会问道：“大师兄，你是不是傻X啊，这么长的命令，光输入就费老大劲了，不嫌麻烦么？”

哈哈，大师兄不仅傻，更喜欢偷懒。每次检查的时候这样子输入肯定麻烦。但在linux下面，有个小窍门，可以极大极大极大地提高你的工作效率，那就是~/.bashrc 文件中的 alias：（注意前面的点 . ）

.bashrc 文件怎么用呢？它在哪里？ 怎么打开？

答：在home目录下，

1） 终端（terminal）里面直接输入 cd 这个命令，会自动返回到home 目录下面：

然后运行命令：vi .bashrc

2） 当然，如果你不想跳转回去，也可以这样：vi ~/.bashrc 波浪号代表home

它怎么用?

http://www.linuxidc.com/Linux/2015-02/113310.htm （参考网址）

1）打开 .bashrc 文件：（大师兄的bashrc文件为例：）

2） 将前面我们用到的命令写进 .bashrc 文件：

注意的部分：

A) alias 和 gr 之间有空格；

B) gr 是大师兄随便想的，你也可以用自己想的其他字母；

C) gr之间没有空格；

D) gr后面紧跟着等号 = ， 中间没有空格；

E) 等号 = 后面紧跟着双引号，且等号= 和双引号之间没有空格；

F) 双引号之间，把我们的命令放进去，一定要确保命令在引号里面；

G)等号 = 后面可以用单引号，也可以用双引号，因为我们的命令中已经有单引号了，这里我们用的双引号，下面两者效果是一样的；

1
2
alias gr="grep 'reached required accuracy' OUTCAR"
alias gr='grep "reached required accuracy" OUTCAR'

注意：上图中两者任选一个即可，不要都放进去，否则会乱套。保存后退出。

3) source 一下 .bashrc 文件。下图中的几个命令效果是一样的！

在这里点 . 和 source 命令的效果一样。

4 运行命令：（敲一下gr，回车即可）

1
2
3
4
5
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13/nsw10$ ls
CHG  CHGCAR  CONTCAR  DOSCAR  EIGENVAL  IBZKPT  INCAR  job_sub  KPOINTS  OSZICAR  OUTCAR  PCDAT  POSCAR  POTCAR  REPORT  slurm-1023075.out  vasprun.xml  WAVECAR  XDATCAR
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13/nsw10$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13/nsw10$ gr
 reached required accuracy - stopping structural energy minimisation

当你看到这个信息的时候，说明优化任务就已经算完了。 在这里gr 命令就是前面 .bashrc 文件中alias 后面的那个命令。如果还不明白，看下面的例子。大师兄把.bashrc 文件中的 gr 替换成了bigbro , source 了一下后，运行 bigbro命令，得到了和前面一样的结果。

1
alias bigbro='grep "reached required accuracy" OUTCAR'
1
2
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex13/nsw10$ bigbro
 reached required accuracy - stopping structural energy minimisation

5 ）注意：

避免你自己定义的命令名字和linux下面自带的命令重复。

1
alias  cd="grep 'reached required accuracy' OUTCAR"

但是可以这样设置：

1
alias  cd="cd && ls "

比如：在linux系统下，使用cp命令复制文件夹的时候，需要用到 cp -r A B 。有时候会忘记加 -r而出错，我们可以这样设置，alias cp='cp -r '来避免出错。

1
2
3
4
5
6
7
8
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW$ cp ex13 ex14
cp: omitting directory `ex13'
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW$ cp ex13 ex14 -r
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW$ vi ~/.bashrc
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW$ .  ~/.bashrc
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW$ cp ex13 ex15
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW$

如果你极其地讨厌一个人，可以在他的服务器下面这样设置！！！ 建议大家报复的时候可以采用这样温柔点的设置，这样，他每敲一次命令就会被你骂一次。（史上最狠报复，只适用于对付那些深深伤害了你的人！）大家只可以感受下这把刀的威力，切不可自作孽!

1
2
3
alias cd='cd  ~'
alias cp='echo "You are a XXX" '
alias ls='echo "You are a XXX" '

扩展练习：

1 分别设置O$_2$的初始键长为1.5和 0.9 Å，运行任务；

2 查看收敛情况；

3 查看结构信息和能量信息；

4 如果有不合理的地方，自己尝试解决。

总结：

本节讲解了

1) 查看任务计算完成的日常做法，

2) 一个在日常计算中极为重要的处理问题的方法，以及简化日常命令的方式。

通过设置.bashrc文件中的alias 命令,我们便可以自己随心所欲地将计算中的长命令、复杂的命令转化为简单的命令,从而极大地提高了我们的工作效率,尤其是适用于那些:使用频繁的但又很复杂的命令。

3) 一个思想：如果我们已知一个结果，那么便可以通过调控不同的参数，然后再与已知结果进行对比，从而分析不同参数在计算中的作用。

这一个思想贯穿了本书的写作过程，对于大家加深对VASP的理解以及练习意义非凡。熟练运用这一种方法，你的洞察力和判断力会得到极大地提升。希望大家能认真思考下本节中提出的处理问题的思想，以及熟练运行alias简化自己的命令。

4） 为了锻炼大家的判断力，请做下题：以最快的速度找到图中的小猫：

LVASPTHW

Ex12 VASP中MAGMOM的简单设定

很多师弟师妹问到：”师兄， 我见上一节计算的时候打开自旋极化的时候，设置MAGMOM怎么没有讲啊？”

为了解答师弟师妹们的疑惑（主要是给师妹讲的），本节介绍一下INCAR中MAGMOM的设置以及易出错的地方。

MAGMOM 的设定

谈到自旋极化计算，避免不了的就是MAGMOM这个参数。通过MAGMOM我们可以指定体系中原子的初始磁矩。通俗点说，就是给VASP指条计算的明路。对于复杂体系来说，合理的初始值可以加快计算速度，并保持计算结果的正确性。但是对于一些简单的磁性体系，我们可以直接使用ISPIN=2, MAGMOM不必进行设置。

我们通过上图可以知道：O$_2$ 分子的基态是三重态，也就是每个原子都有一点单电子，所以每个氧原子的初始的磁矩可以设置为1。但VASP对MAGMOM的要求不一定非得是1，也就是初始值是一个可以模糊也可以精确的数值。因为很多时候，我们不知道确切的磁矩是多少。如果我们已经知道磁矩是多少（比如O$_2$），可以直接写上去，也可以差不多写一个； 但如果不知道的话，就需要根据自己的假设，猜一个初始值了。

下面通过一个问题，先简单看下MAGMOM 应该怎么设置。 对于O$_2$分子的计算，以下几行中，MAGMOM设置错误的是第几行？

1
2
3
4
5
6
7
8
9
MAGMOM = 2 * 1 
MAGMOM = 1 * 2 
MAGMOM = 1  1 
MAGMOM = 2  2
MAGMOM = 3  3 
MAGMOM = 2*1
MAGMOM = 1*2
MAGMOM = 1.5*2
MAGMOM = 2*1.5

有兴趣的可以尝试一下前面中的各个选项。答案是：1，2，7，8 行。当你用了1，2，7，8行中任何一个时，会得到这样的错误信息：

（此时的你，应该知道从哪里找到的这个错误信息，不知道的请查看前面几节的内容）

显示的错误告诉我们MAGMOM设置的不合理。为什么会出错呢？ 我们看下官网MAGMOM的说明：

敲黑板：

  • 首先：MAGMOM是一个实数的排列（real array）

  • 其次：注意红色框框中的部分。默认值是，原子的个数乘以1.0，也就是原子个数在前面，MAGMOM的值在后面。

  • 最后： * 前后没有空格，没有空格，没有空格！！！

因此前面的问题中：

第1行存在1个错误：* 前后有空格

第2行存在2个错误：* 前后有空格，原子数目和MAGMOM值颠倒了

第7和8行存在1个错误：原子数目和MAGMOM值颠倒了

有的人会问，那其他行中，比如第5行设置的 3 3 难道也正确吗？ 正确的磁矩不应该是1 1 吗？

是的，这是正确的，看官网的话：

If one is searching for a spin polarised (ferro- or antiferromagnetic) solution, it is usually safest to start from larger local magnetic moments, because in some cases, the default values might not be sufficiently big. A safe default is usually the experimental magnetic moment multiplied by 1.2 or 1.5.

  • MAGMOM并不要求严格按照我们已知的数据去输入

  • 如果你知道体系的磁矩是多少，初始的时候可以设置的更大些，1.2或者1.5倍。

  • 如果你设置的很大，通常情况下（MAGMOM = 5 5 或者MAGMOM = 2*5，或者更疯狂：MAGMOM = 2*10） ，VASP会自动矫正回来。但是大师兄建议1.5倍足够了。

  • 如果你不知道体系的磁矩是多少，可以根据原子所处的化学环境， 根据成键情况，大体推测有多少个未成对电子，然后将未成对电子数目*1.5即可。

  • VASP的wiki版中：* 前面有空格，这是不对的。

当我们正确解决了POSCAR的结构搭建，懂得了POTCAR和POSCAR的关系，以及如何设置MAGMOM时，就可以计算O$_2$的单点能量了。当然，本例中MAGMOM不用设置，因为默认值就是1，如果你的体系很简单，那么直接ISPIN = 2 就可以了，MAGMOM可以不管。

示例：

大师兄计算了一个复杂大分子在Ni(111)表面上的吸附情况。有56个Ni，17个C，20个H，以及6个O原子。前面已经学到，当体系中含有Ni的时候，要考虑自旋的情况。 Ni(111)表面是一个简单的自旋体系，完全可以使用默认值。因为体相的Ni磁矩为0.56 μB左右。默认值为1.0，已经足够大。为了让大家有一个更加清晰的印象。大师兄专门设置了一个MAGMOM。 如下图：

INCAR 中MAGMOM的设置：

通过该例子，相比大家对MAGMOM的设置以及该注意的地方已经熟练掌握了，当然还有非线性的MAGMOM设置，这个在后面详细描述。

总结

1） 对于简单体系来说，MAGMOM可以采用默认值；

2） MAGMOM设置的时候，初始值不要求与实验值完全一致，一般取大些（1.5倍）比较好。

3） MAGMOM磁矩中*前后没有空格。

LVASPTHW

Ex11 O$_2$单点计算和优化结果分析

前面几节，我们讲解了O分子单点计算和O$_2$的POSCAR和POTCAR的准备。这一节我们主要讲解一下：

1）O$_2$的分子结构分析

2）如何初步进行构型优化计算。

O$_2$ 的单点计算

首先解释下什么是单点计算：顾名思义就是不优化结构，直接算个能量，电子相关的性质。你也可能会听到很多人说静态计算，或者自冾计算。其实都一样的：几何结构计算前后不发生变化。

首先提交O$_2$ 静态计算任务，运行，等待结束。如果你的任务出错，请跟大师兄的输入文件进行对比、改正，直至计算正常结束。如下图：

任务结束后查看OSZICAR:

从OSZICAR最后，得到体系的磁矩为2μB，你应该知道这个磁矩是怎么回事，由哪两个电子贡献。

如果不知道，看下图O$_2$的分子轨道结构：

看完该图，相信大家对于O$_2$的成键方式有了一个更加深刻的印象。我们对比下VASP的输出结果。首先分析α电子的排布情况：

在这里，你会发现能带3和4 是简并的，应该是π(p2x）和π(p2y） 轨道中的α电子。能带 5 对应的应该是σ(2pz) 的电子。在O$_2$分子的电子构型中，两个O原子的2pz轨道以头碰头的方式形成一个σ键，其能量要比2px（2py）以肩并肩方式形成的π键能量要低。但是，能带5的能量（-13.3126）比3和4的（-13.3870）要高些，这与O$_2$的电子构型不一致，表明VASP的单点计算结果是不可靠的。

再看一下β电子的排布情况：

能带4和5对应的是应该是π（2px） 和π（2py） 轨道中的β电子。且总的轨道能量与前面图中一致，这说明VASP对β电子的描述是合理的。

为什么出现这样的情况呢？难道跟前面O原子的计算一样，VASP又算不准啦？

不是的，VASP怎么着也是个老牌的，响当当的计算程序，总不能让人每天指着鼻子骂算不准！

这里的主要原因是：来自于实验的键长值未必就是计算程序所认可的。也就是说，实验值和理论值之间存在偏差，实验的结构不能直接用来计算其性质，只可以作为一个理想的初始值。所以，O$_2$的分子结构需要优化一下。

VASP 优化分子结构

VASP优化分子结构的时候，需要用到一个参数：IBRION。引用官网的话：IBRION determines how the ions are updated and moved. 也就是说IBRION 这个参数决定了结构的优化过程。当你去官网查看的时候（google 搜索VASP IBRION这两个关键词），会发现IBRION有很多值。

想要正确进行计算，你就需要去硬着头皮去了解各个值的含义了，这个过程必须自己去做，只听别人的建议去设置参数，而不自己去主动学习的，你的能力永远不会得到提升！！！

一般来说，优化结构的时候有3个选择：

IBRION=3：你的初始结构很差的时候；

IBRION=2：共轭梯度算法，很可靠的一个选择，一般来说用它基本没什么问题。

IBRION=1：用于小范围内稳定结构的搜索。

如果你的体系遇到结构不收敛的时候，首先检查自己的结构是否合理，也就是物理化学意义是否清晰。如果结构没问题，可以尝试下换下IBRION的参数。

下面，我们在INCAR中加上IBRION参数（IBRION=2），其他输入文件保持不变，重新进行计算：

如果是用上图中的INCAR，你会发现任务很快就算完了。而且只有一步，难道输入的结构就是VASP计算出来的的稳定结构吗？有这种可能，但几率极低。

如果我们仔细查看下OUTCAR中的电子构型，发现它的信息和前面的单点计算一样。这说明，vasp并没有优化，而是又运行了一次单点计算。

为什么呢？ 这是因为另一个参数：NSW。

NSW 控制几何结构优化的步数。也就是VASP进行多少离子步。

官网查看下NSW选项，发现默认值是0，也就是没有进行优化。（默认值，也叫缺省值，英文里面是 Default。 意思是，如果你不输入这个参数，程序将默认使用XXX的数值）

现在原因找到了，继续进行优化任务。问题来了：NSW怎么设置呢？

  • 首先，它必须是大于等于0的整数。
  • 其次，一般来说，简单的体系200步内就可以正常结束。
  • 不知道什么时候收敛，初始结构很差，或者设置了很严格的收敛标准，那么你就要增大一下NSW的取值了，比如NSW=500或者更大。
  • 我们的这个O$_2$例子很简单，设置了NSW=10（你也可以设置为100,200或者500，不会影响计算结果的。）
1
2
3
4
5
6
7
8
System = O
ISMEAR = 0
SIGMA  = 0.01
ISPIN  = 2
MAGMOM = 2*2

IBRION = 2
NSW = 10

计算完成后，打开OSZICAR：

可以看到，结构优化进行了3步便停止了（如果你设置了NSW=1000，那么也是3步后就结束）。其中，每一步内又包含了若干电子步。此时的你应该知道是什么参数控制优化的结束，如果不知道请查看前面Ex09中关于收敛的文章。思考一下： 同样优化200步，设置EDIFF=1E-7和EDIFF=1E-4会有什么区别呢？

查看下OUTCAR：

你会发现，优化过后，OUTCAR中α电子的占据状态调整过来了，β电子的保持不变。这说明计算成功了，优化起作用了。

那么优化过后的结构怎么查看呢？键长又是多少呢？

下面我们要开始真正掌握VASP的输出文件了：CONTCAR。 CONTCAR是VASP的一个输出文件，它包含了VASP计算中最后一步几何优化的结构信息，也就是优化完的结果。它也是文本格式，可以直接打开查看，如图：

怎么才能知道优化完的O-O键长是多少呢？

1）通过坐标直接算：

此时，要注意CONTCAR输出的是Direct坐标，也就是分数坐标，需要转换成笛卡尔坐标。

xy两个方向不用考虑（都是0），z方向的坐标相减即可：键长保留四位就足够了。

(0.136785-( -0.001953))*8.9 = 1.2348 （$\AA$）

2) 使用可视化软件测量：

这样手动算，简单来说还可以，等复杂了就麻烦了，幸运的是，常用的建模软件中：p4vasp，VESTA，VMD，MS，ASE-gui等，都有测量键长的办法。这里简单讲解一下p4vasp 的用法。

p4vasp的安装

这个详见本书的附录2：https://www.bigbrosci.com/2017/11/18/A02/

1） Ubuntu 系统： sudo apt-get install p4vasp

2） Windows系统： 下载p4vasp的安装包，解压后直接打开即可。

下面是Windows下通过p4vasp读取CONTCAR的方法。Linux的也可以这样做，但终端里面：p4v CONTCAR 更高效直接。打开VASP后的界面如下： 注意图中： System：后面的三个问号部分，后面会进行对比：

点击左侧栏中的 Open 选项。下图中，左侧点击进入CONTCAR所在的目录，进入后，在右侧会显示CONTCAR

双击右侧的CONTCAR, 你会发现之前的三个问号？？？部分发生了变化：

显示的 O 表明VASP读取CONTCAR成功。

点击左侧的Show按钮，就可以查看结构了。

可视化界面的基本操作

鼠标左键按住不放，可以3维空间旋转结构；

鼠标中间摁住不放，可以上下，左右移动结构：

鼠标右键摁住不敢，动动鼠标可以缩放结构；

选择合适的观察位置：

1 把鼠标移动到你要选中的原子上，

2 通过空格（键盘上最长的键）用来选择或者取消选择原子：

3 选中两个氧原子后，如下图：

选中之后，点Structure –>Measure 显示键长为1.234774 Å，和我们手动计算的结果一样。

当然，也可以先点击：Structure –>Measure，然后再选择感兴趣的原子。

小结一下：

实验值为1.2075 Å，VASP计算结果为1.2348 Å。两者之间的差值为： 0.0289 Å，偏差为：(1.2348-1.2075)/1.2075 = 2.26%。对于理论和实验之间的偏差，如果小于5%，我们一般可以认为吻合的很好。有时候很多同学揪着VASP的计算结果与实验值的偏差不放，误差已经千分之几了，感觉心里还是不放心，有着一种不完全匹配不罢休的冲动。这大可不必放在心上，如果你的结果偏离实验值千分之几，直接用就可以了。

“To err is human; to describe the error properly is sublime.”
— Cliff Swartz, Physics Today 37 (1999),388.

对于其他的软件程序（MS，VESTA，VMD等），大家下载安装后，百度里面搜一搜教程，基本操作应该很快就能掌握。

扩展练习：

1 IBRION：https://cms.mpi.univie.ac.at/vasp/vasp/IBRION_tag_NFREE_tag.html

2 NSW：http://cms.mpi.univie.ac.at/vasp/guide/node108.html

3 CONTCAR：https://cms.mpi.univie.ac.at/wiki/index.php/CONTCAR

4 从头开始重现本节的所有操作；

5 尝试不同的初始键长，运行vasp，查看输出结果。

总结：

1）学会结合自己所学的化学知识，分析双原子的电子构型；

2）知道IBRION + NSW进行结构优化；

3）知道什么参数控制结构优化的停止，以及单个离子步内电子步数；

4）学会使用可视化软件查看输出的几何结构；

5）知道理论结果和实验值之间没有100%吻合。

LVASPTHW

Ex10 O$_2$结构搭建以及POTCAR注意部分

前面O原子的能量计算已经告一段落，下面我们把体系变得更加复杂些：O${_2}$分子。扩展练习已经做的同学，相信现在已经完成了O${_2}$分子的静态计算。本节我们主要介绍一下:

  • O${_2}$分子模型搭建的细节;
  • POTCAR和POSCAR的对应关系，注意事项。

搭建模型

O${_2}$分子模型的搭建其实很简单，将两个原子在一个直线上连起来就行了。 键长大小可以书本上查， 数据库里找，也可以去网上找些参考资料。大师兄推荐常见常用的数据库:

1）CRC hand book:

  • PDF 电子书，这个网上到处都有，大家自己搜搜就能找见。具体在1403 页(本书有2661页，建议大家查询的时候先把左侧的目录展开，要不然实在是太痛苦!)

2） NIST 数据库: http://cccbdb.nist.gov/exp2x.asp

3) Wikipedia，谷歌和百度:

在谷歌或者百度里面直接搜索O$_2$、键长这两个关键词，就可以得出很多的链接，比如下面两个。

https://zhidao.baidu.com/question/101615129.html

http://www.science.uwaterloo.ca/~cchieh/cact/c120/bondel.html

大师兄你特么在逗我嘛？

氧气原子这么简单的分子，随便一翻，当然到处都能得到结果啦!

哈哈，没有逗你。在这里，大师兄想要告诉你的是：

  • 搜索的时候尽量用谷歌，而不用百度；
  • 尽量搜英文，不要输汉语，这样你得到的信息会相对多一些。

师兄，你又在逗我，谷歌我们访问不了，怎么用啊？

没有逗你，现实是残酷的，没有google，科研会被百度严重拖后腿。虽然我也不想批评国产，谁让它不争气呢？

所以，科学上网的技巧，你要get到。

划重点：

前面介绍了一堆，主要目的有3个：

1) 给大家提供一个数据查询手册(CRC handbook)和数据库(NIST)。当然根据自己的研究方向，还有其他更加专业的数据库，随着学习的深入，后面会讲解。 此外，写论文的时候，如果实在不想找相关的参考文献，可以直接引用CRC hand book这本书（前提是你在这本书里面找到了相关的数据）

2) 查找数据库另一个重要的原因是：如果你的初始结构，比较理想，这会大大加速你的计算过程，还会避免一些意想不到的计算错误（后面的章节会具体介绍到）。

3）除了查找数据库，别人发的文章也是初始结构的重要或者主要来源。

模型搭建

知道了键长信息后，开始结构搭建O$_2$的模型，我们知道它是直线型分子（算大师兄没说），其中一个原子在原点的位置上了，那么我们在三个方向上随便找个坐标就可以了。修改POSCAR如下:

两点要注意：

  • 第7行原子个数：把1改成2 ，表示氧元素有2个原子。

  • 第10行添加第二个氧原子的坐标：

方向 坐标
沿x轴 1.2074 0.0 0.0
沿y轴 0.0 1.2074 0.0
沿z轴 0.0 0.0 1.2074

上图中为沿着z轴方向。注意：

  • xyz三个方向的数值之间可以有1个或者若干个空格，不要用tab分割数字；
  • 在搭建模型的时候，脑子里要有一个立体的概念，原子在三维方向上的排列，移动变化等。

关于POTCAR

大师兄，氧原子多了一个，POTCAR怎么办？是不是每个原子对应一个POTCAR啊?

不是的。POTCAR 是根据POSCAR中的元素顺序创建的(第6行)，与原子数目无关。

第6行是体系中的元素，只有O元素，所以我们的POTCAR还是用之前O原子练习的那个;

第7行是与上一行相对应的体系中元素的原子数。

但是！在你的POSCAR中，如果把O写了两遍，如下图：

此时，第6行中有2个O，且第7行中有两个O的原子数目分别为1。 POTCAR中就要对应的两个O原子的Potentials!

主动出错：

如果我们使用上图的POSCAR（O元素写了2次），且POTCAR中只有1个O原子的Potential。

提交任务后，瞬间完成。这么快就算完肯定不是因为你的服务器多么的牛逼! 而是因为出错了!!

那么我们需要找到错误的原因，怎么去找错误信息呢?

1) 查看OUTCAR

VASP就这么几行，从这里看不出来是哪里错了。

2) 一般来说，VASP计算的时候，会生成两个额外的文件，一个是关于服务器集群计算error的，另一个是VASP运行的out文件。这个out和OUTCAR不太一样，记录着VASP的运行过程和出错信息。看大师兄的文件目录:

1
2
3
4
5
qli@bigbro:~/test$ ls
CHGCAR   DOSCAR   e_Single.142845   INCAR   o_Single.142845  OUTCAR  POSCAR  PROCAR  sub12    WAVECAR  CHG      CONTCAR  EIGENVAL  IBZKPT   KPOINTS  OSZICAR   PCDAT   POTCAR  REPORT  vasprun.xml  XDATCAR
qli@bigbro:~/test$ cat e_Single.142845
qli@bigbro:~/test$

不同的组可能命名不一样，但一般都会有这两个文件。这里的e对应的是服务器的出错信息，o对应的是VASP的out文件。 首先打开查看下服务器出错文件：(图中cat 命令) 发现什么都没有输出，说明服务器没有出错。那我们打开一下 o_Single.142845 文件，如下图

在这个文件里面最后一行，给出了错误的信息：

一般出现这个错误的时候，你就要去检查POSCAR和POTCAR中的元素是否对应了。

解决问题

本例中的错误该怎么解决呢？ 既然POSCAR和POTCAR不一致，解决的话，有2个办法：要么变POSCAR，要么变POTCAR。

1）变POSCAR：

  • 第六行中：把O O换成O;
  • 第七行：氧原子数目改成2，就如刚开始的POSCAR。
  • POTCAR保持不变即可（POTCAR中只有1个O的Potential）:

2）变POTCAR

  • 使其中的元素与POSCAR中的一致，也就是有两个氧元素的Potentials！
1
2
3
4
5
6
7
qli@bigbro:~/test2$ grep TIT POTCAR
   TITEL  = PAW_PBE O 08Apr2002
qli@bigbro:~/test2$ mv POTCAR POTCAR1
qli@bigbro:~/test2$ cat POTCAR1 POTCAR1  > POTCAR
qli@bigbro:~/test2$ grep TIT POTCAR
   TITEL  = PAW_PBE O 08Apr2002
   TITEL  = PAW_PBE O 08Apr2002

这样，再运行就不会出错了。

POTCAR 的制备

前面的演示中，大师兄教给你了POTCAR的制备方法: 使用cat 命令将2个POTCAR1 连在一起然后输出新的POTCAR。如果看不明白，继续看下面的操作：

如果你的体系中含有其他的元素，比如: Fe C H O

那么你就要先准备这四个元素的POTCAR：POTCAR-Fe、POTCAR-C、POTCAR-H、POTCAR-O，

然后运行命令:cat POTCAR-Fe POTCAR-C POTCAR-H POTCAR-O > POTCAR 就可以了。

具体操作：
1
2
3
4
5
6
7
8
9
10
11
12
13
qli@bigbro:~/test2$ ls
POTCAR-C  POTCAR-Fe  POTCAR-H  POTCAR-O 
qli@bigbro:~/test2$ grep TIT POTCAR*
POTCAR-C:   TITEL  = PAW_PBE C 08Apr2002
POTCAR-Fe:   TITEL  = PAW_PBE Fe 06Sep2000
POTCAR-H:   TITEL  = PAW_PBE H 15Jun2001
POTCAR-O:   TITEL  = PAW_PBE O 08Apr2002
qli@bigbro:~/test2$ cat POTCAR-Fe POTCAR-C POTCAR-H POTCAR-O > POTCAR
qli@bigbro:~/test2$ grep TIT POTCAR
POTCAR:   TITEL  = PAW_PBE Fe 06Sep2000
POTCAR:   TITEL  = PAW_PBE C 08Apr2002
POTCAR:   TITEL  = PAW_PBE H 15Jun2001
POTCAR:   TITEL  = PAW_PBE O 08Apr2002

这四个元素POTCAR数据从哪里找，在哪个目录下面？

这就得问问你们组里的师兄师姐们，或者老师了，网上也有很多。

扩展练习:

1 正确运行O$_2$分子的静态计算;

2 查找相关的O$_2$分子的轨道排布，并分析结果的合理性;

3 主动制作错误的POSCAR和POTCAR文件，运行查看结果和错误。

总结

1 整理自己研究方向相关的数据库资料；知道去哪里查询信息；

2 学会怎么根据已知的结构参数搭一些简单的分子模型: CO，H$_2$O， N$_2$， H$_2$ 等;

3 知道去哪里找出错文件;

4 熟练掌握通过 cat命令制备POTCAR的方法;

5 必须掌握POSCAR和POTCAR的对应关系。

LVASPTHW

Ex09 学会收敛点!

前面一节我们学会了正确计算氧原子的能量并分析电子在轨道中的占据情况，这一节稍作总结补充。

复习上节内容

在VASP计算的O原子的电子结构讨论中，来自不莱梅大学的群友，画了一幅画，大家可以看下:

VASP的结果分析完毕后如下:

看画填数，照猫画虎：

由此可见，VASP对于氧原子的描述不是很准确（不理解的请看上一节的详细描述)。而这一点，VASP官网你是找不到的。因此，在计算的时候，对于自己的体系有一个清晰的化学物理印象很重要。VASP的结果也要有针对的去判断。有的时候，即使计算收敛了，但没有体现出任何正确的物理或者化学意义，那就是单纯的数学收敛，结果是不可取的。VASP是一个计算的工具，如何正确分析和判断它的输出才是关键的。这一点我们在后面的学习中也会反复强调。

VASP是怎么判断收敛的?

我们看一下VASP的迭代计算过程：

​

  • 首先，它会猜一个初始的电子密度，然后据此计算体系的势能，求解KS方程，并给出体系的总能量以及对应的电子密度，也就是我们之前说的电子步的优化。

  • 将最新一步的结果与前面进行对比，当前后两者的差值达到我们预设的收敛标准时，计算结束。

  • 这个预设的收敛标准在VASP中通过两个参数来描述：EDIFF 和 EDIFFG。

  • EDIFF 控制电子步（自洽）的收敛标准。在O原子的计算中，由于我们不需要优化，直接进行静态计算，完全由EDIFF控制计算的收敛情况。EDIFFG后面我们再介绍。

从图中可以看出：第19步和18步的能量差为：-0.3484 E-04，停止迭代，开始输出结果。这是因为VASP计算中EDIFF的默认值是 1 x 10-4。

鱼和熊掌的关系

那么：

是不是精度越高，计算越准确呢？ 这是肯定的！

是不是计算都需要这么高的精度呢？ 这肯定不是的！

精度高（鱼：<・)))><<）意味着需要更多的迭代次数，也就是需要更多的计算时间（熊掌）。下图是收敛标准从默认值的1E-4 降低到1E-7后的收敛情况：单迭代次数从19增加到了27，相当于增加了原来1/2的工作量。所以，精度太高，计算量会增加。应了那句流行语：请在wifi下观看，土豪随意。

我们看一下扩展阅读的内容:

第一点：官网说（红色方框中内容）：如果我们把收敛标准设置成0，那么迭代会永远进行下去。大师兄很感兴趣，便设置：EDIFF = 0测试了一番，结果如下:

可以看到，在45步迭代之后，精度收敛到1E-12，VASP便停止了，而不是所谓的always! VASP又把俺们给骗了!

骗归骗，但计算量的增加确是铁打的事实，从1E-4的19步，增加到了45步。

思考另外一个问题：既然精度从默认值1E-4提高到1E-12 (8个数量级)，那么我们算出来的氧原子能量有什么变化呢？看下图:

1
2
3
4
5
6
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex09$ ls
0  4  7
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex09$ ta.sh
0        -1.89224598
4        -1.89218308
7        -1.89224603

图中的0， 4， 7 分别代表EDIFF为 0， 1E-4和 1E-7 的情况，我们可以看出，这三个计算结果中，能量随着EDIFF的减小变化甚微。当然，这也与我们简单的测试体系有关系。倘若计算N多个原子的体系时，不同的精度间的差值会有所增大。上面图中使用了一个小脚本：ta.sh，内容如下：

1
2
3
4
5
6
#!/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

如何看VASP的说明书

那么我们再看一遍VASP官网的说明，注意底部的Mind部分:

大家在浏览VASP官网的时候，凡是带有Mind，Important字样的，它们后面的文字一定要认真阅读，把握其中的含义。因为这都是大家常见的疑难点以及易错的地方。收敛速度跟EDIFF的值成指数关系，在大多数的情况下，1E-4足以胜任，没有必要采用其他的数值，如果大伙感觉不放心，大师兄建议1E-5即可。

几何优化的收敛标准

目前还没有介绍几何优化，先简单介绍下EDIFFG这个参数，大家知道即可，后面学习到了再回来看。

我们优化几何结构的时候，当结构前后变化达到我们设定的要求时，便停止优化。而这个设定的要求，就是通过EDIFFG来控制。

  • 对于优化，我们可以使用力作为收敛标准，此时EDIFFG为负值。一般来说取值在-0.01到-0.05之间(-0.01对于力收敛来说已经是一个很严格的要求了)。

  • 当然，对于较大的体系，我们也可以使用能量作为标准：此时，EDIFFG 为正值，一般取值范围在0.0001-0.001即可。

    注意:

不要把正负值弄错了。大师兄见过一个群友这么设置的: EDIFFG = -0.0001。然后在群里咨询为什么他的计算不收敛。这就好比在没有WiFi的时候，装成土豪看视频，卡里的钱扣完了，视频才看到一半…..前面说到，EDIFFG = -0.01的时候收敛就已经够费劲了， 你再加个0，莫不是要算到天荒地老？

小结一下

EDIFF，EDIFFG 是控制收敛标准的两个参数。前者负责电子自洽过程(单个离子步内)，取值为1E-4或者1E-5即可，没有特殊要求，不建议设置的太低。后者负责结构优化的过程(可正可负)。对于EDIFFG，默认值采用能量收敛，标准为:EDIFFx10。

扩展训练:

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

2 在1的页面中，左侧框中搜索EDIFFG，查找相关参数;

3 继续查询之前了解的相关参数，学会使用改网址;

4 根据前面所学，进行O$_2$分子的静态计算，并分析其轨道结果，判断是否合理。

总结

  • 重温了一下上节的结果，不要完全相信程序的结果，物理，化学意义更为重要;
  • 初步了解VASP的迭代过程，
  • 必须掌握EDIFF在电子步中的作用，取值大小与收敛速度的联系，也就鱼和熊掌的关系；
  • 死死记住EDIFFG代表的含义，取值可正可负及其默认值;
  • 经常浏览VASP官网，尤其是Mind，Note， Important等后面的部分;（眼睛要尖！）
  • 建议大家浏览VASP的Wiki版网页(左下方原文链接)，搜索相关参数，里面的公式较之前模糊的形式，有了很大的改进。
  • 目前VASP官网国内上不去，掌握科学上网的技巧也是做科研的必须技能。
个人经验

Ubuntu18.使用VESTA报错的解决办法

今天在Ubuntu18的版本上用了一下VESTA，发现17.XX版本的问题还依然存在。

当你敲命令：VESTA之后，会得到包含下面信息的错误。

error while loading shared libraries libpng12.so.0

解决步骤：

1）点下面的链接，下载libpng12的deb文件。

https://packages.debian.org/jessie/amd64/libpng12-0/download

2）双击下载的deb文件，安装即可。

3）再用命令打开VESTA就不会出错了。

个人经验

通过关键词判断优化任务收敛

今天QQ群里，有人（小风）问怎么通过reached这样的关键词来判断自己的任务是否收敛。

首先：阿牛哥哥说说了，实现的方法有很多种？这大实话一点都不假！

然后：乐平老师说，问问题的人没有问用什么语言写？这Python高手深藏不漏。

紧接着：杠精大神用bash语言coding了一番，敲出若干行神秘的文字：

1
2
3
4
5
6
7
check=$(grep "reach" vasp.log | tail -1)
if [ "$check" == "" ];then
   echo "$n $i $j cp POSCAR CONTCAR" >> ../conver.log
   cp POSCAR CONTCAR
else
   echo "$n $i $j $check" >> ../conver.log
fi

下面我们就通过学习赋值和if语句来揭开上面的神秘文字。

原理

如果：f(A) = B，且f(C) = B，那么A的某一个性质 = C的某一个性质， 其中f() 是对操作对象在固定范围所实施的一个固定的方法。

也就是说对A和C进行同样的操作，同时都得到了B，那么说明A的某一个性质 = C的某一个性质 。

操作

我们写脚本的目的是要通过B来判断A与C的关系。

A是我们已知的结果；

C是我们未知的结果；

B也就是我们说的赋值的对象。

首先，我们看一个已经收敛的例子：

1
2
3
4
5
6
qli@bigbro:~/ch4$ grep reached OUTCAR
------------------------ aborting loop because EDIFF is reached ----------------------------------------
------------------------ aborting loop because EDIFF is reached ----------------------------------------
------------------------ aborting loop because EDIFF is reached ----------------------------------------
 reached required accuracy - stopping structural energy minimisation

我们知道，如果使用grep reached OUTCAR这个命令，我们会得到2个结果，一个是aborting loop because EDIFF is reached，暂且记为A1， 另一个是reached required accuracy - stopping structural energy minimisation ，暂且记为A2。如果一个优化任务收敛的话，我们是通过A2来判断的。那么现在我们就可以把A2名字改成A了。A有什么特点呢？ 我们随便列2个。

1） A在OUTCAR中只出现一次，且出现在最后一行。

2） A的第一个单词是reached。

现在，我们将reached这个单词赋值为B。也就是让B为收敛冠名，只要满足通过上面2个条件筛选出来的结果是B，就说明计算收敛了。

筛选这个操作怎么进行呢？

  • 利用特点1）： 使用tail 这个命令。
1
2
3
4
5
6
7
8
qli@tekla2:~/ch4$ grep reached OUTCAR
------------------------ aborting loop because EDIFF is reached ----------------------------------------
------------------------ aborting loop because EDIFF is reached ----------------------------------------
------------------------ aborting loop because EDIFF is reached ----------------------------------------
 reached required accuracy - stopping structural energy minimisation
qli@tekla2:~/ch4$ grep reached OUTCAR  |tail -n 1
 reached required accuracy - stopping structural energy minimisation

  • 利用特点2）：使用cut或者awk命令来获取第一个单词。

    1
    2
    3
    4
    5
    6
    7
    qli@tekla2:~/ch4$ grep reached OUTCAR  |tail -n 1
     reached required accuracy - stopping structural energy minimisation
    qli@tekla2:~/ch4$ grep reached OUTCAR  |tail -n 1 | cut -c 2-8
    reached
    qli@tekla2:~/ch4$ grep reached OUTCAR  |tail -n 1 | awk '{print $1}'
    reached

  • 判断C = B ? 

    1
    2
    3
    4
    5
    6
    7
    8
    9
    qli@tekla2:~/ch4$ b='reached'
    qli@tekla2:~/ch4$ c=$(grep reached OUTCAR  |tail -n 1 | awk '{print $1}')
    qli@tekla2:~/ch4$ echo $c
    reached
    qli@tekla2:~/ch4$ c=`grep reached OUTCAR  |tail -n 1 | awk '{print $1}'`
    qli@tekla2:~/ch4$ echo $c
    reached
    qli@tekla2:~/ch4$ if [ c=b ]; then echo bigbro ; fi
    bigbro
LVASPTHW附录

Ubuntu下p4vasp的安装

使用vasp计算，后处理的小程序、小脚本有很多。其中最为出名的非p4vasp莫属了。从搭建模型，调整结构，DOS、能带计算，功函数等等，都可以使用它来完成。举个例子，本人博士四年，一直用它来查看和处理结果，再结合一些其他辅助的小程序和脚本，可以说是极大地提高了自己的工作效率。今天因写脚本需要，浏览了一下p4vasp官网，发现官网对p4vasp的安装解释非常简单。于是便测试了一下。本节主要包含3部分:

  • 第一部分：p4vasp官方的安装流程，相信大家按照这一步都可以完成安装；
  • 第二部分：本人亲自安装的具体流程；
  • 第三部分：傻瓜化安装流程。

第一部分

p4vasp官网： http://www.p4vasp.at/。点击下方的链接，根据操作进行：

http://www.p4vasp.at/#/doc/getstarted

第二部分：动手安装流程：

这一部分内容，是根据前面的官方介绍，亲自操刀实践的部分。

1） 首先说明下，本人的操作系统是：Ubuntu 16.04.3 LTS

2) 下载p4vasp：

目前，对Linux系统来说，最新的貌似只有source版本的了，binary的已经成为了历史。

对Windows用户来说，更可怜，p4vasp都没有更新，版本还是老样子，Windows的获取方法：

Linux用户则可以通过p4vasp官网下载：

3) 解压：

这里将p4vasp软件包解压到Documents目录下，当然，也可以使用官网介绍的解压命令进行。

4) 进入解压的目录，然后按照官网说的，Ubuntu用户，使用 install-ubuntu-dependencies.sh 这个脚本安装一些一些必须的dependencies. 如下图：运行命令后，会提示你是不是下载，输入Y 或者 y ，然后回车。屏幕此时开始疯狂地颤抖，说明正在下载安装。

注意：

  • 解压到什么地方不用要，可以放到桌面上，也可以直接解压在Download里面，本人习惯解压到Documents这个文件下。

  • 这里大师兄先进入的 install 这个目录， 然后运行的脚本。安装完记得再返回去。

5) 安装完成后，返回上一级目录，然后正式安装p4vasp。

等待屏幕颤抖完毕。

看到这里，基本上99.99% 成功安装了。

6) 测试一下： 终端直接输入： p4v ，效果如下：

打开p4vap了，说明安装成功。

7）检查下版本：（点击右上角的help）

最新版的p4vasp已经顺利安装成功。欢呼吧，少年！！

傻瓜化安装

在进行第二部分操作的前4年，本人一直用的这个办法：

1） 终端一个命令搞定：

1
sudo apt-get install p4vasp

效果如下：

输入Y，等待屏幕颤抖结束

2） 运行p4vasp：

3） 检查下版本：该版本为v0.3.29r1，比官网的最新版稍微旧一点。但不影响使用。

小结

对于Ubuntu用户：

  • 建议大家使用官网介绍的办法。安装最新版的p4vasp。
  • 如果遇到自己解决不了的问题，可以使用 sudo apt-get install p4vasp 这个命令一键搞定。
  • 其他linux系统的用户，本人没有尝试过。

Windows用户，目前只能使用本文中链接下载旧版本的了。

Mac 用户：浏览p4vasp的github网站，根据Readme文件自己尝试着安装吧。

https://github.com/orest-d/p4vasp

LVASPTHW附录

通过Ubuntu连接服务器的方法

本教程主要讲述的是如何在Ubuntu16.04操作系统上连接国科智算超算。分3部分：

  • 1）前期准备；
  • 2）连接和挂载超算中心；
  • 3）数据传输测试。

在使用本教程前，首先确认已经收到了超算管理员分配的秘钥。通过Vim或者其他文本编辑工具，可以查看一下这个秘钥。超算中心好比是一个宝藏，这个宝藏的大门上有把锁，而开启这把锁的钥匙，就是管理员给我们的秘钥。

前期准备

下面的工作，一步一步来，相信99.99% 的人都可以顺利完成。

  • 为了方便理解，本人将管理员给的秘钥重新命名了一下。将下面的999改成管理员给你的那个数字。
1
qli@bigbro:~$ cp  ~/Desktop/id_rsa_gkzshpc999.gkzshpc999  .ssh/my_key

经过上面的一步，管理员给的秘钥就被我们命名成：my_key了。这好比是我们在钥匙上贴了个标签，开门的时候直接找这个标签对应的钥匙就可以了。

如果你的Ubuntu系统下没有 .ssh 文件夹，可以自己先建一个，然后再运行上面的命令

1
qli@bigbro:~$ mkdir .ssh

  • 安装sshfs 

    1
    2
    3
    4
    5
    qli@bigbro:~$ sudo apt-get install sshfs 
    [sudo] password for qli: 
    Reading package lists... Done
    Building dependency tree       
    Reading state information... Done

  • 创建挂载国科智算的文件目录：

    1
    qli@bigbro:~$ mkdir gkzs

  • 打开~/.bashrc文件，并添加下面的两行：

    1
    2
    alias gkzs='ssh -i /home/qli/.ssh/my_key  gkzshpc999@59.49.37.9 -p 9236'
    alias mgkzs='sshfs -o IdentityFile=/home/qli/.ssh/my_key  -p 9236 gkzshpc999@59.49.37.9:  /home/qli/gkzs'

    上面我们的两行alias后面的内容主要为：

    • 我们的手（ssh）拿着钥匙（-i /home/qli/.ssh/my_key） 打开大门 （gkzshpc999@59.49.37.9 -p 9236）
    • 我们的手（sshfs）拿着钥匙 （-o IdentityFile=/home/qli/.ssh/my_key） 打开大门 （-p 9236 gkzshpc999@59.49.37.9）后，并将宝藏运回家（/home/qli/gkzs）

  • source 一下 ~/.bashrc 文件：

    1
    qli@bigbro:~$ . ~/.bashrc

连接超算中心

上面的工作完成之后，剩下的就是命令操作的事情了：

  • 连接超算中心：

    1
    2
    3
    4
    5
    6
    qli@bigbro:~$ ls
    Desktop  Documents  Downloads  gkzs  lvliang  Music  Pictures  Public  SCVPN  teklahome  Templates  Videos
    qli@bigbro:~$ gkzs
    Last login: Thu Nov 22 05:07:16 2018 from 80.29.50.15
    [gkzshpc999@login02 ~]$ ls
    perl5

  • 打开新的一个Terminal，我们挂载超算中心到我们的电脑上面，以便传输数据。

1
2
3
4
5
6
7
qli@bigbro:~$ ls
Desktop  Documents  Downloads  gkzs  lvliang  Music  Pictures  Public  SCVPN  teklahome  Templates  Videos
qli@bigbro:~$ ls gkzs/
 
qli@bigbro:~$ mgksz
qli@bigbro:~$ ls gkzs/
perl5

数据传输测试

  • 在连接到服务器的界面：我们创建一个文件：mount_test。

    1
    2
    3
    4
    5
    6
    [gkzshpc999@login02 ~]$ echo  'I love BigBro' >  mount_test
    [gkzshpc999@login02 ~]$ ls
    mount_test  perl5
    [gkzshpc999@login02 ~]$ cat mount_test 
    I love BigBro
    [gkzshpc999@login02 ~]$

  • 在挂载的目录下查看：目录下多出来了刚刚创建的 mount_test文件。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    qli@bigbro:~$ cd  gkzs/
    qli@bigbro:~/gkzs$ ls
    mount_test  perl5
    qli@bigbro:~/gkzs$ cat mount_test 
    I love BigBro
    qli@bigbro:~/gkzs$ 
    qli@bigbro:~/gkzs$ cp ~/Desktop/CONTCAR  .
    qli@bigbro:~/gkzs$ ls
    CONTCAR  mount_test  perl5
    qli@bigbro:~/gkzs$

上面的最后操作中，我们从桌面上复制了一个CONTCAR到挂载的目录下，然后查看下服务器的终端：

1
2
3
4
5
6
7
8
9
[gkzshpc999@login02 ~]$ ls
CONTCAR  mount_test  perl5
[gkzshpc999@login02 ~]$ head -n 5 CONTCAR  
Ru\(0\0\1)                              
   1.00000000000000     
     8.1377999784000004    0.0000000000000000    0.0000000000000000
     4.0688999892000002    7.0475415119999996    0.0000000000000000
     0.0000000000000000    0.0000000000000000   21.5631999968999999
[gkzshpc999@login02 ~]$

完成了上面的操作，下面你就可以将自己电脑上准备的一些计算文件或者文件夹通过命令复制到超算中心，然后就可以提交任务了。

LVASPTHW附录

通过Ubuntu连接国科智算的方法

本教程主要讲述的是如何在Ubuntu16.04操作系统上连接国科智算超算。分3部分：

  • 1）前期准备；
  • 2）连接和挂载超算中心；
  • 3）数据传输测试。

在使用本教程前，首先确认已经收到了超算管理员分配的秘钥。通过Vim或者其他文本编辑工具，可以查看一下这个秘钥。超算中心好比是一个宝藏，这个宝藏的大门上有把锁，而开启这把锁的钥匙，就是管理员给我们的秘钥。

前期准备

下面的工作，一步一步来，相信99.99% 的人都可以顺利完成。

  • 为了方便理解，本人将管理员给的秘钥重新命名了一下。将下面的999改成管理员给你的那个数字。
1
qli@bigbro:~$ cp  ~/Desktop/id_rsa_gkzshpc999.gkzshpc999  .ssh/my_key

经过上面的一步，管理员给的秘钥就被我们命名成：my_key了。这好比是我们在钥匙上贴了个标签，开门的时候直接找这个标签对应的钥匙就可以了。

如果你的Ubuntu系统下没有 .ssh 文件夹，可以自己先建一个，然后再运行上面的命令

1
qli@bigbro:~$ mkdir .ssh

  • 安装sshfs 

    1
    2
    3
    4
    5
    qli@bigbro:~$ sudo apt-get install sshfs 
    [sudo] password for qli: 
    Reading package lists... Done
    Building dependency tree       
    Reading state information... Done

  • 创建挂载国科智算的文件目录：

    1
    qli@bigbro:~$ mkdir gkzs

  • 打开~/.bashrc文件，并添加下面的两行：

    1
    2
    alias gkzs='ssh -i /home/qli/.ssh/my_key  gkzshpc999@59.49.37.9 -p 9236'
    alias mgkzs='sshfs -o IdentityFile=/home/qli/.ssh/my_key  -p 9236 gkzshpc999@59.49.37.9:  /home/qli/gkzs'

    上面我们的两行alias后面的内容主要为：

    • 我们的手（ssh）拿着钥匙（-i /home/qli/.ssh/my_key） 打开大门 （gkzshpc999@59.49.37.9 -p 9236）
    • 我们的手（sshfs）拿着钥匙 （-o IdentityFile=/home/qli/.ssh/my_key） 打开大门 （-p 9236 gkzshpc999@59.49.37.9）后，并将宝藏运回家（/home/qli/gkzs）

  • source 一下 ~/.bashrc 文件：

    1
    qli@bigbro:~$ . ~/.bashrc

连接超算中心

上面的工作完成之后，剩下的就是命令操作的事情了：

  • 连接超算中心：

    1
    2
    3
    4
    5
    6
    qli@bigbro:~$ ls
    Desktop  Documents  Downloads  gkzs  lvliang  Music  Pictures  Public  SCVPN  teklahome  Templates  Videos
    qli@bigbro:~$ gkzs
    Last login: Thu Nov 22 05:07:16 2018 from 80.29.50.15
    [gkzshpc999@login02 ~]$ ls
    perl5

  • 打开新的一个Terminal，我们挂载超算中心到我们的电脑上面，以便传输数据。

1
2
3
4
5
6
7
qli@bigbro:~$ ls
Desktop  Documents  Downloads  gkzs  lvliang  Music  Pictures  Public  SCVPN  teklahome  Templates  Videos
qli@bigbro:~$ ls gkzs/
 
qli@bigbro:~$ mgksz
qli@bigbro:~$ ls gkzs/
perl5

数据传输测试

  • 在连接到服务器的界面：我们创建一个文件：mount_test。

    1
    2
    3
    4
    5
    6
    [gkzshpc999@login02 ~]$ echo  'I love BigBro' >  mount_test
    [gkzshpc999@login02 ~]$ ls
    mount_test  perl5
    [gkzshpc999@login02 ~]$ cat mount_test 
    I love BigBro
    [gkzshpc999@login02 ~]$

  • 在挂载的目录下查看：目录下多出来了刚刚创建的 mount_test文件。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    qli@bigbro:~$ cd  gkzs/
    qli@bigbro:~/gkzs$ ls
    mount_test  perl5
    qli@bigbro:~/gkzs$ cat mount_test 
    I love BigBro
    qli@bigbro:~/gkzs$ 
    qli@bigbro:~/gkzs$ cp ~/Desktop/CONTCAR  .
    qli@bigbro:~/gkzs$ ls
    CONTCAR  mount_test  perl5
    qli@bigbro:~/gkzs$

上面的最后操作中，我们从桌面上复制了一个CONTCAR到挂载的目录下，然后查看下服务器的终端：

1
2
3
4
5
6
7
8
9
[gkzshpc999@login02 ~]$ ls
CONTCAR  mount_test  perl5
[gkzshpc999@login02 ~]$ head -n 5 CONTCAR  
Ru\(0\0\1)                              
   1.00000000000000     
     8.1377999784000004    0.0000000000000000    0.0000000000000000
     4.0688999892000002    7.0475415119999996    0.0000000000000000
     0.0000000000000000    0.0000000000000000   21.5631999968999999
[gkzshpc999@login02 ~]$

完成了上面的操作，下面你就可以将自己电脑上准备的一些计算文件或者文件夹通过命令复制到超算中心，然后就可以提交任务了。

LVASPTHW

Ex08 氧原子能量的正确计算

本节提要

前面的扩展练习中,要求大家结合自己的化学知识,判断计算结果是否正确。相比很多同学已经知道了氧原子含有磁性,我们在计算用需要考虑进去。在本节练习中,大师兄教给你如何分分析氧原子在轨道中的占据情况,判断计算是否合理,并正确进行氧原子能量的计算。

氧原子轨道

氧原子的电子构型是 $1s^2 2s^2 2p^4$ ,轨道占据情况如图:

所以氧原子含有两个单电子,具有磁性,且磁矩为2 $μ_B$。

计算结果

我们看一下前面的计算结果:

OSZICAR

OSZICAR

能量信息:

能量信息

电子占据情况: (OUTCAR 中查看)

电子占据情况

  • Band1 含有两个电子, 对应的 $2s^2$ ;
  • Band 2-4 为三个简并态,每个能带上有 1.33 个电子, 对应的 $2p^4$ 。

但在这里,电子的占据是不正确的:

  • 因为 $2p_x$ 轨道上含有两个电子, 而不是把 $2p_x$ 上的一个电子平均分配到 3 个 p 轨道上;
  • 此外,通过这个电子的占据情况,我们不能正确得到氧原子的磁矩信息。

那么怎么解决这个问题呢?

分两步解决问题

第一步:

由于氧原子具有磁性, 自旋极化计算要进行考虑, 在 VASP 计算中, 需要在 INCAR 中添加一项: ISPIN = 2

ISPIN 的取值为 1 或 2:
​ 1. 代表不考虑自旋极化,是VASP的默认值;
​ 2. 代表打开自旋极化该选项,计算中考虑。

自旋极化什么时候该加呢?

需要考虑自旋极化一般来说有以下几种：（新手记住前面4个即可。）

  • 单原子的计算,
  • O$_2$ 分子(基态为三重态)
  • 自由基相关的计算
  • 含Fe,Co, Ni 的体系
  • 要计算的体系具有磁性：顺磁，铁磁，反铁磁等，要打开自旋极化。
  • 当关注体系的电子性质时，且自己不知道加或者不加的时候，建议加上。

如果你不知道自己的体系是否要考虑自旋极化，可以简单做这样的一个测试：

  • 测试1：不加自旋极化，正常算，得到结果1；
  • 测试2：加上自旋极化，并在INCAR中添加： LORBIT = 11
  • 测试2中OUTCAR的末尾会输出各个原子的磁矩信息：
    • 如果体系中原子的磁矩不为0，那么需要考虑自旋极化，
    • 如果体系所有原子的磁矩都为0，那么就不需要考虑自旋极化，而且此时，你会发现测试1 和 2的结果中：结构和能量是一样的。

添加该选项后,重新计算,结果如下:

INCAR

1
2
3
4
SYSTEM = O
ISMEAR = 0
SIGMA = 0.01
ISPIN = 2

OSZICAR

OSZICAR

可以看出, 加了自旋极化后, 电子迭代步数增加了, 这是因为体系中的电子被强制分成了 $\alpha$ 和 $\beta$ 两种, 并分别进行计算, 因此增加了计算量。

体系的磁矩大小在 $OSZICAR$ 中查看, 对应右下角的 mag= 该项, 从这里我们可以看到, 磁矩信息与我们已知的 2 $\mu_B$ 吻合, 这一点是合理的, 正确的。

OUITCAR 中的能量信息:

OUTCAR

考虑自旋极化后,体系的能量降低了。

我们重点通过查看体系的电子占据情况, 判断该计算是否合理:

OUTCAR

从图中可以看出,结果由两个 spin 组成: spin component 1 和 2:
spin component 1 中含有 4 个电子, 另一个中含有 2 个电子:

OUTCAR

参考该图,可以看出 spin component 1 为自旋向上的 $\alpha$ 电子, 另一个是自旋向下的 $\beta$ 电子。

但是:

1. Spincomponent 1中,能带2-4中3个alpha 电子的能量是一样的。这不正确，因为其中一个2px轨道中的alpha电子已经成对,能量要比py 和 pz 低;描述失败!
2. spincomponent 2 中,另一个beta电子还是平均分配在了三个能带上面。

电子占据依然不合理。

上面的结果是由于体系的高对称性导致的简并所造成的。 $8 \times 8 \times 8$ 的立方体格子, 在这里具有高阶的点群对称性: $O_h$, 见 OUTCAR

OUTCAR

在 GGA 泛函中, 为了获得体系更低的能量,对于原子来说,通常会采用一种 symmetry broken solution 的处理方法。 但是在 VASP 计算中, 体系的对称性则是通过晶胞来获得, 即把这个晶胞当成一个原子来处理, 因此我们需要手动改变晶胞的形状来消除对称性造成的简并结果。 (此处的解释可能有些牵强, 有兴趣的可以查找与 symmetry broken solution 相关的文献)

INCAR, KPOINTS, POTCAR等均保持不变。修改POSCAR, 如下:

POSCAR

1
2
3
4
5
6
7
8
9
O
1.0
7.5 0.0 0.0
0.0 8.0 0.0
0.0 0.0 8.9
O
1
Cartesian
0.0 0.0 0.0

通过改变格子在三个方向的大小, 降低对称性。 大师兄还尝试了三斜的格子, 结果和这个是一致的。 (有兴趣的可以尝试下)

POSCAR

1
2
3
4
5
6
7
8
9
O
1.0
	7.5		7.6		0.0
	0.0		8.2564		8.2
	8.5		0.0		8.9
O
	1
Cartesian
	0.0		0.0		0.0

下面我们看一下改变晶胞为7.5 x 8.0 x 8.9后的计算结果:

OSZICAR

OSZICAR

OSZICAR 中电子迭代步数又增加了, 这是因为对称性的降低增加了计算量。

磁矩为 $2\mu_B$ ,这一点跟前面一样,正确!

体系的对称性: 从 $Oh$ 降低到了 $D{2h}$

OSZICAR

OSZICAR

能量信息:

TOTEN

TOTEN

跟前面能量相比, 可以看出, 对称性降低后, 体系的能量进一步降低了。

电子占据情况:

OUTCAR

TOTEN

分析

  1. 在 $\alpha$ 电子(spin component 1)中, 能带1 是 $2s$ 中的 $\alpha$ 电子; 能带2-3为两个简并轨道, 对应 $p_y$ 和 $p_z$ 电子, 能带 4 位于 $p_x$ 中的 $\alpha$ 电子。 这里能带 4 的能量应该比 2 和 3 要低, 但结果恰恰相反, 描述不合理;
  2. 在 $\beta$ 电子(spin component 2)中, 能带 1 是 $2s$ 的 $\beta$ 电子, 能带 2 为 $p_x$ 中的 $\beta$ 电子, 此时该 $\beta$ 电子占据了一个轨道, 而不是分布在三个 $p$ 轨道上, 结果是合理的。

  • 计算到这里, 已经是 VASP 官网对于单原子能量最为完整和准确的计算了。 虽然上面 $p_x$ 的 $\alpha$ 电子比 $p_y$ 和 $p_z$能量高, 至少电子的轨道电子占据情况是正确的, 且磁矩和简并的错误已经消除。 由此可见, VASP对于单原子的电子占据情况(至少对于 O 原子来说), 虽然取得了部分理想结果, 但还是有不足的地方。

  • 此外,体系的对称性可以通过 ISYM 来控制( 扩展阅读 ISYM ), ISYM = 0 的时候, 不考虑对称性进行计算.

  • 大师兄也尝试过了将 ISYM = 0 应用在 $8 \times 8 \times 8$ (test3)以及 $7.5 \times 8.0 \times 8.9$ (test4) 的格子里进行计算。 对于电子占据情况来说, 均得到了与前面 $7.5 \times 8.0 \times 8.9$ 一致的结果。 但是使用 $7.5 \times 8.0 \times 8.9$ 所得到的能量更低(test4 的能量比 test3 低)。

  • 使用 $7.5 \times 8.0 \times 8.9$ 格子时, ISYM = 0 (test4) 和采用默认值, 得到的能量结果是相同的, 因此, 改变晶胞的对称性在计算单原子能量的时候是必须的, 单独通过 ISYM 这一项不考虑对称性是远远不够的。

扩展阅读

  1. ISYM: http://cms.mpi.univie.ac.at/vasp/guide/node115.html
  2. 重现本节所有练习, 并自己认真分析结果, 查看电子占据情况, 体系的磁性, 以及能量信息。

总结

  1. 知道什么时候使用自旋极化计算,怎么进行自旋极化计算;
  2. 知道为什么自旋极化后计算量增加了;
  3. 知道在哪里查看体系的磁性信息;
  4. 知道为什么对称性降低后计算量增加了;
  5. 会分析单原子中电子轨道的占据情况;判断结果是否合理;
  6. 知道通过改变晶格大小调节体系的对称性从而消除简并对能量的影响;
  7. 知道ISYM这一个参数,在计算单原子能量中,ISYM的作用远远不及改变晶格大小;
  8. 知道vasp也存在自己局限性,不要完全相信。