Multiwfn使用的高效的靜電勢算法的介紹文章已于PCCP期刊發表！
Multiwfn使用的高效的靜電勢算法的介紹文章已于PCCP期刊發表！
文/Sobereva@北京科音 2021-Aug-31
近日，深圳灣實驗室的張鋆和北京科音自然科學研究中心（http://www.keinsci.com）的盧天共同發表的名為Efficient Evaluation of Electrostatic Potential with Computerized Optimized Code的文章刊登于Phys. Chem. Chem. Phys., 23, 20323 (2021)，可訪問此鏈接查看：https://doi.org/10.1039/d1cp02805g。此文專門介紹了一種非常高效的靜電勢計算算法，并已實現在十分流行的波函數分析程序Multiwfn（http://www.shanxitv.org/multiwfn）中，下面簡要說明一下。
靜電勢（electrostatic potential, ESP）是量子化學中的一種極為重要的實空間函數，可以用來解釋和預測分子間靜電相互作用、計算擬合靜電勢電荷用于分子模擬、預測親核和親電反應位點、計算基于靜電勢定義的描述符預測體系凝聚相性質，等等。Multiwfn程序支持非常豐富的靜電勢相關的分析，已被大量計算化學文章所使用。筆者寫過諸多相關博文，連同一些重要文獻匯總在了《靜電勢與平均局部離子化能相關資料合集》（http://bbs.keinsci.com/thread-219-1-1.html）帖子里。
最早，Multiwfn的靜電勢計算代碼是大約2011年時候寫的，使用的是相對容易實現但速度較慢的代碼，當時做稍微大一些體系的靜電勢分析很費勁。后來，Multiwfn支持了調用Gaussian的cubegen程序做靜電勢計算部分，使得靜電勢相關的分析耗時有了極大降低，介紹見《Multiwfn現已可以調用cubegen使靜電勢分析耗時有飛躍式的下降！》（http://www.shanxitv.org/435）。再后來，曾經開發過LIBRETA電子積分庫的張鋆和筆者共同合作，最終在Multiwfn 3.7版程序中理想地實現了一種非常高效的靜電勢計算算法，筆者在《Multiwfn的計算靜電勢的內部代碼速度得到了巨幅提升！》（http://www.shanxitv.org/563）專門做了介紹。這個算法如今被正式發表在前述的Phys. Chem. Chem. Phys.期刊上，并且也是當前版本算靜電勢的時候默認使用的。
Phys. Chem. Chem. Phys.這篇文章里關于算靜電勢的電子積分層面的數學細節在此不做介紹，請感興趣的讀者自行閱讀，這里只把文章中的靜電勢計算速度測試部分說一下。
下面對兩個體系在不同的有代表性的基組下按照《RESP擬合靜電勢電荷的原理以及在Multiwfn中的計算》（http://www.shanxitv.org/441）文中的做法計算RESP原子電荷，這個過程耗時幾乎完全來自于計算大量的擬合點上的靜電勢。其中較小的多巴胺（Dopamine）體系的耗時是在普通8核筆記本上測的，較大的瑞德西韋（Remdesivir）體系的耗時是在雙路E5-2696v3共36核的服務器上測的。表中的orca_vpot是ORCA量子化學程序中自帶的算靜電勢的獨立程序，code1是最知名的商業量子化學程序自帶的速度非常快的能夠計算靜電勢的工具，Multiwfn (old)和(new)分別是最早的Multiwfn內部的靜電勢計算代碼和目前版本默認的代碼，N_GTO是Primitive Gaussian函數數目。
由上表可見，Multiwfn的新靜電勢代碼比最早的代碼快了幾十倍，也顯著勝于ORCA的計算靜電勢的工具，在多數情況下比需要付費的代碼速度還快。
為了證明Multiwfn的新的靜電勢代碼可以算得動很大體系，筆者從MOF-5晶體當中截取了一個重復單元，邊緣用氫恰當飽和，此時化學組成是H120C144O104Zn32，共400個原子。然后通過xtb程序在GFN2-xTB理論下計算得到了molden格式的波函數文件，用于Multiwfn計算表面靜電勢，然后按照《使用Multiwfn+VMD快速地繪制靜電勢著色的分子范德華表面圖和分子間穿透圖》（http://www.shanxitv.org/443）中的ESPpt的做法作圖，結果如下所示。此體系的Gaussian函數數目多達4840個，被計算靜電勢的表面頂點數多達259262個，但即便如此，在36核服務器下花了半小時也算完了。這說明Multiwfn的靜電勢分析、結合VMD對靜電勢繪圖已經可以用于相當大的體系。
Phys. Chem. Chem. Phys.文中對Multiwfn的靜電勢計算的并行效率也做了測試。測試的是對ATP分子在B3LYP/6-311+G(2d,p)下產生的波函數計算RESP電荷。如下圖可見，并行效率很理想，隨著調用的核數增多，相對于單線程運行的速度提升得也很理想。因此Multiwfn的靜電勢分析可以充分利用核數較多的CPU的運算能力。
雖然Multiwfn調用Gaussian的cubegen計算靜電勢的速度也不錯，但相對于直接用Multiwfn當前內部的靜電勢計算代碼有幾個不足：
(1)必須花錢購買程序
(2)只直接支持fch格式。而Multiwfn的靜電勢代碼是特別普適的，用任意Multiwfn支持的格式作為輸入文件都可以，目前包括fch、molden、mwfn、wfn、wfx、gms（不過，Multiwfn倒是也有將mwfn/molden/gms轉換為fch的功能，例如早先寫過的《巨大體系的范德華表面靜電勢圖的快速繪制方法》http://www.shanxitv.org/481里把xtb產生的molden文件轉成了fch再讓Multiwfn調用cubegen算靜電勢）。
(3)Multiwfn有些涉及靜電勢的功能不支持調用cubegen，比如對靜電勢用主功能2做拓撲分析。特別值得一提的是在《巨幅降低Multiwfn結合VMD繪制分子表面靜電勢圖耗時的一個關鍵技巧》（http://www.shanxitv.org/602）中筆者提出的巨幅降低繪制分子表面靜電勢圖的方法不能結合cubegen使用，用此文的方法依靠Multiwfn自己的靜電勢計算代碼的耗時遠低于結合cubegen的時候。
目前只有一個情況讓Multiwfn調用cubegen算靜電勢的耗時比用Multiwfn自己的代碼明顯更低，就是利用主功能5計算均勻分布的靜電勢格點數據的時候（《繪制靜電勢全局極小點+等值面圖展現孤對電子位置的方法》http://www.shanxitv.org/493等文章介紹的分析里需要用到）。大抵cubegen專門為這種情況做了特殊考慮，不是一個一個點單獨計算而是批量計算，期間共享了很多中間數據，通過避免重復計算一些中間數據而大幅節約了時間。所以如果算中、大體系靜電勢格點數據的時候有條件的話依然建議調用cubegen算靜電勢，沒有cubegen的話就盡量用CPU核數較多的機子算，并且格點間距設得略大一些讓被計算的點數不那么多。
Phys. Chem. Chem. Phys.這篇文章提出的靜電勢算法以及Multiwfn中相應的靜電勢計算代碼已經很高效了，但還有進一步改進的余地，可能速度在未來還會有不小的提升，包括
(1)更充分地利用CPU的SIMD指令集
(2)充分利用積分屏蔽來降低巨大體系的耗時
(3)針對收縮型基函數優化代碼（當前是純粹基于primitive函數算的。這倒也帶來一個好處是即便是用wfn、wfx這樣不含基函數信息的輸入文件也可以算靜電勢）
最后提醒一下，如果大家用Multiwfn做靜電勢相關分析時使用的是程序內部代碼計算的靜電勢，即默認的情況，請在發表的文章中對前述Phys. Chem. Chem. Phys., 23, 20323 (2021)文章進行提及和引用，例如可以在Computational Details段落里寫上：The electrostatic potential involved in the analyses was evaluated by Multiwfn based on the highly effective algorithm proposed in Ref.[PCCP文章]。