http://www.shanxitv.org/category/第一性原理 http://www.shanxitv.org/739 2025-03-31T23:18:00+08:00 使用CP2K過程中常用的可視化工具 文/Sobereva@北京科音   2025-Mar-31 0 前言 使用計算化學程序過程中普遍離不開可視化問題，對如今已經非常流行、十分強大、巨快、開源免費的第一性原理程序CP2K自然也不例外。最近有人在思想家公社QQ群里問“cp2k可視化一般搭配哪個程序呢？”，針對這個問題，筆者簡要羅列一下我推薦的在CP2K使用過程中的可視化程序，主要面向初學者。本文只是粗略概述，這些程序在筆者講授的北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/KFP）中都會反復利用和詳細演示具體操作。除本文說的外，還有大量其它可視化程序，比如Avogadro、gabedit、Chemcraft等，我認為對于CP2K用戶來說，用好本文提到的這些就已經完全足夠了。文中提到的Multiwfn可以在官網http://www.shanxitv.org/multiwfn下載，相關信息見《Multiwfn FAQ》（http://www.shanxitv.org/452）。 1 建模階段的可視化 能夠構建用于做第一性原理計算的周期性體系的程序很多，其中我十分推薦GaussView，這是絕大多數內行做量子化學計算的人都會用的可視化程序，其圖形界面設計得很好，所見即所得，大多數情況用這一個程序就可以很好完成CP2K計算前的建模。GaussView可以直接讀取記錄周期性體系原子坐標和晶胞信息的最常見的cif文件。在GaussView里構建好結構后，可以保存成Gaussian程序的輸入文件（gjf），里面有原子坐標信息和晶胞的平移矢量信息（Tv開頭的行），然后用Multiwfn載入這樣的gjf文件后，就可以按照《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）文中說的創建CP2K輸入文件了，整個過程相當方便。 許多人誤以為GaussView只是結合Gaussian使用的，或者至多是用于量子化學程序算孤立體系建模的，這是大錯特錯！由于Gaussian本來就有計算周期性的功能（但由于功能性和效率的原因，極少有人用Gaussian算周期性體系），所以Gaussian的御用圖形界面GaussView也支持周期性體系，可以直接對周期性體系添加/刪除/替換原子、調整鍵長/鍵角/二面角、平移/旋轉片段、把其它體系粘貼到當前體系里、測量幾何參數，等等，對周期性體系和對分子體系操作一樣非常靈活和順手。并且GaussView專門有個PBC editor界面，在里面可以定義晶胞參數、重新定義晶胞、修改周期性、晶軸變換、擴胞、切晶面、判斷空間群、居中體系、刪除晶胞外的原子、把晶胞外的原子卷入晶胞，等等，常用的操作幾乎都提供了。GaussView對于周期性體系建模相當好用，然而其這方面的價值卻被很多人低估，甚至有人明明不了解還胡亂貶損。 點擊下圖紅色箭頭處的按鈕就可以進入前述的PBC editor界面。紅、綠、藍邊框分別是晶胞的a、b、c軸。   還值得一提的是，GaussView的菜單欄的Tools中的Atom selection界面很有用。用滑框選擇工具選擇一批原子成為黃色后，Atom selection界面里就可以看到選中的原子序號，格式和Multiwfn程序里要求的完全一致。例如在Multiwfn里創建幾何優化任務的輸入文件時，若你選擇要凍結某些原子，就可以把序號從這里復制出來并直接粘貼到Multiwfn的窗口中，這樣實現諸如slab邊緣原子的凍結設置巨方便。   Multiwfn還可以載入CP2K的輸入文件、Quantum ESPRESSO的輸入文件、VASP的POSCAR等格式，主功能100的子功能2里選擇保存為gjf或者cif文件后，也可以用GaussView載入并觀看和編輯。順帶一提，在建模方面Multiwfn也提供了很多實用功能，很多還是GaussView沒有的，見《Multiwfn中非常實用的幾何操作和坐標變換功能介紹》（http://www.shanxitv.org/610）。 2 幾何優化結果的可視化 如果要看CP2K做幾何優化得到的最終結構，最簡單的方法是用Multiwfn載入此任務產生的restart文件，之后進入主功能0直接就看到了，例如下面的聚噻吩體系。界面上方的菜單以及圖形界面右側有一大堆選項可以調整效果，諸如視圖的旋轉和縮放、鍵的粗細、成鍵判斷閾值、原子序號是否顯示/字號/顏色、原子球大小、高亮特定原子。菜單欄的Tools列表里還有一些輔助工具，比如測量幾何參數、導出所有內坐標、輸出笛卡爾或分數坐標，等等。   很多時候我們還關心優化的過程，尤其是當體系結構復雜時，不觀看優化軌跡的話，往往都不好判斷優化過程中哪里發生了何種變化，心里沒數。觀看優化軌跡最好的程序是VMD，在http://www.ks.uiuc.edu/Research/vmd/可免費下載，強烈建議用VMD 1.9.3（撰此文時1.9.4還沒發布正式版，其測試版的bug巨多，強烈不建議用）。CP2K的幾何優化任務默認會輸出.xyz文件，這是多幀的軌跡文件，里面記錄了優化過程的每一幀的坐標，不熟悉xyz格式的話看《談談記錄化學體系結構的xyz文件》（http://www.shanxitv.org/477）。把xyz文件載入VMD就可以觀看優化軌跡了。即便任務沒跑完，也可以載入此文件看看最新一步的結構。 xyz文件里不記錄晶胞信息，因此在VMD里沒法顯示出晶胞邊框，給觀看周期性體系造成不便。對于不變胞的優化，可以用Multiwfn載入restart文件，此時屏幕上不僅顯示晶胞信息，還會顯示在VMD里定義晶胞的命令，如下圖所示。  把這命令復制到VMD的命令行窗口運行，之后再輸入pbc box命令，晶胞在VMD里就顯示出來了，如下所示。   如果做的是變胞優化任務，晶胞在優化過程中不是固定不變的，此時建議讓CP2K輸出pdb格式的優化軌跡，里面記錄了每一幀的晶胞信息。此文件載入VMD并顯示盒子后，會看到隨著軌跡的播放，盒子也實時變化。 上述說明不僅適用于優化極小點，對于dimer方法優化過渡態也是一樣的。 VESTA（http://jp-minerals.org/vesta/en/）也是一個很好且免費的主要面向周期性體系的可視化工具，默認情況下的顯示效果往往比VMD甚至還更好些。若想用VESTA看CP2K優化的結果，可以把restart文件載入Multiwfn，主功能100的子功能2里選擇導出cif或VASP的輸入文件（POSCAR），然后就能載入VESTA了。 3 反應路徑的可視化 CP2K做NEB類任務需要對給定的初始結構之間插點。CP2K雖然能夠自動插點（點=image），但沒法在計算前進行預覽。利用《sobNEB：產生CP2K的NEB的插點的方便的工具》（http://www.shanxitv.org/660）里介紹的筆者的sobNEB程序插點的話，可以直接產生traj.xyz軌跡文件，放到VMD里通過播放動畫或者多幀疊加顯示，就可以直觀判斷初始的NEB的image是否合理了。 NEB類任務計算過程中是好多副本一起算的，每個副本都會在計算過程中輸出它所負責的那個image優化過程的xyz軌跡文件。要想觀看最終收斂的NEB軌跡，需要自己把這些副本輸出的軌跡的最后一幀合并在一起作為新的xyz軌跡，放到VMD里就可以觀看最終的反應路徑動畫了。手動做比較麻煩，建議用腳本實現。北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/KFP）里專門給了個shell腳本自動做這件事情，運行后就會在當前目錄下產生[項目名]_traj.xyz文件，里面記錄了NEB最新一步的軌跡（即便任務還沒跑完，也可以看最新的NEB軌跡是什么樣），同時還會產生[項目名]_ene.txt，里面記錄了最新一步的各個點的能量。下面是培訓里的一個實例，對載入的[項目名]_traj.xyz文件做多幀疊加顯示，直觀顯示了Na的遷移軌跡   4 振動分析的可視化 CP2K做完振動分析后往往需要觀看振動矢量了解振動的特征。推薦做法是用《使CP2K計算的振動模式可以被GaussView觀看的程序：MfakeG》（http://www.shanxitv.org/656）里介紹的筆者開發的MfakeG工具，把記錄振動模式信息的.mol后綴的molden文件轉化成仿Gaussian輸出文件，之后載入GaussView就可以用Results - Vibrations選項觀看振動模式了，例如下圖的草酸晶體的例子   還有一種做法是使用免費的可視化程序Jmol（https://jmol.sourceforge.net），也可以載入CP2K產生的molden文件觀看振動模式。不過我還是覺得GaussView用起來舒服。 如果要基于CP2K的振動分析觀看振動光譜，推薦用Multiwfn。Multiwfn具有非常強大的繪制各種光譜的功能，見《使用Multiwfn繪制紅外、拉曼、UV-Vis、ECD、VCD和ROA光譜圖》（http://www.shanxitv.org/224）和《使用Multiwfn繪制NMR譜》（http://www.shanxitv.org/565）。Multiwfn可以載入CP2K做振動分析產生的輸出文件按博文所述繪制紅外光譜，如果要求計算了拉曼還可以繪制拉曼光譜。此外，Multiwfn還可以載入常規TDDFT和XAS-TDDFT任務的輸出文件分別繪制UV-Vis光譜和X光吸收光譜，還可以考慮旋軌耦合效應，例如下圖是CP2K+Multiwfn繪制的NaAlO2晶體的Al的K-edge XAS。Multiwfn還可以載入CP2K的NMR任務的輸出文件繪制NMR譜。這些在北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/KFP）里有非常系統的講解和豐富的例子。   5 動力學軌跡可視化 CP2K的第一性原理動力學（叫從頭算動力學也行）是CP2K的杰出強項之一。和GROMACS、AMBER、NAMD等經典力場動力學程序的情況一樣，VMD也是觀看其動力學軌跡的不二的選擇。例如培訓里有個模擬質子轟擊石墨烯層的例子，VMD載入軌跡后并多幀疊加顯示、根據幀號著色，直觀展現了模擬過程中質子的運動軌跡：   還可以自己寫VMD腳本從CP2K產生的記錄原子速度的xyz文件中把速度信息讀入VMD，并根據速度著色，展現出質子的動能是如何傳遞到石墨烯板上并擴散開來的。下圖黃色是打入的質子，越藍的石墨烯原子的速度越大。   VMD不僅可以觀看動力學過程中結構的變化，還可以繪制等值面圖觀看電子結構的變化，例如下圖這個例子，直觀展現了水合電子是怎么在模擬過程中出現的。   通過寫VMD tcl腳本，還可以實現很復雜的分析，如下面幻燈片里講的高溫下正癸烷的裂解產物分析。PS：做動力學的，不會寫分析腳本的話，稍微做深一點就會寸步難行。   北京科音分子動力學與GROMACS培訓班（http://www.keinsci.com/KGMX）專門用約100頁幻燈片特別完整、詳細講VMD的使用，并且還額外用近100頁幻燈片深入講VMD腳本的編寫，如果想系統學習VMD的話這是極佳的途徑。 筆者偶爾看到有人用OVITO程序可視化CP2K產生的軌跡，我沒用過那個程序，也完全不理解為什么有人不用VMD而用OVITO，明明VMD已經可以完美地滿足一切需求。筆者在答疑時看到有一些CP2K用戶還被OVITO坑了：OVITO根據邊緣原子位置自動確定了盒子邊框，有人以為那就是實際的晶胞邊界，導致對結果產生了嚴重誤解。 6 電子結構的可視化和分析 Multiwfn可以基于CP2K的輸出文件繪制效果非常理想的能帶結構圖，如CrO2體系：   《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》（http://www.shanxitv.org/651）中介紹了怎么用CP2K產生molden文件。Multiwfn將之載入后，就可以在《使用Multiwfn繪制態密度(DOS)圖考察電子結構》（http://www.shanxitv.org/482）講的基礎上舉一反三繪制效果很好的DOS、PDOS、OPDOS、LDOS圖，下圖是WO3體系：   Multiwfn還能基于CP2K的molden文件觀看軌道，做法可參考《使用Multiwfn觀看分子軌道》（http://www.shanxitv.org/269）。還支持對特定k點看軌道：   Multiwfn還能對周期性體系做超級豐富的波函數分析，其中很多都是以圖形方式展現的，比如可以基于《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）和《使用Multiwfn考察周期性體系的芳香性》（http://www.shanxitv.org/722）講的做法計算LOL-pi格點數據，之后可以在Multiwfn里直接觀看等值面；也可以導出成cube文件后，載入VMD或VESTA程序繪制，如下所示，極為生動直觀地展現了一個COF體系的pi電子共軛路徑   再比如下圖是Multiwfn載入CP2K對硅表面計算產生的molden文件，做軌道定域化后直接顯示出軌道等值面圖，直觀展現了體系不同位置的電子結構特征。相關信息見《Multiwfn的軌道定域化功能的使用以及與NBO、AdNDP分析的對比》（http://www.shanxitv.org/380）。   為避免此章太長，Multiwfn可以針對周期性體系做的巨量的可視化分析這里就不一一提及，很多分析我都寫過博文，感興趣的讀者請閱讀： 使用Multiwfn結合CP2K的波函數模擬周期性體系的隧道掃描顯微鏡（STM）圖像 http://www.shanxitv.org/671（http://bbs.keinsci.com/thread-37740-1-1.html） 使用Multiwfn對周期性體系做鍵級分析和NAdO分析考察成鍵特征 http://www.shanxitv.org/719（http://bbs.keinsci.com/thread-47176-1-1.html） 使用Multiwfn結合CP2K做周期性體系的atom-in-molecules (AIM)拓撲分析 http://www.shanxitv.org/717（http://bbs.keinsci.com/thread-46927-1-1.html） 使用Multiwfn結合CP2K對周期性體系做電荷分解分析（CDA） http://www.shanxitv.org/716（http://bbs.keinsci.com/thread-46878-1-1.html） 使用Multiwfn做IGMH分析非常清晰直觀地展現化學體系中的相互作用 http://www.shanxitv.org/621（http://bbs.keinsci.com/thread-28147-1-1.html） 使用IRI方法圖形化考察化學體系中的化學鍵和弱相互作用 http://www.shanxitv.org/598（http://bbs.keinsci.com/thread-23457-1-1.html） 使用Multiwfn結合CP2K通過NCI和IGM方法圖形化考察固體和表面的弱相互作用 http://www.shanxitv.org/588（http://bbs.keinsci.com/thread-21742-1-1.html） 使用Multiwfn繪制分子和固體表面的距離投影圖 http://www.shanxitv.org/589（http://bbs.keinsci.com/thread-21754-1-1.html） 使用Multiwfn圖形化展現原子對色散能的貢獻以及色散密度 http://www.shanxitv.org/705（http://bbs.keinsci.com/thread-44723-1-1.html） 使用Multiwfn做Hirshfeld surface分析直觀展現分子晶體和復合物中的相互作用 http://www.shanxitv.org/701（http://bbs.keinsci.com/thread-43178-1-1.html） 使用CP2K結合Multiwfn對周期性體系模擬UV-Vis光譜和考察電子激發態 http://www.shanxitv.org/634（http://bbs.keinsci.com/thread-28006-1-1.html） 使用Multiwfn計算分子和晶體中孔洞的直徑 http://www.shanxitv.org/643（http://bbs.keinsci.com/thread-30696-1-1.html） 使用Multiwfn計算晶體結構中自由區域的體積、圖形化展現自由區域 http://www.shanxitv.org/617（http://bbs.keinsci.com/thread-25241-1-1.html） 使用CP2K結合Multiwfn繪制密度差圖、平面平均密度差曲線和電荷位移曲線 http://www.shanxitv.org/638（http://bbs.keinsci.com/thread-28225-1-1.html） 此外，用VMD載入CP2K產生的Hartree勢（靜電勢的負值），還可以實現對晶體表面靜電勢的可視化，一個COF體系的例子如下所示。靜電勢在研究靜電主導的相互作用方面有重要意義，參見《靜電勢與平均局部離子化能相關資料合集》（http://bbs.keinsci.com/thread-219-1-1.html）里面的資料。   和靜電勢同等重要的是筆者提出的范德華勢，適用于考察范德華作用主導的相互作用，可以由Multiwfn計算，見《談談范德華勢以及在Multiwfn中的計算、分析和繪制》（http://www.shanxitv.org/551）。 http://www.shanxitv.org/729 2024-11-28T23:00:00+08:00 駁網上流傳的對CP2K缺點的不實描述 Refuting the false description of shortcomings about CP2K circulating on the Internet 文/Sobereva@北京科音  2024-Nov-28 0 前言 CP2K是極其流行的第一性原理程序，在《2024年計算化學公社舉辦的計算化學程序和DFT泛函的流行程度投票結果》（http://www.shanxitv.org/706）的流行度得票統計中甚至已超過了VASP。筆者在北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/KFP）里也積極推廣CP2K、全面介紹其相關理論背景知識和正確的使用方式。筆者之前在網上答疑時，就曾偶爾聽到一些關于CP2K的不實謠言，諸如什么算導體慢、算相互作用能不準之類。今天在群里看到有人發了一頁講CP2K缺點的幻燈片，里面列了四條，根本就沒有一條是對的，令我深感這種幻燈片實在太以訛傳訛了！估計已有很多對CP2K一無所知的人看了這種幻燈片對CP2K產生了嚴重錯誤印象，導致其放棄了其實極其適合他們的CP2K而轉向使用某些又貴又慢的程序。因此筆者忍不住寫本文，對幻燈片里的不實說法進行駁斥和斧正。實際上，如果讀者參加過北京科音CP2K第一性原理計算培訓班，正確、深入、系統學習過CP2K計算相關的理論知識，并結合培訓里的例子用CP2K做過各類計算，就自然而然就知道那幻燈片里對CP2K缺點的說法有多么離譜。我也希望讀者看到有人再給CP2K亂扣屎盆子時轉發本文鏈接。 1 對導體計算慢 那幻燈片里說CP2K“對導體計算較慢，CP2K之所以算的快是因為OT算法，對有帶隙的體系可以迅速SCF迭代收斂。但是此法原理上不適用于導體。所以導體只能用對角化方法”。 以上說法大錯特錯。首先說CP2K為什么快。對于純泛函的DFT計算，CP2K對中、大體系遠比VASP、Quantum ESPRESSO、Abinit等純粹基于平面波的程序快，這在于CP2K是原子中心基函數（具體來說是Gaussian型基函數）展開軌道波函數，同時用傳統平面波程序的方式計算經典庫侖作用。用原子中心基函數使得CP2K的得以利用矩陣的稀疏性巨幅節約算較大體系的時間，而且CP2K的開發者把程序效率視為重中之重，把CP2K的效率優化到了極致。由于相關公式較多，這里直接貼北京科音CP2K第一性原理計算培訓班講CP2K的GPW算法的兩頁ppt進行說明 此外，對于雜化泛函的計算，使用原子中心基函數的CP2K由于可以充分利用積分屏蔽技術大幅節約要計算的雙電子積分的量，使得CP2K在主流的雙路服務器上做雜化泛函計算算到上千原子都不是難事。 上述原因是CP2K做DFT特別快的最關鍵原因，而不是在于CP2K有OT。OT只不過是在這CP2K本來就特別高效的基礎上進一步令大體系每一輪的SCF耗時更低，因為OT避免了算大體系耗時很高的對角化巨大KS矩陣的步驟（體系越大時越重要，對小體系OT和對角化的差異不大）。另外，用OT時SCF往往比對角化更容易收斂。一些相關信息看《CP2K中遇到SCF難收斂時的解決方法》（http://www.shanxitv.org/665）。即便CP2K只用傳統對角化算法做SCF而不用OT，算中、大體系的速度依然吊打VASP等純平面波程序。 算金屬體系時需要開smearing，此時需要算空軌道，由于用OT做SCF過程中無法求解空軌道，因此OT一般來說不能用于金屬體系。金屬體系的KS矩陣不是很稀疏，因此前述的CP2K基于原子中心基函數帶來的優勢沒有非金屬的情況那么明顯。但即便如此，CP2K算不小的以金屬為主的體系也完全不比平面波程序慢！（這里只從gamma點的能量和受力計算來說，至于對稱性的利用與k點的考慮、SCF收斂性等其它方面另談，那是其它層面的事） 上述的北京科音CP2K第一性原理計算培訓班里的豐富的例子中就有不少是金屬為主的體系，比如“Pd(100)晶面對苯分子的吸附”、“Cu(001)表面的銀原子遷移”、“Au(111)表面H2分子解離”、“Ir(001)表面活化甲烷”等，學員親自用CP2K跑過一遍這些例子就會感受到CP2K對這些金屬為主的體系的計算效率也非常高，何來“對導體計算較慢”？ 另外，前述網傳的ppt中說OT“對有帶隙的體系可以迅速SCF迭代收斂”也是錯誤說法。不少情況下OT達到收斂所需的SCF圈數反倒顯著多于對角化。OT也并不保證肯定能收斂，而且也經常遇到OT需要一百輪以上SCF才收斂的收斂緩慢的情況。OT遠沒有那么神乎其神。 2 磁性體系處理比較麻煩 那幻燈片里說CP2K“對于磁性體系處理比較麻煩，因為CP2K計算需要提前指定自旋多重度，不像VASP一樣可以自動尋找合適的磁矩”。 以上說法嚴重不實。為了保證收斂到正確的磁性狀態，CP2K需要用MAGNETIZATION指定初猜波函數中的原子自旋磁矩或者用&BS指定初始的原子電子組態，VASP也需要用MAGMOM來設定初猜波函數中的原子自旋磁矩，顯然在這點上CP2K根本沒有什么額外的更麻煩之處。 CP2K的計算根本就不是必須“提前指定自旋多重度”。北京科音CP2K第一性原理計算培訓班的下面這兩頁幻燈片里我把情況說得非常清楚： 3 k點功能不完善、小體系計算不準確 那幻燈片里說CP2K“K點功能不完善，只有用Gamma點計算比較好用，在添加的K點參數以后一些很好的加速算法就不能用了比如ASPC波函數外推，老版本的CP2K里只有Gamma點計算的功能。這意味著對大體系可以，但是小體系計算不準確” 以上說法極具誤導性。說CP2K的k點（不叫K點）的功能不完善我沒意見，確實尚有不少特性不支持k點，截止到2024版包括OT、NMR、SCCS溶劑模型、DFT+U等等，但若說成“只有用Gamma點計算比較好用”就完全是在開玩笑。CP2K早就能在純泛函計算的時候很好地考慮k點，而且效率也很好，在能量計算、（變胞）優化、CI-NEB產生反應路徑和搜索過渡態、能帶計算等等方面都表現得很好，北京科音CP2K第一性原理計算培訓班里大量例子都是開著k點做的，根本什么問題都沒有。而且如http://bbs.keinsci.com/thread-43683-1-1.html所說，CP2K從2024開始還支持了雜化泛函結合k點計算的RI-HFk算法。CP2K做DFT計算對k點考慮的一個不足（至少截止到撰文時最新的2024版）是不支持利用晶格的對稱性節約k點來降低耗時（而Quantum ESPRESSO、VASP支持），但這也不是什么大問題，因為本來CP2K就已經很快了。 CP2K支持波函數的外推，也就是在結構優化、動力學過程中用之前幾步的波函數外推出這一步的較好的初猜波函數，ASPC是其中一種。用k點時任何外推方法都沒法用，因此上面那段話提到的ASPC不支持k點不假。然而這在實際中僅僅是一個小問題，上面那段話把這問題的重要性完全夸大了。考慮k點時做幾何優化完全可以用USE_PREV_P直接使用上一步SCF收斂的密度矩陣當初猜，從而讓SCF收斂更快、降低幾何優化每一步的耗時。而分子動力學用的盒子尺寸通常不小，否則會有虛假的周期性效應導致動力學行為/現象不合理，因此此時一般都用不著考慮k點。 那幻燈片里說CP2K“小體系計算不準確”純屬無稽之談。CP2K考慮k點計算小晶胞體系不僅效率上足夠好，也不存在任何bug，何來不準確？北京科音CP2K第一性原理計算培訓班講“能量的計算及相關問題”的部分還專門給了個計算例子，令學員清楚地看到CP2K做n*n*n的超胞只考慮gamma點算的能量正好是用n*n*n k點算原胞的n*n*n倍，充分體現出CP2K考慮k點計算小體系的正確性沒有絲毫問題。 4 高斯基組帶來巨大的BSSE 那幻燈片里說CP2K用的“高斯基組帶來巨大的BSSE，基組不完備誤差，不適合計算結合能。VASP使用平面波基組就不存在這個問題”。 以上說法真是**！ BSSE問題我在《談談BSSE校正與Gaussian對它的處理》（http://www.shanxitv.org/46）有簡要介紹，北京科音中級量子化學培訓班（http://www.keinsci.com/KBQC）里“弱相互作用的計算與相關問題”一節有深入介紹。BSSE是計算弱相互作用能時候需要注意的問題，而不說具體基組談BSSE的大小根本毫無意義！！！CP2K能用的基組多了去了，最常用的是MOLOPT系列，按照此順序尺寸依次增大、BSSE依次減小：SZV-MOLOPT-GTH、DZVP-MOLOPT-SR-GTH、DZVP-MOLOPT-GTH、TZVP-MOLOPT-GTH、TZV2P-MOLOPT-GTH、TZV2PX-MOLOPT-GTH。在計算弱相互作用能方面，常用的DZVP-MOLOPT-SR-GTH有嚴重的BSSE問題，而到了TZV2P-MOLOPT-GTH檔次BSSE問題就已經相當小了。下面是MOLOPT基組原文JCP, 127, 114105 (2007)中的不同基組的對比測試，m-開頭的對應MOLOPT系列，可見TZV2P-MOLOPT-GTH算小分子相互作用能BSSE帶來的誤差僅僅有區區不到0.2 kcal/mol，遠遠小于相互作用能自身的數量級（水分子二聚體間相互作用能約5 kcal/mol），難道這能叫“巨大的BSSE”？更何況CP2K還直接支持做counterpoise校正，可以令本來就很小的BSSE問題更是完全忽略不計。雖然算大體系之間弱相互作用會比小體系間的BSSE更大，但只要用到TZV2P-MOLOPT-GTH結合counterpoise檔次依然誤差微乎其微。 有些人對CP2K使用的基組可謂近乎一無所知，就隨便網上看幾個零碎的（往往還是不適合自己情況或者有嚴重誤導性的）例子，拿比如TZVP-GTH、DZVP-MOLOPT-SR-GTH這種BSSE挺大的基組算弱相互作用能，也不知道counterpoise是何物，發現算出來結果差（一般會嚴重高估）就直接賴CP2K算不準，這些人太令人無語了。我不得不說，參加北京科音CP2K第一性原理計算培訓班把這些必備知識好好系統學學再做計算實在太重要了。以不恰當的方式做計算，不管什么程序都算不準，諸如VASP計算時用很小的平面波截斷能、Gaussian計算時用6-31G，結果能不爛么？ 注意前面說的是弱相互作用能，對于涉及到成鍵的相互作用，如化學吸附的能量變化，相同的基組體現出的BSSE問題會更小，并且也不建議用counterpoise校正，見《計算化學鍵鍵能時以counterpoise方式考慮BSSE不僅是多余的甚至是有害的》（http://www.shanxitv.org/381），直接用足夠大的比如TZV2P-MOLOPT-GTH檔次做單點計算就夠了（優化和振動分析一般用DZVP-MOLOPT-SR-GTH即可）。 前述網上的幻燈片里的說法不僅嚴重誤導了CP2K的業余和潛在用戶，還令所有使用Gaussian基函數的程序中槍！諸如Gaussian、ORCA、NWChem、GAMESS-US等等等等，難道這些被廣為使用的量子化學程序算相互作用能都有巨大誤差不成！？ 那幻燈片最后還來了句“VASP使用平面波基組就不存在這個問題”，顯得好像VASP、Quantum ESPRESSO平面波程序相對于CP2K就明顯很適合算相互作用似的。不僅前面說了，這種說法有嚴重誤導性，而且還忽略了純平面波程序在這方面的一個普遍劣勢，即必須用三維周期性。關于這點，挺值得一提的是在Mol. Phys., 117, 1298 (2019)文中，作者專門比較了Quantum ESPRESSO以三維周期性的方式和Q-Chem以孤立體系的方式計算的22種小分子相互作用能，結論是基于平面波計算的時候必須盒子取得很大，即便開偶極校正，也建議至少留15埃真空層，以確保令誤差小到不至于嚴重影響結果。CP2K相較于純平面波程序的一個優點是還可以以非周期性（像量子化學程序一樣）、一維周期性、二維周期性的方式做計算，因此在非周期性方向可以不考慮周期性，相應方向就可以用小得多的真空區（具體看Poission solver的選擇，在培訓班里會很詳細講），避免需要設很大的真空區浪費耗時。 值得一提的是，CP2K里有個專門做純平面波計算的SIRIUS模塊，可以得到和Quantum ESPRESSO極為相近的結果。如果你就是有特殊原因非要做平面波計算也可以用這個模塊。 5 CP2K真正有哪些缺點？ 在我來看CP2K的缺點是以下這些 上面這張幻燈片里關于CP2K的k點的局限性前文已經說過一部分了，對一般應用并不構成問題。如果做的某些計算真是必須考慮k點但CP2K此時不支持的話，實際上擴胞到只需要考慮gamma點也可以解決，由于CP2K算大體系的速度非常快此時一般也都算得動。 上面這張幻燈片里說的第三方輔助工具相對有限，注意這只是目前跟VASP、Quantum ESPRESSO等歷史更長的程序相比，實際上CP2K已有大量輔助工具可以用。比如創建輸入文件可以非常便利地用Multiwfn實現，見《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）；Multiwfn基于CP2K的輸出文件可以做非常豐富的波函數分析以及繪制DOS和能帶圖；計算各種熱力學量CP2K可以結合Shermo，見《使用Shermo結合量子化學程序方便地計算分子的各種熱力學數據》（http://www.shanxitv.org/552）；觀看振動模式可以用MfakeG+GaussView，見http://www.shanxitv.org/soft/MfakeG；計算聲子譜可以CP2K結合Phonopy；做晶體結構預測可以CPK結合USPEX；做非絕熱動力學可以CP2K結合Newton-X；ASE有CP2K的接口，等等。由于CP2K的巨大、無可取代的優點，無疑以后輔助工具還會越來越多！ 上面這張幻燈片里說的算某些磁性的d族金屬表面體系SCF難收斂，典型例子就是諸如Ni(111)表面、Fe(001)表面這種，而且改各種SCF收斂相關設置也不太好解決。碰到這種情況可以用CP2K的SIRIUS模塊基于平面波計算（CP2K用戶做DFT默認用的是Quickstep模塊），或者改用免費的Quantum ESPRESSO程序。注意絕對絕對絕對、千萬千萬千萬不要以為CP2K算磁性體系就有SCF必然難收斂的問題！北京科音CP2K第一性原理計算培訓班里講了很多例子，諸如“Fe2O3的鐵磁性和反鐵磁性狀態的單點計算”、“Fe (bcc)原胞的單點計算”、“反鐵磁性的UO2的DFT+U的計算”、“順磁性物質CuCl2晶胞的優化”等等，都是順利收斂的。 上面這張幻燈片里說的沒有解析Hessian，因此只能基于解析受力做有限差分得到Hessian，是相對于Quantum ESPRESSO、VASP的一個缺點。但在實際中也不是大問題。有對稱性的小晶胞的振動分析結合Phonopy程序可以利用對稱性巨幅節約要做的有限差分的次數，再加上CP2K算其中要涉及的超胞的受力計算本來就很快，所以振動分析也花不了什么時間。而對于沒對稱性的原子較多的大晶胞，多數情況也可以通過恰當的凍結巨幅節約要做的有限差分的次數（如計算表面吸附，基底的不與被吸附分子接觸的區域都可以凍住）。 除了上述以外CP2K還有些次要缺點，例如SCCS溶劑模型下容易遇到SCF難收斂、考慮k點的變胞優化容易中途卡住（但完全可以解決，培訓里講了，計算化學公社論壇首頁google框搜也能找到我的回答）等等。 總的來說，CP2K的缺點相對于其優點是相當次要的，因此我十分推薦把極為高效、功能十分強大、特別流行還開源免費的CP2K作為第一性原理計算的默認選擇、當做主力程序，結合GaussView建模和Multiwfn產生輸入文件和做各種后處理分析（北京科音CP2K第一性原理計算培訓班里有非常全面系統的講解，零散的相關博文見http://www.shanxitv.org/category/CP2K/里的文章）用起來相當絲滑，有額外需求時再利用其它程序作為補充。 6 正確分辨網上關于CP2K的說法 最后我想再強調一下，大家切勿隨便聽信網上的可信度不明的關于CP2K的評價。比如今天看到群里有人說什么“CP2K對磁性系統太難收斂了”、“我的體系沒有磁性，也很慢”。這種說法非常不負責任。但凡是嚴謹的內行，首先就不會在沒有具體前提的情況下說這種話，而如果不是內行，又有什么能力客觀評價CP2K？像是上來就說CP2K磁性體系難收斂的，我要反問：具體是什么磁性體系（是否容易收斂都是case by case的事）？你確認初始原子磁矩合理設置了么？其它程序算這個體系就能容易收斂么、有參照物么？沒有這些最基本前提，就說“CP2K對磁性系統太難收斂了”毫無意義。至于那個說“我的體系沒有磁性，也很慢”的，我要問：慢具體是多慢、什么體系花了多少耗時？你的計算資源如何、什么CPU多少核并行？確保以合理的方式做并行計算了么（別告訴我說你用的是往往很慢的ssmp）？算的什么任務？用的什么計算級別和CUTOFF等直接影響耗時的設定？你監控SCF或幾何優化步數了么、是否發生震蕩？我這里說的這些因素都是嚴重關乎收斂性和耗時的，但凡有人沒有把這些細節準確交代出來就批CP2K，大概率其水平堪憂，因此切勿隨便聽信他們的說法，否則極大概率會誤導你。我在網上答疑時見到過太多外行人自己根本不懂怎么正確使用程序，遇到使用程序不順心的情況就直接賴程序有問題，他們應該自我檢討。 http://www.shanxitv.org/719 2024-07-12T14:52:00+08:00 使用Multiwfn對周期性體系做鍵級分析和NAdO分析考察成鍵特征 Using Multiwfn to perform bond order analysis and NAdO analysis to study bonding characters for periodic systems 文/Sobereva@北京科音  2024-Jul-12 0 前言 鍵級是考察化學鍵特征的一類重要方法，非常強大的波函數分析程序Multiwfn可以方便地計算很多類型的鍵級，在《Multiwfn支持的分析化學鍵的方法一覽》（http://www.shanxitv.org/471）的鍵級部分有具體的說明，Multiwfn手冊4.9節有一些實際計算例子，在“量子化學波函數分析與Multiwfn程序培訓班”（http://www.keinsci.com/workshop/WFN_content.html）里做了特別全面詳細的講解和演示。另外，Multiwfn還支持NAdO方法，可以將原子間的模糊鍵級以軌道形式圖形化展現，對于了解成鍵本質極其有用，介紹和實例看《使用鍵級密度(BOD)和自然適應性軌道(NAdO)圖形化研究化學鍵》（http://www.shanxitv.org/535）。 以上博文都是舉的孤立體系的例子，而如今Multiwfn已經將上述方法擴展到了對周期性體系的分析上。將Multiwfn與CP2K做周期性DFT計算產生的波函數相結合，可以很容易地用上述方法考察周期性體系的成鍵情況，在本文將通過一個層狀的共價有機框架化合物（COF）體系，以及一個Pd表面吸附苯的體系，對具體做法進行演示，希望讀者舉一反三將本文介紹的方法運用到實際研究中。使用本文的方法計算的結果若用于發表，記得需要按照Multiwfn啟動時的提示恰當引用Multiwfn的原文。 讀者務必使用2024-Jun-24及以后更新的Multiwfn版本，否則情況與本文所述不符。Multiwfn可以在主頁http://www.shanxitv.org/multiwfn免費下載。不了解Multiwfn者建議閱讀《Multiwfn入門tips》（http://www.shanxitv.org/167）和《Multiwfn FAQ》（http://www.shanxitv.org/452）。本文要用到非常流行、高效且免費的第一性原理程序CP2K，筆者假定讀者已具備了相關常識，推薦不熟悉者通過“北京科音CP2K第一性原理計算培訓班”（http://www.keinsci.com/workshop/KFP_content.html）系統性學習。本文要使用Multiwfn創建CP2K輸入文件，這在《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）里有簡要介紹。本文用的CP2K是2024.1版。 下面的例子涉及到的所有文件都可以在http://www.shanxitv.org/attach/719/file.rar里得到。 1 單層COF體系 這部分示例的COF體系的cif文件是本文文件包里的COF_16371N2.cif，用GaussView顯示的結構如下，可見晶胞里有兩層。 1.1 準備輸入文件 此體系兩層之間的相互作用主要是pi-pi堆積，對電子結構影響微乎其微，因此只計算一層就夠了，可以節約時間，尤其是能顯著節約NAdO分析過程中計算原子重疊矩陣（AOM）的時間。于是在GaussView里刪除一層，然后保存為COF.gjf。由于當前體系不是導體且晶胞邊長約15埃已經不小了，所以可以不用擴胞就直接做gamma點計算。注意CP2K導出molden文件只支持gamma點的情況，這在《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》（http://www.shanxitv.org/651）里已明確說了。 啟動Multiwfn，然后輸入 COF.gjf cp2k  //進入產生CP2K輸入文件的界面 [回車]  //輸入文件名用默認的COF.inp -2  //要求產生molden文件 -7  //設置周期性 XY  //僅在平行于COF的方向考慮周期性，而垂直于COF的方向用非周期性 -9  //其它設置 16  //要求將Kohn-Sham矩陣導出到.csr文件里。此文件在Multiwfn計算NAdO軌道能量時會用到 0  //返回 0  //產生輸入文件 現在當前目錄下就有了COF.inp，對應PBE/DZVP-MOLOPT-SR-GTH級別的單點計算。用CP2K運行之，當前目錄下會產生記錄波函數信息的COF-MOS-1_0.molden和記錄KS矩陣的COF-KS_SPIN_1-1_0.csr。按照《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》所述，在molden文件開頭插入以下內容定義晶胞信息和當前贗勢基組描述的各元素的價電子數  [Cell]  14.94670000     0.00000000     0.00000000  -7.47335000    12.94422190     0.00000000   0.00000000     0.00000000     8.00000000  [Nval]  C 4  N 5  H 1 之后就可以基于這個molden格式的波函數文件做各種周期性體系的分析了。 1.2 計算Mayer鍵級 首先計算非常常用的Mayer鍵級，它衡量原子間等效的共享電子的對數。啟動Multiwfn，載入COF-MOS-1_0.molden，然后輸入 9  //鍵級分析 1  //Mayer鍵級 Multiwfn首先計算重疊矩陣，然后立刻輸出了Mayer鍵級，只有大于閾值的（由settings.ini里的bndordthres參數控制）才直接顯示了出來。要想看所有原子間的Mayer鍵級可以再選擇y導出鍵級矩陣然后查看相應的矩陣元。為了便于對照，下圖顯示了當前體系各個原子的序號，本文主要關注其中綠色虛線圈住的那些原子。 Multiwfn輸出的一些有代表性的鍵的Mayer鍵級值如下 11(C )   19(C )    1.12454138 11(C )   22(N )    1.48268711 19(C )   48(C )    1.27625346 20(C )   50(C )    1.39369508 18(C )   48(C )    1.24722702 可見這些鍵的Mayer鍵級都不同程度地明顯大于1.0，因此可以推測它們都不僅僅形成了sigma鍵，還有一定pi成份。在后文會計算pi鍵級和使用NAdO分析更進一步確認這一點。上面列出的C-C鍵鍵級可以體現出體系中哪種C-C鍵的強度更強。可見連接萘單元和C3N3單元的C11-C19的鍵是相對最弱的。 1.3 計算模糊鍵級（fuzzy bond order） Mayer鍵級怕彌散函數，而模糊鍵級則沒有這個問題，用于任意基組都可以。Mayer鍵級用于MOLOPT系列基組都是沒問題的，不過這里作為例子，也計算一下模糊鍵級。由于計算模糊鍵級需要先構造AOM，對較大的體系很耗時，所以對于Mayer鍵級適用的情況建議優先用Mayer鍵級。注意對孤立體系Multiwfn計算模糊鍵級默認是基于Becke劃分的，而對于周期性體系計算模糊鍵級默認是基于Hirshfeld劃分的，因此默認設置下周期性體系的計算結果和孤立體系的沒有嚴格的可比性。 在鍵級分析主功能的菜單里選擇7，然后Multiwfn就會開始計算AOM，之后立刻給出模糊鍵級的計算結果，前述的那些鍵的模糊鍵級如下 11(C )   19(C )    1.05202861 11(C )   22(N )    1.50866214 19(C )   48(C )    1.22313834 20(C )   50(C )    1.38494183 18(C )   48(C )    1.17984661 雖然模糊鍵級與Mayer鍵級在定量上有一些差別，但不同的鍵的鍵級的大小順序是完全一致的，沒有結論上的差異。 1.4 計算pi電子貢獻的Mayer鍵級 我在《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）中專門介紹了pi鍵級的計算方法，也即pi占據軌道對Mayer鍵級的貢獻，沒讀過的話務必先閱讀。這一節對單層COF這個周期性體系也做這個計算，看看各個鍵的pi作用程度的差異。由于當前體系是精確平行于XY平面的，所以讓Multiwfn自動指認pi軌道前不需要先做軌道定域化。 啟動Multiwfn并載入COF-MOS-1_0.molden后，依次輸入 100  //其它功能（Part 1） 22  //檢測pi軌道 0  //當前的軌道都是離域形式（如正則分子軌道） 2  //設置其它軌道占據數為0 0  //返回 9  //鍵級分析 1  //Mayer鍵級 前面提到的那些鍵的鍵級計算結果如下，這對應的就是pi鍵級。可見C3N3六元環中的C-N鍵的pi共享電子作用很強，比這里列出的C-C鍵的還要更顯著。 11(C )   19(C )    0.14583409 11(C )   22(N )    0.41775145 19(C )   48(C )    0.28202057 20(C )   50(C )    0.38270736 18(C )   48(C )    0.27414145 使用此做法發表文章時，除了引用Multiwfn原文外還建議同時引用介紹Multiwfn做pi電子結構分析的文章Theor. Chem. Acc., 139, 25 (2020)。 1.5 NAdO分析 利用前述的《使用AdNDP方法以及ELF/LOL、多中心鍵級研究多中心鍵》中介紹的NAdO分析，可以將sigma和pi作用分別以軌道形式呈現出來，非常有價值。例如在我的研究18氮環的ChemPhysChem (2024) https://doi.org/10.1002/cphc.202400377文中就用了NAdO分析非常清楚、嚴格地展示了較短的N-N鍵的sigma+pi作用特征。本節我們用NAdO分析展示一下C11-C19的sigma和pi軌道相互作用情況。 啟動Multiwfn并載入COF-MOS-1_0.molden后，依次輸入 15  //模糊空間分析 3  //計算AOM。對周期性體系默認是基于Hirshfeld原子空間計算的 [回車]  //用默認的0.2 Bohr格點間距。此設置對當前體系耗時不高，但對于明顯更大體系為了節約耗時可以適當用更大的格點間距如0.35 Bohr，但會多多少少犧牲精度 計算完成后當前目錄下出現了記錄AOM的AOM.txt。在Multiwfn里接著輸入 0  //返回 200  //其它功能（Part 2） 20  //BOD和NAdO分析 -1  //要求計算NAdO的能量 2  //從外部文件讀取Fock矩陣 COF-KS_SPIN_1-1_0.csr  //輸入此文件的實際路徑 1  //基于AOM計算原子間的相互作用 [回車]  //讀取當前目錄下的AOM.txt 11,19  //分析11和19號原子間的相互作用 馬上NAdO分析就做完了，屏幕上顯示以下內容。這說明所有NAdO的本征值加和為1.052，和上一節得到的模糊鍵級一致。并且當前只有前兩個NAdO的本征值明顯大于0，它們是最值得考察的，其它的基本可以忽略  Eigenvalues of NAdOs: (sum=   1.05204 )    0.78955   0.19252   0.06784   0.02133   0.00807   0.00651   0.00127    0.00096   0.00031   0.00011   0.00010   0.00005   0.00004   0.00000 ...略 所有的NAdO軌道現在都已經導出到了當前目錄下的NAdOs.mwfn中。現在輸入y讓Multiwfn載入此文件，之后退回到Multiwfn主菜單，進入主功能0查看NAdO軌道。不熟悉Multiwfn看軌道功能的話參考《使用Multiwfn觀看分子軌道》（http://www.shanxitv.org/269）。在圖形界面里分別選擇1號和2號軌道，可分別看到下面的圖像，同時在文本窗口顯示的NAdO本征值和軌道能量也標上了。這里等值面數值用的是0.08。 由上可見，C11-C19同時具有sigma和pi作用特征，前者遠比后者顯著得多。而且由于C19處在pi共軛區域，所以NAdO方法產生的主要描述C11-C19的pi作用的軌道并不完全定域在C11和C19上，也同時一定程度離域到了周圍原子上。通過NAdO軌道能量還可以看到，C11和C19之間共享的sigma電子的能量明顯比pi電子要低，這十分符合一般化學常識。 2 Pd(100)表面吸附苯 本節通過Pd(100)晶面吸附苯的體系，演示一下苯與Pd(100)基底這兩個片段間Mayer鍵級的計算，以及對于像這樣很大的體系怎么盡可能節約整個NAdO分析的時間。此體系的CP2K做結構優化的輸入和輸出文件，以及優化任務跑完后產生的restart文件都在本文文件包里提供了。優化完的結構如下，可見吸附作用很強，以至于Pd-C的作用都令苯產生彎曲了。 由于此例子的molden文件、FOM.txt和NAdOs.mwfn文件太大，本文的文件包里就沒提供了。 2.1 準備輸入文件 用Multiwfn載入本文文件包里的Pd+ben_opt-1.restart讀取優化完的坐標和晶胞信息，然后輸入以下命令創建一個做單點且產生molden文件的任務 cp2k  //產生CP2K輸入文件 Pd+ben_SP.inp  //產生的輸入文件名 -7  //設置周期性 XY -2  //要求產生molden文件 6  //由于Pd基底是導體，開啟smearing 0  //產生輸入文件 為了幫助SCF收斂，把Pd+ben_SP.inp里的ALPHA設為0.15，NBROYDEN設為12。用CP2K計算Pd+ben_SP.inp，得到Pd+ben_SP-MOS-1_0.molden。然后在里面開頭部分加入以下內容。  [Cell]  23.33880000     0.00000000     0.00000000  0.00000000    23.33880000     0.00000000  0.00000000     0.00000000    24.79750184  [Nval]  Pd 18  C 4  H 1 2.2 片段間Mayer鍵級 這一節將苯和Pd(100)基底分別定義為一個片段計算它們之間的Mayer鍵級。啟動Multiwfn，載入Pd+ben_SP-MOS-1_0.molden，然后輸入 9  //鍵級分析 -1  //定義片段 289-300  //苯的部分 1-288  //Pd(100)基底部分 1  //計算鍵級 Multiwfn開始產生重疊矩陣，之后給出了原子間鍵級，以下是其中一部分，可見苯中的C和有的Pd之間的Mayer鍵級較大，即共價作用很顯著。 ...略  # 1412:       246(Pd)  270(Pd)    0.21887832  # 1413:       246(Pd)  271(Pd)    0.21771061  # 1414:       246(Pd)  289(C )    0.08497491  # 1415:       246(Pd)  293(C )    0.09391700  # 1416:       246(Pd)  294(C )    0.55236636  # 1417:       247(Pd)  271(Pd)    0.23611928  # 1418:       247(Pd)  272(Pd)    0.23841126 ...略 老有人問怎么判斷固體表面對分子是物理吸附還是化學吸附，判斷方法很多，如結合能的數量級、原子間距離和原子半徑的關系、電子密度差、IGMH（http://www.shanxitv.org/621）等等，而用Mayer或模糊鍵級是最能說明問題的。以上Mayer鍵級體現的C和Pd之間的明顯的共價作用意味著二者之間形成了明顯的化學鍵作用，顯然Pd(100)對苯是化學吸附。如果鍵級只有比如零點零幾的程度，那一般可以說是物理吸附，但也不排除形成典型離子鍵的可能，這通過原子電荷可以判斷，計算方式可參考《使用Multiwfn對周期性體系計算Hirshfeld(-I)、CM5和MBIS原子電荷》（http://www.shanxitv.org/712）。 最后Multiwfn還給出了片段間鍵級： The bond order between fragment 1 and 2:    4.309784 即苯分子和Pd(100)之間等效共享電子對數多達4.3，顯然太共價了，如果這都不叫化學吸附... 2.3 NAdO分析 這一節我們將Pd(100)與苯之間顯著的共價作用通過NAdO方法以軌道形式直觀展現。本節的做法和1.5節示例的常規的NAdO分析流程有兩個明顯的不同： (1)對當前這樣原子數又多、軌道數又多的大體系計算所有原子的AOM是相當耗時間的事情，而且導出、載入AOM.txt的耗時巨高而且巨占硬盤（正比于占據軌道數的平方乘以原子數）。為了極大程度節約時間，可以使用Multiwfn提供的一個特殊的策略，也就是在主功能15里用選項33計算感興趣的兩個片段間的片段重疊矩陣（FOM）并導出為FOM.txt，然后對兩個片段之間做NAdO分析時直接從中讀取要用的FOM，這樣不僅省得導出和載入巨大的AOM.txt，而且產生FOM的過程中不需要計算片段里沒有涉及的原子的AOM，當片段里的原子數占整個體系的原子數比例較少的時候這可以大幅節約時間。對于當前的例子，苯可以作為片段1，與苯相互作用最密切的那部分Pd原子適合作為片段2，這比起把所有Pd都定義為片段2能節約巨量耗時。但用戶很難判斷哪些Pd原子與苯作用密切，此時可以利用Multiwfn特意提供的一種定義片段2的方法，即片段1以外的原子中，與片段1任何一個原子之間Mayer鍵級大于指定閾值的原子都被定義為片段2。利用這個做法可以使得與苯的軌道相互作用最顯著的Pd原子被定義為片段2，其數目遠少于總的Pd數目。 (2)在計算FOM和做NAdO分析之前，需要先用主功能6的子功能38令電子根據軌道能量由低到高重新排列，使得軌道占據數都為整數。這是因為當前CP2K計算用了smearing，導致費米能級附近的軌道處于分數占據狀態，這不是單行列式波函數應具有的特征，而Multiwfn的NAdO分析目前只支持單行列式波函數。重新排列后，當前的波函數就成了標準的單行列式波函數形式，從而可以做NAdO分析。另一方面，這么重排占據數后在計算AOM和FOM時可以只對占據軌道之間計算，也只有這部分是NAdO分析所實際需要的，只計算占據軌道之間的AOM/FOM的耗時、對內存的占用遠遠低于考慮所有軌道的情況， 啟動Multiwfn，載入Pd+ben_SP-MOS-1_0.molden，然后輸入 6  //修改波函數 38  //按照軌道能量由低到高重排占據數 -1  //返回 15  //模糊空間分析 33  //計算FOM 3  //定義兩個片段并計算FOM，不屬于片段1的原子若與片段1的任意原子的Mayer鍵級大于指定值則被定義為片段2 289-300  //苯的部分 0.001  //Mayer鍵級的閾值 之后看到以下信息，告訴了你哪些原子被作為了片段2，并且可見片段1和2之間的Mayer鍵級為4.114，和上一節看到的苯與整個Pd(100)表面的片段間Mayer鍵級4.310相差很小，這說明當前自動定義的片段2是合理的。如果你把片段2的原子在Multiwfn主功能0里用菜單欄的Other settings - Set atom highlighting功能高亮顯示，會看到它們都是離苯不遠的原子。  Atoms in fragment 2:  76,77,80-83,86,87,97,101-103,106-108,112,195-197,200-202,205-207,220,221,224-22  6,230,231,241,245-247,250-252,256,265,266,269-272,274-277,280,281  Mayer bond order between fragments 1 and 2:   4.114 之后輸入0.35，Multiwfn就開始用0.35 Bohr間距的立方格點計算AOM，然后構造出FOM并導出到當前目錄下的FOM.txt。在雙路7R32機子上96核并行計算，整個過程花了26分鐘（耗時與格點間距的三次方呈反比，耗時太高的話可以用比如0.5 Bohr格點間距，誤差也還可以接受，此時只需要不到10分鐘）。然后輸入 0  //返回 200  //其它功能（Part 2） 20  //BOD/NAdO分析 4  //直接從當前目錄下的FOM.txt中讀取NAdO分析要用的片段1和2的FOM 之后開始了NAdO的計算，結果如下  Eigenvalues of NAdOs: (sum=   6.95512 )    0.79476   0.77393   0.71745   0.63686   0.62698   0.35965   0.22325    0.16243   0.14986   0.14655   0.14545   0.13818   0.12618   0.12058    0.11388   0.11078   0.10617   0.09335   0.09007   0.08979   0.08423    0.08224   0.07942   0.07434   0.07391   0.07171   0.06475   0.05751    0.05507   0.05478   0.05246   0.04964   0.04669   0.04030   0.03771 ...略 當前顯示的所有NAdO本征的加和6.955對應于苯和Pd(100)之間在Hirshfeld原子空間劃分下算的模糊鍵級，和前面看到的片段間Mayer鍵級4.31有一定差異，這很正常，畢竟對原子空間的定義截然不同。由以上數據可見有不少NAdO軌道對模糊鍵級的貢獻都顯著，尤其是前六個。現在輸入y載入新產生的NAdOs.mwfn，在主功能0里觀看NAdO軌道，其中本征值最大的6個NAdO軌道如下所示，等值面數值用的0.03 從上圖中可以看到這些軌道在苯與Pd(100)之間都是相位相同方式疊加的，必然都是起到明顯成鍵作用的，所以本征值都不小。從軌道圖形上可以看出這些軌道來自于苯上C原子的垂直于苯環的p軌道與Pd的原子軌道混合，而且有的圖里直接就能清楚看出Pd用的是d原子軌道。例如NAdO(2)中和苯接觸的Pd明顯用的是dz2軌道，從Pd區域的軌道等值面形狀就能看出這一點。感興趣的話可以進一步按照《談談軌道成份的計算方法》（http://www.shanxitv.org/131）介紹的用SCPA方法做一下軌道成分分析。在Multiwfn主菜單里輸入 8  //軌道成分分析 3  //SCPA方法 2  //2號NAdO軌道 Multiwfn馬上就輸出了軌道成份，下面列出的是Multiwfn返回的各個原子上各個角動量基函數產生的貢獻。明顯可以看到Pd主要都是用D基函數，而C用的都是P基函數。因此對NAdO(2)的軌道成份的分析體現了Pd的d原子軌道與C的p原子軌道的顯著混合對Pd吸附苯有關鍵性貢獻。可見以這種方式討論成鍵特征能分析得巨清楚透徹！ Composition of each shell Shell  1370 Type: D    in atom  196(Pd) :     0.50345 % Shell  1405 Type: D    in atom  201(Pd) :     0.88530 % Shell  1440 Type: D    in atom  206(Pd) :     0.50251 % Shell  1538 Type: D    in atom  220(Pd) :     0.80192 % Shell  1545 Type: D    in atom  221(Pd) :     0.76004 % Shell  1570 Type: S    in atom  225(Pd) :     0.53572 % Shell  1573 Type: D    in atom  225(Pd) :    10.40186 % Shell  1577 Type: S    in atom  226(Pd) :     0.56745 % Shell  1580 Type: D    in atom  226(Pd) :    10.21554 % Shell  1608 Type: D    in atom  230(Pd) :     0.80255 % Shell  1615 Type: D    in atom  231(Pd) :     0.76058 % Shell  1713 Type: D    in atom  245(Pd) :     0.53415 % Shell  1717 Type: S    in atom  246(Pd) :     0.96339 % Shell  1720 Type: D    in atom  246(Pd) :    21.79157 % Shell  1748 Type: D    in atom  250(Pd) :     0.53493 % Shell  1752 Type: S    in atom  251(Pd) :     0.96056 % Shell  1755 Type: D    in atom  251(Pd) :    21.80729 % Shell  2019 Type: P    in atom  289(C ) :     1.48815 % Shell  2024 Type: P    in atom  290(C ) :     1.47584 % Shell  2029 Type: P    in atom  291(C ) :     6.18539 % Shell  2034 Type: P    in atom  292(C ) :     1.44475 % Shell  2039 Type: P    in atom  293(C ) :     1.41651 % Shell  2044 Type: P    in atom  294(C ) :     6.20361 % 感興趣的讀者還可以用CDA方法從苯的分子軌道與Pd表面的晶體軌道的混合角度來考察Pd對苯吸附造成的電子轉移和軌道相互作用的細節。詳見《使用Multiwfn結合CP2K對周期性體系做電荷分解分析（CDA）》（http://www.shanxitv.org/716）。 3 總結 本文通過COF二維層狀體系和Pd(100)表面吸附苯作為例子，演示了如何用Multiwfn程序計算原子間和片段間的Mayer鍵級和模糊鍵級，還演示了如何計算pi鍵級單獨考察pi電子的貢獻，以及講解了如何做NAdO分析以軌道形式清楚直觀地考察原子間或片段間的共價作用的內在特征。可見使用Multiwfn可以把成鍵特征情況展現得超級清楚透徹。本文的做法對其它情況的周期性體系，諸如三維原子晶體、過渡態結構等等，也都是完全適用的。更多細節和相關知識請參看本文提到的相關博文。 Multiwfn對周期性體系能算的鍵級不止本文示例的這些。對周期性體系Multiwfn還能算Mulliken鍵級和Wiberg鍵級，還能把Mulliken和Mayer鍵級分解成不同軌道的貢獻，后者稱為軌道占據數擾動的Mayer鍵級。具體介紹見《Multiwfn支持的分析化學鍵的方法一覽》（http://www.shanxitv.org/471），做法見Multiwfn手冊4.8節的相關例子，本文就不特意演示了。除了本文介紹的這些外Multiwfn還有很多其它考察成鍵的方法對周期性體系都可以用，比如AIM拓撲分析，見《使用Multiwfn結合CP2K做周期性體系的atom-in-molecules (AIM)拓撲分析》（http://www.shanxitv.org/717）。 使用Multiwfn計算的結果若用于發表，記得需要按照Multiwfn啟動時的提示恰當引用Multiwfn的原文。給別人代算的話也必須明確告知對方這一點。 http://www.shanxitv.org/717 2024-07-03T11:37:00+08:00 使用Multiwfn做周期性體系的atom-in-molecules (AIM)拓撲分析 Using Multiwfn to perform atom-in-molecules (AIM) topology analysis for periodic systems 文/Sobereva@北京科音   2024-Jul-3 0 前言 Bader提出的非常知名的atom-in-molecules (AIM)理論在分析孤立體系和周期性體系的電子結構方面有非常廣泛的應用，相關資料匯總見《AIM學習資料和重要文獻合集》（http://bbs.keinsci.com/thread-362-1-1.html）。AIM分析框架中最常見的分析手段之一是電子密度的拓撲分析，也常被稱為AIM拓撲分析，例如根據這種拓撲分析得到的鍵臨界點（BCP）的位置可以用來討論成鍵和弱相互作用特征，參考《Multiwfn支持的分析化學鍵的方法一覽》（http://www.shanxitv.org/471）和《Multiwfn支持的弱相互作用的分析方法概覽》（http://www.shanxitv.org/252）里的相關介紹。非常強大的波函數分析程序Multiwfn（http://www.shanxitv.org/multiwfn）很早以前就已經支持了結合Gaussian、ORCA等量子化學程序產生的波函數文件對孤立體系做AIM拓撲分析，見《使用Multiwfn做拓撲分析以及計算孤對電子角度》（http://www.shanxitv.org/108）和Multiwfn手冊4.2節的大量例子。如今Multiwfn也已支持了對周期性體系做AIM拓撲分析，本文的目的是對此進行演示，以使得讀者可以輕松地舉一反三充分將這種手段應用于實際研究中。使用Multiwfn做AIM拓撲分析在發表文章時記得需要按照Multiwfn啟動時的提示恰當引用Multiwfn的原文。 讀者務必使用2024-Jun-16及以后更新的Multiwfn版本，否則情況與本文所述不符。Multiwfn可以在主頁http://www.shanxitv.org/multiwfn免費下載。不了解Multiwfn者建議閱讀《Multiwfn入門tips》（http://www.shanxitv.org/167）和《Multiwfn FAQ》（http://www.shanxitv.org/452）。本文要用到非常流行、高效且免費的第一性原理程序CP2K，不熟悉者推薦通過“北京科音CP2K第一性原理計算培訓班”（http://www.keinsci.com/workshop/KFP_content.html）系統性學習。本文要使用Multiwfn創建CP2K輸入文件，這在《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）里有簡要介紹。本文用的CP2K是2024.1版。 下面的例子涉及到的所有文件都可以在http://www.shanxitv.org/attach/717/file.rar里得到。 1 實例1：NaCl晶體 這一節對NaCl做AIM拓撲分析，要使用CP2K產生的molden文件作為Multiwfn的輸入文件。molden文件的產生方法在《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》（http://www.shanxitv.org/651）里有詳細說明。由于考慮k點時CP2K沒法產生molden文件，因此必須使用足夠大的超胞模型并只考慮gamma點。本文文件包里NaCl.cif是NaCl單胞的晶體結構文件，基于它我們創建CP2K算單點任務的超胞的輸入文件，同時得到molden文件。 用Multiwfn載入NaCl.cif，然后輸入 cp2k  //創建CP2K輸入文件的界面 [回車]  //產生的輸入文件名用默認的NaCl.inp -2  //要求產生molden文件 -11  //幾何操作界面 19  //構造超胞 2  //第1個方向平移復制為原先的2倍 2  //第2個方向平移復制為原先的2倍 2  //第3個方向平移復制為原先的2倍 -10  //返回 0  //產生輸入文件 現在當前目錄下就有了NaCl.inp，對應的是PBE/DZVP-MOLOPT-SR-GTH計算級別，對此體系是合適的。當前晶胞邊長為11.28埃，對于絕緣體來說只考慮gamma點基本夠用。用CP2K運行NaCl.inp，得到NaCl-MOS-1_0.molden。用文本編輯器打開NaCl.inp，將以下內容插入文件開頭來定義晶胞信息從而使得Multiwfn在分析時能考慮周期性，并且定義Na和Cl的價電子數分別為9和7（不了解為什么這么設置的人看前述的《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》）。  [Cell] 11.28112000     0.00000000     0.00000000  0.00000000    11.28112000     0.00000000  0.00000000     0.00000000    11.28112000  [Nval]  Na 9  Cl 7 下面開始做AIM拓撲分析。啟動Multiwfn，載入NaCl-MOS-1_0.molden，然后輸入 2  //拓撲分析（默認是對電子密度做拓撲分析） 2  //以每個原子核位置為初猜搜索臨界點 當前從屏幕上會看到找到了64個(3,-3)臨界點，分別對應64個原子核位置的電子密度極大點。之所以雖然當前用的是贗勢基組，在原子核位置還能找到電子密度極大點，是因為Multiwfn自動對這些原子指認了合適的EDF信息來描述內核電子，詳見《在贗勢下做波函數分析的一些說明》（http://www.shanxitv.org/156）。 然后再輸入 3  //以每兩個原子間連線中點作為搜索臨界點的初猜位置 4  //以每三個原子的中心位置作為搜索臨界點的初猜位置 0  //觀看結果 此時文本窗口顯示了各種臨界點的數目，如下所示，并且提示Poincare-Hopf關系已經滿足了，因此大概率所有臨界點都已經找全了。  The number of critical points of each type:  (3,-3):    64,   (3,-1):   384,   (3,+1):   384,   (3,+3):    64  Poincare-Hopf verification:    64  -   384  +   384  -    64  =     0  Fine, Poincare-Hopf relationship is satisfied, all CPs may have been found 在Multiwfn的圖形窗口可以看到下圖，棕紅色大球是Na，綠色大球是Cl，橙色、黃色、綠色小球分別是電子密度的(3,-1)、(3,+1)和(3,+3)臨界點，而(3,-3)臨界點都被原子球擋住了。在Multiwfn手冊3.14節對Multiwfn的拓撲分析功能有完整、非常詳細的說明，建議閱讀。 下面再把AIM拓撲分析里常涉及的鍵徑產生出來。在圖形界面右上角點Return關閉之，然后選8產生(3,-3)和(3,-1)之間的拓撲路徑，也即鍵徑。在筆者的i9-13980HX CPU上24核并行很快就產生完了。之后再選0進入圖形界面觀看，看到下圖，可見所有鍵徑（棕色細線）都產生了出來，既連接相鄰的Na和Cl，也連接每一對臨近的兩個Cl。 默認情況下Multiwfn把處于晶胞邊緣的所有鏡像原子和臨界點也顯示了出來，如果不想顯示它們的話（也即只顯示晶胞中唯一的那些），在圖形界面頂端依次選擇Other settings - Toggle showing all boundary atoms和Toggle showing all boundary CPs and paths。在文本窗口也會看到當前設置的狀態。 下面來考察一下Na和Cl之間的BCP的屬性。用選項0的圖形界面右邊的選項把其它臨界點關閉顯示而只顯示(3,-1)臨界點，選擇CP labels復選框要求顯示臨界點序號，關閉鍵徑顯示并要求顯示原子，然后就能看到下圖，這是當前圖像的一個邊緣區域（點擊save picture按鈕可在當前目錄下得到高分辨率圖像文件，下圖是剪切出的一部分） 可見189、512、238等序號的臨界點都對應Na-Cl的BCP，考察哪個都可以，都是等價的。點擊Return按鈕關閉圖形窗口，然后選擇7進入臨界點屬性計算界面，然后再輸入臨界點序號189，此時就看到了這個臨界點處各種函數的值，如下所示。這些函數在Multiwfn手冊2.6節都有簡要介紹，且在“量子化學波函數分析與Multiwfn程序培訓班”（http://www.keinsci.com/workshop/WFN_content.html）里有十分全面深入詳細的講解。 CP type: (3,-1)  Density of all electrons:  0.1144132659E-01  Density of Alpha electrons:  0.5720663296E-02  Density of Beta electrons:  0.5720663296E-02  Spin density of electrons:  0.0000000000E+00  Lagrangian kinetic energy G(r):  0.1168592222E-01  G(r) in X,Y,Z:  0.7296255010E-03  0.1022667133E-01  0.7296253875E-03  Hamiltonian kinetic energy K(r): -0.3025433987E-02  Potential energy density V(r): -0.8660488230E-02  Energy density E(r) or H(r):  0.3025433987E-02  Laplacian of electron density:  0.5884598270E-01  Electron localization function (ELF):  0.1993303578E-01  Localized orbital locator (LOL):  0.1249064363E+00 ...略 如我在全面提到的《Multiwfn支持的分析化學鍵的方法一覽》文中所述，BCP位置上電子密度拉普拉斯函數和能量密度都為正是離子鍵的典型特征。當前BCP這兩個值分別為0.0588 a.u.和0.00302 a.u.，都為正值，完全體現出Na-Cl是離子鍵的實事。此外，當前BCP的電子密度僅為0.01144 a.u.、ELF僅為0.0199，BCP處很小的電子密度和ELF值也是典型的離子鍵的普遍共性。 我在《使用Multiwfn+VMD快速地繪制高質量AIM拓撲分析圖（含視頻演示）》（http://www.shanxitv.org/445）介紹了一種極好的將Multiwfn和VMD聯用快速、理想地繪制AIM拓撲分析圖的做法，比在Multiwfn里直接顯示的效果更好，而且視角可以自由旋轉，沒看過的話務必先看一下。對于當前體系，也可以像此文一樣用Multiwfn文件包里的examples\scripts\AIM.bat和AIM.txt來作圖，得到的圖像如下。 上圖效果明顯不好，顯得殘缺不全，這是因為此時顯示的原子、臨界點、鍵徑都只有晶胞里唯一的那些，而處于棱和面上的邊界原子、臨界點、鍵徑的周期鏡像都沒有顯示出來。為了避免這個問題，對于當前體系有原子、臨界點、鍵徑處于棱和面上的情況，應當改用examples\scripts\目錄下的AIM+boundary.bat和AIM+boundary.txt來繪制，其用法與AIM.bat和AIM.txt的組合完全一致。只不過前者傳遞給Multiwfn的指令要求Multiwfn導出mol.pdb、CPs.pdb和paths.pdb文件時把邊界的鏡像也都納入進去。利用AIM+boundary.bat和AIM+boundary.txt來繪制的結果如下，可見此時的效果就很好了，和前面在Multiwfn里看到的完全相同。注意此時這三個pdb文件里只有序號靠前的那些才是晶胞中唯一的那些。 2 層狀COF 下面再舉個例子，用Multiwfn對一個層狀的共價有機框架化合物（COF）體系做AIM拓撲分析。此體系的cif文件是本文文件包里的COF-12000N2.cif，其結構圖如下所示。 由于此體系的晶胞尺寸不小，計算時可以直接對原胞做只考慮gamma點的計算來得到molden文件。如《實驗測定分子結構的方法以及將實驗結構用于量子化學計算需要注意的問題》（http://www.shanxitv.org/569）所說，由于X光衍射測的晶體結構中氫的原子位置通常不準確，最好做分析前先凍結重原子而優化一下氫原子的坐標，但此例為了省事就跳過這一步了。如果你要用AIM拓撲分析考察氫鍵，由于氫的位置是重中之重，那顯然是絕對不能略過這一步的。 當前體系用PBE結合CP2K的GTH贗勢基組計算完全可以，而此例作為演示，這回使用PBE結合pob-TZVP-rev2全電子基組計算。啟動Multiwfn，載入COF-12000N2.cif，然后輸入 cp2k  //創建CP2K輸入文件的界面 [回車]  //產生的輸入文件名用默認的COF-12000N2.inp -2  //要求產生molden文件 2  //選擇基組 15  //pob-TZVP-rev2 4  //開OT，對當前用的基組通常比對角化明顯更容易收斂 0  //產生輸入文件 用CP2K計算COF-12000N2.inp，然后將以下晶胞設置信息插入到計算得到的COF-12000N2-MOS-1_0.molden文件的開頭  [Cell]  22.55600000     0.00000000     0.00000000 -11.27800000    19.53406901     0.00000000   0.00000000     0.00000000     6.80000000 啟動Multiwfn，載入COF-12000N2-MOS-1_0.molden，然后輸入 2  //拓撲分析（默認是對電子密度做拓撲分析） 2  //以每個原子核位置為初猜位置搜索臨界點 3  //以每兩個原子間連線中點作為搜索臨界點的初猜位置 4  //以每三個原子的中心位置作為搜索臨界點的初猜位置 8  //產生(3,-3)到(3,-1)的拓撲路徑 0  //觀看結果 現在看到下圖。可見臨界點和鍵徑都成功地得到了，不僅有化學鍵對應的BCP和鍵徑，COF內氫鍵N-H...O作用以及層之間相互作用相對應的BCP和鍵徑都有。而且不僅有晶胞內兩層COF之間相互作用的BCP和鍵徑，當前晶胞里的COF和相鄰鏡像晶胞里的COF相互作用對應的BCP和鍵徑也都有，顯然周期性完全充分考慮了。 雖然Multiwfn的文本窗口提示當前并未符合Poincare-Hopf關系，但由于感興趣的臨界點和鍵徑都有了，所以就沒必要進一步搜索了。 再利用前述名為AIM+boundary.bat的Windows批處理腳本結合AIM+boundary.txt里記錄的指令，用Multiwfn+VMD繪制體系結構+臨界點+鍵徑圖。然后再在VMD的文本窗口里輸入pbc box命令顯示盒子，VMD main界面里選擇Display - Orthographic用正交視角。此時VMD顯示的圖像如下，可見非常理想！ 順帶一提，對于當前體系要想完整展現層間相互作用的話，最理想的做法是用我提出的IGMH，見《使用Multiwfn做IGMH分析非常清晰直觀地展現化學體系中的相互作用》（http://www.shanxitv.org/621）和《一篇最全面介紹各種弱相互作用可視化分析方法的文章已發表！》（http://www.shanxitv.org/667）。 3 對Ag晶胞基于電子密度格點數據做AIM拓撲分析 Multiwfn不僅可以基于CP2K產生的波函數文件做AIM拓撲分析，還可以基于從cube等文件里載入的電子密度的格點數據做拓撲分析，此時各個位置的電子密度、梯度和Hessian矩陣基于格點數據以B-spline方式插值得到。由于這個特性，Multiwfn做周期性體系的AIM拓撲分析不限于CP2K用戶，諸如從VASP產生的CHGCAR文件里讀取電子密度格點數據做AIM拓撲分析也可以。Multiwfn支持哪些格點數據文件格式見手冊2.5節的說明。而且CP2K用戶還可以在考慮k點的情況下對原胞計算得到電子密度cube文件結合Multiwfn做AIM拓撲分析，而不必在晶胞較小時非得弄成超胞模型再算。但基于電子密度格點數據做AIM拓撲分析的關鍵不足是對臨界點只能得到電子密度及其導數信息以及一些與之直接相關的量，而依賴于波函數才能計算的如動能密度、ELF等函數就沒法計算了。另外，如果電子密度格點數據是在贗勢下計算產生的，基于它插值出的電子密度函數中，原子核處將沒有電子密度極大點，也因此無法產生鍵徑。 此例演示用Multiwfn對Ag的單胞基于電子密度cube文件做AIM拓撲分析的過程。使用的級別是PBE/TZVPP-MOLOPT-PBE-GTH-q19，這個基組是CP2K 2024.1的data目錄下的BASIS_MOLOPT_UZH基組文件里的。之所以不用常用的BASIS_MOLOPT里的Ag的基組，在于那里面Ag的基組都是-q11的，即只描述11個價電子。使用贗勢的情況下做AIM拓撲分析時，被描述的價電子越多，即使用越小核贗勢，結果越準確、越可能接近全電子計算的結果。 啟動Multiwfn，載入本文文件包里的Ag單胞的晶體結構文件Ag.cif，然后輸入 cp2k  //創建CP2K輸入文件的界面 [回車]  //產生的輸入文件名用默認的Ag.inp -3  //要求產生cube文件 1  //電子密度 6  //開smearing 8  //設置k點 10,10,10 0  //產生輸入文件 把Ag.inp里BASIS_SET_FILE_NAME設為BASIS_MOLOPT_UZH，POTENTIAL_FILE_NAME設為POTENTIAL_UZH。Ag的BASIS_SET設為TZVPP-MOLOPT-PBE-GTH-q19，POTENTIAL設為GTH-PBE-q19。然后用CP2K計算，之后會得到記錄晶胞內價電子密度的Ag-ELECTRON_DENSITY-1_0.cube。如果不了解cube格式的話，看《Gaussian型cube文件簡介及讀、寫方法和簡單應用》（http://www.shanxitv.org/125）。用文本編輯器打開這個cube文件，在第一行里寫上box2cell。這樣Multiwfn載入這個cube文件時，一旦發現第一行的內容包含box2cell，就會把格點數據的盒子信息轉化為晶胞信息，從而在計算的時候能夠考慮周期性（也可以不修改cube文件，而是載入cube文件后用隱藏的主功能1000里的子功能18來實現盒子信息向晶胞信息的轉化）。 把Multiwfn的settings.ini里的iuserfunc設為-3，使得user-defined function對應基于載入的格點數據通過B-spline插值的函數。然后啟動Multiwfn，載入Ag-ELECTRON_DENSITY-1_0.cube，之后輸入 2  //拓撲分析 -11  //修改被分析的函數 100  //User-defined function，當前對應于插值方式產生的價層電子密度 3  //將每兩個原子之間的中點作為搜索臨界點的初猜 然后進入選項0，并把原子標簽顯示出來、把鍵設細，就會看到下圖。可見臨近原子之間的BCP都已經找到了 下面再把更多的臨界點也找出來。當前情況不適合選擇選項4、5分別用3、4個原子中心作為初猜點來搜索臨界點。選項3用每兩個原子之間中點當做初猜點時是考慮循環相鄰晶胞中的鏡像原子的，而選項4、5則只會循環當前晶胞里那些唯一的原子，因此此例用選項4、5搜索后還是會遺漏許多重要的臨界點。當前最適合的方式是選擇選項6里的子選項-1，此時會在每個原子附近特定半徑內撒一批初猜點來搜索臨界點。這么選過之后，會發現找到了一大堆新的臨界點，進入選項0后看到的圖如下所示。下圖左側把所有臨界點都顯示了出來，可以看到離原子核較近的區域臨界點有不少，而且各種類型的都有。這是因為當前用了贗勢，由于電子密度不是從原子核開始向四周單調下降的，因此會在價層區域出現各種各樣的臨界點，這些不是我們感興趣的。把原子球設大掩蓋這些沒什么意義的臨界點，并且旋轉體系使得視角與晶軸平行后，就看到了下面的右圖。可見在原子間相互作用區域找到的臨界點的分布整齊、對稱，因此可以認為有意義的臨界點都找到了。 有一個概念叫平面性指數，在“量子化學波函數分析與Multiwfn程序培訓班”（http://www.keinsci.com/workshop/WFN_content.html）里有一頁幻燈片做了介紹，如下所示。我們這里來對Ag晶體計算一下這個平面性指數，看看是否和文章里說的情況吻合。注意這里說的平面性指數和《使用Multiwfn定量化和圖形化考察分子的平面性（planarity）》（http://www.shanxitv.org/618）里介紹的我提出的衡量體系幾何結構平面性的參數是兩碼事。 首先我們得確定要考察的籠臨界點也即(3,+3)的臨界點編號。下圖左邊只顯示了籠臨界點，可見當前體系有兩類籠臨界點，一種是諸如晶胞中心的那個（15號），它處于八面體中心。另一種是處于四面體中心的，比如63號臨界點。下圖右邊只顯示了所有BCP，此體系中所有的BCP都是等價的，我們考慮其中任意一個，比如1號。利用當前拓撲分析界面里的選項7依次對這三個臨界點的屬性進行考察，可以得到的這三個位置的基于格點數據插值產生的電子密度，也即User-defined real space function后面的值，結果如下（注意不要讀Multiwfn直接輸出的Density of electrons后面的值，那個是用孤立狀態原子的電子密度疊加得到的準分子密度） CP(15): 0.1572533472E-01 a.u. CP(63): 0.2212960297E-01 a.u. CP(1): 0.2706526462E-01 a.u. 可見兩類籠臨界點里電子密度最小值是0.01572 a.u.，BCP處的電子密度為0.02706 a.u.，因此平面性指數為0.01572/0.02706=0.58，在0.5左右，符合幻燈片里說的“其它金屬及合金”的情況。 4 總結&其它 本文通過三個有代表性的例子詳細介紹了在Multiwfn中對周期性體系做電子密度的拓撲分析，也即一般說的AIM拓撲分析的完整過程，所有要注意的細節都做了解釋。通過例子可見，Multiwfn做周期性體系的拓撲分析的過程很簡單，并不比分析孤立體系更復雜。讀者如果實際操作一遍，還會感受到Multiwfn做這些分析的速度相當快。AIM拓撲分析對于考察周期性體系的電子結構有重要意義，推薦讀者在實際研究中運用，使文章增光添彩。 很值得一提的是，Multiwfn的拓撲分析功能極其普適，絕不僅限于電子密度的拓撲分析。對于Multiwfn自身可以直接計算的函數，以及載入的格點數據里記錄的函數，Multiwfn都可以做拓撲分析，對分子體系的例子參看比如《使用Multiwfn對靜電勢和范德華勢做拓撲分析精確得到極小點位置和數值》（http://www.shanxitv.org/645）對靜電勢、范德華勢的拓撲分析，以及《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）里面對ELF-pi函數的拓撲分析。對周期性體系，Multiwfn也完全可以以相同的方式對任意函數做拓撲分析。 使用本文的做法做AIM拓撲分析在發表文章時必須按照Multiwfn啟動時的提示恰當引用Multiwfn，如果用于給別人代算也要明確告知對方這一點。 http://www.shanxitv.org/716 2024-07-01T22:17:00+08:00 使用Multiwfn結合CP2K對周期性體系做電荷分解分析（CDA） Using Multiwfn in combination with CP2K to perform charge decomposition analysis (CDA) for periodic systems 文/Sobereva@北京科音  2024-Jul-1 1 前言 筆者之前寫過《使用Multiwfn做電荷分解分析(CDA)、繪制軌道相互作用圖》（http://www.shanxitv.org/166）介紹了怎么用Multiwfn程序基于Gaussian等量子化學程序產生的波函數文件對分子體系做電荷分解分析（CDA）。這種方法可以從片段軌道間相互作用的角度深入了解片段間是如何轉移電子的。由于CDA的很高價值，Multiwfn做CDA分析已經得到非常廣泛的使用。如果讀者沒看過此http://www.shanxitv.org/166的話一定要先仔細看一下，否則無法理解下文的內容。 從2024-Jun-5更新的Multiwfn開始，CDA模塊已經支持了基于CP2K產生的molden文件做周期性體系的CDA分析，下面將通過一個簡單的例子進行演示。Multiwfn可以在主頁http://www.shanxitv.org/multiwfn免費下載，不了解此程序者參看《Multiwfn FAQ》（http://www.shanxitv.org/452）和《Multiwfn入門tips》（http://www.shanxitv.org/167）。我假定讀者已經具備了CP2K的基本使用知識，不了解者推薦通過北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）完整系統地學習。 本文的例子是計算NaCl板吸附CO分子的體系，用CDA方法考察CO與NaCl板之間電子轉移情況。這個體系優化后的結構如下，C和O分別是1和2號原子，3-110號原子是NaCl板。對這個體系，之前我在《使用CP2K結合Multiwfn繪制密度差圖、平面平均密度差曲線和電荷位移曲線》（http://www.shanxitv.org/638）中通過密度差圖的方式定性考察了電子轉移，在《使用Multiwfn對周期性體系計算Hirshfeld(-I)、CM5和MBIS原子電荷》（http://www.shanxitv.org/712）里通過片段電荷的方式定量考察了電子轉移。而下面用CDA分析，可以以軌道角度在明顯更深層次理解電子轉移細節。 本文例子涉及的所有文件都在http://www.shanxitv.org/attach/716/file.7z里。CP2K用的版本是2024.1，Multiwfn是2024-Jun-5更新的版本，讀者絕對不要用Multiwfn的更老版本。 2 產生輸入文件 首先我們要用CP2K計算NaCl+CO、NaCl、CO各自的molden格式的波函數文件。如果你不知道怎么用CP2K產生molden文件的話，先仔細閱讀《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》（http://www.shanxitv.org/651），我假定讀者已經看過本文了。這里有幾個要點： (1)整體和片段的molden文件里的原子坐標必須對應 (2)整體和片段的計算必須使用嚴格相同的級別和設置，用的盒子也應當一致 (3)必須要求CP2K計算所有空軌道 (4)整體的molden文件里需要加入[Cell]字段定義晶胞信息，而片段的molden文件里加不加入不影響CDA結果。[Nval]字段可加入可不加入，也不影響CDA結果 本文文件包里的NaCl_CO.cif是以前我在PBE-D3(BJ)/DZVP-MOLOPT-SR-GTH下優化的NaCl板吸附CO的結構文件。用Multiwfn載入它，然后輸入 cp2k  //創建CP2K輸入文件 NaCl_CO.inp  //產生的輸入文件名 -2  //要求產生molden文件 -9  //其它設置 12  //設置計算的空軌道數 -1  //計算所有空軌道 0  //返回 0  //產生輸入文件 用CP2K運行NaCl_CO.inp，得到NaCl_CO-MOS-1_0.molden。然后將以下字段插入到molden文件的開頭 [Cell] 16.92168000     0.00000000     0.00000000  0.00000000    16.92168000     0.00000000  0.00000000     0.00000000    25.00000000 將NaCl_CO.inp復制為NaCl.inp，再將NaCl.inp中&COORD里的CO部分刪掉、PROJECT NaCl_CO改名為PROJECT NaCl，然后用CP2K運行之，得到NaCl-MOS-1_0.molden。將上面的[Cell]加入其中（這和CDA分析無關，加入這個的目的是為了之后在Multiwfn中能正確觀看NaCl板的周期性的晶體軌道）。 將NaCl_CO.inp復制為CO.inp，再將CO.inp中&COORD里的NaCl部分刪掉、PROJECT NaCl_CO改名為PROJECT CO，然后用CP2K運行之，得到CO-MOS-1_0.molden。這個可以不用刻意加上[Cell]，因為當前體系中CO離盒子邊界很遠，顯示軌道時是否考慮周期性對看到的圖像沒影響。 接下來就可以進行CDA分析了。 3 進行CDA分析 啟動Multiwfn，載入NaCl_CO-MOS-1_0.molden，然后依次輸入 16  //CDA分析 2  //兩個片段 CO-MOS-1_0.molden  //CO的波函數文件。注意要按整個體系里原子序號順序載入片段 NaCl-MOS-1_0.molden  //NaCl板的波函數文件 現在屏幕上輸出了一大堆復合物軌道的CDA結果，同時還輸出了ECDA結果： The net electrons obtained by frag. 2 = CT( 1-> 2) - CT( 2-> 1) =    0.1580 這個ECDA結果告訴你CO（片段1）往NaCl板（片段2）轉移了0.158個電子，這和《使用Multiwfn對周期性體系計算Hirshfeld(-I)、CM5和MBIS原子電荷》（http://www.shanxitv.org/712）里通過Mulliken方法計算的CO的片段電荷嚴格相同。 當前CDA分析輸出的復合物軌道太多，一共437個（對應整個體系的占據軌道數），而絕大多數對d、b、r項的貢獻微乎其微或完全為0，沒有考察的意義。為了簡化輸出，在當前界面里輸入 -3  //設置CDA輸出時d、b、r項的閾值，它們之一的絕對值大于這個閾值的復合物軌道才會輸出 0.005 0  //顯示CDA和ECDA結果 當前結果如下    Orb.      Occ.          d           b        d - b          r     272    2.000000    0.006220   -0.000012    0.006232    0.000488     275    2.000000    0.125862   -0.001874    0.127736    0.051155     297    2.000000    0.001098    0.000305    0.000792   -0.005373     318    2.000000    0.001053    0.000326    0.000727   -0.005720     341    2.000000    0.000781    0.000200    0.000581   -0.006708     352    2.000000    0.002291    0.000731    0.001560   -0.015705     360    2.000000    0.001101    0.000298    0.000803   -0.009297     414    2.000000    0.000635    0.000214    0.000421   -0.005229     427    2.000000    0.000739    0.000219    0.000520   -0.005731 ------------------------------------------------------------------- Sum:     874.000000    0.153805    0.017732    0.136074   -0.037831 Note: The "Sum" includes all terms including those not printed above CDA的總結果是CO向NaCl轉移了d=0.154個電子，NaCl向CO反饋了b=0.018個電子，NaCl凈獲得d-b=0.136個電子。b項基本可以忽略不計，而d的主要貢獻是由275號復合物軌道造成的。這個復合物軌道的貢獻顯然值得進一步探究，看看是由CO和NaCl板的哪些片段軌道混合所造成的。為此，輸入 6  //將特定復合物軌道對CDA項的貢獻進行分解 275  //要考察的復合物軌道序號 0.005  //輸出閾值 現在看到如下結果 Occupation number of orbital   275 of the complex:  2.00000000 FragA Orb(Occ.)  FragB Orb(Occ.)      d           b        d - b          r    5( 2.0000)     271( 2.0000)    0.000000    0.000000    0.000000    0.005183    5( 2.0000)     289( 2.0000)    0.000000    0.000000    0.000000    0.005642    5( 2.0000)     336( 2.0000)    0.000000    0.000000    0.000000    0.005127    5( 2.0000)     341( 2.0000)    0.000000    0.000000    0.000000    0.008960    5( 2.0000)     435( 0.0000)    0.009358    0.000000    0.009358    0.000000    5( 2.0000)     439( 0.0000)    0.011934    0.000000    0.011934    0.000000    5( 2.0000)     443( 0.0000)    0.010929    0.000000    0.010929    0.000000    5( 2.0000)     447( 0.0000)    0.008032    0.000000    0.008032    0.000000    5( 2.0000)     453( 0.0000)    0.009177    0.000000    0.009177    0.000000    5( 2.0000)     465( 0.0000)    0.010755    0.000000    0.010755    0.000000    5( 2.0000)     469( 0.0000)    0.013735    0.000000    0.013735    0.000000    5( 2.0000)     473( 0.0000)    0.011282    0.000000    0.011282    0.000000 Sum of above terms:               0.085201    0.000000    0.085201    0.024911 可以看到CO向NaCl轉移電子靠的都是CO的5號占據軌道和NaCl空軌道相互作用。NaCl板接收CO轉移來的電子分散在其大量空軌道上，比如NaCl板的443號空軌道接受了0.011個電子。在當前的輸出閾值下，上面輸出的這些項的所有的d項加和為0.085，明顯和275號復合物軌道的總d值0.153有很大差距，說明CO還向上面沒輸出的很多其它的NaCl板的空軌道轉移了電子。 下圖把275號復合物軌道、CO的5號占據軌道，以及NaCl接收CO貢獻來的電子最多的兩個軌道（439、469）展示在了一起，便于大家直觀了解情況，等值面數值用的是0.02。復合物和片段的軌道可以用Multiwfn載入相應的波函數文件后用主功能0按照《使用Multiwfn觀看分子軌道》（http://www.shanxitv.org/269）說的觀看。從此圖可以清晰直觀地認識到CO的孤對電子往吸附CO的那個Na位點的空軌道上轉移了電子。還有大量也具有這種特征的NaCl板的空軌道也接收了CO的孤對電子，這里沒全都畫出來。 下面再看一下上圖中275號復合物軌道是怎么由片段軌道構成的。按0返回到CDA菜單，然后進入選項2，輸入275，看到以下信息  Note: Only the fragment orbitals with contribution >  1.0 % will be shown below, the threshold can be changed by "compthresCDA" in settings.ini  Occupation number of orbital   275 of the complex:  2.00000000  Orbital     5 of fragment  1, Occ: 2.00000    Contribution:   87.42 %  Sum of values shown above:     87.42 % 這告訴你275號復合物軌道由CO的5號占據軌道貢獻了87.4%。默認情況下只有貢獻大于1%的片段軌道才會輸出，雖然NaCl板的一大堆非占據軌道都通過與CO的5號軌道混合而參與了這個復合物軌道的構成，但貢獻都小于輸出閾值，所以上面沒看到。從當前信息可以認識到275號復合物軌道主要體現CO的5號軌道特征，這是為什么上面圖中275號復合物軌道和CO的5號軌道看起來十分相似。 4 總結 此文通過一個簡單的表面吸附的例子演示了怎么用Multiwfn的CDA功能分析周期性體系的電荷轉移本質。討論電子轉移常見的套路就是繪制電子密度差、算算片段電荷，而從本文的例子看到Multiwfn可以從片段軌道相互作用對電子轉移的本質進行深入的考察，掌握了本文介紹的方法明顯可以給分析增光添彩。 使用Multiwfn做本文介紹的分析請記得在發表文章時引用Multiwfn程序啟動時提示的Multiwfn原文，以及引用介紹CDA/GCDA方法的文章（Multiwfn用的實質上是我對CDA廣義化后的GCDA形式），這在進入Multiwfn的CDA功能時屏幕上明確提示了。 http://www.shanxitv.org/712 2024-06-08T15:43:00+08:00 使用Multiwfn對周期性體系計算Hirshfeld(-I)、CM5和MBIS原子電荷 Using Multiwfn to calculate Hirshfeld(-I), CM5 and MBIS atomic charges for periodic systems 文/Sobereva@北京科音 First release: 2024-Jun-8    Last update: 2024-Dec-27 1 前言 原子電荷（atomic charge）即原子所帶的凈電荷，對于討論化學體系中原子的帶電狀態有重要意義，作為原子的一種描述符它也和原子的很多屬性有密切聯系。強烈建議閱讀《一篇深入淺出、完整全面介紹原子電荷的綜述文章已發表！》（http://www.shanxitv.org/714）里提到的筆者寫的原子電荷的綜述以全面了解原子電荷。 強大的波函數分析程序Multiwfn（http://www.shanxitv.org/multiwfn）支持非常豐富的原子電荷計算方法，其中多數已經支持周期性體系。Multiwfn支持的方法中有一類是基于模糊式原子空間劃分，也就是每個原子對應一個平滑的權重函數，權重函數與體系的電子密度的乘積的全空間積分就是這個原子帶的電子數（原子的布居數），原子的核電荷數減去它就是原子電荷。Multiwfn支持的這類方法包括Hirshfeld、Hirshfeld-I、MBIS、Becke，以及對Hirshfeld電荷做后校正得到的ADCH和CM5電荷。從2024-May-25更新的Multiwfn版本開始，Hirshfeld、CM5、Hirshfeld-I、MBIS都已經支持了周期性體系，本文將具體介紹怎么用Multiwfn與非常流行、高效且免費的第一性原理程序CP2K相結合非常方便快速地計算這些原子電荷。Multiwfn計算這些電荷的功能是普適的，并不限于結合CP2K，也可以基于其它第一性原理程序產生的體系的電子密度的cube文件來算，還支持VASP產生的記錄電子密度的CHGCAR文件。 注：Multiwfn更老的一些版本也支持基于周期性波函數計算Hirshfeld和CM5電荷，但那時候用的算法對于周期性體系效率極低，本文介紹的做法與之有天壤之別。 下文第2節簡要介紹一下Hirshfeld、Hirshfeld-I、CM5和MBIS原子電荷的基本特點，第3節介紹用Multiwfn+CP2K計算它們的具體方法，第4節給出具體例子，同時還會介紹怎么計算片段電荷。使用Multiwfn計算原子電荷在發表文章時記得需要按Multiwfn啟動時的提示恰當引用Multiwfn。 Multiwfn可以在官網http://www.shanxitv.org/multiwfn免費下載，如果對Multiwfn不了解的話建議看《Multiwfn FAQ》（http://www.shanxitv.org/452）和《Multiwfn入門tips》（http://www.shanxitv.org/167）。如果你用CP2K創建給Multiwfn算原子電荷用的輸入文件的話，筆者假定你對CP2K的使用已經有基本且正確的了解，不具備這些知識的話非常建議通過北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）完整、系統地學習一遍。本文的例子利用到Multiwfn創建CP2K的輸入文件，相關介紹見《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）。 2 Hirshfeld、Hirshfeld-I、CM5和MBIS等原子電荷簡介 Hirshfeld、Hirshfeld-I、CM5和MBIS原子電荷在Multiwfn手冊3.9節的相應小節里都有很詳細的介紹，因此這里不介紹細節，只是簡單說一下它們的特點。 Theor. Chim. Acta (Berl.), 44, 129 (1977)中提出的Hirshfeld電荷是被廣為使用的原子電荷，計算方式簡單，物理意義清楚，在量子化學研究孤立體系和第一性原理研究周期性體系的領域中都用得很多。Hirshfeld電荷為人詬病的地方是原子電荷數值整體嚴重偏小，對偶極矩、靜電勢重現性很差，這在我寫的《原子電荷計算方法的對比》（http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml）里的對比測試中有充分體現。 JCP, 126, 144111 (2007)提出的Hirshfeld-I是在Hirshfeld基礎上的改進，通過迭代方式不斷更新原子權重函數直到收斂，這使得原子權重函數所描述的原子空間可以響應實際化學環境。其原子電荷數值顯著大于Hirshfeld，對靜電勢的重現性也有明顯改進。代價是Hirshfeld-I需要做迭代，往往幾十輪，計算耗時明顯高于Hirshfeld，而且需要事先提供體系中各種元素的各個氧化態的原子徑向密度信息，實現起來較麻煩（為了實現Hirshfeld-I，我計算了360多個原子的徑向密度文件并在Multiwfn中自帶）。 務必注意，雖然CP2K程序自己也支持Hirshfeld-I電荷的計算，但我測試發現其結果和一般意義的Hirshfeld-I電荷不同，我檢查了其源代碼也確認了其實現方式明顯不是基于Hirshfeld-I的標準定義，因此如果你要得到一般意義的Hirshfeld-I電荷，必須用Multiwfn來算。CP2K算的Hirshfeld電荷在結合SHAPE_FUNCTION DENSITY選項時的結果雖然不離譜，但結果和Multiwfn算的也存在不可忽視的差異，應當以Multiwfn的結果為準。 JCTC, 8, 527 (2012)中提出的CM5電荷和我在J. Theor. Comput. Chem., 11, 163 (2012)中提出的ADCH電荷一樣都是對Hirshfeld電荷的后校正，對大部分情況都能使得它的原子電荷變得更大、更符合化學直覺，對靜電勢的重現性也明顯變得更好。ADCH沒有CM5那么強的經驗性、不牽扯一堆經驗參數，因而原理更為理想、更推薦使用。CM5最初是用于分子體系的，如今在應用于晶體方面也已經有了一定探索，例如JCTC, 16, 5884 (2020)將CM5與其它一些原子電荷對比考察了它們對一些固體的電荷分布的描述。CM5還被用于一些計算化學領域的方法，比如《比SMD算溶解自由能更好的溶劑模型uESE的使用》（http://www.shanxitv.org/593）里介紹的uESE溶劑模型依賴于CM5電荷，CM5乘以1.2后適合結合OPLS-AA力場做動力學模擬，見《計算適用于OPLS-AA力場做模擬的1.2*CM5原子電荷的懶人腳本》（http://www.shanxitv.org/585）。 ADCH電荷在分子體系的研究方面已被廣為使用，之前有人問我是否會把ADCH電荷擴展到周期性體系。我目前沒這個打算，因為ADCH電荷的思想對周期性體系往往不適用。ADCH方法將原子偶極矩展開為周圍原子的校正電荷，但對于諸如NaCl這樣的體系，每個原子所處的環境是中心對稱的，因此原子偶極矩為0、校正電荷都為0，故ADCH電荷會與Hirshfeld電荷完全相同。 MBIS (Minimal Basis Iterative Stockholder)電荷于JCTC, 12, 3894 (2016)中提出。類似于Hirshfeld-I，此方法也是通過迭代方式（因此耗時較高）優化原子空間直到收斂。MBIS比Hirshfeld-I的一大好處是不需要利用各個元素各個氧化態的原子徑向密度信息，因此實現起來簡單得多，也更優雅。原文里的測試體現出MBIS比Hirshfeld-I整體更具有優勢。但MBIS經常收斂較慢，對于某些體系甚至需要200多輪才能充分收斂。目前Multiwfn里MBIS最高支持到Rn元素。 根據我的測試，對于大多數固體體系，原子電荷大小的關系是MBIS和Hirshfeld-I最大、誰大不一定，CM5大小介中，Hirshfeld電荷最小。如果你對數值的絕對大小不那么在乎，只是想體現相對大小關系的話，最“樸素”且很便宜的Hirshfeld也可以用。但如果對定量數值關心的話，建議用其它的。CM5、Hirshfeld-I和MBIS用于晶體哪個最理想沒有定論，都可以算一算看看，如果傾向于得到較大數值的結果可以優先考慮Hirshfeld-I和MBIS。如果其中有的原子電荷明顯不合理那就棄掉換別的方法，比如O的電荷如果算出來-2.34那就不應該用，因為顯然O最多只能帶-2電荷。 特別要強調的是，絕對不要把原子電荷與氧化態搞混，也切勿用原子電荷判斷氧化態！氧化態怎么正確地判斷看《使用Multiwfn結合CP2K計算晶體中原子的氧化態》（http://www.shanxitv.org/711）和《使用Multiwfn通過LOBA方法計算氧化態》（http://www.shanxitv.org/362）。 還有很多其它原子電荷計算方法，諸如Mulliken電荷、Lowdin電荷、SCPA電荷、AIM電荷（被個別文章稱為Bader電荷）等，Multiwfn也都支持將它們用于周期性體系，這不屬于本文的范疇。Mulliken電荷結合CP2K里常用的MOLOPT系列和pob系列基組用于周期性體系也說得過去。雖然它是基于希爾伯特空間劃分的原子電荷中最原始、思想最樸素的一種，但實際結果一般還算定性正確，而且計算耗時極低。被用得很多的AIM電荷實際上很糟糕，往往出現違背化學常識的結果，這在前述的《原子電荷計算方法的對比》里已經充分體現了，不建議用。 3 Hirshfeld、Hirshfeld-I、CM5和MBIS原子電荷在Multiwfn中的計算方法 3.1 CP2K用戶的情況 首先說對于CP2K用戶怎么計算。Hirshfeld、Hirshfeld-I、CM5和MBIS原子電荷都是完全基于電子密度計算的，不直接牽扯到波函數，在Multiwfn中計算它們有兩種可以用的輸入文件： (1)用CP2K產生的體系的電子密度cube文件作為輸入文件。這種情況在CP2K計算時可以考慮k點，是我最推薦的，耗時非常低。如果不了解cube文件的格式的話，參看《Gaussian型cube文件簡介及讀、寫方法和簡單應用》（http://www.shanxitv.org/125）。注意CP2K有bug（起碼直到我撰文時最新的CP2K 2024.1為止），產生的cube文件里的原子有效電荷都為0，這會導致Multiwfn無法正確計算原子電荷。這里說的原子的有效核電荷是指當前基組描述的原子的電子數，對于贗勢基組它對應價電子數。為此，要么用文本編輯器手動改cube文件中每個原子的有效核電荷，要么把cube文件第一行寫成各個元素的有效核電荷，例如N 5 B 3 Ti 12，這代表當前體系中的N、B、Ti的有效核電荷數分別為5、3、12 (2)用CP2K產生的記錄了體系波函數的molden文件作為輸入文件。詳見《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》（http://www.shanxitv.org/651）。這個方法的缺點在于沒法考慮k點，對于小晶胞體系必須先擴胞到足夠大，從而在計算時只考慮gamma點，這顯然使得在CP2K中和Multiwfn中的計算耗時都會比用(1)的方式高很多 用Multiwfn載入上述輸入文件其中一種后，進入主功能7，然后選擇相應選項即可計算Hirshfeld、Hirshfeld-I、CM5和MBIS原子電荷。對于molden文件當輸入文件時，Multiwfn還會讓你輸入計算用的格點間距，格點間距越小耗時越高而計算精度越高，通常用0.2 Bohr就可以，是精度和耗時的較好權衡。而使用cube文件做為輸入文件時，計算原子電荷用的格點間距與此文件的格點間距直接對應，CP2K的CUTOFF設得越大則產生的cube文件里的格點間距會越小。 由于Multiwfn計算上述電荷用的是均勻分布的格點，不可能準確積分較重原子的內核區域電子密度，因此上述輸入文件應當是CP2K做GPW計算得到的，也需要用贗勢基組，此時cube文件記錄的電子密度或者molden文件里記錄的波函數只對應于價電子。 順帶一提，如果你要計算Mulliken電荷，用Multiwfn載入上述的molden文件，然后進主功能7，依次選擇5、1即可得迅速到結果。Multiwfn中計算AIM電荷的方法筆者以后會另文說明。 3.2 其它程序用戶的情況 如果你是VASP用戶，可以用記錄晶胞中價電子密度的CHGCAR文件給Multiwfn用于計算前述原子電荷，等同于3.1節說的第2類輸入文件。只要文件名里包含CHGCAR字樣就會被Multiwfn當做是CHGCAR文件來載入。注意為了讓Multiwfn得知各類原子的有效核電荷數，CHGCAR文件的第一行需要以Nval為開頭，后面依次寫上各類元素的名字和有效核電荷，并以空格分隔。比如寫為Nval Na 1 Cl 7，這代表Na和Cl的有效核電荷數分別為1和7。一個具體例子是http://www.shanxitv.org/attach/712/CHGCAR_Si8.rar。 注意VASP的PAW計算可能在原子核附近區域產生負值的電子密度。從2024-Dec-27更新的Multiwfn開始可以對這種情況給出正確的Hirshfeld、Hirshfeld-I、CM5電荷，但這種情況的電子密度原理上不兼容MBIS計算，所以算MBIS電荷時應當不用PAW。 如果你是其它第一性原理程序的用戶，若程序可以產生記錄晶胞中價電子密度的cube文件（即等同于3.1節說的第2類輸入文件），也都可以用于給Multiwfn計算前述的原子電荷。 之后在Multiwfn中的計算過程和CP2K用戶沒任何區別。 4 實例 下面給出Multiwfn結合CP2K計算Hirshfeld、Hirshfeld-I、CM5和MBIS原子電荷的具體例子。相關的文件都可以在http://www.shanxitv.org/attach/712/file.zip中得到（4.2、4.3和4.5節的cube/molden文件除外，因為文件過大）。 4.1 TiO2晶體 此例對TiO2的金紅石晶型計算前述的原子電荷。它的實驗晶體結構的cif文件是本文文件包里的TiO2-Rutile.cif，結構圖如下所示 先用CP2K計算出這個晶胞的電子密度的cube文件，用于之后Multiwfn做原子電荷計算。首先用Multiwfn創建CP2K輸入文件，啟動Multiwfn并載入TiO2-Rutile.cif，然后輸入 cp2k  //創建CP2K的輸入文件 [回車]  //產生的文件名用默認的TiO2-Rutile.inp -3  //產生cube文件 1  //電子密度 8  //設置k點 4,4,6  //與各個晶胞邊長相乘皆為18埃左右，夠用了 0  //產生輸入文件，計算級別為默認的PBE/DZVP-MOLOPT-SR-GTH 用CP2K運行剛得到的TiO2-Rutile.inp，馬上就算完了，得到了TiO2-Rutile-ELECTRON_DENSITY-1_0.cube文件，在本文的文件包里已經提供了。用文本編輯器打開此文件，把第一行改為Ti 12 O 6并保存。這代表Ti和O的有效核電荷分別為12和6，直接對應于輸入文件里它們分別用的DZVP-MOLOPT-SR-GTH-q12和DZVP-MOLOPT-SR-GTH-q6基組名里面的q后面的值。 啟動Multiwfn，輸入 TiO2-Rutile-ELECTRON_DENSITY-1_0.cube  //寫此文件的實際路徑 y  //從第一行載入各個元素有效核電荷 7  //布居分析與原子電荷計算 1  //Hirshfeld電荷 立刻看到如下Hirshfeld電荷計算結果。由于格點間距不是無窮小，因此等價的原子算出來的原子電荷也會有輕微的數值差異，可以自己手動取平均。如果用更大的CUTOFF，則cube文件的格點間距會更小，不等價性會更低。 Final atomic charges: Atom    1(Ti):     0.59349459 Atom    2(O ):    -0.29461903 Atom    3(O ):    -0.29472601 Atom    4(Ti):     0.58363165 Atom    5(O ):    -0.29385236 Atom    6(O ):    -0.29385236 現在Multiwfn還問你是否把原子電荷導出為chg文件，這里選n不導出。 接下來計算CM5電荷。在Multiwfn主功能7里選擇16，馬上就得到了結果                       ---------- CM5 charges ---------- Atom:    1Ti  CM5 charge:    1.348393  Hirshfeld charge:    0.593495 Atom:    2O   CM5 charge:   -0.672093  Hirshfeld charge:   -0.294619 Atom:    3O   CM5 charge:   -0.672191  Hirshfeld charge:   -0.294726 Atom:    4Ti  CM5 charge:    1.338580  Hirshfeld charge:    0.583632 Atom:    5O   CM5 charge:   -0.671306  Hirshfeld charge:   -0.293852 Atom:    6O   CM5 charge:   -0.671306  Hirshfeld charge:   -0.293852 Summing up all CM5 charges:     0.00007649 可見CM5電荷的大小明顯大于Hirshfeld電荷，看起來更合理，而Hirshfeld電荷值明顯偏小。 下面再演示計算Hirshfeld-I電荷。把Multiwfn自帶的examples目錄下的atmrad目錄挪到當前目錄下，這樣Multiwfn計算Hirshfeld-I電荷時就會自動用此目錄下的各種元素各個價態的徑向電子密度。之后在主功能7里選15，再選1用默認的設置進行Hirshfeld-I電荷計算，經過16輪迭代得到了結果： Final atomic charges: Atom    1(Ti):     3.01789464 Atom    2(O ):    -1.50888824 Atom    3(O ):    -1.50886131 Atom    4(Ti):     3.01787620 Atom    5(O ):    -1.50897240 Atom    6(O ):    -1.50897240 可見Hirshfeld-I電荷的大小不僅顯著大于Hirshfeld電荷，比CM5電荷還要大很多。雖然數值很大，但完全在合理范圍內，畢竟都沒有大過原子的氧化態。 接下來再計算MBIS電荷。在主功能7里選擇20，再選1用默認的設置進行MBIS電荷計算，經過19輪迭代得到了結果： Final atomic charges: Atom    1(Ti):     1.73473636 Atom    2(O ):    -0.86733169 Atom    3(O ):    -0.86731253 Atom    4(Ti):     1.73467041 Atom    5(O ):    -0.86734303 Atom    6(O ):    -0.86734303 對于當前例子，MBIS的數值介于CM5和Hirshfeld-I之間。但也有很多反例，MBIS也經常比Hirshfeld-I電荷還大不少，后文有例子。 以上原子電荷的計算順序是隨意的，不要誤以為必須按上述順序計算，你想算哪個就算哪個即可。 4.2 AlN晶體 這一節計算AlN晶體的原子電荷。完全可以按照和上一節相同的做法基于電子密度的cube文件來計算，但本節演示一下基于超胞的molden文件的計算過程。AlN原胞的cif文件是本文文件包里的AlN.cif，將之載入Multiwfn，然后輸入以下命令創建CP2K輸入文件 cp2k  //創建CP2K輸入文件 AlN_443.inp  //產生的輸入文件名 -11  //進入幾何操作界面 19  //構造超胞 4  //a方向擴胞成原本的4倍 4  //b方向擴胞成原本的4倍 3  //c方向擴胞成原本的3倍 -10  //返回 -2  //要求產生molden文件 0  //產生輸入文件 用CP2K運行剛產生的AlN_443.inp，得到AlN_443-MOS-1_0.molden。在此文件開頭插入以下內容定義晶胞和各個元素的價電子數 [Cell] 12.44400000     0.00000000     0.00000000 -6.22200000    10.77682012     0.00000000  0.00000000     0.00000000    14.93400000 [Nval] Al 3 N 5 啟動Multiwfn，載入AlN_443-MOS-1_0.molden，然后輸入 7  //布居分析與原子電荷計算 16  //CM5電荷 [回車]  //使用均勻分布的格點計算 [回車]  //使用默認的0.2 Bohr格點間距。格點間距越小結果越精確但計算越耗時 在i9-13980HX CPU上24核并行，50秒就算完了，結果為                       ---------- CM5 charges ---------- Atom:    1Al  CM5 charge:    0.467945  Hirshfeld charge:    0.302813 Atom:    2N   CM5 charge:   -0.467945  Hirshfeld charge:   -0.302812 Atom:    3Al  CM5 charge:    0.467945  Hirshfeld charge:    0.302812 Atom:    4N   CM5 charge:   -0.467946  Hirshfeld charge:   -0.302814 ...略 Al的電負性很小，N的電負性很大，因此二者的原子電荷理應相差懸殊。當前的Hirshfeld電荷明顯偏小。CM5電荷的數值雖然更大，但還是顯得有些偏低。 下面再計算Hirshfeld-I電荷。在主功能7里輸入 15  //Hirshfeld-I 1  //開始計算 [回車]  //使用默認的0.2 Bohr格點間距 經過37輪收斂，結果如下。可見Hirshfeld-I電荷的數量級比較合理，體現出AlN中成鍵的極強離子性  Atom    1(Al):     1.80317384  Atom    2(N ):    -1.80317397  Atom    3(Al):     1.80317257  Atom    4(N ):    -1.80317539 ...略 最后再計算下MBIS電荷。在主功能7里輸入 20  //MBIS 1  //開始計算 [回車]  //使用默認的0.2 Bohr格點間距 對此體系MBIS收斂很慢，經過190輪才達到收斂，在i9-13980HX上24核并行算了5分鐘，在雙路7R32 96核并行下算了不到一分鐘。結果如下。可見對此體系，MBIS電荷的數值很大，和氧化態幾乎正好是一樣的。所以，如前所述，若你就是想算出來比較大的原子電荷，可以先考慮MBIS，但也有一些例外，后面會提到。 Atom    1(Al):     2.99599978 Atom    2(N ):    -2.99601488 Atom    3(Al):     2.99600129 Atom    4(N ):    -2.99599935 4.3 MOF-5晶體 MOF-5晶體含有Zn、O、C、H。其cif文件是本文文件包里的MOF-5.cif。由于此體系的晶胞邊長約26埃，因此計算電子密度格點數據時不用考慮k點。按照3.1節的做法，基于PBE/DZVP-MOLOPT-SR-GTH級別的電子密度cube文件（第一行改為Zn 12 O 6 C 4 H 1），用Multiwfn計算各種原子電荷，結果如下。此體系中O有兩類，分別給出，C則處于較多不同的化學環境，其電荷就不給出了。 Hirshfeld： Zn=0.372 O=-0.332/-0.200 H=0.056 CM5： Zn=0.644 O=-0.584/-0.305 H=0.116 Hirshfeld-I（36輪收斂，2*7R32機子上96核并行4分多鐘算完） Zn=1.568 O=-1.917/-0.746 H=0.139 MBIS（39輪收斂，2*7R32機子上96核并行3分多鐘算完） Zn=1.546 O=-1.662/-0.782 H=0.190 對這個體系，MBIS和Hirshfeld-I相符得很好，而且Zn的電荷大小很合理。Hirshfeld電荷再次數值最小、CM5比之大一些。 4.4 其它：MoS2、BaTiO3、MgF2、NaCl 注意MBIS電荷并非總是較大。例如《使用Multiwfn結合CP2K計算晶體中原子的氧化態》（http://www.shanxitv.org/711）里的單層MoS2的例子，在PBE/DZVP-MOLOPT-SR-GTH級別的波函數下，Mo的電荷的計算結果是Hirshfeld=0.330，CM5=0.827，Hirshfeld-I=1.987，MBIS=-0.059，Mulliken=-0.771。對這個體系來說CM5和Hirshfeld-I顯得靠譜，MBIS電荷接近0有點說不通，Mulliken電荷明顯有誤導性。 BaTiO3在PBE/DZVP-MOLOPT-SR-GTH級別下的原子電荷計算結果如下。Hirshfeld還是嚴重偏小。Hirshfeld-I大得離譜，Ba的原子電荷都超過氧化態+2了。相較之下CM5和MBIS的結果比較靠譜，正好二者算的O的原子電荷差不多，差異主要在Ba和Ti的電荷的相對大小上。 Atom:    1Ba  CM5 charge:    1.103253  Hirshfeld charge:    0.404357 Atom:    2Ti  CM5 charge:    1.200467  Hirshfeld charge:    0.537212 Atom:    3O   CM5 charge:   -0.767963  Hirshfeld charge:   -0.310211 Atom:    4O   CM5 charge:   -0.767892  Hirshfeld charge:   -0.315693 Atom:    5O   CM5 charge:   -0.767892  Hirshfeld charge:   -0.315693 Hirshfeld-I Atom    1(Ba):     2.37169239 Atom    2(Ti):     3.02878350 Atom    3(O ):    -1.78799416 Atom    4(O ):    -1.80625420 Atom    5(O ):    -1.80625434 MBIS： Atom    1(Ba):     0.57391985 Atom    2(Ti):     1.69416825 Atom    3(O ):    -0.74678425 Atom    4(O ):    -0.76066530 Atom    5(O ):    -0.76066537 再看兩個典型離子化合物的情況，都在PBE/DZVP-MOLOPT-SR-GTH級別下計算。對NaCl計算的Na的原子電荷為：Hirshfeld=0.217，CM5=0.443，Hirshfeld-I：1.053，MBIS：1.043，Mulliken：0.596 對MgF2計算的Mg的原子電荷為： Hirshfeld=0.414，CM5=0.828，Hirshfeld-I：2.048，MBIS：1.888，Mulliken：1.439 可見對于強離子性化合物，Hirshfeld-I和MBIS電荷都在氧化態附近，并且由于方法定義的缺陷，往往電荷的大小還輕微大于氧化態的大小（但如果用大核贗勢，對Na和Mg只描述最外層的3s價電子，則不會有這個問題）。其它原子電荷的大小關系為Mulliken>CM5>Hirshfeld。 以MgF2為例，這里順帶一提CP2K算的Hirshfeld-I電荷明顯不對，千萬不要用。在CP2K輸入文件的&DFT/&PRINT里加上以下內容要求輸出Hirshfeld-I電荷后，CP2K直接給出的Mg的Hirshfeld-I電荷為0.522，和Multiwfn算的相距甚遠。 &HIRSHFELD   SHAPE_FUNCTION DENSITY   SELF_CONSISTENT T &END HIRSHFELD 4.5 NaCl板吸附CO的片段電荷 NaCl板吸附CO是《使用CP2K結合Multiwfn繪制密度差圖、平面平均密度差曲線和電荷位移曲線》（http://www.shanxitv.org/638）里的一個例子體系，這里基于這個模型計算一下原子電荷和CO片段的電荷。CP2K優化這個體系后產生的NaCl_CO-1.restart文件在本文的文件包里提供了，將之載入Multiwfn，然后輸入 cp2k  //創建CP2K輸入文件 NaCl_CO.inp  //產生的輸入文件名 -7  //設置周期性 XY -3  //要求產生cube文件 1  //電子密度 -2  //產生molden文件（用于之后算Mulliken電荷的目的） 0  //產生輸入文件 計算完畢后，把NaCl_CO-ELECTRON_DENSITY-1_0.cube文件第一行改為Na 9 Cl 7 C 4 O 6。先用Hirshfeld方法各個原子的電荷和CO的片段電荷。啟動Multiwfn，載入此cube文件，然后輸入 7  //布居分析與原子電荷計算 -1  //定義片段 109,110  //CO的原子序號 1  //Hirshfeld電荷 結果為 ...略  Atom  109(C ):     0.13622804  Atom  110(O ):    -0.03514734 Fragment charge:    0.10108070 Fragment population:    9.89891930 即CO的Hirshfeld片段電荷為0.101。 之后再選擇CM5方法，結果為 ...略  Atom:  109C   CM5 charge:    0.147306  Hirshfeld charge:    0.136228  Atom:  110O   CM5 charge:   -0.081972  Hirshfeld charge:   -0.035147  Summing up all CM5 charges:     0.00510030  Fragment charge:    0.06533361 之后再選擇Hirshfeld-I方法，結果為  Atom  109(C ):     0.09790054  Atom  110(O ):    -0.08744136  Fragment charge:    0.01045918 之后再選擇MBIS方法，結果為 ...略  Atom  109(C ):     0.19642171  Atom  110(O ):    -0.15937966  Fragment charge:    0.03704205 再來算Mulliken電荷。用Multiwfn載入加入了恰當的[Cell]和[Nval]信息的CP2K計算后產生的NaCl_CO-MOS-1_0.molden文件，輸入 7  //布居分析與原子電荷計算 -1  //定義片段 109,110  //CO的原子序號 5  //Mulliken電荷 1  //在屏幕上輸出結果 結果為 ...略 Atom   109(C )    Population:  3.60671601    Net charge:  0.39328399 Atom   110(O )    Population:  6.23484007    Net charge: -0.23484007 Total net charge:   -0.00000086 Fragment charge:    0.15844392 可見不同方法計算的結果在定量上有一定差異，但都有共性，即CO整體帶少量正電，即電子從CO往NaCl板上轉移了一些。并且C和O分別帶正電和負電，且原子電荷的絕對值都是C大于O。 5 總結 本文簡要介紹了四種知名的基于模糊式原子空間劃分計算原子電荷的方法的特點，包括Hirshfeld、Hirshfeld-I、CM5和MBIS，并介紹了它們在Multiwfn中的計算方法，并且給了很多具體例子，以使得讀者能了解如何使用Multiwfn結合CP2K或其它第一性原理程序很容易地計算它們。同時還說明了怎么直接計算片段電荷，這是定量討論體系中兩部分間電子轉移的最嚴格的方式。推薦將這些原子電荷應用于實際問題的研究當中以考察固體與表面體系的電荷分布情況。記得屆時應恰當引用Multiwfn啟動時提示的Multiwfn原文以及相應的原子電荷的原文。 通過本文的一些體系的對比可以明顯看到，Hirshfeld電荷原理上最簡單，但明顯缺點是整體嚴重偏小，CM5通常比它大不少，而且電荷的穩定性不錯。雖然結合化學直覺來看，CM5電荷的絕對大小仍然往往能偏小，但至少也不會明顯離譜，總比Hirshfeld電荷更建議使用。Hirshfeld-I和MBIS電荷普遍比較大，誰更大不一定，至少幾乎總比CM5還大不少，從絕對大小來看往往比CM5更合理，但穩健性比CM5弱一些，對極個別體系可能會表現得不符合常識。對于離子性特征很強的化合物，Hirshfeld-I和MBIS電荷往往接近氧化態，由于方法并不完美，有時候原子電荷甚至比氧化態還要大一點點。當然，沒有任何原子電荷計算方法是絕對完美的，本文介紹的這些原子電荷在實際當中都可以用，可以根據它們的特點和實際結果決定選用，在橫向對比時方法必須統一。 順帶一提，對于計算固體表面或者多孔物質的原子電荷以用于基于力場的動力學模擬的目的，我目前最推薦CP2K算REPEAT電荷，北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）里講“電子結構的分析”的部分做了專門的介紹。但REPEAT電荷需要在真空區定義擬合點，因而無法用于諸如本文所示的AlN、TiO2等致密的固體，而且對于遠離真空區的原子也沒法擬合出質量靠譜的原子電荷。所以REPEAT電荷有特定用處，但不是像本文介紹的那些原子電荷一樣屬于普適性的原子電荷計算方法。 http://www.shanxitv.org/711 2024-05-23T16:05:00+08:00 使用Multiwfn結合CP2K計算晶體中原子的氧化態 Using Multiwfn in combination with CP2K to calculate oxidation state for atoms in crystals 文/Sobereva@北京科音  2024-May-23 1 前言 波函數分析程序Multiwfn可以基于電子波函數計算原子和片段的氧化態，這在筆者之前寫的《使用Multiwfn通過LOBA方法計算氧化態》（http://www.shanxitv.org/362）一文中已做了詳細介紹，沒讀過者一定要先閱讀此文，否則無法理解下文的內容。此文介紹的做法已經被廣泛用來準確地考察分子中的氧化態。如今Multiwfn已經支持將這種方法用于CP2K計算的周期性體系的波函數，從而考察晶體中原子的氧化態，這非常有實際價值，具體做法將在本文介紹。注意氧化態和原子電荷完全是兩碼事，周期性體系的原子電荷的計算看《使用Multiwfn對周期性體系計算Hirshfeld(-I)、CM5和MBIS原子電荷》（http://www.shanxitv.org/712）。 Multiwfn可以在官網http://www.shanxitv.org/multiwfn免費下載，本文的讀者必須使用2024年5月23日及以后更新的版本，否則沒有本文介紹的特性。對Multiwfn不了解者看《Multiwfn FAQ》（http://www.shanxitv.org/452）和《Multiwfn入門tips》（http://www.shanxitv.org/167）。 2 用法 在Multiwfn中對周期性體系做LOBA或mLOBA分析氧化態的具體操作流程和對孤立體系做別無二致，都是先用主功能19產生定域化軌道（介紹見《Multiwfn的軌道定域化功能的使用以及與NBO、AdNDP分析的對比》http://www.shanxitv.org/380），然后進入主功能8（軌道成分分析），選擇LOBA/mLOBA分析，之后輸入一個LOBA方法判斷定域化軌道電子歸屬的閾值，或者輸入m讓Multiwfn按照mLOBA的方式判斷歸屬即可。定域化軌道的電子歸屬是什么意思，在前述的《使用Multiwfn通過LOBA方法計算氧化態》里已經專門說過了，這里不再累述。之后各個原子的氧化態就會立馬輸出出來。用戶也可以自定義片段來得到片段的氧化態。 對周期性體系做LOBA/mLOBA分析的輸入文件必須是CP2K產生的molden文件，并且需要自己在里面加入盒子信息。具體方法在《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》（http://www.shanxitv.org/651）里有詳細介紹。由于考慮k點時不能產生molden文件，因此原本晶胞較小的話需要擴胞成足夠大的超胞，從而對其計算時只需要考慮gamma點。Multiwfn不支持Quantum ESPRESSO、VASP等其它第一性原理程序產生的波函數，也沒有任何支持的必要，因為對于絕大多數情況，用Multiwfn創建CP2K單點輸入文件并用CP2K計算是非常簡單且快速的事，可參考《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）。 對于當前目的，CP2K計算時用贗勢基組也行，用全電子基組也行，沒有什么特別的講究。CP2K計算的時候不要開smearing，否則軌道占據數不再精確為整數、不再是標準的單行列式波函數，此時也就沒法做軌道定域化。如果有特殊原因非開smearing不可，可以在Multiwfn載入molden文件后進入主功能6，選擇選項38，此時Multiwfn會按照軌道能量由低到高重排電子確定占據數，軌道占據數就都為整數了，之后退回到主菜單，再照常按前述步驟做LOBA/mLOBA分析即可。 3 例子 下面給出一系列通過Multiwfn使用LOBA/mLOBA計算固體的氧化態的例子。CP2K的輸入輸出文件在http://www.shanxitv.org/attach/711/file.zip里都給了。只有3.1節例子對應的molden文件我在里面直接給了，如果所有例子的molden文件都給的話壓縮包就太大了。CP2K輸入文件里的各種設置我這里就不細說了，在北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）里有非常全面系統的講解。 3.1 單層氮化硼 氮化硼（BN）是知名的二維材料，這里完整介紹單層BN的氧化態計算的流程。本文文件包里的BN_bulk.cif是氮化硼晶體的結構文件，用GaussView打開它，可看到晶胞里有兩層BN。刪除其中的一層，然后保存為BN_ml.gjf。之后啟動Multiwfn，載入BN_ml.gjf，輸入 cp2k   //創建CP2K輸入文件 [回車]  //產生的輸入文件用默認的名字BN_ml.inp -2  //要求導出molden文件 -11  //進入幾何操作界面 19  //構造超胞 6  //第1個平行于層的方向擴胞成原先的6倍 6  //第2個平行于層的方向擴胞成原先的6倍 1   //垂直于層的方向不擴胞 -10  //返回 -7  //設置周期性 XY 0  //產生輸入文件（計算級別為默認的PBE/DZVP-MOLOPT-SR-GTH） 用CP2K運行Multiwfn產生的BN_ml.inp，我這里用雙路7R32 96核機子4秒鐘就算完了。現在當前目錄下有了BN_ml-MOS-1_0.molden，按照《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》（http://www.shanxitv.org/651）所述，把以下信息插入進去，以定義晶胞信息和原子價電子數信息。改好的文件已經在本文的文件包里提供了，是BN_ml-MOS-1_0.molden。  [Cell] 14.98944000     0.00000000     0.00000000 -7.49472000    12.98123580     0.00000000  0.00000000     0.00000000     8.00000000  [Nval] B 3 N 5 啟動Multiwfn，載入BN_ml-MOS-1_0.molden，然后輸入 19  //軌道定域化 1  //只定域化占據軌道，用默認的Pipek-Mezey方法 8  //軌道成分分析 100  //LOBA/mLOBA分析 輸入一個閾值，比如50，馬上看到各個原子的以LOBA方法計算的氧化態，如下所示。可見B的氧化態為3、N為-3，整個體系所有原子氧化態加和為0，非常符合常識。這里用的閾值50%對于成鍵方式比較典型的情況是適合的，它的意思是如果有某原子對某個定域化軌道貢獻超過50%（Multiwfn以Hirshfeld方法計算貢獻值），則這個軌道的電子就完全歸屬于這個原子。  Oxidation state of atom   1(B ) :  3  Oxidation state of atom   2(N ) : -3  Oxidation state of atom   3(B ) :  3  Oxidation state of atom   4(N ) : -3  Oxidation state of atom   5(B ) :  3  Oxidation state of atom   6(N ) : -3 ...略  The sum of oxidation states:     0 我提出的modified LOBA (mLOBA)方法是對LOBA的重要改進，它將每個定域化軌道的電子完全歸屬到對它貢獻最大的那個原子上，可以確保所有原子的氧化態加和總為0，而且也沒有閾值選擇的任意性（有的時候LOBA的結果對閾值敏感），比LOBA在原理上明顯更好。在當前界面里輸入m即得到mLOBA方法算的氧化態的結果，也是B和N的氧化態分別為3和-3。 3.2 BaTiO3 BaTiO3晶體的cif文件是本文文件包里的BaTiO3.cif。用Multiwfn創建其4*4*4超胞在PBE/DZVP-MOLOPT-SR-GTH下做單點的CP2K計算的輸入文件。當前的超胞在每個方向有16埃，足夠大了。CP2K計算后得到molden文件，在里面加入恰當的[Cell]和[Nval]字段后，啟動Multiwfn并將之載入，然后輸入 19  //軌道定域化 -9  //設置產生定域化軌道后自動計算軌道成份的方法 0  //不輸出軌道成份。這些信息不是當前感興趣的，當前體系又大，不輸出可以節約不少時間 1  //只定域化占據軌道，用默認的Pipek-Mezey方法 由于當前體系大、軌道數很多，因此軌道定域化過程耗時較高。在EPYC 7R32機子上花了約8分鐘。之后輸入 8  //軌道成分分析 100  //LOBA/mLOBA分析 m  //顯示mLOBA方法計算的原子氧化態 馬上看到以下結果。可見Ba、Ti、O的氧化態都和常識完全一致  Oxidation state of atom   1(Ba) :  2  Oxidation state of atom   2(Ti) :  4  Oxidation state of atom   3(O ) : -2  Oxidation state of atom   4(O ) : -2  Oxidation state of atom   5(O ) : -2 ...略 3.3 MOF-5 MOF-5是知名的金屬有機框架化合物，介紹見https://en.wikipedia.org/wiki/MOF-5。其中含有Zn，這里算算Zn的氧化態是多少。MOF-5晶體的cif文件在本文的文件包里提供了，是立方晶胞，邊長25.86埃。由于邊長已經大到可以只考慮gamma點做計算，因此CP2K計算時不需要擴胞，直接對原胞算單點得到molden文件即可。相應的CP2K輸入文件已提供在了本文的文件包里。氧化態計算過程同前，7R32機子上幾分鐘可以算完整個過程。氧化態結果為Zn=2、O=-2、H=1、C=2/0/2。可見Zn、O、H的氧化態都完全符合預期。當前體系中C的氧化態有三種，和所處的化學環境差異有關。 3.4 單層MoS2 MoS2是知名的二維材料體系，它的體相結構文件是本文文件包里的MoS2_bulk.cif。完全效仿3.1節氮化硼的例子，取單層并基于5*5*1的超胞模型做它的氧化態的計算。如果以LOBA方法計算氧化態并輸入50作為閾值，結果為  Oxidation state of atom   1(Mo) :  6  Oxidation state of atom   2(S ) : -2  Oxidation state of atom   3(Mo) :  6  Oxidation state of atom   4(S ) : -2  Oxidation state of atom   5(S ) : -2  Oxidation state of atom   6(S ) : -2  Oxidation state of atom   7(Mo) :  6 ...略  The sum of oxidation states:   50 雖然Mo的價電子有6個，看似Mo的氧化態為6也不是不可能，但實際上明顯是錯的，因為當前所有原子的氧化態加和為100。考慮到S的氧化態不可能比-2更負，明顯當前的結果高估了Mo的氧化態。即便嘗試各種閾值，也始終得不到能令LOBA方法算的氧化態加和為0的情況。因此LOBA對此體系完全失敗。 輸入m使用mLOBA計算氧化態，結果為  Oxidation state of atom   1(Mo) :  4  Oxidation state of atom   2(S ) : -2  Oxidation state of atom   3(S ) : -2  Oxidation state of atom   4(Mo) :  4  Oxidation state of atom   5(S ) : -2  Oxidation state of atom   6(S ) : -2  Oxidation state of atom   7(Mo) :  2  Oxidation state of atom   8(S ) : -2  Oxidation state of atom   9(S ) : -2 ...略  The sum of oxidation states:     0 可見mLOBA方法算的所有原子氧化態加和為期望的0，但是Mo有的氧化態是2有的是4，和當前體系中所有Mo等價的事實不符，原因后面會說。解決方法是輸入-1定義片段，然后輸入所有Mo對應的序號。序號的快速查詢方法是在Multiwfn主功能0的菜單里選擇Tools - Get atom indices of a given element，然后輸入Mo，點OK，就會返回下面的內容，將之粘貼到Multiwfn定義片段的窗口里即可。 1,4,7,10,13,16,19,22,25,28,31,34,37,40,43,46,49,52,55,58,61,64,67,70,73 片段定義好后，再次輸入m，可以看到Oxidation state of the fragment: 100，即這個片段總的氧化態為100。由于當前體系里有25個完全等價的Mo，因此Mo的氧化態為100/25=4，是完全正確的值。 為什么mLOBA算出來的明明等價的Mo會有不同的氧化態？這可以通過考察定域化軌道的特征結合mLOBA的原理予以了解。默認情況下，Multiwfn的主功能19做完軌道定域化后會自動輸出各個定域化軌道的成份，可以看到除了單中心、雙中心的定域化軌道外，還有大量三中心的：  More delocalized LMOs: (Three largest contributions are printed)   301:   64(Mo) 24.8%   46(Mo) 24.8%   61(Mo) 24.8%   302:   46(Mo) 24.8%   49(Mo) 24.8%   31(Mo) 24.8%   303:   16(Mo) 24.8%   34(Mo) 24.8%   31(Mo) 24.8% ...略 從軌道成份上明顯可知這些三中心軌道都是等價的。軌道定域化做完后按照《使用Multiwfn觀看分子軌道》（http://www.shanxitv.org/269）說的，在Multiwfn的主功能0里隨便看一個三中心軌道，如下所示 可見這個軌道體現的是三個Mo之間的三中心作用。按照mLOBA的思想，這個軌道上的兩個電子會完全歸屬到對它貢獻最大的原子上，但由于三個Mo的貢獻本質上相同，實際的貢獻量由于數值計算誤差會有極其輕微的差異，所以這個軌道上的電子歸屬到哪個Mo是有隨機性的，這是為什么Mo的氧化態有的為2有的為4。可見，只要把軌道定域化和mLOBA的思想都理解了，遇到可疑的結果自然就能通過具體分析搞明白是怎么回事、知道怎么恰當解決。 4 總結 本文詳細介紹了怎么用Multiwfn基于CP2K對固體做DFT計算得到的波函數通過LOBA和mLOBA方法得到原子的氧化態，給了很多例子。可見本文的方法操作簡單，結果非常合理，這對于討論原子在固體中的狀態、對體系進行分類很有價值。作為練習，讀者可以對NaCl、MgO、AlN、SiO2晶體也都算算氧化態。值得一提的，如本文的例子所體現的，如有了我后來提出的mLOBA方法，原本的LOBA就沒什么使用價值了，我建議總是默認用mLOBA。 使用本文的方法計算在發表結果時務必記得按照Multiwfn程序啟動時的提示對Multiwfn進行正確的引用。 http://www.shanxitv.org/706 2024-05-05T14:31:00+08:00 2024年計算化學公社舉辦的計算化學程序和DFT泛函的流行程度投票結果 Results of the Computational Chemistry Commune 2024 poll on the popularity of computational chemistry programs and DFT functionals 文/Sobereva@北京科音  2024-May-5 0 前言 2024年4月4號，在北京科音建立的人氣最高、專業性最強的綜合性計算化學論壇“計算化學公社”（http://bbs.keinsci.com）上開展了為期一個月的五項投票： 你最常用的做量子化學計算的DFT泛函投票（http://bbs.keinsci.com/thread-44387-1-1.html） 你最常用的做第一性原理計算的DFT泛函投票（http://bbs.keinsci.com/thread-44391-1-1.html） 你最常用的量子化學程序投票（http://bbs.keinsci.com/thread-44388-1-1.html） 你最常用的分子動力學程序投票（http://bbs.keinsci.com/thread-44389-1-1.html） 你最常用的第一性原理程序投票（http://bbs.keinsci.com/thread-44390-1-1.html） 現對投票結果進行總結和簡單評論。未來預計每三年重新開展一次投票。要強調的是，這個投票只是體現流行程度，和方法/程序的好壞并沒直接關系。本投票結果主要反映中國的計算化學領域的專業、內行群體的情況，不反映業余/外行群體的情況。本次投票的結果也有助于業余/外行研究者正確認清什么才是主流，減少被他人忽悠、聽信歪曲說辭誤入歧途的概率。 上一次的投票于2021年舉行，當時的結果和很多相關評論見下文： 2021年計算化學公社論壇“你最常用的計算化學程序和DFT泛函”投票結果統計 http://www.shanxitv.org/599（http://bbs.keinsci.com/thread-23482-1-1.html） 1 你最常用的做量子化學計算的DFT泛函投票 本次可投的泛函有43種，帶不帶色散校正算同一種泛函。在2021年的投票條目基礎上有所增減，特別是增加了雙雜化泛函，明顯幾乎不會有人用的泛函沒納入可投范圍。投票范疇僅限做量子化學計算的情況，不包含第一性原理計算的情況。關系特別近的，比如M05-2X和M06-2X、wB97X和wB97XD和wB97X-D3、SCAN和r2SCAN、revDSD-PBEP86-D3(BJ)和DSD-PBEP86-D3(BJ)等等當做同一個泛函來計。此次投票者共713人。本投票每個人最多選6項，且所投的泛函必須占平時全部研究工作的5%以上。按照得票率（票數除以總投票人數）繪制的圖如下。為了避免條目過多，只把得票前30名的列出。此圖中諸如某泛函對應50%就代表有50%的人平時較多使用此泛函，后文的統計圖同理。 總的來說本年度的投票結果和2021年時沒太大變化，前六名的順次沒有改變，還是依次為B3LYP、(M05/M06)-2X、PBE0、wB97X-(/D/D3)、PBE、CAM-B3LYP老幾樣，各自有各自的用處（參看我對2021年投票結果的評論http://www.shanxitv.org/599，這里不再贅述）。它們的得票率的相對比例也基本沒變，也就是量化領域里用處相對有限的PBE的比例有一定降低，以后肯定還得跌。2021年時候的第7名M06雖然得票率如今還是5%左右，但排名已下滑到第10，被wB97M-V、SCAN/r2-SCAN、PWPB95-D3(BJ)所超過。M06在我來看用處著實不大，雖然計算過渡金屬配合物體系有一定用處，但PBE0-D3(BJ)/D4多數情況是更好的選擇，而且M06還有后繼者MN15可用。wB97M-V的得票率從2018年的3.1%提升到了如今的10%，已經算是增幅很快了，再過3年統計時肯定還會增高。此泛函在國內量化研究者中一定程度的流行，很大程度在于在計算化學公社論壇和思想家公社QQ群的討論中時常被提及、在《簡談量子化學計算中DFT泛函的選擇》（http://www.shanxitv.org/272）博文中和我在北京科音基礎（中級）量子化學培訓班（http://www.keinsci.com/workshop/KBQC_content.html）中的推薦、被免費的ORCA程序支持。提出時間不算很長的純泛函SCAN及其改進版r2SCAN現在得票率能到6%著實不容易，2021年時得票率還不到1%，這主要在于有越來越多的程序已經支持此泛函，而且綜合素質整體強于更早的經典泛函PBE，因而在純泛函范疇內能有重要的位置。 2021年投票的時候沒納入雙雜化泛函，這次得票率超過1%的雙雜化泛函的排名順序是PWPB95-D3(BJ)(5.9%) > (rev)DSD-PBEP86-D3(BJ)(3.1%) > B2PLYP-D3(BJ) (2.7%) > ωB97X-2-D3(BJ) (2.0%)。PWPB95-D3(BJ)和(rev)DSD-PBEP86-D3(BJ)能在國內用戶中變得流行和上述wB97M-V的情況很類似。本身這倆泛函各有長處，有流行開來的必然性。PWPB95-D3(BJ)比較robust，算過渡金屬配合物能量問題較好，能在ORCA里用；而revDSD-PBEP86-D3(BJ)算主族體系反應能、勢壘以及弱相互作用能都是雙雜化里頂尖的，而且在Gaussian里通過《Gaussian中非內置的理論方法和泛函的用法》（http://www.shanxitv.org/344）我介紹的方法能直接用。此外，ORCA中DSD-PBEP86適合算TDDFT和NMR目的也是其加分項。這倆泛函現在流行度能超過經典且最早提出的雙雜化泛函B2PLYP是其應得的。 BLYP這次的排名降幅很大，從第10名已跌到第22名，本身這個泛函如今鮮有用武之地。BLYP一般也就算算主族體系，在Gaussian里用這個明顯不如用B3LYP，耗時也持平，而以前在ORCA里用這個搭配def2-SVP結合RIJ加速做有機體系幾何優化速度效率高是個用處，以前我在《大體系弱相互作用計算的解決之道》（http://www.shanxitv.org/214）里也提過，但如今也不如改用B97-3c來跑。 2 你最常用的做第一性原理計算的DFT泛函投票 可投的泛函有26種，帶不帶色散校正算同一種泛函。此投票在2021年沒有，是本次新加的。此次投票者共442人。本投票每個人最多選6項，且所投的泛函必須占平時全部研究工作的5%以上。結果如下，零票的未顯示 96年提出的PBE至今仍穩居第1的位置，如同B3LYP在量子化學領域的地位，而且和第二名相差更懸殊。PBE能有這樣的地位是必然的，PBE提出年代早、被程序支持得極為廣泛，計算便宜，有豐富的專門為其搞的贗勢，幾何優化和分子動力學目的大多數時候夠用（加色散校正后又拓寬了其普適性），而且在基態能量相關問題方面依然有使用價值而且沒被已流行的其它純泛函甩開特別多（這和B3LYP在量子化學領域的情況截然不同，B3LYP算能量早過時了，很難再發得出去文章，見http://bbs.keinsci.com/thread-12773-1-1.html）。B3LYP在這次投票里得了第2，令我有點意外，大概率是很多人不好好看投票帖子的說明，誤把量子化學研究用的泛函也在此進行投票了。PBE0能排第3不意外，需要一個HF成分適當的雜化泛函做TDDFT/TDDFPT算激發態、算強相關體系的問題時經常用得著。HSE06流行度排得上第4主要來自于其計算帶隙、能帶方面公認很好，以及晶胞參數優化方面表現不錯。PBEsol是優化晶體結構、晶胞參數的好把式，而且還是便宜的純泛函，能排到第5很正常。M06-2X能排第6令我有點意外，一方面必定是有人誤當成量子化學計算的情況來投，另一方面是計算主族晶體/液體相關的化學反應、吸附的相關能量問題必定有一些人在用。SCAN/r2SCAN已經有一定流行度，由于在文獻中出現頻率越來越高，在未來的流行度必定也會逐漸提升，但流行度超越PBE在可預見的未來還不太可能，畢竟對于大量PBE就已經表現得夠用的范疇，作為更貴但沒帶來顯著優勢的meta-GGA的SCAN/r2SCAN不可能顯著侵占這方面的市場。第一性原理領域里用BLYP的人我不很理解是什么心態。revPBE和與之相似的RPBE能有一定流行度在于算化學吸附方面不錯，被不少人用。第一性原理方面的泛函投票中CAM-B3LYP顯得遠不如量子化學領域里來得流行，估計用這個的大部分都是CP2K用戶用來算激發態和UV-Vis譜方面，只占投票的少量群體。算化學吸附很好的BEEF-vdW和算物理吸附很好的optB88-vdW能有一定票數不算意外。純泛函中矮子里拔將軍算帶隙往往可以接受的HLE17在本次投票中獲得了一點流行度，略意外的是算帶隙整體更好些的純泛函mBJ反倒在這次投票中顯得無人問津，可能是前者在CP2K里能直接用而后者不能所致。作為PBE后繼提出來的知名的TPSS流行度那么低有點令我意外，倒也確實實際用處不太大，現在又有了理論上以及實際整體表現得更好的r2SCAN。PW91雖然在文獻里出現得很多，但這次得票率相當低，畢竟實際中有PBE就基本沒有更老的PW91能派上用場的時候。B97M-rV能有少量票數，主要在于算熱力學量方面在純泛函里是表現得較突出的。 3 你最常用的量子化學程序投票 可投程序有29種，投票者共679人。本投票每個人最多選三項，且所投的程序必須占平時全部研究工作的10%以上。按照得票率繪制的圖如下，只顯示了得票前20名的 前三位和2021年投票的結果一樣，還是Gaussian > ORCA > xtb，Gaussian依然是約90%的量子化學研究者日常必用的程序，穩穩占據主導位置。而ORCA和xtb的得票率比2021年時都有了約10%增長，這是意料之中的。實際上這三個程序也是我自己最常用的：xtb拿來快速預優化和結合molclus（http://www.keinsci.com/research/molclus.html）做構象搜索的初篩，Gaussian做基于DFT的opt、freq、掃描、IRC等涉及幾何變化的任務以及算一些屬性（NMR、超極化率等），ORCA算高精度單點。這三個程序的流行度遠遠甩開了其它程序，一方面是它們比較容易安裝和使用，一方面各有各的獨特優勢，有被大量使用的剛性原因。而且它們把大部分量子化學計算的應用領域都給覆蓋了，對于日常應用性研究來說其它程序難以有和它們競爭的顯著空間。Dmol3、ADF、Q-Chem、Turbomole這四個商業程序日子愈發不好過。在量化計算方面沒有長處還巨貴的Dmol3的用戶越來越少，從2021年的6.2%已經進一步萎縮到了4.3%，以后肯定還會明顯進一步萎縮。ADF從2021年時候的僅僅2.2%進一步萎縮到了1.5%。Q-Chem從2021年的3.0%萎縮到了1.0%。Turbomole從2021年的1.6%萎縮到了1.0%。以GPU加速為賣點的TeraChem更不幸，2021年時候還有1人投票，今年變成了0人。 4 你最常用的第一性原理程序投票 可投程序有25種，投票者共565人。本投票每個人最多選三項，且所投的程序必須占平時全部研究工作的10%以上。按照得票率繪制的圖如下（0票的沒顯示） 根據這次投票的結果可見，至少在計算化學公社論壇里，CP2K的流行程度已經趕超了VASP。這令我90%程度感到意外，但也有10%程度感覺是在情理之中。由于Multiwfn在2021年開始提供了極其易用的創建CP2K輸入文件的功能（http://www.shanxitv.org/587），我后來又對此功能反復打磨并在Multiwfn中提供了對CP2K繪制DOS、能帶、STM、電子激發、成鍵分析等許多功能，再加上2023年3月、2023年10月和2024年3月開辦了三期北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）非常全面系統講解了CP2K的使用，無疑令中國的CP2K用戶猛增。但即便我已預料CP2K的得票率肯定會遠高于2021年時候的33.3%，但也沒預料到這次居然能達到60.9%，直接翻了將近一倍，甚至把一直霸占流行度榜首的VASP給擠下去了。由于免費且十分高效的CP2K的用戶還在不斷激增，而且CP2K更新很快，不斷完善和發展新功能，Multiwfn在未來還會給其提供更多的相關分析處理功能，CP2K的位置在以后無疑還會更加牢固。相比之下，VASP的流行度從2021年投票時候的65.8%降到了現在的58.1%。由于現在有非常強大的競爭者CP2K，而且CP2K不具備的一些功能還有免費的Quantum ESPRESSO能用來平替VASP，售價較昂貴且算大體系速度遠不如CP2K的VASP在未來的前景不樂觀。以前很多人一提到第一性原理計算就言必稱VASP，以后就不再是如此了。除了CP2K的流行度猛增外，其它程序的流行度都普遍出現了下降，如CASTEP和Dmol3分別從2021年的19.0%和9.3%下降到了13.8%和6.4%。Wien2k今年更是連一票都沒有，而2021年時還有3票。那些程序流行度百分比減少，一方面是它們的票數大多數確實有一定減少，用戶發現有更理想或免費的程序可用，另一方面原因應當是有很多通過CP2K程序開始入手第一性原理計算的人加入了投票，他們只對CP2K的得票率有貢獻而間接地拉低了其它程序的得票率。 5 你最常用的分子動力學程序投票 可投程序有18種，投票者共551人。本投票每個人最多選三項，且所投的程序必須占平時全部研究工作的10%以上。按照得票率繪制的圖如下 GROMACS依舊是用戶數的龍頭老大，而且流行度從2021年投票時的69.3%還進一步略微提升到了71.3%，得票數大約等于其它所有程序用戶數之和，和曾經的情況一樣。第2、3位依然分別是Lammps和AMBER。Lammps和OpenMM得票率略漲，而AMBER、Forcite和NAMD的流行度都有較多降幅，GULP、DL_POLY和CHARMM更是快跌沒了。 http://www.shanxitv.org/690 2023-11-07T17:53:00+08:00 CP2K做雜化泛函計算的關鍵要點和簡單例子 Key points and simple examples of hybrid functional calculations with CP2K 文/Sobereva@北京科音  First release: 2023-Nov-7   Last update: 2024-Mar-9 0 前言 對周期性體系，雜化泛函的計算耗時眾所周知遠高于純泛函。著名的第一性原理程序CP2K不光純泛函的計算超級快，對大體系做雜化泛函計算也非常快（盡管還是遠比純泛函貴），遠遠快于Quantum ESPRESSO (QE)、VASP等純粹基于平面波的程序。CP2K結合像樣的基組做雜化泛函的單點計算，在一般雙路服務器上算到上千原子都沒問題。然而，使用CP2K做高效的雜化泛函計算有很多非常關鍵性的要點必須知道，遠遠不是像Gaussian程序那樣光寫個雜化泛函名字就了事的。然而，筆者經常在網上回答CP2K問題時看到很多人完全不具備這些最基本常識，胡用瞎用CP2K做雜化泛函計算，導致耗時巨高、始終出不來結果、任務崩潰、SCF完全不收斂等各種問題。因此筆者覺得很有必要寫個小文將CP2K做雜化泛函計算的最關鍵常識和要點挨個提一下。Multiwfn程序（http://www.shanxitv.org/multiwfn）具有極其便利的創建CP2K輸入文件的功能，在本文的最后會還用Multiwfn產生C60晶體的CP2K雜化泛函計算的輸入文件作為實際例子。 本文的內容至少對CP2K 2024.1，以及2023-Nov-7及以后更新的Multiwfn是適用的。 筆者在北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）的“能量的計算及相關問題”部分里講雜化泛函計算遠比本文全面、深入得多，并列舉諸多經驗技巧和例子，還詳細介紹HSE06、CAM-B3LYP、wB97X等各種范圍分離形式的雜化泛函以及revDSD-PBEP86-D3(BJ)等各種雙雜化泛函的用法和要點，使用雜化泛函結合k點計算的要點，歡迎想學習更多知識者關注和參加。 1 CP2K雜化泛函計算的特性 --------- 2024-Feb-18更新：以下之前寫的關于k點的說法已經過時。從CP2K 2024.1版開始支持了雜化泛函計算時考慮k點，稱為RI-HFXk算法，并且可以無縫結合ADMM降低耗時。此時需要以RI近似計算HF交換部分。不僅能算能量還能做優化（包括變胞）和能帶計算。從第3屆的“北京科音CP2K第一性原理計算培訓班”開始已經加入了RI-HFXk的全面深入的介紹和具體例子，參考http://bbs.keinsci.com/thread-43683-1-1.html。 首先要知道CP2K做雜化泛函計算的一個關鍵不足是不支持考慮k點，因此如果被計算的是周期性體系，晶胞必須大到只考慮gamma點就足夠。如果你要計算的體系的晶胞原本比較小，顯然必須先擴胞到足夠大才行。如果不知道擴到多大夠用，嚴格來說需要做你感興趣的屬性的計算結果隨擴胞倍數的收斂性測試。雖然CP2K能夠對含有非常多原子的超胞做得動雜化泛函計算，但耗時會明顯高于QE等程序對原胞在考慮k點的情況下做雜化泛函計算。 最能充分體現CP2K雜化泛函計算效率優勢的情況是算一些較大分子晶體、MOF、COF等一些大晶胞體系，以及做第一性原理動力學等情況。對小晶胞體系如上所述經過擴胞后也能算，只不過不體現CP2K的長處。 CP2K顯然沒法用雜化泛函算能帶，因為不支持k點，這種目的只能用QE等程序代替。必須用CP2K算的話，應當用CP2K支持的純泛函里算帶隙相對比較好的比如HLE17。 --------- 由于雜化泛函遠比普通泛函昂貴，尤其是用于幾何優化、NEB、動力學等涉及到結構變化的任務。因此如果你算的問題用純泛函就足以描述得不錯、用雜化泛函不會帶來明顯好處的話，顯然應該用純泛函。記住，非必要甭考慮雜化泛函，純泛函的耗時僅有雜化泛函的一個微不足道的零頭。對幾何優化、振動分析、NEB等高耗時任務用純泛函，算勢壘、反應能的時候用更好的雜化泛函是可以接受且常見的策略。 周期性體系的雜化泛函計算不僅耗時遠高于純泛函計算，對內存要求更是遠遠高于純泛函計算。經常用CP2K做雜化泛函任務的服務器的內存容量絕對是多多益善，選擇服務器的時候要注意這一點。也別天真地試圖拿個人計算機將就著跑CP2K的雜化泛函計算。 CP2K做雜化泛函計算的前提是編譯時必須帶著libint電子積分庫，否則沒法算雜化泛函中HF交換項中涉及的雙電子積分。CP2K的編譯方法見http://www.shanxitv.org/586。 雜化泛函用的贗勢基組和贗勢就用常見的給PBE等純泛函優化的就可以。從CP2K 9.1開始，data目錄下的BASIS_MOLOPT_UZH和POTENTIAL_UZH文件里還分別給出了專門給PBE0優化的贗勢基組和贗勢，如C的DZVP-MOLOPT-PBE0-GTH-q4和GTH-PBE0-q4，原理上來說它們用于雜化泛函計算會更理想一些（但實際未必有什么優勢）。 2 庫侖截斷 雜化泛函計算之所以比純泛函昂貴得多，在于計算雜化泛函的HF交換項時要計算大量的雙電子積分（two-electron integral, ERI），尤其是對于周期性體系來說，不光要計算晶胞內基函數之間的ERI，還要計算中心晶胞與周圍鏡像晶胞基函數之間的ERI。再加上由于只能算gamma點而晶胞必然較大，原子數、基函數通常非常多，無疑要算的ERI是巨量的。 由于周期邊界條件，原理上來說要計算的ERI是無窮的，顯然這會導致計算無法進行。因此實際計算時必須使用庫侖截斷（Coulomb truncation），也就是將ERI中描述電子庫侖相互作用的1/r算符截斷在某個距離（截斷半徑），超過截斷半徑時1/r視為0，相應地ERI就不需要計算了，這個策略使得ERI的計算數目是有限的。如果你計算的是非周期性體系，由于本身ERI數目就是有限的，因此雖然也可以用庫侖截斷，但這就不是必須的了。 雜化泛函的庫侖截斷的半徑設置顯然既影響耗時也影響結果。截斷半徑越小，要計算的ERI越少，耗時就越低，但結果越偏離泛函原本的精確結果。顯然不能為了一味地節約時間而把截斷半徑設得太小，否則結果會很垃圾。6埃通常是耗時和精度的較好權衡，設得明顯更小會令精度有明顯損失，設得更大并不會令精度有明顯提升卻會增加耗時不少。如果對精度要求很高也可以用8埃，設得更大就完全沒必要了。 要注意像6、8埃這種程度的截斷半徑其實還沒有大到令絕對能量充分收斂，但能量在相互求差時，由截斷半徑帶來的系統性誤差能很大程度抵消。可以認為用不同的截斷半徑相當于用不同的理論方法，直接對比能量的情況必須保持截斷半徑的設置相統一！ 值得一提的是，不一定截斷半徑小結果就一定差，有可能截斷半徑造成的對ERI忽略的誤差和泛函本身的誤差相抵消導致恰好對某些體系某些問題的計算比用更大截斷半徑時結果還更好。 截斷半徑應當小于晶胞最短尺寸的一半以免導致得到缺乏物理意義的結果。這里所謂的最短尺寸，是(001)晶面間距、(010)晶面間距、(100)晶面間距三者中的那個最小值。也因此為了能夠使用6埃截斷半徑，晶胞最短尺寸應當大于12埃，如果晶胞不夠大則應當擴胞（即便拋開庫侖截斷問題不談，單從只能計算gamma點這點來說，晶胞也不能過小）。 Multiwfn創建的雜化泛函輸入文件中，若周期性設成了NONE，則不使用庫侖截斷，若是周期性，Multiwfn會自動加上庫侖截斷的設置。如果晶胞最短尺寸夠大，Multiwfn會自動把截斷半徑設為6埃，如果不夠大，則會根據當前晶胞最短尺寸自動設成能設的最大半徑。 以下是庫侖截斷的設置方式，加入到&DFT/&XC里面。POTENTIAL_TYPE TRUNCATED代表用庫侖截斷的1/r算符，CUTOFF_RADIUS設置截斷半徑（埃）。T_C_G_DATA t_c_g.dat是默認的可以不寫，這指定的是庫侖截斷用的gamma函數的數據文件t_c_g.dat，在CP2K的data目錄下。 &HF   &INTERACTION_POTENTIAL     POTENTIAL_TYPE TRUNCATED     CUTOFF_RADIUS 6.0     T_C_G_DATA t_c_g.dat   &END INTERACTION_POTENTIAL &END HF 3 積分屏蔽 積分屏蔽（integral screening）是指盡可能在不明顯犧牲精度的前提下減少要計算的ERI的量，這對于CP2K雜化泛函計算的耗時和內存使用量有至關重要的影響！上一節的庫侖截斷本身就是一種積分屏蔽策略，CP2K里還有另外兩種非常重要的積分屏蔽方法，也是普遍應用在Gaussian等量子化學程序中的： (1)利用Schwarz不等式做積分屏蔽：利用Schwarz不等式關系，可以根據數目較少的二指標ERI估計出數目甚巨的三、四指標ERI中哪些數值肯定小于閾值從而可以忽略掉，這可以大幅減少ERI要算的數目（原本總共要算的ERI中占比最大的就是四個基函數序號不同的那些ERI）。閾值通過&DFT/&XC/&HF/&SCREENING中EPS_SCHWARZ設置，此值越大，被忽略的ERI數目越多、計算精度越差。默認的1E-10雖然足夠保證精度，但此時太貴，沒必要。通常設為1E-6就可以了，是精度和耗時的較好權衡，想要精度更好點用1E-8。 (2)結合密度矩陣做積分屏蔽：ERI對能量的貢獻取決于它與特定密度矩陣元的乘積。因此哪怕ERI不是很小，但與之相乘的密度矩陣元很小，那么這樣的ERI也完全不需要計算。CP2K定義了特定方法快速估計一個ERI與相應密度矩陣元乘積值的上限，如果小于上述的EPS_SCHWARZ，則此ERI就會忽略掉。這種依賴于密度矩陣的積分屏蔽默認是關閉的，如果開啟則把&DFT/&XC/&HF/&SCREENING里的SCREEN_ON_INITIAL_P設為T。由于這種積分屏蔽能令耗時巨幅降低，因此我強烈建議總是開啟，也因此Multiwfn產生的雜化泛函的輸入文件里都自動把SCREEN_ON_INITIAL_P設成T。 基于密度矩陣做積分屏蔽有一個極其關鍵的要點是初猜的密度矩陣質量不能太差！然而很多人對此一無所知！默認情況下初始的密度矩陣是CP2K自動做初猜產生的，由于它和SCF收斂時的密度矩陣相差很大、質量很糙，因此此時積分屏蔽通常做得極爛，不僅可以忽略的ERI可能沒忽略，還會把一些重要的ERI給忽略掉，由此可能導致KS矩陣質量很爛、SCF極難收斂。而且忽略哪些ERI是在SCF過程一開始決定好的，即便之后隨著SCF迭代密度矩陣質量不斷變好，也不會由此使得積分屏蔽逐步變得合理，這導致即便最后SCF收斂了，得到的能量的精度也可能很爛。因此當SCREEN_ON_INITIAL_P設T時，必須讀取一個基本合理的波函數當初猜。通常是在雜化泛函計算前先用相同的基組做一個便宜的純泛函計算，比如PBE0計算前先用PBE做一個單點計算，得到記錄了PBE收斂的波函數信息的wfn文件，然后PBE0計算時用SCF_GUESS RESTART，并把WFN_RESTART_FILE_NAME指定成那個wfn文件，以讀取其作為初猜。這不僅使得基于密度矩陣的積分屏蔽可以合理進行，還有一個額外的好處是由于初猜波函數質量較好（至少比默認自動產生的好得多），使得雜化泛函計算時達到SCF收斂所需的迭代次數會比較少，這進一步節約了時間（注意即便拋開ERI計算耗時不談，雜化泛函在構造KS矩陣時花的時間也明顯高于純泛函，因為涉及到密度矩陣元與ERI的大量的乘加操作）。 一般的量子化學程序雖然也基于密度矩陣做積分屏蔽，但不是非得像CP2K這樣得讀取一個像樣的初猜波函數，這是因為它們會根據當前最新的密度矩陣判斷這一輪SCF用的ERI哪些是能被忽略的。CP2K的做法和它們不同，主要是因為CP2K傾向于讓用戶做incore形式的SCF而非一般量化程序一般用的direct形式的SCF，見后文。 當SCREEN_ON_INITIAL_P設T時，對于雜化泛函的幾何優化、分子動力學等任務的續算，要注意光有.restart文件還不行，還必須提供上一步的.wfn文件當初猜波函數，否則也是會由于基于自動初猜的密度矩陣做積分屏蔽非常惡心而令任務無法正常跑下去。 如果你嫌麻煩而不想先做純泛函計算產生收斂的波函數再做雜化泛函計算，那就不要寫SCREEN_ON_INITIAL_P T。對于小體系、小基組特別是非周期性的計算，不基于密度矩陣做積分屏蔽時耗時也不高。而對于計算量較大的雜化泛函計算，則強烈建議不要偷這個懶而多花巨量不必要的時間。 4 ADMM CP2K還支持其開發者獨創的ADMM (Auxiliary Density Matrix Methods)方法減少要算的ERI數目從而進一步顯著加快雜化泛函的計算，其原理見以下北京科音CP2K第一性原理計算培訓班的ppt。 ppt里提到的主基組就是指平時實際用的基組，諸如DZVP-MOLOPT-SR-GTH。對于不同的主基組，有一些專門構造的與之匹配的ADMM輔助基組，比如MOLOPT系列基組適合搭配的是CP2K的data目錄下的BASIS_ADMM_UZH文件里的那些輔助基組，從小到大依次是admm-dz、admm-dzp、admm-tzp、admm-tz2p。admm-dzp是2-zeta結合極化函數的輔助基組，通常算是保底質量。DZVP-MOLOPT-SR-GTH主基組適合搭配admm-dzp，更好的主基組如TZVP-MOLOPT-GTH結合admm-dzp也可以接受，若結合admm-tzp精度會更好些，這都比起不用ADMM能巨幅節約計算時間和內存。輔助基組的大小與主基組越接近，ADMM給雜化泛函HF交換項帶來的誤差越小，但需要計算越多ERI因而越耗時。如果輔助基組和主基組設為相同的，那么ADMM不會造成任何誤差，相應地也不會節約任何耗時（反倒還會增加一些耗時）。 data目錄下的BASIS_ADMM文件中的FIT3、pFIT3等輔助基組都比較過時了，對元素涵蓋也遠不如BASIS_ADMM_UZH里的完整，因此如今不再建議使用（它們出現在不少官方例子文件和文獻中，如今不要效仿）。目前版本的Multiwfn產生利用ADMM的雜化泛函的輸入文件時默認設的是admm-dzp。 注意前述這些輔助基組都是對GTH贗勢基組用的，千萬別做全電子計算時也用它們。要想全電子計算時用ADMM，可以用(aug-)pcseg系列基組作為主基組，在《在線基組和贗勢數據庫一覽》（http://www.shanxitv.org/309）中介紹的BSE基組數據庫中有和它匹配的(aug-)admm輔助基組，諸如pcseg-2適合搭配admm-2。pcseg是一類構造得很好的對DFT計算非常劃算的基組，《談談量子化學中基組的選擇》（http://www.shanxitv.org/336）里簡要提及了，在北京科音中級量子化學培訓班（http://www.keinsci.com/workshop/KBQC_content.html）里講基組的地方做了具體介紹。 用MOLOPT系列基組做周期性體系的雜化泛函計算幾乎是一定要結合ADMM的，否則極難算得動！因為MOLOPT是完全廣義收縮基組，而CP2K利用的計算源高斯函數（PGTF）之間ERI的libint庫又沒有專門考慮廣義收縮的情況，因此造成等效的PGTF數目甚巨，ERI的數目又與PGTF數目的四次方近似成正比，數目顯然多得恐怖。如果你用片段收縮基組，比如def2-SVP、6-31G**、pcseg-1、pob-DZVP-rev2之類去做雜化泛函計算，由于需要考慮的PGTF數目遠比檔次與之相仿佛的DZVP-MOLOPT-SR-GTH要少，因此ADMM方法不是非用不可，但不用的話耗時會顯著高于DZVP-MOLOPT-SR-GTH搭配admm-dzp的情況。 不能有的原子用ADMM加速計算而有的原子不用。因此如果你用混合基組，其中有的原子用的主基組有ADMM輔助基組，有的沒有，你又想用ADMM，那么沒有ADMM輔助基組的那些原子的輔助基組就只能設定成其主基組。雖然這些原子沒法被ADMM技術加速，但至少其它那部分原子能被加速。 前面提到的各種節約耗時的策略是彼此完全兼容的，同時利用可以最大程度節約時間。即周期性體系雜化泛函計算時一般要用庫侖截斷、利用Schwarz不等式做積分屏蔽、結合密度矩陣做積分屏蔽、開啟ADMM。所有這些設定都會對能量產生影響，因此寫文章的時候最好注明，橫向對比時應統一。 啟用ADMM需要在&KIND里加上BASIS_SET AUX_FIT [輔助基組名]，同時增加一行BASIS_SET_FILE_NAME [輔助基組文件名]。并且在&DFT里加入以下內容 &AUXILIARY_DENSITY_MATRIX_METHOD &END AUXILIARY_DENSITY_MATRIX_METHOD 如果用的是對角化而非OT，還需要在以上字段里加入ADMM_PURIFICATION_METHOD NONE要求不做純化。 5 ERI的儲存和內存使用量的控制 SCF有兩種常見形式，incore SCF是指SCF計算一開始將所有要用到的ERI全都算出來并儲存在內存里，之后SCF每一輪在構造KS矩陣時會直接從內存里讀取ERI而不再重新計算。direct SCF是指SCF每一輪都重新計算要用到的ERI，這也叫on-the-fly方式計算ERI。顯然，incore比direct總耗時低得多，但代價是機子的必須內存非常大，因為要存的ERI通常非常多。在CP2K的雜化泛函計算中，如果你給的內存足夠存得下所有ERI，那么只有SCF第一輪的時候由于要計算所有ERI所以耗時很高，而之后SCF每一輪的耗時就很低了。如果給的內存不夠儲存所有ERI，那么能在內存里存多少ERI就存多少，存的那部分在之后每一輪SCF中會直接讀取，存不下的ERI在每一輪SCF中會重新算，相當于介于incore和direct之間的情況。顯然，做大體系雜化泛函計算時應當盡量在大內存機子上做、分配給CP2K盡可能多的內存，以讓盡可能多的ERI能存到內存中、盡可能減少每輪SCF過程中要on-the-fly計算的ERI的量來節約耗時。 常有人問為什么他做雜化泛函計算時老是卡在SCF剛開始的位置不動，這是因為SCF第一輪肯定是要計算所有需要考慮的ERI，耗時顯然往往很高，對大體系、大基組的情況花一、兩個小時都是常事，看上去好像一直停在這里似的。如果你給的內存較少而導致大部分ERI都得on-the-fly方式去算，那么之后SCF每一輪都得花將近這么多時間，整個任務跑完可能得半天甚至一天時間。如果PRINT_LEVEL設的是MEDIUM，在SCF第一輪算完后就會輸出Number of sph. ERI's calculated on the fly和Number of sph. ERI's stored in-core，前者是之后每一輪SCF都需要重新算的ERI數，后者是已存到內存中而之后不需要SCF每一輪重新算的ERI數。顯然前者相對于后者越小，之后每一輪SCF比第一輪時間少得越多。 CP2K用戶通常默認用popt版本，每個MPI進程在計算HF交換項時最大內存使用量由&DFT/&XC/&HF/&MEMORY里的MAX_MEMORY控制，單位為MB，扣除掉儲存零碎數據的占用量外，其余都用來儲存ERI。MAX_MEMORY與MPI并行進程數的乘積必須小于當前機子空余物理內存量，顯然應當盡量給大以盡可能減少on-the-fly方式算的ERI量，但也不要太頂著頭分配，因為進入HF交換項計算模塊之前CP2K還會占一定量的內存，故要有適當余量。在SCF第一輪不斷往內存里寫入ERI、內存占用量不斷增大的過程中，若空余內存已用完時CP2K還在繼續往里存ERI，CP2K就會馬上崩潰，甚至還可能導致計算機暫時失去響應。所以MAX_MEMORY應當設得恰到好處。 如果空余物理內存很有限，為了能讓MAX_MEMORY能設得比較大，可以減小MPI進程數。為了不因此浪費CPU運算能力，此時可以用psmp版，此時每個MPI進程下屬的OpenMP線程都可以共享這個MPI進程儲存的ERI。 6 雜化泛函的定義方式 這一節簡要說一下雜化泛函怎么正確在輸入文件里定義。對于非周期性體系，比如用PBE0，直接在&DFT里寫上以下內容即可。此時不用庫侖截斷而用完整的1/r算符，EPS_SCHWARZ為默認的1E-10、SCREEN_ON_INITIAL_P為默認的F。 &XC   &XC_FUNCTIONAL PBE0   &END XC_FUNCTIONAL &END XC 顯然對周期性體系絕對不能簡單寫成上面那樣，至少也得自己去指定為使用庫侖截斷。結合前面的講解，周期性體系PBE0計算的泛函定義部分一般寫為下面這樣。PBE0泛函由75%的PBE交換項+25%的HF交換項+100%的PBE相關項構成，所以需要恰當地定義&PBE里面的參數，并在&HF里用FRACTION指定HF成份。 &XC   &XC_FUNCTIONAL     &PBE       SCALE_X 0.75       SCALE_C 1.0     &END PBE   &END XC_FUNCTIONAL   &HF     FRACTION 0.25     &SCREENING       EPS_SCHWARZ 1E-6       SCREEN_ON_INITIAL_P T     &END SCREENING     &INTERACTION_POTENTIAL       POTENTIAL_TYPE TRUNCATED       CUTOFF_RADIUS 6.0     &END INTERACTION_POTENTIAL     &MEMORY       MAX_MEMORY 3000     &END MEMORY   &END HF &END XC 有的泛函在&XC_FUNCTIONAL里有特定的字段，可以省得自己去定義細節參數，但HF成份依然必須自己定義，否則相當于只做了殘缺的純泛函計算，結果毫無意義。比如HF成分為10%的TPSSh泛函可以寫為 &XC   &XC_FUNCTIONAL        &HYB_MGGA_XC_TPSSH     &END HYB_MGGA_XC_TPSSH   &END XC_FUNCTIONAL   &HF     FRACTION 0.1     ...同上   &END HF &END XC 更多的泛函定義方式可以直接看Multiwfn產生的相應雜化泛函的輸入文件里的寫法，都是嚴格正確的。 7 一個常見警告：The Kohn Sham matrix is not 100% occupied 做雜化泛函計算，特別是用了ADMM時，很容易在SCF開始后看到The Kohn Sham matrix is not 100% occupied警告。具體原因在https://www.cp2k.org/faq:hfx_eps_warning有介紹，這里就不多說了。對于實際計算來說，看到這個警告時，先把&QS里的EPS_PGF_ORB設為比默認明顯更小的1E-12，如果計算正常，即便還有這個警告也可直接無視。EPS_PGF_ORB設得越小，越無法利用重疊矩陣和KS矩陣的稀疏性節約時間，但這對耗時的影響相對于雜化泛函計算的總耗時來說微不足道。目前版本的Multiwfn產生的雜化泛函的輸入文件里自動就把EPS_PGF_ORB設成了1E-12。 對于晶胞較小的情況，本身矩陣稀疏性就差，開發者建議上來就在&QS里設MIN_PAIR_LIST_RADIUS -1從而完全不利用矩陣的稀疏性節約時間，此時絕對不會再出現那個警告，但比起用EPS_PGF_ORB 1E-12時耗時會高一些。 8 簡單例子：C60分子晶體 最后，這里以C60分子晶體的單點任務作為例子演示一個標準的周期性雜化泛函的計算。如果你還不了解Multiwfn創建CP2K輸入文件的功能，建議先閱讀《CP2K第一性原理程序在CentOS中的簡易安裝方法》（http://www.shanxitv.org/586）。 此體系的cif文件和下面涉及的輸入文件都可以在http://www.shanxitv.org/attach/690/file.zip中找到。此體系結構如下所示，一共240個C原子。晶胞邊長是14.07埃，對于分子晶體來說不擴胞的情況下只考慮gamma點可以接受。 啟動Multiwfn，輸入C60.cif的路徑載入之，然后輸入 cp2k PBE.inp   //要產生的PBE計算的輸入文件名 0   //基于默認設置產生輸入文件 現在當前目錄下就有了PBE.inp，對應PBE/DZVP-MOLOPT-SR-GTH單點任務。然后接著在Multiwfn里輸入 cp2k PBE0.inp   //要產生的PBE0計算的輸入文件名 1   //選擇理論方法 -6   //PBE0結合ADMM 0   //產生輸入文件 由PBE0.inp內容可見，Multiwfn自動指定的計算設置完全合適，對應的是PBE0/DZVP-MOLOPT-SR-GTH結合admm-dzp輔助基組、6埃庫侖截斷、EPS_SCHWARZ 1E-6、SCREEN_ON_INITIAL_P T、開啟OT。 手動把PBE0.inp中的WFN_RESTART_FILE_NAME前頭的#去掉并設為PBE-RESTART.wfn，把SCF_GUESS RESTART開頭的#去掉，從而使得PBE0計算時讀取PBE收斂的波函數當初猜。然后根據你的機子的空余物理內存量和要用的MPI進程數設置合適的MAX_MEMORY。 把PBE.inp和PBE0.inp放到同一個目錄下，先運行PBE.inp，算完后當前目錄下就出現了PBE-RESTART.wfn，然后再運行PBE0.inp。在《淘寶上購買的雙路EPYC 7R32 96核服務器的使用感受和雜談》（http://www.shanxitv.org/653）介紹的筆者的2*7R32 512GB雙路服務器上用96核并行、MAX_MEMORY設5000，用CP2K 2023.2 popt版，PBE計算花了8秒，PBE0計算花了178秒。PBE0任務的SCF迭代過程如下所示，可見第一輪耗時較高，由于當前給的內存足夠儲存所有ERI，因而ERI沒有on-the-fly計算的，所以之后每一輪SCF的耗時遠遠低于第一輪。 9 總結 本文把強大的CP2K高效地做雜化泛函計算的各方面要點進行了介紹，以避免初學者對各種關鍵信息毫不知情而憑感覺胡用瞎用。本文的內容也充分體現出用Multiwfn產生CP2K輸入文件非常容易，令研究者從復雜的輸入文件編寫中充分解放，但這絕不意味著CP2K的使用者就可以對關鍵性技術細節毫不知情，否則會各種踩坑、白浪費時間、算出無意義的結果。CP2K用戶必須掌握的知識、要領遠多于大部分其它計算程序。筆者開設的北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）可以使得學員完整系統透徹掌握全部這些知識。 http://www.shanxitv.org/678 2023-07-28T19:53:00+08:00 談談重復不出來計算化學文獻里的數據的可能原因 On the possible reasons why data in computational chemistry literature cannot be reproduced 文/Sobereva@北京科音   2023-Aug-1 0 前言 經常有做計算化學的人在網上問為什么他沒法重復出來文獻里的計算數據。導致重復不出來的原因實在太多了，然而初學者提問此類問題時幾乎總是含糊其辭，提供的有效信息甚少，導致別人根本無從回復，或者需要別人羅列一大堆可能導致差異的原因。我感到每次回復此類問題太費事，于是在這里專門寫個文章說說這事。遇到這種問題的讀者仔細看完此文后，應該能料想出可能是什么原因導致自己和文獻的結果存在差異。就算還搞不明白，至少也應該能明白在網上問別人的時候應當提供什么信息，知道怎么把導致差異的可能原因的范圍盡可能縮小，免得別人打很多字羅列諸多當前不需要考慮的可能性，或者別人嫌你提供的有效信息完全不夠而根本不想回復。 這里首先要強調一點的是，如果你的目的就是為了精確重復文獻數據，那你就要盡最大可能使用和文獻完全相同的方式去做計算，而不要管人家的做法合不合理，哪怕你明知道人家的計算方式不當。下文說的內容都是為了讓你能盡可能準確重復出文獻里的數據，而不是講怎么合理地計算。如果你的實際目的是想獲得相同體系相同研究問題的合理的結果，那你只要確保你自己的計算是合理的就行了，沒必要非得把別人的數據重復出來，說不定人家的計算本來就有問題。如果最終驗證發現別人的計算有硬傷并因此導致結論錯誤，或者憑基本知識就能斷定別人的計算有問題，那么可以發一篇comment文章指正以免以訛傳訛，比如《我對一篇存在大量錯誤的J.Mol.Model.期刊上的18碳環研究文章的comment》（http://www.shanxitv.org/584）里介紹的文章。 1 結構不一樣 自己計算用的結構和文獻不一致，是導致數據無法重現的最常見原因之一。如果你用的結構就是文獻正文或補充材料里給的坐標，那就可以直接排除這個因素了。注意有些文獻的SI里鍵長、笛卡爾坐標單位是Bohr而不是埃，千萬別弄錯了。 如果你計算的基本結構特征都和文獻里不符，也即算的都不是同一個東西，能重現出文獻里的數據就怪了。很多體系有不同的形態，比如酚酞在酸性、中性、堿性環境下主要的存在形態明顯不一樣，氨基酸在真空和水環境下的質子化態明顯不一樣，姜黃素在水和有機溶劑中一種是烯醇式一種是酮式結構，等等。你必須確保你算的和文獻里算的確實是同一種狀態的結構，如果文獻里給了截圖，一定要仔細對照。 做計算要有基本的化學常識，并且要認真仔細，免得搭出來沒意義的結構去做計算，自然和文獻里對不上。比如文獻里給出Lewis結構式時通常是隱去碳上的氫的，如果計算的時候真不畫氫，或者少畫了個別的氫，自然計算毫無意義。再比如，文獻里算一個羧基配位的過渡金屬配合物，如果你在計算時羧基上還留著氫，明顯也是錯的。再比如，文獻里畫了一個帶一個小配體的過渡金屬配合物的Lewis式，若你算的是水溶液中的狀態并忽略了配位水分子，也自然是錯的。再比如，文獻里給出一個聯苯或者碗烯的二維結構式，看上去是在一個平面上，若你去計算其基態時真的擺成平面來優化，得到的也是大概率是錯誤的結構（有破壞平面的虛頻的結構）。還要注意一些結構的構建不要想當然，比如SiO2表面在水中時表面的O主要是以羥基及其去質子化形態存在的，比例和pH有關，如果都當成是烏禿禿的氧就錯了。結構構建的關鍵性信息在文獻里一般都會有具體說明，但也可能一些文獻中在描述上有疏漏。 有些體系在計算的時候結構構建方式不唯一，比如《要善用簇模型，不要盲目用ONIOM算蛋白質-小分子相互作用問題》（http://www.shanxitv.org/597）、《18碳環（cyclo[18]carbon）與石墨烯的相互作用：基于簇模型的研究一例》（http://www.shanxitv.org/615）和《使用量子化學程序基于簇模型計算金屬表面吸附問題》（http://www.shanxitv.org/540）里提到的簇模型，怎么構建簇明顯會影響結果，特別是當簇用得較小的時候。若要重復文獻里的數據，就要用盡可能相同的模型，哪怕其模型很昂貴。 很多分子有不止一種構象（勢能面極小點），柔性大分子甚至有一萬種以上的構象。幾何結構優化一般會優化到與初始結構最接近的勢能面極小點。而不同構象下計算的屬性可能有不小差異，因此你必須確保你用的構象和文獻里一致，也即結構優化使用的初猜結構要盡可能像文獻里給出的。如果文獻沒明確給出三維結構信息，按理說文獻應當用的是能量最低構象（確切來說是自由能最低構象）做的計算，但也有可能作者忽略了構象搜索并得到了不合理的結果。如果你對構象搜索一點概念都沒有，建議看下文了解些常識 《使用molclus程序做團簇構型搜索和分子構象搜索》（http://bbs.keinsci.com/thread-577-1-1.html） 《gentor：掃描方式做分子構象搜索的便捷工具》（http://bbs.keinsci.com/thread-2388-1-1.html） 《將Confab或Frog2與Molclus聯用對有機體系做構象搜索》（http://bbs.keinsci.com/thread-20063-1-1.html） 《使用Molclus結合xtb做的動力學模擬對瑞德西韋(Remdesivir)做構象搜索》（http://bbs.keinsci.com/thread-16255-1-1.html） 也不是所有情況都適合用能量最低構象來計算，比如某化學反應可能是從某個能量相對較高的構象發生的，可以對反應的IRC的反應物一端做幾何優化得到之。 團簇體系的構型搜索和柔性分子的構象搜索同等重要。哪怕水二聚體這么簡單的兩個小分子構成的團簇都都有不同的極小點構型，并且對應的結合能明顯不同。因此計算這種體系也必須確保和文獻用的構型相一致。如果文中沒明確交代，一般默認其用的是能量最低構型結構做的計算，但也可能作者忽略了這個問題，導致誤用了能量不是最低的構型并因此得到了錯誤的結果。如果你對構型搜索沒有概念，建議看此文：《genmer：生成團簇初始構型結合molclus做團簇結構搜索的超便捷工具》（http://bbs.keinsci.com/thread-2369-1-1.html）。 有的時候兩種極小點構型或構象可能結構相差很小。這種情況需要把可視化程序里的視角和文獻里的圖片擺得盡可能一致再去觀察，免得算的其實不是同一個結構。 上面說的是計算極小點結構的情況，自不必說，重復文獻里過渡態的計算結果時也必須嚴格確保和文獻里用的結構一致。 對于研究固體表面的情況，要注意文獻里的結構是怎么來的。表面是從體相直接切出來的未經優化的（unrelaxed），還是經過優化的（relaxed），還是使用的考慮表面重構的（reconstructed），這需要區分清楚。文獻中計算用的是具體哪個晶面也必須要搞清楚。固體表面吸附問題也有類似構象/構型搜索的勢能面存在多種極小點問題，一個小分子可能在表面上不同位點以不同方式吸附，實際算的結構必須和文獻相一致。文獻里算的吸附是物理吸附還是化學吸附也得搞清楚。對于吸附過程有勢壘的情況，初猜結構中你把被吸附物和表面擺得相距很近，一般優化完了的是形成新鍵的化學吸附，如果擺得較遠，一般得到的是弱相互作用維持的物理吸附，二者的結構、相互作用能相差甚大。研究表面吸附還有覆蓋率的問題，要注意文獻中對所用結構的相關描述，研究低覆蓋率和飽和覆蓋率時用的模型明顯不同。 用slab模型研究固體表面還要注意厚度問題。slab太薄的話，得到的功函數、吸附能等性質都不準確。對于二維層狀材料的研究，還要注意層數對計算的性質的影響。這些方面都應當和文獻嚴格一致。對表面類體系優化時還經常牽扯到對一部分原子凍結使之固定為體相的結構，之后振動分析也相應地牽扯凍結，凍結方式必須和文獻一致。 做某些體系的分子動力學模擬也同樣要注意初始結構的差異，否則現象可能和文獻不符。比如有文獻研究表明水中的蛋白質由于疏水效應會自發進入口徑足夠大的碳納米管，如果你一開始把蛋白質放得離碳納米管老遠，由于質量頗大的蛋白質本身整體運動就較慢，再加上碳納米管的隨意旋轉（沒凍結的話），可能文獻里說的那種蛋白質自發進入管內的現象跑挺長時間也出現不了。再比如想模擬結冰過程，哪怕一開始把溫度設得顯著低于冰點，但沒有在水中加入局部的冰的結構的話，跑很長時間也不會看到冰的結構的出現。至于在模擬的時間尺度內一定會出現的現象，初始結構就不重要了。比如水和乙醇的互溶，無論初始結構里把二者混合均勻，還是分成兩相，最后肯定都是互溶的。 做周期性體系的模擬還要注意盒子尺寸與文獻要一致。用第一性原理計算固體時，假設k點數是固定的，盒子尺寸的不同會影響精度。做動力學模擬時，盒子尺寸還間接影響原子運動行為以及模擬得到的屬性。 用Quantum ESPRESSO等平面波程序時，以及CP2K這種將平面波當輔助基函數的程序時，如果算的是某些維度為非周期性體系，還需要注意真空層的厚度，要盡可能和文獻一致。因為真空層可能對結果有不可忽視的影響，比如影響非周期性方向上與相鄰鏡像的作用是否可以忽略。而且比如CP2K里用MT的Psolver時還有非周期性方向上盒子尺寸大于相應方向有電子密度明顯分布的兩倍以上的要求，顯然真空區不能小了。 假設文獻用的是實驗測定的坐標直接進行性質的計算，還要注意實驗來源必須一致。比如文獻或晶體數據庫里對同一種晶體結構可能提供了不同溫度下測定的坐標，會有一定程度的不同（體現熱膨脹）。顯然壓力也直接影響測定的結構。哪怕是相同的測試環境，不同來源的晶體結構也會有差異，特別是氫的位置的確定，無疑這對氫鍵等問題的計算結果會有直接影響。做量子化學/第一性原理研究的人一般都會至少對氫的坐標重新進行優化來refine之，需要注意文獻里是怎么考慮的。蛋白質的高解析度和低解析度測出來的晶體結構里殘基側鏈的位置也可能相差甚大，甚至明顯影響到做酶催化、分子對接之類問題的結果。同樣的屬性，以不同類型實驗測定的結果注定也會有差異，比如氣相電子衍射和單晶衍射測的一個是氣相結構一個是晶體結構，對于柔性分子差異可能很大，這點參考《實驗測定分子結構的方法以及將實驗結構用于量子化學計算需要注意的問題》（http://www.shanxitv.org/569）。 上面列舉了諸多結構層面影響計算結果的可能性，當然不可能列舉得很完備，實際情況多種多樣。總之希望讀者能意識到嚴格保證和文獻里用的結構的一致對于重現數據有多么重要。 2 計算設置不一樣 為了精確重復文獻數據，要盡可能仔細閱讀文章及其補充材料，了解清楚文獻計算細節，使得自己的計算設置能最大程度與文獻相一致。本節里談到的很多可能導致差異的因素在《談談不同量子化學程序計算結果的差異問題》（http://www.shanxitv.org/573）里都有詳細說明，務必注意閱讀。 量子化學和第一性原理計算用的理論方法，以及理論方法中涉及的設定和參數（如范圍分離泛函的w參數、色散校正的經驗參數、CASSCF和多參考方法的活性空間的選擇、TDDFT計算是否用了TDA近似、DFT+U參數、CASPT2計算用的IPEA shift參數、DLPNO方法用的閾值、GW@DFT和DFT-SAPT具體基于什么泛函、sTDA計算用的經驗參數，等等）必須和文獻精確一致。有的理論方法在同一個程序里甚至有兩種定義，比如B3LYP在ORCA里寫成B3LYP還是B3LYP/G是不同的，要搞清楚文獻實際用的哪個（通常是前者）。還要注意理論方法名和關鍵詞的差異，比如里這里提及的：《Gaussian16里用PBE0關鍵詞計算的實際上是PBE0-DH雙雜化泛函》（http://bbs.keinsci.com/thread-13660-1-1.html）。有的文獻的寫法還可能存在含糊性，比如可能文獻就說對某個開殼層體系用了二階微擾理論，但實際上二階微擾理論用于開殼層的實現形式有一堆，比如UMP2、PUMP2、ROMP2、ZAPT2、ROMP、OPT2等等，如果文中給了引用的文獻，要根據文獻判斷具體是哪種。 對計算結果可能有明顯影響的各種條件必須一致，比如： ? 溶劑效應的考慮。如果用的是隱式溶劑模型，那么用的具體是什么方法？描述的是哪種溶劑？用CP2K的SCCS模型時α、β、γ參數是怎么設的？描述混合溶劑的時是怎么確定溶劑參數的？如果用的是顯式溶劑模型，那么溶劑分子排布是怎么考慮的？有沒有同時用隱式溶劑模型？ ? 計算熱力學校正量的模型。比如《使用Shermo結合量子化學程序方便地計算分子的各種熱力學數據》（http://www.shanxitv.org/552）里介紹的Shermo程序默認用的以及ORCA用的是Grimme的QRRHO模型，而Gaussian用的是較原始的不適合含低頻體系的RRHO模型 ? 計算熱力學量時是否考慮了平動和轉動，顯然算周期性體系時不應當考慮（在Shermo程序中這由settings.ini里的imode決定） ? 計算轉動對熱力學量貢獻時對轉動對稱數的考慮，看Shermo手冊附錄部分了解相關常識 ? 計算熱力學校正量、振動分析用的原子質量 ? k點數目和分布 ? 平面波相關計算的截斷能 ? CP2K中，做雜化泛函時用的庫侖截斷半徑、DFT+U計算用的布居計算方法以及考慮的角動量、Poisson solver的選擇 ? 對相對論效應的考慮 ? 對外場/外勢的考慮 ? DFT積分格點精度。參看《密度泛函計算中的格點積分方法》（http://www.shanxitv.org/69） ? 是否用了凍核、怎么凍的 ? 雙電子積分的積分屏蔽設置 ? 是否用了RI、DLPNO等加速計算方式 ? 是否考慮了偶極校正 ? QM/MM計算劃分的位置、MM與QM區之間的相互作用的考慮、link原子的設置 等等等等 基組必須和文獻嚴格一致，并注意上述http://www.shanxitv.org/573博文里說的許多細節問題。注意一個基組/贗勢名可能有不止一個指代，比如Stuttgart贗勢，有不同方式考慮相對論效應的版本、不同價電子數和氧化態的版本，Gaussian內置的和Stuttgart官網上的贗勢和贗勢基組定義往往又有所不同。再比如，GAMESS-US用戶寫CCD和CCT關鍵詞分別用的是cc-pV(D+d)Z和cc-pV(T+d)Z，而不是更常見的cc-pVDZ和cc-pVTZ。要弄清楚作者實際用的什么，仔細看文中的描述，并注意作者引的基組原文。 不光一般意義上的基組（叫主基組或者叫軌道基組）要和文獻計算時的一致，輔助基組也應當一致。輔助基組種類很多，比如RIJ用的/J輔助基組，RIJK用的/JK輔助基組，電子相關計算用的/C輔助基組，CP2K里加快雜化泛函計算用的ADMM輔助基組，CP2K里做LRIGPW計算用的輔助基組，等等。同一個主基組能夠或者適合搭配的某種目的的輔助基組可能不止一種，而且有的時候有含糊性和任意性，比如ORCA里用ma-def2-TZVP主基組時，可能有人用autoaux關鍵詞自動產生輔助基組，可能有的人用標準的def2/J輔助基組。 對于基于經典力場的模擬，力場參數顯然必須和文獻嚴格一致。如果參數是從文獻里拷來然后用在自己計算中的，注意單位轉換，以及1-4相互作用規則和LJ參數的混合規則等細節必須和文章相同。文獻給出的LJ參數有的是R0有的是σ需要區分。原子電荷也必須和文獻保持一致，如果原子電荷是作者自己算的，要搞清楚計算原子電荷的方法，并且如果是基于波函數算的，那么波函數是在什么幾何結構、溶劑環境、計算級別下得到的要弄清楚。這點很重要，比如對于RESP電荷，氣相下計算的和水環境下計算的，Hartree-Fock和B3LYP計算的，相差大了去了，參看這些討論：《RESP擬合靜電勢電荷的原理以及在Multiwfn中的計算》（http://www.shanxitv.org/441）、《RESP2原子電荷的思想以及在Multiwfn中的計算》（http://www.shanxitv.org/531）。 做動力學模擬的情況，必須確保模擬相關設置和文獻一致。比如步長、步數、坐標/速度/受力/能量保存頻率、參考溫度和壓力、控溫和控壓方法、熱浴和壓浴的時間常數、可壓縮系數、各向同性還是各向異性控壓、初速度怎么產生的、范德華和靜電作用的計算方式（如簡單截斷、twin-range、Ewald、PME、SPME等，以及其中涉及的截斷半徑，是否用了switching function等處理）、水模型是剛性還是柔性的、是否用了鍵/鍵角約束、用沒用能量/壓力色散校正，等等。 如果你用的計算程序和文獻里不同，或者和文獻里的相同但版本不同，那么由于默認設置不同，以及算法實現的不同，可能對結果帶來很大差異，一定要仔細閱讀《談談不同量子化學程序計算結果的差異問題》（http://www.shanxitv.org/573）。 在量子化學波函數分析領域也有很多細節設置會影響分析結果。這里以非常流行的波函數分析程序Multiwfn（http://www.shanxitv.org/multiwfn）為例： ? 按照《使用Multiwfn做電子密度、ELF、靜電勢、密度差等函數的盆分析》（http://www.shanxitv.org/179）、《使用Multiwfn的定量分子表面分析功能預測反應位點、分析分子間相互作用》（http://www.shanxitv.org/159）、《使用Multiwfn對靜電勢和范德華勢做拓撲分析精確得到極小點位置和數值》（http://www.shanxitv.org/645）等文章介紹的方法進行分析時，格點間距數值越小精度越高，但也越耗時 ? 按照《使用Multiwfn繪制態密度(DOS)圖考察電子結構》（http://www.shanxitv.org/482）繪制PDOS和OPDOS圖時，在Multiwfn里選擇的計算軌道成份的方法會直接影響結果 ? 按照《使用Multiwfn模擬掃描隧道顯微鏡(STM)圖像》（http://www.shanxitv.org/549）和《使用Multiwfn結合CP2K的波函數模擬周期性體系的隧道掃描顯微鏡（STM）圖像》（http://www.shanxitv.org/671）繪制STM圖時，偏壓的數值和常高模式掃描的平面的選擇都會明顯影響結果 ? 計算《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）和《使用AdNDP方法以及ELF/LOL、多中心鍵級研究多中心鍵》（http://www.shanxitv.org/138）里介紹的多中心鍵級時，基于原本的基函數和基于自然原子軌道（NAO）計算會得到不同的結果 ? 計算Hirshfeld原子電荷時需要用到原子自由狀態下的密度。用自行提供的原子波函數文件來計算此密度，還是直接用Multiwfn內置的此密度，會有所不同。 本節提到的這些計算細節是理應充分體現在文章的computational details里或者補充材料里的。但由于一些文獻的作者為了省篇幅，或者用的很多設置和程序默認的一致，或者由于作者計算水平比較外行、意識不到很多細節對計算非常重要，就在文中沒特意說。如果有些要點不知道作者具體怎么考慮的，自己把各種可能性都試試就知道了，文章沒特意說的設置就姑且用程序默認的。 3 計算的電子結構狀態不一樣 凈電荷和自旋多重度的設置必須和文獻里的相同，否則相當于算的是其它電子態，自然算出來的各種性質都不同。比如計算基態的氧氣發生的反應、計算富勒烯的原子化能，如果氧氣和碳原子誤當成了單重態來算，顯然結果是錯的。此外，稍微了解點兩態反應性的概念就會知道，不少體系在不同自旋態下反應機理都明顯不一樣，如果自旋多重度或凈電荷設錯了，自然反應物/產物/中間體/過渡態等無法得到和文獻里一致的結果，而且相差懸殊。 如果文獻沒注明、從文獻中看不出相關信息，那么凈電荷一般為0。如果被研究的體系的自旋多重度不那么顯而易見，或者研究中涉及到了不止一個自旋態，文中一般都會說明自旋多重度。如果確實沒說，那就得查閱文獻或者數據庫，或者自己判斷、測試了。 很容易被忽略的是開殼層單重態或者叫自旋極化單重態。用KS-DFT等方法來計算的話應當做對稱破缺計算，參看《談談片段組合波函數與自旋極化單重態》（http://www.shanxitv.org/82）。比如整體為單重態的反鐵磁性耦合體系、雙自由基體系、共價鍵被拉長的均裂狀態就是最典型的例子。還有一些特殊的體系也是如此，比如18碳環的D18h點群的過渡態結構、乙烯二面角扭轉成90度的結構、六并苯。重復文獻里的這些體系的計算時別誤當成了閉殼層單重態來計算（文中如果也是做了對稱破缺計算，一般都會明確說。如果沒說，沒準文獻也誤當成了閉殼層計算）。 對于第一性原理計算磁性材料時必須謹慎，很多情況下不是光設個整體的自旋多重度就完事的，而必須恰當設置初始磁矩才能收斂到真正要考察的磁性狀態。為了重復文獻（或材料數據庫）中的數據，若它們給了收斂的波函數對應的各個原子的自旋磁矩/自旋布居信息的話，最好根據其來設置初猜波函數里的原子自旋磁矩，以盡最大可能收斂到和文獻相同的波函數狀態。 要注意波函數穩定性問題。即便凈電荷、自旋多重度的設置是正確的，也不代表SCF收斂后得到的波函數就是實際想研究的態、和文獻里算的狀態相一致。如果文獻里研究的是基態，你算的結果和文獻對不上，而且體系又不是簡單有機體系而是電子結構特殊的體系，那么應當做一下波函數穩定性檢測，比如Gaussian里可以用stable關鍵詞，發現不穩定的話可以用stable=opt嘗試搞出穩定波函數。不穩定的波函數可能意味著這是某個激發態，還可能意味著此波函數是沒有任何物理意義的虛假的狀態。 4 計算的東西或方式不一樣 重復不出文獻的數據時還要想清楚到底你算的東西、計算的方式是否和文獻一致、有可比性。下面列舉一些典型的需要注意的情況。 比如通過能量差E(AB)-E(A)-E(B)考察倆分子A和B間相互作用強度，復合物結構肯定是要優化的，而單分子是直接從復合物優化后的結構里直接取出來，還是也做優化，這在不同文中的做法就不一樣了。單分子不優化情況下得到的能量差一般稱為（復合物結構下的）相互作用能，單分子也優化情況下得到的能量差一般稱為結合能，二者之間相差各個單體的變形能之和。對于諸如《透徹認識氫鍵本質、簡單可靠地估計氫鍵強度：一篇2019年JCC上的重要研究文章介紹》（http://www.shanxitv.org/513）里面列舉的那些剛性小分子，相互作用能和結合能差異甚微，但如果考察的是兩個容易變形的分子間的作用強度，比如兩個氨基酸分子的結合，那么相互作用能和結合能就必須明確區分了。J. Chem. Phys., 158, 244106 (2023)還專門討論了什么情況變形能會較大。對于鍵能，也是類似地有不同計算方式。所以研究此類問題的時候一定要搞清楚文獻具體是怎么算的。 再比如做SAPT能量分解計算，會得到很多耦合項，比如交換與色散作用的耦合項，這是歸到色散項里，還是歸到交換互斥項里，不同文章的做法可能不一樣，需要注意文章怎么說的。再比如計算反應的自由能問題，是計算標準反應自由能，還是計算特定濃度下的反應自由能，這是不一樣，相差濃度校正項，要搞清楚文獻實際算的是什么。再比如電子親和能、電離能、激發能，都有垂直和絕熱之分，文獻里計算哪種都有，也必須確保和文獻算的東西對應。 文獻里有時候對熱力學量的標注比較含糊，要分清楚文獻算的是什么溫度下什么量。比如E，文獻里可能指電子能量，也可能指內能，需要結合語境判斷。 模擬ECD譜需要轉子強度，有長度形式和速度形式兩種，諸如Gaussian都會輸出，在Multiwfn里按照《使用Multiwfn繪制紅外、拉曼、UV-Vis、ECD、VCD和ROA光譜圖》（http://www.shanxitv.org/224）里說的繪制ECD時可以選采用哪種，一般推薦用速度形式。如果你用的和文獻不一致，會導致峰型存在一定差異。 基組的CBS外推的方法不唯一，有不同公式，有的情況外推方式里也涉及到參數，參看《談談能量的基組外推》（http://www.shanxitv.org/172）。不光是能量外推，還有比如GW計算的準粒子能量隨基組的外推、幾何結構隨基組的外推等。一些文獻對高級別理論方法結合大基組還用一些近似方式估計，比如一種常見的是CCSD(T)/大基組 = CCSD(T)/小基組 + (MP2/大基組-MP2/小基組)，這在《各種后HF方法精度簡單橫測》（http://www.shanxitv.org/378）里專門說了。涉及到這些高精度數據估算時，需注意保證和文獻的做法一致。 還有一些量本來就沒有唯一定義，去重復文獻的數據更是得弄清楚文獻里用什么定義、什么方法算的。比如自由體積，我在《使用Multiwfn計算晶體結構中自由區域的體積、圖形化展現自由區域》（http://www.shanxitv.org/617）里介紹了一種合理的計算方式，但顯然這不是唯一方式。再比如分子半徑，我在《談談分子半徑的計算和分子形狀的描述》（http://www.shanxitv.org/190）就已經列舉了多種計算方式。還有諸如原子電荷，計算方法更是多了去了，如ADCH、RESP、Hirshfeld、NPA等，簡要介紹見《原子電荷計算方法的對比》（http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml）。 為了正確理解文獻具體算的是什么，必須把背景知識搞清楚，典型的就是gap這個概念，有HOMO-LUMO gap、band gap、optical gap、fundamental gap等，參見《正確地認識分子的能隙(gap)、HOMO和LUMO》（http://www.shanxitv.org/543）。如果沒搞清楚概念和定義方式，可能重復不出文獻的數據，或者實際上文獻的計算方式或用詞不當，你卻意識不到這一點。 某些問題的計算中可能涉及到一些實驗值，要搞清楚作者用的什么值。比如用量子化學程序計算含碳物質的生成焓，通過構造熱力學循環可知其中需要利用實驗的（或第一性原理計算的）石墨的升華焓，要搞清楚這個量是從哪來的，具體是多少。再比如， 量子化學計算某分子的標準氧化還原勢，其中利用到標準氫電極的電勢，其數值有不同來源，從4.05到4.44 V都有報道，Phys. Chem. Chem. Phys., 16, 15068 (2014)建議對理論計算的情況用4.28 V。 還有一些研究涉及到一些額外的參數需要注意，典型的就是頻率校正因子，見《談談諧振頻率校正因子》（http://www.shanxitv.org/221）。是否使用頻率校正因子，以及對某一特定級別使用什么目的的、哪個來源的校正因子，都會影響結果。 周期性體系的計算中，軌道能量的絕對數值是沒意義的，不同程序用的能級零點也不一樣，需要關心的只有相對數值，比如帶隙。所以不要問為什么DOS圖的形狀和文獻一致但橫坐標和文獻不一致。這也正是為什么文獻在繪制DOS圖時往往一般把價帶頂或者費米能級位置設為0點，因為這樣才便于公平對比。 還有些量在文獻里用的習俗和你用的可能不同。比如算結合能問題，有少數文章把結合過程的能量降低量作為正值記錄。再比如NBO給出的E2超共軛作用能，程序輸出時把正值作為超共軛作用造成的體系能量降低，但有的文獻報道E2值的時候取的是負值。 如以上列舉的情況可見，算某個東西的時候必須仔細閱讀文獻搞清楚作者是怎么算的、用的什么定義、算的到底是什么。 5 數據的可重現性本身就差 有些數據很容易重現，比如在RHF/6-31G**級別下優化水分子的基態結構，不管用什么主流量化程序、怎么優化，得到的都是相同結果，初始結構很隨意。有些數據本身就是難以精確重現，典型的就是凝聚相體系的分子動力學。由于分子動力學過程有混沌效應，而且實際計算中往往不可避免地引入隨機因素，會導致文獻中的動力學模擬的定量結果幾乎無法嚴格精確重現。關于這點我專門寫過文章，參看《數值誤差對計算化學結果重現性的影響》（http://www.shanxitv.org/88）。更具體來說也看重現文獻里的什么動力學模擬數據。比如小分子的液態性質相對來說是容易重復的。比如重復密度值，文獻里是801.92 g/L，你算出來801.88 g/L，這就可以認為是很好重復出來了。如果你算出來是810.21 g/L，而且模擬流程和設置是合理的，那就是有明顯不可忽視的差異了，已經超過了隨機性的影響了，需要根據前面所說的來考慮導致與文獻存在差異的可能原因。也有的量的重現相對比較困難，比如磷脂膜里插入了一個分子，去計算它的側向擴散系數，這就屬于相對比較難算準的，因為本身分子在磷脂膜這種擁擠的環境里擴散就困難，而且體系里就這么一個分子，沒法通過同類分子的平均來減小統計誤差。就算你自己跑的時間非常長以保證統計精度絕對足夠高，但文獻里給出的數據未必統計精度就夠高。像這種問題，對文獻里的數據的重現精度就別要求太高了。 還有一種情況是動力學模擬的現象在有必然性的同時本身也有一定隨機性，跑出什么現象有運氣成分。比如距離蛋白質的口袋一定距離處放一個小分子配體，在水環境下做動力學模擬，考察2 ns模擬時間內配體能否和蛋白質結合，初速度隨機生成。像這種模擬的結果就很有隨機性，如果跑10次模擬，可能有的模擬剛過幾百ps小分子就進入口袋了，也可能有的模擬開始后小分子恰好往蛋白質相反方向游離，跑到2 ns還沒進入口袋。對于這種問題的研究必須做大量的模擬來得到統計結果，較嚴謹的文獻普遍都會這么考察這種現象明顯有隨機性的問題，你若要重復文獻也必須以相同的方式去模擬。如果模擬的樣本數不夠多的話，也不要指望結果特別相符。比如前述問題，可能文獻里跑10次發現有5次進入口袋，你模擬10次發現有3次進入口袋，這未必是模擬設置和文獻有區別，而純粹是隨機性所致，要以正確的態度來看待這種情況。 順帶一提，量子化學、第一性原理做靜態的計算的可重復性比起分子動力學模擬要好得多得多，但也不排除極個別情況前者可能受到隨機性的明顯影響。比如對勢能面很復雜的大分子進行優化，特別是QM/MM計算時可能牽扯到幾千原子的情況的優化，有可能由于微小的隨機性因素（如使用了并行計算所帶來的）導致兩次優化分別收斂到了勢能面上可能相距挺遠的兩個極小點，且有肉眼可察覺的差異。 6 對相符程度要求過高 有人其實可能已經合理重復出了文獻中的數據，但由于他對重復的精度要求過高，導致誤以為沒重復出來。 頭腦里要有收斂限的概念。比如文獻里優化的結構二面角是123.64度，你優化出來的是123.69，這就可以認為是重現出來了，因為這點差異通常都已經小于計算程序的幾何優化的默認收斂限了，差異基本來自初始結構的任意性。再比如，文獻里DFT計算出的反應能是-8.12 kJ/mol，你算出來的是-8.18 kJ/mol，也完全可以接受了。如果本來你用的計算程序和文獻里都不同，在盡可能令設置和文獻相同的情況下，你算出來比如是-8.32 kJ/mol也可以算重復出了文獻，各種雞毛蒜皮不值得一提的因素足矣帶來這種程度的差異。但如果你算出來是比如-9.8 kJ/mol，那就明顯不是不值得一提的差異了。再比如TDDFT計算電子激發能、DFT計算HOMO能級，和文獻里相差0.02 eV，這都可以算一致，但如果相差到0.1 eV，這就不算一致了，畢竟TDDFT算激發能的精度都在0.3 eV左右（泛函選擇得當的前提下），當差異達到接近方法本身的誤差的數量級時明顯就算顯著差異了。 重復文獻數據的時候用的收斂限應當足夠嚴，起碼令收斂限顯著小于你期望的對文獻數據的重現精度。雖然文獻作者用的收斂限未必夠嚴，因此光是你把收斂限設嚴未必有意義，但至少先做自己這里能做的事。 7 自己的計算存在硬傷 重復不出文獻數據有可能是自己計算存在硬傷導致的。初學者缺乏常識、經驗和敏銳性，計算很容易存在硬傷，例如： ? 關鍵詞沒寫對或設置不當，導致數據沒意義，比如：Gaussian里想要用PBE0泛函但寫的關鍵詞是PBE0（此時做的是PBE0-DH泛函計算）、用贗勢基組時只定義了基組部分而沒定義贗勢、自定義基組時漏定義了某些原子或元素、Gaussian里用IOp自定義泛函時用opt freq關鍵詞（IOp沒法自動傳遞到freq這一步導致振動分析不是在與優化相同級別勢能面的極小點計算的）、Gaussian里用后HF方法算偶極矩時沒寫density結果得到的是HF級別的、計算UV-Vis光譜時算的態數太少導致光譜不夠完整、CP2K計算時用雜化泛函并且基于密度矩陣做積分屏蔽時沒讀取純泛函波函數當初猜（《解決CP2K中SCF不收斂的方法》http://www.shanxitv.org/665文中提到了）、CP2K算磁性體系時沒寫UKS結果當成了閉殼層算、Gaussian里按照特定方向加電場計算時不知道要寫恰當用nosymm導致電場加的方向和期望的不符（nosymm用處見《談談Gaussian中的對稱性與nosymm關鍵詞的使用》http://www.shanxitv.org/297）、用sobtop（http://www.shanxitv.org/soft/Sobtop）直接產生拓撲文件后沒有給里面添加原子電荷、用sobtop產生拓撲文件時誤把可旋轉二面角以諧振勢對待、用Multiwfn做空穴-電子分析時沒按照《使用Multiwfn做空穴-電子分析全面考察電子激發特征》（http://www.shanxitv.org/434）說的在Gaussian輸入文件里寫IOp(9/40=4)，等等 ? 缺乏基本知識瞎算，比如：計算小晶胞體系不知道考慮k點、算極化率時不知道要帶彌散函數、B3LYP算色散主導的弱相互作用卻不加色散校正、CP2K里用MOLOPT基組計算需要CUTOFF很高的Na的時候用的CUTOFF不夠、界面體系用了各向同性控壓而不知道應該用半各向同性控壓、算氣-液界面體系在垂直于界面方向開了控壓結果導致真空區消失、用Multiwfn基于CP2K產生的molden文件分析前沒在里面添加盒子信息（參見《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》http://www.shanxitv.org/651）、做Mulliken布居分析用的基組帶了不該帶的彌散函數、算溶劑環境下溶質自由能的時候忘了考慮標準態濃度轉換問題（參看《談談隱式溶劑模型下溶解自由能和體系自由能的計算》http://www.shanxitv.org/327）、試圖用后HF方法計算分子軌道（卻渾然不知程序雖然給出了結果但實際上是HF部分產生的）、用背景電荷表現環境分子的勢場用的是對靜電勢重現性極爛的QEq電荷（相關討論看《基于背景電荷計算分子在晶體環境中的吸收光譜http://www.shanxitv.org/579），等等 ? 計算流程/方式不對，比如：算高精度能量和性質之前結構沒經過幾何優化、直接把IRC兩端的結構當做反應物和產物的準確結構用于計算能量和屬性數據、振動分析和IRC計算之前結構沒有在嚴格相同級別下優化過、對動力學軌跡計算某個分子的RMSD衡量內部結構變化的時候不知道先要消除整體平動和轉動、算電子發射光譜之前沒優化激發態、模擬拉曼光譜時沒按照《使用Multiwfn計算特定方向的UV-Vis吸收光譜》（http://www.shanxitv.org/648）里說的將拉曼活性轉化成拉曼強度 ? 數據讀取錯誤。比如Gaussian里做了雙雜化泛函計算，結果讀成了其中SCF Done后面的中間量（怎么正確讀看《談談該從Gaussian輸出文件中的什么地方讀電子能量》http://www.shanxitv.org/488）。再比如做幾何優化，Gaussian會對初始和最終結構都做布居分析，本來想讀最終結構下的，結果誤讀成第一次輸出的。再比如做開殼層體系的NBO分析時，NBO會對總密度矩陣、alpha密度矩陣和beta密度矩陣的分析結果依次輸出，如果讀錯地方就糟糕了 ? 對計算數據質量缺乏把控。如果數據質量太糟糕，自然可能無法重現文獻數據。各個部分的計算都要保證數據質量。比如算穩定狀態必須確保波函數穩定、確保無虛頻。雖然做波函數穩定性測試和振動分析會增加耗時，但必要時必須做，否則可能得到錯誤的結果而無法重復文獻。不要輕易使用比如Gaussian里opt=loose這種關鍵詞放寬收斂限，否則很可能收斂限太松導致某些情況對文獻數據重復性太差，尤其是此時做振動分析得到的頻率的準確度很可能很糟糕。用取消SCF收斂情況檢查的IOp(5/13=1)這種危險關鍵詞的時候必須謹慎檢查最終SCF收斂情況，千萬別最后收斂精度很爛卻直接使用其數據。再比如對于通過分子動力學計算物質平衡狀態屬性時，必須保證采樣充分使統計誤差足夠小。再比如考察某種磷脂膜體系的特征時，需要模擬時間足夠長以使得體系充分達到平衡，本身這類體系達到平衡所需的時間就偏長，動輒需要100 ns以上 ? 沒恰當考慮對稱性問題。對于有對稱性的體系，在Gaussian等支持對稱性的量子化學程序里，要盡可能利用對稱性，這樣計算又好又快，內行的研究者的文章里對于能利用對稱性的情況幾乎都會利用。不利用對稱性甚至可能還會得到和文獻不符的結論，比如文獻報道某個體系是平面的，但你一開始擺的結構是彎曲的，由于優化的收斂不可能無窮精確，最后你得到的結構可能輕微偏離平面，與文獻結論不一致。如果體系本身沒有那么高對稱性，你一開始當成高對稱性結構來計算，最后可能錯誤地得到了高對稱性結構，也和文獻不符（這種情況只要做一下振動分析看有沒有虛頻便知）。 ? 其它低級錯誤：單位轉換因子沒用對或者用的不準確（比如百毒搜出來的大量轉換因子都是錯的）、該做構型/構象搜索的情況沒做或者做得不當、該考慮構型/構象權重平均的情況沒考慮、優化表面吸附問題時把被吸附分子放得太遠導致還沒優化到吸附狀態就滿足了收斂判據、處理數據用的Excel表格或者腳本存在問題、科學計數法記錄的數據沒讀對（如漏掉了指數部分），等等 8 文獻自身有問題 千萬不要只懷疑自己而不懷疑文獻，盡信文獻不如無文獻，很多文獻里（包括IF很高的期刊的文章里）的數據是有錯誤的，甚至是非常明顯的錯誤，上一節提到的各種犯錯的例子也都出現在很多文獻里。比如《我對一篇存在大量錯誤的J.Mol.Model.期刊上的18碳環研究文章的comment》（http://www.shanxitv.org/584）里我就指出一篇文章里從頭到尾的一大堆錯誤，諸如作者計算C18的時候大概率由于用的是非平面的初始結構，導致優化出的C18不是精確平面的，然后他們發現C18的極小點結構竟然有ECD信號，明顯這是虛假的結果。 還有些文獻里對數據的解釋、標注是錯的，是他們對數據理解有問題。比如直接把HOMO、LUMO能量的負值說成是氧化、還原勢是很多涉及電化學的文章里常見的事，顯然你如果以正確的方式來算，也即計算溶液下的電極反應的自由能變，結果肯定和文獻對不上。再比如有不少文獻用的費米能級是錯的，對有gap的體系直接就把價帶頂當做費米能級，這明顯不符合正確的確定有限溫度下費米能級的方式（正確定義參看wiki相應條目和Multiwfn手冊3.300.9節）。 還有些文獻里的計算方式從描述上就知道是錯的。比如有的文獻寫G=Eele+ZPE-TS，有本科化學知識的人都知道這明擺著不對，H(T)和U(0)之間的差異上哪去了？就算是作為近似不考慮這部分、假設這部分在求差的時候能極大抵消，至少也應該把等號寫成約等號。這種明顯不合理的數據就別去重復了，非要重復的話那就暫時性降智故意用錯誤的公式。 很多文獻作者也對數據很不負責任。比如之前我在IF很高的文獻里看到補充材料中的Gaussian輸入文件里居然全帶著IOp(5/13=1)，這文獻的數據的可靠性沒法令人不放心（雖然作者可能檢查了最終收斂情況，但誰也不敢說作者確實這么做了）。現在有很多無良代算機構，沒一丁點搞理論計算的人的基本操守，為了能給出令做實驗的人滿意的與實驗相符的計算數據來交差，極有可能在算不出期望的數據的情況下捏造假數據，這樣的數據更別指望能被別人重復出來了。甚至有的代算方連計算都不計算，直接憑空胡編亂造數據，建議讀者看看《談談我對計算化學代算的看法》（http://www.shanxitv.org/505）里提及的情況。 有很多文獻對計算的關鍵要點、直接影響數據可重復性的因素交代得非常粗糙、簡略、語焉不詳，甚至根本不提，這無疑給讀者帶來了很多誤解、給別人重復計算造成了障礙。重復這些文章的數據只能是考慮各種可能情況去試了。另外，還有可能文章里用的參數、設置本來就不慎寫錯了，作者寫的程序本身存在bug等等。 對于已經恰當考慮了前面說的所有因素，竭盡全力去重復文獻，但還是重復不出來（且確實是有不可忽視的差異而不是在復現精度上過度較真），如果沒有絕對必要去重復文獻，那就別去重復了，極有可能文獻的數據就是有問題的，或者沒交代關鍵性計算細節。這種情況只要你能保證自己的計算是合理的就夠了。但如果有特殊原因就是非要重復文獻不可，顯然該做的不是在思想家公社QQ群或者計算化學公社論壇（http://bbs.keinsci.com）等公開場所提問，因為找誰都不如直接找作者，世界上還有誰比作者更清楚數據是怎么來的的人么？世界上還有誰比作者更有回復你的義務？發郵件給通訊作者，向作者索要相應數據的輸入輸出文件，或者向作者提供你復現文獻計算的輸入輸出文件問作者為什么沒重復出來即可。當然很有可能一直沒收到作者的回復，可能性有這些： (1)作者看漏了你的郵件，或者郵件被誤歸入垃圾郵件，或者正忙著什么事沒來及回復你結果后來又忘了。可以過一個禮拜（且可以換個郵箱）再發一次試試。如果通訊作者有多個，也可以給另外的人發 (2)措辭不當，不夠禮貌和誠懇，而作者脾氣又大或者又比較忙，又覺得回復你也沒什么好處，就不回了。可以過一陣重新措辭再發一次。為了讓作者有動力花時間給你回復，最好表示你會在未來的文章中會引用他的文獻。 (3)文獻中的數據是作者找別人代算的，代算方又是非常不負責任的，相關文件也不提供、計算細節也不交代清楚，甚至代算之后就再也聯系不上了，因而作者也無能為力。 (4)做計算的是學生，學生畢業后找不到人了，導師又不清楚計算細節、沒原始文件。 (5)文獻中的計算本身就有嚴重錯誤或是假數據，可能收到你的郵件后作者才意識到，又或者可能之前就已經明知道用的是假數據，總之都不希望你知道他的數據存在問題，免得面子上掛不住，甚至怕文章最后被撤稿。 久久精品国产99久久香蕉