VASPKIT校正氣體分子自由能

目前至少80%的文獻中用VASP計算氣體分子並校正自由能的計算方法都是不準確的（或者說忽略了某些項的貢獻）！！因為VASP並不自帶自由能校正方法，而氣體分子自由能校正的過程又非常的繁瑣。用VASPKIT可以非常方便的讀取頻率計算的OUTCAR文件，自動做自由校正。

使用VASPKIT請引用：

Wang, V.; Xu, N.; Liu, J.-C.; Tang, G.; Geng, W. T. VASPKIT: A Pre- and Post-Processing Program for VASP code arXiv e-prints [Online], 2019. https://ui.adsabs.harvard.edu/abs/2019arXiv190808269W (accessed August 29, 2019).

理論計算催化中的反應自由能變化是非常重要的內容，但是主流第一性原理計算程序(VASP, QE, CP2K等)都沒有內置直接計算分子自由能的模塊 (注: 計算化學程序如Gaussian有很完善的自由能計算模塊)。導致大量計算催化文獻中的自由能校正方法不統一，不準確，甚至出現嚴重的定量錯誤 (比如：比如NRR反應往往需要用實驗值校正N2分子的自由能，某些文章計算出的總包反應能量和實驗值偏差巨大，竟然審稿人看不出來問題，最後錯誤的結果發表)。

無論是氣相化學反應，還是表面、材料的化學變化，都是在一定反應條件下進行的，但是直接用第一性原理程序計算出來的都是指在0K下忽略零點振動能(ZPE)的能量(俗稱電子能量)。自由能校正中首先要校正ZPE，然後計算配分函數q，然後校正體系內能(U)，焓(H)，熵(S)，ZPE。由此可以引入，反應溫度、壓力對能量的影響。隨後可以進一步考慮 pH、電極電勢等等外界因素對化學勢的影響，討論結果才和真實的實驗情況更為接近。

VASPKIT校正分子自由能方法

對氣體分子做熱力學量的校正需要計算振動、平動、轉動、電子激發對配分函數的貢獻，得到配分函數以後帶入到內能、焓、熵和配分函數的關係式中得到內能、焓、熵的校正值，再計算自由能的校正。這些過程涉及到較複雜的計算，好在VASPKIT已經把這所有的校正公式都編寫好了，用的時候只需要氣體分子的：

1，所有振動頻率，
2，氣體分壓，
3，溫度，
4，分子的基態多重度。

第一步：在一個大晶胞中優化分子，比如O2分子，

使用1 1 1的Gamma Center 的k點，
ISPIN = 2 (氧氣分子是三重態)，
ISYM = 0 (防止不合理的對稱性導致錯誤的電子佔據)，
ISMEAR = 0；SIGMA = 0.01 (對於分子體系可以設置很小的sigma確保精度)
EDIFF = 1E-7 (保證精度高一些，因為頻率計算和自由能校的精度對這兩個參數極為敏感)
EDIFFG = -0.001
晶胞邊長至少保證>10 Angstrom。

第二步：計算振動頻率。

修改以下INCAR參數：
IBRION = 5 或 6
NFREE = 2
POTIM = 0.015 (默認值)
EDIFF = 1E-7
計算完成以後grep cm OUTCAR，查看所有頻率：

得到3N個頻率，其中3N-6（直線分子3N-5）是屬于振動的。

  1 f  =   46.979964 THz   295.183821 2PiTHz 1567.082879 cm-1   194.293584 meV
   2 f  =    1.952595 THz    12.268518 2PiTHz   65.131565 cm-1     8.075288 meV
   3 f  =    1.120777 THz     7.042052 2PiTHz   37.385108 cm-1     4.635164 meV
   4 f/i=    0.006984 THz     0.043884 2PiTHz    0.232971 cm-1     0.028885 meV
   5 f/i=    1.046106 THz     6.572876 2PiTHz   34.894327 cm-1     4.326347 meV
   6 f/i=    1.294104 THz     8.131095 2PiTHz   43.166663 cm-1     5.351986 meV
因為線性分子的振動自由度為(3N-5)。發現有三個虛頻和兩個非常小的頻率，這是平動和轉動在振動自由度上的投影，直接忽略即可(VASPKIT自動判斷並忽略)。忽略最小的5個(線型分子)或6個(非線型分子)振動頻率，並不是直接忽略了平動和轉動的貢獻。而是通過平動和轉動的配分函數另外計算其對熱力學量的貢獻。其中平動熵是氣體分子熵的主要貢獻。
第三步：在頻率計算的文件夾裡，直接運行VASPKIT 502功能，依次輸入(1)溫度，(2)壓力，(3)自旋多重度。VASPKIT會自動判斷分子的對稱性和轉動慣量，用於計算轉動配分函數。
502
 +- Warm Tips -+
  -->> (01) Reading Structural Parameters from CONTCAR File...
  -->> (02) Analyzing Molecular Symmetry Information...
 Molecular Symmetry is: D(inf)h
 Linear molecule found!
 +---+
 Please input Temperature(K)!
298.15
 Please input Pressure(Atm)!
1
 Please input Spin multiplicity!--(Number of Unpaired electron + 1)
3
 -->>
  -->> (03) Extracting frequencies from OUTCAR...
  -->> (04) Reading OUTCAR File...
  -->> (05) Calculating Thermal Corrections...
 U(T) = ZPE + Delta_U(0->T)
 H(T) = U(T) + PV = ZPE + Delta_U(0->T) + PV
 G(T) = H(T) - TS = ZPE + Delta_U(0->T) + PV - TS

 Zero-point energy E_ZPE   :       2.240 kcal/mol    0.097144 eV
 Thermal correction to U(T):       3.724 kcal/mol    0.161475 eV
 Thermal correction to H(T):       4.316 kcal/mol    0.187167 eV
 Thermal correction to G(T):     -10.302 kcal/mol   -0.446723 eV
 Entropy S                 :     205.141 J/(mol*K)   0.002126 eV
可見最後的輸出格式和高斯程序的自由能校正模塊非常相似， 依次輸出，ZPE，內能，焓，自由能的校正值，最後還有熵。但是kcal/mol或eV。氧氣分子的自由能校正值是 -0.447 eV，而查閱JANAF-NIST熱力學數據表的數據大約是-0.45 eV，非常接近。
用G(T)加上電子能量(結構優化OSZICAR中最後一個E0，或者頻率計算中第一個E0的能量)值就是自由能校正之後的氧氣分子的能量。
用 -0.447 eV + 結構優化最後的電子能量(-9.855 eV) = -10.302 eV
注1：固體材料的自由能校正需要用phonopy計算聲子，然後做各種熱力學量的校正，這篇文章只針對氣體分子
注2：一般文獻中的表達，不把ZPE包含到內能校正中，所以G(T)可以寫為如下內容：
其中：H(T) = U(T) + PV
好多文獻中的校正方法沒有考慮U(T)，用的公式是G(T) = EDFT + ZPE + TS + PV, 是因為這些文獻的作者認為U(T)的貢獻很小可以忽略，實際上，對於氧氣分子這部分貢獻有大約有H(298.15 K) = 0.161475 - 0.097144 = 0.064331 eV。查Janaf-NIST熱力學表得到的 H(298.15 K) = 8.683 kJ/mol = 0.08999 eV, 由此可見對於精度要求稍高的計算，忽略H(T)是不合理的做法。
配分函數(Partition function)
統計熱力學中最關鍵的量是配分函數Q，
物理意義是體系的平均熱可及的量子態數，由它可以得到各種熱力學性質。計算中只能得到單個分子的性質，所以就需要把系綜的配分函數和分子的配分函數聯繫起來。含N 個分子的系綜的配分函數Q 與單個分子的配分函數q之間有如下關係：
藉由分子配分函數將有大量組成成分（通常為分子）系統中微觀物理狀態（例如：動能、勢能）與宏觀物理量統計規律 （例如：壓力、體積、溫度、熱力學函數、狀態方程等）連結起來的科學。如氣體分子系統中的壓力、體積、溫度。一個分子的能量可以認為是由分子的整體運動的能量(平動能)，以及分子內部運動的能量(轉動能，振動能，電子激發能)之和：
所有熱力學量都可以分解成這些組分的加和：
亥姆霍次自由能：
內能：
熵：
轉動、振動、電子躍遷對U的貢獻等價於對H的貢獻，對Cv的貢獻等價於對Cp的貢獻。(只差一個PV或P)，例如：