使用Multiwfn計算雙描述符勢預測反應位點

使用Multiwfn計算雙描述符勢預測反應位點
Using Multiwfn to calculate dual descriptor potential to predict reactive sites

文/Sobereva@北京科音   2024-Jun-1

1 前言

雙描述符（dual descriptor）是概念密度泛函理論框架下定義的非常流行的預測親電和親核反應位點的方法，見《概念密度泛函綜述和重要文獻合集》（http://bbs.keinsci.com/thread-384-1-1.html）里面的相關資料和我寫過的各種相關博文的匯總，以及筆者的論文《親電取代反應中活性位點預測方法的比較》（http://www.whxb.pku.edu.cn/CN/abstract/abstract28694.shtml）和《Comparative study on the methods for predicting the reactive site of nucleophilic reaction》（https://doi.org/10.1007/s11426-015-5494-7）。最近，雙描述符勢（dual descriptor potential）在J. Math. Chem., 62, 1094 (2024)中被提出，文中指出它在預測反應位點方面比雙描述符具有一定優勢，原理上也更嚴格。從2024-May-13更新的Multiwfn開始，雙描述符勢的計算已被納入到了《使用Multiwfn超級方便地計算出概念密度泛函理論中定義的各種量》（http://www.shanxitv.org/484）所介紹的Multiwfn做概念密度泛函理論分析的功能（主功能22）當中。下文就簡單介紹一下雙描述符勢的概念，并給出很典型的計算例子。

2 原理

在說雙描述符勢之前先說福井勢（Fukui potential）。兩個分子之間反應傾向于在什么位點發生，取決于反應初期反應物之間的相互作用，在這方面靜電勢與福井勢之間存在明顯的互補關系。例如當一個親電物質正在與一個親核物質反應，親核物質對親電物質施加的外勢可以寫為下式，它對于親電物質與它的反應產生了引導作用。式中N-phile下標代表親核物質，ΔN是反應物之間轉移的電子數，V_ESP是靜電勢，f-是福井函數，積分部分對應福井勢。

反應初期的反應物的外勢分布特征決定了優先的反應位點。一般的化學反應總是伴隨著或多或少的電子轉移，而靜電勢則完全沒有考慮電子轉移對外勢的影響，因此靜電勢對反應位點的預測往往很失敗，這在前面提到的我的論文里的對比中都有體現。電子轉移對外勢的影響藉由福井勢體現了出來，而且當電子轉移越顯著（ΔN越大），福井勢的重要性就越顯著、對反應位點越起到決定性作用。

對應于福井函數有三種形式，f-、f+和f0，相應地福井勢也有三種形式，Vf-、Vf+和Vf0，如下所示

根據上面的介紹可知，用福井勢討論反應位點比用福井函數明顯更嚴格，因為前者才是與反應初期的外勢直接相聯系的，而福井函數只是間接地聯系。而之所以人們在討論反應位點時經常用福井函數而很少見到用福井勢，主要在于二者的分布特征通常相似、能得到的結論相同，而計算福井勢則需要計算復雜的積分，相當于計算靜電勢的成本，這遠比計算福井函數昂貴。

對f+和f-福井函數求差可得到雙描述符（Δf）。類似地，雙描述符勢可以對Δf做庫侖積分得到，也相當于Vf+和Vf-之間求差，如下所示

和Δf一樣，雙描述符勢的最正和最負的區域分別對應優先發生親核和親電反應的位點。雙描述符勢的計算代價顯著高于Δf，其好處是函數分布特征比Δf簡單得多，明顯更易于考察反應位點，而且理論上還更為嚴格，因為能與外勢直接掛鉤。

注意雖然基于Δf可以定義簡縮雙描述符（condensed dual descriptor）便于在不同原子間定量對比，但雙描述符勢并沒法搞成簡縮的形式。

3 實例

下面來看怎么用Multiwfn非常方便地計算雙描述符勢，并且通過此例能展現出雙描述符勢相對于Δf的實際好處。讀者務必先仔細閱讀《使用Multiwfn超級方便地計算出概念密度泛函理論中定義的各種量》（http://www.shanxitv.org/484）了解Multiwfn的主功能22的基本用法，這里假定讀者已經讀過此文。計算福井勢和雙描述符勢的功能從2024-May-13及以后版本的Multiwfn中才支持，切勿用更老的版本。Multiwfn可以在其官網http://www.shanxitv.org/multiwfn免費下載。

這個例子是用雙描述符勢預測吲哚（indole）的親電反應位點。下圖是網上找的，可見吲哚的親電反應優先發生在3號位的碳原子上，我們來看看雙描述符勢能否正確預測。下面的例子涉及到的所有文件都可以在http://www.shanxitv.org/attach/708/file.zip里得到。

首先用GaussView之類程序畫出來吲哚的結構，然后用Gaussian等量子化學程序做幾何優化。Gaussian在常用的B3LYP/6-31G*級別做優化的輸入輸出文件在本文的文件包里都提供了。之后啟動Multiwfn，然后輸入
indole.out  //Gaussian優化吲哚任務的輸出文件，Multiwfn會從中載入最后一幀結構
22  //概念密度泛函理論計算
1  //生成N、N-1和N+1態的.wfn文件（N代表原本狀態的電子數）
[回車]  //用默認的B3LYP/6-31G*級別計算
[回車]  //凈電荷和自旋多重度(0 1)、(-1 2)、(1 2)分別用于N、N-1和N+1態
此時當前目錄下已經產生了N.gjf、N+1.gjf和N-1.gjf，可以手動用Gaussian進行計算。如果你把Multiwfn的settings.ini里的gaupath設為了當前實際運行Gaussian的命令，這里也可以輸入y直接讓Multiwfn調用Gaussian去自動計算這些文件。算完后當前目錄下就有了N.wfn、N+1.wfn和N-1.wfn。之后輸入
9  //計算福井勢和雙描述符勢的格點數據
1  //由于當前體系小，而且福井勢和雙描述符勢的計算耗時較高，故這里使用低質量格點

在筆者的i9-13980HX機子上總共20多秒鐘算完。然后在當前的菜單里可以看到直接顯示各種形式的福井勢和雙描述符勢的選項，以及將它們導出成cub文件的選項。選擇4直接顯示雙描述符勢，并把等值面數值設為0.04（單位是a.u.），看到的圖如下所示，藍色代表與當前等值面數值符號相反，也即-0.04 a.u.的等值面。可見雙描述符勢非常清楚地顯示出C7（當于3位的碳）具有最負的福井勢，說明這是親電反應最優先發生的位點，和已知的實驗結論完全一致！

再來看看雙描述符的圖像是什么樣。輸入以下命令
0  //返回
3  //計算福井函數和雙描述符的格點數據
2  //中等質量格點
4  //顯示雙描述符的等值面圖
將等值面數值調整到易于分析的0.005，結果如下。可見雙描述符也體現出C7是優先發生親電反應的位點，因為藍色等值面在C7附近最大。然而雙描述符的等值面明顯顯得很零碎，有很多節面，遠不如雙描述符勢的等值面對反應位點展現得清晰。這在于雙描述符勢是對雙描述符求庫侖積分，在積分時把雙描述符復雜的函數分布特征給消除了。

在量子化學波函數分析與Multiwfn程序培訓班（http://www.keinsci.com/workshop/WFN_content.html）里專門講“反應位點的預測與反應活性分析”的部分里還給了很多其它的例子，這里再列舉其中一個，展現親核位點的：

4 總結

本文對最近提出的雙描述符勢的概念進行了簡要介紹，并演示了用Multiwfn計算雙描述符勢的具體例子。雙描述符勢比常用的雙描述符雖然計算起來更昂貴，但物理意義更嚴格、對反應位點展現得更清楚，因此當實際發現雙描述符預測反應位點不準時，或者其等值面過于復雜而難以預測反應位點時，或者就是想得到更可靠的反應位點預測結果時，都建議用雙描述符勢代替雙描述符。在前面提到的雙描述符勢的原文里還有更多的應用例子，感興趣者可以閱讀。

使用本文的方法計算在發表結果時務必記得按照Multiwfn程序啟動時的提示對Multiwfn進行正確的引用。