使用Multiwfn結合CP2K的波函數模擬周期性體系的隧道掃描顯微鏡（STM）圖像
使用Multiwfn結合CP2K的波函數模擬周期性體系的隧道掃描顯微鏡（STM）圖像
文/Sobereva@北京科音 2023-Jun-4
1 前言
強大的波函數分析程序Multiwfn（http://www.shanxitv.org/multiwfn）有非常方便好用的模擬隧道掃描顯微鏡（STM）圖像的功能，模擬的原理以及Multiwfn中此功能的用法我在《使用Multiwfn模擬掃描隧道顯微鏡(STM)圖像》（http://www.shanxitv.org/549）里有專門說明，但那篇博文是以孤立體系為例。Multiwfn支持對CP2K產生的周期性體系波函數進行大量的分析，在北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）里我做了很多介紹，其中包括模擬周期性體系的STM圖。本文的目的是舉一個完整的例子，能令甚至之前沒用過CP2K的人都能輕易使用Multiwfn模擬周期性體系的STM圖。沒讀過http://www.shanxitv.org/549的讀者必須要先仔細閱讀此文。不了解Multiwfn的話建議看《Multiwfn FAQ》（http://www.shanxitv.org/452）和《Multiwfn入門tips》（http://www.shanxitv.org/167）。值得一提的是，CP2K自己也有模擬STM的功能，但遠沒有Multiwfn來得好用和方便，所以這里就不多說了。
2D Mater., 6, 015005 (2019)文中給出了幾層厚度的黑磷薄片表面的0.7 V偏壓下的常電流模式的STM圖，如下所示，本文將對黑磷計算同樣條件下的STM圖并與之對照。
本文為了節約時間只用單層黑磷，而且圖省事直接從三維晶體中截出來一層結構直接做單點計算得到波函數。嚴格來說應該用多層模型，并且固定最下層黑磷并令上層的自發弛豫。
本文的例子涉及的文件都可以在http://www.shanxitv.org/attach/671/file.rar中得到。讀者務必使用2023-Jun-4及以后更新的Multiwfn版本，否則情況與本文所述不同。
2 得到單層黑磷的molden文件
如《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》（http://www.shanxitv.org/651）所述，用Multiwfn分析CP2K的波函數必須讓CP2K產生記錄波函數信息的molden文件。這一節我們就對單層黑磷做一個單點計算來得到此文件。
本文文件包里的Phosphorus-black.cif是黑磷的晶體結構，用GaussView打開，只保留中間的單層部分而刪除其它原子，如下所示，然后保存為Phosphorus-black.gjf。
Multiwfn有極其便利的創建CP2K輸入文件的功能，見《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587），這里我們用Multiwfn創建產生molden文件的單點任務文件。由于CP2K沒法產生考慮k點時的molden文件，因此必須構造超胞用gamma點計算。當前的晶胞參數是a=3.31埃、b=4.38埃，在a方向復制為7倍、b方向復制為5倍后兩個方向都有>20埃，此時只考慮gamma點是可以接受的。當前模擬STM用的偏壓為正，此時能量最低的一個或多個空軌道對STM會產生貢獻，因此必須要求CP2K求解出一些空軌道才行，幾十個就夠了。模擬STM用的偏壓越大應當計算的空軌道越多，因為越高的空軌道會被涉及。
啟動Multiwfn，載入Phosphorus-black.gjf，然后輸入
cp2k //產生CP2K輸入文件
SP.inp //將要產生的CP2K輸入文件名
-7 //設置周期性
XY //XY二維周期性
-2 //產生molden文件
-11 //進入結構編輯界面
19 //產生超胞
7 //a方向復制的倍數
5 //b方向復制的倍數
1 //c方向（垂直于表面的方向）不變
-10 //返回
-9 //其它設置
12 //計算空軌道
20 //算最低20個空軌道
0 //返回
0 //產生輸入文件
此時SP.inp就產生在當前目錄下了，在本文的文件包里提供了，計算級別是Multiwfn默認的PBE/DZVP-MOLOPT-SR-GTH。用CP2K運行SP.inp，在我的2*EPYC 7R32雙路服務器上不到10秒鐘就算完了。之后按照《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）所述，把晶胞信息和價電子信息手動加入.molden文件，即在開頭插入以下字段。本文文件包里的SP-MOS-1_0.molden是已經改好的。
[Cell]
23.17000000 0.00000000 0.00000000
0.00000000 21.90000000 0.00000000
0.00000000 0.00000000 12.11600000
[Nval]
P 5
3 用Multiwfn模擬黑磷常距離模式的STM
首先模擬一下常距離模式的STM。啟動Multiwfn，載入SP-MOS-1_0.molden，然后輸入
300 //其它功能 Part 3
4 //模擬STM
2 //設置偏壓
0.7 //0.7V，和文獻里的相同
0 //開始計算
默認計算的是Z最大值原子上方0.7埃的XY平面，X和Y方向計算范圍和當前晶胞的X和Y方向的跨度一致。瞬間就算完了，在新出現的作圖設置菜單里輸入0繪圖，圖像如下所示
可以再做一些調整，只讓最上層的P用粗體字標注，下層的用細體字顯示以作區分，并且用和文獻里相似的配色，并且調整坐標軸讓刻度整齊。為此，把圖像關閉后，接著輸入
4 //修改顯示原子標簽的距離閾值
0.71 A //在圖上顯示距離作圖平面0.71埃以內的原子的標簽
y //將距離作圖平面更遠的原子的標簽以細體字顯示
9 //修改色彩變化方式
13 //黑-橙-黃
7 //修改標簽間隔
3,3,0.0005
8 //修改色彩刻度
0,0.002 //下限和上限
1 //保存圖像
此時的圖像如下所示
此圖的效果已經很好了。和實驗STM圖像展現出的信息近似一致。
4 用Multiwfn模擬黑磷常電流模式的STM
這次模擬常電流模式的STM，這與文獻里的STM圖對應。啟動Multiwfn，載入SP-MOS-1_0.molden，然后輸入
300 //其它功能 Part 3
4 //模擬STM
2 //設置偏壓
0.7 //0.7V，和文獻里的相同
1 //切換模式為常電流模式
0 //開始計算
3 //繪制STM圖像
0.0009 //常電流值（反復嘗試，取一個圖像效果較好、和實驗圖像較接近的即可）
4 //修改顯示原子標簽的距離閾值
2.6 A //在圖上顯示距離被計算區域頂端2.6埃以內的原子的標簽（計算的區域頂端的z坐標是計算開始之前界面里7 Set range in Z direction選項后面顯示的第二個值。當前設2.6埃可以只把頂層磷原子納入進來）
y //將距離更遠的原子的標簽以細體字顯示
9 //修改色彩變化方式
13 //黑-橙-黃
7 //修改標簽間隔
3,3,0.1
1 //保存圖像
現在看到下圖，和常高模式的STM圖展現出的信息基本一致
5 其它
本文的例子充分體現出用Multiwfn+CP2K繪制周期性體系表面的STM真是又快又方便，鼓勵大家應用在實際當中。請別忘了恰當引用Multiwfn。作為練習，大家可以繪制一下單層MoS2和石墨炔的STM圖。
對于無gap的體系，如金屬表面（包括吸附小分子后），繪制STM有一點需要注意。這種情況一般都是開smearing的，此時費米能級附近的軌道占據數不為整數，導出的molden文件里的軌道占據數也是如此，此時Multiwfn沒法直接繪制STM（Multiwfn會以為是記錄自然軌道的多組態波函數）。解決方法是先進入主功能300的子功能9（計算費米能級功能），此功能會將電子按照軌道能量由低到高以整數方式填充來修改軌道占據數，相當于0 K的占據方式。此外，在此功能里輸入當前的溫度，程序會給你費米能級值，將之記下來。退出此功能后，就可以照常進入STM繪制界面，屆時再把費米能級值輸入進選項3 Set Fermi level即可。另外值得一提的是，用Multiwfn創建CP2K輸入文件時開啟smearing的話，產生的輸入文件里自動就會有計算空軌道的設置，就不用自己再單獨設置了。