解析gromacs的restraint、constraint和freeze

聲明：此文最初發布于Sobereva的百度博客，后在網上被一堆博客胡亂轉載，卻不注明真正來源。此文真正原文在此！鄙視胡亂轉載但是不注明出處的人！

gromacs里面人為固定坐標有三大類方式：restraint、constraint和freeze，有各自的特點。

限制(restraint)

包括位置限制(也就是常說的正式MD前首先要進行的限制性動力學)，角度限制，二面角限制，方向限制，距離限制。限制的特點是可以讓被限制的東西在一定范圍內運動，而非徹底固定住，實際上就是施加一個諧振勢來限制其移動。

位置限制就不必說了，posre.itp里面[ position_restraints ]定義的就是，默認相對于最初位置進行限制。

角度限制包括兩類，一類是兩對兒原子間角度的限制，用[ angle_restraints ]來指定，例如：
[ angle_restraints ]
; i   j    k    l    type   theta0     fc     multiplicity
651 1211 1683 1211    1      67.0     1500         1
說明限制651-1211與1683-1211原子對兒之間的夾角在67度附近，力常數為1500。type無用。
角度限制另一類包括一對兒原子與z軸夾角角度的限制，用[ angle_restraints_z ]來指定，例如：
[ angle_restraints_z ]
; #1 #2  type  theta0    fc     multiplicity
22  45     1     90      500         2
說明限制22與45號原子的連線與z軸垂直，力常數是500，多重性是2，使得90度夾角時候限制勢能最低，0和180度時最高。type無用。

二面角限制使用[ dihedral_restraints ]段落來定義，實際上improper項就是用二面角限制方式限制的。例如
[ dihedral_restraints ]
;   i    j     k    l    type  label  phi  dphi  kfac  power
    5    7     9    15     1      1  180     0     1      2
被限制的是5,7,9,15原子組成的二面角，type總是1，label沒用，phi是參考角，dphi是超過參考角多少度開始使用限制勢，power沒用。kfac乘上mdp中的dihre_fc將作為限制勢力常數。
最后在mdp中加入例如：
dihre               =  simple
dihre_fc            =  100
dihre_tau           =  0.0
nstdihreout         =  50

距離限制在[ distance_restraints ]里面定義，比如
[ distance_restraints ]
; i j type index type low up1 up2 fac
10 16 1 0 1 0.0 0.3 0.4 1.0
type總是1，index是計算的順序，如果幾個項index都一樣，比如10-28和10-46，就同時計算。勢能圖見gmx3.3手冊p60，low,up1,up2分別指圖上的r0,r1,r2，可見原子間距離在low至up1區間內是不受限制的，這種方式可以達成NMR限制。fac是指這個因子乘上mdp中disre_fc作為限制勢力常數。
也可以定義兩個原子在[ bonds ]里用bond type 6，就是個和普通鍵一樣的諧振子勢，但是這兩個原子間被認為沒有鍵連。

應當注意以上限制方式中[]段落應當緊接著寫在被限制分子的.itp后面（或者說對應的[ moleculetype ]后面），這樣程序才知道其中的原子編號指的是哪種分子中的原子。

約束(constraint)

用shake或lincs方法固定住原子之間的相對位置，固定的幾何關系是絕對不變的。分為兩種類型，type 1和type 2的效果其實一樣，但是type 1被認為將這兩個原子鍵連了，而type 2沒有鍵連。一般被鍵連的原子間都不計算相互作用(在[ moleculetype ]里定義nrexcl來控制相隔幾個鍵內的原子間不計算相互作用，也就是相當于被列進[ exclusions ])，比如type 1時一般就不計算這兩個原子間的非鍵作用了，而type 2時仍然計算。

用哪種約束算法用.mdp里的constraint_algorithm設定，默認lincs。.mdp里的比如constraints = all-bonds也是應用這種約束方法，也就是約束住所有[ bonds ]項，原本[ bonds ]的設定就被覆蓋了，即不體現成鍵效果（來回振動），只體現約束效果（距離固定不變）。

自定義約束項可以在拓撲文件中這么寫:
[ constraints ]
1 8 2 0.153    //原子1 原子2 type 約束的距離(0.153nm)。
之后1 2原子距離會固定保持在0.153nm，即便1 2原子在[ bonds ]中定義是成鍵的也是如此。
注意，.mdp里如果寫constraint=none，只是說批量約束hbonds、all-bonds等沒有了，在拓撲文件中自定義的constraint約束項仍然生效。

順便一提，在默認情況下，gromacs的水的結構是被settle算法約束住的，也就是所說的rigid水，settle就是專用于水的約束算法，在spce.itp里面[settles]即是定義，不計算分子內氫、氧彼此之間的相互作用。而如果.mdp里面設define=-DFLEXIBLE，也就是激活了#ifdef FLEXIBLE后面那段，則不用settle算法約束結構，而是照常按照諧振勢的鍵長鍵角項控制水的結構。-DFLEXIBLE開不開對溶質沒什么影響，對計算速度也沒什么影響(有一丁點減慢可忽略)。除了能量最小化外，都不要用-DFLEXIBLE，因為擬合水、力場參數的時候都一般是認為水是剛性(rigid water)情況下進行的(比如SPC水和使用SPC水的GROMOS力場)，即約束住鍵長鍵角，開了FLEXIBLE實際上參數還得再作調整。而在gmx4中，在能量最小化時自動是限制方法(restraint)，就沒必要設-DFLEXIBLE了，從此可忘掉它。.mdp中constraints以及constraint_algorithm設定皆與水的約束無關，除非define=-DFLEXIBLE，否則水都是用settle算法約束住。

凍住(freeze)

徹底凍在最初坐標，一點也不動。但仍然計算與其它原子間的相互作用，所以并不會省時間，除非是寫進[ exclusions ]。

使用時在.mdp里面定義比如：
freezegrps = protein   設蛋白凍住
freezedim = Y Y Y    設凍住的方向，分別對應X Y Z軸。


同時使用freezegrps, constraints, pressure coupling可能會有問題，不要組合使用。

注意絕對坐標與相對坐標的限制方式的不同，比如想保持一個結構的剛性，用位置限制或freeze是原地不動，而通過約束方法對分子內的結構進行約束，結構仍然可以有整體平動和轉動。應根據實際情況選擇。