VASP_phonopy学习笔记

phonopy课程

课程的主要内容:

  1. phonopy软件入门
  2. 扩散系数前置因子的计算
  3. 不同温度下热力学性质计算(自由能、热容、振动熵)
  4. 温度驱动的相变预测
  5. 材料的热膨胀性能预测
  6. 拉曼/红外光谱计算与模拟

1. phonopy软件入门

1.1 phonopy软件的安装

推荐使用conda安装

  1. 下载anaconda并且安装

  2. 一条命令安装phonopy

    conda install -c conda-forge phonopy
    conda install hdf5=1.8.18   (可选)
  3. 如果想新建环境,可以使用下面的几条命令

    conda create --name phonopy python=3
    conda activate phonopy
    conda install -c conda-forge phonopy h5py
    conda deactivate

1.2 phonopy计算声子谱的小例子

1. 首先需要对原胞进行高精度的结构优化

INCAR_OPT

ISTART = 1
LREAL = .FALSE.
PREC = Accurate
LWAVE = .False.
LCHARG = .False.
ADDGRID= .TRUE.
ENCUT = 400

ISMEAR = 0
SIGMA = 0.05
NELM = 60
NELMIN = 6
EDIFF = 1E-08  # 甚至为1E-12

NSW = 100
IBRION = 2
ISIF = 3
EDIFFG = -0.001  # 甚至更精确

KPOINTS_OPT

KPOINTS
0
G
7 7 7
0 0 0

PORCAR_OPT

GaAs
1.0
        5.6532001495         0.0000000000         0.0000000000
        0.0000000000         5.6532001495         0.0000000000
        0.0000000000         0.0000000000         5.6532001495
   As   Ga
    4    4
Direct
     0.250000000         0.250000000         0.250000000
     0.750000000         0.750000000         0.250000000
     0.750000000         0.250000000         0.750000000
     0.250000000         0.750000000         0.750000000
     0.000000000         0.000000000         0.000000000
     0.000000000         0.500000000         0.500000000
     0.500000000         0.000000000         0.500000000
     0.500000000         0.500000000         0.000000000

优化完之后的结构

GaAs                                    
   1.00000000000000     
     5.7447626994798862    0.0000000000000000    0.0000000000000000
     0.0000000000000000    5.7447626994798862    0.0000000000000000
    -0.0000000000000000    0.0000000000000000    5.7447626994798862
   As   Ga
     4     4
Direct
  0.2500000000000000  0.2500000000000000  0.2500000000000000
  0.7500000000000000  0.7500000000000000  0.2500000000000000
  0.7500000000000000  0.2500000000000000  0.7500000000000000
  0.2500000000000000  0.7500000000000000  0.7500000000000000
 -0.0000000000000000 -0.0000000000000000  0.0000000000000000
  0.0000000000000000  0.5000000000000000  0.5000000000000000
  0.5000000000000000  0.0000000000000000  0.5000000000000000
  0.5000000000000000  0.5000000000000000 -0.0000000000000000

2. 其次通过phonopy对优化完的结构进行扩胞

新建一个文件夹phonon_cal,将高精度结构优化中的CONTCAR拷贝到该文件夹中,并且重命名为POSCAR_UNIT,并且使用phonopy进行扩胞

phonopy -d --dim="2 2 2" -c POSCAR_UNIT

扩胞之后生成的文件

(phonopy) [sunzhenlu@m01 phonon_cal]$ ll
total 44
phonopy_disp.yaml   # 位置偏移信息,还挺详细的
POSCAR-001          # 在POSCAR_UNIT的基础上增加了偏移
POSCAR-002
POSCAR_UNIT         # 原有的高精度优化的结构文件
SPOSCAR             # 将各偏移合并到一起

根据计算使用的不同文件,可以分为两种方法

  • DFPT方法(直接使用SPOSCAR进行计算)
  • Finite displacement方法(将POSCAR-***拷贝到分别的目录中做自洽计算)

官方说明:DFPT方法是更为准确的,但是在很多情况下,比如原胞的原子已经很多,DFPT方法是不现实的,这时候只能采用Finite displacement方法。

实际上两者的计算量是基本相同的。

Finite displacement方法的好处是可以多任务并行。

3.1 DFPT方法

INCAR_DFPT

ENCUT = 400
NELM = 100
ALGO = Normal
PREC = Accurate
ISMEAR = 0
SIGMA = 0.05
EDIFF = 1E-8
AMIX = 0.1
BMIX = 0.01

IBRION = 8
NSW = 1

band.conf

ATOM_NAME = As Ga
DIM = 2 2 2
FORCE_CONSTANTS = READ
BAND = 0.5 0.0 0.0  0.5 0.5 0.5  0.5 0.5 0.0  0.0 0.0 0.0  0.5 0.5 0.5

第四行是高对称点,可以通过下面方法得到:

将原胞文件拖到Materials Studio中,然后顶栏Tools -> Brillouin Zone Path -> create


之后就是提交计算任务

如果因为晶格的原因报错了,应该在声子谱计算之前,利用phonopy重新检查晶格对称性

phonopy --symmetry

计算完成之后进行下面的处理:

  1. 收集力常数
phonopy --fc vasprun.xml
  1. 得到band.yaml文件
phonopy --dim="2 2 2" -c POSCAR-UNIT -p band.conf
  1. 转格式
bandplot --gnuplot band.yaml > band.dat
  1. band.dat文件导入origin中画图
grep - band.dat # 检查是否存在虚频

3.2 Finite displacement方法

INCAR_Finite_displacement

ENCUT = 400
NELM = 100
ALGO = Normal
PREC = Accurate
ISMEAR = 0
SIGMA = 0.05
EDIFF = 1E-8    #精度同样要保证
AMIX = 0.1
BMIX = 0.01

band.conf

ATOM_NAME = As Ga
DIM = 2 2 2
# FORCE_CONSTANTS = READ
BAND = 0.5 0.0 0.0  0.5 0.5 0.5  0.5 0.5 0.0  0.0 0.0 0.0  0.5 0.5 0.5

提交计算任务,针对每一个POSCAR-***做一个静态计算

(phonopy) [sunzhenlu@m01 finite]$ tree
.
|-- 001
|   |-- INCAR
|   |-- KPOINTS
|   |-- POSCAR      # 原来的POSCAR-001
|   |-- POTCAR
|   `-- vasp.pbs
|-- 002
|   |-- INCAR
|   |-- KPOINTS
|   |-- POSCAR      # 原来的POSCAR-001
|   |-- POTCAR
|   `-- vasp.pbs
|-- phonopy_disp.yaml
|-- POSCAR-001
|-- POSCAR-002
|-- POSCAR_UNIT
`-- SPOSCAR

计算完成之后进行下面的处理:

  1. 收集力常数

    将每个子文件中计算得到的vasprun.xml文件,拷贝到上一级目录并且重命名为vasprun.xml-***

phonopy -f vasprun.xml-*
  1. 得到band.yaml文件
phonopy --dim="2 2 2" -c POSCAR-UNIT -p band.conf
  1. 转格式
bandplot --gnuplot band.yaml > band.dat
  1. band.dat文件导入origin中画图
grep - band.dat # 检查是否存在虚频

2. 扩散系数前置因子的计算

在DFPT或者Finite displacement方法计算声子谱的基础上,得到mesh.yaml文件

  1. 新建一个文件夹,将FORCE_CONSTANTSPOSCARmesh.conf放入该文件夹中(POSCAR为原胞)

  2. 得到mesh.yaml文件

    phonopy mesh.conf -c POSCAR
    
    mesh.conf的内容
    ATOM_NAME = Hf Si O
    DIM = 2 2 2
    MP = 13 13 13
    FORCE_CONSTANTS = READ
  3. mesh.yaml中找到gamma点的频率

    phonon:
    - q-position: [    0.0000000,    0.0000000,    0.0000000 ]
     distance_from_gamma:  0.000000000
     weight: 1    
     band:
     - # 1
       frequency:    -0.0036593413
     - # 2
       frequency:     0.0036958358
     - # 3
       frequency:     0.0036958358
     - # 4
       frequency:     4.0238095850
     - # 5
       frequency:     4.5901171620
     - # 6
       frequency:     4.8708104657
    .
    .
    .
     - # 36
       frequency:    30.4419382088
    
  4. 利用公式求解

D=D_{0} \exp \left(-\frac{E_{a}}{k_{\mathrm{B}} T}\right)=g a^{2} \tilde{v} \exp \left(-\frac{E_{a}}{k_{\mathrm{B}} T}\right)
\tilde{v}=\frac{\prod_{i=1}^{3 N-3} v_{i}}{\prod_{i=1}^{3 N-4} v_{i}^{\prime}}

一篇参考文献:On the vineyard formula for the pre-exponential factor in the Arrhenius law

where νi are the eigenfrequencies of vibrations of the cluster in the state corresponding to the minimum potential energy Epot for all 3N – 6 normal coordi nates, and ν‘i are the frequencies of vibrations at the saddle point corresponding to the maximum of Epot for one normal coordinate and to the minimum for all the other coordinates (since one of the 3N – 6 frequencies is imaginary, it is not included in the denominator of expression (2) for A). In the simulation of a macro scopic sample by an N atom supercell with periodic boundary conditions, the numbers of nonzero real fre quencies νi and in expression (2) are equal to 3N – 3 and 3N – 4, respectively.

其中νi是在与所有3N – 6正态坐标的最小势能Epot相对应的状态下,团簇的振动的本征频率,ν'i是在鞍点处的振动频率,对应于一个点的Epot最大值 法向坐标,并且对于所有其他坐标都应为最小值(由于3N – 6频率之一是虚数,因此不包含在A的表达式(2)的分母中)。 在一个具有周期性边界条件的N原子超晶胞对宏观样品的模拟中,非零实际频率νi和表达式(2)的数量分别等于3N – 3和3N – 4。

不是我的专业,如果公式错误,请联系我

3. 不同温度下热力学性质计算(自由能、热容、振动熵)

  • 简谐振子的本征值
E_{n}=\left(\frac{1}{2}+n\right) \hbar \omega_{k} \quad n=0,1,2,3 \ldots
  • 声子内能
E=\sum_{\mathbf{q} \nu} \hbar \omega(\mathbf{q} \nu)\left[\frac{1}{2}+\frac{1}{\exp \left(\hbar \omega(\mathbf{q} \nu) / k_{\mathrm{B}} T\right)-1}\right]
  • 等容热容
C_{V} =\left(\frac{\partial E}{\partial T}\right)_{V}=\sum_{\mathbf{q} \nu} k_{\mathrm{B}}\left(\frac{\hbar \omega(\mathbf{q} \nu)}{k_{\mathrm{B}} T}\right)^{2} \frac{\exp \left(\hbar \omega(\mathbf{q} \nu) / k_{\mathrm{B}} T\right)}{\left[\exp \left(\hbar \omega(\mathbf{q} \nu) / k_{\mathrm{B}} T\right)-1\right]^{2}}
  • 亥姆霍兹自由能
F =-k_{\mathrm{B}} T \ln Z =\varphi+\frac{1}{2} \sum_{\mathbf{q} \nu} \hbar \omega(\mathbf{q} \nu)+k_{\mathrm{B}} T \sum_{\mathbf{q} \nu} \ln \left[1-\exp \left(-\hbar \omega(\mathbf{q} \nu) / k_{\mathrm{B}} T\right)\right]
  • 振动熵
S =-\frac{\partial F}{\partial T}=\frac{1}{2 T} \sum_{\mathbf{q} \nu} \hbar \omega(\mathbf{q} \nu) \operatorname{coth}\left(\hbar \omega(\mathbf{q} \nu) / 2 k_{\mathrm{B}} T\right)-k_{\mathrm{B}} \sum_{\mathbf{q} \nu} \ln \left[2 \sinh \left(\hbar \omega(\mathbf{q} \nu) / 2 k_{\mathrm{B}} T\right)\right]

在DFPT或者Finite displacement方法计算声子谱的基础上,得到mesh.yaml文件

  1. 新建一个文件夹,将FORCE_CONSTANTSPOSCARmesh.conf放入该文件夹中(POSCAR为原胞)

  2. 得到mesh.yaml文件

    phonopy mesh.conf -c POSCAR
    
    mesh.conf的内容
    ATOM_NAME = Hf Si O
    DIM = 2 2 2
    MP = 13 13 13
    FORCE_CONSTANTS = READ
  3. 得到thermal_properties.yaml文件

    phonopy -c POSCAR mesh.conf  -t
  4. 得到thermal_properties.dat文件

    phonopy-propplot --gnuplot thermal_properties.yaml >thermal_properties.dat
  5. 将dat文件导入到Origin中画图

dat文件的5列数据分别为:
temperature[K]
Free energy [kJ/mol]
Entropy [J/K/mol]
Heat capacity [J/K/mol]
Energy [kJ/mol]

4. 温度驱动的相变预测

我也听不懂,抱歉

5. 材料的热膨胀性能预测

分为一下几个步骤

  1. 高精度的原胞优化
  2. 在此基础上施加形变(梯度可以设置为千分之一)
  3. 优化多个结构
  4. 使用DFPT方法计算每个结构的声子谱(如果结构有问题,使用phonopy --symmetry
  5. 在每个文件夹中得到thermal_properties.yaml
  6. 运行脚本得到e-v.dat
  7. 运行phonopy-qha的相关命令得到各种dat文件
  8. Over

下面是能用到的各种命令

扩胞用
phonopy -d --dim="2 2 2" -c POSCAR_UNIT
查找晶格对称性
phonopy --symmetry
收集力常数
phonopy --fc vasprun.xml
得到band.yalm文件(好像不需要这一步)
phonopy --dim="2 2 2" -c POSCAR-UNIT -p band.conf
得到thermal_properties.yaml文件
phonopy -c POSCAR mesh.conf  -t

下面是用到的两个conf文件的格式,需要根据自己的体系进行修改

band.conf

ATOM_NAME = Si
DIM = 3 3 3
FORCE_CONSTANTS = READ
BAND = 0.0 0.0 0.0  0.0 0.5 0.0  0.0 0.5 0.5 0.0 0.0 0.5 0.0 0.0 0.0

mesh.conf

ATOM_NAME = Si
DIM = 3 3 3
MP = 13 13 13
FORCE_CONSTANTS = READ

下面是产生e-v.dat文件的脚本,需要根据自己的目录名进行修改

for i in -0.005 -0.004 -0.003 -0.002 -0.001 0.001 0.002 0.003 0.004 0.005
do
 cd AC_Uni$i;
 pwd;
 volume=`grep volume vasprun.xml | tail -1| awk '{printf("%20.13f", $3)}'`;
 energy=`grep 'E0' OSZICAR | tail -1 | awk '{printf "%12.4f \t", $5}'`;
 echo ${volume} ${energy} >> ../e-v.dat;
 cd ../; #Exit the Strain folder
done

生成文件的名称和含义

(phonopy) [sunzhenlu@m01 file]$ ll
total 272
-rw-rw-r-- 1 sunzhenlu sunzhenlu  4700 Jul 30 16:53 bulk_modulus-temperature.dat
-rw-rw-r-- 1 sunzhenlu sunzhenlu  4200 Jul 30 16:53 Cp-temperature.dat
-rw-rw-r-- 1 sunzhenlu sunzhenlu  4200 Jul 30 16:53 Cp-temperature_polyfit.dat
-rw-rw-r-- 1 sunzhenlu sunzhenlu 35744 Jul 30 16:53 Cv-volume.dat
-rw-rw-r-- 1 sunzhenlu sunzhenlu  4200 Jul 30 16:53 dsdv-temperature.dat
-rw-rw-r-- 1 sunzhenlu sunzhenlu 35623 Jul 30 16:53 entropy-volume.dat
-rw-rw-r-- 1 sunzhenlu sunzhenlu  4700 Jul 30 16:53 gibbs-temperature.dat
-rw-rw-r-- 1 sunzhenlu sunzhenlu  4700 Jul 30 16:53 gruneisen-temperature.dat
-rw-rw-r-- 1 sunzhenlu sunzhenlu 32037 Jul 30 16:53 helmholtz-volume.dat
-rw-rw-r-- 1 sunzhenlu sunzhenlu  5200 Jul 30 16:53 volume-temperature.dat

含义:嗯~ o( ̄▽ ̄)o

官网有

6. 拉曼/红外光谱计算与模拟

尴尬\~决定不学了~

暂无评论

发送评论 编辑评论


				
上一篇
下一篇