1、前言
在第一性原理计算过程中考察体系稳定性是经常遇到的一个问题，也是大部分审稿人容易提问的地方；通常，使用声子谱研究体系的动力学稳定性，使用分子动力学研究体系的热稳定性；另外，可以借助于一些其它方法进一步说明体系的稳定性，例如研究钙钛矿体系时通过离子半径判断稳定性，研究异质结时通过计算结合能判断稳定性等。
2、声子谱计算方法
声子谱的计算主要有两种方法，一种是直接法，另一种是微扰密度泛函方法(DFPT)。
直接法，或称frozen-phonon方法，是通过在优化后的平衡结构中引入原子位移，计算作用在原子上的Hellmann-Feynman力，进而由动力学矩阵算出声子色散曲线。直接法的缺陷在于它要求声子波矢与原胞边界supersize正交，或者原胞足够大使得Hellmann-Feynman力在原胞外可以忽略不计。这使得对于复杂系统，如对称性高的晶体、合金、超晶格等材料需要采用超原胞。超原胞的采用使计算量急剧增加，极大的限制了该方法的使用。
1987年，Baroni、Giannozzi和Testa提出了DFPT方法。DFPT通过计算系统能量对外场微扰的响应来求出晶格动力学性质。该方法最大的优势在于它不限定微扰的波矢与原胞边界正交，不需要超原胞也可以对任意波矢求解。Castep、Vasp等采用的是一种linear response theory 的方法(或者称为 density perturbation functional theory，DFPT)，直接计算出原子的移动而导致的势场变化，再进一步构造出动力学矩阵。进而计算出声子谱。
2.1 计算环境的搭建
计算所用的软件为VASP与phonopy code。
2.2 phonopy code编译
(1) 搭建Anaconda环境，下载链接(https://www.anaconda
.com)。
(2) 安装Anaconda
bash Anaconda3-2019.07-Linux-x86_64.sh
(3) 写入环境变量
echo "export PATH=/……/anaconda3/bin:$PATH" >> ~/.bashrcsource ~/.bashrc
Note： /……/表示自己服务器路径
(3) 下载phonopy code，下载链接(https://pypi.org/project/phonopy/)
(4) 编译phonopy code
tar zxvf phonopy-2.2.0.tar.gzcd phonopy-2.2.0python3 setup.py install --user
(5) 写入phonopy的环境变量
echo "export PATH=/……/phonopy-2.2.0/build/scripts-3.7:$PATH" >> ~/.bashrcsource ~/.bashrc
Note：/……/表示自己服务器路径
(6) 检查phonopy code是否安装完成
$ phonopy
2.3 声子谱计算步骤
(1) 结构优化
声子谱的计算需要对原胞结构做高精度的充分优化，否则很容易出现虚频；下面以二维InSe为例，详细说明声子谱的计算步骤。
INCAR
SYSTEM = 2D_InSe ISTART = 0 NWRITE = 2 PREC = Accurate ENCUT = 500 GGA = PE NSW = 200 ISIF = 3 ISYM = 2 NBLOCK = 1 KBLOCK = 1 IBRION = 2 NELM = 80 EDIFF = 1E-08 EDIFFG = -0.001 ALGO = Normal LDIAG = .TRUE. LREAL = .FALSE. ISMEAR = 0 SIGMA = 0.02 ICHARG = 2 LPLANE = .TRUE. NPAR = 4 LSCALU = .FALSE. NSIM = 4 LWAVE = .FALSE. LCHARG = .FALSE. ICORELEVEL = 1
POSCAR
2D_InSe1.0 4.0836000443 0.0000000000 0.0000000000 -2.0418000221 3.5
3650
13772 0.0000000000
0.0000000000 0.0000000000 25.3775005341 In Se 2 2Direct 0.666670026 0.333330010 0.589979992
0.666670026 0.333330010 0.478789968 0.333329978 0.666669998 0.428439986 0.333329978 0.666669998 0.640340007
KPOINTS
Monkhorst Pack of Gamma centered0Gamma 13 13 1 0.0 0.0 0.0
POTCAR
cat In_d/POTCAR Se/POTCAR > POTCAR
OPTCELL
100110000
不懂OPTCELL的设置请看刘博博客：
http://blog.wangruixing.cn/2019/05/05/constr/
(2) 通过phonopy code建立超胞
a) 将优化后的结构文件拷贝为POSCAR
cp CONTCAR POSCAR
b) 建立超胞
$ phonopy -d --dim="4 4 1"$ lsphonopy_disp.yaml POSCAR-002 POSCAR-005 POSCAR-008 POSCAR-011 POSCAR-014 SPOSCARPOSCAR POSCAR-003 POSCAR-006 POSCAR-009 POSCAR-012 POSCAR-015POSCAR-001 POSCAR-004 POSCAR-007 POSCAR-010 POSCAR-013 POSCAR-016
Note: 建立超胞之前，一定要在MS中找好对称性，否则会增加大量的计算量
c) POSCAR和SPOSCAR的重命名
cp POSCAR POSCAR-unitcellcp SPOSCAR POSCAR
(3) 计算力学Hessian矩阵
力学Hessian矩阵可由VASP计算得到，该矩阵保存于VASP的输出文件vasprun.xml中；计算过程中KPOINTS文件根据服务器计算能力，可适当增加，POTCAR文件无需改变，INCAR文件设置如下：
INCAR
SYSTEM = 2D_InSe ISTART = 0 NWRITE = 2 IBRION = 8 NSW = 1 IALGO = 38 NELM = 200 EDIFF = 1E-07 EDIFFG = -0.001 ISMEAR = 0 SIGMA = 0.02 ENCUT = 500 PREC = Accurate LREAL = .FALSE. LWAVE = .FALSE. LCHARG = .FALSE. ADDGRID = .TRUE
(4) 数据处理
a) 根据vasprun.xml文件生成力学文件FORCE_CONSTRAINS
phonopy --fc vasprun.xml
b) 编辑band.conf文件，该文件给出了高对称点路径的信息
ATOM_NAME = In SeDIM = 4 4 1BAND = 0.5 0.0 0.0 0.0 0.0 0.0 0.333333 0.333333 0.0 0.5 0.0 0.0FORCE_CONSTANTS = READ
Note：
line 1：元素名称，顺序与POSCAR保持一致
line 2：建立超胞的大小
line 3：高对称点路径，每组高对称点之间用两个空格隔开
line 4：表示读取力常数文件FORCE_CONSTANTS
c) 生成band.yaml文件
phonopy --dim="4 4 1" -c POSCAR-unitcell band.conf
Note：4 4 1为建立超胞的大小
d) 得到声子谱数据文件PBAND.dat文件
phonopy-bandplot --gnuplot> PBAND.dat
Note: 高对称点标注说明：phonopy软件默认在两个高对称点之间打点51个，且在PBAND.dat中每组高对称点数据以一个空行分割，据此即可得到完整的声子谱图像
e) 计算结果
文献结果：
3、后记
如果计算完的声子谱其他地方都好，就是在Γ点有一点点虚频，那么这个材料很可能是稳定的，只是你优化做得不够好，进一步提高优化的精度可消除这一点点虚频。
对于二维材料，如果在Γ点出现很小的虚频，基本可以认为这个材料是稳定的，大部分二维材料都会有此现象；尤其是VASP结合phonopy code计算二维材料的声子谱在Γ点更是容易出现虚频；使用Quantum-ESPRESSO的PWSCF和PH模块计算声子谱对于内存的需求较小一些，且对于二维材料的声子谱计算更友好一些，后期会详细介绍Quantum-ESPRESSO计算声子谱的具体步骤。
请大家关注贺勇的个人博客网站：
https://yh-phys.github.io
,
