Gromacs的pull code介紹(3.3.x+4.0.x)

Gromacs 3.3.x:

Pull Code有三種類型：AFM pulling,Umbrella sampling,constraint
AFM是給原子一個外力，實際上是一個彈簧點在恒速遠離，彈簧點拉動一組原子，其力取決于彈簧點與這組原子質心的距離和自己定義的力常數。類似amber的jar=1。
比如A-------B=C   A是參考組，B是pull組，C是彈簧點。模擬開始后AC距離根據constraint_rate1的設定恒速增加，C的移動拉動B，B受C的力等于BC的距離乘以afm_k1。
Constraint是保持兩組原子質心間距離不變(用SHAKE算法實現)，監視它們之間的限制力
Umbrella是給一組原子的質心加一個諧振勢
輸入文件有-pi pull.ppa里面包含pull的控制代碼。-pn index.ndx指定哪寫原子組成哪些組。
輸出文件有-po pullout.ppa是形式化的輸入文件的拷貝，也就是各種參數都寫齊了的輸入的.ppa，類似mdout.mdp。-pdo pull.pdo記錄了pull的過程中的力、位置等的變化信息

定義組的方法：
com 每一步計算一次設定的組內原子的質心
com_t0 每一步計算一次設定的組內原子的質心，但是對原子穿越盒子跑到另一側導致質心位置突然變動這種現象進行校正
dynamic 只有組內的在某個柱形空間范圍內的才被當回事
dynamic_t0 同com_t0與com的關系，做校正而已

一般參數：
verbose = 默認為no，如果設yes，會往.pdo里面多寫一些東西，并通過stderr多輸出一些東西。
Skip steps = pdo里面每多少步輸出一次結果，默認為1
runtype = 可設afm,contraint,umbrella
ngroups = 總共有多少組（不含參考組），最大可設20個
group_1 = LT 設定第一個組
group_2 =
....
reference_group = LTWD 設定參考組。定義的那些group_x與這個組的距離會被計算。比如得到兩個粒子之間的constraint force，則一個定義為group_1一個定義為這個組。如果用AFM或umbrella并且不定義這個組，則外力是從絕對坐標而來。
如果用了參考組，則pull組都是相對于參考組的坐標的，而體系質心不會因為pull了某個東西而平移，字面上看是參考組不動，實際上與pull組相反方向走。如果不定義參考組，pull組在絕對坐標中pull，如果此時設了comm_mode = none，就會導致體系質心也發生移動，若開了消平動，質心基本還是原處。
weights_1 = pull組1的權重
weights_2 = pull組2的權重
reference_weights = 參考組的權重，權重對決定質心位置和相關原子受力有影響，見手冊
reftype = 參考組的類型，可以是com,com_t0,dynamic,dynamic_t0
reflag = 1 參考組的位置由許多步取平均得到
pulldim = N N Y 往哪個方向pull，這里是Z方向。默認都是Y。設了幾個Y，在.pdo里面就會有幾個維度，相對于一個Y時除了時間以外就有幾倍列數的數據。不管afm_dir1怎么設的，這里是N的方向就被忽略。

Dynamic參考組選項：
r = 0 用的是dynamic/dynamic_t0類型參考組時，這個設柱形的半徑
rc = 0 在r至rc的范圍內的粒子，權重為1至0。建議r=1，rc=1.5
update = 1 多少步更新一次dynamic參考組

當runtype=constraint時的參數：
constraint_direction = 限制的方向，默認為0.0 0.0 0.0也就是沒有方向，因此兩個組之間的距離被限制住。
constraint_distance1/constraint_distance2 = 最初的限制距離，默認是0，也就是根據初始坐標決定距離
constraint_rate1 constraint_rate2 = 距離限制的變化速率，為nm/ps，默認為0
constraint_tolerance = 距離限制的容忍限，也就是SHAKE的tolerance，默認1e-6

當runtype=afm時的參數：
例：不設參考組，afm_k1 = 1 -1.0 0，pulldim = N Y N，afm_rate1 = 0.5，不設afm_init1，則pull組會一直朝自己的y的負方向走，不主動向x方向走即便afm_dir1設了這個方向，雖然彈簧點的絕對坐標應該是0 -0.1 0，0 -0.2 0，0 -0.3 0這樣變化從而也應該帶動pull組的x和z接近0，但因為這兩個方向都是N所以不會。但如果pulldim都設Y，則pull不僅y隨彈簧點的y增大而增大，由于彈簧點的x和z一直都是0，x和z也逐漸趨近于0，pull組會追著從絕對坐標0點開始走的彈簧點。
afm_rate1 afm_rate2 = 對于每個組的彈簧移動的速率(nm/ps)
afm_k1 afm_k2 = 對于每個組的彈簧力常數(KJmol^-1*nm^-2)
afm_dir1 afm_dir2 = 彈簧點移動的向量（相對于彈簧點當前位置），如果設了參考組，則參考組動時，彈簧點也隨之動，保持相對位置不變。
afm_init1 afm_init2 = 描述彈簧最初位置相對于參考組的向量。如果想讓開始時的力為0，則彈簧初始位置必須設為pull組相對于參考組的值。

當runtype=umbrella時的參數：
umbrella給pull組施加諧振勢來維持pull組相對于參考組的位置。如果參考組移動了，pull組也會隨之移動。如果沒有定義參考組，pull組就處于絕對坐標中，參考組被設為[0 0 0]。
例如：不設參考組，pos1設0.0 5.0 0.0，當.pdo中的距離很小后，檢查pull組的y值，應當很接近于5.0。相當于把pull組往某個坐標點拉。但是不能設得過大、離譜，否則預計與結果不符甚至被無視了。
例如：設參考組，pos1設0.0 5.0 0.0，當.pdo中的距離很小后，檢查pull組的y值和參考組的y值，其差應當很接近于5.0。動的時候參考組和pull組一起在對稱地動，并非只動pull組。
k1 k2 = 指定每個pull組的力常數(kJ*mol^-1*nm^-2)
pos1 pos2 = 指定pull組被限制到的位置，此值是相對于參考組的。


注：總的來說，umbrella和afm實際上是一樣的，區別只是umbrella的彈簧點是不相對移動的（相當于afm_rate1=0，沒有afm_dir1選項），afm的彈簧點是相對移動的。
pull組始終是向著彈簧點坐標受力。pulldim控制這個力的哪些分量被表達出來。彈簧點的t=0位置(原始位置)都是這兩項的矢量和：參考組(reference_group)+相對參考組的初始位置(afm_init1或pos1)。彈簧點是從原始位置開始根據afm_dir移動。如果參考組不是固定的坐標而是一個組，則這個組的位置變化會加到彈簧點上。
==========
輸出結果
時間都是ps，位置都是nm
對于Constraint，如10.5 -81.2，分別代表時間和限制力，負數說明pull組與參考組相互吸引。
對于AFM，如10.5 3.5 3.63 3.6，分別代表時間，參考組的位置，pull組的位置，彈簧的位置。只有pulldim設為Y的維度的結果才會被輸出。pull組與參考組距離之差的變化與afm_rate相等
對于Umbrella，如10.5 0.002，分別代表時間和pull組相對于其限定的位置。


其它：如果跑起來沒多久就出現一堆東西，提示write pdb之類，說明體系崩了，就把力常數設小點。并且檢查是否拉伸過程合理，比如拉伸過程中是否出現過度碰撞、擠壓。

Gromacs 4:

gromacs4中，pull代碼進行了重寫，參數定義從pull文件改為了.mdp中定義。拉伸的三種方式也有所變化，把afm和umbrella合并為umbrella（因為本質一樣），新增加了常力拉伸(constant_force)，即外加力大小一直不變。
pull_vec1和pull_init1的格式是例如：1.0 2.0 3.0。
主要參數對應關系：3.3 vs 4.0
runtype=pull
pulldim=pull_dim
reference_group=pull_group0 (下文中參考組指的就是pull_group0)
group_1=pull_group1  group_2=pull_group2......
pos1/afm_init1=pull_init1 (與pull_geometry有關)
afm_dir1=pull_vec1 (與pull_geometry有關)
afm_rate1=pull_rate1
ngroups=pull_ngroups (二者數目皆不包含group0)
k1/afm_k1=pull_k1

4.0多了個pull_geometry選項。包含position、distance、direction、cylinder四種模式，這里只討論前三種。
======================
使用pull=umbrella結合這三種模式的情況：

position模式(重點)：等價于3.3.3的afm，即彈簧點位置是參考組位置+pull_init+time*pull_rate*pull_vec，朝著彈簧點拉。如果將pull_rate1設為0，則相當于彈簧點位置不變，就等價于3.3.3版本中的umbrella。

distance模式：始終順著連接參考組與pull組向量的方向拉。用這種模式時彈簧點只會在參考組與pull組連線方向上移動。pull_init1只能設一個數而不是坐標，這個數就是彈簧點t=0時在這個連線方向上與參考組的距離。pull_rate1是彈簧點在連線方向上的移動速度，逐漸離近還是逐漸遠離可以通過此值的正負號來控制。pull_vec1對彈簧點運動方向無任何影響。

direction模式：定義一個過pull組在t=0時坐標的向量，pull_vec1方向。pull_init只能設一個數而不是坐標。設參考組在這個向量上投影點是P，則彈簧點的位置是P+pull_init*uni(pull_vec)+time*pull_rate1*pull_vec1。uni()代表求單位向量。

上述都是假設參考組是絕對坐標或者固定住參考組，如果參考組在拉伸過程中位置變化，彈簧點位置也相應變化。
======================
使用pull=constraint結合這三種模式的情況和umbrella一樣，區別只是將pull組位置用約束算法固定在彈簧點位置上，而不是被彈簧點拉著走。比較常用的是distance模式，這個模式可使參考組與pull組之間的距離固定保持不變(即pull_init1)或者刻意讓其逐漸改變(即pull_init1+time*pull_rate1)，而它們之間的相對角度方位可隨意變化。
======================
使用pull=constant_force結合這三種模式的情況：
注意，拉力大小始終不變，為pull_k1。無論哪種模式皆沒有彈簧點的概念，故pull_rate1、pull_init1對結果皆無影響。

position模式：pull組由t=0時的位置向-1,-1,-1方向拉（不是指朝著(-1,-1,-1)坐標）。pull_group0、pull_vec1的設定不對結果產生影響。這種模式不適合constant_force。

distance模式：如果參考組設的是一個組，則無論何時，pull組始終朝著參考組方向拉，若拉伸過程中經過了參考組，則拉力方向反過來。如果參考組是絕對坐標(0,0,0)，則拉力方向始終是t=0時pull組朝著參考組的方向，即便pull組可能已經經過了參考組。pull_vec1的設定不對結果產生影響。

direction模式(重點)：將pull組（將其假想為坐標原點）向pull_vec1的反方向拉。pull_group0的設定不對結果產生影響。

======================
上述六種組合中，假定參考組位置固定，若不進行位置固定，則group0與group1是對稱運動的，質心位置不變（包括constant_force+direction這樣不直接涉及到group0的模式也是如此），兩個組相當于等價。

pull_dim控制哪個方向上的力真正生效，上述假定此項設為了Y Y Y。無論是哪種模式、怎么設定拉伸拉伸，最終在pull組上真正生效的力是pull_dim允許方向上的分量。但如果是constant_force，則在pull_dim方向上的總受力仍是pull_k1。

總結：前面說的三種pull_geometry模式，總的來說position本意是相對于參考組位置設定一個點，以其為依據進行拉伸操作。用distance時，拉伸的方向是根據pull組與參考組連線的瞬時方向。direction是直接定義pull組的拉伸方向。
常用的組合模式是umbrella+position可固定勢阱位置；constraint+distance可約束pull組與參考組的距離；constant_force+direction方便地直接將pull組往某個方向拉。

另外，pull_start如果設為yes，則最終的pull_init相當于你設的pull_init加上被pull組與參考組之間質心的距離矢量。



mdrun時利用-pf可以指定輸出的包含pull組的時間和對應時刻受力的xvg文件名（輸出的是外加的pull力，而非這個組的凈受力。另外，若為position模式會輸出受力xyz分量，distance和direction只輸出總受力）。-px可以指定輸出的包含pull組的坐標和對應時刻受力的xvg文件名，第一列是時間，第2、3、4列是pull_group0的xyz絕對坐標，5、6、7列是pull_group1相對于pull_group0的xyz坐標。它們的輸出頻率使用pull_nstfout和pull_nstxout選項控制。

同時pull N個組：
設定方法同pull一個組，寫pull_ngroups=N，增加定義pull_group2、pull_init2、pull_k2、pull_rate2；pull_group3、pull_init3...... 受力-時間 和 坐標-時間的.xvg文件內也會輸出這些組。