上一节我们介绍了Phonopy的安装和如何用Phonopy计算声子谱，这一节我们来学习如何用Phonopy计算声子态密度，自由能，和热膨胀系数。话不多说，赶紧来学习吧。

用phonopy算声子态密度

在计算声子谱的基础上，如果要获得声子态密度可以准备mesh.conf:

ATOM_NAME = Bi O Se #元素符号

DIM = 5 7 1 #扩胞倍数

MP = 71 71 1 #取点密度

声子态密度可以由phonopy直接生成，执行命令

phonopy --dim='2 4 1' -p POSCAR-unitcell mesh.conf

可得到声子态密度。

用phonopy算自由能

自由能是热力学中的重要概念，用于衡量体系对外做功的上限。常用的自由能分两种：

一种叫做亥姆霍兹自由能，另一种叫吉布斯自由能。在等温过程中体系对外做功的最大值不会超过体系亥姆霍兹自由能的减小值，而在等温等压过程中体系对外做功的最大值不会超过体系吉布斯自由能的减小值。

自由能的高低还可以被用于判断体系是否会发生相变，我们在常温常压下进行实验时，吉布斯自由能往往更加有用。但如果体系压强和体积的乘积远小于体系的亥姆霍兹自由能时，我们就可以近似用亥姆霍兹自由能来代替体系的吉布斯自由能。

亥姆霍兹自由能可以通过phonopy+VASP很方便的计算出来。现给出其计算过程。

首先建立超胞，然后输入

phonopy -d --dim="1 1 1"

建立文件夹force_001，里面放入四个文件INCAR，POSCAR，POTCAR和KPOINTS

INCAR如下：

SYSTEM=C1Mo1

IBRION =-1

IALGO=38

EDIFF=1E-7

ISMEAR=0

SIGMA=0.01

ENCUT=500

PREC=Accurate

LWAVE=.FALSE.

LCHARG=.FALSE

ADDGRID=.TRUE

POSCAR用POSCAR-001，POTCAR和KPOINTS根据结构信息设置。依次建立force_002、force_003……重复上述过程，POSCAR分别用POSCAR-002、POSCAR-003……将force系列的文件统一运行，运行完毕后输入命令：

phonopy –f ./force_001/vasprun.xml ./force_002/vasprun.xml……

运行

phonopy –t --tmax=2000 mesh.conf

其中mesh.conf为

DIM= 1 1 1

ATOM_NAME = C Mo

MP = 5 5 7

通过上面的步骤算出来的只是晶格振动那一部分的自由能，要想得到整个材料的自由能还需要把电子那部分自由能加上去。

用phonopy算热膨胀系数

物体由于温度改变而有胀缩现象。其变化能力以等压(p一定)下，单位温度变化所导致的长度量值的变化，即热膨胀系数表示。各物体的热膨胀系数不同，一般金属的热膨胀系数单位为1/℃（摄氏）。线胀系数是指固态物质当温度改变摄氏度1℃时，其某一方向上的长度的变化和它在20℃（即标准实验室环境）时的长度的比值。

使用软件：VASP，phonopy。下载phonopy官方文件中Al-QHA的example

官方例子中只将缩放系数1改变为0.95，而分数坐标没有发生变化，可以认为在第一次晶格优化以后，通过改变缩放系数改变晶胞体积，然后进行离子弛豫。

因此计算材料的热膨胀系数步骤为：

1.先进行高精度的晶格优化。

2.修改缩放系数为0.95-1.05新建10个POSCAR.

3.对11个POSCAR晶体进行声子谱计算

4.thermo.conf文件后处理生成11个对应的thermal_properties.yaml

5.QHA计算热膨胀系数

thermo.conf文件

DIM = 2 2 2

ATOM_NAME = sio

MP = 31 31 31

TPROP = .TRUE.

TMAX = 2000

TSTEP = 2

FORCE_CONSTANTS = READ

下面我们写e-v.dat文件

首先找到之前计算的那11个POSCAR的vasprun.xml结果文件，然后在vasprun.xml文件中搜索最后一个e_wo_entrp。

下载phonopy官方文件中Si-QHA的example，解压vasprun_xmls.tar.lzma我们可以得到算例中用到的vasprun.xml文件。

首先看这两行

<i type="int" name="IBRION"> -1</i>

<i type="int" name="NSW"> 0</i>

我们知道这个不是声子谱计算的结果文件应该是原子弛豫后自洽计算的结果文件，我们找到1.05倍缩放系数对应的哪个结果文件的最后一个e_wo_entrp。<i name="e_wo_entrp"> -340.22390268</i>，而官方例子e-v.dat中的值不同。

我们知道在用phonopy声子谱计算的时候会进行扩胞，从解压得到的vasprun.xml文件可知例子中使用了2x2x2的扩胞。

对于单个POSCAR-5来说，对应的能量是对应值的1/8。-340.22390268除以8等于-42.5279879和官方-42.527932非常接近。可能是计算的时候保留的有效位不同造成的。同理试算POSCAR--5的对应值-337.05806784/8 = -42.1322585和-42.132246非常接近。

注意这里官方文档的POSCAR--5和对应的thermal_properties.yaml都是8个原子，而vasprun.xml是扩胞后的数据。

写完e-v.dat文件后输入

phonopy-qha e-v.dat thermal_properties-{1..10}.yaml

便可以得到对应的热膨胀系数文件thermal_expansion.dat和等压热容文件Cp-temperature_polyfit.dat。