在第一性原理計算過程中考察體系穩定性是經常遇到的一個問題，也是大部分審稿人容易提問的地方；通常，使用聲子譜研究體系的動力學穩定性，使用分子動力學研究體系的熱穩定性；另外，可以藉助於一些其它方法進一步說明體系的穩定性，例如研究鈣鈦礦體系時通過離子半徑判斷穩定性，研究異質結時通過計算結合能判斷穩定性等。
聲子譜的計算主要有兩種方法，一種是直接法，另一種是微擾密度泛函方法(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 ~/.bashrcNote： /……/表示自己伺服器路徑
(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 ~/.bashrcNote：/……/表示自己伺服器路徑
(6) 檢查phonopy code是否安裝完成
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.5365013772 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
不懂OPTCELL的設置請看劉博博客：
http://blog.wangruixing.cn/2019/05/05/constr/
(2) 通過phonopy code建立超胞
a) 將優化後的結構文件拷貝為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
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.confNote：4 4 1為建立超胞的大小
d) 得到聲子譜數據文件PBAND.dat文件
phonopy-bandplot --gnuplot> PBAND.dat
Note: 高對稱點標註說明：phonopy軟體默認在兩個高對稱點之間打點51個，且在PBAND.dat中每組高對稱點數據以一個空行分割，據此即可得到完整的聲子譜圖像
e) 計算結果
文獻結果：
如果計算完的聲子譜其他地方都好，就是在Γ點有一點點虛頻，那麼這個材料很可能是穩定的，只是你優化做得不夠好，進一步提高優化的精度可消除這一點點虛頻。
對於二維材料，如果在Γ點出現很小的虛頻，基本可以認為這個材料是穩定的，大部分二維材料都會有此現象；尤其是VASP結合phonopy code計算二維材料的聲子譜在Γ點更是容易出現虛頻；使用Quantum-ESPRESSO的PWSCF和PH模塊計算聲子譜對於內存的需求較小一些，且對於二維材料的聲子譜計算更友好一些，後期會詳細介紹Quantum-ESPRESSO計算聲子譜的具體步驟。
請大家關注賀勇的個人博客網站：
https://yh-phys.github.io