Sobtop:A tool of generating forcefield parameters and GROMACS topology file
Latest version:1.0(dev3.1) Last update: 2022-Aug-9
Index of this page:
Developer Citation Download Introduction Examples Skill FAQ TODO list Update history
Developer
Dr. Tian Lu (Contact: sobereva[at]sina.com. Beijing Kein Research Center for Natural Sciences 北京科音自然科學研究中心)
If you encountered any difficulty in using Sobtop, or you found bug, or you have any suggestion on improving Sobtop, please post on Molecular Modeling board of http://bbs.keinsci.com or send E-mail to me.
Sobtop can be freely used for both academic and commerical purposes. Sobtop must be properly cited if the work will be published. You may also use Sobtop to do calculations for others (i.e. 代算), while you should explicitly tell them that Sobtop needs to be cited in their publication.
Citation
If Sobtop is utilized in your work, it must be cited as follows in main text (rather than SI):
Tian Lu, Sobtop, Version [當前版本], http://www.shanxitv.org/soft/Sobtop (accessed on 日 月 年)
Sobtop在進一步完善后會寫詳細的手冊并發布正式版,然后專門寫成介紹paper發在學術期刊上,屆時請引用期刊上的文章(記得提交proof前關注本頁面看這里寫的最新的引用方式)。
Download
Executable file: sobtop_1.0(dev3.1).zip ("sobtop.exe" and "sobtop" are executable files of Sobtop in Windows and Linux, respectively)
Introduction of Sobtop
What is Sobtop?
GROMACS is the most popular molecular dynamics code. To simulate a molecule by GROMACS, topology file must be created first for every kind of species involved in the simulation system. Although there have been numerous codes that can generate topology file of GROMACS, all of them have severe limitations in different aspects, for example: (1) Difficult to use and install (2) Only support specific kind of system, such as organic molecules only, while molecules containing transition metals and boron can hardly be treated (3) Cannot automatically determine missing bonded parameters based on Hessian matrix (4) Too slow and unable to deal with large systems such as polymer (5) Only support isolated system (6) Not flexible enough to create topology for complicated systems.
Sobtop aims at eliminating all limitations in existing topology file generators, providing an extremely easy-to-use, universal and robust code to generate forcefield parameters and GROMACS topology file for any system!!! I sincerely hope Sobtop could be a perfect and final solution of generation of GROMACS topology file.
(Meaning of "Sobtop": Sobereva topology)
Features of Sobtop
- Freely available.
- Very easy to use, the steps for generating GROMACS topology file are minimized as much as possible. Preparing any running environment is not needed (in contrast to many Python scripts. Sobtop is written in Fortran).
- The program is interactive, users do not need to memorize any keyword or argument. User can also write shell script to invoke Sobtop to generate topology file for a batch of systems.
- Very fast. Even macromolecules containing thousands of atoms can be processed.
- Can deal with any kind of system (organic molecules, inorganic molecules, transition metal coordinates, covalent atomic clusters, polymers, inorganic materials, etc.) containing any element.
- Both isolated and periodic systems are supported. The x2top utility in GROMACS can be completely replaced.
- GAFF, AMBER and UFF atomic types can be automatically determined (for the former two cases, .mol2 file should be used, in other cases .pdb file can also be used). Users can also provide a file containing rules of determination of atomic types according to chemical environment.
- Prebuilt forcefield parameters and definition of atomic types are completely open to users, thus they can very conveniently modify and extent the parameters by editing force field files.
- Force constants of bond, angle and dihedrals can be directly derived based on Cartesian Hessian matrix calculated by Gaussian, ORCA and xtb (for isolated systems) or CP2K (for periodic systems). Four methods are supported for this purpose: Seminario, mSeminario, m2Seminario, DRIH. These methods can be automatically employed during generating .itp file, and minimum, maximum and average of each kind of derived bond/angle/dihredral types are outputted to a text file so that users can then easily put them into forcefield parameter file. In addition, a specific function is also provided in Sobtop to allow users manually calculate bond/angle/dihredral parameters by directly inputting atomic indices.
- Bonded parameters can be assigned in flexible ways. For example, one can request Sobtop to derive parameters for bonds and angles based on Hessian (thus more accurate than GAFF parameters), while employ torsion parameters of GAFF to account for rotatability of dihedrals. Another example, one can specify a rigid region, all bonded terms involved in this part are represented by harmonic potential, whose parameters are derived based on Hessian, while other parts employ GAFF parameters and thus dihedral flexibility can be represented.
- Atomic charges can be directly loaded from .chg file exported by atomic charge calculation (e.g. RESP, ADCH) of Multiwfn. EEM and Gasteiger atomic charges (via invoking Multiwfn) and MMFF94 atomic charges (by invoking OpenBabel) can also be automatically assigned.
- Very detailed comments are given in the generated topology file, which greatly facilitates users to check and manually modify parameters.
- .gro file can be exported, which exactly corresponding to the generated .itp file.
- .rtp file can be generated, which can be used by pdb2gmx tool.
- Some utility functions: Printing current redundant internal coordinates, exporting Wilson B matrix and related matrices, generating and exporting Hessian matrix in redundant internal coordinate, and so on.
Examples
全面包含各類體系例子的手冊目前還沒時間寫。這里僅給一些簡單的例子,以及一些要點的說明,大家請舉一反三。基本上看過一遍后就能對Sobtop的各方面主要特征有個全面的了解。強烈建議按順序閱讀,因為這些例子和相關講解是循序漸進的,也絕對不要斷章取義。使用Sobtop的過程中要注意看屏幕上的提示,信息非常豐富、細致、易懂。
例子列表:
例1:產生SF6的GROMACS的拓撲和gro文件:力常數全基于Hessian計算得到
例2:產生苯甲酸甲酯的GROMACS的拓撲和gro文件:參數完全基于GAFF的
例3:產生苯甲酸甲酯的拓撲和gro文件:bond和angle參數基于Hessian得到,其它部分用GAFF的
例4:產生過渡金屬配合物的拓撲文件
例5:產生聚合度為40的氯丁二烯聚合物的拓撲文件
例6:產生跨越盒子的無限長氯丁二烯聚合物的拓撲文件
例7:產生碳納米管的GROMACS拓撲文件
例8:產生金剛石晶體的拓撲文件
例9:產生二氧化硅晶體的拓撲文件
例10:產生金屬有機框架化合物MOF-5的拓撲文件
例11:產生FeCl2二維材料的拓撲文件
例1:產生SF6的GROMACS的拓撲和gro文件:力常數全基于Hessian計算得到
注:絕對不要以為什么時候都得像本例這樣提供含有Hessian的文件并讓Sobtop計算力常數。如果你要對普通有機體系產生拓撲文件,應該效仿的是例2直接用GAFF的參數,那是最省事的。老有人只看了本例,居然連例2都不看就去胡搞。也記得把FAQ10和FAQ11看了,搞清楚什么時候才有必要基于Hessian算力常數。
SF6是典型的非有機體系,那些基于有機分子力場的拓撲文件產生工具如acpype、antechamber、Ligpargen、CGenFF、ATB等碰到這種體系就完全無能為力了。碰到這種體系一般都建議讓Sobtop基于Hessian矩陣計算出來力場參數。對于[H3O]+、Sb4O6、P4、[PF6]-、[AuCl4]-、B2H6等等亂七八糟的無機剛性體系一律都可以用本節的做法非常容易地得到拓撲文件。
首先用Gaussian對SF6做優化+振動分析并保留.fch/fchk文件,或使用ORCA做優化+振動分析并保留.hess文件,里面都包含了要用的笛卡爾Hessian矩陣和優化后的坐標信息。然后把優化后的結構保存成mol2文件,這用GaussView就可以。注意此文件里的鍵連關系必須和實際要產生的力場項對應,在GaussView里可以肉眼檢查,不符的話就調整。只要連接關系正確就行,而成鍵方式隨意。
在Windows下可雙擊Sobtop.exe圖標啟動。在Linux下可進入Sobtop目錄后運行chmod +x *加上可執行權限,再輸入./sobtop啟動(不要在其它目錄下啟動,否則無法正常運行)。啟動Sobtop后依次輸入以下命令,//后是注釋
examples\SF6\SF6.mol2
2 //產生gro文件
[回車] //默認的gro文件產生路徑
-1 //設置產生力常數的方法
4 //DRIH方法
1 //產生GROMACS拓撲文件
2 //指認GAFF原子類型,沒法指認的自動用UFF原子類型
2 //所有力常數都通過DRIH方法得到
examples\SF6\SF6.fchk //Gaussian的freq任務得到的,用于給Sobtop提供Hessian矩陣。這一步用ORCA的freq任務產生的.hess文件、xtb的--ohess任務產生的名為hessian的文件、CP2K振動分析的輸出文件亦可
[回車] //用默認的top文件產生路徑
[回車] //用默認的itp文件產生路徑
當前目錄下已經有了SF6.itp、SF6.top和SF6.gro。此例原子類型都是GAFF的,平衡鍵長來自SF6.fchk,力常數使用DRIH方法基于Hessian產生。SF6.itp里包含的是SF6分子的拓撲信息,包括定義此分子的[ moleculetype ]字段,以及定義其中涉及到的原子類型的[ atomtypes ]字段。由[ moleculetype ]下面的信息可見Sobtop自動設的分子名就是結構文件去除后綴后的名字。SF6.top是主top文件,包含了GAFF/AMBER力場下模擬適用的[ defaults ]字段、整個體系名[ system ]字段,以及分子數目字段[ molecules ]。
如果mol2文件里記錄了原子電荷的話,它們會被讀取并寫入到產生的itp里。由于當前mol2文件里沒有記錄原子電荷信息,所以所有原子電荷都為0。因此要自行用Multiwfn算RESP或RESP2(0.5)原子電荷添加進itp中的[ atoms ]部分,溶液中的模擬用RESP2(0.5)最為推薦。具體做法參考下文,可以手動用Multiwfn基于Gaussian產生的fch文件或ORCA產生的molden文件算,也可以直接用下面的懶人腳本一行命令便利地得到。
RESP擬合靜電勢電荷的原理以及在Multiwfn中的計算
http://www.shanxitv.org/441(http://bbs.keinsci.com/thread-10880-1-1.html)
計算RESP原子電荷的超級懶人腳本(一行命令就算出結果)
http://www.shanxitv.org/476(http://bbs.keinsci.com/thread-12858-1-1.html)
RESP2原子電荷的思想以及在Multiwfn中的計算(內含懶人腳本)
http://www.shanxitv.org/531(http://bbs.keinsci.com/thread-16190-1-1.html)
ORCA結合Multiwfn計算RESP、RESP2和1.2*CM5原子電荷的懶人腳本
http://www.shanxitv.org/637(http://bbs.keinsci.com/thread-28178-1-1.html)
順帶一提,也可以先用Multiwfn計算出包含原子電荷的chg文件,在產生拓撲文件之前讀進Sobtop,這樣chg文件里的原子電荷直接就出現在產生的itp的[ atoms ]里。具體做法是啟動Sobtop后,在主菜單里選7 Assign atomic charges,再選10 Load atomic charges from .chg file of Multiwfn,輸入.chg文件的路徑,然后回到主菜單,按照前面的例子操作即可。圖省事的話,也可以在Sobtop主菜單里直接輸入.chg文件的路徑,也會從里面載入原子電荷。
練習:examples\exercise\COBH3.mol2是用Gaussian在B3LYP/def2-TZVP級別下優化好的COBH3分子對應的mol2文件,請讀者對它效仿本節的步驟產生拓撲文件。
小知識:如何檢驗拓撲文件的合理性
如果想檢驗拓撲文件是否合理、其中的參數是否合適,可以寫個對應真空、常溫條件下短時間動力學的mdp,與前文得到的itp、top、gro相結合跑模擬。通過在VMD里看軌跡動畫,或者每隔一定幀數疊加顯示,靠肉眼判斷結構變化和波動是否合理。通常憑直覺就可以大致判斷合理性,如果有基于DFT的從頭算動力學的軌跡作參照那就更容易檢驗,可以對比基于力場跑的軌跡的RMSF、角度分布等參數與從頭算動力學跑的偏差。
在examples目錄下給了個md_vac.mdp,是控溫在300 K、用步長1 fs跑真空下100 ps動力學的mdp。沒有用PBC條件,只兼容GROMACS 2018及以前的版本。可以用以下命令進行動力學模擬
gmx grompp -f md_vac.mdp -c SF6.gro -p SF6.top -o md.tpr -maxwarn 1
gmx mdrun -v -deffnm md -ntmpi 1
也要注意并不是真空下模擬下看起來合理,在溶劑環境下模擬也就一定合理。因為溶劑環境下涉及到分子與溶劑的非鍵相互作用,以及溶劑與分子的碰撞,情況明顯更為復雜。如果構建的拓撲文件的實際目的是用于溶劑下的模擬,顯然還是要以溶劑中模擬的軌跡作為合理性的最終判據。
文獻里衡量力場參數準確性的一種參見做法是基于力場先做充分的能量極小化再做振動分析,和量子化學算的頻率對比看偏離程度。GROMACS的雙精度版本是可以做振動分析的,雖然有時候有bug。雖然這樣檢驗力場參數確實顯得比較嚴格,但其實對于實際模擬來說我覺得用處不大。因為對剛性體系只要結構能在動力學過程中維持好就行了,頻率相符程度差一些也無所謂,只不過是結構波動情況可能沒那么精準而已,但這不是一般關心的。而對于涉及到可旋轉鍵的柔性分子,就算極小點位置的頻率算得不錯,也不代表動力學行為就一定很理想,因為實際的動力學牽扯到勢能面上較廣闊區域,而頻率的計算精度僅體現極小點位置的局部勢能面曲率的合理性而已。
小知識:基于Hessian產生力常數的方法
Sobtop支持多種基于Hessian產生力常數(k)的方法,對bond、angle、dihedral都可以產生力常數。顯然這種情況下的二面角是視為剛性的而不可旋轉,此時Sobtop也不會試圖產生維持pi共軛區域平面性的improper項。默認方法可以通過settings.ini里的k_method設置。
Sobtop支持的各種基于Hessian產生力常數的方法簡要說明如下:
? Seminario方法:原文是Int. J. Quantum Chem.: Quantum Chemistry Symposium, 30, 1271 (1996)。已經過時了,算的鍵角力常數嚴重高估。總是不建議使用。不過鍵角力常數高估其實對實際MD模擬來說不算大問題,只不過是高估了鍵角的剛性而已。
? Modified Seminario (mSeminario)方法:原文是J. Chem. Theory Comput., 14, 274 (2018)。是對Seminario方法產生鍵角力常數的改進,解決了高估的問題。但mSeminario和Seminario方法在原理上不嚴格,導致的一個問題是對于某些體系對稱的鍵角給出來的力常數可能明顯不同,可能因此跑MD的時候分子對稱性被破壞,比如對CLi4、SF6就都有這個問題。用戶可自行檢查一下itp里[ angles ]部分等價鍵角的力常數值。另外,此方法產生的i-j和j-i鍵的力常數可能有輕微不同,即序號順序會影響結果,這也是違背邏輯的,但由于此問題很輕微所以不用太在意。
? Diagonal elements of redundant internal coordinate Hessian (DRIH)方法:先將讀入的笛卡爾Hessian矩陣通過Wilson B矩陣轉化成冗余內坐標下Hessian矩陣,其中每個冗余內坐標直接對應一個bonded項,然后取對應的對角元當做力常數。此方法原理上比mSeminario更嚴格。對于有高對稱的剛性小體系如CLi4、Sb4O7、SF6,用DRIH多數情況比mSeminario更好(從基于力場計算的頻率相對于量子化學算的頻率的吻合來說),也完全沒有其空間對稱等價的鍵角力常數不相同(實際中由于數值精度原因,可能還是有微小不同,可以無視)、i-j和j-i鍵力常數輕微不同的問題。另外用DRIH方法時itp文件里會有以[ angles ]字段記錄的鍵-鍵耦合項,因此對鍵伸縮描述精度比mSeminario好得多,雖然這對于實際模擬來說不是重點。但對有些體系,比如苯、二茂鐵,用DRIH時我發現二面角力常數偏小,導致模擬過程中結構波動偏大。而且DRIH結合GAFF力場的二面角項使用時,對量子化學算的中、低頻率的重現性往往比mSeminario結合GAFF的二面角項時差一些。注意:使用DRIH方法產生力常數時,如果體系實際是有點群對稱性的,一定要讓初猜結構就滿足點群對稱性,通常這可以令最終優化出的結構也滿足點群對稱性。如果偏離實際對稱性,有可能令DRIH方法給出的力常數不合理。
? The second modified Seminario (m2Seminario)方法:我提出。改進了mSeminario方法中用的不太合理的Hessian矩陣非對角塊的投影方法,使得等價的鍵角的力常數完全相同,而且i-j和j-i的力常數也完全相同。雖然此方法比mSeminario原理上更嚴格,但對于個別體系的個別力常數,mSeminario方法由于內在的“錯錯得對”效應,算的力常數反倒比此方法更好。例如此方法算的水的鍵角力常數明顯偏小,而mSeminario算的則更合理。另外,結合GAFF力場的二面角項使用時,此方法對量子化學算的頻率的重現性往往稍微比mSeminario差一些(但也有反例)。
總的來說,沒有絕對完美的方法。一般對于純剛性,尤其是有對稱性的體系,我建議先考慮DRIH或m2Seminario。對于含有可旋轉鍵的有機分子,建議mSeminario計算鍵+鍵角,結合GAFF的二面角項(會順帶產生pair和improper項)。如果你不怎么講究,有機分子全用內置的GAFF參數也可以,如后面的苯甲酸甲酯的例子,但如果體系有成鍵方式比較特殊、在GAFF力場中沒充分考慮的部分,可能此時很不理想。
如果你實在不確定哪種最好,可以都試試(Seminario方法除外),然后按照前述的“檢驗拓撲文件的合理性”里的做法予以檢驗。
值得一提的是,力常數越大,動力學模擬過程中相應變量的波動程度就會越小。Sobtop的主界面上有個-2 Set scale factor for derived force constants選項,可以給基于Hessian算出來的不同類型力常數乘上特定的系數。比如模擬中你嫌分子骨架整體波動有點太厲害,就可以用這個選項把二面角力常數要乘的系數設成一個明顯大于1的值,然后再模擬試試。
例2:產生苯甲酸甲酯的GROMACS的拓撲和gro文件:參數完全基于GAFF的
本節的做法是產生普通有機分子拓撲文件最簡單快速的做法,沒有之一。
啟動Sobtop,然后輸入
examples\Methyl_benzoate\Methyl_benzoate.mol2
2 //產生gro文件
[回車] //默認的gro文件產生路徑
1 //產生GROMACS拓撲文件
2 //指認GAFF原子類型,沒法指認的自動用UFF原子類型
4 //從力場庫文件中取合適的參數,缺的就猜一個(此模式對于bond和angle項,平衡值設為當前結構文件里的值,力常數分別用差不多是c3-c3和c3-c3-c3的值,但取了個盡量整的數。對于二面角項則將旋轉勢壘設為0)
[回車] //用默認的top文件產生路徑
[回車] //用默認的itp文件產生路徑
當前目錄下已經有了Methyl_benzoate.itp、Methyl_benzoate.top和Methyl_benzoate.gro。其中LJ參數和成鍵相關參數都是GAFF力場的。在屏幕上沒有提示有任何bond、angle、dihedral參數缺失。雖然提示ca-ca-ca-c的improper參數缺失,但Sobtop自動用X -X -ca-ha的improper參數來代替的做法通常都是合理的,足夠維持局部平面性。
最后,記得用Multiwfn算出來RESP或RESP2(0.5)原子電荷添加進去,在上一節明確說了。
練習:examples\exercise\7-helicene.mol2是7-螺烯的mol2文件,大家對此分子也按本例的流程產生拓撲文件,然后在真空下跑跑動力學,看看結構維持得怎么樣。
小知識:關于缺參數的處理
至于怎么搞缺的參數,看具體缺什么:
如果是缺可旋轉二面角參數,最好是構建個簡單的模型體系(關注的二面角所處的化學環境要和實際分子中一致,無關區域應盡可能簡化以避免礙事和浪費計算時間),用量子化學程序做二面角的柔性掃描(看比如《詳談使用Gaussian做勢能面掃描》http://www.shanxitv.org/474),然后靠Origin、自己寫的程序或者現成工具(如http://rotprof.lncc.br)擬合成周期勢。顯然這個過程略麻煩,如果這個二面角不重要的話,也可以自己去bonded_param.dat的$DIHEDRAL里憑感覺找個能湊合用的參數。另外,還可以用Google(never baidu)搜A-B-C-D parameter,其中A/B/C/D替換為實際缺的原子類型名,說不定有的文章里提供了參數。如果是缺剛性二面角參數,讓Sobtop基于Hessian直接算出來就完了,因此上面選4的地方應該選3,此時Sobtop會要求你輸入.fch/fchk或.hess的路徑,產生itp文件過程中有缺的bonded參數時就會自動算。
如果缺的是bond和angle參數,一般都應當讓Sobtop基于Hessian算出來。如果你希望圖省事,不怎么要求準確度,那就干脆直接用Sobtop給你猜的,平衡鍵長/鍵角就是當前結構文件里的值,而力常數是一個足夠大的值。
小知識:關于指認原子類型
在指認原子類型的那一步里,從屏幕上可見Sobtop支持的指認原子類型方式相當靈活豐富。除前述例子用到的外,還有其它選項:自動指認AMBER原子類型、給沒指認的原子類型自動用UFF的、手動設置原子類型、根據assign_AT.dat里的規則判斷原子類型、從外部文件讀取各個原子的原子類型、導出當前的原子類型等。在這一步里,屏幕上會顯示當前所有原子的原子類型,用戶需要用界面上的選項把所有未定義類型的原子的類型都定義了,才能選0進入下一步。而前面的例子用的選項2算是快捷選項(先指認GAFF原子類型,沒指認成功的自動用UFF原子類型),選完了直接就進入了下一步。
小知識:Sobtop的力場庫文件
由于Sobtop的力場庫文件完全透明、格式清晰、注釋詳細,因此Sobtop極為靈活、普適,相當于一個通用的框架。由于用戶可以自己修改力場庫文件,故Sobtop能自動提供的參數絕不僅限于GAFF和AMBER力場的。不過,Sobtop目前只支持GAFF/AMBER的力場形式,其它力場如果用了特殊形式,就不可能被支持,例如CHARMM力場有Urey-Bradley項因此不支持,GROMOS力場的bond和angle項是特殊形式因此不支持。
在Sobtop目錄下有LJ_param.dat文件,里面包含所有原子類型的LJ參數和原子類型含義的注釋。itp文件里的LJ參數就是直接取自這里的。此文件里默認包含各種UFF力場、GAFF力場(原子類型為小寫)和AMBER19SB+parmbsc0+OL3力場(原子類型為大寫)的原子類型。用戶可以自己用文本編輯器編輯此文件,修改參數或添加新原子類型,詳見文件中的注釋部分。添加新原子類型時原子類型名不要超過10個字符,注意區分大小寫,也不要和已有的沖突。
在Sobtop目錄下有bonded_param.dat文件,里面包含了所有預置(prebuilt)的bond、angle、dihedral、improper參數(統稱為bonded參數),分別在不同字段里記錄,各字段以$開頭,包括$BOND、$ANGLE、$DIHEDRAL、$HARM_DIH、$IMPROPER。Sobtop會根據被指認的原子類型和成鍵關系,在此文件里搜索對應的參數。搜索是從前往后搜索,取最后一次匹配的,所以二面角參數中涉及到X的通配項一定要出現在專一項之前,以確保最終用的盡可能是專一項。此文件里的注釋極其詳細,用戶可以很方便地自行修改里面的值、添加新的值。此文件里默認包含了所有GAFF和AMBER19SB+parmbsc0+OL3力場定義的bonded參數。
應知道二面角可以分為兩類,一種是可旋轉(rotatable)二面角或者叫柔性(flexible)二面角,用周期性的扭轉勢(torsion potential)描述;另一種是剛性二面角,用諧振勢(harmonic potential)描述。bonded_param.dat文件里$DIHEDRAL字段記錄的是前者的參數,$HARM_DIH記錄的是后者的參數。Sobtop先搜索前者后搜索后者,因此如果某個二面角參數在兩處都定義了,那么會取$HARM_DIH里的,并因此作為剛性二面角對待。
bonded_param.dat的$IMPROPER字段不僅決定了improper參數是什么,也決定了Sobtop會生成哪些improper項。此文件里A-B-C-D項對應于C在中央,A、B、D在其周圍的情況。例如此文件里有一行是X -c3-n -c3,則原子類型為n的原子若與三個原子成鍵,其中兩個是c3,一個是任意類型,就會產生一個improper項,并且用這里提供的參數。順帶一提,在GAFF/AMBER力場里,improper項在函數形式上和普通二面角項并無差異,也是周期勢函數。如果想當諧振勢對待,應當在$HARM_DIH里進行定義。
小知識:關于pair項的產生
Sobtop對于柔性二面角會產生對應的pair項,由此GROMACS會計算1-4作用,它和dihedral參數對應的扭轉勢相組合定義實際二面角扭轉時感受到的勢。那些缺參數的二面角也會產生對應的pair項,因此如果不補參數的話,這些二面角也并非完全是自由旋轉的。那些被諧振勢描述的剛性二面角就不會產生pair項了,本來也沒任何必要,而且基于Hessian計算的剛性二面角參數本身就已經把二面角改變時感受到的勢完全體現了。
注意pair項和柔性二面角并不是一一對應關系。比如環丁烷,四個碳C1,C2,C3,C4按順序排列,此體系里有C1-C2-C3-C4二面角,用GAFF參數時是當做柔性二面角對待的,雖然好像應該產生(C1,C4)的pair項,但由于C1和C4實際上是直接成鍵的,因此它們是1-2而不是1-4關系,因此得到的itp里不會對它產生pair項。
例3:產生苯甲酸甲酯的拓撲和gro文件:bond和angle參數基于Hessian得到,其它部分用GAFF的
GAFF力場里雖然定義了bond和angle參數,質量整體也不錯,但相對于基于較好質量Hessian矩陣直接計算的bond和angle參數來說還是有一定差距,尤其是對于成鍵方式特殊的體系,畢竟GAFF里的那些參數不是針對當前體系專門獲得的。而且涉及到特殊成鍵、有機體系里不常見的元素時,GAFF里可能根本沒有相應的bond和angle項。Sobtop允許bond和angle參數基于Hessian矩陣計算出來,而其它部分(dihedral和improper)從力場文件里查詢得到,其中若有缺的再基于Hessian計算得到。這樣做bond和angle項的參數質量又好又不會出現缺失項,又能兼顧柔性二面角的可旋轉特征。
啟動Sobtop,然后輸入
examples\Methyl_benzoate\Methyl_benzoate.mol2
2 //產生gro文件
[回車] //默認的gro文件產生路徑
1 //產生GROMACS拓撲文件
2 //指認GAFF原子類型,沒法指認的自動用UFF原子類型
7 //bond和angle參數基于Hessian計算得到,其它的用力場文件里的,其中若有缺失的基于Hessian計算得到
examples\Methyl_benzoate\Methyl_benzoate.fchk
[回車] //用默認的top文件產生路徑
[回車] //用默認的itp文件產生路徑
從屏幕上可以看到哪些參數是靠Hessian算出來的。現在當前目錄下已經有了Methyl_benzoate.itp、Methyl_benzoate.top和Methyl_benzoate.gro。最后,記得用Multiwfn算出來RESP或RESP2(0.5)原子電荷添加進去。
順帶一提,如果你還想要更好精度的拓撲文件,建議對可旋轉的二面角做勢能面掃描再擬合得到參數,然后將當前itp里的相應參數替換掉,但這就比較費時費事了。僅當你的研究對二面角旋轉勢壘、極小點位置準確度要求很高的時候才值得這樣做。
例4:產生過渡金屬配合物的拓撲文件
此例專門有個很詳細的文章《使用Sobtop超級方便地創建二茂鐵的GROMACS的拓撲文件》 http://www.shanxitv.org/635。對任何過渡金屬配合物都可以按照此例來搞。雙核甚至多核金屬配合物也都可以這么搞,沒有特殊性。
此例一個關鍵點是Fe這樣GAFF不支持的元素,要求Sobtop自動用UFF的原子類型,其原子類型名以UF_開頭,后面跟元素名。UFF力場里實際上每種元素不區分原子類型,而且UFF力場對周期表里幾乎所有元素都有定義,因此當一個原子遇上沒有專門的原子類型的時候都可以用UFF原子類型來湊合。Fe的UFF原子類型在Sobtop里就叫UF_Fe,可見在LJ_param.dat文件中有其LJ參數的定義。
有人覺得對過渡金屬是不是用Merz等人搞得很全面的離子的LJ參數會更好,我的意見是否定的,在于:Merz的金屬離子參數對于重現離子在溶液中的不同屬性(自由能、離子-氧距離、第一配位層的配位數)分別優化了不同的LJ參數,而且往往數值相差還不小,沒法說用哪種合適。而UFF力場每種元素就只有一種LJ參數,因此沒有這種含糊性。本來也沒證據證明在配合物中用某一種Merz的離子參數就有什么優勢。另外,過渡金屬配合物中金屬的數目相對于配體原子來說是非常少的,它對于配合物與其它分子間的范德華作用貢獻很小,因此就算LJ參數不是特別理想也不至于有明顯實際問題。順帶一提,UFF和Merz的過渡金屬參數都是離子的,如果是計算金屬單質,都不適用,必須用一些文章里提供的專門的LJ參數。另外,如果體系中過渡金屬原子數目非常多,比如FeCl2、MoS2二維材料,可能得自行專門去優化過渡金屬的LJ參數才能得到很理想的結果(UFF的LJ參數還是略糙的,而且作為普適性的參數肯定不如考慮了具體化學環境專門優化出來的參數好)。
練習:examples\exercise\ZPT.mol2是吡硫鎓鋅二聚體,里面有兩個Zn,每個Zn和兩個硫和兩個氧配位。請讀者按照本節的做法對它產生拓撲文件,并檢驗在真空下做動力學是否能正常維持結構。
Hint:搞柔性與剛性混合體系的拓撲文件
如果你的體系一部分是剛性,一部分柔性,例如過渡金屬卟啉配合物上面還接了一段柔性有機基團(顯然柔性部分必須考慮二面角的可旋轉性),在Sobtop問你怎么產生bonded參數的時候應當選擇選項5,然后輸入一批原子序號將之作為剛性部分。然后剛性部分涉及到的bonded項會自動基于Hessian計算得到,柔性部分會從力場文件里取對應的參數,柔性部分若有缺失參數的話也會基于Hessian計算得到。
更具體來說,用選項5時Sobtop是用這樣的規則:
bond項:只要有哪怕一個原子在選中的范圍,就使用基于Hessian計算的參數
angle項:只要有哪怕一個原子在選中的范圍,就使用基于Hessian計算的參數
dihedral項:只要二面角中間兩個原子都在選中的范圍,就基于Hessian計算參數。因此,如果中間兩個原子中有一個不在選中的范圍,這個二面角就是可旋轉的(如果力場文件里$DIHEDRAL里有定義的話)。
其它情況直接檢索從力場文件里取參數。
例5:產生聚合度為40的氯丁二烯聚合物的拓撲文件
這個例子產生聚合度為40,末端用氫封閉的氯丁二烯聚合物的拓撲文件,結構文件是examples\polymer\Neoprene_40.mol2。GAFF對于描述有機聚合物是完全沒問題的,本例直接指定GAFF原子類型并且用相應的bonded參數。對此體系為了圖省事,就直接用MMFF94原子電荷,它只適合有機體系,且只依賴于原子間連接關系經驗性地得到原子電荷(也因此不受輸入文件里聚合物構象的影響,沒有構象依賴性在某種程度上算是一個優點),可以瞬間對巨大體系算出來。雖然其對靜電勢重現性肯定不及RESP電荷理想,因此對于動力學模擬的目的并不算多理想,但要求不高的話對有機類體系勉強能湊合用,起碼比QEq、Gasteiger、Mulliken、NPA之類強多了,也明顯比不用原子電荷要好。如果RESP電荷是10分的話,MMFF94電荷可以給個7分左右。Sobtop可以調用免費的OpenBabel程序(http://openbabel.org/wiki/)計算MMFF94原子電荷,啟動Sobtop之前要把sobtop.ini里的OpenBabel_cmd設為OpenBabel可執行文件的實際路徑(對于Windows版,強烈不建議把OpenBabel裝到默認的帶空格的路徑下!否則可能調用不成功)。
啟動Sobtop,然后輸入
examples\polymer\Neoprene_40.mol2
7 //設置原子電荷
5 //MMFF94電荷。由屏幕上提示可見,Sobtop調用了OpenBabel,一瞬間402個原子的MMFF94電荷就產生出來了
0 //返回
2 //產生gro文件
[回車] //默認的gro文件產生路徑
1 //產生GROMACS拓撲文件
2 //指認GAFF原子類型,沒法指認的自動用UFF原子類型
4 //用預置的力場參數,缺失的自動猜
[回車] //用默認的top文件產生路徑
[回車] //用默認的itp文件產生路徑
現在當前目錄下就有了Neoprene_40.itp、Neoprene_40.top和Neoprene_40.gro了。可以嘗試做一下真空下的模擬,應當會看到聚合物鏈不斷扭動,這很合理。
從Sobtop窗口中提示可看到(在Neoprene_40.itp中搜missing也可以發現),c3-c2-cl參數是缺失的,即缺了“sp3雜化碳-sp2雜化碳-Cl”這樣的鍵角參數。當前用Sobtop自動猜的參數就跑得很理想,所以沒什么必要去補它,本來這種剛性鍵角的力常數大點小點對動力學行為都不會有什么影響。只有可旋轉的二面角參數才是真正關鍵的,因為它們直接影響分子的柔性和構象變化,而當前itp里這些參數都沒缺。另外,從Sobtop窗口中的提示還可見c3-ha-c2-c2和c2-c3-c2-cl的improper項是缺失的,sp2雜化的碳理應有improper項維持它所在的局部三角區域的平面性。對缺失的improper項Sobtop都是自動用X -X -ca-ha型improper參數湊合,這是沒什么問題的,從實際模擬結果會看到局部平面性維持得是合理的。
如果你特別講究的話,就是想把那個缺的c3-c2-cl力常數準確得到,可以構造一個三個氯丁二烯單元組成的分子,兩邊加氫,然后用Gaussian或ORCA做優化和振動分析。之后啟動Sobtop載入其mol2文件,選5 Derive rigid parameter for specific bond/angle/dihedral,再載入fch/fchk或hess文件,然后輸入對應c3-c2-cl鍵角的那三個原子的序號,Sobtop馬上就會輸出相應的平衡鍵角和力常數值。之后在bonded_param.dat的$ANGLE字段里加入c3-c2-cl項及對應的參數。下次再重復前面的步驟產生聚合物的拓撲文件時就會用這些參數了。
Hint:關于給聚合物用RESP原子電荷
如果能給聚合物所有原子都用RESP電荷那當然是最理想的。但十分不建議對整個聚合物直接算RESP電荷,一方面在于優化、算電荷的耗時會很高,尤其是原子數達到500以上很難算得動(對于2022年的主流幾十核服務器來說);另一方面是在于這樣算的原子電荷構象依賴性較大,用聚合物的不同構象來算的時候有些原子電荷可能差異不小;還有一方面在于沒法令每個中間的重復單元里等價的原子的電荷是相同的,且每個重復單元的凈電荷也不恰好為整數,這不夠優雅(做過蛋白質動力學模擬的人都知道,每種殘基里相同原子名的原子電荷都是相同的,每個殘基凈電荷為整數,聚合物也應當如此)。
一般來說,正確得到聚合物RESP電荷的做法是構造一個寡聚物(至少三個單元,五個更好,一般沒必要更多),兩邊進行飽和,并且要用平直的構象(但如果優化后卷起來了,導致其它部分的原子妨礙感興趣部分原子的RESP電荷的擬合,應當考慮做限制性優化,適當固定一些關鍵的二面角以避免卷曲)。之后用量子化學程序在DFT級別下做幾何優化,通常B3LYP/6-31G*就可以。基于優化后的結構再用B3LYP/def2-TZVPP算個單點任務得到更好質量的波函數用于給Multiwfn算RESP電荷。電荷算三次,第一次算RESP電荷的時候把首端單元的總電荷約束為0;第二次將最中間的單元的總電荷約束為0;第三次將尾端單元的總電荷約束為0。這樣三次計算后,首端、中間、尾端單元的原子電荷就都有了。上述過程的具體操作可以效仿《RESP擬合靜電勢電荷的原理以及在Multiwfn中的計算》(http://www.shanxitv.org/441)里獲得氨基酸殘基原子電荷的例子。之后,對實際的聚合物,首端和尾端單元的原子電荷套用前面得到的,中間每個單元都套用前面得到的中間單元的。至于怎么令這個套用原子電荷的過程盡可能方便,請大家結合實際情況開動腦筋,不要懶得去寫腳本(筆者以后會提供個完美、便利的解決方案)。顯然,整個聚合物里各單元中的原子排列順序必須得和寡聚物模型中相一致,否則就亂套了。
上面說的得到聚合物單元RESP電荷的做法是對于均聚物來說的,對于共聚物,請在理解上面文字的原理的前提上當意即妙考慮怎么構造寡聚物模型和獲得不同單元的原子電荷。
例6:產生跨越盒子的無限長氯丁二烯聚合物的拓撲文件
這個例子示例一下怎么得到跨越盒子無限延展的物質的拓撲文件,還是拿氯丁二烯聚合物作為例子。從此例會看到,對這種體系,只要mol2里記錄了跨盒子的鍵即可,在產生拓撲文件的過程上沒有絲毫特殊性。
examples\polymer\Neoprene_10_cross.mol2是一個聚合物為10,兩頭跨越盒子,因此無限延展的聚合物結構。用文本編輯器打開它的話,會看到@<TRIPOS>BOND字段里有4 1 73 1,這是說此體系中第4個鍵是1和73號原子之間的鍵,形式鍵級為1。1和73之間的鍵正是跨盒子的鍵(如果你用GaussView看的話,會看到有一個超長的從體系最左端的碳連到最右端的碳的鍵,這是因為GaussView沒按照周期性來顯示鍵)。
Sobtop構造拓撲文件的流程和構造孤立體系的流程沒有絲毫差異。啟動Sobtop,載入這個mol2文件后,屏幕上會看到
Atom X Y Z Charge Bonds to
1C -1.90300 0.55280 0.00000 0.00000 2H 3H 4C 73C
...略
可見確實1和73之間跨盒子的鍵被識別了。之后還是讓Sobtop自動指認GAFF原子類型,從屏幕的提示中會看到1和73都正確地識別成了sp3雜化的碳(c3原子類型)。然后照常產生itp和top文件,會看到itp里的[ bond ]項里也有1-73的,[ angles ]和[ dihedrals ]也都有涉及1-73的,顯然拓撲文件是很合理的。
在GROMACS里做此體系的動力學模擬的時候別忘了在mdp里要加上periodic-molecules = yes,這樣GROMACS才能考慮跨盒子的力場項,否則當前體系從拓撲關系上會被當做一個孤立的環狀體系看待。
對當前體系也可以照常讓Sobtop調用OpenBabel指認MMFF94原子電荷。如果要用RESP電荷算當前體系,顯然每個聚合物單元都應當套用模型體系算出來的中間單元的RESP電荷。
例7:產生碳納米管的GROMACS拓撲文件
值得一提的是,Sobtop不僅支持mol2也支持pdb作為輸入文件,但pdb無法給sobtop提供連接關系信息,因此sobtop會根據原子間距離自動猜,只要兩個原子距離小于它們的共價半徑和的1.15倍就被視為成鍵(閾值可以通過sobtop.ini里的bondcrit改)。注意載入后應看看屏幕上顯示的連接關系對不對,另外也可以用Sobtop主界面里的選項4導出記錄了當前判斷出的連接關系的gjf文件,然后載入GaussView用肉眼檢查。用pdb的話Sobtop沒法自動指認GAFF原子類型,需要以其它方式指認。
此例產生碳納米管拓撲文件,結構文件是examples\CNT\CNT.pdb。由結構可見,位于中間的碳形成三個C-C鍵,位于邊緣的碳形成兩個C-C鍵。這個體系適合所有碳原子都用ca(GAFF里的芳香碳原子)。如果以mol2作為輸入文件且讓Sobtop自動判斷GAFF原子類型的話會發現中間的碳判斷成了ca,而邊緣的碳判斷成了c1(GAFF里sp雜化的碳),明顯不合適。另外,由于當前體系一大堆環連在一起,讓Sobtop自動判斷GAFF原子類型時耗時也較高。由于以上原因,再加上本來此體系就只有一種原子類型,所以此例最適合手動直接設置原子類型。
啟動Sobtop,然后輸入
examples\CNT\CNT.pdb
1 //產生GROMACS拓撲文件
6 //手動設置原子類型
[回車] //選擇所有原子
ca //設為ca原子類型。然后屏幕上提示所有原子的原子類型已經被指認了
0 //進入下一步
4 //用預置的力場參數,缺失的自動猜
[回車] //用默認的top文件產生路徑
[回車] //用默認的itp文件產生路徑
現在當前目錄下就有了CNT.itp和CNT.top了。此例就不產生gro文件了,因為GROMACS的grompp也能基于pdb文件產生tpr文件。
實際上GAFF、OPLS-AA、GROMOS等力場里都沒有專供碳納米管(以及富勒烯、石墨烯)的力場參數,從itp里的注釋中可見當前的二面角參數是取的GAFF里的通配項X-ca-ca-X,故只算是湊合用。如果比如你想讓碳納米管剛性更強,一種做法是直接替換itp里當前的二面角的勢壘高度值(當前是15.167)為更大的值,也可以比如在bonded_param.dat的$DIHEDRAL字段末尾加上ca-ca-ca-ca項,在里面定義專供碳納米管、石墨烯等體系用的純ca間的二面角的參數。另外,從當前itp文件里[ dihedrals ] ; impropers下面的部分還可以看到用于維持平面性的improper項,實際上GAFF沒有定義ca-ca-ca-ca的improper參數,從注釋上可知當前Sobtop自動套用了X -X -ca-ha的improper項,大家也可以在bonded_param.dat的$IMPROPER里加上相應的專一項。
例8:產生金剛石晶體的拓撲文件
這個例子產生金剛石晶體的拓撲文件,主要有三個目的,一是讓大家了解怎么寫assign_AT.dat文件自定義判斷原子類型的規則,另一方面也讓大家看到用pdb作為輸入文件的時候Sobtop也可以直接判斷跨盒子的鍵,本節還示例怎么自定義新的原子類型。看過這個例子,大家肯定就再也不想用GROMACS里的x2top了,因為Sobtop靈活強大太多了。
assign_AT.dat文件類似于x2top的n2t文件,具體語法和例子在此文件里的注釋中寫得都很明白。Sobtop可以根據此文件里定義的判斷規則,根據原子在結構文件中所處的實際環境指定原子類型和原子電荷。可以用的判斷條件包括元素、成鍵的總數、必須滿足的與其它元素原子成鍵的形式鍵級、與其它元素原子的距離范圍。相比之下,x2top程序的n2t文件只能通過元素間距離在某個值附近來指定原子類型,遠沒有Sobtop靈活。
此例用的結構文件是examples\diamond_3x3x3\diamond_3x3x3.pdb,是Multiwfn基于金剛石原胞構造出的3*3*3超胞的pdb文件,注意此文件開頭的CRYST1字段定義了此超胞的盒子信息,這個信息必須有,否則Sobtop沒法根據坐標判斷出跨盒子的鍵。當前這個體系其實用sp3雜化的碳作為原子類型就可以,可以直接手動在Sobtop里把所有原子都定義成c3原子類型,但這里作為示例我們通過定義assign_AT.dat來實現。當前體系特征非常簡單,每個原子都形成四個鍵,故只需要把規則定義為形成四個鍵的原子就被設為某種原子類型即可。在assign_AT.dat中加入以下內容
$Ctest
nbond 4
代表成鍵數為4的原子的原子類型都判斷為Ctest。原子類型后面沒寫原子電荷,因此此類原子的原子電荷默認為0。
由于Ctest原子類型是我們自己創造的,故還要在LJ_param.dat中加入它的LJ參數和注釋(分號后的信息):
Ctest 1.9255 0.105 ;My carbon
此體系涉及Ctest原子類型間的bond、angle、dihedral項,因此在bonded_param.dat的$BOND里加入自定義參數:
Ctest-Ctest 300.9 1.5375
在$ANGLE里加入:
Ctest-Ctest-Ctest 62.9 111.51
在$DIHEDRAL里加入:
Ctest-Ctest-Ctest-Ctest 1 0.18 0.0 -3.
Ctest-Ctest-Ctest-Ctest 1 0.25 180.0 -2.
Ctest-Ctest-Ctest-Ctest 1 0.20 180.0 1.
PS:上面這些參數值其實都是GAFF里c3間的參數,本節借用過來作為例子。
啟動Sobtop,依次輸入
examples\diamond_3x3x3\diamond_3x3x3.pdb //載入后從屏幕上可見所有原子都形成四個鍵,確實Sobtop根據原子坐標和晶胞信息猜的成鍵都是對的
1 //產生GROMACS拓撲文件
5 //根據assign_AT.dat里的規則指認原子類型。之后從屏幕上的信息可見所有原子都被指認為了Ctest原子類型
0 //進入下一步
4 //用預置的力場參數,缺失的自動猜
[回車] //用默認的top文件產生路徑
[回車] //用默認的itp文件產生路徑
現在當前目錄下就有了diamond_3x3x3.itp和diamond_3x3x3.top,與diamond_3x3x3.pdb相結合就可以做動力學模擬了。還是別忘了mdp里要加上periodic-molecules = yes。
實際上,這一節的例子如果把$Ctest設為$c3的話,就不需要自己修改LJ_param.dat和bonded_param.dat了,因為本身GAFF力場里就有c3-c3、c3-c3-c3、c3-c3-c3-c3的參數。
效仿本例,大家也可以產生比如碳化硅、氮化硼等原子晶體的拓撲文件,只要在assign_AT.dat里根據情況定義多個判斷原子類型的規則,并在bonded_param.dat里添加涉及的bonded參數即可。諸如硼這種元素,既可以在assign_AT.dat里定義你專門命名的原子類型并在LJ_param.dat里定義LJ參數(B的LJ參數可借用里面的UF_B的,即UFF力場的硼的);也可以在assign_AT.dat里定義$UF_B項,即設置UFF力場的B的原子類型判斷規則,這樣LJ_param.dat里就不用再新加一個原子類型了。由于碳化硅、氮化硼這些體系是化合物而非單質,因此別忘了應當在assign_AT.dat里也把相應原子類型要用的原子電荷寫上去。原子電荷用Multiwfn基于CP2K產生的molden文件計算CM5電荷通常就不錯(可參看http://www.shanxitv.org/588了解CP2K的molden文件怎么產生。Multiwfn載入后用主功能7的子功能16就能計算CM5電荷)。
對于MOF、COF、黑磷、MoS2、MXene等更復雜的周期性體系,也都可以在這一節的做法基礎上舉一反三。對于里面涉及到的缺的bonded參數,平衡鍵長/鍵角/二面角可以直接用晶體結構里的值,結合自己猜的足夠大的力常數(只要力常數夠大,就能維持剛性。對這些體系通常我們對其實際的結構波動程度不感興趣,故不用太較真,只要能維持晶體結構穩定就行了)。也可以截取不同的簇并飽和邊緣,用Gaussian或ORCA做振動分析,然后通過Sobtop得到里面涉及的一些項的力常數(顯然,被求力常數的這些內坐標應當處于簇模型的靠中間部分,以避免邊界效應)。也可以用CP2K直接對周期性體系做振動分析,輸出文件里有Hessian矩陣,Sobtop可以靠它來一次性得到所有力常數。這類體系,由于直接與真空區接觸,原子電荷可以CP2K算REPEAT電荷是不錯的選擇,輸入文件可以直接由Multiwfn創建,見《使用Multiwfn非常便利地創建CP2K程序的輸入文件》(http://www.shanxitv.org/587)。
例9:產生二氧化硅晶體的拓撲文件
這一節用Sobtop產生二氧化硅晶體的拓撲文件,本例要把Chem. Mater., 26, 2647 (2014)一文中給出的體相二氧化硅的原子電荷、LJ 12-6參數、bond和angle項參數補進來,見此文表3和表4。此文作者認為對這種體系沒必要考慮dihedral項,因此沒給dihedral參數。
用VMD自帶的Extensions - Modeling - Inorganic Builder可以構建二氧化硅晶體結構。本例用的examples\SiO2_3x3x3\SiO2_3x3x3.pdb是通過此工具構造的二氧化硅3*3*3超胞。
根據文獻里的原子電荷和LJ參數,在LJ_param.dat里寫入以下信息,新增加兩種原子類型
Si_bulk 2.075 0.093 ;Si in bulk SiO2
O_bulk 1.735 0.054 ;O in bulk SiO2
注意:Chem. Mater., 26, 2647 (2014)文獻里公式1的形式是錯的,此式里的σ實際上應當為r0(勢能曲線極小點的距離),因此上面LJ_param.dat里定義的(r0)/2項是此文里給的σ的一半(而如果你看到的某篇文獻里給的是一般意義的σ參數,即勢能曲線為0處的距離,那么在LJ_param.dat里應當按照(r0)/2 = 2(1/6)*σ/2的關系來定義)。
在bonded_param.dat里的$BOND中添加
Si_bulk-O_bulk 285 1.65
在$ANGLE里添加
O_bulk-Si_bulk-O_bulk 100 109.5
Si_bulk-O_bulk-Si_bulk 100 149.0
在assign_AT.dat里添加以下內容。1.1和-0.55分別是體相二氧化硅里硅和氧的原子電荷。這里寫element Si要求元素必須是硅,這是為了能和上個例子里形成四個鍵的碳進行區分。
$Si_bulk 1.1
nbond 4
element Si
$O_bulk -0.55
nbond 2
element O
啟動Sobtop,然后輸入
examples\SiO2_3x3x3\SiO2_3x3x3.pdb
1 //產生拓撲文件
5 //根據assign_AT.dat里的規則指認原子類型。之后從屏幕上的信息可見所有原子都被正確指認了原子類型
0 //進入下一步
4 //用預置的力場參數,缺失的自動猜
[回車] //用默認的top文件產生路徑
[回車] //用默認的itp文件產生路徑
現在當前目錄下就有了SiO2_3x3x3.itp和SiO2_3x3x3.top。由于我們沒有添加二面角參數,故當前itp里[ dihedrals ]中的參數都是缺失的,Sobtop自動對它們填的參數里勢壘為0,因此等價于不設這些二面角項。你如果把整個[ dihedrals ]字段全都刪掉也可以,不影響結果。如果你把sobtop.ini里的iskipgendih設為1,那么所有二面角項都將不會產生(對很大的體系,這能大幅節約產生和寫入二面角項的耗時!)。[ pairs ]不要刪掉,因為Chem. Mater.這篇文章里在給出參數的時候是已經考慮了計算1-4作用項時優化的參數。
之后將SiO2_3x3x3.itp、SiO2_3x3x3.top和SiO2_3x3x3.pdb相結合,就可以試試跑模擬了。在examples目錄下給了個md_PBC.mdp,是在NPT下跑100 ps常溫下動力學的文件,由于當前晶體尺寸很小,為了cutoff不超過晶胞最短邊長的一半,所以cutoff設得很小。運行以下命令即可執行測試。觀看軌跡會發現晶體結構和晶胞尺寸維持得特別好,極度理想!
gmx grompp -f md_PBC.mdp -c SiO2_3x3x3.pdb -p SiO2_3x3x3.top -o md.tpr -maxwarn 2
gmx mdrun -v -deffnm md
例10:產生金屬有機框架化合物MOF-5的拓撲文件
此例對金屬有機框架化合物MOF-5的單胞產生拓撲文件,此體系里有Zn-O配位鍵。examples\MOF-5\MOF-5.cif是MOF-5的晶體文件。用GaussView打開它,確認里面成鍵關系(特別是Zn-O配位鍵)都是正確的,然后保存為mol2文件。
注1:此例用mol2而非pdb格式原因有二:(1)Sobtop會從mol2文件里讀取連接關系,并且連接關系可以在GaussView里可視化檢驗正確性。若用pdb格式的話,Sobtop根據原子間距離和元素半徑自動猜的連接關系和你期望的未必一致,特別是對MOF-5這樣鍵的類型略多的情況檢驗起來略麻煩 (2)用mol2格式時才能讓Sobtop自動指認GAFF原子類型。 如果你打算自己寫合適的assign_AT.dat來定義原子類型的判斷規則,那么像之前的例子一樣用pdb格式也完全可以
注2:如果你要讓Sobtop自動猜參數缺失的參數,應當按照http://www.shanxitv.org/soft/Sobtop/#FAQ14里說的,手動給GaussView保存出來的mol2文件里加入晶胞信息,否則Sobtop不會將之視為周期性體系,猜的參數可能完全錯誤!
啟動Sobtop,輸入
examples\MOF-5\MOF-5.mol2
1 //產生拓撲文件
2 //指認GAFF類型,缺失的用UFF的。從屏幕提示上會看到Zn被指認為了UFF原子類型,其它的都是GAFF原子類型
4
//用預置的力場參數,缺失的自動猜
[回車] //用默認的top文件產生路徑
[回車] //用默認的itp文件產生路徑
現在當前目錄下就有了MOF-5.itp和MOF-5.top。從屏幕上會看到缺失的bond項有UF_Zn-os,缺失的angle項有:
UF_Zn-os-UF_Zn
UF_Zn-os-ce
os-UF_Zn-os
os-ce-os
缺失的dihedral項有:
UF_Zn-os-UF_Zn-os
UF_Zn-os-ce-os
UF_Zn-os-ce-ca
os-UF_Zn-os-ce
os-ce-ca-ca
缺失的improper項有ce-ca-ca-ca,這個不用管,原因同前。
對于不是特別講究的情況,bond和angle項缺就缺了,如上可見這些都是涉及Zn的,而這些項的剛性也都很強,就用Sobtop默認猜的參數就能維持其狀態。缺失的dihedral項里前四個都是和Zn有關的,你去看MOF-5的結構圖的話,就會認識到這些二面角項沒有也沒什么關系,因為光靠angle項也足矣維持每4個Zn所在的局部單元的剛性和立體結構,從實際模擬得到的動力學軌跡上也會發現這一點。
缺的os-ce-ca-ca這個二面角完全無視掉不太好。如果直接用Sobtop默認賦予它的參數(旋轉勢壘為0)的話,模擬過程中苯環部分會過于自由、快速地旋轉,而實際中旋轉不可能容易到這種程度。當然構建個模型體系做二面角掃描再去擬合參數是最好的,但步驟略繁瑣。這里要求不高,就試圖直接用力場庫文件里其它的二面角參數湊合借用一下。這個二面角中間部分其中一個是ca,對應于苯環上的碳,這個是必須要存在的,而且不好替換成別的原子類型。因此要在bonded_param.dat的$DIHEDRAL部分里搜索-ca-,可以找到好多項,再要求另一個中間部分原子的元素也必須是碳,篩選出的就只剩下面這些了,都是通配項:
X -c -ca-X 4 4.000 180.000 2.000
X -c1-ca-X 2 0.000 180.000 2.000
X -c2-ca-X 4 2.800 180.000 2.000
X -c3-ca-X 6 0.000 0.000 2.000
X -ca-ca-X 4 14.500 180.000 2.000
X -ca-cp-X 4 14.500 180.000 2.000
X -ca-cq-X 4 14.500 180.000 2.000
去看LJ_param.dat里原子類型的定義可知,我們缺的二面角中間部分涉及的ce原子類型對應的是共軛區域內部的sp2雜化的碳原子(但不是在芳環中),在當前體系中它來自于羧基碳,上面的這些通配項里中間原子與之最接近的就是c2了(sp2雜化的碳,只不過沒限定是處于共軛部分)。因此我們把X -c2-ca-X 4 2.800 180.000 2.000復制并修改為
os-ce-ca-ca 4 2.800 180.000 2.000 ;Same as X -c2-ca-X
并插入到$DIHEDRAL里面。我習慣上插入到字段開頭,便于修改。
之后再按照前面的過程重新生成一次MOF-5的拓撲文件,我們添加的os-ce-ca-ca就會被利用了,產生拓撲文件時屏幕上也不會提示缺os-ce-ca-ca的參數了。
之后我們跑一下MOF-5的動力學。給grompp用的結構文件可以這么得到:啟動Multiwfn,載入MOF-5.cif,進入主功能100的子功能2,選擇保存為pdb文件,命名為MOF-5.pdb(注意不能用GaussView載入cif文件后保存pdb文件,因為這樣的pdb里沒有盒子信息。也不能用Sobtop將mol2文件轉成gro文件直接用于動力學模擬,是因為GaussView保存的mol2里沒盒子信息,因此得到的gro文件里也沒盒子信息,除非你自己修改其最后一行把晶胞信息如實填進去)。
之后將上面得到的MOF-5.pdb、MOF-5.itp、MOF-5.top以及examples\md_PBC.mdp相結合,就可以跑動力學模擬了:
gmx grompp -f md_PBC.mdp -c MOF-5.pdb -p MOF-5.top -o md.tpr -maxwarn 3
gmx mdrun -v -deffnm md
從模擬軌跡中會看到MOF-5的形態維持得很好,該剛性的部分一直保持剛性,苯環也可以適度地發生旋轉。
以上沒有涉及原子電荷,顯然如果你要研究比如MOF-5吸附其它分子之類的模擬,原子電荷肯定得有。CP2K程序對周期性體系可以算REPEAT電荷,對靜電勢重現性較好,故很適合這種目的。算REPEAT電荷的CP2K輸入文件靠Multiwfn能很容易地生成,看《使用Multiwfn非常便利地創建CP2K程序的輸入文件》(http://www.shanxitv.org/587)。examples\MOF-5\CP2K_REPEAT.in是這個任務的輸入文件,同名的.out是輸出文件,CP2K_REPEAT-RESP_CHARGES.resp是此任務產生的記錄各個原子REPEAT電荷的文件,把里面的電荷拷到MOF-5.itp里即可。
基于簇模型和Hessian矩陣獲得缺失的和Zn有關的參數
如果你做的模擬要用于發表文章,對于以上那些和Zn有關的缺失的參數,我建議還是自己搞出來,這樣更嚴格,這其實很容易。因為根據MOF-5的結構可看出,和Zn有關的缺失的二面角參數用諧振勢描述就夠了,因此都不涉及到費事的勢能面掃描。具體做法是把MOF-5復制成超胞,再把其中一個連接單元截出來并且把邊緣用氫飽和,然后做優化和振動分析得到fch文件。這個模型的Gaussian輸入文件是examples\MOF-5\MOF-5_fragment.gjf,產生的fch文件比較大,我就不提供了。將優化完的結構保存為MOF-5_fragment.mol2。
作為例子,這里獲得缺失的UF_Zn-os-UF_Zn的angle參數。當前的模型是Td點群對稱的,隨便找對應這個鍵角類型的三個原子,從MOF-5_fragment.gjf中可見Zn11-O15-Zn13就適合。啟動Sobtop,然后依次輸入
MOF-5_fragment.mol2
-1 //修改產生力常數的方法
3 //m2Seminario。當前體系有對稱性,此方法此時比mSeminario更好,用后者算等價的鍵角可能得到明顯不同的結果。此體系也不太適合DRIH,會發現有些鍵角力常數明顯太小
5 //對特定鍵、鍵角、二面角產生剛性參數
MOF-5_fragment.fchk
11,15,13 //對應Zn11-O15-Zn13
屏幕上明確提示了加入$ANGLE里要用的參數,應當在里面加入UF_Zn-os-UF_Zn 22.34 109.47。
此體系的缺失的二面角如UF_Zn-os-UF_Zn-os不太好弄成諧振勢形式加入$HARM_DIH,因為此體系里這種二面角有0、120、-120度三種情況,因此比較適合弄成有三個極小點的周期勢函數,但對這個體系不太好做旋轉掃描。可以用當前功能得到其諧振力常數,然后自己設計一個周期勢函數,使其在0、120、-120度極小點處的力常數與m2Seminario得到的大致吻合。
如果你只需要模擬一個MOF-5單胞的話,其實補全參數最簡單的方法是直接用CP2K做優化和振動分析,基于Hessian自動算出來所有缺的力場項,參考下面FeCl2的例子。
例11:產生FeCl2二維材料的拓撲文件
此例產生FeCl2層狀二維材料的拓撲文件。這一節分為兩部分,第一部分演示懶人法,直接基于晶體結構,靠Sobtop猜的參數;第二部分做法較嚴格,基于CP2K對周期性結構做振動分析得到Hessian,讓Sobtop嚴格地計算出所有力常數,但CP2K里振動分析的耗時很高。如果你只為了在模擬過程中維持FeCl2的剛性結構,用第一種做法就行了。FeCl2的原胞很小,僅有一個Fe和兩個Cl,這一節我們對FeCl2的5*5*1的超胞結構來產生拓撲文件。
1 懶人法:直接用猜的參數
首先構造超胞的pdb文件。啟動Multiwfn,載入Sobtop目錄下的examples\FeCl2\FeCl2.cif,然后輸入
300 //主功能300
7 //幾何操作
19 //構造超胞
5 //第一個晶格矢方向復制成原先的5倍
5 //第二個晶格矢方向復制成原先的5倍
1 //第三個晶格矢(垂直于界面方向)保持不變
-2 //保存pdb文件
FeCl2_5x5x1.pdb //要產生的文件名
啟動Sobtop,載入FeCl2_5x5x1.pdb,從屏幕上的信息可見每個Fe都是與6個Cl相連,每個Cl都與3個Fe相連,自動猜的連接關系符合實際情況。然后輸入
1 //產生拓撲文件
6 //手動設置原子類型(由于當前用的不是mol2文件,故接下來沒法選2先自動用GAFF原子類型然后用UFF的補全,而需要手動設置原子類型)
Cl //選擇所有氯原子
cl //設為GAFF力場的cl原子類型
1 //對未設定類型的原子用UFF原子類型
0 //進入下一步
4
//用預置的力場參數,缺失的自動猜
[回車] //用默認的top文件產生路徑
[回車] //用默認的itp文件產生路徑
屏幕上提示的缺失力場項如下,很正常
Missing bond types (equilibrium length comes from current geometry, while force
constant is arbitrary guessed):
1 cl-UF_Fe
Missing angle types (equilibrium angle comes from current geometry, while force
constant is arbitrary guessed):
1 cl-UF_Fe-cl
2 UF_Fe-cl-UF_Fe
Missing dihedral types (rotational barriers are set to zero):
1 cl-UF_Fe-cl-UF_Fe
前述的FeCl2_5x5x1.pdb里有晶胞信息,Z方向真空層足夠大,不做修改就可以直接用于跑動力學。現將FeCl2_5x5x1.pdb、FeCl2_5x5x1.itp、FeCl2_5x5x1.top以及examples\md_PBC.mdp相結合,跑動力學模擬:
gmx grompp -f md_PBC.mdp -c FeCl2_5x5x1.pdb -p FeCl2_5x5x1.top -o md.tpr -maxwarn 3
gmx mdrun -v -deffnm md
gmx trjconv -s md.tpr -f md.xtc -o md_new.xtc -ur tric -pbc atom,選system
gmx trjconv -s md.tpr -f md.gro -o md_new.gro -ur tric -pbc atom,選system
上面用trjconv處理是為了讓新得到的xtc和gro文件里原子坐標的周期性記錄方式滿足當前FeCl2的三斜晶胞的情況。基于md_new.gro/xtc觀看軌跡,可見FeCl2二維結構維持得非常好,剛性很強。
2 較嚴格方法:基于CP2K周期性DFT計算的Hessian得到參數
首先用Multiwfn載入FeCl2.cif,參考《使用Multiwfn非常便利地創建CP2K程序的輸入文件》(http://www.shanxitv.org/587)創建一個對FeCl2原胞的變胞幾何優化任務的CP2K輸入文件,用5*5*1 k點,優化過程中保持晶胞夾角和空間群不變,輸入文件是examples\FeCl2\exact\opt.inp。跑完后,用Multiwfn載入優化任務產生的restart文件,在創建CP2K輸入文件的界面里擴成5*5*1的超胞,再創建CP2K振動分析的輸入文件,見examples\FeCl2\exact\freq.inp。之后退回到Multiwfn主菜單(如果你已經把Multiwfn關了也沒事,此時載入freq.inp即可),用主功能100的子功能2將當前超胞結構保存成pdb文件opted_5x5x1.pdb。用CP2K跑freq.inp,從其輸出文件freq.out中搜Hessian的話會看到里面確實記錄了笛卡爾Hessian矩陣。
啟動Sobtop,然后輸入
opted_5x5x1.pdb
-1 //修改基于Hessian計算力常數的方法
3 //m2Seminario。當前體系里有大量等價的鍵、鍵角、二面角,用默認的mSeminario的話得到的力常數會嚴重違背等價性。DRIH對于當前體系不理想故不考慮
1 //產生拓撲文件
6 //手動設置原子類型
Cl //選擇所有氯原子
cl //設為GAFF力場的cl原子類型
1 //對未設定類型的原子用UFF原子類型
0 //進入下一步
2 //所有參數皆基于Hessian計算產生
examples\FeCl2\exact\freq.out //CP2K振動分析輸出文件
[回車] //用默認的top文件產生路徑
[回車] //用默認的itp文件產生路徑
現在當前目錄下就有了opted_5x5x1.itp和opted_5x5x1.top。當前目錄下還產生了FFderiv.txt,里面記錄了基于Hessian計算的每種類型的bond、angle、dihedral項的計算次數、所得參數的平均值/最小值/最大值。雖然當前用了m2Seminario方法計算力常數,但由于CP2K靠有限差分算的Hessian有一定數值誤差,因此等價的力場項的力常數仍可能多多少少有點差異(但至少遠小于用mSeminario的情況)。
順帶一提,對于Sobtop對當前體系算出來的最小值和最大值相差很小的bonded參數,可以考慮將它的平均幾何參數和平均力常數添加到bonded_param.dat里,以供今后對類似體系直接使用。例如在FFderiv.txt里可見的如下所示的UF_Fe-cl-UF_Fe項就是這種情況,對當前體系總共算了150次,其鍵角和力常數的min和max值都比較相近,因此可以挪到bonded_param.dat里。
2 Type: UF_Fe-cl-UF_Fe Times of evaluations: 150
Avg, min, max of k: 675.30 664.43 695.20 kcal/mol/rad^2
Avg, min, max of a0: 94.548 94.507 94.576 degree
但FFderiv.txt里比如cl-UF_Fe-cl型的鍵角項的最小鍵角為85.437,最大鍵角為179.990,差異太大,因此不適合將此類鍵角項的平均參數放到bonded_param.dat里。
為檢驗拓撲文件合理性,用下面的命令跑動力學
gmx grompp -f md_PBC.mdp -c opted_5x5x1.pdb -p opted_5x5x1.top -o md.tpr -maxwarn 3
gmx mdrun -v -deffnm md
gmx trjconv -s md.tpr -f md.xtc -o md_new.xtc -ur tric -pbc atom,選system
gmx trjconv -s md.tpr -f md.gro -o md_new.gro -ur tric -pbc atom,選system
從結果可見FeCl2的結構很好地維持住了。
別忘了,當FeCl2板和其它分子一起模擬時,需要給FeCl2賦予原子電荷,REPEAT電荷就不錯。
例12:產生非標準殘基的rtp文件字段
待寫
Skill
技巧1:通過命令行方式使用Sobtop
Sobtop通常是以交互式方式使用,但也可以通過命令行方式使用,這可以很方便地將Sobtop嵌入shell腳本里,實現對大批量分子一鍵產生拓撲和gro文件。相關知識非常建議閱讀《詳談Multiwfn的命令行方式運行和批量運行的方法》(http://www.shanxitv.org/612)。下面僅舉個最簡單的例子,通過命令行方式運行Sobtop產生examples\7-helicene.mol2的itp、top、gro,假設用戶用的是Linux環境。
先創建一個文本文件,比如叫genGAFF.txt,內容是在Sobtop里要輸入的所有命令,將之放到Sobtop目錄下
2
1
2
4
0
然后進入Sobtop目錄,輸入./Sobtop examples/7-helicene.mol2 < genGAFF.txt > 7-helicene.txt。之后當前目錄下就有了7-helicene.itp、7-helicene.top、7-helicene.gro了。7-helicene.txt是Sobtop在屏幕上輸出的那些信息,供檢查用。
FAQ
FAQ 1:我是初學者,我用Sobtop產生完了拓撲文件,然后怎么用GROMACS做模擬啊?
這個不是Sobtop的問題,而是GROMACS本身使用的問題,問出這種問題說明沒有最基本的GROMACS常識性知識。怎么做模擬這種問題絕對不是三言兩語就能說明白的,不懂分子動力學和分子力場的基礎知識、不懂GROMACS的各種常見子程序的用法和拓撲文件的規則、不懂對于各類問題模擬的流程和應注意的要點,根本不可能開展模擬。零基礎的話光是照著網上零七八碎的教程(其中往往還有錯,或者是過時的)去摸索效率極低,還容易被坑,還學得很不系統。我非常建議參加北京科音分子動力學與GROMACS培訓班(http://www.keinsci.com/workshop/KGMX_content.html)特別系統性地學一遍,就全都明白了,在理解和使用Sobtop上也不會再遇到任何障礙。沒機會到現場參加也沒關系,往屆的培訓資料都可以隨時購買,頁面里明確說了。
FAQ 2:為什么Sobtop直接支持的是GAFF和AMBER力場,會支持CGenFF、OPLS-AA、GROMOS么?
Sobtop支持的函數形式是GAFF/AMBER力場的,而且自帶力場庫包含的是GAFF/AMBER的參數,主要原因在于Sobereva比較喜歡GAFF/AMBER力場,因為:
(1)力場函數簡單且標準
(2)流行程度極高
(3)被分子模擬程序支持廣泛,在GROMACS中完美支持
(4)原子類型數目適中,而且名字非常簡潔清楚,因此很便于管理
(5)基于RESP這種標準化的原子電荷計算方法
(6)GAFF(普通有機體系)、AMBER(蛋白質和核酸)、GLYCAM(糖類)、Lipid/Slipids(磷脂)、Merz(各類單原子離子)等高質量力場彼此兼容,組合使用可以描述無機以外的幾乎各種體系
CGenFF、OPLS-AA、GROMOS力場能描述的體系用GAFF/AMBER也都能描述,GAFF/AMBER不能描述的那些力場也基本不能描述,因此我不覺得有任何必要去直接支持GAFF/AMBER以外的力場。
由于前面說了,Sobtop的力場庫文件是完全開放的,所以如果其它力場里的力場項的函數形式與GAFF/AMBER的相同,也可以把它們的相應類型參數自己寫到Sobtop的力場庫里使用。Sobtop能創建的拓撲文件涉及的力場絕不僅限于GAFF/AMBER,諸如例子里的二氧化硅、MOF-5、FeCl2體系都是完全脫離了GAFF/AMBER力場,過渡金屬配合物的例子里對過渡金屬也用了GAFF/AMBER以外的原子類型。
AMBER和GAFF的關系:我在答疑時經常發現有不少人搞不清楚GAFF和AMBER的關系,居然試圖用AMBER描述小分子。這里順帶強調一下二者的聯系。AMBER力場提出的目的是用于蛋白質和核酸的模擬,里面的原子類型、參數主要也是面向這類生物大分子的。雖然AMBER也不是不能描述有機小分子,但是很容易遇到缺適合的原子類型、缺需要的bonded參數的情況。后來提出的GAFF的全稱是“廣義化的AMBER力場”,開發的目的專門是用來描述各種雜七雜八的有機小分子,并且它和AMBER力場形式完全兼容,也有不少參數是直接從AMBER繼承過來的,最適合搭配的都是RESP原子電荷(搭配后來的RESP2(0.5)也很好)。GAFF的原子類型都是小寫,而AMBER是大寫,所以二者一起用的時候不會有沖突。GAFF的原子類型比AMBER能更充分涵蓋五花八門的有機小分子的情況,bonded參數也更豐富,顯然對于小分子體系沒有任何理由用AMBER而不用GAFF。而對于蛋白質和核酸,雖然比AMBER更普適的GAFF也能描述,但明顯不如AMBER這樣的對它們專一性的力場更合適。所以,對于諸如蛋白質+配體這種類型體系的模擬,蛋白質要用AMBER,配體要用GAFF。由于GAFF比較普適,所以只要不是像蛋白質、核酸這樣有專門特別適合的力場的情況,都可以用GAFF。但也有很多情況有些元素不被GAFF支持,有些處于非典型電子結構狀態的原子在GAFF里沒有合適的原子類型能夠描述,或者特殊的bonded參數在GAFF里沒有,這時候就要考慮讓Sobtop用UFF原子類型進行補充,以及用Sobtop基于Hessian矩陣計算力常數或者通過從文獻搜索/自行擬合等其它方式補充bonded參數。
FAQ 3:Sobtop會支持AMBER、Lammps等程序么?
Sobtop沒打算支持其它程序。我不用Lammps。AMBER我如今也很少用,因為GROMACS更快更靈活還完全免費,常用功能都有,還有gmx_mmpbsa、PLUMED等很多工具能擴展其功能。而且AMBER的拓撲文件相比于GROMACS過于復雜。
如果你是其它動力學程序的用戶,也可以先借助Sobtop產生GROMACS的拓撲文件,然后再轉成其它程序的。比如GROMACS拓撲文件轉成AMBER、NAMD和CHARMM的可以用ParmED,轉成Lammps的可以用GRO2LAM或InterMol。
FAQ 4:離子液體的拓撲文件怎么產生?陰陽離子合在一起處理還是分開處理?
我強烈建議(這也是GROMACS的良好使用習慣),對體系中每種通過bond項(用于描述有顯著共價特征的化學鍵,也包括配位鍵)連接到一起的片段,單獨用Sobtop產生拓撲文件,并計算其原子電荷。比如離子液體,由于陰陽離子間是非共價作用,因此就應當對陽離子和陰離子分別照常產生拓撲文件,然后把二者的itp都include到主top文件里,并且把二者中的[ atomtypes ]也都合并、去重后挪到主top的[ defaults ]字段的后面去。
然而,在使用Packmol等程序構造離子液體凝聚相模型時,應當把陰陽離子對作為一個整體來插入,而不建議陰、陽離子單獨插入,否則產生的結構文件里可能恰好有的地方陽離子非常密集、有的地方陰離子非常密集,顯然模擬一開始,帶相同電荷的離子間會產生強烈的靜電互斥,可能造成動力學不穩定甚至崩潰。而陰陽離子作為一對來插入的話,至少保證了每個離子旁邊有個相反電荷的離子,不至于出現局部靜電互斥特別大的情況。別忘了top里的[ molecules ]字段中離子出現順序要和結構文件里相一致。
再順帶一說,有些零基礎的GROMACS初學者,居然試圖用Sobtop產生整個模擬體系(含一大堆分子)的拓撲文件,然后直接用GROMACS跑,這種做法是極端野蠻粗暴不可理喻的。這么做可能在Sobtop里耗時很高,而且所有分子的拓撲信息都攪合在一起作為一個[ moleculetype ]出現時會特別難以管理和修改。
FAQ 5:Sobtop可以產生生物分子的拓撲文件么?
可以是可以,畢竟Sobtop是普適的,對這類體系可以用AMBER原子類型。但是對這類體系我不建議優先考慮Sobtop,因為GROMACS里自帶的pdb2gmx是專門干這個的,它對蛋白質、核酸會自動修改質子化態、加氫、處理末端殘基,這些對生物分子專一性的處理在Sobtop里是沒有考慮的,而且pdb2gmx直接就會根據rtp文件里的定義去給氨基酸/核苷酸殘基設置力場原文里定義的原子電荷。
Sobtop倒是可以幫助你去產生非標準殘基的rtp文件。具體來說就是給標準殘基兩側用ACE(乙酰)和NME(甲胺)封閉成模型分子,用量子化學程序做優化和振動分析。基于得到的波函數文件,用Multiwfn對非標準殘基計算RESP電荷(將中間的殘基的凈電荷約束為整數)并得到chg文件。然后在Sobtop里載入模型分子的mol2文件,從chg文件里載入原子電荷,然后選產生rtp文件,其間讓Sobtop自動指認AMBER原子類型,并要求從預置力場庫文件里取合適的參數,缺的讓Sobtop基于Hessian算出來。然后把得到的rtp里的字段恰當處理從分子狀態變成為殘基狀態后,挪到GROMACS力場目錄下的rtp里。
FAQ 6:Sobtop會支持二面角扭轉勢的擬合么?
早晚會。其實這個功能實現的復雜度完全看做到什么程度。最簡單的情況,只是擬合一個二面角勢函數而且不考慮優化phase參數,是非常簡單的事。但如果同時考慮擬合多個二面角,還要優化phase參數,就是個非線性全局優化問題,這就很復雜了。
FAQ 7:我希望Sobtop能像MCPB.py那樣產生金屬蛋白的拓撲信息
做成MCPB.py那樣原理上是不可能的,因為GROMACS產生蛋白質拓撲文件是靠pdb2gmx,其邏輯和AMBER的leap程序截然不同。
實際上目前的Sobtop就能對付金屬蛋白體系,只不過手動操作較多。比如蛋白質中有個鐵硫簇和幾個殘基結合,那就把鐵硫簇以及與之配位的殘基摳出來并飽和邊緣作為一個簇模型,凍結邊緣原子,用量子化學程序做優化和振動分析,參考《要善用簇模型,不要盲目用ONIOM算蛋白質-小分子相互作用問題》(http://www.shanxitv.org/597)里的討論。然后對這個簇模型照常用Sobtop產生拓撲文件、計算RESP電荷,其間要求Sobtop完全基于Hessian矩陣獲得力常數。然后,用Sobtop對鐵硫簇部分產生參數為空的itp文件。之后,用pdb2gmx對純蛋白質部分產生[ moleculetype ],將之和鐵硫簇部分的[ moleculetype ]都引入主top里,把鐵硫簇和配位殘基的原子電荷替換為前面計算的,把鐵硫簇部分的參數替換為簇模型中鐵硫簇部分的。然后在top里加入[ intermolecular_interactions ]字段,這個字段里可以設置跨越不同[ moleculetype ]的bonded項,其中加入[ bonds ]、[ angles ]、[ dihedrals ]字段,在里面手動寫入鐵硫簇和配體殘基兩部分之間的bond、angle、dihedral項和之前對簇模型計算的相應的力場參數。注意這里寫的原子序號是全局序號,即整個模擬的結構文件里對應的序號。
以上過程手動操作的地方較多,以后筆者可能會寫輔助工具,給出個盡量自動化的方案,減少手動操作量。
FAQ 8:為什么我基于Sobtop產生的拓撲文件,在動力學模擬過程中GROMACS報錯/崩潰?
mdrun跑計算期間不穩定乃至崩潰時會在屏幕上顯示LINCS warning、segmentation fault、1-4 interaction not within cut-off、water can not be settled、Listed nonbonded interaction between particles ... larger than the table limit其中的一種或多種提示,還可能產生一大堆step****.pdb文件。為什么出現、怎么解決和Sobtop沒直接關系,純粹是GROMACS使用層面上的問題。既然出現崩潰,自然得試圖搞清楚為什么崩潰,必須得基于GROMACS和分子動力學的常識性知識一點點排查,包括初始結構合理性、拓撲文件里的項和參數的合理性、mdp設置的合理性,都合理了自然就不崩潰了。下面就全面地談一下,以下文字無論你是否用Sobtop產生拓撲文件都是適用的。PS:如果你欠缺GROMACS和分子模擬基本理論知識的話,光看以下文字可能也很難理解和解決問題,這種情況我強烈建議參加北京科音分子動力學與GROMACS培訓班(http://www.keinsci.com/workshop/KGMX_content.html)從頭非常系統性地學一遍,問題就全都明白了,也不會因為犯各種低級錯誤導致模擬老也沒法順利跑起來。
檢查mdp:從mdp角度排查來說,首先確保設置的合理性。如果你覺得原本是合理的,在合理的范疇內可嘗試不同的參數,比如原先模擬的溫度較高就先試試低溫、原本加的電場較大就把電場改小點,等等。還找不到原因就把能去掉的設置一點點去掉,比如用了控壓就把控壓去掉、用了凍結就把凍結去掉(或改用位置限制基本等效地實現)、用了約束就把約束去掉、用了電場/pull/外力設置就將它們去掉、用了能量-壓力色散校正就把它去掉,等等,修改設置后重新跑,反復對比找原因,看是否是哪些設置造成的問題。
動力學步長太大是導致動力學不穩定的常見因素,如果原先用的2 fs步長結合氫有關的鍵的約束,應當改成1 fs乃至0.5 fs不加約束再試,如果此時能順利跑下去,等跑一段時間體系弛豫了、相對平穩了,續跑時候可再切換回2 fs步長+約束照常跑。順帶一提,用約束時絕對不要用某些誤人子弟的教程里的constraints = all-bonds。
產生初速度、控溫方式是特別要注意的,溫度和原子速度分布關系密切,若有原子速度過大可能造成崩潰。如果你用了gen_vel = yes結合默認的gen-temp = 300要求程序產生對應300 K時Maxwell分布的原子初速度,建議把溫度設低點,比如50或100 K,否則可能有個別原子初速度太大造成崩潰。關于控溫,我一般都建議設置退火,從較低溫度(0或100K)經過100~200 ps線性、平緩地升溫到實際研究的溫度。如果不設退火,而且熱浴時間常數tau-t比較小的話,可能由于溫度變化太快導致出現有的原子速度瞬間過大而崩潰。
檢查結構文件:從結構文件(給grompp用的gro或pdb)角度來說,一定要注意結構文件里的分子順序和top文件里的[ molecules ]里記錄的順序必須一致,且結構文件里各個分子里原子的順序和拓撲文件里相應的[ moleculetype ]中的[ atoms ]里的順序必須一致,否則用的參數就亂套了,就算不崩潰結果也毫無意義。
還要確保動力學初始結構中沒有明顯不合理的原子間接觸。如果之前沒做能量極小化的話記得先做能量極小化以消除不合理的接觸、過大的斥力,否則動力學一開始由于有的原子受力太大導致之后速度太大,很容易崩潰。如果之前做過能量極小化,而mdp里emtol設得較大,比如500,可以設得更小點(如100或50),讓能量極小化做得更充分一些再試。也注意nsteps的設置,這個控制的是能量極小化的步數上限,如果設得不夠大的話,可能斥力還沒充分釋放時能量極小化任務就自動停掉了。
檢查拓撲文件:從拓撲文件角度來說,要注意檢查[ bonds ]、[ angles ]、[ dihedrals ]、[ pairs ]字段,是否該有的項都有,而且里面的參數都是合理的。把分子結構連同原子序號顯示出來,對照著一個一個看,基于分子力場和GROMACS拓撲文件的基本知識來判斷合理性并不困難。如果參數是你自己照著文獻里的值填進去的,或者從其它力場/程序里挪過來的,注意原本的地方用的力場函數形式是否和當前相同、單位轉換是否弄對了。還要注意檢查[ atoms ]里的原子電荷是否都是合理的。如果里面涉及到了你額外補充的原子類型,也檢查[ atomtypes ]里的范德華參數的合理性,注意函數形式、單位等問題。
注意看grompp的提示:記得在用grompp產生tpr文件時,一定要仔細看屏幕上的NOTE和Warning,理解是什么含義,能忽略的就忽略(出現N個warning時可以grompp加-maxwarn N來忽略),如果是可能導致動力學出問題的因素則要想辦法解決,不能簡單地無視。
仔細觀察軌跡:通過在VMD里看動力學軌跡,經常從肉眼上就能發現結構變化不合理的地方,進而判斷出崩潰的原因。比如拓撲文件里少了個bond項,自然會導致動力學過程中該維持住的地方散架;再比如某個鍵角參數的平衡值或力常數嚴重不合理,很容易造成結構嚴重扭曲,這顯然從軌跡里是能看出來的。如果動力學模擬剛一開始就崩潰,為了便于找原因,可以把nstxout設為1(或很小的值),這樣每一步的結構都會輸出到軌跡文件里便于細致分析結構變化,這對于一些特殊情況很有用。比如曾經我用某個程序產生帶磺酸基的分子的拓撲文件,看似拓撲文件沒任何毛病,但跑起來就是會崩。后來一幀幀地觀看模擬開始后的結構變化,發現是磺酸基上的氫在一瞬間跑到距離氧很近的地方去了。分析發現是由于這個基團上的氫帶的電荷很正,氧的電荷很負,它們彼此間有1-4靜電作用項因此有很強的靜電吸引,再加上氫的質量很小、很容易運動,導致了上述看到的問題。像這種情況,可以把1-4作用項去了(從[ pairs ]中刪除對應的),也可以把氫改成質量更大的重氫(在[ atoms ]里改原子質量)減緩其運動,或者用很小的模擬步長。顯然,從拓撲文件和參數中找出原因必須對拓撲文件里的各種項、力場的函數形式都理解得很清楚才行。
當懷疑一個分子的拓撲文件有問題時,應當先對對分子在真空下進行模擬,只有一個分子時比較容易觀看;而如果帶著其它分子一起跑,或者模擬的是某種分子的凝聚相,顯然就不方便觀看結構了。但也有時候在真空下模擬時沒問題,僅當其它分子存在,因而涉及到分子間相互作用時問題才會暴露出來。
簡化體系:如果你模擬的是復雜體系,遇到動力學過程報錯時應想到到不斷簡化體系進行測試,體系越簡單越容易找出報錯的原因。比如,模擬一個蛋白質+配體+納米管+水的體系遇到報錯,對這樣復雜的體系顯然很難直接排查出問題出在哪里,因此應當分別嘗試蛋白質+水、納米管+水、配體+水能否模擬。如果發現只有配體+水不能模擬,那再看看真空中跑一個配體分子能不能模擬。如果還不能,那自然就應當著重分析是不是配體的拓撲文件存在問題。這種把復雜問題簡化是遇到計算報錯問題最基本的邏輯。
并行計算方式:mdrun通常默認會使用多個thread-MPI來并行計算,此時會使用域分解策略進行計算。但對于個別情況,可能由于不當的域分解策略造成崩潰。可以給mdrun加上-ntmpi N和-ntomp M,分別指定thread-MPI并行進程數和每個下屬的OpenMP并行線程數,反復調節N和M。對于非GPU加速的情況,N和M的乘積應等于CPU的物理核心數以最大程度利用CPU運算性能。遇到崩潰問題,建議先嘗試只寫-ntmpi 1,這樣的話相當于不會利用域分解,可以由此判斷是否是域分解造成的崩潰問題。
該說的上面都說了,沒什么可補充的了。以后問我為什么GROMACS動力學崩潰時,不要就說個崩潰了,沒有信息量的問題不可能回答。先盡自己最大可能按上面說的自己排查、分析、解決,如果花費兩天時間還解決不了,屆時把問題描述詳細并發到計算化學公社論壇(http://bbs.keinsci.com)的分子模擬板塊,必須上傳mdp、拓撲和結構文件(別就傳個mdp文件,別人都不知道你跑的結構是什么。而且前面說了,mdp里設置不合理僅僅是可能導致出現問題的眾多原因之中的一個),我總會看到,在有空的時候會回答。
FAQ 9:為什么Sobtop調用OpenBabel計算MMFF94原子電荷時失敗?
首先要知道,MMFF94是面向普通有機體系構造的力場,此力場定義的計算原子電荷的方法也只適合普通有機體系。如果你的體系里有比如金屬、硼等普通有機體系里不存在的元素,或者你的體系里有一些很特殊的成鍵方式,顯然就沒法成功計算MMFF94原子電荷了。而且就算能計算出來,結果也可能不靠譜(至少應當結合理論化學常識判斷一下)。
如果安裝OpenBabel的方式不當,也可能由于OpenBabel的功能不正常而得不到MMFF94原子電荷。切勿用一些古靈精怪的方法安裝OpenBabel。如果你用的是Windows,應當去OpenBabel官網上下載Windows版安裝包并安裝,不建議用conda等特殊方式安裝(除非你安裝完了后通過手動用命令行方式使用OpenBabel時能確保功能是正常的),而且安裝時不要裝到默認的帶空格的Program Files目錄下。
mol2文件格式不標準也可能導致OpenBabel沒法正常載入而失敗,應當自己檢查mol2文件內容。
FAQ 10:使用Sobtop必須提供Hessian矩陣么?
當 然 不 是!!!我很難理解為什么經常有人有這種誤解。明明本頁面里的一些例子諸如“產生苯甲酸甲酯的GROMACS的拓撲和gro文件:參數完全基于GAFF的”就根本沒有利用用戶提供的Hessian矩陣。看教程一定要一個字一個字認真看完整,絕對不要斷章取義。教程里寫得很清楚,Sobtop有bonded_param.dat力場庫文件,如果你要求Sobtop用的參數全都直接來自于力場庫文件,當然就用不著提供Hessian。僅當你要求Sobtop基于Hessian計算力常數的情況才必須提供振動分析任務產生的Hessian。
FAQ 11:什么時候才需要Sobtop基于Hessian算力常數?
對那些適合諧振勢描述的bonded項(鍵長、鍵角、不可周期性旋轉的剛性二面角),而且GAFF力場里沒有適合參數,而且文獻里也找不到合適的參數的情況,才有必要自己算力常數。最典型的就是諸如《使用Sobtop超級方便地創建二茂鐵的GROMACS的拓撲文件》 (http://www.shanxitv.org/635)所示的過渡金屬配合物,通常涉及過渡金屬配位鍵的相關力場項都適合當做剛性的,因而可以用諧振勢描述,但GAFF里卻沒有適用的參數,文獻里也很難找到能用的參數,所以才明顯有必要讓Sobtop基于Hessian計算它們的力常數并寫到拓撲文件里。
我之前看到有人計算柔性有機分子(比如多巴胺)居然也讓Sobtop計算力常數,這明顯是費力不討好。這種分子有柔性的可旋轉的二面角,若計算力常數而把這樣的二面角用諧振勢來描述,二面角還怎么旋轉?分子的柔性特征還怎么體現?而且GAFF對于普通有機分子就已經描述得挺好,直接讓Sobtop按照苯甲酸甲酯的例子指認GAFF的參數就完了,這還免得你去用量子化學程序事先做結構優化再做振動分析。
FAQ 12:基于Hessian算的力常數產生周期性體系的拓撲文件時,振動分析必須用超胞么?
做周期性體系動力學模擬時,往往涉及大量的重復單元,甚至能達到好幾千原子,而直接對這么大體系用CP2K做振動分析得到力常數耗時極高,或者根本不可能算得動。
對于這種情況,你可以取一個合適大小的晶胞(如果單胞就已經夠大了,比如MOF-5,用單胞即可。如果單胞很小,必須先平移復制到足夠大,以確保得到的力常數靠譜,參考FeCl2的例子的做法),對它做振動分析,并用Sobtop得到體系中涉及的各種力常數。然后將這些力常數連同平衡幾何參數寫入到bonded_param.dat里,并且在assign_AT.dat里定義原子類型判斷規則。之后當前體系的任意大的超胞結構都可以迅用Sobtop速產生拓撲文件了,產生的時候讓Sobtop根據assign_AT.dat的規則判斷原子類型,并基于原子類型直接從力場庫文件里取合適的參數。
對原子電荷也是如此,絕對不是說實際要模擬多大的超胞就必須對多大的超胞體系用CP2K計算。也是取一個算得動大小的模型體系,算出來原子電荷后寫入assign_AT.dat的原子類型判斷規則里即可。之后Sobtop按照此文件里的規則指認原子類型時,原子電荷也就一并指認了。
FAQ 13:為什么Sobtop載入mol2文件時元素不識別?為什么識別出的元素和實際不符?
這都是因為你提供的mol2文件不規矩所致。給Sobtop用的mol2文件不能太不規矩,然而有些來源、有些程序產生的mol2文件相當不規矩。用文本編輯器打開.mol2文件,可看到在@<TRIPOS>ATOM字段下面記錄了原子的信息,比如
7 H7 -4.0563 -0.4333 -0.0000 H
第一個數字是原子序號。H7是這個原子的原子名。然后是原子的XYZ坐標。再往后是atom type(這和實際Sobtop指認的原子類型是兩碼事)。當前atom type是H,因此此原子會被Sobtop解析為是氫原子。在mol2標準格式中,atom type也可以是諸如N.am、S.O2、O.3等形式,可見共性都是圓點前面是元素名,因此Sobtop如果發現atom type中有圓點的時候會根據圓點前的字符判斷元素。
根據以上規則,通過檢查你的mol2文件,自然就能明白為什么有的原子沒有被Sobtop正確識別了。原子信息不規矩的話自然需要自己手動修改,或者換其它程序產生mol2文件。
FAQ 14:為什么sobtop產生的周期性體系的拓撲文件里面有平衡鍵長很大(如幾nm)的bond項?
這是因為你用GaussView產生的mol2當做sobtop一開始的輸入文件,而且要求sobtop自動猜或者基于Hessian矩陣計算bonded參數。即便GaussView載入的是cif這樣含晶胞信息的文件,它產生的mol2文件也不包含晶胞信息,因此sobtop會以為當前體系是孤立體系,那些本來是跨越晶胞的鍵就成了孤立體系內長得離譜的鍵了。有兩個解決辦法:
(1)用記錄了晶胞信息的pdb文件作為輸入文件,sobtop會根據原子間距離(考慮周期性)和元素半徑判斷連接關系。此做法的缺點是沒法保證sobtop判斷的連接關系和GaussView圖形界面里看到的鍵完全一致
(2)在mol2文件里寫入晶胞信息,使得sobtop把此體系視為周期性體系。在mol2文件末尾加入一行@<TRIPOS>CRYSIN,在下一行寫晶胞的a、b、c三個邊長(埃)以及alpha、beta、gamma夾角(度),每個值之間以逗號分隔。例如:
@<TRIPOS>CRYSIN
3.785,3.785,9.514,90,90,90
FAQ 15:為什么sobtop對巨大體系指認原子類型花費時間極長?
如果體系非常大,比如幾千個原子,指認原子類型那一步不要選“Assign GAFF atom types...”或“Assign AMBER atom types...”,否則程序通過分析化學環境來判斷GAFF或AMBER原子類型花的時間極多,甚至好幾個小時也完成不了。這種情況應當用其它方式指認原子類型,比如按照AT_assign.dat文件定義的規則指認原子類型,或者手動通過輸入原子序號指認原子類型。
FAQ 16:我對很大的固體表面用sobtop產生拓撲文件,但生成二面角的過程耗時極高,或者產生出的itp文件尺寸極大,怎么解決?
要求不產生二面角項就行了,即把sobtop.ini里的iskipgendih設為1。很大的體系二面角項巨多,考慮二面角項的話必定導致產生拓撲文件耗時非常高,以及產生出的itp文件巨大。對于研究固體表面類型的體系通常二面角項沒有什么用處,因為靠bond和angle項就能維持住結構。固體表面大多是呈高度剛性的,或者其可變形特征對于實際研究的問題可忽略不計,這種情況就連bond和angle項都可以不必考慮(如果手動從itp里刪掉的話可以減小itp的體積),動力學模擬時把原子坐標凍結或者加上位置限制勢固定就夠了。
TODO list
Supporting fitting dihedral parameters based on scanning curve obtained by quantum chemistry calculation.
Designing an convenient and ideal protocol to prepare topology file of metalloprotein.
Designing an convenient protocol to prepare .rtp file of non-standard residue.
Facilitating assignment of RESP charges for polymers
Update History
Version 1.0(dev3.1): First release: 2022-May-18, Last update: 2022-Aug-9
Version 1.0(dev3): 2022-Mar-26
Version 1.0(dev2): 2022-Mar-14
Version 1.0(dev): 2022-Feb-16