DynaPhoPy计算高温非谐声子谱

前言

在诸多关于材料的理论计算中为了探究材料的热学等性质会计算该材料的声子谱，为了证明这个材料在理论上具有稳定性，会通过各种方式获得一个没有虚频的声子谱。根据笔者和诸多计算材料的研究者的经验，部分结构的材料，尤其是根据实验结果测试获得的材料，往往在低温下使用简谐近似计算得到的声子谱存在一个很难去除虚频的声子谱，但是带有温度的考虑非谐效应的声子谱时则会得到一个无虚频的结果。同时，为了探究材料在高温时所具有的声子运输性质，研究非谐性效应，对于理论上和应用上都具有非常大的意义。

本篇文章主要介绍使用VASP，phonopy，DynaPhoPy软件计算考虑非简谐效应的声子谱，并和简写近似的声子谱进行对比。

参考链接：

有限温度声子谱之Dynaphopy实现 - 哔哩哔哩

Dynaphopy的使用 - 知乎

GitHub - wenlibin02/DynaPhoPy: Phonon Analysis from Molecular Dynamics Data （官网）

非谐声子模拟方法(Normal Mode Decomposition) - 知乎

Dynaphopy     （使用说明）

计算流程

首先我们将需要计算的模型充分弛豫进行结构优化，并通过phonopy扩包，使用SPOSCAR计算得到其声子谱。

这个过程中主要获得0K时的声子谱数据和力常数文件FORCE_CONSTANTS（DFPT方法，如果使用有限位移法计算的声子谱，得到的力常数文件为FORCE_SETS，可通过DynaPhoPy作者所写的脚本转换为FORCE_CONSTANTS，脚本在DynaPhoPy软件安装包中的example文件夹内）。

使用计算声子谱时所扩包得到的SPOSCAR进行第一性分子动力学模拟（AIMD），主要为得到OUTCAR和XDATCAR。

这里我的设置为300K，INCAR参数如下：

Global Parameters

  ISTART =  0            (Read existing wavefunction; if there)

  # ISPIN =  2           (Spin polarised DFT)

  # ICHARG =  11         (Non-self-consistent: GGA/LDA band structures)

  LREAL  = Auto          (Projection operators: automatic)

  ENCUT  =  600       (Cut-off energy for plane wave basis set, in eV)

  PREC   =  Normal       (Precision level)

  LWAVE  = .FALSE.        (Write WAVECAR or not)

  LCHARG = .FALSE.        (Write CHGCAR or not)

  ADDGRID= .TRUE.        (Increase grid; helps GGA convergence)

Electronic Relaxation

  ISMEAR =  -1            (fermi smearing for considering electron temperature)

  SIGMA  =  0.02585201664 (electron temperature setting kb*T)

  EDIFF  =  1E-05

  IALGO = 38

  NELMIN = 4

Molecular Dynamics

  IBRION =  0            (Activate MD)

  NSW    =  20000         (Max ionic steps)

  EDIFFG = -1E-03        (Ionic convergence; eV/A)

  POTIM  =  1            (Timestep in fs)

  SMASS  =  0            (MD Algorithm: -3-microcanonical ensemble; 0-canonical ensemble)

  TEBEG  =  300     (Start temperature K)

  TEEND  =  300   (Start temperature K)

  MDALGO =  2         (Andersen Thermostat)

  ISYM   =  0          (Symmetry: 0=none; 2=GGA; 3=hybrids)

  ISIF = 2

这里尤其要注意NSW的设置，尽量设置的大一点，依据笔者经验，NSW/1000的数值应该大于原胞的总原子数，具体取多少可以自行测试。

然后设置一个input文件，这个文件名可以自定义，主要内容如下：

STRUCTURE FILE POSCAR

POSCAR-afterlr #结构优化后的POSCAR

FORCE CONSTANTS

FORCE_CONSTANTS_443 #力常数文件

PRIMITIVE MATRIX #原胞与晶胞的基矢转换矩阵

1.0 0.0 0.0

0.0 1.0 0.0

0.0 0.0 1.0

SUPERCELL MATRIX #扩胞系数

4 0 0

0 4 0

0 0 3

MESH PHONOPY #网格 默认40

40 40 40

BANDS #高对称点

  0.3333333333   0.3333333333   0.0000000000       0.0000000000   0.0000000000   0.0000000000

  0.0000000000   0.0000000000   0.0000000000        0.5000000000   0.0000000000   0.0000000000

   0.5000000000   0.0000000000   0.0000000000        0.3333333333   0.3333333333   0.0000000000

   0.3333333333   0.3333333333   0.0000000000      0.3333333333   0.3333333333   0.5000000000

  0.3333333333   0.3333333333   0.5000000000       0.0000000000   0.0000000000   0.5000000000

   0.0000000000   0.0000000000   0.5000000000     0.5000000000   0.0000000000   0.5000000000

   0.5000000000   0.0000000000   0.5000000000     0.3333333333   0.3333333333   0.5000000000

如果你可以跳过终端连接软件（比如Xshell等）直接操作服务器或电脑主机，意思就是可以看到窗口弹出的图像显示，那么可以使用下面这个命令进如操作窗口，并获得程序弹出的数据和图像。

dynaphopy INPUT OUTCAR -i

这时我们为了获得考虑了非谐性的声子谱的数据，需要执行以下命令：

dynaphopy input_file TRAJECTORY -sfc FORCE_CONSTANTS

input_file就是刚才我们所写的input文件的名称，TRAJECTORY文件在VASP计算中则使用OUTCAR或XDATCAR替换，在lammps的计算中则使用其轨迹文件。此时，则写入了考虑非谐性计算的力常数FORCE_CONSTANTS。（可自定义文件名称，但建议使用这个）

然后我们使用phonopy依据这个导出数据并绘图。计算声子谱图像和数据导出的命令因个人习惯而定。