有限位移法计算硅的声子谱与声子态密度
计算声子谱和态密度，目前有两种常用的方法：有限位移法（Finite Displacement） 和 DFPT（密度泛函微扰理论）。
- 有限位移法：简单粗暴。把超胞里的某一个原子人为推开一点点，用VASP计算其他所有原子受到的力，从而推导出力常数。必须用超胞。
- DFPT（密度泛函微扰理论）：更高级的量子力学方法，直接在倒空间算。理论上是原胞即可，但为了兼容Phonopy，有时也需要结合超胞。计算量较大，但精度高。
如果要使用VASP计算整个声子谱或者声子态密度，则需要借助辅助软件，如Phonopy、PHON、Phonon等。这里以Phonopy为例，计算硅的声子谱和声子态密度。
上述软件都是用超元胞方法来计算声子色散关系，但首先第一步结构优化是一样的。
第一步 结构优化
同使用VASP计算能带一样，需要进行结构优化（OPT）。
需要额外注意的几点：
- 要进行高精度的结构优化，即参数
ISIF = 3。 - 由于声子频率对原子受力的精确度非常敏感，后续所有计算都要保持相同的截断能（ENCUT）和精度（EDIFF），确保力常数计算的准确性。
文件参数示例
以下是我进行OPT时所用的文件参数：
1. INCAR
SYSTEM = FCC Si
PREC = A
ENCUT = 350
ISMEAR = 0; SIGMA = 0.05
ADDGRID = .T.
LCHARG = .F.
LWAVE = .F.
EDIFF = 1e-7
NSW = 50
ISIF = 3
EDIFFG = 0.02
IBRION = 2
POTIM = 0.2下面带注释的版本，我不小心用了Tab键，会引入一些看不见的符号，运行VASP时可能会报错，所以复制的时候请复制上面的代码。
# 常规参数
SYSTEM = FCC Si
PREC = A # 即accuracy，测试下来我用Normal画出来的也不赖
ENCUT = 350 # 可根据你的体系适当调
ISMEAR = 0; SIGMA = 0.05
ADDGRID = .T. # 待会马上在下面讲这个参数
LCHARG = .F. # 不输出电荷密度，节约空间
LWAVE = .F. # 不输出波函数，节约空间
# KSPACING = 0.15 # 可用这个自动生成K点，也可手动写KPOINTS，这里为了加深理解不使用这个参数，手写KPOINTS
# 结构优化参数
# 电子步
EDIFF = 1e-7 # 电子步精度
# 离子步
NSW = 50 # 离子步步数，因为FCC Si很容易优化，所以这里的步数偏少
ISIF = 3
EDIFFG = 0.02 # 离子步精度
IBRION = 2 # 共轭梯度算法，进行结构优化
POTIM = 0.2 # 稍微减小IBRION步长，让优化更稳其中的参数 ADDGRID
这是一个用于提升计算精度的参数。简单来说，它会在VASP计算中启用一个额外的、更密集的网格，专门用来更精确地处理原子核附近的电荷密度。这在计算声子（晶格振动）时很重要，因为声子频率对原子受力的精确度非常敏感。
为什么需要额外的网格？
要理解这个参数，可以先了解VASP是如何处理电荷密度的：
- 常规网格：VASP会使用一个基础的网格来“描绘”整个晶胞内的电子波函数和电荷密度。这个网格的疏密程度由
PREC（精度）和NGX/NGY/NGZ等参数控制。 - “补丁”网格（PAW方法）：由于原子核附近的电子密度变化极其剧烈，为了精确描述这个区域，VASP（采用PAW方法时）会在原子核周围加上一个“补丁”式的球内网格来进行计算。这个“补丁”区域的电荷被称为增广电荷（augmentation charge）。
- ADDGRID的作用：
ADDGRID = .T.会在这个“补丁”区域启用一个更密集的辅助网格来更精确地计算增广电荷。根据VASP手册，这个额外网格的格点数是常规网格的8倍。
对声子计算的意义
在声子计算中，每个原子的受力是由电子密度分布决定的。如果原子核附近的电荷分布计算不精确，力就会出现微小误差，进而导致声子频率出现误差。
- 对于硅（Si）这类半导体，开启
ADDGRID通常能显著提高力的计算精度，有助于获得更平滑、更准确的声子色散曲线，并减少虚假的“虚频”。 - VASP官方也确认，
ADDGRID=.TRUE.有时可以进一步提高力的计算精度，尤其是在需要高精度力（如声子计算）的场景下。
2. KPOINTS
K-POINTS
0
Gamma
10 10 10
0 0 03. POSCAR
FCC Si
5.43
0.5 0.5 0.0
0.0 0.5 0.5
0.5 0.0 0.5
Si
2
c
0 0 0
1.35 1.35 1.354. POTCAR
还需要Si的POTCAR。
有了上面这四个文件，直接丢给VASP进行OPT即可。
第二步 生成超胞与位移文件
重命名
首先，将前面OPT得到的 CONTCAR 文件重命名为 POSCAR-unitcell（名字可以随意，接下来马上会讲）。
在结构优化后，需要使用Phonopy的命令把初基原胞扩展成超胞。超胞在空间中重复原来原胞的结构，形成一个更大的‘原胞’。
对于Si，以2×2×2超胞为例已经足够，可以使用以下Linux命令：
phonopy -d --dim="2 2 2" -c POSCAR-unitcell（注意：最新版输入该指令可能会报错，不过终端会提示你该怎么做，在接下来“注意”部分会讲到为什么）
指令参数的意义
1. -d 这是 --dim 的简写，全称是 --dimension（维度/倍数）。 作用：启用Phonopy生成超胞并同步生成位移文件。
2. --dim='2 2 2' 这里双引号和单引号均可。 含义：指定超胞大小，这里 2 2 2 表示沿着晶格的 a、b、c 三个基矢方向，各扩展 2 倍。 对于本次例子的FCC Si，原胞有2个硅原子，扩完之后超胞里就有 \(2 \times 2 \times 2 \times 2 = 16\) 个原子。
3. -c POSCAR-unitcell 含义：-c 是 --cell 的缩写, 表示读取目录中名为 POSCAR-unitcell 的文件作为初始原胞。 所以前面说到，对 CONTCAR 重命名是可以任意的，只要这里跟着变就可以了。若这里不填，phonopy就会默认搜寻目录中名为 POSCAR 的文件。
输出的文件
POSCAR-unitcell：备份的原胞。SPOSCAR：没有位移的完美超胞。POSCAR-001,POSCAR-002…：根据晶体对称性，Phonopy自动判断需要移动哪些不等价的原子，每一个文件代表一种位移构型（比如硅只需移一次，生成一个POSCAR-001；如果对称性低，可能会生成几十个）。
注意
前面提到，输入 1. 的命令可能会报错，这是因为 Phonopy 是 v4 或更高版本，其命令结构发生了重大变化（本人深受其害）。
在 v4 中，Phonopy 将命令拆分为两个部分：负责“准备”的 phonopy-init 和负责“计算”的 phonopy。--dim 被移动到 phonopy-init 模块。
所以你要把命令换成：
phonopy-init -d --dim="2 2 2" -c POSCAR-unitcell除此之外，还有以下这些变化：
| 操作类型 | v3 命令 (旧) | v4 命令 (新) |
|---|---|---|
| 扩胞/生成位移 | phonopy -d --dim ... |
phonopy-init -d --dim ... |
| 提取力常数 | phonopy -f vasprun.xml-* |
phonopy-init -f vasprun.xml-* |
| DFPT 计算 | phonopy --fc vasprun.xml |
phonopy-init --fc vasprun.xml |
| 绘制能带图 | phonopy -p band.conf |
phonopy -p band.conf (不变) |
第三步 自洽计算
本步的目的是优化电子的电荷密度分布，计算的过程和前面的SCF基本一致。
需要额外注意的点
- 若你前面第二步生成了多个位移文件
POSCAR-XXX，你就要进行多次SCF计算。最好为每个POSCAR-00x建立一个独立的文件夹（例如disp001,disp002…）。 - 别忘了保持相同的截断能（ENCUT）和精度（EDIFF）。
输入文件参数
1. INCAR
# 常规参数
SYSTEM = FCC Si
PREC = A # 即accuracy，测试下来我用Normal画出来的也不赖
ENCUT = 350 # 可根据你的体系适当调
ISMEAR = 0; SIGMA = 0.05
ADDGRID = .T. # 提升受力精度
LCHARG = .F. # 不输出电荷密度，节约空间
LWAVE = .F. # 不输出波函数，节约空间
# KSPACING = 0.15 # 可用这个自动生成K点，也可手动写KPOINTS，这里为了加深理解不使用这个参数，手写KPOINTS
# 结构优化参数
# 电子步
EDIFF = 1e-7 # 电子步精度
# 离子步
NSW = 0 # 离子步步数，因为FCC Si很容易优化，所以这里的步数偏少
# ISIF = 3
# EDIFFG = 0.02 # 离子步精度
IBRION = -1 # 固定离子位置，计算单点能
POTIM = 0.2 # 稍微减小IBRION步长，让优化更稳同样的，这里给出没有注释的干净版本：
SYSTEM = FCC Si
PREC = A
ENCUT = 350
ISMEAR = 0; SIGMA = 0.05
ADDGRID = .T.
LCHARG = .F.
LWAVE = .F.
EDIFF = 1e-7
NSW = 0
IBRION = -12. KPOINTS
K-POINTS
0
Gamma
5 5 5
0 0 0注意这里的K点，在各个方向减半。这是因为扩胞后，超胞的晶格常数是原来原胞的两倍，即正空间中各个基矢的长度变为原来的两倍，那么倒空间中的基矢的长度变为原来的 \(\frac{1}{2}\)，为了保持相同的采样密度（k点之间的间距保持不变），这里必须这样做。
可以参考前面章节
输入文件的参数 -> POSCAR -> 对POSCAR的理解 -> 扩胞下的情况
前面的INCAR里提到了参数 KSPACING = 0.15，在不写KPOINTS的情况下，这个参数直接限定k点之间的间距，自动生成 KPOINTS。所以使用这个就不用自己更改 KPOINTS 文件了，会方便得多。
对于声子谱计算，推荐的做法是：
- 统一使用
KSPACING：在结构优化（原胞）和声子计算（超胞）的所有VASP步骤中，都在INCAR里设置相同的KSPACING值，并确保计算目录中不存在KPOINTS文件。 - 做好收敛测试：虽然
KSPACING很方便，但0.15这个值是否适合你的体系，仍需要通过测试来判断。一个可靠的方法是，用KSPACING生成K点后，在OUTCAR里查看VASP实际使用的K点网格，然后手动在KPOINTS文件里设置一个更密的网格（例如，把生成的4 4 4手动改为6 6 6），对比两次计算的总能或力是否收敛。如果变化很小，说明原来的KSPACING已经足够。
3. POSCAR 将第二步生成的位移文件 POSCAR-XXX，改名为 POSCAR。
4. POTCAR 同上。
第四步 提取力常数
这一步是把VASP算出来的力，转化为Phonopy认识的矩阵。
对于有限位移法： 在包含所有 disp00x 文件夹的总目录下运行（若有多个位移文件，就这样子）：
phonopy -f disp001/vasprun.xml disp002/vasprun.xml ... 新版本使用（这里只有一个位移文件，只生成了一个 vasprun.xml）：
phonopy-init -f vasprun.xml这会生成 FORCE_SETS 文件。
计算逻辑：
- VASP本身不知道什么是“力常数”。它只是在你把原子人为推开（位移）后，算出每个原子此时此刻受到的力（Forces）。
-f告诉 Phonopy：去读 VASP 算出来的这些受力。- Phonopy拿到受力后，会结合扩胞时生成的位移说明书（
phonopy_disp.yaml），通过数学公式计算（组装）出力常数。 - 输出文件：生成
FORCE_SETS文件。
第五步 准备配置文件与绘图
这里需要准备两个文件 band.conf 和 mesh.conf，来分别计算声子色散关系和声子态密度。
声子能带图配置（band.conf）
ATOM_NAME = Si # 你的元素名
DIM = 2 2 2 # 必须和你扩胞的倍数一模一样
BAND_LABELS = G X W K G L
FC_SYMMETRY = .TRUE. # 开启对称性加速Phonopy会自动找同一目录下的 FORCE_SETS。
声子能带图绘制
准备好 band.conf 后，确保目录下有 band.conf, FORCE_SETS，执行：
若不画图，只生成 band.yaml 文件，使用：
phonopy --config band.conf --band "0 0 0 0.5 0 0.5 0.5 0.25 0.75 0.375 0.375 0.75 0 0 0 0.5 0.5 0.5"指令的意义
--config：读取目录中的band.conf文件，进行配置。--band：表示开启绘制能带模式，其后面的数字表示，按band.conf中BAND_LABELS = G X W K G L顺序的k点坐标。
若画图，则在指令后加上 -p 指令，表示启用绘图，即：
phonopy --config band.conf --band "0 0 0 0.5 0 0.5 0.5 0.25 0.75 0.375 0.375 0.75 0 0 0 0.5 0.5 0.5" -p声子态密度配置（mesh.conf）
ATOM_NAME = Si
DIM = 2 2 2
MP = 30 30 30 # 在布里渊区划分的网格密度，越大越平滑其中 MP 表示声子态密度计算所需的 k 点。
同样的，若不画图，只生成 mesh.yaml 文件：
phonopy --config mesh.conf --mesh "30 30 30"指令的意义
--mesh：表示开启绘制网格样本模式，后面的数字表示空间中k点的分割数目。
投影态密度（PDOS）
若你想看投影态密度（PDOS），则在声子态密度配置（mesh.conf）中加上 PDOS 参数：
ATOM_NAME = Si
DIM = 2 2 2
MP = 40 40 40
PDOS = 1 2 # 表示分别计算第1个和第2个原子的投影态密度使用其它软件绘图
上一节我们知道输出文件中有：band.yaml 和 mesh.yaml。通常其它绘图软件不会识别这两个文件，因此我们要把它们转为 .dat 格式：
phonopy-bandplot --gnuplot band.yaml > band.datband.dat 的第一列是 X 轴（高对称点路径的距离），后面几列是 Y 轴（不同声子支的频率）。注意：不同高对称点之间会用空行隔开。
phonopy-dosplot --gnuplot mesh.yaml > mesh.datmesh.dat 第一列是频率（Y轴），第二列是态密度（X轴）。
DFPT法计算硅的声子谱与声子态密度
DFPT（密度泛函微扰理论）：更高级的量子力学方法，直接在倒空间算。理论上使用原胞即可，但为了兼容Phonopy，有时也需要结合超胞。计算量较大，但精度高。
前两步与有限位移法相同，这里从第三步开始讲起。
第三步 自洽计算
输入文件参数：
与有限位移法相比，主要不同在于 IBRION = -8 以及 NSW = 1。
# 常规参数
SYSTEM = FCC Si
PREC = A # 即accuracy，测试下来我用Normal画出来的也不赖
ENCUT = 350 # 可根据你的体系适当调
ISMEAR = 0; SIGMA = 0.05
ADDGRID = .T.
LCHARG = .F. # 不输出电荷密度，节约空间
LWAVE = .F. # 不输出波函数，节约空间
# KSPACING = 0.15 # 可用这个自动生成K点，也可手动写KPOINTS，这里为了加深理解不使用这个参数，手写KPOINTS
# 结构优化参数
# 电子步
EDIFF = 1e-7 # 电子步精度
# 离子步
NSW = 1 # 只有NSW不为0，VASP才会去检测IBRION的值
# ISIF = 3
# EDIFFG = 0.02 # 离子步精度
IBRION = -8 # 使用DFPT方法计算声子谱
SYMPREC=1E-6 # 如果不加，会遇到 "non-integer element in rotation matrix" 报错由此得到 vasprun.xml 文件。
第四步 提取力常数
这里在含 vasprun.xml 文件的目录使用：
phonopy --fc vasprun.xml这会生成 FORCE_CONSTANTS 文件。 (注：如果在v4版本中报错，请替换为 phonopy-init --fc vasprun.xml)
第五步 准备配置文件与绘图
同有限位移法，需要准备 band.conf 与 mesh.conf。
区别在于，要加上一个指令：FORCE_CONSTANTS = READ：
band.conf:
ATOM_NAME = Si # 你的元素名
DIM = 2 2 2 # 必须和你扩胞的倍数一模一样
BAND_LABELS = G X W K G L
FORCE_CONSTANTS = READ # 告诉Phonopy直接读取上一步生成的文件
FC_SYMMETRY = .TRUE. # 开启对称性加速mesh.conf:
ATOM_NAME = Si
DIM = 2 2 2
MP = 30 30 30
FORCE_CONSTANTS = READ
发表回复