Phono3py是一个主要用python写的声子-声子相互作用相关性质的模拟包,可以基于有限位移算法实现三阶力常数和晶格热导率的计算过程,同时输出包括声速,格林奈森常数,声子寿命和累积晶格热导率等参量。
相关介绍和安装请参考往期推荐。
phono3py快速安装教程
QE+phono3py 热导率计算
VASP+phono3py:快速计算晶格热导率
理论到实践:VASP+Phono3py计算Phonon Lifetime
phonopy 和 phono3py 安装教程Phonopy-Spectroscopy计算材料红外和Raman光谱
在phonopy和phono3py的计算过程中,会将过程文件和结果文件写入到hdf5的格式文件中,以此来方便程序读入和写入,phonopy有关的诸多输出量可以通过额外设置获得文本文件,并用于数据处理,phono3py则不能自主输出文本数据。但是在phono3py介绍网址上有如何读取hdf5文件的说明。
http://phonopy.github.io/phono3py/hdf5_howto.html
本文则以此来介绍读取phono3py计算完晶格热导率后的数据读出。
生成有限位移文件后,并通过VASP计算获得力常数fc3.hdf5和fc2.hdf5后,进行热导率计算,采用11×11×11的网格,
phono3py-load --mesh 11 11 11 --br
有关计算结果会被写入到 kappa-m111111.hdf5中,首先通过pip安装 ipython
pip install ipython
随后打开ipython编辑器
ipython

随后输入(In [1]:是ipython的输入提示,只需要输入冒号后面的部分就行)
In [1]: import h5pyIn [2]: f = h5py.File("kappa-m111111.hdf5")In [3]: list(f)
则输出kappa-m111111.hdf5中包含的数据的title
Out[3]:['frequency','gamma','group_velocity','gv_by_gv','heat_capacity','kappa','kappa_unit_conversion','mesh','mode_kappa','qpoint','temperature','weight']
获得热导率张量
In [4]: f['kappa'].shapeOut[4]: (101, 6)
In [5]: f['kappa'][:]Out[5]:array([[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,0.00000000e+00, 0.00000000e+00, 0.00000000e+00],[ 2.11702476e+05, 2.11702476e+05, 2.11702476e+05,6.64531043e-13, 6.92618921e-13, -1.34727352e-12],[ 3.85304024e+04, 3.85304024e+04, 3.85304024e+04,3.52531412e-13, 3.72706406e-13, -7.07290889e-13],...,[ 2.95769356e+01, 2.95769356e+01, 2.95769356e+01,3.01803322e-16, 3.21661793e-16, -6.05271364e-16],[ 2.92709650e+01, 2.92709650e+01, 2.92709650e+01,2.98674274e-16, 3.18330655e-16, -5.98999091e-16],[ 2.89713297e+01, 2.89713297e+01, 2.89713297e+01,2.95610215e-16, 3.15068595e-16, -5.92857003e-16]])
获得温度的张量
In [6]: f['temperature'][:]Out[6]:array([ 0., 10., 20., 30., 40., 50., 60., 70.,80., 90., 100., 110., 120., 130., 140., 150.,160., 170., 180., 190., 200., 210., 220., 230.,240., 250., 260., 270., 280., 290., 300., 310.,320., 330., 340., 350., 360., 370., 380., 390.,400., 410., 420., 430., 440., 450., 460., 470.,480., 490., 500., 510., 520., 530., 540., 550.,560., 570., 580., 590., 600., 610., 620., 630.,640., 650., 660., 670., 680., 690., 700., 710.,720., 730., 740., 750., 760., 770., 780., 790.,800., 810., 820., 830., 840., 850., 860., 870.,880., 890., 900., 910., 920., 930., 940., 950.,960., 970., 980., 990., 1000.])
获得特定温度热导率和模式热导率(需提前定义)
In [7]: f['kappa'][30]Out[7]:array([ 1.09089896e+02, 1.09089896e+02, 1.09089896e+02,1.12480528e-15, 1.19318349e-15, -2.25126057e-15])In [8]: f['mode_kappa'][30, :, :, :].sum(axis=0).sum(axis=0) / weight.sum()Out[8]:array([ 1.09089896e+02, 1.09089896e+02, 1.09089896e+02,1.12480528e-15, 1.19318349e-15, -2.25126057e-15])
计算声子寿命
In [9]: g = f['gamma'][30]In [10]: import numpy as npIn [11]: g = np.where(g > 0, g, -1)In [12]: lifetime = np.where(g > 0, 1.0 / (2 * 2 * np.pi * g), 0)
声子群速和热容等数据处理方式同理。
具体修改请查看h5py使用方法和hdf5格式文件读取以及ipython命令使用方法。