使用CP2K過程中常用的可視化工具
使用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）。