跳至主要內容

频率分析

Northword大约 8 分钟VASP反应路径

频率分析

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+16N+1 步。N 为体系中的振动的原子数，这是因为：

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

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

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

    官网原文如下，还要查阅 IBRION 和 NFREE 的相关内容。

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.`

例如，乙醇分子CH3CH2OHCH_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），坐标是分数坐标系。

四个频率的换算：

E=hc/λν=c/λν~=1/λT=1/νin which,E=energy (eV)λ=wavelength (m)λ~ = wavenumber (m−1)T=period (s)ν=frequency (s−1 or Hz)h=Planck’s constant = 4.135667516×10−15 eV⋅sc=speed of light = 299792458 m/s \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}

此外，1THz=1012Hz,1cm−1=100m−11 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/if/i指虚频。

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

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

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

零点能

零点能ZPE=1/2hν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 中的E0E_0）: E0E_0
  2. 频率计算后得到分子的零点能: ZPEZPE
  3. 零点能校正之后分子的能量为：EZPE=E0+ZPEE_{ZPE}=E_0+ZPE

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

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

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

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

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

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

不考虑零点能的反应能垒 (EaE_a) 和反应热 (ΔE\Delta E)：

Ea=E0(TS)–E0(IS)ΔE=E0(FS)–E0(IS) E_a = E_0(\textrm{TS}) – E_0(\textrm{IS}) \\ \Delta E = E_0(\textrm{FS}) – E_0(\textrm{IS})

考虑零点能校正：

Ea′=EZPE(TS)–EZPE(IS)=E0(TS)+ZPE(TS)–E0(IS)–ZPE(IS)=Ea+ZPE(TS)–ZPE(IS)ΔE′=ΔE+ZPE(FS)–ZPE(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 个。哪里出了问题