使用Multiwfn研究分子動力學中的弱相互作用
后記：針對aNCI方法用于蛋白質-配體相互作用的分析，筆者后來專門寫了一篇文章《使用Multiwfn做aNCI分析圖形化考察動態過程中的蛋白-配體間的相互作用》（http://www.shanxitv.org/591）。本文的讀者千萬別忘了看此文！
使用Multiwfn研究分子動力學中的弱相互作用
文/Sobereva @北京科音
First release: 2013-May-18 Last update: 2021-Jul-21
1 前言
曾經筆者在《使用Multiwfn圖形化研究弱相互作用》（http://www.shanxitv.org/68）中詳細介紹了楊偉濤課題組在JACS,132,6498 (2010)提出的研究弱相互作用的約化密度梯度(RDG)方法，也叫NCI方法。原先的RDG方法只考慮單個結構下的弱相互作用，然而由于熱運動，實際體系的結構是不斷波動的，其中的弱相互作用也因此不可能只通過單一的結構全面表現出來。在JCTC,9,2226 (2013)中，楊偉濤課題組提出了averaged RDG (aRDG)方法（也叫aNCI方法），它是原先RDG方法的擴展，使得RDG分析方法可以用于分析多幀結構，尤為適合結合分子動力學模擬技術來研究平衡的動態環境中的弱相互作用。雖然aRDG原則上說也可以結合諸如蒙特卡羅模擬，但下文討論的只限于aRDG方法對全原子分子動力學軌跡的分析。在《一篇最全面介紹各種弱相互作用可視化分析方法的文章已發表！》（http://www.shanxitv.org/667）介紹的筆者的綜述文章中有對aRDG非常充分詳細的介紹并給了諸多應用例子，強烈建議大家閱讀。
已經有大量文章使用了Multiwfn的aRDG分析發表了研究論文，可以參考，比如
ACS Omega (2020) DOI: 10.1021/acsomega.0c02863
Carbohyd. Res., 496, 108134 (2020) DOI: 10.1016/j.carres.2020.108134
Russ. J. Phys. Chem. A, 94, 1356 (2020) DOI: 10.1134/S0036024420070195
J. Mol. Liq., (2020) DOI: 10.1016/j.molliq.2020.113305
Chin. Phys. B, 27, 083103 (2018) DOI: 10.1088/1674-1056/27/8/083103
PLoS ONE, 13, e0196651 (2018) DOI: 10.1371/journal.pone.0196651
Sci. Rep., 6, 21763 (2016) DOI: 10.1038/srep21763
Sci. Rep., 5, 7572 (2014) DOI: 10.1038/srep07572
2 原理
如果不熟悉RDG方法，請務必先閱讀《使用Multiwfn圖形化研究弱相互作用》。aRDG方法和原先的RDG方法在原理上十分類似，同樣是通過RDG函數的等值面來展現弱相互作用區域，并且將sign(lambda2)rho函數用不同的顏色投影到RDG等值面上來展現出弱相互作用類型。但是在aRDG方法中，計算RDG和sign(lambda2)rho函數所用的電子密度、電子密度的梯度以及電子密度的Hessian矩陣都是來自對多幀結構取平均得到的。也就是說，比如要研究的分子動力學軌跡有1000幀，那么aRDG所用的密度就是對這1000幀都依次計算密度然后取平均得到。因此aRDG中的RDG函數和sign(lambda2)rho函數都表現的是整條軌跡中平均的RDG和平均的sign(lambda2)rho。
除了這點區別外，aRDG方法還定義了一個熱波動指數(thermal fluctuation index, TFI)，TFI=std(rho)/avg(rho)。其中avg(rho)就是上面說的整個軌跡中的平均密度，std(rho)就是整個軌跡中密度的方差，一個函數的方差越大就表明函數的波動越強。熱波動指數的用處是展現不同區域弱相互作用的穩定程度，分析的時候將熱波動指數通過不同顏色映射到平均的RDG等值面上，通過考察顏色，就能判斷不同區域弱相互作用形成的穩定與否。
3 在Multiwfn中的使用
aRDG方法可通過Multiwfn的主功能20里的子功能3來實現。Multiwfn可在http://www.shanxitv.org/multiwfn上免費下載。
由于aRDG分析用的是平均的密度及其梯度和Hessian，所以分析結果會隨著考慮的幀數增加而逐漸趨于收斂。考慮的幀數越多，結果越精確，圖像也越細膩。如果只考慮一幀的結構，那么和原先的RDG分析就沒區別了。建議至少考慮500幀結構，如果想作出漂亮的圖，應該用1000幀或更多。
使用Multiwfn的aRDG分析功能時，需要在一開始載入.xyz軌跡文件。.xyz格式很簡單，.xyz軌跡文件只不過是簡單地將一般的只記錄單個結構的.xyz文件合并在一起而已，所以算是個比較通用的記錄軌跡的格式。用戶可以用任何自己最常用的分子動力學程序來跑出軌跡，然后載入VMD程序，之后直接轉換為.xyz格式軌跡文件。
在軌跡中，感興趣的分子的坐標應當保持固定不變，這可以在動力學程序中將其坐標凍結來實現。并且感興趣的分子最好是在整個體系的靠中間的位置。計算格點數據時設定的計算區域應當能夠容納這個感興趣的分子，且四周有一定延展區域，以避免重要區域的RDG等值面被截斷。
由于需要考慮幾百、上千幀，所以計算那么多幀的電子密度、梯度及Hessian的格點數據是不可能直接基于波函數的，因為無論是計算那么多次波函數還是計算那么多格點數據都是極其耗時的。aRDG原文用的是QM/MM方法，感興趣的區域靠QM來描述并基于波函數獲得密度，而外環境部分靠分子力學描述。但QM/MM動力學在實際操作上并不方便，對很多用戶來說都有不小門檻，而且依然較耗時，且輸出的格式與動力學程序的依賴性也太強因此不夠普適。所以，Multiwfn中的aRDG分析一律使用promolecular近似，也就是基于自由原子密度的疊加來近似產生電子密度及其梯度和Hessian。這個近似顯著簡化了計算過程，大大節約了計算量，而且分析結果依然比較可靠。
4 實例：標況下純水體系中水分子之間的平均弱相互作用
限于時間和精力，此文只給出一個應用aRDG的例子，也就是分析常溫常壓下的純水體系中水分子之間的弱相互作用。對于其它體系的研究請根據本例舉一反三。本節涉及到的所有文件包括軌跡都可以在此處下載：http://www.shanxitv.org/attach/186/aRDG_files.rar。
4.1 生成軌跡
用任何程序來獲得動力學軌跡都是可以的，本文用gromacs 4.5版來獲得。如果讀者習慣使用其它動力學程序，下文的gromacs的具體指令就不用管了，只要弄明白每一步的目的就行了。注意從GROMACS 5.0開始，操作步驟和本文的過程大不一樣。
首先建立一個空盒子文件emptybox.gro，盒子各邊長都為2.5nm，內容如下：
WAT
0
2.5 2.5 2.5
盒子邊長當然可以更長，只不過計算要更費時間。2.5nm算是底限了，再小的話就難以保證盒子長度是非鍵相互作用cut-off距離的兩倍以上了。
然后運行以下命令往盒子里面填充水得到純水盒子，其中將含有511個水分子
genbox -cp emptybox.gro -cs spc216.gro -o water.gro
運行此命令，并且選GROMOS96 53a6 force field，得到水盒子體系的top文件。水分子用的是SPC/E模型。
pdb2gmx -f water.gro -o water.gro -p water.top -water spce
接下來跑100ps NPT動力學，令水盒子在298.15K 1atm下平衡。范德華作用的cut-off為1.2nm。靜電相互作用用PME方法計算，實空間cut-off距離為1.2nm。步長0.001ps。
grompp -f pr.mdp -c water.gro -p water.top -o water-pr.tpr
mdrun -v -deffnm water-pr
現在用VMD打開water-pr.gro，從中選一個比較靠近盒子中央的水分子，我們選取101號水分子，下圖中將它用vdW方式突出顯示
注意101號水分子的三個原子編號是301、302和303，氧原子是301號原子。
這個水分子的坐標在接下來的動力學過程中將保持凍結。為此，先生成index文件，輸入以下命令
make_ndx -f water-pr.gro
ri 101 // 之所以是ri而不是r是因為ri是從1開始編號的，r是從0開始編號的
q
得到的index.ndx里的group r_101便是101號水分子。
運行以下命令來做298.15K下的1ns的平衡動力學，每隔1ps保存一幀結構，共得到1000幀。注意用NVT系綜，因為如果允許控壓的話，會對盒子進行scale從而影響101號水的坐標。由于之前已經在1atm下平衡了，水的可壓縮性又很低，所以改用NVT模擬依然合理。101號水的坐標靠freezegrps = r_101關鍵詞來凍結住。
grompp -f md.mdp -c water-pr.gro -p water.top -o water-md.tpr -n index.ndx
mdrun -v -deffnm water-md
將water-pr.gro和剛生成的water-md.xtc載入VMD，核對一下是否101號水確實固定住了。然后在VMD main窗口中點擊與這個體系對應的項，并在其上點擊右鍵，選Save coordinate，將File type設為xyz，Selected atoms選all，Frames框當中的First和Last分別寫1和1000，然后點Save按鈕，將軌跡保存為wat.xyz。
特別注意：有非常重要的一點是，.xyz文件對各個原子應當記錄的是元素名，而不應當是原子名。但是，按照上述方式得到的.xyz軌跡里記錄的卻是原子名，如果你用文本編輯器打開wat.xyz就會看到里面每個水包含OW、HW1、HW2，這都是原子名。一般情況下，做本文的分析之前，應當手動把這些原子名替換為對應的元素名，這樣才能100%確保Multiwfn能正常分析。但是，由于周期表里沒有叫做OW、HW1、HW2的元素，因此Multiwfn在讀取的時候會進一步試圖通過首字母判斷元素，因此Multiwfn是可以恰當地將這三種原子判斷為氧原子和氫原子的，因此接下來的分析也沒有問題。但如果Multiwfn載入.xyz文件后，屏幕上一開始提示的Formula后面的化學組成里出現了體系里不存在的元素，那么就說明必須把相應原子名替換為元素名了。比如.xyz里如果有原子類型為CA的碳，如果被Multiwfn載入，就會被當做鈣離子對待，肯定結果將和正常情況明顯不符。
4.2 用Multiwfn生成格點數據
啟動Multiwfn，輸入以下命令
wat.xyz
20 // 主功能20
3 // aRDG分析
1,1000 // 分析的幀號范圍是1~1000
7 // 通過原子編號自定義盒子中心，并自定義各方向盒子尺寸和格點數
301,301 // 輸入兩個原子編號，它們的中點即為盒子中心。這里輸入的兩個原子編號都是301（被固定的水的氧原子編號），也就是說將301號原子作為盒子中心
80,80,80 // 三個方向格點數都設為80
4.5,4.5,4.5 // 盒子在xyz方向的前后延展距離都為4.5 Bohr，因此盒子尺寸是9*9*9 Bohr。盒子設太大的話會使格點間距較大，導致圖像粗糙；如果太小，則會截斷RDG等值面。所以設置的大小應該根據實際體系憑借直覺來定，如果不確定的話可以先試幾個值，且每次都只考慮很少的幀數（如前20幀）來節約時間
之后，Multiwfn開始計算每一幀的電子密度及其梯度和Hessian，并取平均，之后用平均量計算出平均的RDG和sign(lambda2)rho。整個過程需要花一陣子時間，在目前一般的四核機子上耗時約半個小時到一個小時。如果是更大的動力學體系則更費時間，所以可以先暫時干別的去。
計算完畢后可以選選項1來查看平均RDG vs 平均sign(lambda2)rho散點圖（分析方法和《使用Multiwfn圖形化研究弱相互作用》介紹的一致，這里就不談了），也可以用相應選項將其保存到圖形文件中或者把數據點輸出到文本文件里。此體系的散點圖如下
選擇6將平均RDG和平均的sign(lambda2)rho分別導出到當前目錄下的avgRDG.cub和avgsl2r.cub。
由于我們要考察弱相互作用的穩定性，所以選7來計算熱波動指數并將之輸出到當前目錄下的thermflu.cub里。計算熱波動指數需要重新計算一遍每一幀的密度，所以也是比較耗時的，和計算平均RDG及sign(lambda2)rho過程耗時基本一樣。
4.3 分析結果
本文使用的是VMD 1.9來觀看結果。在本節開頭提到的壓縮包中提供了avgRDG.vmd和avgRDG_TFI.vmd，這是VMD作圖腳本。將它拷到VMD程序目錄下，并且將avgRDG.cub、avgsl2r.cub和thermflu.cub也都拷到VMD程序目錄下。
啟動VMD，在控制臺窗口輸入source avgRDG.vmd，就可以顯示出平均RDG等值面，isovalue為0.25。同時，平均sign(lambda2)rho函數也通過不同顏色投影到了等值面上。默認情況下所有分子都以CPK模式顯示了出來，這給分析101號水分子周圍的弱相互作用帶來了不便，所以進入Graphics-Representation，選Style為CPK的那一欄，將Selected Atoms里的內容改為101號水分子對應的三個原子，即改為serial 301 302 303然后點回車，此時就只有101號水分子顯示出來了。經過視角的調整，結果如下所示（在VMD里點t之后可以平移視角，點r可以恢復旋轉視角模式。點c再點擊一個原子可以將視角旋轉中心設在相應原子上）
期望的結果在圖上都展現了，等值面也沒有被截斷，說明盒子設得合適。不過不幸的是，在格點數據計算區域的邊緣，有一大堆零零碎碎亂七八糟的等值面，這和我們感興趣的弱相互作用無關卻十分礙眼。利用Multiwfn強大的格點數據處理功能，可以將這些離特定的一些原子（此例即101號水分子的三個原子）比較遠的區域的等值面給消掉。
啟動Multiwfn并輸入：
avgRDG.cub // 即剛才生成的平均RDG的格點數據
13 // 主功能13，用于處理格點數據
13 // 設定距離特定原子比較遠的格點的數值
1.5 // 如果格點距離特定原子的距離超過相應原子范德華半徑的1.5倍，則格點的數值將被設為下面用戶輸入的值。范德華半徑倍數的設定是依靠經驗的，可以多試幾次并檢驗結果以達到滿意效果
100 // 新設定的值為100，這個值是隨意的，只要是個明顯大于作RDG等值面圖時的isovalue的值就可以將這些點屏蔽掉
2 // 讀入特定原子的編號
301-303 // 101號水分子的三個原子的編號
0 // 將當前已經被改過的格點數據保存到新的cube文件中
avgRDG.cub // 新輸出的格點數據文件名
然后將新得到的avgRDG.cub移動到VMD目錄下覆蓋原先的avgRDG.cub（可以將原先的avgRDG.cub先備份一下）。然后重新用avgRDG.vmd腳本作圖并做前述的調整，結果如下所示（給出了兩個視角）
這次的圖非常清楚。默認的色彩刻度是-0.25~0.25，對應“藍色-綠色-紅色”的變化。越藍的地方說明靜電、氫鍵效應越強，越紅的地方說明位阻效應越明顯，而綠色表明相應位置平均密度值較低，對應于很弱的作用，即色散作用。在兩個O-H鍵對著的地方都有兩個藍色小圓片，表現出在動態環境下，O-H能與其它水形成很顯著的氫鍵。圖中有個顯著的，胖次形狀的綠色等值面，這清楚地展現了水分子傾向于以這些方向和其它水形成范德華作用。在氧的上方有一大坨等值面，其中兩側是藍色，中間偏紅，其藍色區域體現出氧原子傾向于以這樣的方位作為氫鍵受體，而紅色則在一定程度體現出實際環境中氧容易以這樣的方向和其它的水接觸而表現出互斥效應。氧上方的那一大坨等值面可以看出是由兩層值面相連構成的，如果將isovalue設低一些，如0.18，那么兩層等值面就分開了，如下所示
之所以會有兩層等值面原因尚不明，但可能分別與第一配位層、第二配位層分子的作用有聯系。最外層的等值面的形狀偏離對稱性相對顯著一些，而且等值面上投影的顏色有些“斑駁”，這說明對于這樣靠外區域的弱相互作用，只考慮1000幀還不足以使結果充分收斂。
關掉之前的等值面圖，然后運行source avgRDG_TFI.vmd，并作適當調整，結果如下
默認色彩刻度范圍是0~1.5，還是對應藍-綠-紅的變化。越偏藍的地方，說明熱波動指數越小，弱相互作用在動力學過程中越穩定；越紅的地方熱波動指數越大，弱相互作用在動力學過程中就越不穩定。可見O-H對著的小圓片顏色稍微發藍，這表明水分子靠O-H作為氫鍵給體與其它水分子構成的氫鍵在動力學過程中是很穩固的，很不容易被破壞。而胖次形狀的等值面整體呈紅色，說明水分子間的范德華作用相對于氫鍵是明顯不易維持穩定的，這容易理解，畢竟范德華作用是很弱的。氧上方的等值面的內層部分主體都是綠色，表現出水分子以氧作為受體與其它水分子構成的氫鍵還算比較穩定。圖中最上方的一大片紅色等值面表現出這部分弱相互作用很難保持得住，波動性極大，因此將之忽略了也無妨（aRDG原文中干脆直接把這個區域給截掉了，以避免說不清楚）。
最后，我們來看看如果我們只考慮動力學軌跡中的一幀得到的結果會是如何。下面的圖顯示的是隨機選取的一幀得到的RDG等值面填色圖
可見，如果只考慮一幀，表現出的弱相互作用和平均的弱相互作用有很大差距。圖中只顯示了兩個氫鍵區域，和零星的范德華作用區域，從這樣的單幀的圖上我們顯然是沒法直接了解到整個動力學過程中這個水分子是怎么和其它水分子作用的。
5 總結&其它
aRDG方法給研究動態體系中的分子間相互作用增加了一個很重要的新途徑。以往研究這類問題往往就是計算徑向分布函數、計算原子間距離/角度數值分布之類，是間接表征手段，遠沒有aRDG方法直觀、表現的信息豐富。aRDG不僅可以研究動力學過程中小分子與其它小分子的平均弱相互作用，還可以研究諸如受體-蛋白間、小肽-小肽之間、分子與固體表面間的平均弱相互作用。aRDG原文中還研究了過渡態結構下體系與溶劑的平均弱相互作用。讀者別忘了閱讀《使用Multiwfn做aNCI分析圖形化考察動態過程中的蛋白-配體間的相互作用》（http://www.shanxitv.org/591），里面有Multiwfn做aRDG分析更多的信息。
相對于單一結構，動力學過程中包含了很多可以影響弱相互作用的因素，特別是溫度。比如研究不同溫度下平均RDG、平均sign(lambda2)rho以及熱穩定指數的差異，對于研究溫度對弱相互作用的影響十分有益。其它因素，諸如電場、濃度、壓力等等因素對平均弱相互作用的影響也可以嘗試通過aRDG分析來研究討論。
有人問做aRDG分析需要跑多少ns的動力學才夠用，這種問法很不妥，提供的信息不充分。做aRDG分析的關鍵在于：(1)要讓被統計的粒子采樣充分，避免由于模擬時間不夠長導致粒子本來會出現的地方在有限的模擬時間內沒有出現 (2)要有足夠多的幀數描述粒子的分布，幀數不夠的話得到的aRDG等值面會非常粗糙，棱角、窟窿明顯。滿足(1)需要軌跡足夠長，多長合適應當根據實際體系運動特征判斷，諸如像本文這樣的體系，被統計的分子在中心分子附近運動得非常自由，而且被統計的分子是大量的（環境分子全都是被統計的分子），因此跑1 ns這樣不算長的時間就能采樣比較充分。滿足(2)需要幀數夠多，幀數需要多少不光在于跑的動力學總時間長度，還在于軌跡保存頻率，顯然同樣長模擬時間內軌跡保存頻率越高幀數就越多。另外，如果被統計的粒子運動范圍比較廣，也需要越多的幀數，以使得粒子在各個位置的出現都有足夠多的幀描述；反之如果粒子分布范圍非常有限，那么很少的幀數就足矣采樣充分（極限情況下粒子只出現在一個位置，此時只需要一幀就夠了，故只需要做RDG而沒必要做aRDG了）。如果你不知道幀數夠不夠，姑且先跑跑一段動力學做aRDG試試，如果等值面還不理想，再續跑來得到更多幀。
最后，我把筆者講授的分子動力學與GROMACS培訓班（http://www.keinsci.com/workshop/KGMX_content.html）里示例的苯酚在水環境中的aRDG分析也展示出來，大家用GROMACS可以很容易跑出來。
