使用Sobtop超級方便地創建二茂鐵的GROMACS的拓撲文件
使用Sobtop超級方便地創建二茂鐵的GROMACS的拓撲文件
文/Sobereva@北京科音 2022-Feb-17
0 前言
Sobtop是筆者開發的主要用于產生各種類型體系的GROMACS拓撲文件的程序，可以在主頁http://www.shanxitv.org/soft/Sobtop免費下載。過渡金屬配合物體系的拓撲文件在以往是很難產生的一類情況，因為里面牽扯到GAFF、OPLS-AA、GROMOS等主流有機分子力場不支持的元素，里面涉及到的鍵/鍵角/二面角參數在那些力場里也沒有，所以acpype、Ligpargen、ATB等程序根本沒法搞。而這類體系，用Sobtop則可以超級容易地產生拓撲文件。AmberTools里的MCPB.py程序雖然能產生這類體系的拓撲文件，但是過程麻煩至極，在http://bbs.keinsci.com/thread-27796-1-1.html里有人演示了怎么借助MCPB.py產生二茂鐵的GROMACS的拓撲文件，估計大多數人看了后會打退堂鼓。Sobtop的詳情以后會另文專門系統性地介紹，本文僅簡單演示一下怎么用Sobtop快捷地產生二茂鐵的拓撲文件，通過對比大家可以充分認識到用Sobtop有多么方便。
本文牽扯到的所有輸入輸出文件都可以在http://www.shanxitv.org/attach/635/file.rar里找到。Sobtop用的是2022-Feb-16發布的預覽版1.0(dev)。Multiwfn用的是2022-Feb-17更新的版本，Multiwfn可以在主頁http://www.shanxitv.org/multiwfn免費下載。
1 用Gaussian做優化和振動分析
首先用Gaussian對二茂鐵做優化和振動分析。要注意，二茂鐵的兩個茂環間僅在低溫晶體中是交錯式的，而在真空中極小點結構是重疊式的（參看Materials, 8, 7723 (2015)和https://en.wikipedia.org/wiki/Ferrocene#Determining_the_structure中的相關信息），因此本例用重疊式（D5h點群）的初始結構。記得Gaussian計算前最好用GaussView做對稱化，這樣優化得又精準速度又快。
關于計算級別，本例用的TPSSh/def2-TZVP級別對于過渡金屬配合物體系通常是很理想的選擇。關于泛函的選擇看《簡談量子化學計算中DFT泛函的選擇》（http://www.shanxitv.org/272）。如果想明顯更便宜，Fe可以用SDD，配體可以用6-311G*。
此例的Gaussian輸入文件如下
%chk=Ferrocene.chk
#p tpssh/def2tzvp opt freq
[空行]
Generated by Multiwfn
[空行]
0 1
C 0.71056100 -0.97800331 1.65426528
C -0.00000000 1.20887857 1.65426528
C -0.71056100 -0.97800331 1.65426528
H 2.17848562 0.70783289 1.62195576
H -2.17848562 0.70783289 1.62195576
Fe 0.00000000 0.00000000 -0.00000000
C 1.14971184 0.37356402 1.65426528
C -1.14971184 0.37356402 1.65426528
H 1.34637816 -1.85313055 1.62195576
H -0.00000000 2.29059533 1.62195576
H -1.34637816 -1.85313055 1.62195576
C 1.14971184 0.37356402 -1.65426528
C -1.14971184 0.37356402 -1.65426528
C 0.71056100 -0.97800331 -1.65426528
C 0.00000000 1.20887857 -1.65426528
C -0.71056100 -0.97800331 -1.65426528
H 2.17848562 0.70783289 -1.62195576
H -2.17848562 0.70783289 -1.62195576
H 1.34637816 -1.85313055 -1.62195576
H 0.00000000 2.29059533 -1.62195576
H -1.34637816 -1.85313055 -1.62195576
[空行]
[空行]
算完了之后，用formchk把當前目錄下的Ferrocene.chk轉換為Ferrocene.fchk。不知道formchk是什么和怎么用的話看《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）里相關部分。
還需要得到二茂鐵的mol2文件，并且要求其中記錄的鍵連關系和期望的一致，這決定了Sobtop給出的拓撲文件里都有哪些力場項。很多程序都可以實現這個，此例我們用常用的GaussView打開Ferrocene.fchk，如果看到Fe和各個C之間沒顯示鍵的話就手動連上（設幾重鍵都無所謂，Sobtop只看有沒有鍵）。然后保存為Ferrocene.mol2。
2 計算RESP電荷
RESP電荷很適合用于動力學模擬的目的，我們要讓拓撲文件里的原子電荷對應RESP電荷。計算RESP電荷最強大、方便、靈活的程序是Multiwfn，見《RESP擬合靜電勢電荷的原理以及在Multiwfn中的計算》（http://www.shanxitv.org/441）。要注意對二茂鐵體系不宜直接按照常規方式算RESP電荷，因為擬合點的分布不會恰好滿足D5h對稱性，因此得到的等價原子的原子電荷會有微小的偏差。為了讓所有等價的原子的電荷都嚴格相同，在Multiwfn中計算RESP電荷時應當設置等價性。
啟動Multiwfn，然后輸入
Ferrocene.fchk
7 //計算原子電荷
18 //計算RESP電荷
5 //等價性約束設置
11 //根據點群自動判斷等價性
a //根據整個體系的結構判斷點群
y //當前屏幕上顯示的點群確實是D5h，列出的原子等價性沒錯，所以這里輸入y確認。然后當前目錄下有了等價性設置文件eqvcons_PG.txt
q //返回
1 //從文本文件里讀取等價性設置
eqvcons_PG.txt //剛才產生的文件
2 //開始一步式RESP電荷擬合
[回車] //程序提示Fe缺半徑，按回車表示用推薦的半徑
Multiwfn算RESP電荷效率很高，RESP電荷很快就產生了，從屏幕上給出的結果看確實原子電荷合理，而且滿足等價性。按y把產生的原子電荷導出為當前目錄下的Ferrocene.chg。
3 用Sobtop產生拓撲文件
啟動Sobtop，然后依次輸入
Ferrocene.mol2
7 //指定原子電荷
10 //從Multiwfn的chg文件中讀取原子電荷
Ferrocene.chg
0 //返回
2 //產生GROMACS的.gro文件
[回車] //此時當前目錄下產生了Ferrocene.gro
1 //產生GROMACS的拓撲文件
2 //優先用GAFF原子類型，缺的自動用UFF力場的。原子類型決定了拓撲文件里的LJ參數
2 //所有成鍵相關參數基于振動分析產生的fchk文件里的Hessian矩陣生成
Ferrocene.fchk //如果此文件和Ferrocene.mol2在相同目錄下，此處直接按回車可以快捷載入
[回車] //在當前目錄下產生Ferrocene.top
[回車] //在當前目錄下產生Ferrocene.itp
現在拓撲文件在當前目錄下出現了，Sobtop可以關了。現在可以打開itp文件看看內容，會看到沒有任何問題，文件整齊干凈，注釋也特別清楚。
4 測試拓撲文件合理性
產生拓撲文件之后，我一般建議在真空下對這個體系跑一下動力學，根據模擬過程中結構變化情況大致判斷拓撲文件是否合理。做這個模擬要用的Ferrocene.gro、Ferrocene.top、Ferrocene.itp前面我們都已經產生了。做這個模擬用的md.mdp在本文的文件包里也給了，可見是控溫在298.15K下模擬100 ps，1 fs一步，每100 fs往xtc里寫入一幀，沒有用PBC。筆者此例用的是GROMACS 2018.8。這個mdp不兼容GROMACS較新版本，較新版本用戶需要在一個足夠大的空盒子里在考慮PBC的情況下跑NVT模擬。
把前述的gro、itp、top、mdp都放到當前目錄下，運行以下命令
gmx grompp -f md.mdp -c Ferrocene.gro -p Ferrocene.top -o md.tpr
gmx mdrun -v -deffnm md
然后用VMD載入軌跡，用Extensions - Analysis - RMSD Trajectory Tool工具對system選區做align消除平動轉動后，每10幀疊加顯示一次，得到下圖
可見模擬過程中二茂鐵很好地維持了剛性，說明拓撲文件合理。
5 總結&其它
本例充分體現了使用Sobtop結合Multiwfn程序產生過渡金屬配合物體系拓撲文件的便利，把復雜度降低到了極限。此文只以極為簡單的體系做了示例，對更復雜的情況，諸如多核金屬配合物、除了剛性部分還帶有可旋轉鍵的柔性部分的配合物，Sobtop也都能很容易地產生拓撲文件。凡是體系全為剛性的情況都可以直接效仿本文的做法；對于既有剛性配位區域也有柔性基團的配合物，在Sobtop問你怎么產生成成鍵相關(bonded)參數時，建議選擇相應選項，讓bond和angle參數用特定方法基于Hessian產生，而其它的用預置力場參數，這樣可以使得GAFF的參數用于描述柔性部分二面角可旋轉特征。
實際中二茂鐵的茂環的旋轉勢壘是非常低的，常溫下是可以旋轉的，但本文得到拓撲文件對應的是純剛性二茂鐵。本文暫不考慮茂環間可相對旋轉的問題，這也不是什么重要問題。
本文的做法比起MCPB.py實在方便百倍，而且MCPB.py用的Seminario方法產生力常數也不如本文的方法合理（當前版本Sobtop默認用的是modified Seminario，產生的鍵角力常數遠好于Seminario）。
做優化和振動分析用免費的ORCA亦可，文中用fchk的地方改成用ORCA振動分析產生的.hess文件即可。
對于體系沒有對稱性的情況，算RESP電荷時不需要本文這樣做對稱等價設置，直接在Multiwfn的RESP界面里選擇1用常規的兩步式擬合即可。更多細節看《RESP擬合靜電勢電荷的原理以及在Multiwfn中的計算》（http://www.shanxitv.org/441）。
本文的例子直接從chg文件中載入了RESP電荷。也可以在Sobtop里先不設置電荷，等itp產生出來之后再手動把要用的電荷寫入到[ atoms ]部分。
使用Sobtop產生拓撲文件發表文章時記得按照Sobtop主頁http://www.shanxitv.org/soft/Sobtop的說明進行引用。引用方式在未來會有所變化，請記得發文章之前看一下主頁上最新的引用說明。