gromacs決定力場參數的方式

Q:
43a1力場中，有的鍵長、鍵角、二面角都是空白的，只有原子的信息，這種情況是不是constraint啊？

比如這樣：

   51    91     2         
   52    53     2    gb_4

A:
[bonds]里面指定了原子間連接關系，不管后面力場參數是不是空的，都說明原子是相連的。此時，如果在mdp里面設了constraint=all-bonds，所有[bonds]里面的鍵就會被約束住(不是按照距離遠近判斷，而是看[bonds]里面怎么連接的)。比如51 91 2后面雖然沒寫力場參數，此時也會被約束住。如果寫了力場參數，則力場參數視為無效，也按約束來處理，用gmxdump -s察看.tpr，看不到它們處于成鍵項(比如810 type=366 (G96BONDS) 796 797)，而會看到處于約束項(比如672 type=429 (CONSTR) 796 797)。也就是說constraint的優先級最高。


不同的力場決定成鍵參數的方式不同。對于設constraint=none的情況下有兩種情況：
對于GROMOS96力場來說，用什么參數只取決于grompp時讀入的所要研究體系的.top(或者被include的.itp)，里面寫了哪些力場參數就發揮哪些效果，沒有寫的力場參數就當作沒有(被設為0)。比如[bonds]里寫51 91 2，雖然指定了連接關系，但沒寫力場參數，又沒有設約束，等跑起來之后就可能會逐漸遠離。

對于OPLSAA力場來說，哪些鍵用哪個力場參數是在grompp時根據原子類型來決定的(類似amber中leap程序的方式)，而不像GROMOS96那樣在.rtp里面事先定義。OPLSAA力場在所要研究體系的.top(或者被include的.itp)里面只定義有鍵結項(本文中鍵結項泛指鍵長、鍵角、二面角)并不直接定義參數，例如：
[ bonds ]
   1     2     1
......

[ pairs ]
    1     9     1
......

[ angles ]
    2     1     3     1
......

[ dihedrals ]
    2     1     4     5     3
......
這和GROMOS96力場處理方式截然不同。雖然[bonds]里1 2 1沒定義參數，但是用gmxdump -s看tpr卻發現連接1和2原子的鍵存在，模擬過程中也保持一定距離。這說明對于OPLSAA力場，grompp從.top(或者被include的.itp)獲知的僅僅是有哪些鍵結項，這些鍵結項涉及的原子會通過開頭[ atoms ]段里面的定義將原子序號轉化為原子的type(如opls_136)，再由ffoplsaanb.itp的定義轉化為bond_type(如CT)，再從ffopsabon.itp里面找到對應參數，放進.tpr里面，在模擬過程中生效。


對于[bonds]里的內容處理方法簡單來說就是，設了constraint=all-bonds，不管什么力場，[bonds]里寫了哪些項就約束哪些。如果constraint=none，用GROMOS96力場時有項有參數的就起作用，有項沒參數就等于沒寫。OPLSAA里面有項沒參數也起作用。

PS：關于自定義constraints的情況可參考《解析gromacs的restraint、constraint和freeze》（http://www.shanxitv.org/10）。

善用gmxdump -s察看.tpr是最可靠的方法。

---

在這里我將問題擴展一下，說一下gromacs決定鍵結參數的機制。實際上OPLSAA力場與GROMOS96力場都是遵循同樣的方法來處理，但是如上所述，決定鍵結參數有很大區別，并非因為程序對它們處理機制的差異，而是力場參數文件的差異。下面都不考慮約束問題。

例如用某個力場，研究的體系的.top里面的一個普通的成鍵定義：
[ bonds ]
   1     2     1    ga_1                //注意第二個1代表函數類型，不是原子號
首先gmx會用C預處理程序(一般就是cpp)將ga_1還原為實際的參數，轉換關系定義在ff????bon.itp里，在grompp時加上-pp，從得到的processed.top中就會看到ga_1被還原了。當然如果ga_1位置直接寫的就是參數，就不必轉換了。這樣1和2之間的鍵參數就定下來了。

另一種情況，雖然定義了成鍵項，但沒寫參數：
[ bonds ]
   1     2     1
那么程序會首先將根據原子序號，在[ atoms ]里面查找到對應的原子類型，比如轉換為CH2 N，然后按照原子類型到ff????bon.itp里面的[ bondtypes ]里面找原子類型相對應的參數，找到了就定下來了，如果沒找到，在grompp的時候就會出Warning提示這個鍵的參數被設為0，相當于不存在。

GROMOS96與OPLSAA定義力場參數的差異就在于它們的ff????bon.itp的內容。
比如GROMOS96的ffG43a1bon.itp里面[ bondtypes ]幾乎沒有什么內容，只是二硫鍵和血紅素用的幾個鍵的定義。所以必須在[ bonds ]里面直接寫出參數(這步由pdb2gmx完成，主要將.rtp里面的相應殘基的鍵/角/二面角項和參數搬過來)，否則在ffG43a1bon.itp的[ bondtypes ]里也不會找到參數，那個成鍵項等于沒有被定義。
而比如OPLSAA的ffoplsaabon.itp，[ bondtypes ]里面包含了全部OPLSAA力場的成鍵項，所以就不必在體系.top的[ bonds ]里面直接寫參數，grompp時在這里一般都能找到對應的。在ffoplsaa.rtp里面定義的除了原子電荷外，主要就是殘基的連接關系，也就是[ bonds ]，沒有參數，在[ angles ]和[ dihedrals ]里面幾乎沒有內容或者根本沒有。這種情況下pdb2gmx的主要作用除了把.rtp的對應內容直接搬到.top里面以外，另一個主要工作就是根據連接關系自動補全角和二面角項，所以.top里面看到的鍵/角/二面角項是全的，如上所述沒有寫參數。

從大體來說，實際上GROMOS96力場是有些特殊的，它直接在rtp中事先定義好了全部鍵結項和參數，top也包含全部鍵結項和參數，grompp調用的參數直接取自top。而OPLSAA則算是比較普通的情況，rtp只定義連接關系，在top中更進一步把全部的鍵結項表達出來，該用什么參數等grompp時再去從力場參數庫中搜。amber的leap也是如此，載入結構載入的只是連接關系，各個鍵結項用的參數在saveamberparm的時候對應去搜。

決定angle、dihedral的參數與決定bond的方式相同。另外如果比如[ bonds ]里面已經寫了參數，同時在[ bondtypes ]里面也能找到原子類型對應的參數，則優先使用在[ bonds ]里面直接寫的。

這里再說一下對二硫鍵等特殊鍵的處理方式。pdb2gmx讀取結構時，也載入specbond.dat，specbond.dat內如如下
6
CYS    SG    1    CYS    SG    1    0.2    CYS2    CYS2
......
這就說明如果發現在同一個鏈上有兩個CYS上的SG間的距離在0.2nm附近(注意不是僅僅指的在0.2nm以內)，就認為這兩個SG成鍵，并且殘基名由CYS變成CYS2。這時這個半胱氨酸的參數就用的.rtp里面CYS2殘基的參數，加氫時也使用.hdb里面CYS2殘基的加氫方式而非默認的CYSH，也就是硫上面不加氫。在.top里面可以看到出現了這兩個硫原子之間的成鍵項以及相關的角、二面角項，但是沒有參數(無論是GROMOS還是OPLSAA)，這就說明要從比如[ bondtypes ]里面讀對應的參數，這就是為什么比如ffG43a1bon.itp的[ bondtypes ]里面有S S 2 gb_33這樣的二硫鍵參數定義。硫原子的原子名是SG但是原子類型是S，所以二硫鍵和這個參數項對應。對于血紅素的一些特殊成鍵也是這么處理的。
也可以自行添加、修改specbond.dat的內容來處理類似的特殊情況。

---

這里也說說非鍵作用和1-4作用，下文所指的非鍵作用不包括1-4作用。

非鍵作用：

對于VDW參數的確定，GROMOS96力場在ffG43a1nb.itp里面的[ nonbond_params ]中提供了各種各樣原子類型的組合，計算i、j兩個原子間的VDW作用時會根據原子類型對應獲得C(ij)參數(包括C6和C12項)。若計算兩個原子間VDW作用時，發現[ nonbond_params ]里面沒有對應的，就通過組合方法確定，從[ atomtypes ]中得到原子i的C(i)與原子j的C(j)，由于ffG43a1.itp里面comb-rule設的是1，會根據(C(i)*C(j))^1/2公式來獲得C(ij)。由于OPLSAA并沒有[ nonbond_params ]段，故所有C(ij)項都是如上算出來的(OPLSAA的comb-rule是3，詳見手冊)。

對于靜電作用，只依賴于原子電荷，沒什么特殊的。

要注意，在拓撲文件中的[ moleculetype ]段中，nrexcl決定了相隔多少個鍵以內不計算非鍵相互作用，例如nrexcl=3時，與某個原子相隔5個及以內的原子間都不計算非鍵作用。一般為3，即表明直接bonded的兩個原子、構成angle項兩端的兩個原子，以及構成二面角項兩端的兩個原子都不計算非鍵作用。


1-4作用：

1-4作用項雖然看起來也屬于非鍵作用，常被稱為1-4非鍵作用，但gromacs對它的處理與非鍵作用完全不同。1-4作用包括庫侖1-4和1-4VDW相互作用，哪些原子間有1-4相互作用都定義在體系.top里面[ pairs ]，里面每一項就是一個1-4作用項，如果沒寫就說明這兩個原子間不擁有1-4相互作用（即便從成鍵關系上屬于1-4關系的原子也不算）。如果兩個原子包含在1-4作用項中，則不計算它們的非鍵作用，相當于從非鍵列表中排除了。nrexcl對1-4作用完全無效，根據nrexcl的設定，即便某兩個原子間不計算非鍵作用，若存在于1-4作用項中，也會照常計算1-4作用。即nrexcl只管非鍵作用。

1-4作用項中的庫侖1-4作用，就是根據1-4原子間原子電荷照常計算靜電作用然后乘以FudgeQQ(在比如ffoplsaa.itp的[ defaults ]中定義)。

1-4作用項中的1-4 VDW作用，所用C(ij)參數與非鍵VDW中的不同，需要額外的參數，但是GROMOS96和OPLSAA的體系.top在[ pairs ]里默認都沒寫參數，這是因為：

GROMOS96力場，明確定義了各種原子類型組合的1-4VDW參數，全部儲存在比如ffG43a1nb.itp里面的[ pairtypes ]當中，體系.top里面[ pairs ]當中定義的每一個1-4VDW項的參數都從這里對照原子類型自動提取。在ffG43a1.itp當中會看到gen-pairs設為了no，意思是如果讀不到對應的1-4VDW參數，就出現warning并將參數用0值代替。

OPLSAA力場，各種原子類型之間的1-4VDW作用都是普通VDW作用乘以FudgeLJ因子0.5計算得到，也就是進行削弱(GROMOS96定義的1-4VDW參數也是被削弱的)，在ffoplsaa.itp里面可以看到gen-pairs被設為yes，意思是讀不到的1-4VDW參數都通過上述方法計算生成。實際上ffoplsaanb.itp里面根本就沒有[ pairtypes ]，所以如果不做人為修改，1-4VDW參數都是通過計算生成。

如果直接在[ pairs ]當中的項的后面給出參數，這些項就直接用給的參數當作1-4VDW的參數。如果直接給的參數為0，或者沒直接給參數、在ff????nb.itp里的[ pairtypes ]也沒有對應項、而且把gen-pairs設為了no，則參數值會默認為0，這兩種情況下1-4VDW作用效果皆為0。但此時1-4靜電力仍會照常計算，除非把FudgeQQ也設為0，則徹底不計算1-4作用了。

不計算1-4作用實際上有兩種情況：一種是不寫[ pairs ]的情況，1-4作用完全忽視掉，程序根本不算。另一種是寫了[ pairs ]，程序也算了，但是因為上面一段所述的參數和設定原因，導致算出來的1-4作用為0，相當于沒算，這種情況下仍然會花費計算時間，在mdrun的結尾會看到1-4作用計算的開銷。


PS: 基于GROMOS87的ffgmx力場，決定VDW參數方式同GROMOS96，決定鍵結參數方式同OPLSAA(但OPLSAA多了一步原子名的轉換過程，即type到bond_type)。