3.7.3 晶格的振动自由能

声子计算不单单可以给出声子谱来说明体系的动力学稳定性，还可以从声子态密度数据中获得晶格振动的自由能。前面我们提到过，对于一个结构单纯做DFT静态自洽计算得到的能量，是不包含晶格振动的自由能的。要获得晶体的振动自由能，比较主流的方法之一就是算声子，其晶格振动的自由能与声子态密度的关系表达式如下（具体推导见固体物理教材）：

$$F_{vib}\left(V,T\right)={\sum }_{q,j}\left[\frac{\hslash {\omega }_{qj}\left(V\right)}{2}+k_{B}T\ln \left(1-\exp \left(-\frac{\hslash {\omega }_{qj}\left(V\right)}{k_{B}T}\right)\right)\right]$$

求和写成积分形式则是：

$$F_{vib}\left(V,T\right)={\int }_0^{+\infty }\left(\frac{\hslash \omega \left(V\right)}{2}+k_{B}T\ln \left(1-\exp \left(-\frac{\hslash \omega \left(V\right)}{k_{B}T}\right)\right)\right)\rho \left(\omega \right)d\omega$$

其中$\rho \left(\omega \right)$就是声子态密度。晶体的总自由能可以写成：

$$F(V,T)=F_{elec}(V,T)+F_{vib}\left(V,T\right)$$

其中$F_{elec}(V,T)$可以直接从OUTCAR的能量中的Free energy部分获得，即OSZICAR中的F（注意，不是E0）。DFT直接做relax可以获得0K下的体积（没有原子热振动效应），如果假设体积不随温度膨胀的话，那么用该体积算一次声子，得到声子态密度；并且假设不考虑电子熵，再结合上面的公式就可以得到体系的自由能随温度的变化。这一“体积不随温度变化”的假设就是简谐近似，Harmonic Approximate（HA）。更严格一点的是准谐近似Qusia-Harmonic Approximate（QHA）；QHA考虑了晶格热膨胀效应，但是忽略了温度对力常数（频率）的修正。QHA常用于计算固态物质的热膨胀系数、或是高温高压下的相图。

现在我们用phonopy来获得给定体积V下，体系的晶格振动的自由能数据：

只需要在前面的KPATH.phonopy基础上，加上如下参数：

TPROP = .TRUE.
TMIN = 0
TMAX = 1000
TSTEP = 100

然后运行phonopy -c xxx KPATH.phonopy那一串命令即可。此时band相关的参数不起作用并且运行过程中屏幕会输出当前体积下晶格振动自由能随温度的变化，如下图。第二列F就是晶格振动自由能，E是晶格振动的内能，S是晶格振动熵，F=E-TS。0K下晶格振动内能不为0，因为谐振子存在零点能，具体内容请复习固体物理。并且这里的mol是以原胞为单元的，并非每原子，可以从thermal_properties.yaml文件前几行中了解这一点。

#      T [K]      F [kJ/mol]    S [J/K/mol]  C_v [J/K/mol]     E [kJ/mol]
       0.000      79.4251899      0.0000000      0.0000000     79.4251899
     100.000      79.2614195      7.1856905     23.1953433     79.9799886
     200.000      76.8294562     45.4642063     94.9035257     85.9222974
     300.000      69.8479486     94.2311910    144.4295012     98.1173059
     400.000      58.0944499    139.9832679    172.3314374    114.0877570
     500.000      42.0339920    180.3250134    188.4559003    132.1964988
     600.000      22.1977920    215.6263310    198.3321427    151.5735906
     700.000      -0.9509022    246.7132646    204.7302454    171.7483830
     800.000     -27.0301131    274.3526452    209.0788710    192.4520031
     900.000     -55.7273834    299.1660850    212.1551338    213.5220931
    1000.000     -86.7855026    321.6411732    214.4049800    234.8556705

单位换算，1KJ/mol~0.0103643eV/f.u.

注意，计算晶格振动自由能的前提是体系的声子谱和态密度总没有虚频。否则得到的自由能是有问题的。当然，phonopy提供了一种灵活但不保证准确性的处理方式，即把虚频设置为正值，通过PRETEND_REAL参数实现，具体见phonopy官网的参数说明。

https://phonopy.github.io/phonopy/setting-tags.html#thermal-properties-tag

上面获得晶格振动自由能是建立在简谐近似的基础，即认为晶体没有热膨胀。对于大部分体系，热膨胀效应不可忽略，可以通过QHA来研究。QHA的教学不属于基础计算教学内容，请自行从phonopy官网学习。

提示：在做QHA过程中时，需要拟合F-V数据。如果不加振动自由能，F-V数据很光滑；如果加入振动自由能，F-V数据可能没那么光滑；可能原因是声子计算精度不够，考虑INCAR中加入ADDGRID = .TRUE. ; LREAL = .FALSE. ; PREC = Accurate三行参数。