频率分析

Northword2021年1月21日
  • VASP
  • 反应路径
大约 9 分钟

频率分析

Todo.... 频率分析是个啥,理论知识...

频率分析作用

  1. 确定结构是否稳定;

  2. 看振动方式和大小,用来和实验对比,棋博士最新的文章就是一个非常好的例子;

  3. 反应热,反应能垒,吸附能等的零点能矫正;

  4. 确认过渡态(有一个振动的虚频)

  5. 热力学中计算 entropy,用于计算化学势,微观动力学中的指前因子和反应能垒。

步骤

结构优化

在常规结构优化基础上进行下一步。

待解决的问题

乙醇结构优化中,若指定 EDIFF=1E-4(第一次)或 1E-6(第二次),EDIFFG=-2E-2,POTIM 默认(0.015),计算无法收敛(每个离子步都算满了 60 个电子步,19-20 个离子步后报错),提示如下:

ZBRENT: fatal error in bracketing
     please rerun with smaller EDIFF, or copy CONTCAR
     to POSCAR and continue

但(第三次)EDIFF 和 EDIFFG 默认,POTIM=0.05 时可在 7 步收敛(每个离子步仍是 60 电子步)。(此时查 OUTCAR 有 EDIFF=0.1E-3,EDIFFG=0.1E-2)(POTIM 默认 0.5)

初始结构为 www.chemspider.com 下载,20 20 20 的 cell,K 点 gamma 111。原因待测试。(其余参数 ISMEAR=0,SIGMA=0.01,IBRION=2,NSW=100,未给出均为默认)(测试一下①EDIFF=1E-4,EDIFFG=-2E-2,POTIM=0.05;②EDIFF=0.1E-3,EDIFFG=-0.1E-2,POTIM 默认)

频率计算

IBRION = 5       # Use 5 for Freq calculation
NSW    = 1
NFREE  = 2       # Do not use NFREE=1
POTIM  = 0.02
EDIFF  = 1E-6

# NCORE = 4   # comment this line
  • IBRION 的值改成 5
  • POTIM 用一个更小的值,我们这里用的 0.02,默认值是 0.015
  • NSW 设置成 1,这个可以直接不管,继续采用优化时的 NSW 值,因为你设置成 1, 2, 3, 4, 5, …, 1000 都不会影响计算;但不能不设置(因为默认值是 0,这时算个单点后任务便停止了。)
  • NFREE=2 添加这一个参数,表明原子在某一方向上正反两个方向移动;
  • 此外,EDIFF 也要设置一个严格的值(频率计算时,默认值为 1E-6,足够了!下一节会讲到)

结果分析

步数

当设置了 NFREE=2 且所有原子弛豫的时候,频率计算需要 $6N+1$ 步。N 为体系中的振动的原子数,这是因为:

  1. 第一个离子步是个频率计算前的单点计算。

  2. N 个原子,每个原子在 x、y、z 三个方向均有一个自由度,共 3N。

  3. 设置 NFREE=2,也就是在每个方向上 +POTIM–POTIM 都移动并算一下,就有了$3N \times 2 = 6N$步。

    官网原文如下,还要查阅 IBRIONNFREE 的相关内容。

The parameter NFREE determines how many displacements are used for each direction and ion, and POTIM determines the step size. The step size is defaulted to 0.015 ? (starting from VASP.5.1), if too large values are supplied in the input file. Expertise shows that this is a very reasonable compromise.

NFREE=2 uses central differences, i.e., each ion is displaced by a small positive and negative displacement, ±POTIM, along each of the cartesian directions.`

例如,乙醇分子$CH_3CH_2OH$,含有 9 个原子,其振动频率计算应有 55 步。 Ex24 乙醇分子振动频率计算(二) | LVTHWopen in new window

这一过程在 stdout 里也有较为明显的表示:

   1 F= -.10036430E+02 E0= -.10036285E+02  d E =-.289628E-03
 Finite differences POTIM= 0.02000 DOF=  27
 bond charge predicted

   2 F= -.78734041E+01 E0= -.78734041E+01  d E =-.373678E-15
 Finite differences progress:
  Degree of freedom:   1/ 27
  Displacement:        1/  2
  Total:               1/ 54
 bond charge predicted

   3 F= -.67069872E+01 E0= -.67026196E+01  d E =-.873513E-02
 Finite differences progress:
  Degree of freedom:   1/ 27
  Displacement:        2/  2
  Total:               2/ 54
 bond charge predicted

   4 F= -.67462590E+01 E0= -.67409236E+01  d E =-.106707E-01
 Finite differences progress:
  Degree of freedom:   2/ 27
  Displacement:        1/  2
  Total:               3/ 54
 bond charge predicted

······

  55 F= -.98834544E+01 E0= -.98834523E+01  d E =-.431696E-05
 Finite differences progress:
  Degree of freedom:  27/ 27
  Displacement:        2/  2
  Total:              54/ 54
 Finite differences POTIM=  2.000000000000000E-002

振动频率可视化

使用 p4vasp 或 jmol。 Ex25 乙醇分子振动频率计算(三) | LVTHWopen in new window

OUTCAR 中的信息

 Finite differences progress:
  Degree of freedom:  27/ 27
  Displacement:        2/  2
  Total:              54/ 54

 SECOND DERIVATIVES (NOT SYMMETRIZED)
 ------------------------------------
               1X          1Y          1Z          2X          2Y          2Z          3X          3Y          3Z          4X          4Y          4Z          5X          5Y          5Z          6X          6Y          6Z          7X          7Y          7Z          8X          8Y          8Z          9X          9Y          9Z
  1X    -0.796290   -0.233038    0.000000    1.493917   -0.390431    0.000000   11.997934    0.713060    0.000000   -0.502744   -0.458102   -1.112604   -0.502744   -0.458102    1.112604   -9.852689   -0.544908    0.000000    1.558071   -1.815756    2.790667    1.558071   -1.815756   -2.790667   -4.953526    5.003033    0.000000
  1Y     0.375968    0.109966    0.000000   -0.221500   -0.104078    0.000000   -7.444189    0.714797    0.000000    0.061864   -0.018602   -0.070447    0.061864   -0.018602    0.070447    6.078079   -0.495659    0.000000    0.526607   -0.086673    0.245169    0.526607   -0.086673   -0.245169    0.034699   -0.014475    0.000000
  1Z     5.808229   -2.434202   -0.224578    9.890712   -0.191510   -0.835513 -196.271299   22.155307   -2.396373    4.894997   -1.437502    1.327563    3.314884   -1.079152   -0.255396  159.821373  -14.859397    2.205595    2.708645   -0.577444   -1.520109    5.153556   -1.595313    1.097139    4.678904    0.019213    0.601672
  ······
  9Z    -2.239638   -0.224936    0.640338   -2.848485   -0.518046    0.053144   92.115496   -4.125652   -4.274775   -5.253060    0.192219    0.259185   -3.326232   -0.413466   -0.479808  -74.488855    6.180803    1.430131   -0.656228   -1.103483    0.253415   -0.780173    0.700478    3.745246   -2.522825   -0.687916   -1.626877

 
 Eigenvectors and eigenvalues of the dynamical matrix
 ----------------------------------------------------

   1 f  =  201.746767 THz  1267.612322 2PiTHz 6729.547573 cm-1   834.357861 meV
             X         Y         Z           dx          dy          dz
      8.658517  9.304797  0.000000    -0.012567    0.004216    0.089194
      0.120860  9.184513  0.000000    -0.000391    0.000152    0.006966
      9.405744 18.133372  0.000000     0.241400   -0.053508   -0.052395
     12.979810 18.660074 17.736563    -0.057237    0.039130    0.448592
     12.979810 18.660074  2.263437     0.019013    0.039514    0.091267
      8.203737 18.242837  0.000000    -0.670436   -0.088927    0.146993
     16.776380  8.675669  2.229820     0.035997   -0.046760   -0.002012
     16.776380  8.675669 17.770180     0.062323   -0.406063    0.024770
     14.098762 10.462995  0.000000    -0.025516    0.172042   -0.171822

   2 f  =   47.211040 THz   296.635710 2PiTHz 1574.790721 cm-1   195.249235 meV
             X         Y         Z           dx          dy          dz
      8.658517  9.304797  0.000000     0.042971    0.024710    0.035150
      ······
   27 f/i=  203.242065 THz  1277.007557 2PiTHz 6779.425348 cm-1   840.541919 meV
             X         Y         Z           dx          dy          dz
      8.658517  9.304797  0.000000    -0.002381   -0.002248   -0.090493
      ······

 Finite differences POTIM=  2.000000000000000E-002
  LATTYP: Found a simple cubic cell.
 ALAT       =    20.0000000000

频率相关的信息会被输出到 OUTCAR 的这两个部分,

第一部分:二阶导,没啥用

第二部分:特征值和特征向量,主要看这个

1 f 行(line20)是四个频率单位的数值。下面几行是每个原子的坐标(X、Y、Z)及其在 x y z 方向上的振动大小(dx、dy、dz),坐标是分数坐标系。

四个频率的换算:

$$ \begin{aligned} & E = hc/\lambda \ & \nu = c / \lambda \ & \tilde{\nu} = 1 / \lambda \ & T = 1 / \nu \ \end{aligned}\

\begin{aligned} \text{in which,} \ & E=\text{energy}\space (eV) \ & \lambda=\text{wavelength (m)} \ & \tilde{\lambda}\text{ = wavenumber }(m^{−1}) \ & T=\text{period (s)} \ & \nu=\text{frequency }(s^{−1}\space or\space Hz) \ & h=\text{Planck’s constant = } 4.135667516×10^{−15} ~eV \cdot s \ & c=\text{speed of light = 299792458 m/s} \ \end{aligned} $$

此外,$1 THz = 1012 Hz, \quad 1 cm^{-1} = 100 m^{-1}$ 。

还可以用 http://halas.rice.edu/conversionsopen in new window 在线转换单位。

频率提取:

[2020223055092@mu02 freq]$ grep cm-1 OUTCAR 
   1 f  =  201.746767 THz  1267.612322 2PiTHz 6729.547573 cm-1   834.357861 meV
   2 f  =   47.211040 THz   296.635710 2PiTHz 1574.790721 cm-1   195.249235 meV
   3 f  =   35.921110 THz   225.698994 2PiTHz 1198.199235 cm-1   148.557825 meV
   4 f  =   30.557648 THz   191.999365 2PiTHz 1019.293390 cm-1   126.376319 meV
   5 f  =   28.299918 THz   177.813630 2PiTHz  943.983631 cm-1   117.039096 meV
   6 f  =   24.737229 THz   155.428593 2PiTHz  825.145113 cm-1   102.304992 meV
   7 f  =   20.159900 THz   126.668391 2PiTHz  672.461876 cm-1    83.374677 meV
   8 f  =   17.283332 THz   108.594381 2PiTHz  576.509899 cm-1    71.478143 meV
   9 f  =   16.416363 THz   103.147049 2PiTHz  547.590902 cm-1    67.892643 meV
  10 f  =   12.378931 THz    77.779114 2PiTHz  412.916663 cm-1    51.195160 meV
  11 f  =    7.042735 THz    44.250808 2PiTHz  234.920339 cm-1    29.126420 meV
  12 f  =    6.004684 THz    37.728545 2PiTHz  200.294706 cm-1    24.833387 meV
  13 f  =    3.621816 THz    22.756539 2PiTHz  120.810763 cm-1    14.978631 meV
  14 f  =    1.485344 THz     9.332691 2PiTHz   49.545738 cm-1     6.142891 meV
  15 f/i=    0.608073 THz     3.820638 2PiTHz   20.283146 cm-1     2.514790 meV
  16 f/i=    2.581155 THz    16.217876 2PiTHz   86.098066 cm-1    10.674804 meV
  17 f/i=    4.872529 THz    30.615003 2PiTHz  162.530067 cm-1    20.151167 meV
  18 f/i=    6.118100 THz    38.441159 2PiTHz  204.077858 cm-1    25.302439 meV
  19 f/i=    8.804759 THz    55.321930 2PiTHz  293.695124 cm-1    36.413568 meV
  20 f/i=   10.508365 THz    66.026003 2PiTHz  350.521305 cm-1    43.459119 meV
  21 f/i=   15.745766 THz    98.933566 2PiTHz  525.222205 cm-1    65.119277 meV
  22 f/i=   18.917161 THz   118.860028 2PiTHz  631.008549 cm-1    78.235117 meV
  23 f/i=   21.091439 THz   132.521422 2PiTHz  703.534668 cm-1    87.227213 meV
  24 f/i=   24.556339 THz   154.292030 2PiTHz  819.111283 cm-1   101.556892 meV
  25 f/i=   29.978804 THz   188.362378 2PiTHz  999.985217 cm-1   123.982410 meV
  26 f/i=   35.952889 THz   225.898667 2PiTHz 1199.259268 cm-1   148.689252 meV
  27 f/i=  203.242065 THz  1277.007557 2PiTHz 6779.425348 cm-1   840.541919 meV

共 27 个振动模式,最后$f/i$指虚频。

前面我们提到过,虚频可以判断结构是否稳定。那这里,我们计算出的乙醇分子结构肯定不稳定喽?不一定。

因为频率计算和软件的数值积分有关(我也不清楚数值积分怎么进行的);

计算过程中我们的设置对频率计算影响很大,KPOINTS, ENCUT, EDIFF, POTIM 等都会影响计算的精度;综合这些因素,对于分子的振动频率来说(注意:声子谱不适用)一般低于 $100 cm^{−1}$ 的频率可以忽略。严格点可以降到 $50 cm^{−1}$,也就是说:如果你在计算中发现有个 $50 cm^{−1}$ 左右的虚频,完全可以不考虑。

零点能

零点能$ZPE=1/2 h \nu$

# 所有振动的能量之和 (所有的hv之和,单位meV)
grep 'f  =' OUTCAR | awk '{print $10}' | paste -sd+ | bc

# 零点能(eV)  将以下两行写脚本(meV转换eV除以1000,然后1/2,等于上式结果除以2000)
hv_sum=$(grep"f  =" OUTCAR | awk '{print  $10}'| paste -sd+ | bc)
echo "scale =6; $hv_sum/2000" | bc

零点能校正:

  1. 结构优化之后得到分子的能量(OSZICAR 中的$E_0$): $E_0$
  2. 频率计算后得到分子的零点能: $ZPE$
  3. 零点能校正之后分子的能量为:$E_{ZPE}=E_0+ZPE$

过渡态和反应热的零点能校正:

对一个反应:IS --> TS --> FS

  1. 优化反应物 IS 和产物 FS 的结构,获得能量:$E_0(\textrm{IS})$, $E_0(\textrm{FS})$;

  2. 对反应物和产物进行频率计算,获得各自的零点能:$\textrm{ZPE(IS)}, \textrm{ZPE(FS)}$。

  3. 搜索过渡态,获得结构和能量 $E_0(\textrm{TS})$;

  4. 过渡态频率分析,获得零点能 $\textrm{ZPE(TS)}$。

不考虑零点能的反应能垒 ($E_a$) 和反应热 ($\Delta E$):

$$

E_a = E_0(\textrm{TS}) – E_0(\textrm{IS}) \

\Delta E = E_0(\textrm{FS}) – E_0(\textrm{IS}) $$

考虑零点能校正:

$$ \begin{aligned}

E_{a}^{'} &= E_{\textrm{ZPE}}(\textrm{TS}) – E_{\textrm{ZPE}}(\textrm{IS}) \

 &= E_0(\textrm{TS}) + \textrm{ZPE(TS}) – E_0(\textrm{IS}) – \textrm{ZPE(IS}) \\
 &= E_a + \textrm{ZPE(TS)} – \textrm{ZPE(IS)}\\

\Delta E^{'} &= \Delta E + \textrm{ZPE(FS)} – \textrm{ZPE(IS)} \end{aligned} $$

零点能校正的情况:

频率计算时放开哪些原子看体系,看关注哪些部分。在过渡态中,IS、FS、TS 固定和放开的要一致。

影响频率计算的因素

Ex 27 乙醇分子的振动频率计算(五) | LVTHWopen in new window

EDIFFG,增强收敛标准对虚频并没有什么好的效果。

ENCUT,对零点能影响很小,增大截断能可以减小虚频,但并不是算频率就要增大截断能。

PREC,

POTIM

POINTS

备注

获取虚频

grep 'f/i'  OUTCAR | awk '{print $1 "\t " $2 "\t" $7 "\t" $8 "\t " $9 "\t" $10 "\t" $11}'

获取时间

grep Elapsed */OUTCAR | sort -n

待解决的问题

LVTHW 算出来 3 个虚频,我算的 13 个。哪里出了问题