美文网首页生物信息杂谈分子模拟
#分子模拟#Gromacs 5.1.2做拉伸动力学的几点笔记

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

作者: 生信杂谈 | 来源:发表于2017-08-19 15:05 被阅读335次

    因为后期实验需要用到,最近几天在学习拉伸动力学模拟,虽然拉伸动力学教程为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的问题,修改一下即可

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

    相关文章

      网友评论

        本文标题:#分子模拟#Gromacs 5.1.2做拉伸动力学的几点笔记

        本文链接:https://www.haomeiwen.com/subject/looorxtx.html