DFT-D4色散校正的簡介與使用

DFT-D4色散校正的簡介與使用

First release: 2019-Feb-13 Last update: 2019-Aug-10

之前筆者在《DFT-D色散校正的使用》（http://www.shanxitv.org/210）中詳細介紹了DFT-D3的函數形式以及使用方法，還有其它相關的文章：《亂談DFT-D》（http://www.shanxitv.org/83）、《談談“計算時是否需要加DFT-D3色散校正？”》（http://www.shanxitv.org/413）。大家如果對DFT-D3尚不了解的話請先閱讀以上幾篇。本文筆者簡要介紹一下DFT-D3的后繼者DFT-D4以及其實際使用方法。由于筆者精力有限，此文在原理方面就不寫得很詳細了。

1 關于DFT-D4方法

DFT-D3色散校正如今已極度流行，不光可以使那些計算弱相互作用差的泛函在計算弱相互作用精度方面有質的改進，對于計算很多中、大體系，特別是涉及過渡金屬反應的熱力學數據往往也有明顯改進。DFT-D3最大的一個軟肋就是沒法考慮體系的實際電子結構。同種元素的原子，在不同化學體系中的電子結構顯然是不同的，這從原子電荷上就可以直接體現出來。特別是過渡金屬，在不同配合物里往往以不同氧化態形式存在，雖然其原子電荷的差異沒有形式上的氧化態的差異那么夸張（參見《使用Multiwfn通過LOBA方法計算氧化態》http://www.shanxitv.org/362里的一些論述），但是差異還是相當顯著的。DFT-D方法計算色散校正能是基于原子的色散系數計算的，這從《DFT-D色散校正的使用》一文的公式里可以看得很清楚。實際上原子的色散系數是顯著依賴于原子的電子結構的，特別是原子電荷這個直接表征原子在分子中帶的凈電荷的量對色散系數的影響是明顯不可忽視的。但是由于DFT-D3沒有考慮這點，而只純粹依賴于體系的幾何結構，當體系中存在顯著帶電的原子的時候，比如過渡金屬配合物、離子化合物，DFT-D3在原理上就不是特別理想了（盡管從實際表現上來看，對于牽扯過渡金屬的體系，用DFT-D3比不用一般還是明顯有改進的）。

色散校正的方法很多，有的是像DFT-D3一樣不考慮實際電子結構，如OBS、DFT-ulg；也有的是能體現實際電子結構的，如XDM、dDSC、TS、LRD、vdW-DF/VV10、MBD@rsSCS等。那些能體現實際電子結構的做法普遍比形式超級簡單的DFT-D3更昂貴，編程實現往往更費勁，大多沒有解析梯度，被主流程序支持得也比較少，比如筆者撰文時最新的G16 B.01一個都不直接支持，而ORCA 4.1.1版只支持VV10（也被叫做DFT-NL）。

為了解決DFT-D3的上述不足，Grimme在JCP, 147, 034112 (2017)中提出了DFT-D4，關鍵性的改進就是讓原子色散系數與原子電荷掛鉤，從而能更好地反應實際情況。2017年這篇文章實際上是DFT-D4的最初版，發表之后，Grimme一直也沒公開提供能做DFT-D4計算的程序，令人覺得很莫名其妙，不禁猜測是不是DFT-D4還存在什么問題。后來在2018年7月，Grimme往ChemRxiv上上傳了介紹DFT-D4最終版文章的v1版，在2019年1月25日又上傳了這篇文章的v2版，參見https://doi.org/10.26434/chemrxiv.7430216.v2。我估計過不多久，這篇介紹DFT-D4最終版本的重量級文章就會正式在期刊上發表了。目前計算DFT-D4的程序已經公開了。
后記：DFT-D4最終版的正式文章已發表于J. Chem. Phys. 150, 154122 (2019)。

根據這篇預印版的DFT-D4最終形式的文章，可以看到最終形式的DFT-D4會有以下特征：
·和DFT-D3(BJ)一樣使用BJ阻尼函數。目前DFT-D4已經對所有主流泛函都擬合了阻尼參數
·帶有和DFT-D3一樣的Axilrod–Teller–Muto (ATM)形式的三體校正項。DFT-D4也可以結合MBD形式的三體校正項，但更為昂貴，而且實測的精度并沒有ATM那么好
·色散系數依賴于electronegativity equilibration (EEQ)方法計算的原子電荷。EEQ這個方法和已經很流行的EEM方法差不多，都是基于電負性均衡原理的思想，依賴于和元素有關的經驗參數，可以瞬間計算很大體系的原子電荷。EEM的一些介紹可以看Multiwfn程序手冊的3.9.15節
·有解析梯度（Hessian目前只有半數值的）。這是VV10之類其它依賴于電子結構的色散校正方法所不具備的
·最高支持到86號元素（Rn）
·和DFT-D3一樣基本上都是零耗時
·根據Grimme的測試，DFT-D4的精度能達到VV10這樣更復雜的色散校正方法的精度，甚至有時更好

對于那些不存在顯著帶電原子的情況，D4和D3精度半斤八兩，至少不輸于D3。而對于有顯著帶電原子的情況，D4比D3的改進還是挺明顯的，除了上述提到的Grimme的D4相關文章里有體現外，還可以看他的Acc. Chem. Res., 52, 258 (2019)，里面專門展示了DFT-D3、D4在金屬有機體系上的效果。

DFT-D4原理上也可以結合其它原子電荷。2017年那篇JCP里初代DFT-D4結合的是GFN-xTB方法（一種半經驗形式的DFT方法，類似DFTB）計算時順帶產生的Mulliken電荷，結合這種電荷時被稱為DFT-D4(TB)，但測試發現這種形式不僅算原子電荷的過程更麻煩（需要依賴于專做GFN-xTB的xtb程序），增加了耗時，結果還不如結合EEQ電荷的最終版DFT-D4好。

Grimme的文章也表示，DFT-D這種色散校正形式到了DFT-D4這一代就基本沒有油水進一步可榨了。由于最終形式的DFT-D4計算又快效果又好，應該在未來會很快普及、被大量量化程序直接支持。而且筆者也認為，有了DFT-D4，VV10等絕大部分以往提出的能夠展現實際電子結構的色散校正方法就可以徹底退出歷史舞臺了（但有個別這些方法由于仗著更高級的計算形式，精度可能做到比DFT-D4更高一些，所以還不是完全沒存在意義）。

當然，未來肯定有很多DFT-D4相關的測試文章接踵而至，也不排除DFT-D4方法的定義又出現一些改變，以上說法只是筆者撰文時的情況和看法。

2 Grimme的dftd4程序

2.1 簡介

Grimme寫的獨立的做DFT-D4的程序dftd4已經發布在了https://www.chemie.uni-bonn.de/pctc/mulliken-center/software/dftd4，可以直接下載源代碼包。下文說的是筆者2019年2月13日在此處下載到的dftd4.2.0.tar.xz的情況，以后版本可能有些地方會有變動。

這個DFT-D4程序用法、用處和DFT-D3很相似。由于此程序里面沒有包含xtb程序，因此做的DFT-D4只能是基于EEQ電荷的（如上所述，這也正是標準的DFT-D4的情況）。此程序的三體校正項支持前述的ATM和MBD形式，因此確切來說，此程序可以算的DFT-D4具體形式包括DFT-D4(EEQ)-ATM和DFT-D4(EEQ)-MBD。

dftd4得先編譯才能運行。為便于讀者使用，按照本節做法編譯好的Linux版和Windows版我直接提供在此：http://www.shanxitv.org/soft/dftd4.2.0_binary.rar。其中有.exe的是Win32版（libiomp5md.dll是其運行時需要的動態庫文件，應放到與之相同的目錄），沒后綴的是Linux版，都是可執行文件，直接就能運行。

2.2 最簡單粗暴的編譯方法

下面說一下如何以最簡單的方式編譯。

Windows 32bit版：筆者機子里裝了Visual Studio 2017以及Intel Parallel Studio XE 2019 Update 1里面的Intel Fortran compiler和Intel MKL庫。啟動Visual Studio 2017，新建一個項目，把dftd4壓縮包里的include目錄下的.f文件和source目錄下的.f90文件都放到項目的文件夾里，并且加入項目里。把編譯的配置從默認Debug改為Release。之后在項目屬性里，在Fortran標簽頁下把Libaries里的Runtime Library設成Multithreaded，Use Intel Math Kernel Library設Sequential。Language標簽頁里Process OpenMP Directives設Generate Parallel Code。之后如果你愿意的話還可以去設對CPU和指令集的優化選項，這里就不用了。然后對項目直接編譯，項目目錄下的Release目錄里就產生了dftd4.exe可執行文件，直接就能用了，而且傳給別人別人也能用。

Linux版：筆者的機子里已經裝了Intel Fortran compiler 19.0.1.144和MKL數學庫（是按照《VASP最簡單的安裝方法》http://www.shanxitv.org/455里面的做法裝的Intel Parallel Studio XE 2019 Update 1里面的）。將dftd4壓縮包解壓后，將include目錄下的.f文件拷入source目錄，再進入source目錄，直接運行以下命令即在當前目錄下產生了dftd4可執行文件
ifort *.f90 -qopenmp -qopenmp-link=static -mkl -static-intel -stand f08 -check bounds -o dftd4
其中那些帶static的選項是為了讓編譯出的可執行文件不依賴于動態庫，從而在別人的機子里也確保能運行。

2.3 標準但麻煩的編譯方法

dftd4基于meson+ninja的代碼生成系統組合，如果你想形式上最“優雅”地編譯（其實沒任何好處），需要系統里有meson和ninja。這倆玩意兒比較新，可能很多人不熟悉，這倆程序其實就相當于常見的cmake和make，只不過有一些額外的特點而已。筆者認為DFT-D4這么簡單的程序用mason+ninja純粹是莫名其妙，也給用戶添麻煩，一般系統里都不自帶它們，還得現裝。下面我們來通過meson+ninja編譯DFT-D4，機子里已經裝了上述Intel Fortran編譯器和MKL庫。筆者也嘗試了用此系統自帶的gfortran 4.8.5，發現編譯不過去。

先運行以下命令以安裝meson、lapack靜態庫（順帶會安裝blas靜態庫）、ninja。由于yum這么安裝之后ninja的可執行文件名是ninja-build，所以得做個符號鏈接成為ninja
yum install lapack-devel
yum install meson
yum install ninja-build
ln -s /usr/bin/ninja-build /usr/bin/ninja

將下載到的DFT-D4程序的壓縮包解壓，進入此目錄，直接將以下內容輸入命令行窗口，就會通過ifort來編譯
FC=ifort meson setup build && ninja -C build
然后在當前目錄下的build目錄下就出現了名為dftd4的可執行程序，可以直接用了。編譯過程中終端里會提示什么不識別-pipe之類選項的警告，不用管它，只要程序能用就行了。如果想運行得更快，而且你的機子里有MKL庫的話，可以把編譯時依賴的lapack給替換為MKL庫，把dftd4目錄下自帶的meson.build替換為我提供的這個即可：http://www.shanxitv.org/attach/464/meson.build。

2.4 基本用法

在含有dftd4的目錄下直接用./dftd4就可以顯示出此程序可以接的選項，一看就懂，就不多說了。比如下面的命令是對H2O.xyz這個文件里的體系基于BLYP泛函的參數計算標準形式的DFT-D4，即DFT-D4(EEQ)-ATM形式的色散校正能
./dftd4 H2O.xyz --func blyp

瞬間看到以下信息
           -------------------------------------------------
          |                Calculation Setup                |
           -------------------------------------------------
   coordinate file      : H2O.xyz
   number of atoms      :      3
   charge               :      0
   non-additivity corr. : ATM
   charge model         : EEQ
   functional           : b-lyp
   omp threads          :      1
   memory needed (est.) :    0.0 Mb

O      q(ref)        CN(ref)   cov. CN(ref)     α(AIM,ref)
     0.0000000      0.0000000      0.0000000      5.1967092
    -0.3547147      0.9924594      0.8041679      6.1223924
    -0.5988153      1.9886928      1.6112394      6.7491941
     0.0000000      0.9985672      0.9798440      5.1987934
H      q(ref)        CN(ref)   cov. CN(ref)     α(AIM,ref)
     0.0000000      0.0000000      0.0000000      5.0540161
     0.0000000      0.9117922      0.8942243      2.7207580

           -------------------------------------------------
          |               Damping Parameters                |
           -------------------------------------------------
   s6                   :     1.0000
   s8                   :     2.3408
   a1                   :     0.4449
   a2                   :     4.0933

           -------------------------------------------------
          |                     Results                     |
           -------------------------------------------------
Edisp /kcal,au:    -0.2980 -0.00047485

dispersion energy written to file .EDISP

即色散校正能為-0.00047485 Hartree。如果你用cat .EDISP顯示剛產生的.EDISP隱藏文件，會看到更高精度的數據：-0.00047485257046358403。

dftd4另外常用的選項是--chrg，用于設定體系凈電荷；--grad，用于計算色散校正能的梯度；--hess，用于計算色散校正能的半數值的Hessian。

DFT-D4程序支持的泛函在dftd4源代碼包里dfuncpar.f90里可以看到。在這個文件里搜get_d4eeqbjatm，凡是下面有名字的都可以用。值得一提的是，截止到2019-May-8，目前官網上的dftd4的這個參數庫文件里仍沒有常用的M06-2X泛函的參數，因此DFT-D4沒法結合M06-2X使用。

3 在量子化學程序中使用DFT-D4

后記：從ORCA 4.2開始，支持且只支持DFT-D4(EEQ)-ATM形式的DFT-D4，直接寫個D4關鍵詞即可使用。下面的關于ORCA的文字可以忽略了。

ORCA從4.1開始支持了DFT-D4，實際上是調用程序目錄里otool_dftd4可執行文件算的，這個文件其實就是專供ORCA的dftd4。只要在ORCA輸入文件里簡簡單單地寫上D4關鍵詞，計算時就會自動計算出D4校正能并加到總能量上。ORCA也支持DFT-D4下的優化任務。

對于筆者撰文時最新的ORCA 4.1.1，大家會發現ORCA輸出的DFT-D4校正能和直接用Grimme的dftd4程序給出的不同（選的都是相同泛函的情況下）。這是因為起碼對于ORCA 4.1.1來說，默認算的是DFT-D4(TB)形式的校正能，因此計算時會自動調用ORCA自帶的otool_xtb（專供ORCA的xtb程序）做GFN-xTB計算。由于DFT-D4(EEQ)-ATM形式更便宜且結果更好，我預感ORCA以后早晚會把D4關鍵詞的默認情況改成DFT-D4(EEQ)-ATM形式。

Molpro從2018.1開始支持了DFT-D4。日落西山的Turbomole從7.3也支持了DFT-D4。

不知Gaussian什么時候會正式支持DFT-D4，估計是早晚的事。對于尚不支持DFT-D4的Gaussian，如果想在Gaussian計算中使用的話，可以利用《將Gaussian與ORCA聯用搜索過渡態、產生IRC、做振動分析》（http://www.shanxitv.org/422）和《將Gaussian與Grimme的xtb程序聯用搜索過渡態、產生IRC、做振動分析》（http://www.shanxitv.org/421）里介紹的Gaussian的external關鍵詞自己寫個接口讓Gaussian可以在計算能量、受力、Hessian的時候調用DFT-D4，實現起來很簡單。