VASP-phonopy计算声子谱

主要参考：
Jun’s Blog
Bilibili-VASP+phonopy计算声子-中山大学王伟良

背景

在学校的高性能计算平台上计算。

phonopy下载与安装

# 注意conda环境，可以参考我上一篇内容安装anaconda，记得换源
pip install phonopy

声子谱计算

(以下个人见解，可能很多错漏)

声子简单来说就是描述晶格振动模式的量子化的准粒子。
它是没有实体对应的，类似于你将多个不同周期的正弦函数每一个都看成一个“实体”，把它们叫做“正弦子”之类的东西。
类似的概念也有“能级(能带)”(能带是离域化的能级)。
能级是电子的振动模式，声子是晶格的振动模式；能级谱/能带是电子振动模式的可能，声子谱是晶格振动模式的可能。

实操

声子谱计算有两种方法：

1.有限位移法

对原胞进行高精度的优化

主要影响参数：EDIFF,EDIFFG,ISIF
一个INCAR文件的参考：

PREC = Accurate
ENCUT = 500
IBRION = 2  # 结构优化算法
ISIF = 3  # 3是既优化原子坐标也优化晶格
NSW = 100  # 优化步数
NELMIN = 5  # 最小电子自洽步数
EDIFF = 1.0e-08  # 自洽能量收敛判据
EDIFFG = -1.0e-08  # 结构优化能量收敛判据
IALGO = 38  # 算法选择
ISMEAR = 0  # Gaussian smearing
SIGMA = 0.1  # smearing 宽度
LREAL = .FALSE # 以下三个都是控制输出文件的
LWAVE = .FALSE
LCHARG = .FALSE
NCORE = 4  # 并行核数

在OUTCAR中看到total force优化到1.0e-04就可以尝试下一步了。
如果优化了很多步，total force 都不减小只上下跳动，那么更改ISIF的参数，改为ISIF = 1。

扩胞

phonopy -d --dim='3 3 3' -c CONTCAR

–dim='3 3 3’相当于三个方向上都有三个晶胞，三维空间中有27个晶胞。
一般而言控制在100个原子左右较好(老师傅说的)
扩胞后，会出现 POSCAR-001~POSCAR-011 这样文件名的文件，将它们分别放进子文件夹中。
这些东西是 N 个位移后的 POSCAR ，N的数量取决于结构的对称性，对称性越好 N 越少。
自洽计算各个位移之后的 POSCAR ，可以得到二阶力常数，同时也得到声子谱、声子态密度、等容热容等热力学性质的原始数据。
可以参考使用shell的循环命令：

for i in {001..011};do mkdir $i && cp POSCAR-$i $i\POSCAR;done

# 一个简单的循环语句
# {001..011}类似一个列表，类似于python的range(1,12)
# 注意{001..011}大括号内没有空格
# do后面是每次循环要干的事包括 mkdir 与 cp 两句命令
# 中间的&&是连接两个命令的，只有前面mkdir执行成功才会执行cp命令

for i in {001..011};do cp INCAR-Phono $i/INCAR; cp POTCAR $i/POTCAR; cp SubLa.sh $i/SubLa.sh; cp KPOINTS-Phono $i/KPOINTS ;done
# 批量准备文件

for i in {002..011};do cd $i; sbatch SubLa.sh; cd ..;echo $i is done;done
# 批量提交任务,各位的服务器并不一定是sbatch的，在上面sbatch SubLa.sh 改成你自己的提交任务的脚本文件。

更改INCAR文件：
INCAR之中有以下几点需要注意：

1. 磁性的数量，扩胞之后原子数量更改，磁性因此要跟着更改。
2. IBRION项改为IBRION = -1
3. NSW 项改为 NSW = 0 (单点计算)
4. 其他的更改项可以等提交出错了之后再更改……比如我设置了NBANDS项，提交出错之后，把这一项删了让VASP自己设置就好了。
5. 应该尝试减少KPOINTS网格密度(将会大幅度减小计算时间以及出错的概率，同时也能减少 vasprun.xml 的大小从而降低 phonopy 的运算时间)
   将INCAR、POTCAR、KPOINTS、提交脚本都放进/001里，然后尝试提交一下。
   如果尝试可以了，就把/001中的 INCAR 改名为 INCAR-Phono、KPOINTS 改名为 KPOINTS-Phono复制出来。(不改名也可以，把我上面第二个循环里面的相应的文件名改成你自己的)

提交任务

在各个文件夹中提交一次vasp任务，计算完成之后会生成一个 vasprun.xml 文件。

ls */vasprun.xml

# 该命令可以列出所有子文件夹中的vasprun.xml
# 如果没有的话可能是你处于子文件夹中

收集结果

phonopy -f 001/vasprun.xml 002/vasp.xml

运行完成之后即生成了一个 FORCE_SETS 文件
准备phonopy的输入文件，文件名可以随便取，比如叫band.conf

ATOM_NAME = Na Cl
DIM = 2 2 2
#PRIMITIVE_AXIS = 0.0 0.5 0.5   0.5 0.0 0.5   0.5 0.5 0.0
FORCE_SETS = READ
EIGENVECTORS = .TRUE.
BAND_CONNECTION = .TRUE.
BAND_LABELS = Gamma M K
BAND = 0.0 0.0 0.0  0.5 0.0 0.0  0.5 0.5 0.0

EIGENVECTORS = .TRUE.输出文件包含原子振动信息
BAND = 0.0 0.0 0.0 0.5 0.0 0.0 0.5 0.5 0.0是画声子谱的时候，是从 Gamma(0.0 0.0 0.0) 画到 M 再画到 K

phonopy -c POSCAR-unitcell -p -s band.conf

运行完成之后会生成一个 band.pdf，即可看到图。

结果与计算优化

先看声子谱

原则上声子谱不应该有负频率。
如果原点有负，则需要提高原胞结构优化的精度。
如果其他地方有负频率，需要扩大超胞。

原子振动信息

在 band.yaml 中可以看到原子的振动信息。

2 密度泛函微扰理论

高精度优化原胞

同2

构建超胞

与上面的命令相同。
在使用上面的命令之后，除了POSCAR-001等等，还会出现一个SPOSCAR文件。

使用VASP进行微扰计算

INCAR 参数调整:

1. NSW = 1
2. IBRION = 8,微扰计算
3. NPAR, NCORE, KPAR 关闭以免报错
   KPOINTS 调整 k 点密度

计算力常数

phonopy --fc calc-force/vasprun.xml

准备后处理文件band_conf

phonopy -p -s band_conf

TODO

Phonopy准备的位移文件的依据（为什么我的计算会有11个需要计算的东西？）（是布里渊区吗？）
使用晶胞与原胞看Phonopy的位移如何