3.4 非谐声子

发表于 2023-05-09 更新于 2024-10-29 分类于 DFT学习

介绍

ALAMODE中非谐声子计算部分采用了自洽声子计算方法。

具体教程：https://alamode.readthedocs.io/en/latest/anphondir/formalism_anphon.html#self-consistent-phonon-scph-calculation

https://www.bilibili.com/video/BV1mv411M7W7/?spm_id_from=333.999.0.0&vd_source=17084cc867ff3bb86398ff1a79b443f2

TDEP 是其采用技术Temperature dependent effective potential的缩写。
下面是开发者的网站，但是并不能下载到源代码。http://ollehellman.github.io/page/index.html好消息是第一性原理计算软件ABINIT中集成了TDEP方法，官网user guide>a-TDEP。需要研究的可以自行研究。ABINIT官网：https://www.abinit.org/

Dynaphopy 是采用了一种基于AIMD（第一性原理分子动力学）的杂化方法的python程序包：http://abelcarreras.github.io/DynaPhoPy/index.html
文中给出了参考文献，需要的自行下载。

ALAMODE

安装

https://alamode.readthedocs.io/en/latest/index.html#

ALAMODE计算Si简谐声子谱

1、利用ALM计算displacement patterns

优化后的超胞结构重命名POSCAR.orig，准备输入文件si_alm.in：

&general
  PREFIX = si222
  MODE = suggest
  NAT = 64; NKD = 1
  KD = Si
/

&interaction
  NORDER = 1  # 1: harmonic, 2: cubic, ..
/

&cell
  20.406 # factor in Bohr unit
  1.0 0.0 0.0 # a1
  0.0 1.0 0.0 # a2
  0.0 0.0 1.0 # a3
/

&cutoff 
  Si-Si None
/

&position
  1 0.0000000000000000 0.0000000000000000 0.0000000000000000   
  1 0.0000000000000000 0.0000000000000000 0.5000000000000000
  1 0.0000000000000000 0.2500000000000000 0.2500000000000000
  1 0.0000000000000000 0.2500000000000000 0.7500000000000000
  1 0.0000000000000000 0.5000000000000000 0.0000000000000000
  1 0.0000000000000000 0.5000000000000000 0.5000000000000000
  1 0.0000000000000000 0.7500000000000000 0.2500000000000000

晶格常数替代为自己的超胞结构。

运行ALM: alm si_alm.in > si_alm.log1，将会生成si222.pattern_HARMONIC文件，包含坐标信息。对于Si只有一个displacement patterns。

2、计算原子间力

在此之前要确定位移的大小，一般简谐0.01-0.04，0.01大多可以。运行python displace.py --VASP=POSCAR.orig --mag=0.01 -pf si222.pattern_HARMONIC

之后对每个结构提交计算，脚本：

#!/bin/bash

# Assume we have 20 displaced configurations for QE [disp01.pw.in,..., disp20.pw.in].

for ((i=1;i<=20;i++))
do
    num=`echo $i | awk '{printf("%02d",$1)}'`
    mkdir ${num}
    cd ${num}
    cp ../disp${num}.pw.in .
    pw.x < disp${num}.pw.in > disp${num}.pw.out
    cd ../
done

计算完成后，提取信息：python extract.py --VASP=POSCAR.orig vasprun*.xml > DFSET_harmonic，生成DFSET_harmonic文件。

3、拟合力常数

修改1中的 si_alm.in ：

&general
  PREFIX = si222
  MODE = optimize   # <-- here
  NAT = 64; NKD = 1
  KD = Si
/

&optimize
  DFSET = DFSET_harmonic
/

运行alm si_alm.in > si_alm.log2，生成si222.fcs and si222.xml

4、计算声子谱和DOS

准备si_phband.in文件，结构中用晶胞结构，注意&cell中单位是Bohr，具体可设置如下

&general
  PREFIX = si
  MODE   = phonons
  FCSXML = si.xml

  NKD = 1; KD = Si
  MASS = 28.085
/

&cell
  10.184247326051628 
  0.0000000000000000    0.5071343999939496    0.5071343999939496
  0.5071343999939496    0.0000000000000000    0.5071343999939496
  0.5071343999939496    0.5071343999939496    0.0000000000000000
/

&kpoint
  1   # KPMODE = 1: line mode
  G 0.0 0.0 0.0 X 0.5 0.5 0.0 51
  X 0.5 0.5 1.0 G 0.0 0.0 0.0 51
  G 0.0 0.0 0.0 L 0.5 0.5 0.5 51
/

执行anphon si_phband.in > si_phband.log会产生si.bands文件，可以用plotband.py脚本画图。执行python plotband.py si.bands得到

准备si_phdos.in文件，结构中用晶胞结构，注意&cell中单位是Bohr，具体可设置如下

&general
  PREFIX = si
  MODE   = phonons
  FCSXML = si.xml

  NKD = 1; KD = Si
  MASS = 28.085
/

&cell
  10.184247326051628 
  0.0000000000000000    0.5071343999939496    0.5071343999939496
  0.5071343999939496    0.0000000000000000    0.5071343999939496
  0.5071343999939496    0.5071343999939496    0.0000000000000000
/

&kpoint
  2   # KPMODE = 2: uniform mesh mode
  20 20 20
/

执行anphon si_phdos.in > si_phdos.log会产生si.dos文件，可以用plotdos.py脚本画图。执行python plotdos.py –emax 550 –nokey si.dos得到

5. 估算对热导率的立方IFCs(非简谐IFCs)

这里计算非简谐IFCs，我们把si_alm1.in复制为si_alm3.in，将NORDER = 1改为NORDER = 2.我们还要注意这里考虑计算三阶力常数，所以&cutoff下的参数我们要再指定个截止半径，一般要求略大于第二近邻的距离（可以从si_alm1.log查看，如下图）。当然你接着用None也行，这样可能大大增加参数的数量，从而增加计算成本。

运行alm si_alm3.in > si_alm3.log，可以看到多了si_cubic.相关文件。运行python displace.py –OpenMX=Si.dat –mag=0.04 -pf si_cubic.pattern_ANHARM3，可见产生了disp01.dat到disp11.dat11个OpenMX输入文件，我们准备shell脚本批量提交，shell脚本命名为pbs.sh，具体为

#!/bin/bash

for ((i=1;i<=11;i++))
do
    num=`echo $i | awk '{printf("%02d",$1)}'`
    mkdir ${num}
    cd ${num}
    cp ../disp${num}.dat .
    cp disp${num}.dat Si.dat
    cp ../openmx.pbs .
    qsub openmx.pbs
    cd ../
done

总结

简谐计算步骤

输入文件：优化后的超胞结构重命名POSCAR.orig，准备输入文件si_alm.in（其中mode：suggest）；输出：si222.pattern_HARMONIC文件，包含坐标信息。
displace.py对上面的.pattern文件处理，得到随机位移结构；对随机位移机构进行力的计算，生成不同vasprun.xml；extract.py处理上面vasprun.xml，生成DFSET_harmonic文件。
力常数拟合输入文件：DFSET_harmonic文件，修改1中的 si_alm.in，其中mode改成opt, 添加下面；生成si222.fcs and si222.xml。
1
2
3
&optimize
DFSET = DFSET_harmonic
/
后处理输入文件：准备phband.in文件，修改MODE = phonons
FCSXML = si.xml；修改kmesh格式可得到声子谱、态密度、热性质等。

非谐计算

准备简谐计算得到的si222.fcs and si222.xml，AIMD得到的vasprun.xml，超胞结构。
displace.py 对AIMD的vasprun.xml处理，得到随机结构；对随机位移机构进行力的计算，生成不同vasprun.xml；extract.py处理上面vasprun.xml，生成DFSET_AIMD_random文件。
Cross validation (CV)：alm_cv.in文件中添加简谐力常数 FC2XML = ../harmonic.xml，添加高阶力常数部分标签和cv相关标签，mode = cv；得到*cvset文件；python3 cvscore.py *cvset > si222.cvscore得到Minimum cvscore
拟合力常数输入文件：alm_opt.in中mode=opt；修改cv相关标签；生成高阶力常数cv.fcs and cv.xml.在这一步alm 的输入文件&general 部分加入 FC3_SHENGBTE=1、 FC4_SHENGBTE=1可输出ShengBTE格式的输入文件。

SCPH计算：scph.in文件中修改MODE=SCPH和温度标签；添加&scph模块，见上面例子。anphon scph.in > scph.log提交，查看收敛：grep “conv” scph.log。生成输出文件：.scph_dymat、.scph_dfc2、.scph_bands。对.scph_dfc2处理得到有效二阶力常数文件，用于后续计算。

>dfc2
DFC2--a generator of renormalized harmonic FCs from SCPH outputs.
XML file containing original Fc2 : STo222.xml 
Output xml filename with anharmonic correction : ST0222_SCPH2-2_300K.xml
FC2 correction file fromSCPH calculation :STo_scph2-2.scph_dfc2
Target temperature : 300


New XML file STo222_SCPH2-2_300K.xml was created successfully.

Alamode计算高温声子色散

Anharmonic Phonons (self-consistent phonon theory)
主要使用alamode软件。版本：alamode-1.1.0

输入文件

vasprun.xml #分子动力学输出文件
Dfile  #包含原子的力和位移
extract_disp.py #选取结构，并计算位移
extract_force.py #选取结构，并获得受力
PPOSCAR #初始超胞
alm.in #提取力常数的输入文件
anphono.in #计算声子谱的输入文件

从分子动力学数据中选取一定数量的结构，并提取它们的力和计算与初始状态相比的位移，输出到Dfile文件当中

python extract_disp.py
#运行后根据提示输入开始的步数和间隔步数
#例如总共1000步，从第一步开始提取，每间隔五步提取一次
#得到disp.dat

python extract_force.py

方法同上得到force.dat
将disp.dat和force.dat文件按如下形式合并到一起

disp1x  disp1y  disp1z   force1x  force1y  force1z
       ......
       ......
dispNx  dispNy  dispNz   forceNx  forceNy  forceNz

获取力常数

提取四阶力常数的alm.in文件如下（使用LASSO方法）：

二阶力常数可以通过有限位移法提取(详见官网手册)

&general
  PREFIX = ThO2_ENET
  MODE = optimize
  NAT = 81; NKD = 2
  KD = Th O
#  PRINTSYM=1
/

&interaction
  NORDER = 5  # 1: harmonic, 2: cubic, ..
  NBODY = 2 3 3 2 2
/

&cell
  1 # factor in Bohr unit
  0.0000000000000   15.919498944104   15.919498944104
  15.919498944104   0.0000000000000   15.919498944104
  15.919498944104   15.919498944104   0.0000000000000
/

&optimize
  LMODEL= enet
  FC2XML = ThO2_2nd.xml
  DFSET=ThO2_dffile

  CV = 0
  L1_RATIO = 1.0
  L1_ALPHA = 1.0e-06
  CV_MINALPHA = 1.0e-6
  CV_MAXALPHA = 0.02
  CV_NALPHA = 100

  STANDARDIZE = 1
  CONV_TOL = 1.0e-9
/

&cutoff 
  *-* None None 12.0 12.0 12.0
/

&position
1  0.0000000000000000  0.0000000000000000  0.0000000000000000
1  0.3333333333333333  0.0000000000000000  0.0000000000000000
..................

注意晶格常数单位的转换 $1 Bohr= 0.529177208 Angstrom$

在终端运行：

bash
alm alm.in > alm.log
得到Nb.xml

Anharmonic Phonon

anphono.in输入文件如下：

&general
 PREFIX =scph
 MODE = SCPH
 FCSXML =Nb.xml #forceconstant file

 NKD = 1; KD = Nb
 MASS = 92.90638

 #NONANALYTIC = 0; BORNINFO = PbTe.born

 TMIN = 0; TMAX = 300; DT = 100   #Temperature and step
/

&scph
 KMESH_SCPH        = 6 6 6
 KMESH_INTERPOLATE = 6 6 6

 SELF_OFFDIAG = 0
 RESTART_SCPH = 0

 MIXALPHA = 0.1
 MAXITER = 500
 TOL_SCPH = 1.0e-10
/

&cell   ## primitive cell is best
  6.2507550151332
  -0.5  0.5  0.5
   0.5 -0.5  0.5
   0.5  0.5 -0.5
/

&kpoint
  1  # KPMODE = 1: get phonon band  2: get others properties
  G 0.0 0.0 0.0 H 0.5 -0.5 0.5 201
  H 0.5 -0.5 0.5 P 0.25 0.25 0.25 201
  P 0.25 0.25 0.25 G 0.0 0.0 0.0 201
  G 0.0 0.0 0.0  N 0.0 0.0 0.5 201
/

在终端运行：

anphon anphono.in > anphono.log
得到scph.scph_bands，里面有不同温度下的声子谱

进入安装目录下的tools文件夹：make Makefile，编译df2c。

[scf1355@ln151%nc-l ah]$ ~/soft/alamode/tools/dfc2 
 DFC2 -- a generator of renormalized harmonic FCs from SCPH outputs.
 XML file containing original FC2 : crte2-442.xml
 Output xml filename with anharmonic correction : crte2-new.xml
 FC2 correction file from SCPH calculation : crte2-442.scph_dfc2
 Target temperature : 300

输出crte2-new.xml.利用https://github.com/ttadano/alamode/discussions/192 中提到的脚本转换成phonopy格式。此时的力常数文件由于对称性数量不匹配，利用phonopy读取并写出FULL_FORCE_CONSTANTS。

Alamode_to_ShengBTE

Alamode计算热导率步骤

所有的计算都是基于密度泛函软件 VASP 来进行的。
首先利用第一性原理的密度泛函对晶体结构进行优化，弛豫出稳定的晶格结构，然后利用密度泛函微扰理论计算出波恩有效电荷和静态介电张量，以便用于后续求解动力学矩阵的非解析部分。
随后为了获得更加精准的二阶力常数，通过差分位移的方法生成微扰结构，后利用静态密度泛函理论计算出力和位移信息，再利用最小二乘法得到所需要的二阶力常数。
对于高阶力常数，本文利用分子动力学模拟得到 80 个准随机结构，对每个准随机结构的原子施加随机方向的微小位移，以便得到更加随机的结构，最后利用最小绝对收缩和选择算子的方法得到高阶力常数。
结合上述已经得到的介电张量、波恩有效电荷、二阶力常数和三阶力常数作为输入，最后通过ShengBTE声子玻尔兹曼输运方程求得简谐近似下晶格热导率、声子群速度、声子散射率等热输运系数。
结合介电张量、波恩有效电荷、二阶力常数和四阶力常数作为输入进行自洽声子理论的计算，这将得到与温度相关的非谐动力学矩阵，对其进行傅里叶变换得到与温度相关的有效二阶力常数，用有效二阶力常数代替上述的二阶力常数，在此基础上求解声子玻尔兹曼输运方程求得自洽声子理论下晶格热导率、声子群速度、声子散射率等热输运系数。
在简谐近似下存在虚频，与实验观察到的稳定性不一致。因此，传统的基于简谐原子间力常数的玻尔兹曼输运方程解的结果不再有效，并且非谐效应在足够大的位移振幅内变得显著。在 ALAMODE 软件包中，通过有限位移法在 2×2×2 超胞中获得简谐原子间力常数 (IFCs) ，位移大小设置为 0.01Å。
随后，通过压缩感知晶格动力学方法训练三阶和四阶力常数。具体来说，本文首先模拟了在 300 K、时间步长为 2fs 的 4000 步的从头算分子动力学，得到了 80 个快照。为了得到准随机结构，本文将到得到快照中的原子向随机方向移动 0.1 Å。随后，本文将获得的准随机结构在 4×4×4 的 k 点网格中进行了静态 DFT 计算，以计算 Hellmann–Feynman 力。最终，利用得到的 80 个准随机结构的位移和力组成的数据集通过最小绝对收缩和选择算子 (least absolute shrinkage and selection operator，简称 LASSO)方法训练出非谐力常数。
在计算中本文考虑了超胞内的所有原子间最近邻相互作用来计算简谐和三阶力常数，四阶力常数只考虑了第五最近邻内的相互作用，由于五阶和六阶力常数对结果的影响很小，因此只考虑了次近邻相互作用。
本文根据 ALAMODE 软件包中采用的 SCP 理论，利用获得的简谐和非谐力常数来计算非谐声子频率。
声子玻尔兹曼输运方程的解是通过 ShengBTE 的修订版 FourPhonon 软件包来实现的，并建立了具有高斯涂抹 (scale broad 参数为 0.05) 的 18×18×18 的 q 点网格来模拟相关积分和声子波矢量。此外，在计算中得到了大约 2.3×105 个三声子过程和 1.2×109 个四声子过程。因此，对于三声子散射，本文采用迭代方案来求解玻尔兹曼输运方程，而通过单模弛豫时间近似(SMRTA) 来处理四声子散射过程，这是因为计算成本巨大。

步骤分析

上面第3步，利用alm来算大概分为三步，mode设为suggest产生随机位移，vasp求力和能量，mode设为opt拟合二阶力常数即可。alm模块。
上面第4步，利用AIMD产生随机位移，vasp求力和能量，压缩感知晶格动力学方法训练三阶和四阶力常数mode设为cv，随后mode设为opt拟合高阶力常数，alm模块。最后mode设为SCPH获得变温声子谱和有效二阶力常数,anphon模块。后续运行dfc2模块生成指定温度的有效二阶力常数，用https://github.com/ttadano/alamode/discussions/192 中提到的脚本转换成phonopy格式。
上面第5步，为了产生ShengBTE的力常数输入文件，在拟合高阶力常数时，alm 的输入文件&general 部分加入 FC3_SHENGBTE=1、 FC4_SHENGBTE=1即可。如果简谐声子谱有虚频，上面的有效二阶力常数也需要，即第7步所述。
上面第6步见ShengBTE教程3.2节。
Alamode可以得到的结果包括格林艾森常数、声子群速度、三声子相空间、热导率原子贡献率和元素、变温声子谱和态密度、MSD、可视化振动模式、QHA等,anphon模块。

Alamode 连接 shengBTE

首先，github上安装最新版本的develop branches 或 feature/shengbte_4ph branches

[3rd 4th IFCs]

拟合力常数时，alm 的输入文件&general 部分加入 FC3_SHENGBTE=1、 FC4_SHENGBTE=1即可。

&general

PREFIX = BN-1000

MODE = opt

NAT = 128

NKD = 2; KD = B N

TOLERANCE = 1.0e-3

FC3_SHENGBTE=1

FC4_SHENGBTE=1

Note：输出的文件名为：FORCE_CONSTANT_3RD、FORCE_CONSTANT_4TH，和ShengBTE的文件名差一个S，FORCE_CONSTANTS_3RD、FORCE_CONSTANTS_4TH

[2nd IFCs]

首先通过QE的DFPT计算一个不含SCP的二阶力常数，BN.ifc2。

做SCPH计算，会生成计算温度点的二阶力常数SCP等效修正部分，例如BN.scph_dfc2（仅考虑四声子散射对重整化的影响），若进一步考虑三声子散射的重整化作用，文件名为BN.scph+bubble(0)_dfc2。

alamode/tools/scph_to_qefc.py 可以将SCP的修正部分与QE计算的二阶力常数相加，得到指定温度重整化后的二阶力常数。

Usage:

python scph_to_qefc.py original_QEfc2 scph_fc2_correction temperature > scp_fc2

python scph_to_qefc.py BN.ifc2 BN.scph_dfc2 1000 > espresso.ifc2

Authors and references

Weizhe Yuan yuanwz@stu.hit.edn.cn

W. Yuan et al., Physical Review B 110, 165412 (2024).

注意的问题

The accuracy of the anharmonic force constants can be checked by increasing the size of the training dataset. If the physical properties, such as thermal conductivity and free energy, change significantly with the size of the training dataset, it indicates that the training is insufficient. If the change is minor, your force constants may be OK.
报错：Warning: The following force constant doesn’t exist in the original file:
You can neglect this warning as long as the correction term (the last column) is small.
This warning is raised when the absolute value of the correction is larger than 1.0e-10 and the corresponding force constant does not exist in the original FCSXML file. In your case, the correction term is only slightly above the tolerance. So, it should not cause any problems hopefully.
When the correction term is much larger than 1.0e-10, this warning indicates a possible error in the inputs. Likely errors are the input structure is not fully symmetrized, or the numerical digits of lattice vectors and fractional coordinates are insufficient.
the given KMESH_INTERPOLATE is not commensurate with the supercell size of the original FCSXML.
In the current implementation of alamode, only the former choice,
renormalized 2FCs + conventional (old) 3FCs,
is supported. Another possible choice may be
renormalized 2FCs + renormalized 3FCs,
but it is not supported within alamode yet.
When MODE = RTA, the anphon code only reads the 2FCs and 3FCs from given XML files. By contrast, when MODE = SCPH, the code additionally reads the 4FCs because they are necessary to compute the renormalized 2FCs based on the self-consistent phonon theory.
Therefore, the second combination of
renormalized 2FCs + high order FCs using SCPH
is not possible.
I performed the SCPH calculation. It is not converged at low temperatures (0K, 100K, 200K), while it achieved convergence at high temperatures (300~1200K).A potential solution is to use a smaller MIXALPHA and a smaller DT.
I would like to use second-order, third-order, and fourth-order force constants to obtain anharmonic potentials and calculate the potential energy as a function of atomic displacement. My goal is to generate a plot similar to the one in the paper “Orbitally driven giant phonon anharmonicity in SnSe” (Nature Phys. 2015).Please use https://github.com/ttadano/alamode/blob/develop/tools/calcpes.py.
now that I have obtained the renormalized force constant file xxx.xml at a finite temperature, I now want to convert the renormalized force constant file directly to FORCE_CONSTANS
How to convert a second-order force constant matrix obtained by almode to FORCE_CONSTANTS in phonopy format？

python convert.py alamode.xml > FORCE_CONSTANTS

Could you give me some suggestion? How shoud I covert this with phonopy?

FULL_FORCE_CONSTANTS = .TRUE.
WRITE_FORCE_CONSTANTS = .TRUE.
Note: the old FORCE_CONSTANTS will be overwrited

How small the fitting error should be?

It depends on the Taylor expansion potential and displacement magnitude you choose.
In the standard harmonic calculation where –mag=0.01 is used in displace.py and the all harmonic interactions are considered (cutoff = None), the fitting error is usually less than 5% (~1–2% in most cases).
In the calculation of third-order force constants with –mag=0.04, the fitting error should be small as well. Indeed, in many cases, we obtain much smaller fitting errors (< 1%) than the harmonic case.
In the temperature-dependent effective potential method, we try to fit the harmonic potential to the displacement-force datasets sampled by ab initio molecular dynamics at finite temperature. Therefore, the fitting error tends to be much larger (> 10%).

dynaphopy

简谐声子谱的计算。
dynaphopy的非谐修正是在phonopy的力常数基础上进行的，需要先进行0k声子谱的计算。
结构优化的INCAR，力的收敛标准是E-3，大部分的文献也是在这个精度范围内。结构优化的结构是单胞，10个原子。
这是phonopy计算的声子谱，计算的结构为222的超胞，80个原子，在gmma点和X点有很强的虚频。这个虚频即使提高收敛精度和扩大超胞也是无法消除的。
AIMD计算
dynaphopy通过读取AIMD的计算结构进行拟合（说法不准确），需要预先准备AIMD的outcar或xdatcar。
AIMD的计算设置，基本参考dynaphopy的example中的设置，详细的参数涵义，都有教程介绍。
这里只是采用了3000步的计算，步长是2fs，能量的收敛精度是E-6，不算太高。同样是采用80个原子的超胞。大部分做非谐计算的文献，所采用的都是在80-160原子数范围内的超胞，步长多是1-2fs，总步数大概是3000-6000，总的时长是5ps-10ps。从下面的非谐声子谱结果来看，这个设置对于非谐计算是合适的？
dynaphopy非谐声子谱计算。
dynaphopy的官网如下：Dynaphopy (abelcarreras.github.io)。
将结构优化的原始晶胞POSCAR，phonopy输出的力常数，和AIMD的OUTCAR放在一个文件夹内。
设置dynaphopy的输入文件input如下：分别是输入结构文件，力常数文件，原胞和超胞与POSCAR的关系设置，band是声子谱的高对称点路径，设置方式和vasp计算PBE能带一致。
然后输入dynaphopy input OUTCAR -i，进入下的界面，选择6，进入非谐声子谱计算，再选择1，可以同时输出简谐声子谱和非简谐声子谱，此时会对peak进行拟合。拟合结束后，会输出声子谱。
一开始使用默认设置，输出的声子谱，可以看到高温声子谱虚频加重了。从@get-it 大佬处了解到，dynaphopy默认使用最大熵方法拟合，效果可能不如fftw有效。于是小弟将拟合方式改为fftw,通过指令dynaphopy input OUTCAR -i -psm 3,采用快速傅里叶变换拟合。得到300k下的声子谱，可以看到，原来gamma点和X点的虚频，已经被成功消除了。而实验上，这个材料在室温下肯定是稳定的。

当然，如果想要得到更准确的声子谱，需要提高优化精度，采用更大的超胞计算力常数，AIMD需要更多的步数和更小的步长，比如dynaphopy的例子是采用0.7 fs,跑200000步。如果只是想得到一个比较合理的声子谱，适当降低精度是没问题的。比如这个计算，在24核的机器上，结构优化用时不到3小时，dfpt计算1小时，AIMD计算14小时，差不多用时一天，能够得到一个比较合适的声子谱。如果提高精度，计算用时可能会提升数倍，但是声子谱的提升可能并不会太明显。