使用Multiwfn模擬掃描隧道顯微鏡(STM)圖像
注：筆者后來另寫了《使用Multiwfn結合CP2K的波函數模擬周期性體系的隧道掃描顯微鏡（STM）圖像》（http://www.shanxitv.org/671），專門介紹怎么用Multiwfn結合CP2K對固體表面體系繪制STM，是對本文的重要補充，請讀者別忘了閱讀。
使用Multiwfn模擬掃描隧道顯微鏡(STM)圖像
文/Sobereva@北京科音
First release: 2020-Apr-28 Last update: 2020-May-1
0 前言
掃描隧道顯微鏡(scanning tunneling microscopy, STM)是非常常用的在原子級別對固體表面或者被吸附在固體表面的分子進行成像的方法。從2020-Apr-27更新的Multiwfn開始，Multiwfn支持了模擬分子體系STM圖像的功能，易學易用、設置靈活、計算速度快、作圖效果好。本文將對其原理和使用方法進行介紹，并給出具體例子。
如果對Multiwfn不熟悉，建議參看《Multiwfn FAQ》（http://www.shanxitv.org/452）和《Multiwfn入門tips》（http://www.shanxitv.org/167），此程序可以在http://www.shanxitv.org/multiwfn免費下載。STM圖和態密度(DOS)有密切關系，Multiwfn另有強大的功能專門用來繪制DOS，見《使用Multiwfn繪制態密度(DOS)圖考察電子結構》（http://www.shanxitv.org/482）。
1 STM的原理和模擬方法
STM在wiki上有清楚的介紹：https://en.wikipedia.org/wiki/Scanning_tunneling_microscope，也有相關專著比如Introduction to Scanning Tunneling Microscopy (2ed, Julian Chen, 2008)。本節只是介紹一下最關鍵的STM的相關知識，以讓讀者能正確理解和使用Multiwfn的STM模擬功能。在Multiwfn手冊3.300.4節還有更細致的介紹。
STM的示意圖如下
如上所示，被測樣品平行于XY平面上。STM針尖(tip)在被測物原子上方進行掃描，范圍如上圖粉框所示。在樣品與針尖足夠近，且在二者之間施加偏壓（bias voltage，以下簡寫為V）時，由于隧道效應，樣品原子與針尖之間會形成隧道電流（tunneling current，以下簡寫為I）。針尖在不同(x,y,z)的位置時I都不同，STM本質上就是根據I(x,y,z)函數來進行成像。當V<0，電子從樣品向針尖流動，當V>0，電子從針尖向樣品流動。
STM有兩種模式：
(1)常高模式（constant height mode）：針尖保持在特定的z位置，對x和y做二維掃描，這樣得到的圖像就體現出不同(x,y)處I的差異。對于比較簡單的情況，針尖與樣品原子距離越近I就越大，因此通過這樣的圖就可以判斷樣品表面的原子排布情況。
(2)常電流模式（constant current mode）：針尖也是對x和y做掃描，但對每個位置，都調節針尖高度找到I等于特定值的位置（稱為z'）。常電流模式的STM圖就相當于z'(x,y)的平面圖。對于簡單情況，原子出現的地方z'會較大（針尖在較高處時就已經產生了特定大小的電流），反之z'較小（針尖需要放得較低才能產生特定大小的電流）。
對于被固體吸附的分子體系，影響STM圖特征的不僅是原子位置，還有其電子結構和偏壓。
模擬STM的方法很多，比如JACS, 121, 5392 (1999)、JPCB, 101, 5996 (1997)、JCP, 104, 2410 (1996)等，但是真正流行的只有Tersoff-Hamann (TF)模型。此方法在Phys. Rev. Lett., 50, 1998 (1983)提出，后來在Phys. Rev. B, 31, 805 (1985)進行了擴展討論。TF模型是基于Bardeen計算隧道電流公式通過推導和簡化得到的近似的模型。原本的TF模型局限性很大，對應于低溫和小偏壓極限，而且只適用于周期性體系。TF模型當初考慮了兩種情況，一種是將針尖近似為無限小的點，還有一種是將針尖視為是特定半徑的球形，一般用的TF模型都是前者，相當于把針尖的真實特征徹底忽略掉了（本來真實特征也難以確定），此時得到的STM圖對應于無窮高解析度的情況。TF模型的表達式傳遞出的一個關鍵性信息是某個位置的I正比于此處的局部態密度（local density-of-states, LDOS），故實際在用TF模型模擬STM時其實就是等價于計算和繪制LDOS，而對于I的絕對大小并不關心（這本來也沒法準確計算得到）。
對于分子、團簇等孤立體系，TF模型寫為下式（雖然這并非是TF模型原文里提出的，但寫文章時往往都這么叫）
其中r是針尖所在的位置矢量，E_F是費米能級，φ是軌道波函數，e是基元電荷（elementary charge），V是前述的偏壓。i循環能量在E_F+eV到E_F范圍的占據分子軌道，a循環能量在E_F到E_F+eV范圍的非占據分子軌道。費米能級對于孤立體系而言沒有確切定義，但在習俗上往往采用HOMO和LUMO能量的平均值，下文也這么考慮。通過上式可見，計算用來模擬STM圖的I實質上就是計算特定能量區間內的前線軌道的概率密度的總和，這被研究STM的文章叫做LDOS（注意此語境下的LDOS和《使用Multiwfn繪制態密度(DOS)圖考察電子結構》里說的LDOS不完全一樣），因此使用TF模型近似模擬STM時說到單位的時候應當寫成LDOS的單位，在Multiwfn里用的是a.u.。
2 Multiwfn模擬STM圖的實例：18碳環
關于Multiwfn模擬STM圖像功能的細節看Multiwfn手冊4.300.4節的介紹，這里直接通過18碳環這個體系作為例子講解怎么模擬，使用的是2020-Apr-27更新的Multiwfn。
18碳環這個體系在《一篇最全面、系統的研究新穎獨特的18碳環的理論文章》（http://www.shanxitv.org/524）介紹的筆者論文中做了非常系統、全面的研究和分析，歡迎閱讀。本例得到的圖像已被筆者用于在Carbon期刊上發表的論文An sp-hybridized all-carboatomic ring, cyclo[18]carbon: Bonding character, electron delocalization, and aromaticity（https://doi.org/10.1016/j.carbon.2020.04.099）的補充材料里，非常推薦一讀。
2.1 準備輸入文件
用Multiwfn模擬STM圖像用的輸入文件必須包含軌道波函數信息，可以用比如wfn、wfx、fch、molden、mwfn等等，詳見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）。如果你是Gaussian用戶，建議用fch文件，因為它不僅包含占據軌道也包含空軌道信息，因此偏壓<0和>0的情況都可以繪制，而如果用比如wfn這樣只包含占據軌道信息的輸入文件，則只能繪制偏壓<0的情況。18碳環體系用Gaussian在wB97XD/def2-TZVP下優化并產生的fch文件可以在這里下載：http://www.shanxitv.org/multiwfn/extrafiles/C18.zip，計算時用的關鍵詞是# wB97XD/def2TZVP opt。
模擬STM圖像用的輸入文件當中分子應當基本平行于XY平面以比較符合分子被固體表面吸附時的狀態，這并不要求體系非得是平面的。從上面鏈接給出的C18.fchk文件中可見這個環狀分子已經精確處在XY平面上，因此可以直接用于模擬。具體來說，此分子是處在Z=0的XY平面，對于純平面體系建議都擺成這樣以便于分析討論（Multiwfn對原子的Z坐標并沒有嚴格要求，所有原子的Z坐標比如都為2.5埃也可以，只要恰當解釋結果即可）。
有的時候Gaussian在計算純平面體系時會把體系擺到比如YZ平面去，此時沒法直接作圖。解決辦法之一是修改Gaussian輸入文件里的坐標，讓體系處于Z=0的XY平面，然后加上nosymm關鍵詞避免被擺到標準朝向下，詳見《談談Gaussian中的對稱性與nosymm關鍵詞的使用》（http://www.shanxitv.org/297）。但更好的解決辦法是利用Multiwfn的旋轉結構和波函數的功能處理一下，這不花任何額外的時間，詳見Multiwfn手冊3.300.4節模擬菲的STM圖的例子。
STM繪制功能只支持能給出分子軌道的理論方法，一般就是用KS-DFT，顯然后HF方法、CASSCF等都是不支持的，雙雜化泛函也不行。限制性閉殼層、限制性開殼層、非限制性開殼層波函數都支持。
2.2 模擬常高模式的STM圖像
啟動Multiwfn，載入C18.fchk，然后輸入
300 // 其它功能(Part 3)
4 // 模擬STM圖的功能
此時從屏幕上的提示會看到，當前費米能級被設為了HOMO和LUMO能量的平均值，并且默認的偏壓被設為了費米能級與HOMO能量的差值。這兩個值具體是多少可以從當前STM界面的選項上的提示看到，分別是-3.374 eV和-5.075 V，后者為負值說明在這樣的設定下電子是通過隧道效應從C18往探針轉移的。當前考慮的占據軌道的能量的下限為E_F + eV = -5.075 eV + e*(-3.374 V) = -8.449 eV，這正對應于18碳環的HOMO能級。也就是說，在默認設置下直接繪制STM圖，則這個圖將直接對應HOMO的概率密度乘上其簡并度。
18碳環這個體系的最高一批占據軌道信息如下所示，這在Multiwfn的主功能0里通過Orbital info. - Show occupied orbitals功能可以直接輸出（詳見《使用Multiwfn觀看分子軌道》http://www.shanxitv.org/269）
Orb: 48 Ene(au/eV): -0.405113 -11.0237 Occ: 2.000000 Type:A+B
Orb: 49 Ene(au/eV): -0.402139 -10.9428 Occ: 2.000000 Type:A+B
Orb: 50 Ene(au/eV): -0.402139 -10.9428 Occ: 2.000000 Type:A+B
Orb: 51 Ene(au/eV): -0.314238 -8.5509 Occ: 2.000000 Type:A+B
Orb: 52 Ene(au/eV): -0.314238 -8.5509 Occ: 2.000000 Type:A+B
Orb: 53 Ene(au/eV): -0.310478 -8.4485 Occ: 2.000000 Type:A+B
Orb: 54 Ene(au/eV): -0.310478 -8.4485 Occ: 2.000000 Type:A+B
可見，此體系最高四個占據的分子軌道能量非常接近，其中53、54是簡并的in-plane pi軌道，51、52是簡并的out-plane pi軌道。我們當前模擬STM圖希望將這四個軌道的貢獻都納入，為此我們需要恰當設置偏壓。比如偏壓設為-3.5 V時，考慮的占據軌道能量的下限就是-5.075+(-3.5)=-8.575 eV，它大于MO50的能量但是小于MO51的能量，故此時只有序號>=51的那四個占據軌道會對STM圖產生貢獻。
默認情況下，STM圖是常高模式，并且繪制的XY平面的Z值比Z值最大的原子的Z值高0.7埃，因此當前體系默認情況下繪制的常高STM圖對應于Z=0.7埃的位置，本例就直接用默認值了。通過屏幕上的選項還可以修改格點數（越大圖越精細，但耗時越高）、修改X和Y的計算范圍，本例都用默認值。默認的圖像的X、Y范圍是對邊界原子延展3 Bohr來確定的，這通常是適合的。
我們接著輸入以下命令
2 // 修改偏壓
-3.5 // 偏壓值(V)
0 // 開始計算
由于Multiwfn代碼效率高，而且常高STM圖像的模擬只需要計算一個平面上的數據，所以瞬間就算完了，在屏幕上可以看到以下信息
Lower limit of MO energy considered in the calculation: -8.575 eV
Upper limit of MO energy considered in the calculation: -5.075 eV
The MOs taken into account in the current STM simulation:
MO 51 Occ= 2.000 Energy= -8.5509 eV Type: Alpha&Beta
MO 52 Occ= 2.000 Energy= -8.5509 eV Type: Alpha&Beta
MO 53 Occ= 2.000 Energy= -8.4485 eV Type: Alpha&Beta
MO 54 Occ= 2.000 Energy= -8.4485 eV Type: Alpha&Beta
Totally 4 MOs are taken into account
Grid spacings in X and Y are 0.100278 0.099212 Bohr
Calculating, please wait...
Maximal value (LDOS) is 0.017536 a.u.
以上信息顯示了考慮的軌道能量的上下限，并提示了哪些軌道被納入了考慮。可見我們期望的那四個軌道都已納入考慮了，說明偏壓設置合理。計算完畢后還提示了平面數據中最大值（即最大的LDOS）是0.0175 a.u.。
此時進入了一個新菜單，用來將計算完的LDOS平面數據繪制成STM圖像，里面有大量選項用來調整作圖效果，還可以把平面數據導出成文本文件以便通過第三方工具如Sigmaplot重新作圖。如果你經常用Multiwfn的主功能4（平面圖繪制功能），對這些選項肯定都比較熟悉。由于這些選項一看名字就知道是干什么的，這里我就不一一解釋了，不清楚的自己玩玩便知。我們就直接選擇選項0來繪圖，馬上看到下圖
上圖的色彩刻度條對應的是LDOS值，實際中的STM圖像展現的各個位置的隧道電流大小也正比于它。默認設置下色彩刻度下限為0，上限為此平面中的LDOS最大值。上圖體現出在碳原子和部分C-C鍵上方0.7埃的位置有相對明顯的隧道電流，因此顏色很白，而在另外一些C-C鍵上方0.7埃處隧道電流則可以忽略不計，因此是黑色的。如果你仔細讀過前述我的18碳環的研究文章就會知道此體系里有兩類C-C鍵，一類較短，一類較長，上圖體現出前者上方的隧道電流比后者的更強，這是因為較短C-C鍵的pi軌道作用相對顯著得多。
我們可以根據實際需要調節作圖設置以得到更滿意的圖，比如我們輸入以下命令做一些修改：
5 // 顯示化學鍵
14 // 用棕色
2 // 修改圖像類型
2 // 填色圖+等值線
7 // 修改坐標軸的刻度間隔
1.5,1.5,0.002
-3 //其它作圖設置
2 // 修改小數點位數
1 // X軸小數點位數
1 // Y軸小數點位數
3 // Z軸小數點位數
8 // 修改原子標簽大小
70
0 // 返回上一級菜單
0 // 重新繪圖
此時看到的圖像如下，效果比之前的更好了。
2.3 模擬常電流模式的STM圖像
下面我們對18碳環模擬常電流模式的STM圖像。進入STM圖像繪制功能后，先選擇選項1將繪圖模式切換為常電流模式，并且同前例一樣把偏壓設為-3.5 V。常電流模式要對一個三維區域中均勻分布的各個點計算隧道電流，相當于要計算LDOS的三維格點數據。默認的X、Y、Z范圍以及各個方向的格點數從屏幕上的選項的文字中都可以看到。格點數越多，計算耗時越高，而圖像越平滑，默認用的150*150*80的格點算是圖像效果和耗時的一個權衡。默認的X、Y范圍和常高模式一樣，而Z范圍默認是Z最大的原子上方的0.7~2.5埃區域，通常默認設置就是比較合適的。
現在直接選擇0開始計算隧道電流（即計算LDOS），耗時明顯高于只考慮一個平面的常高模式，不過在一般4核機子上對這個體系也就不到一分鐘就算完了。之后屏幕上提示LDOS格點數據最大值是0.0175 a.u.。接下來會看到幾個選項，選1可以進入圖形界面觀看隧道電流的等值面，LDOS=0.003 a.u.的等值面如下所示。筆者也選了show data range復選框，使實際計算的區域用藍色矩形展現了出來
從上圖可見，計算的X、Y、Z范圍的設置是比較恰當的。等值面在C11-C12這樣較短的C-C鍵上方的位置相對較高（即維持特定的隧道電流值時針尖位置較高），而在諸如C12-C13這樣較長的C-C上方比較低，這充分體現出這兩種C-C鍵特征的差異，也體現出分子的電子結構對常電流模式STM圖的影響。然后我們點Return按鈕關閉圖形窗口。
當前我們還可以用選項2將LDOS格點數據導出成為STM.cub，這樣就可以用VMD繪制成等值面了，效果會更好。利用《在VMD里將cube文件瞬間繪制成效果極佳的等值面圖的方法》（http://www.shanxitv.org/483）文中的做法可以輕松得到下面的圖像
接下來，我們繪制常電流模式STM的平面圖。選擇3 Calculate and visualize constant current STM image，輸入常電流的值。由于格點數據最大值是0.0175，顯然輸入的數值必須大于0且小于0.0175。這里我們就輸入0.003。然后程序就基于LDOS三維格點數據計算出平面數據，并輸出最小值和最大值：
Minimal Z is 0.700000 Angstrom
Maximal Z is 1.288940 Angstrom
即曰，在對應LDOS=0.003 a.u.的常電流STM平面圖中，Z最小值為0.7埃，這也正是格點數據的Z范圍的最小值；最大值為1.2889埃，這沒達到格點數據Z范圍最大值（2.5埃），說明在當前電流值下我們計算的格點數據的空間范圍是合適的，即Z的上限足夠大。
然后我們進入到了和常高模式時一樣的繪圖界面，有豐富的選項控制作圖。在選擇選項0進行繪制之前我們像上例一樣對作圖設置進行修改，最后作圖效果如下
在這張圖中，色彩刻度對應的是探針的Z位置（埃），此圖傳遞出的信息和前面的等值面圖一樣，也是在維持電流不變的情況下，探針在較短的C-C鍵上方的位置比較長的C-C鍵明顯高得多。這張圖和前面2.2節繪制的常高模式的STM圖的基本特征也非常相似。
如果你還想考察其它常電流下的STM圖，就退出作圖界面，然后再次選擇3 Calculate and visualize constant current STM image并輸入其它常電流的值即可。
值得一提的是，在首次于凝聚相中觀測到18碳環的文章Science, 365, 1299 (2019)的補充材料的中，作者給出了吸附在NaCl上的18碳環的隧道電流為0.7 pA、偏壓為0.1 V的實驗STM圖像，但那個圖模糊不清，而且圖像受基底NaCl影響比較明顯。
4 總結&其它
本文介紹了STM實驗的基本概念和STM圖像的模擬方法，并且以目前正熱門的18碳環體系為例演示了怎么繪制此分子的STM圖像，可見Multiwfn繪制STM圖像又方便效果又好。大家在以量子化學方式研究新奇特的分子的時候，都可以用此文的方法作圖放在文章里，以供實驗化學家參考、與實驗圖像進行對照。
本文的做法是針對孤立體系。如果你是通過第一性原理程序研究固體表面的STM，則本文的方法不適用，不過很多第一性原理程序都直接提供了STM模擬功能，比如Abinit、WIEN2k、CASTEP、OpenMX、GPAW、exciting等。在量子化學范疇中據我所知目前直接支持STM模擬的程序只有xtb一個，但沒什么實際意義（參數沒法調節，沒法直接產生圖像，手冊里也沒有任何介紹，而且根據其源代碼我發現它用的方法不對），而Multiwfn基于xtb產生的molden文件可以如上文一樣直接作圖。