#分子模拟#Gromacs 5.1.2做拉伸动力学的几点笔记

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

英文版已经升级到5.1最新了，可能用5.1版本不会遇到这些问题
首先第一条就遇到了一点小问题

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

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

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

; 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写了一个，如下：

#-*-coding:utf-8-*-
#-*-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环境包弄一下。

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

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

#!/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的问题，修改一下即可

更多原创精彩视频敬请关注生信杂谈: