因为后期实验需要用到,最近几天在学习拉伸动力学模拟,虽然拉伸动力学教程为5.X版本,但是自5.1版本对于pull模块还是有较大的变化,这里和大家分享一下其中遇到的问题和我自己的解决方法。

英文版已经升级到5.1最新了,可能用5.1版本不会遇到这些问题

首先第一条就遇到了一点小问题

gmx pdb2gmx -f input.pdb -ignh -ter -o complex.gro

其提示选择什么水模型,我选择的SPC模型,个人觉得加-water tip3p更佳?因为不知道最后采用的什么水模型。也没要TIP3P,TIP4P之类的选项让人没有底。

在进行构型生成的时候,会有一系列的错误,我们首先来看教程中的

1
2
3
4
5
6
7
8
9
10
11
12
13
; Pull code
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 ; simple distance increase
pull_coord1_groups = 1 2
pull_coord1_dim = N N Y
pull_coord1_rate = 0.01 ; 0.01 nm per ps = 10 nm per ns
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2
pull_start = yes ; define initial COM distance > 0

第一个为pull_start提示的warning,现在该参数在新版中已经没有,需要改成pull_coord1_start (英文版本已经改过来了)
参考的为:https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2016-March/104109.html

NOTE 1提示的没有设置特定的cutoff-scheme的值,cutoff-schemegromacs中只有group方案和Verlet方案,group方案较老,可能新的版本已经不存在了,同时速度也较差,具体可以查看。http://www.gromacs.org/Documentation/Cut-off_schemes

nstlist可以设置超过10,因为在Verlet算法中nstlist不影响精确度

leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1这个问题有人说无关紧要,有人说可以
http://gromacs.org_gmx-users.maillist.sys.kth.narkive.com/cXo18Y84/dr-lemkul-s-umbrella-sampling-tutorial-grompp-note-on-leap-frog-nose-hoover

nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy中nstomm选项是用来去除质心的整体运动的,如果不使用的话,体系整体会有平动速度,沿某一方向平移,去除后就没有这种现象了,轨迹看起来也正常一些。一个办法是不去管他(我的方法),另外一个可以nstcalcenergy设置为10

distances.pl脚本李老师翻译的版本中gmx distance-oaxy命令已经过时,可以使用英文网站上的perl脚本,但是我使用还是会有阅读不全的问题,好像是阻塞的缘故,但是自己不会用perl就重新用python写了一个,如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#-*-coding:utf-8-*-
# usage:python distances.py [filenum=501]
# example:python distances.py
import os,re,sys
try:
num=int(sys.argv[1])
except:
num=501
for i in range(num):
cmd="gmx distance -s pull.tpr -f conf%s.gro -n index.ndx -oall dist%s.xvg -select \'com of group \"Chain_A\" plus com of group \"Chain_B\"\' " %(str(i),str(i))
os.system(cmd)
output=open('summary_distances.dat','a')
for i in range(num):
file='dist'+str(i)+'.xvg'
if os.path.exists(file):
input=open(file,'r')
for line in input:
if re.match(r'^[#@]',line):
pass
else:
num,distance=line.split()
#print(distance)
output.write(str(i)+"\t")
output.write(distance)
output.write('\n')
input.close()
output.close()

setupUmbrella.pypython版本为2.x,像我这种用3版本的就尴尬了,又不想重新写,就用conda环境包弄一下。

1
2
conda create --name python2 python=2.7
source activate python

 

run-umbrella.sh中需要前面需要加gmx,如下:

1
2
3
4
5
6
7
#!/bin/bash
# Short equilibration
gmx grompp -f npt_umbrella.mdp -c confXXX.gro -p topol.top -n index.ndx -o nptXXX.tpr
gmx mdrun -deffnm nptXXX
# Umbrella run
gmx grompp -f md_umbrella.mdp -c nptXXX.gro -t nptXXX.cpt -p topol.top -n index.ndx -o umbrellaXXX.tpr
gmx mdrun -deffnm umbrellaXXX -pf pullf-umbrellaXXX.xvg -px pullx-umbrellaXXX.xvg

后面也是pull_coord1_start的问题,修改一下即可