計算分子動力學軌跡中兩個環平面間的距離和夾角

計算分子動力學軌跡中兩個環平面間的距離和夾角

文/Sobereva@北京科音  2021-Mar-3

1 前言

在《全面探究18碳環獨特的分子間相互作用與pi-pi堆積特征》（http://www.shanxitv.org/572）曾介紹了筆者做的關于18碳環弱相互作用的研究，其中跑了18碳環二聚體的從頭算動力學。當時為了統計模擬過程中兩個碳環的間距和相對角度隨模擬時間的變化，筆者在Multiwfn里寫過一個子程序來實現。今天有人在思想家公社QQ群里問怎么獲得環平面間夾角隨模擬時間的變化，我決定干脆把Multiwfn里的這個隱藏功能給大家簡單介紹一下，順便又增加了計算環平面與笛卡爾平面間夾角的功能，可能有人用得著。請讀者使用2021-Mar-3及以后更新的Multiwfn，否則會與本文不符。Multiwfn可以在http://www.shanxitv.org/multiwfn免費下載。

使用本文的功能需要用xyz軌跡文件作為輸入，格式見《談談記錄化學體系結構的xyz文件》（http://www.shanxitv.org/477）。如果用ORCA或CP2K等程序跑動力學，直接就有了xyz軌跡；如果用GROMACS、AMBER等程序跑，可以用VMD載入后再保存為xyz軌跡。

2 計算兩個環平面間夾角隨時間的變化

這里就用18碳環二聚體在100 K下的動力學軌跡作為演示。此軌跡是我用ORCA跑的，方法參見《使用ORCA做從頭算動力學(AIMD)的簡單例子》（http://www.shanxitv.org/576），一共4001幀，文件可以在這里下載：http://www.shanxitv.org/attach/590/C18dimer.rar。

啟動Multiwfn，然后輸入
C18dimer.xyz
1000 //隱藏的主功能
201 //隱藏的子功能
1  //獲得片段間幾何中心距離和夾角
4001  //軌跡總幀數
1-18  //第一個片段里的原子序號，即第一個18碳環
19-36  //第二個片段里的原子序號，即第二個18碳環

瞬間就算完了。在當前目錄下出現了distangle.txt，每一列的內容是
第1列：幀號
第2列：兩個片段的幾何中心間的距離（埃）
第3列：第一個片段平面的法矢量與兩個片段幾何中心連線之間的夾角，姑且可以叫錯位角，體現相對錯位程度，示意圖見下
第4列：兩個片段的法矢量之間的夾角，即一般意義上的平面間夾角

將distangle.txt中的數據放到一起繪圖，得到下圖

注意程序給出的夾角的范圍是[0,90]。Multiwfn是通過最小二乘法來得到與指定的片段結構盡可能吻合的平面，然后再利用其法矢量算夾角。這比起一些程序通過三個原子來定義平面合理多了，尤其是諸如18碳環這種環的柔性較高、在實際模擬過程中往往顯著偏離平面的體系。

3 計算環平面與笛卡爾平面的夾角隨時間的變化

例如有個軌跡文件叫Mio_Akiyama.xyz，想考察其中1,4-9,15,18這些原子構成的片段與XY、YZ、XZ平面的夾角隨時間的變化，就在啟動Multiwfn后輸入
Mio_Akiyama.xyz
1000 //隱藏的主功能
201 //隱藏的子功能
2  //獲得片段與笛卡爾平面間的夾角
4001  //軌跡總幀數
1,4-9,15,18  //片段里的原子序號
然后在當前目錄下出現了angle.txt。第一列是幀號，第2、3、4列分別是這個片段擬合出的平面與XY、YZ、XZ平面間的夾角，范圍是[0,90]。