Castep 与 vasp+phonopy 声子谱计算
chempeng / 2018-03-06
本文分别介绍 Materials Studio 中的 Castep 模块和 vasp+phonopy 使用 Finite difference method 和 DFPT 方法进行声子谱计算。
Castep
以 Ge 和 Fe 为例，介绍 Castep 模块的声子谱计算方法。
线性响应或密度泛函微扰理论（DFPT)
线性响应或密度泛函微扰理论是点阵动力学从头开始计算中最受欢迎的方法之一，尽管如此，这种方法的应用已经扩充到对振动属性的研究。线性响应提供了一种分析方法用于计算给定混乱的二级派生的整体能量。可以计算出许多属性，主要依赖于混乱的种类。在离子位置的混乱可以引起动力矩阵和声子；在磁场中引起NMR效应；在单位晶格矢量中产生弹性常数；在电场中引起非传导性效应等。
此部分以锗为例，使用 castep 进行线性响应声子谱计算。包含以下内容：
- 优化 Ge 单胞结构
- 计算声子谱及态密度
- 显示声子谱及态密度
- 显示热力学属性
1. 优化 Ge 单胞结构
导入锗的结构，在菜单栏中选择 File | Import。遵循下列路径 structures/metals/pure metals 选中 Ge.xsd。将其转换为原胞，使得计算速度加快。从菜单栏中选择Build | Symmetry | Primitive Cell。接下来优化锗的几何机构。
菜单栏中选择 Modules | CASTEP | Calculation。
- 在Setup标签上，把 Task 从 Energy 改为 Geometry Optimization，把 Functional 改为 LDA ,不选 Metal 。
- 选中 Electronic 标签，把 Energy cutoff 设置为 Ultra-fine，把 SCF tolerance 设置为 Ultra-fine，把 k-point set 设置为 Fine 以及把 Pseudo-potentials 设置为 Norm-conserving。
- 选中 Job Control 标签。选择你想要在其上运行工作的 Gateway location。把 Runtime optimization 设置为 Speed。
- 按下Run按钮开始运行。
工作提交后开始运行。它大概需要2分钟时间，这主要取决于你的电脑的速度。结果放在一个被称为 Ge CASTEP GeomOpt 的新文件夹中。
2. 计算声子谱及态密度
为了计算声子散射和声子的能态密度，在从 CASTEP Calculation 对话框的 Properties 标签选定适当的属性后，我们必须完成一个单点能量计算。
- 确定 Ge CASTEP GeomOpt 文件夹中的 Ge.xsd 文件示激活文档。
- 选中 CASTEP Calculation 对话框中的 Setup 标签，把 Task 设置为 Energy。
- 在 CASTEP Calculation 对话框中选定 Properties 标签。选择 Phonon，勾选 Both，按下 More… 按钮，勾选 Use Interpolation，并且 把Dispersion，Density of states 设置为 fine。
- 选中 Job Control 标签。选择你想要在其上运行工作的 Gateway location，提交作业。
在 Ge CASTEP GeomOpt 文件夹中创建了一个名为 Ge CASTEP Energy 的新文件。当能量计算完成后，两个新文件 Ge_PhononDisp.castep 和 Ge_PhononDOS.castep 放在此文件夹中。
3. 显示声子谱及态密度
声子散射曲线显示出声子能量沿着布里渊区高对称性方向如何依赖于 q 向量。此信息可以从单晶的中子散射实验中获得。只有为数不多的物质可以获得这样的信息，所以用来确定建模方法是否正确的理论偏差曲线对于论证从头开始计算方法的预测性能力是非常有用的。在一定情形下，它可能测量态密度而不是声子散射。而且与声子的态密度有直接关系的电子—声子交感作用可以通过隧道实验直接测量。所以能够从第一原理计算出声子的态密度是非常重要的。
Materials Studio 可以从任何 .phonon CASTEP 输出文件中产生声子散射图和态密度图。这些文件隐藏在 Project Explorer 里，但是一个 .phonon 文件会和每一个拥有 PhonDisp 或 PhonDOS 后缀的 .castep 文件一起产生。接下来使用先前的计算结果来创建声子散射图。
从 Materials Studio 的菜单栏中选择 Modules | CASTEP | Analysis。从属性列表中选择 Phonon dispersion。确定 Results file 选择框中显示的是 Ge_PhononDisp.castep。从 Units 下拉列表中选择 cm-1。从 Graph style 下拉列表中选择 Line。按下 View 按钮。
在结果文档中创建了一个新的图形文档 Ge Phonon Dispersion.xcd。它应和下图相似：
声子散射的实验图如下所示：
预测的频率可从Ge_PhononDisp.castep文件中得到。
在 Project Explorer 中双击 Ge_PhononDisp.castep 。按 CRTL+F 键，搜索 Vibrational Frequencies。结果文件中的部分内容如下所示：
==============================================================================
+ Vibrational Frequencies +
+ ----------------------- +
+ +
+ Performing frequency calculation at 40 wavevectors (q-pts) +
+ ========================================================================== +
+ +
+ -------------------------------------------------------------------------- +
+ q-pt= 1 ( 0.500000 0.250000 0.750000) 0.0487804878 +
+ -------------------------------------------------------------------------- +
+ Acoustic sum rule correction < 0.010115 cm-1 applied +
+ N Frequency irrep. +
+ (cm-1) +
+ +
+ 1 114.755657 a +
+ 2 114.755657 a +
+ 3 204.802542 b +
+ 4 204.802542 b +
+ 5 274.989890 a +
+ 6 274.989890 a +
+ .......................................................................... +
+ Character table from group theory analysis of eigenvectors +
+ Point Group = 32, Oh +
+ .......................................................................... +
+ Rep Mul | E 2 2 m -4 +
+ | -------------------- +
+ a 2 | 2 0 0 0 -1 +
+ b 1 | 2 0 0 0 1 +
每一个q点和每一个分支（纵向光波或声波(LO/LA),横向光波或声波(TO/TA)的频率以 cm-1 给出，同时也给出了 q 点在倒易空间中的位置。高对称性点 G, L 和 X 在倒易空间中的位置各自为(0 0 0), (0.5 0.5 0.5) 和 (0.5 0 0.5)。这些点和 q 点12, 6 以及19相对应。预测的频率（cm-1）和实验的频率（cm-1）如下：
总体来说，计算的精度是可以接受的。通过对更好的 SCF k 点格子的计算，我们可以获得更加另人满意的实验结果。
现在创建声子态密度图： 从 Materials Studio 菜单栏中选择 Modules | CASTEP | Analysis 。从属性列表中选择 Phonon density of states。确定 Results file 选择框中显示的是 Ge_PhononDOS.castep。把 DOS display 设置为 Full。按 View 按钮。创建了一个新的图形文档 Ge_PhonDOS Phonon DOS.xcd，如下图：
4. 显示热力学属性
在CASTEP中的声子计算可以用来评价近似准谐波晶体的焓，熵，自由能，格子的热容对于温度的依赖性。可以用这些结果和实验数据（如，热容的测量）相比较以预测不同的结构经过修正后的相稳定性和相转变。
所有与能量相关的属性均画在同一种曲线图中，并且零点能量的计算值也包括在内。热容被独自画在图表的右侧。现在使用声子计算的结果创建热力学属性图表。
从 Materials Studio 菜单栏中选择 Modules | CASTEP | Analysis 。从属性列表中选择 Thermodynamic properties。确定 Results file 选择框中显示的是 Ge_PhononDOS.castep。勾选上 Debye temperature 图，按下 View 按钮。在结果文件夹中创建了两个新的图形文档 Ge Thermodynamic Properties.xcd 和 Ge Debye Temperature.xcd。所示图形如下：
没有非谐性的实验结果表明在高温极限的 Debye 温度是395(3)K。模拟的 Debye 温度是396 K，很好的与实验值相符。考虑到所完成的计算的级别（given the level of calculation performed），这个结果还是可以接受的。使用更好的k点格子可以提高精确度。通过实验，推导出 Debye 温度在低温极限时的数值为374K，而 CASTEP 的预测值为460K。
我们也应该注意到类似的谐波近似是声子计算的基础，但是当温度高于 Debye 温度的三分之一时，类似的谐波近似是无效的。结果，在温度高于 Debye 温度时，非谐性影响出现，不能够仅仅依靠量的改变来解释。
尽管如此，总的来说，实验图表和由 CASTEP 产生的图表在质量上是非常相似的。当 Debye 最低温度达到255K时，在实验图表上会骤然降低25K; CASTEP 显示出温度骤降出现在相同的位置并且预测出 Debye 温度的最小值大约为270K。
有限位移法 （Finite difference method）
在点阵动力学计算中有两种主要的方法：密度泛函微扰理论(DFPT)和有限位移方法。第一种方法通常比较快速并且更加精确，但是它的执行是有问题的，并且受到一系列的限制。目前，在 CASTEP 中 DFPT 只能用于具有 fixed occupancies (绝缘体)、没有自旋极化并且只能对 norm-conserving 赝势使用。这样，对广泛的材料种类 (包括磁性材料和金属)，声子计算只能通过采用有限位移算法来执行。
此部分以铁为例，使用 castep 进行有限位移声子谱计算。包含以下内容：
- 优化 Fe 晶胞结构
- 计算声子谱及态密度
- 显示声子谱及态密度
- 显示热力学属性
1. 优化 Fe 晶胞结构
导入铁的结构，菜单栏选择 File | Import…，定位到 Structures/metals/pure metals 并选择 Fe.msi。通过将晶胞转化为原胞通常能减少计算时间。 从菜单栏选择 Select Build | Symmetry | Primitive Cell。从工具栏选择CASTEP工具，然后选择 Calculation，或者从菜单栏选择 Modules | CASTEP | Calculation。
在 Setup 选项卡上，将 Task 从 Energy 改为 Geometry Optimization，设置 Quality 为 Medium，设置 Functional 为 LDA。选中 Spin polarized 复选框，取消选择 Use formal spin as initial。设置 Initial spin 值为2。
LDA/CA-PZ 局域交换关联泛函被认为是可获得的最准确的描述之一，将 initial spin 值设为2是因为我们正在模拟铁磁性的 Fe 晶体。
几何优化的默认值不包含晶胞的优化。
- 点击 More… 按钮，在 CASTEP Geometry Optimization 对话框上，选中 Optimize cell，关闭对话框。
- 选择 Electronic 选项卡并设置 Pseudopotentials 值为 Ultrasoft。
- 选择 Job Control 选项卡，为 job 选择 Gateway 位置并设置 Runtime optimization 为 Memory。
- 击More…按钮，打开 CASTEP Job Control Options 对话框，在 Live updates 区域取消选择所有选项，关闭对话框。
- 点击 Run 按钮启动 job。
2. 计算声子谱及态密度
为了计算声子散射和声子态密度，必须执行单点能量计算，并为计算选择适当的性能。
- 确保 Fe CASTEP GeomOpt 目录中的 Fe.xsd 是活动文档。
- 在 CASTEP Calculation 对话框上选择 Setup 选项卡，并设置 Task 为 Energy 。
- 在 Properties 选项卡上，选择 Phonons 复选框，选择 Both 选项以选中态密度和散射。取消选择 Calculate LO-TO splitting，将 Method 选为 Finite displacement。
- 点击 More… 按钮，显示 CASTEP Phonon Properties Setup 对话框。确保 Method 为 Finite displacement。设置 Supercell defined by cutoff radius 值为3.6 Å。将 Dispersion 和 Density of states 的 Quality 值都设置为 Fine。关闭 CASTEP Phonon Properties Setup 对话框。
- 选择 Job Control 选项卡并为计算选择 Gateway。
- 点击 More… 按钮，打开 CASTEP Job Control Options 对话框并选中所有 Live update s选项，关闭对话框。
- 点击Run按钮。
注意：纵向光学横向光学(LO)劈裂不能对金属进行计算，因为它们在Γ-point是相同的。有限位移方案被设计用于金属和自旋极化系统(以及为那些具有高效超软势能的系统)。对铁磁性Fe这是理想的计算声子势能的方法。
Cutoff radius 的选择对有限位移运算是至关重要的参数。当使用较大的 cutoff radius 值时精度较高，因为这时考虑了更多的近邻。然而，随着该值的增加，计算时间增加的非常迅速。出于实际原因，在本教程中，对该参数选择了较小的值。声子频率的收敛作为 cutoff radius 的函数在执行有意义的实验计算时应该被研究。
Job 被提交并开始运行。一个叫做 Fe CASTEP Energy 的新目录被建立在 Fe CASTEP GeomOpt 目录中。当能量计算结束时，新的结果文件被放置在该目录中，包括 Fe_PhonDisp.castep 和 Fe_PhonDOS.castep。
3. 显示声子谱及态密度
生成声子谱
- 从 Materials Studio 菜单栏选择 Modules | CASTEP | Analysis，从性能列表中选择 Phonon dispersion，确保 Results file 选择器显示 Fe_PhononDisp.castep。
- 从 Units 下拉列表选择 cm-1，从 Graph style 下拉列表选择 Line。
- 点击View按钮。
生成声子态密度曲线。
- 从 Materials Studio 菜单栏选择 Modules | CASTEP | Analysis ，从性能列表中 选择 Phonon density of states，确保 Results file 选择器显示 Fe_PhonDOS.castep。
- 设置 DOS display 为 Full，点击 More… 按钮，打开 CASTEP Phonon DOS Analysis Options 对话框，选择 Interpolation 作为 Integration method，设置 Accuracy level 为 Fine。点击 OK 按钮。
- 在 CASTEP Analysis 对话框上点击 View 按钮。
4. 显示热力学属性
现在使用声子计算的结果创建热力学属性图表。
- 从 Materials Studio 菜单栏中选择 Modules | CASTEP | Analysis 。从属性列表中选择 Thermodynamic properties。确定 Results file 选择框中显示的是 Fe_PhononDOS.castep。
- 勾选上 Debye temperature 图，按下 View 按钮。
在结果文件夹中创建了两个新的图形文档 Fe Thermodynamic Properties.xcd 和 Fe Debye Temperature.xcd。
Vasp + Phonopy
- Phonopy 安装
- DFPT
- Finite difference method
1. Phonopy 安装
Phonopy 是一款由 Python 编译，基于第一性原理进行声子计算的开源软件。常规安装需要各种依赖与相应配置，繁琐复杂。安装 Anaconda 则可大大减少其复杂程度。
- Anaconda 安装
- Phonopy 安装
Anaconda 安装
下载 Anaconda 清华镜像 （国内下载速度快），本文以 Anaconda3-5.0.1-Linux-x86_64.sh 为例进行安装。
bash Anaconda3-5.0.1-Linux-x86_64.sh
按照提示进行安装，当提示是否添加添加环境变量时，输入 yes。或输入：
echo 'export PATH="~/anaconda3/bin:$PATH"' >> ~/.bashrc
source ~/.bashrc
终端输入 conda 查看是否安装成功，显示
[yourname@localhost ~]$ conda
usage: conda [-h] [-V] command ...
conda is a tool for managing and deploying applications, environments and packages.
..........
Phonopy 安装
至官网下载 Phonopy，本文以 phonopy-1.11.2.tar.gz 版本为例，提供百度云下载，解压
tar xzvf phonopy-1.11.2.tar.gz
进入安装目录，安装
python setup.py install
输入 phonopy 查看是否安装成功，显示
[yourname@localhost ~]$ phonopy
_
_ __ | |__ ___ _ __ ___ _ __ _ _
| '_ \| '_ \ / _ \| '_ \ / _ \ | '_ \| | | |
| |_) | | | | (_) | | | | (_) || |_) | |_| |
| .__/|_| |_|\___/|_| |_|\___(_) .__/ \__, |
|_| |_| |___/
1.11.2
Crystal structure file of POSCAR (default file name) could not be found.
___ _ __ _ __ ___ _ __
/ _ \ '__| '__/ _ \| '__|
| __/ | | | | (_) | |
\___|_| |_| \___/|_|
2. DFPT
此部分以官网算例 NaCl 介绍，做简单翻译补充。
1. 准备原胞，命名为 POSCAR-unitcell
Na Cl
1.00000000000000
5.6903014761756712 0.0000000000000000 0.0000000000000000
0.0000000000000000 5.6903014761756712 0.0000000000000000
0.0000000000000000 0.0000000000000000 5.6903014761756712
4 4
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 0.0000000000000000 0.5000000000000000
0.5000000000000000 0.5000000000000000 0.0000000000000000
0.5000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.5000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.5000000000000000
2. 扩胞
phonopy -d --dim="2 2 2" -c POSCAR-unitcell
3. vasp 计算力常数
以上一步生成的 SPOSCAR 为 POSCAR，准备 KPOINTS、POTCAR 和 INCAR，进行 vasp 计算。
mv SPOSCAR POSCAR
INCAR
PREC = Accurate
IBRION = 8
NSW = 1
ENCUT = 500
EDIFF = 1.0e-08
ISMEAR = 0
SIGMA = 0.01
LREAL = .FALSE.
ADDGRID= .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
4. 收集力常数
用 vasp 计算得到的 vasprun.xml 生成力常数文件 FORCE_CONSTANTS。
phonopy --fc vasprun.xml
5. 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
BAND = 0.0 0.0 0.0 0.5 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.5 0.5 0.5
FORCE_CONSTANTS = READ
$ phonopy --dim="2 2 2" -c POSCAR-unitcell band.conf #band plot
NOTE:
- 后处理时加参数
--factor=521.471可将频率单位转换为 cm-1。 bandplot --gnuplot band.yaml > band.dat可以把声子的结构转化成 dat 文件用 oringin 来画图。- 如果计算超内存，使用命令
ulimit -s unlimited。
3. Finite difference method
同以官网算例介绍
1. 准备原胞，扩胞
Si O
1.00000000000000
4.0444481195475825 0.0000000000000000 0.0000000000000000
0.0000000000000000 4.0444481195475825 0.0000000000000000
0.0000000000000000 0.0000000000000000 2.6211679887202837
2 4
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.5000000000000000 0.5000000000000000 0.5000000000000000
0.3030082337326380 0.3030082337326380 0.0000000000000000
0.6969917662673621 0.6969917662673621 0.0000000000000000
0.1969917662673621 0.8030082337326379 0.5000000000000000
0.8030082337326379 0.1969917662673621 0.5000000000000000
创建一系列有限位移的超胞
phonopy -d --dim='2 2 3'
产生 disp.yaml phonopy_disp.yaml POSCAR POSCAR-001 POSCAR-002 POSCAR-003 SPOSCAR
2. vasp 计算力常数
准备输入文件KPOINTS、POTCAR、POSCAR和INCAR，使用脚本批量计算，generation.sh (根据所产生 POSCAR 数量修改）
#!/bin/bash
for i in 1 2 3; do
mkdir $i
cp POSCAR-00$i $i/POSCAR
cp INCAR KPOINTS POTCAR $i/
done
INCAR
PREC = Accurate
IBRION = -1
ENCUT = 500
EDIFF = 1.0e-08
ISMEAR = 0
SIGMA = 0.01
LREAL = .FALSE.
LWAVE = .FALSE.
LCHARG = .FALSE.
3. 收集力常数
phonopy -f 1/vasprun.xml 2/vasprun.xml 3/vasprun.xml
4. phonopy 后处理
声子态密度，热性质（自由能，熵，等压热容）的计算需要准备 mesh.conf 文件，投影态密度的计算准备 pdos.conf 文件，声子谱计算准备 band.conf 文件
# mesh.conf
ATOM_NAME = Si O
DIM = 2 2 3
MP = 8 8 8
$ phonopy -p mesh.conf #dos plot
$ phonopy -t mesh.conf #thermal properties print
$ phonopy -t -p mesh.conf #thermal properties print and plot
# pdos.conf
ATOM_NAME = Si O
DIM = 2 2 3
MP = 8 8 8
PDOS = 1 2, 3 4 5 6
$ phonopy -p pdos.conf #pdos plot
# band.conf
ATOM_NAME = Si O
DIM = 2 2 3
BAND = 0.5 0.5 0.5 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5 0.0
$ phonopy -p band.conf #band plot
参考资料：