gmx 拉伸分子动力学及PMF计算

参考资料

前期工作

建模；
添加溶剂和平衡离子；
能量最小化和平衡；

SMD

mdp: pull code

pull的代码基于md而非steep等，以下为Xilock的pull部分code：

;其他代码如run parameters, output parameters, bond parameters, cut off, couple等
;.....
;Pull code
;根据需要确定是否需要freeze住无关部分
;freezegrps      = RUT
;freezedim       = Y Y Y
;或者用-DPOSRES固定住
;define      = -DPOSRES_B
;拉伸，通过速率和模拟时间来控制拉伸距离
pull            = yes
pull-ngroups    = 2
pull-ncoords    = 1
pull-group1-name = Chain_A
pull-group2-name = Chain_B
pull-coord1-type     = umbrella      ; harmonic biasing force
pull-coord1-geometry = distance
pull-coord1-groups = 1 2   ;group2相对于group1做SMD
pull-coord1-dim      = N N Y
pull-coord1-start    =yes
pull_coord1_rate        = 0.01          ; 0.01 nm per ps = 10 nm per ns
pull_coord1_k           = 1000          ; kJ mol^-1 nm^-2
;pull-nstxout    = 0 ;不输出pullx文件
;pull-nstfout    = 0 ;不输出pullf文件

后面做production MD时只需将速度改为0即可，留意freeze和define的修改与使用；
若pull后的构象没有速度，记得打开gen_vel = yes和gen_temp = 310来生成初始速度；

牵引方式与牵引几何的组合

I:pull=umbrella

pull-geometry=distance:始终沿参考组质心与牵引组质心连线的方向进行牵引. 弹簧点只会沿两质心的连线方向移动. 弹簧点t=0时刻在此连线方向上距参考点的距离由pull-coord*-init设定, 弹簧点在连线方向上的移动速率由pull-coord*-rate设定, 逐渐靠近还是逐渐远离可以通过此值的正负号来控制.在这种情况下,pull-coord*-vec的设定对弹簧点的运动方向无任何影响, 也不会对结果产生影响.
pull-geometry=direction:定义矢量pull-coord*-vec, 为t=0时刻过牵引组质心的矢量. 设参考组在这个向量上投影点为P, 则弹簧点在t时刻的位置为P+(pull-coord*-init+tpull-coord*-rate)uni(pull-coord*-vec).其中uni()代表求单位矢量.. 上面几种情况都是假定参考组是绝对坐标或者被固定住, 如果参考组在牵引过程中位置会发生变化, 弹簧点位置也会相应变化, 以维持间距均匀增长.

II:pull=constraint 与pull=umbrella一样, 区别在于使用约束算法将牵引组位置固定在弹簧点位置上, 而不会被弹簧点拉着走.

pull-geometry=distance:可使参考组与牵引组之间的距离保持固定不变(即pull-coord*-init)或者刻意让其逐渐改变(即pull-coord*-init+time*pull-coord*-rate), 而它们之间的相对角度方位可随意变化.

III:pull=constant-force 这种情况下拉力大小始终不变, 由pull-coord*-k设定. 无论使用哪种几何都没有弹簧点的概念, 因此pull-coord*-rate,pull-coord*-init对结果无影响.

pull-geometry=distance:如果设定的参考组为原子组, 则无论何时, 始终朝参考组方向拉伸牵引组. 若牵引过程中越过了参考组, 拉力会反向. 如果参考组为绝对坐标, 则拉力方向始终是t=0时刻牵引组朝着参考组的方向, 即便牵引组可能已经越过了参考组.pull-coord*-vec的设定对弹簧点的运动方向无任何影响, 也不会对结果产生影响.
pull-geometry=direction:将牵引组(将其假想为坐标原点)向pull-coord*-vec的反方向拉.pull-group0-name的设定不会对结果产生影响.

从轨迹中抽取构型

用下述命令处理拉伸轨迹，得到conf0.gro，conf1.gro等文件：

echo 0 | gmx trjconv -s pull.tpr -f pull.xtc -o conf.gro -sep

计算质心距离

使用jalemkul@vt.edu的脚本perl distances.pl

在bash环境下运行该脚本，如git的bash，cmd下可能报错。
使用前注意修改构象个数和group名称。

perl distances.pl

附Xilock修订后的distances.pl:

#!/usr/bin/perl -w

use strict;

# loop g_dist command - measure distance in each frame, write to a file
# Total_Num is the total number of confXXX.gro
our $Total_Num = 250;
for (my $i=0; $i<=$Total_Num; $i++) {
    print "Processing configuration $i...\n";
    # Note1: in Win, replace ' by "
    # Note2: \ is used to transfer; & should be delect sometimes
    # Note3: group name in "" such as group "RUT" or groug with numbe such as group 25
    # Note4: group name should be exclusive in ndx file, do not have repeated group name!!
    # Note5: -oxyz will output distance in x, y, and z anxials respectively, we used z in the following.
    system("gmx distance -s pull.tpr -f conf${i}.gro -n index_distance.ndx -oxyz dist${i}.xvg -select \'com of group \"RUT\" plus com of group \"Lysozyme\"\' &>/dev/null");
}

# write output to single file
open(OUT, ">>summary_distances.dat");

 for (my $j=0; $j<=$Total_Num; $j++) {
    open(IN, "<dist${j}.xvg");
    my @array = <IN>;

    my $distance;

    foreach $_ (@array) {
        if ($_ =~ /[#@]/) {
#             do nothing, it's a comment or formatting line
        } else {
            my @line = split(" ", $_);
#            third one is in z anxial
            $distance = $line[3];
        }
    }

    close(IN);
    print OUT "$j\t$distance\n";
}

close(OUT);

# clean up	
print "Cleaning up...\n";

 for (my $k=0; $k<=$Total_Num; $k++) {
    unlink "dist${k}.xvg";
}

exit;

注：主要有三处要修改的内容

总构象数Total_Num；
gmx distance代码；
要输出的距离方向，如$line[1]为x轴，$line[3]为z轴；

提取固定间距构象并生成gmx命令的sh文件

使用Mike Harms的脚本setupUmbrella.py

python setupUmbrella.py summary_distances.dat 0.2 run-umbrella.sh > caught-output.txt

注意脚本中要修改的几个地方：

line150之后有三处print函数后要加括号()，如最后一行改为：print (““.join(out));
若使用python3，则将line61的keys.sort()改为:sorted(keys)
若只输出最后一个configarition则可能是排序有问题，将sorted(keys)改为keys = sorted(keys)即可。（如在python3.6中）

run-umbrella.sh代码示例：

#!/bin/bash

# NVT
gmx grompp -f mdp/nvt_10ns.mdp -c confXXX.gro -p topol.top -n index_ions.ndx -o nvt_10ns_XXX.tpr
gmx mdrun -v -deffnm nvt_10ns_XXX -pf pullf-XXX.xvg -px pullx-XXX.xvg

得到一堆sh文件后，可以用cat *.sh > MergeSh.bsh得到组合的批处理脚本，去掉多余的部分即可。
或者在bsh里使用循环： merge.bsh

#!/bin/bash
for i in 6 14 20 27 34
do
	gmx grompp -f mdp/nvt_10ns.mdp -c conf$i -p topol.top -n index_ions.ndx -o nvt_10ns_$i.tpr
	gmx mdrun -v -deffnm nvt_10ns_$i -pf pullf$i.xvg -px pullx$i.xvg
done

MD Run

记得指定用-pf和-px分别指定输出的pullf和pullx文件命，如：

gmx mdrun -deffnm umbrella0 -pf pullf-umbrella0.xvg -px pullx-umbrella0.xvg

否则因为输出文件同名而会被覆盖.

PMF

建立两个文件（可以用ls *.tpr > tpr-files.dat等），tpr-files.dat和pullf-files.dat，分别包括tpr和pullf.xvg文件名：如tpr-files.dat:

umbrella0.tpr
umbrella1.tpr
...
umbrella22.tpr

或如pullf-files.dat:

pullf0.tpr
pullf1.tpr
...
pullf22.tpr

然后用下述指令计算PMF，注意选择从平衡后开始：

gmx wham -it tpr-files.dat -if pullf-files.dat -b xxx -o -hist -unit kCal

手机版“神探玺洛克”请扫码