以氢原子和氧原子为例的VASP单原子能量的计算
type
status
date
slug
summary
tags
category
icon
password
在学习VASP计算的过程中，单原子的计算是不错的起点。
氢原子由一个带有一个正电荷单位的质子组成的原子核加带有一个单位负电荷的电子组成。
以H原子为例，H原子的基态波函数可以用薛定谔方程精确求解得到E=-13.6eV，如果是单原子也可以假定E=0，那么VASP的计算结果如何呢？
实际上VASP计算结果既不是-13.6eV也不是0，这是因为在VASP的POTCAR文件中使用的原子的量能都是参考能量，也就是相对于一个参考的基准上的能量，计算结果也是一个相对值。这个参考基准并不是H原子的基态能量-13.6eV，因为参考基准的计算并不允许有自旋极化，也就是赝势都不是自旋极化的，但是薛定谔方程求解的H原子的基态波函数是自旋极化的；另一个原因是因为我们算体系的能量是在DFT的框架下算的，而-13.6是严格的求解薛定谔方程得到的，在DFT的框架下电子的这个自能项跟交换项实际上是不能完全抵消的，这也导致了参考态的能量跟氢原子真实的基态-13.6有一定的差距。实际上有参考基准下的VASP计算的结果显示为-0.xxeV，也是一个很小的数值，这个数值应该是在参考态的基准上加上这个计算结果。
VASP四个核心必要的输入文件：
INCAR/POSCAR/KPOINTS/POTCAR，必须全部大写（linux环境下文件名是区分大小写的）
POTCAR文件
1、第1行表示赝势用的是PAW_PBE
2、第2行表示1个价电子。
3、6行EATOM = 12.4884 eV, 0.9179 Ry，指的就是参考基准的能量，这里的12.4884eV是绝对值，实际上是负数的值。
VASP为元素周期表中的几乎所有元素提供了相应的势函数文件。当计算的体系中含有多种元素时，我们需要根据元素在POSCAR文件中出现的顺序，将各元素的势函数文件按顺序拼接起来，从而生成对应的POTCAR文件。
在本节的例子中，考虑到精度和效率的综合因素，我们选择使用PAW_PBE势文件。
INCAR 文件
- VASP的核心控制文件，主控文件；
- 告诉VASP需要干什么和怎么干；
- 格式使用tag=value的形式，tag是关键词，value是赋值；
- 不同的tag之间没有先后顺序之分，重复出现的使用第一个。
1、变量或者说标签的赋值格式如下，SYSTEM = XXXX
SYSTEM是标签变量名，后面加空格，再加=，再加空格，再加自定义字符串，这里的赋值仅作为我们自己使用，对计算没有任何作用，自定义知道这个计算是计算什么的即可，英文无特殊字符。
2、ISMEAR决定了如何设置每个轨道的局部占位率 $$f_{nk}$$
ISMEAR和SIGMA决定做布里渊区积分时如何计算分布函数，默认值1和0.2；
ISMEAR = 0 表示使用gaussian展宽，一般用于半导体和绝缘体，同时孤立体系的SIGMA一般设置较小值，如0.05 eV；
ISMEAR = 1或2 表示采用一阶或者二阶的Methfessel-Paxton方法，一般用于金属体系（避免用于半导体和绝缘体），同时可设置较大的SIGMA值，如0.2eV，以保证VASP计算的熵值小于1.0meV/atom；
ISMEAR = -5 表示四面体积分方法，一般适用于高精度的总能量和态密度计算，且不需设置SIGMA值。但四面体积分法对k点非常少（如只有1或2个)的情况以及一维体系不适用。
综上分析，采用默认值并非计算小分子的优选策略，推荐使用ISMEAR=0、SIGMA=0.05。
3、ISPIN = 2 表示打开自旋极化计算，因为只有一个H原子，不是H2，是有自旋极化的，如果是H2则两个原子磁矩抵消，那么就要赋值ISPIN=1。
这里关于该变量设置的VASP的wiki详细介绍：
https://www.vasp.at/wiki/index.php/ISMEAR
POSCAR文件
- 用于描述待计算对象的结构信息，包括模拟盒子（原胞或超胞）的尺寸和形状，以及原子的坐标；
- 在进行分子动力学计算时，POSCAR文件还可能包含初始速度等信息。
- 值得注意的是，VASP要求所有的计算对象必须为周期性体系。
因此，对于本节的小分子气体体系，我们可以将这些小分子置于一个较大(例如，尺寸为10Å×10Å×10Å)的超胞中，以确保镜像原子间的相互作用可以忽略不计。
1、如果在INCAR文件中没有设置的标签均采用系统默认值。
2、第1行是注释行，完全自定义，不要有中文和特殊字符，自己看的懂做什么计算即可。
3、第2行为元胞大小的缩放因子，通常设置为1。
4、第3至5行为元胞的三个基矢，它们乘第二行的缩放因子得到元胞真正的大小。晶格参数xyz三组坐标。晶格参数设置得足够大，这样相邻细胞中的原子之间就不会发生明显的相互作用。
5、第6行和第7行分别为元胞中元素类型和对应的原子个数，若含有多种元素，不同元素及其对应原子数间需用空格分开。
6、第8行为原子坐标的类型，可选择“Direct”或“Cartesian”。VASP其实只识别该行的第一个字母，若为“D"或“d"即是分数坐标，为“C"或“c"则是直角坐标。cart表示此处使用笛卡尔坐标系。
7、从第9行开始，每一行对应一个原子的坐标，0 0 0 表示在原点位置，笛卡尔坐标系的原点。其元素类型的顺序需要和第六行元素类型顺序一致。
8、POSCAR文件可以导入OVITO软件进行三维可视化。
KPOINTS文件
- 用于描述在倒易空间布里渊区k点的分布。
- 在VASP中k点有三种生成模式，分别为自动产生k点、手动输人所有k点坐标和线性模式。
对于本节的小分子算例，无须考虑Bloch定理的影响，只需自动产生一个k点，即 $$γ$$点。
一个k点就足以描述孤立的原子或分子，因为波矢量k连接相邻的晶格点，即从一个晶格点平移到相邻的晶胞。因此，包含更多的k点可以更准确地描述相邻晶胞中原子之间的相互作用。但如果我们在POSCAR文件中确实已将盒子设置得足够大以描述孤立的原子，那么这种相互作用应该是零。
1、第1行为注释；
2、第2行o，表示VASP自动生成k点；
3、第3行为产生k点的方式，有Gamma和Monkhorst-Pack两种选择，其中Gamma方式确保k点一定包含倒空间原点（0，0，0)，而选择Monkhorst-Pack方式则k点不一定包含原点，Monkhorst-Pack方式不适用于六边形结构；
4、第4行为由空格分开的三个整数，表示沿三个倒易基矢方向的k点的数目，切记：孤立体系必须取1 1 1；
5、第5行表示对所有k点进行平移，取0即表示不平移。
检查文件：
文件在计算服务器上准备完毕后，检查必须的文件和提交脚本：
使用cat和head检验核心文件和提交脚本：
注意sub.vasp中的#不是注释。来自管理员的解答。
提交：
氧原子的能量计算跟H很多一致的地方，重点关注不一样的地方。
POTCAR文件
1、第1行表示赝势用的是PAW_PBE
2、第2行表示vasp计算考虑6个价电子。
3、第4行：氧的价电子构型（vasp赝势的参考态？）
4、第12行：ZVAL = 6表示有6个价电子
INCAR 文件
这里与H原子计算不同的地方是考虑ISPIN 取值=1还是2？我们先假定=1，也就是没有自旋极化，因为考虑到s2p4电子分布，自旋相互抵消，磁矩=0，那么可以令ISPIN=1，作为非自旋极化计算，此时我们看一下计算的结果：
OUTCAR在ISPIN=1的情况下：
因为是孤立系统，所以此处的能带可以等同于能级和原子轨道。
能级1的能量是-23.8438，为S轨道，occupation共占用2个电子，是整数没问题。
能级2,3,4能量是-8.9040，为简并轨道P轨道，但是占用的是1.333，这与实际氧原子的基态电子排布不符。并且按照此计算的氧原子的磁矩与实际氧原子的磁矩也不符。
现在我们把ISPIN=2重新计算，注意在相同文件夹中重新计算vasp需要删除WAVECAR文件。
重新计算的OUTCAR文件：
新的计算结果中有spin component 1和2，可以理解为自旋方向向上和向下的两个分量：
spin component 1中有4个电子，可设定为自旋向上。
spin component 2中有2个电子，对应为自旋向下，上下差2个电子，计算的磁矩=2 $$μ_B$$即两个玻尔磁子，此时符合实际情况。但本次计算中占据数依然有非整数的问题。
问题出在一个vasp计算的细节问题上，当计算一个原子的基态波函数的时候需要打破对称性，因为氧原子的基态波函数是不满足立方晶格对应的旋转对称性的，我们在书写POSCAR文件时，通常是强加了对称性上去，比如下面的8*8*8的晶胞。那么计算结果得到的基态能量实际上是高于氧原子的基态能量的，那么正确的POSCAR应该打破对称性。
POSCAR文件
采用非对称的晶胞打破对称性后重新计算的总能量应该小于有对称性的基态能量计算结果。
新的计算结果OUTCAR如下：
此时可以得到如下：电子数：自旋向上4个，2s低能轨道1个，2p高能轨道3个；自旋向下的2个，2s低能轨道1个，2p高能轨道1个。所以出现右侧的示意图，注意示意图表示了简并的轨道，而实际轨道的能量有微小的差距。
KPOINTS文件
与H原子的计算完全一致。
Loading...