談談自旋密度、自旋布居以及在Multiwfn中的繪制和計算
談談自旋密度、自旋布居以及在Multiwfn中的繪制和計算
文/Sobereva @北京科音
First release: 2016-Nov-9 Last update: 2021-Aug-30
1 自旋密度
1.1 基本概念
對于開殼層體系，Alpha電子密度分布和Beta電子密度分布是不同的，為了考察未配對電子（或者說單電子）在三維空間中的分布情況，定義了自旋密度(Spin density)的概念：顯然，對于三維空間中某個點，自旋密度若為正，說明此處Alpha電子比Beta電子多；為負則說明此處Beta電子比Alpha電子多。對于閉殼層體系，由于Alpha和Beta電子完全匹配，二者密度處處相同，所以沒有自旋密度。對自旋密度進行全空間積分，結果是Alpha與Beta電子數的差值。
初學者經常問一個含糊不清、令人難以回答的（或者說需要打好多字才能嚴謹地回復的）問題：“怎么計算自旋密度？”。自旋密度是個三維實空間函數，具體要計算哪個位置的自旋密度？到底打算怎么考察它？究竟想怎么表征其分布？必須說清楚！
有不同方式的來分析、表征自旋密度的分布。一般通過等值面圖來考察。下面我們來看一些開殼層體系的自旋密度等值面圖，看看它能說明什么。圖中綠色代表正值，藍色為負值。等值面數值調節為了能盡可能清楚體現自旋密度分布特征的情況。圖片皆為筆者開發的Multiwfn基于Gaussian09產生的波函數文件繪制的。
這是三重態卡賓(CH2)的自旋密度圖。可見沒成對的Alpha單電子主要繞著碳分布，并且集中分布在它沒成鍵的區域（對比甲烷結構可知本來能再形成兩個共價鍵的）。
這是NO的自旋密度圖，可見單電子是在垂直于分子軸的p軌道上。回憶N2有兩個pi成鍵軌道和一個sigma成鍵軌道，NO比N2多一個電子，故這個電子應該呆在pi反鍵軌道上（sigma反鍵軌道能量太高）。從圖上也確實看出了pi反鍵軌道的特征。
下圖是把乙烷拉開一定距離，使得C-C鍵充分斷開，然后通過非限制性開殼層計算得到的自旋密度。可見左邊的甲基帶著一坨Alpha單電子，右邊甲基帶著一坨Beta單電子。當兩個甲基距離足夠近的時候，這兩坨自旋彼此相反的單電子就會配對形成C-C鍵，變成閉殼層狀態。類似這樣把共價鍵拉斷后的電子結構特征和雙自由基體系是相同的，有關計算上的討論見《談談片段組合波函數與自旋極化單重態》（http://www.shanxitv.org/82）
把乙烯的兩個亞甲基扭成彼此成90度時的自旋密度如下。此構型下顯然乙烯的pi鍵已經破壞了，兩個碳各懸著一坨單電子在p軌道上。
這是雙核過渡金屬配合物，呈現反鐵磁性耦合特征。此例兩個過渡金屬各帶一部分單電子，且兩個過渡金屬帶的單電子自旋相反。
下圖是甲基自由基進攻乙烯得到丙烯自由基反應的過渡態結構下的自旋密度。可見在此結構下，單電子已經不完全在甲基上了，而是有一部分已經出現在了乙烯上。而且還可以看到乙烯上出現了一塊beta密度占主導的區域，可想而知隨著反應進行，這個碳的beta單電子就可以和甲基上alpha單電子配對形成新的共價鍵了。
這是HOO自由基的自旋密度，可以視為把雙氧水弄走一個氫。可見自旋密度最大的正是丟了氫的氧原子。這個氧單電子最多，最活潑，如果有另一個自由基和HOO碰上，肯定也是最先在這個地方成鍵。
這是Li@Calix[4]pyrrole的自旋密度。此體系是把帶著一個單電子的Li原子塞到閉殼層的Calix[4]pyrrole的籠中。這個體系中單電子是怎么分布的？憑化學直覺不好說，而從自旋密度圖上立刻就知道Li原本的單電子分布位置發生了很大變化，不是再繞著Li核了，而是鼓出去并彌散在了很大一片區域。
1.2 與自旋密度相關的問題
使用CASSCF、TDDFT等方法是沒法得到被計算的態的自旋密度的，比如哪怕你用TDDFT算三重態激發態也是得不到其自旋密度的，這和理論方法的形式有關系。用這些理論方法時要考察單電子分布只能用Multiwfn計算odd electron density，它和自旋密度一樣可以圖形化或定量展現單電子分布，而且對任意能產生密度矩陣的方法都適用，但它不像自旋密度那樣可以通過正、負體現單電子是alpha還是beta自旋。詳見《使用Multiwfn計算odd electron density考察激發態單電子分布》（http://www.shanxitv.org/583）。
還有一個與自旋密度密切相關的函數叫做自旋極化參數函數ζ，灑家最早看到此函數是在Density functional theory of atoms and molecules (Parr,Weitao Yang)一書中，其定義為
可見和自旋密度的差異是多了分母項，分母就是總電子密度。這個函數體現了不同位置處未成對電子占總電子的比例。在某個位置若此函數值為0，說明此處沒有單電子，若在某個位置此函數值達到上限1.0，說明這個地方的電子全是單電子。我們看看丁烷雙自由基的自旋密度和自旋極化參數函數的等值面圖：
可見這兩種函數都體現了Alpha和Beta單電子分別集中分布在丁烷雙自由基的兩端。但是自旋密度等值面明顯離核比較近，因為離核較遠處無論是alpha還是beta電子都很少了。
下面也順帶提一下自旋密度與概念密度泛函的一點聯系。福井函數是預測化學反應位點極為重要的實空間函數，簡要介紹和應用可參見筆者的《親電取代反應中活性位點預測方法的比較》（http://www.whxb.pku.edu.cn/CN/abstract/abstract28694.shtml）。福井函數用于預測親核反應位點時的具體形式是f+，定義為：
類似地，福井函數用于預測親電反應位點時的形式f- = ρ(N)-ρ(N-1)可以近似計算為N-1電子態的自旋密度，對應于Parr函數中的P-。
另外提一下，對于限制性開殼層(RO)形式的計算，明確區分了雙占據和單占據軌道，因此此時的自旋密度可以對所有單占據軌道波函數求模方再加和得到，因此容易討論軌道對自旋密度的貢獻。但是RO下的自旋密度明顯沒有非限制性開殼層(U)得到的準確，一些相關討論可參看Pople等人的Int. J. Quantum Chem., 56, 303。
利用原子核處的自旋密度可以計算原子對超精細耦合常數的Fermi contact項的貢獻，參見J. Phys. Chem. A, 101, 3174 (1997)。
1.3 繪制自旋密度
下面老夫介紹下如何用Multiwfn程序繪制自旋密度，本文用的是3.3.9版，以后版本可能選項會有所不同，以所用版本實際屏幕提示和手冊為準。對于最常用的Gaussian程序來說用.wfn/.wfx或.fch文件作為Multiwfn的輸入文件即可，產生這些文件的做法在Multiwfn手冊第四章開頭有明確說明。Multiwfn可以在http://www.shanxitv.org/multiwfn免費下載，使用入門參見《Multiwfn入門tips》（http://www.shanxitv.org/167）。下文用的例子是Multiwfn自帶的示例文件三重態甲酰胺的.wfn文件。1.3.1 等值面圖
這一節我們繪制三重態甲酰胺的自旋密度等值面圖。啟動Multiwfn，依次輸入examples\formamide-m3.wfn
5 //計算格點數據
5 //自旋密度
2 //用中等質量格點，總共算約51200個點。對于較大體系，建議選3用高質量格點（對于很大體系則建議選4并輸入較大格點數），否則等值面可能不平滑
此時Multiwfn開始計算自旋密度格點數據，算完后屏幕上看到的Summing up positive value and multiply differential element是全空間積分值，結果和2.0很接近，即對應三重態時Alpha比Beta多的兩個電子。
輸入-1可以進入圖形界面觀看自旋密度等值面，各個選項按鈕功能不言自明不再累述。選圖形界面菜單中的Set lighting里的Enable all，并將Isosurface style設為Use solid face+mesh后看到的圖形如下，可見單電子主要在氧上，其次是碳上。
如果想保存圖片，可以點RETURN關閉圖形窗口后選1。如果想把自旋密度格點數據導出成cube文件(.cub)以便在VMD等第三方可視化程序中顯示等值面，選2。
將Multiwfn與免費的VMD程序（http://www.ks.uiuc.edu/Research/vmd/）相結合可以繪制出更漂亮的自旋密度等值面、視角可以更自由調節，過程見筆者錄制的演示視頻“使用Multiwfn結合VMD繪制自旋密度等值圖”（https://www.bilibili.com/video/av26312131）。更好的、步驟明顯省事得多，而且還可以達到絕佳的繪制效果的做法是直接用筆者寫的VMD腳本，見《在VMD里將cube文件瞬間繪制成效果極佳的等值面圖的方法》http://www.shanxitv.org/483），其繪制過程之簡單、效果之好完爆任何一個其它的可視化程序！
如果想繪制前述的自旋極化參數函數的等值面，把Multiwfn目錄下的settings.ini里的ipolarpara設為1，保存文件并重啟Multiwfn。之后，實空間函數5就不再是自旋密度而是自旋極化參數函數了，因此再按照前述步驟繪制出的等值面圖就是自旋極化參數了。默認的等值面數值不適合考察這個函數，大家需要手動調節等值面數值直到便于觀看，比如調為0.7。繪制自旋極化函數的平面圖、曲線圖也是設好ipolarpara后再按照后文繪制自旋密度的步驟操作即可。
1.3.2 平面圖
這里我們繪制三重態甲酰胺的N3-C1-O6三個原子定義的平面上的自旋密度填色圖+等值線圖，這能十分充分地展現指定平面上的自旋密度分布細節信息。啟動Multiwfn并輸入examples\formamide-m3.wfn
4 //繪制平面圖
5 //自旋密度
1 //填色圖（想繪制什么類型的平面圖就選什么）
[按回車用默認格點數]
4 //用三個原子定義作圖平面
3,1,6
當前看到的圖像效果還不太理想。我們右擊圖像關閉之，選1，然后輸入色彩刻度下限和上限，即-0.1,0.06。我們輸入的色彩刻度上限比默認的上限數值要低一些，這會使得不同區域數值的大小能更好地通過色彩區分開。然后選2，讓等值線圖也同時顯示在填色圖上以便于定量考察。之后選-1重新顯示圖像，我們看到效果極好。圖中越紅的地方自旋密度越正，白色區域代表數值超過了當前色彩刻度上限的區域
想保存圖像就關了圖像后選0。
1.3.3 曲線圖
最后，我們將三重態甲酰胺的碳-氧鍵上的自旋密度繪制成為曲線圖。啟動Multiwfn輸入examples\formamide-m3.wfn
3 //繪制曲線圖
5 //自旋密度
1 //通過兩個原子核定義作直線來作圖
1,6
馬上看到曲線圖。圖中虛線是數值為0的位置，X軸上的紅點是原子核位置，左邊的是C1，右邊的是O6，可見自旋密度在原子核處還是很高的。
Multiwfn提供了很多選項對作圖效果進行各種調整、改進，在屏幕上提示得極為簡單易懂，手冊相應章節也有介紹，這里就不多說了。而且平面圖、曲線圖繪制后，在菜單中也都可以看到有選項用來導出原始數據，可以很容易地再導入到sigmaplot、origin等程序里重新繪圖。
2 自旋布居
2.1 基本概念和計算原理
自旋密度是三維函數，每個點一個數值。而我們往往想討論某個片段、某個原子、某個原子軌道帶多少單電子，這就需要做自旋布居分析(Spin population analysis)，有很多具體算法，分析結果定量不同但一般定性一致。其中，Mulliken、SCPA、Bickelhaupt、NPA等方法可以得到基函數、原子軌道、原子、分子片段的自旋布居數，而Hirshfeld、Becke、Voronoi、AIM方法只能得到原子和分子片段的自旋布居數。自旋布居數的定義是Alpha布居數減Beta布居數，比如某個原子的自旋布居數是0.3，就是說明這個原子帶的Alpha電子比它帶的Beta電子多0.3個。顯然正的布居數對應帶Alpha單電子，負的布居數對應帶Beta單電子。這里很簡單地提一下計算原理，詳細介紹和對比分析參見筆者的《分子軌道成分的計算》（化學學報,69,2393，http://sioc-journal.cn/Jwk_hxxb/CN/abstract/abstract340458.shtml）。Mulliken、SCPA、Stout-Politzer、Bickelhaupt這幾種方法原理類似，它們做自旋布居分析時可以認為是先通過某種方式計算出每個基函數的自旋布居，然后把每個原子的所有基函數的自旋布居加和就得到了原子的自旋布居，而原子的自旋布居再進一步加和就得到了分子片段的自旋布居。如果想得到原子軌道的自旋布居，就得弄清楚基函數和原子軌道的對應關系，比如6-31G*是每個價層原子軌道用兩個基函數來描述，故把那兩個基函數的自旋布居相加就得到的對應的原子軌道的自旋布居。這幾種方法分析速度都很快，Stout-Politzer和Bickelhaupt不推薦用，Mulliken和SCPA的結果通常合理，但關鍵要注意的是用的基組死活都不能有彌散函數，否則結果根本沒法用！
Hirshfeld、Becke、Voronoi和AIM方法是把分子空間以不同方式劃分成一個個原子空間，然后在空間中對自旋密度進行積分來得到原子自旋布居數。之后可以再加和成分子片段的自旋布居數。AIM計算耗時很高，也沒什么特別的好處，不建議用。Voronoi的物理意義不強，不建議用。對于只需要原子/片段的自旋布居數而不需要原子軌道的自旋布居數的時候，筆者十分推薦用Hirshfeld或Becke方法進行計算，原理十分簡單清晰，而且可靠性高，也不怕有彌散函數。
想必看過以上文字的人早已分清自旋密度與自旋布居的關系了。令筆者經常渾身難受的是老看到有人問“怎么算原子的自旋密度”，明擺著這些人根本沒搞懂最基本的概念。原子怎么算自旋密度？到底算哪個點的自旋密度？是原子核位置的還是原子核附近具體某個點的？如果要考察原子帶多少單電子，要計算的是“原子自旋布居”，那絕對不叫“原子自旋密度”！關于這點，筆者以前在《量子化學中的一些常見不良寫法和用詞》（http://www.shanxitv.org/298）就專門明確強調過了。
2.2 在Multiwfn中計算自旋布居
在Multiwfn中可以實現基于Mulliken、SCPA、Stout-Politzer、Bickelhaupt、Hirshfeld、Hirshfeld-I、Becke、Voronoi、AIM的自旋布居分析，下面就說說其中常用方法的具體計算流程。
2.2.1 Mulliken/SCPA/Lowdin等
實際上Gaussian在做開殼層體系計算的時候末尾直接就會輸出原子的自旋布居數，如Mulliken charges and spin densities:
1 2
1 C 0.162279 0.781603
2 H 0.139270 -0.012974
3 N -0.708466 0.140471
4 H 0.333784 0.035523
5 H 0.344143 0.001055
6 O -0.271011 1.054322
這里第二列就是原子自旋布居。可見Gaussian開發者多糊涂，居然把spin populations寫成了spin densities！而且，居然Exploring Chemistry With Electronic Structure Methods第三版里面作者在討論相關問題的時候也用spin density這個詞來說原子，真是嚴重誤人子弟！
在Multiwfn中，可以很方便地通過Mulliken、SCPA、Lowdin等方式做自旋布居分析，比Gaussian輸出的詳細多了。用Mulliken、SCPA、Lowdin分析時不能用.wfn/.wfx文件，必須含有基函數信息的文件，比如.fch、.molden等，詳見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）。這里還是以三重態甲酰胺為例。
啟動Multiwfn，輸入
formamide-m3.fch
7 //布居分析
5 //Mulliken分析（如果用SCPA，選7；用Lowdin，選6）
1
程序依次輸出以下內容：
Population of basis functions：基函數的布居數
Population of shells：基函數殼層的布居數
Population of each type of angular moment orbitals：每種角動量基函數的布居數
Population of atoms：原子的布居數
輸出中的Spin pop.那一列就是自旋布居，是Alpha pop.和Beta pop.之差。
原子布居輸出結果和Gaussian一致，我們來看看輸出的各個原子各角動量基函數的布居情況
Atom Type Alpha pop. Beta pop. Total pop. Spin pop.
1(C ) s 1.64782 1.51428 3.16210 0.13354
p 1.63876 0.99055 2.62931 0.64821
d 0.02309 0.02323 0.04631 -0.00014
2(H ) s 0.42388 0.43685 0.86073 -0.01297
3(N ) s 1.77929 1.78117 3.56046 -0.00188
p 2.12848 1.99169 4.12017 0.13680
d 0.01669 0.01114 0.02783 0.00556
4(H ) s 0.35087 0.31535 0.66622 0.03552
5(H ) s 0.32846 0.32740 0.65586 0.00106
6(O ) s 1.99532 1.98659 3.98191 0.00873
p 2.66628 1.62313 4.28940 1.04315
d 0.00107 -0.00137 -0.00030 0.00244
Total s 6.52563 6.36164 12.88727 0.16399
p 6.43352 4.60537 11.03888 1.82815
d 0.04085 0.03299 0.07384 0.00786
從數據中可見，體系總共的兩個單電子，其中1.828個都在體系中原子的p軌道上。具體來說，O6的p軌道貢獻了1.04個，C1的p軌和s軌分別貢獻了0.648和0.133個，在N3上只有極少量單電子。
各個原子軌道上自旋布居是多少？其實這也很容易得到，仔細看《利用布居分析判斷基函數與原子軌道的對應關系》（http://www.shanxitv.org/418）。我們一旦按照此文提到的過程判斷出了基函數與原子軌道的對應關系，就可以把Multiwfn輸出的基函數的自旋布居加和成原子軌道的自旋布居了。
經常有人問體系的磁矩怎么算，這里順帶說一下。分子的磁矩來自于電子自旋以及電子的軌道運動，前者是分子磁矩最主要的貢獻者。若體系的alpha電子數減beta電子數為n，并且將電子的自旋g因子近似為2（精確值是2.0023193），則分子的電子自旋磁矩的z分量（μZ）就為-n*μ_B，這里μ_B=e*h_bar/(2m_e)是玻爾磁子。因此，μZ可以分解成原子軌道、原子、片段或者三維空間中各個點的貢獻。比如某d原子軌道上自旋布居為0.4，就說這個d原子軌道對μZ貢獻是-0.4μ_B；若某點自旋密度為x，就可以說這個點的單電子對μZ的貢獻是-x*μ_B。所以，對于上面三重態甲酰胺的例子，μZ=-2*μ_B，我們可以說C1對體系μZ的貢獻是-0.78*μ_B，其中p軌道貢獻-0.64*μ_B，其它的是s軌道貢獻的。另外，電子磁矩的總大小|μS|的計算方式是2*μ_B*sqrt(S*(S+1))，此處體系的電子自旋量子數S等于alpha減beta電子數再除以2。|μS|顯然沒法分解為體系不同部分的貢獻。單重態體系S=0，故電子自旋磁矩總大小為0。而對于單重態雙自由基，總的電子自旋磁矩也為0，但可以視為有局部電子自旋磁矩。
如果在做布居分析之前先用主功能7里面的選項-1定義片段，那么做如上分析時，還會輸出片段的電荷以及片段的自旋布居。
2.2.2 Hirshfeld/Becke
在Multiwfn中使用Hirshfeld或Becke方法計算原子自旋布居時用.wfn、.fch、.molden等各種Multiwfn支持的含有波函數信息的文件皆可。還是用三重態甲酰胺的例子，啟動Multiwfn后輸入formamide-m3.fch
15 //模糊空間分析。默認是用Becke方式劃分原子空間，如果想用Hirshfeld方式，選-1再選3
1 //在各個原子空間內積分指定實空間函數
5 //自旋密度
結果如下
Atomic space Value % of sum % of sum abs
1(C ) 0.69370502 34.685254 34.685254
2(H ) 0.01750339 0.875170 0.875170
3(N ) 0.19357310 9.678656 9.678656
4(H ) 0.03433842 1.716921 1.716921
5(H ) 0.00590491 0.295246 0.295246
6(O ) 1.05497500 52.748754 52.748754
Value這一列就是原子自旋布居，和上一節Mulliken分析結果定性一致。% of sum這一列是相應行的數值占所有數值加和的百分比，可見O6貢獻了單電子的一半以上。
如果在主功能15里面先選-5，把某個片段里的原子定義為計算時被考慮的原子，那么再按照以上操作后，跟著輸出的Summing up above values的數值就是這個片段的自旋布居。
2.2.3 AIM
如前所述，筆者不推薦用AIM方法做自旋布居分析，不過姑且也說一下基本流程。也是用.wfn/.wfx/.fch/.molden作為輸入文件皆可。啟動Multiwfn后輸入
formamide-m3.fch
17 //盆分析
1 //產生盆
1 //用電子密度零通量面作為劃分盆的依據
2 //中等質量格點
7 //在盆中對指定實空間函數進行積分
5 //自旋密度
結果為
Atom Basin Integral(a.u.) Vol(Bohr^3) Vol(rho>0.001)
1 (C ) 4 0.67507872 376.036 74.955
2 (H ) 5 0.02216879 587.073 45.094
3 (N ) 2 0.20244172 553.973 115.300
4 (H ) 6 0.03091178 362.484 28.203
5 (H ) 1 0.00400347 312.129 27.089
6 (O ) 3 1.06545431 754.626 118.934
Sum of above integrals: 2.00005879
Integral下面的數值就是相應原子的自旋布居，和前面其它方法算得也相仿佛。這里顯示總積分值很接近2.0，說明結果的積分精度是沒問題的，如果偏離實際值很多，就說明積分格點設置不合理，需要在設定格點的那一步做適當調整。
3 其它
原子的自旋布居還可以通過對原子著色的方式直觀展現，見《使用Multiwfn+VMD以原子著色方式表現原子電荷、自旋布居、電荷轉移、簡縮福井函數》（http://www.shanxitv.org/425）。
如果你想了解自旋密度是怎么由軌道構成的話，那么一定要仔細看《用于非限制性開殼層波函數的雙正交化方法的原理與應用》（http://www.shanxitv.org/448）。簡單來說，對于限制性開殼層(RO)計算，自旋密度就是各個單占據軌道自旋密度的加和，自旋布居也相當于是這些單占據軌道中被考察對象所占的成份的加和。但是RO計算耗時高、局限性大、難收斂、結果差，因此一般我們都是用非限制性開殼層(U)形式計算開殼層體系，但此時所有軌道都會對自旋密度有所貢獻，這就很難從軌道角度考察自旋密度的構成了。不過，將這些非限制性軌道用Multiwfn做雙正交化變換后，就只有極少幾條軌道會對自旋密度有所貢獻，討論起來就方便多了。
自旋密度還可以通過自旋自然軌道的方式考察，這是一種比較獨特的考察方法，見《在Multiwfn中基于fch產生自然軌道的方法與激發態波函數、自旋自然軌道分析實例》（http://www.shanxitv.org/403）。