使用Multiwfn計算Bond length/order alternation (BLA/BOA)和考察鍵長、鍵級、鍵角、二面角隨鍵序號的變化
使用Multiwfn計算Bond length/order alternation (BLA/BOA)和考察鍵長、鍵級、鍵角、二面角隨鍵序號的變化
文/Sobereva@北京科音
First release: 2019-Aug-11 Last update: 2019-Nov-30
1 基本知識
共軛寡聚物是量子化學經常研究的對象,如10聚乙烯、12聚噻吩等。這類體系都有共軛長鏈,pi電子能夠在上面長程離域。這種共軛鏈上的鍵長通常呈現交替變化特征,即長-短-長-短...,相應地,鍵級的變化也有類似的交替波動規律。為展現這點,可以繪制鍵長或鍵級隨共軛鏈上的鍵的序號變化的折線圖,例如ACS Cent. Sci., 2, 309 (2016)文中的圖2。還專門有人定義了Bond length alternation (BLA)和Bond order alternation (BOA)來定量化描述這種交替特征的大小。按照Handbook of Thiophene-based Materials: Applications in Organic Electronics and Photonics公式7.6的定義,BLA一般化的定義可表達為:[偶數鍵的平均鍵長] - [奇數鍵的平均鍵長]。假設共軛鏈的始端到末端的原子序號順序是3-5-6-9-10-12,那么BLA就是[R(5-6)+R(9-10)]/2 - [R(3-5)+R(6-9)+R(10-12)]/3。
將BLA里的鍵長改為鍵級,就對應于BOA。相對于BLA從幾何層面上體現鍵的交替特征,BOA是從電子結構層面上體現這點。由于對于同類鍵,鍵越長鍵級通常越小,因此BLA和BOA一般呈負相關(BLA為正時BOA一般為負,BLA越正則BOA傾向于越負)。如《Multiwfn支持的分析化學鍵的方法一覽》(http://www.shanxitv.org/471)的鍵級部分所介紹的,鍵級有很多不同的定義。Multiwfn在計算BOA的時候用的是Mayer鍵級,因為其計算速度很快,普適性強,數值大小和形式鍵級接近,有較好化學意義,適合用于作為BOA中的鍵級定義。原理上,用筆者在J. Phys. Chem. A, 117, 3100 (2013)提出的拉普拉斯鍵級(LBO)來考察BOA也完全可以,只不過耗時會高一些。
BLA和BOA與共軛寡聚物的gap、(超)極化率等很多屬性都有十分密切的關系,如BLA越大gap越大,還專門有人擬合了關系式。外電場、取代基也會影響BLA/BOA值。更多相關內容這里就不細說了,在這里有較全面的介紹:http://photonicswiki.org/index.php?title=Structure-Property_Relationships。很多此類體系的研究文章中也都討論了BLA或BOA,如J. Chem. Phys., 136, 094904 (2012)、Syn. Met., 141, 171 (2004)等等。
順帶一提,在寡聚物的共軛鏈上,除了鍵的特征有規律性的波動,原子的特征也有規律性的波動,如筆者的RSC Adv., 3, 25881 (2013)一文中圖5考察了納米線的共軛鏈上原子電荷和原子的pi布居的變化。
研究共軛鏈的特征時候還經常考察鍵角、二面角隨共軛路徑的變化。路徑上的各二面角偏離平面的程度整體越低,往往暗示共軛越強。
2 Multiwfn的BLA/BOA的計算功能
Multiwfn是功能非常全面的波函數分析程序,可以從官網http://www.shanxitv.org/multiwfn免費下載,相關知識見《Multiwfn FAQ》(http://www.shanxitv.org/452)。從2019-Aug-11更新的Multiwfn開始,BLA和BOA計算功能作為Multiwfn的主功能200的子功能18出現。Multiwfn還順便給出選定的共軛鏈上各個鍵的鍵長和鍵級,便于用戶之后畫出“鍵長-鍵序號”和“鍵級-鍵序號”的折線圖。此功能用的鍵級是Mayer鍵級,注意由于Mayer鍵級怕彌散函數,所以用的波函數文件在產生的時候不能帶彌散函數。
如果使用含有基函數信息的文件比如.fch、.molden作為輸入文件,那么鍵長和鍵級數據都會被給出來;如果使用只含分子結構信息的格式如.xyz、.pdb、.mol、.mol2作為輸入文件也可以,但此時鍵級數據和BOA值不給出來。更多相關信息,以及怎么產生這些文件見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》(http://www.shanxitv.org/379)。
雖然手動提取鍵長、計算鍵級,然后再手算BLA/BOA也不難,但是當考察的體系很多,特別是聚合度很高的話,手動處理非常費時費力,還容易弄錯,而Multiwfn的這個功能大大簡化了這個過程。
此外,在這個功能中用戶還可以要求Multiwfn將各個鍵角、二面角的數值沿著路徑的變化輸出出來,從而便于繪制“鍵角-鍵角序號”和“二面角-二面角序號”的折線圖,從后面的例子就會看到這種圖對展現結構特征很有價值。
3 例子
3.1 考察五聚噻吩的BLA、BOA以及鍵長/鍵級隨鍵序號的變化
下面以五聚噻吩為例,介紹一下如何在Multiwfn中計算BLA、BOA以及繪制鍵長/鍵級隨鍵序號的變化。此體系以下簡稱TP5,用Gaussian 16在PBE0/6-31G*下優化并得到的fchk文件可以在此下載:http://www.shanxitv.org/attach/501/TP5.zip。值得一提的是,這個體系看似好像應該是平面結構,但實際上如TP5.fchk所示,其無虛頻的結構是稍微彎折的。
在計算之前,顯然必須先確定出來從共軛鏈的始端到末端的原子序號順序。雖然可以在可視化程序里通過觀看分子,把序號一個一個記錄下來,但這樣很費勁。利用Multiwfn算BLA/BOA則不需要這么麻煩。用戶只需要告訴Multiwfn這條鏈上都有哪些原子,以及始端和末端的原子序號是什么就行了,Multiwfn會根據這些信息以及原子間連接關系自動識別出鏈上的原子序號順序。
獲取鏈上原子序號最省事的做法是把體系載入GaussView,選Builder面板上的Select Atoms by Brush按鈕(對于gview 6.0.16版是倒數第四行最右邊那個。gview 6之前的版本沒有這個按鈕),然后按住鼠標左鍵,讓光標經過共軛鏈上的每一個原子來選中它們,之后窗口里的圖像應該是下面這樣
注:由于C-S、C-H鍵不參與共軛,所以H、S原子不納入考慮。由于當前體系里所有碳原子都參與共軛,所以你在gview的Atom List Editor里直接把所有碳原子選上效果也是等價的。
然后進入gview的Tools - Atom selection,會看到當前被選中的原子序號范圍是:1-4,9-10,12,14,16-17,19,21,23-24,26,28,30-31,33,35。從上圖還可以看到,選中的鏈的始端的原子序號是1,末端的原子序號是35。
現在可以開始計算了。啟動Multiwfn,然后輸入
TP5.fchk
200 //主功能200
18 //計算BLA/BOA和考察鍵長/鍵級隨鍵序號的變化
1-4,9-10,12,14,16-17,19,21,23-24,26,28,30-31,33,35 //把前面得到的鏈當中的原子序號范圍直接粘貼到窗口里
1,35 //鏈的始端和末端原子序號
屏幕上首先顯示出了自動判斷的共軛鏈上的原子序列:
Sequence of the atoms in the chain from the beginning side to the ending side
1 2 3 4 9 10 12 14 16
17 19 21 23 24 26 28 30 31
33 35
你可以對照分子結構圖檢驗一下判斷的對不對,不對的話后面的結果都沒意義,當前識別出來的是對的。原子間連接關系是自動根據鍵長和共價半徑判斷的,細節這里提了《談談原子間是否成鍵的判斷問題》(http://www.shanxitv.org/414)。如果你載入的是.mol或.mol2文件,則連接關系是直接從里面讀取而不是自動判斷的。
接下來輸出的是共軛鏈上從始端到末端的各個鍵的序號、構成這個鍵的相應原子,以及這個鍵的鍵長和Mayer鍵級:
Bond Atom1 Atom2 Length (Angstrom) Mayer bond order
1 1 2 1.3678 1.6303
2 2 3 1.4232 1.2924
3 3 4 1.3795 1.5445
4 4 9 1.4468 1.0985
5 9 10 1.3799 1.5273
6 10 12 1.4153 1.2976
7 12 14 1.3811 1.5164
8 14 16 1.4427 1.1098
9 16 17 1.3814 1.5139
10 17 19 1.4143 1.3017
11 19 21 1.3814 1.5140
12 21 23 1.4427 1.1097
13 23 24 1.3811 1.5164
14 24 26 1.4153 1.2976
15 26 28 1.3799 1.5273
16 28 30 1.4468 1.0985
17 30 31 1.3795 1.5445
18 31 33 1.4232 1.2924
19 33 35 1.3678 1.6303
如屏幕上的提示所示,上述數據還被輸出到了當前目錄下的bondalter.txt中。
之后輸出的是統計數據,包括序號為偶數和奇數的鍵的數目,這兩類鍵的平均鍵長、平均鍵級,以及相應的BLA和BOA值。
The number of even bonds: 9
The number of odd bonds: 10
Average length of even bonds: 1.4300 Angstrom
Average length of odd bonds: 1.3779 Angstrom
Bond length alternation (BLA): 0.0521 Angstrom
Average bond order of even bonds: 1.2109
Average bond order of odd bonds: 1.5465
Bond order alternation (BOA): -0.3356
如果想繪制“鍵長 vs. 鍵序號”和“鍵級 vs. 鍵序號”的折線圖,就把bondalter.txt拖入到Origin里繪制即可。下圖是繪制后的效果,兩個折線圖作為兩個不同的layer繪制在了一起,此圖的Origin .opj文件在前述的TP5.zip文件包里也提供了。
之后程序還問你是否把鍵角、二面角在路徑上的變化輸出出來,此例我們不考察這個,因此輸入n。對于其它共軛寡聚物的考察方法與本文類似,請大家自行嘗試。
后來有讀者問,在這個例子里,gview的Atom selection里已經顯示了1-4,9-10,12,14,16-17,19,21,23-24,26,28,30-31,33,35這個序列,干嘛Multiwfn還要再自己判斷一次序列?對當前這個體系,確實Multiwfn這一步沒必要,因為原子序號本身就是順著的;但是對于其它體系,很多情況下,你要考察的鏈上的原子序號順序不是順著的,比如有時可能是6,5,2,3,4,9,10這樣交錯的。Multiwfn自動判斷鏈上的原子順序的好處是,你在給Multiwfn提供鏈上的原子序號時不需要考慮提供的順序,你按照實際順序寫成6,5,2,3,4,9,10也行,寫成2-6,9-10也行,寫成9-10,2-6也行,怎么寫都無所謂,只要序號范圍對即可。這就省得你一邊看結構一邊一個個記錄鏈上的原子序號了,對于鏈比較長的情況省事多了,你想考察哪條鏈只要在gview里用前述功能一劃來選中,在Atom selection里拷出來序號范圍就行了。
3.2 考察18碳環上的鍵角、二面角變化
在《談談18碳環的幾何結構和電子結構》(http://www.shanxitv.org/515)中筆者專門介紹了18碳環體系。這個體系在勢能面極小點情況下鍵角都是160度,二面角都是0度,但是在實際環境中,由于有熱運動,體系結構會發生明顯變形。筆者對這個體系在200 K下做了從頭算分子動力學,隨便截取了其中一幀,結構如下所示。當前這個例子就用Multiwfn來考察一下此結構中的鍵角、二面角沿著環是怎么變化的。
啟動Multiwfn然后輸入
examples\C18_MD_1.xyz
200
18
1-18 //即當前體系里所有原子的序號
1,1 //對于路徑是個封閉的環的情況,這一步輸入的兩個原子序號必須相同。輸入環上哪個原子序號都可以,這里我們將此體系的1號原子作為起始原子
此時屏幕上輸出了一共18個鍵長。之后選擇y來把鍵角、二面角在環上的變化輸出出來,看到以下信息:
Note The unit of printed values is degree
Atoms: 1 2 3 Angle: 164.841
Atoms: 2 3 4 Angle: 154.976
Atoms: 3 4 5 Angle: 155.200
Atoms: 4 5 6 Angle: 158.564
Atoms: 5 6 7 Angle: 167.583
Atoms: 6 7 8 Angle: 161.476
Atoms: 7 8 9 Angle: 155.205
Atoms: 8 9 10 Angle: 157.061
Atoms: 9 10 11 Angle: 161.160
Atoms: 10 11 12 Angle: 157.517
Atoms: 11 12 13 Angle: 162.158
Atoms: 12 13 14 Angle: 165.757
Atoms: 13 14 15 Angle: 155.054
Atoms: 14 15 16 Angle: 158.404
Atoms: 15 16 17 Angle: 159.228
Atoms: 16 17 18 Angle: 158.207
Atoms: 17 18 1 Angle: 162.906
Atoms: 18 1 2 Angle: 158.520
Atoms: 1 2 3 4 Dihedral: 28.69, deviation to planar: 28.69
Atoms: 2 3 4 5 Dihedral: 25.51, deviation to planar: 25.51
Atoms: 3 4 5 6 Dihedral: 22.98, deviation to planar: 22.98
Atoms: 4 5 6 7 Dihedral: 18.44, deviation to planar: 18.44
Atoms: 5 6 7 8 Dihedral: 27.47, deviation to planar: 27.47
Atoms: 6 7 8 9 Dihedral: 12.77, deviation to planar: 12.77
Atoms: 7 8 9 10 Dihedral: 5.86, deviation to planar: 5.86
Atoms: 8 9 10 11 Dihedral: 8.16, deviation to planar: 8.16
Atoms: 9 10 11 12 Dihedral: 3.88, deviation to planar: 3.88
Atoms: 10 11 12 13 Dihedral: 15.47, deviation to planar: 15.47
Atoms: 11 12 13 14 Dihedral: 16.28, deviation to planar: 16.28
Atoms: 12 13 14 15 Dihedral: 7.30, deviation to planar: 7.30
Atoms: 13 14 15 16 Dihedral: 3.10, deviation to planar: 3.10
Atoms: 14 15 16 17 Dihedral: 8.71, deviation to planar: 8.71
Atoms: 15 16 17 18 Dihedral: 12.94, deviation to planar: 12.94
Atoms: 16 17 18 1 Dihedral: 18.89, deviation to planar: 18.89
Atoms: 17 18 1 2 Dihedral: 16.91, deviation to planar: 16.91
Atoms: 18 1 2 3 Dihedral: 18.13, deviation to planar: 18.13
可見,總共18個鍵角,以及18個二面角都給出了。其中deviation to planar指的是偏離平面的程度。Multiwfn給出的二面角范圍是0~180,對于0~90的情況,deviation to planar就等于二面角;而對于90~180的情況,deviation to planar等于180減去這個角度。
之后如果想作圖的話,按照Multiwfn手冊5.5節說的做法把相應的列從屏幕上拷出來,可以拷到Origin、Excel之類的窗口里,之后做成折線圖就是下面的樣子
此圖非常直觀地體現出,在當前結構下,鍵角普遍明顯偏離18碳環極小點結構對應的160度,說明發生了顯著的環平面內的變形;從二面角的圖上可以看到數值普遍明顯大于0,說明18碳環也發生了顯著的平面外扭曲。可見這個分析對于展現某些體系的某些環的幾何特征非常有價值而且非常方便。
在《一篇文章深入揭示外電場對18碳環的超強調控作用》(http://www.shanxitv.org/570)介紹的筆者的研究論文中,筆者還用此文介紹的方法考察了不同電場下18碳環結構的鍵長、鍵角交替曲線,如下所示。這個研究的題材很有意思,建議大家讀讀。