天玑智研_学术论坛_如何显示振动模？学习这篇文章就够了

今天我们将更新一篇用VESTA, VMD, v_sim画振动模的文章。首先我们会介绍如何让振动模在VESTA中显示，并让它在VMD中动起来。这里将用到一个脚本，用于提取DFPT计算得到的Gamma点振动模。

最后我们会介绍如何让整个布里渊区中任意k点的振动模显示出来，并且也让它动起来（在VMD或v_sim中）。

步入正题，我们需要准备四个文件，POSCAR，OUTCAR，poscar.vesta，还有就是本文提到的脚本modes_to_vesta.py。

该脚本来自于github，下载方式见文尾：

https://github.com/Stanford-MCTG/VASP-plot-modes

POSCAR ：用于计算的VASP结构文件。

OUTCAR：用IBRION=7，8计算得到的DFPT OUTCAR文件。这里注意计算是还需要加入参数 NWRITE=3，用于写出“Eigenvectors after division by SQRT(mass)”。

poscar.vesta，是通过vesta打开POSCAR，另存为.vesta格式得到的文件。它是作为后面写出振动模的模板。

脚本安装和运行：

拷贝modes_to_vesta.py到bin文件夹下：

$ cp modes_to_vesta.py ~/bin

给执行权限：

$ chmod +x modes_to_vesta.py

得到python2或python3位置：

$ which python3/usr/bin/python3

打开脚本，修改第一行：

$ vi modes_to_vesta.py

第一行，输入python的位置。

#!/usr/bin/python3

运行脚本

$ modes_to_vesta.py poscar.vesta

原理：

以MoS2最高振动模为例，vasp计算OUTCAR中的本征矢和本征值如下：

1 f = 13.733531 THz 86.290322 2PiTHz 458.101281 cm-1 56.797341 meV X Y Z dx dy dz 1.591824 0.919040 7.110820 -0.000000 -0.000000 -0.064613 -0.000000 1.838080 5.548443 -0.000001 -0.000001 0.096326 -0.000000 1.838080 8.673147 0.000001 0.000001 0.097036

程序中取了每个原子的dx，dy和dz，然后乘以scaling_factor作为VESTA中该原子的矢量。

所以，458.10 cm-1 振动模给出的矢量就是：

1 -0.00000 -0.00000 -2.58452

2 -0.00004 -0.00004 3.85304

3 0.00004 0.00004 3.88144

注意：程序中给出的scaling_factor是40，可以修改。

振动模动起来

振动模动起来主要是用freqmov.py这个脚本。

freqmov.py：该脚本主要用于把频率计算的结果转化为动画。

该脚本来自于github，下载脚本见文尾：

https://github.com/Angel-Jia/VASP-script

用法：

(freqmov.py 按照上一节方法，给执行权限，放到bin文件夹下，另外需要拷贝VASP.py到bin文件夹下，在freqmov.py 和VASP.py 中加入# coding=utf-8，可以让python 读取中文不报错)

$ freqmov.py freqfile 30 0.6

首先要制作freqfile文件,其格式为（这里只显6行）：

30 f/i= 43.595023 THz 273.915607 2PiTHz 1454.173389 cm-1 180.294587 meV X Y Z dx dy dz -1.334507 2.369499 6.388597 -0.015951 0.008861 -0.024862 2.621978 4.667179 6.419931 -0.008397 -0.017126 -0.024545 1.361682 2.515564 6.433748 0.010846 0.019286 -0.030780 2.704039 0.032985 6.442084 0.010923 -0.007642 -0.018801

其中：1、2行必须有，之后是每个原子的坐标及其振动方向，原子的数量必须与本体系原子数量一致。例如30个原子的体系，就必须有32行数据，至于是不是虚频则无所谓。

其中，

30是动画的帧数，数字越大动画越细腻，30足够了；

0.6是系数，用于调节振幅。OUTCAR中给出的振幅很大，需要适当缩减，0.6即可；OUTCAR可以采用IBRION = 5, 6 或IBRION = 7, 8计算得到。

最后会产生一个freqfile.xyz文件，用VMD打开即可播放。

freqfile文件怎么得到呢？freqfile文件可用下列方法提取：

$ grep "30 f/i" OUTCAR -A 5 > freqfile

"30 f/i" 为字符串30 f/i，有时候可能是"30 f"，依据不同振动模而定，具体见OUTCAR。

-A 5表示输出grep查找字符串所在行信息（这里是“30 f/i”）以及其后5行。也可以自己打开OUTCAR，找到所需信息直接存为文本格式。

最后下载VMD导入xyz文件即可看到振动模了。

VMD => Graphics => Representations => Drawing Method => CPK

任意k点的振动模

这里，我们将介绍两种方法，来实现任意k点下振动模动起来的目标。

1. 老方法

首先是我们之前讲到的，绘制振动模实际上是把本征矢量乘以一个系数得到的一个矢量。那么我们想绘制哪个点，只要知道他的本征矢量就可以。

在phonopy中我们可以通过BAND选项来设置任意k点，然后通过EIGENVECTORS=TRUE.这个参数来实现显示这些点的本征矢量。这样我们就可以按照之前的方法来绘制振动模了。

具体设置位置可以放在band.conf文件中。

具体本征矢的位置在band.yaml里,格式为:

- q-position: [ 0.0000000, 0.0000000, 0.0000000 ]

distance: 0.0000000

band:

- # 1

frequency: -0.0370089502

eigenvector:

- # atom 1

- [ -0.00000000000000, 0.00000000000000 ]

- [ 0.00000000000000, 0.00000000000000 ]

- [ -0.63347564223346, 0.00000000000000 ]

- # atom 2

- [ -0.00000000000000, 0.00000000000000 ]

- [ 0.00000000000000, 0.00000000000000 ]

- [ -0.77376263201120, 0.00000000000000 ]

第一行为q点，

第二行是距离第一个点的距离,用于绘制能带时的横轴，

第三行，band标志，

第四行1，代表第一个能带，

第五行为该能带的频率值，

第六行为本征矢名称，

第七行表示第一个原子，

第八行到第十行是本征矢的矩阵。这里写的是实部加虚部的形式。一般只需要实部部分就好。

根据这些信息自己写一个之前给的freqmov.py 的输入文件，即可让振动模动起来了。也可以直接写一个OUTCAR类似的输出文件。让modes_to_vesta.py来读取。代码能力强的话可以直接写一个导出VESTA格式的脚本出来，就能用VESTA打开了。

2. 新方法

phonopy里有这样一个参数ANIME，它可以用于写出动画格式文件。（在band.conf加入该参数）

ANIME = 0.5 0.5 0

让ANIME等于某个q点，就会写出这个q点下所有振动模的信息。写出的文件为anime.ascii。

我们只要安装一个可视化软件就可以了。phonopy推荐的可视化软件为v_sim。

这里用的是ubuntu的虚拟机。一条命令就可以安装。

sudo apt install v-sim

下载安装包自己编译可以访问如下教程：

https://blog.csdn.net/CSDN296218029/article/details/101851656

结果：运行命令

$ v_sim anime.ascii

选择如下窗口：

undefined

选择振动模，然后点击右下角绿色三角形：

动图如下：

v_sim 是一个强大的软件，大家研究研究还能玩出新花样。快快跟着教程练习吧。

如何显示振动模？学习这篇文章就够了

评论