Gromacs

GROMACS(全称:英语:GROningen MAchine,全称格罗宁根华讯模拟体系)是一套分子动力学模拟程序包,主要用来模拟研究蛋白质、脂质、核酸等生物分子的性质。它起初由荷兰格罗宁根大学生物化学系开发,目前由来自世界各地的大学和研究机构维护。平时用GROMACS比较多,故总结一些常用的命令、技巧,分享一些资源,方便查找。

1. gromacs 的安装

2. gromacs xtc文件的输出

grommp.mdp控制文件中加入以下语句就会生成traj_com.xtc文件

nstxtcout                = 10

模拟轨迹单精度文件(.xtc)。单精度轨迹文件,文件较.trr.trj小很多,为常用分析文件。包含模拟系统中原子坐标,模拟时间,和模拟盒子信息。

3. MDAnalysis的使用

3.1. 读取轨迹文件

不废话,上代码:

#!/usr/bin/env  python
import MDAnalysis as mda

top = mda.topology.TPRParser.TPRParser("topol.tpr").parse()
u = mda.Universe(top, "traj_comp.xtc")
print(" trajectory loaded!")

MDAnalysis的轨迹读入很有意思,它会根据文件格式自动分析,如果你有top文信息可以从topol.tpr文件中分析得到,然后一起传入mda.Universe,会得到完整的体系的信息。

4. 轨迹文件坐标展开成真实坐标

#!/bin/sh
# usage: sh convert.sh frame_42.pdb
echo "0" > tmp.dat
pdb=$1
atom_pdb=$(echo $pdb |cut -d . -f1)"_atom.pdb"
atom_whole_pdb=$(echo $atom_pdb |cut -d . -f1)"_whole.pdb"
gmx_mpi trjconv -f $pdb -pbc atom -o $atom_pdb  < tmp.dat > /dev/null 2>&1
gmx_mpi trjconv -f $atom_pdb -pbc whole -o $atom_whole_pdb  < tmp.dat > /dev/null 2>&1
rm tmp.dat
rm $atom_pdb

其中输入的0可以用EOF完成。<< EOF EOF的作用是在命令执行过程中用户自定义输入,它类似于起到一个临时文件的作用,只是比使用文件更方便灵活。
这里<< 是一个特殊的符号,是标准输入来自命令行的一对分隔号的中间内容。

#!/bin/sh
# usage: sh convert.sh frame_42.pdb
#!/bin/sh
pdb=$1
atom_pdb=$(echo $pdb |cut -d . -f1)"_atom.pdb"
atom_whole_pdb=$(echo $atom_pdb |cut -d . -f1)"_whole.pdb"
gmx_mpi trjconv -f $pdb -pbc atom -o $atom_pdb   > /dev/null 2>&1 < /dev/null 2>&1 << EOF
0
EOF
rm $atom_pdb

cat 与EOF的结合tips

cat << EOF >test.sh 
> red
> yellow
> blue
> EOF

即先向临时文件写入几行,然后cat到test.sh文件中。

Note 1: EOF是END Of File的缩写,表示自定义终止符。所以EOF可以换成别的字符。
Note 2: 发现EOF需要顶格写,呵呵哒

5. 延长作业的命令

gmx convert-tpr -s topol.tpr -f traj.trr -e ener.edr -o (一个新的文件夹)/topol.tpr -time  (时间) -until (时间)
在新的文件夹下gmx mdrun -s topol.tpr &

6. coarse grained AA gromacs trajectory

目标:将所有残基的质心提取出来
首先将所有原子的index按照residue分成N个groups,对每个groups输出com即可,用如下脚本得到mapping.ndx文件:

#!/usr/bin/env python
import sys
import MDAnalysis as mda
import numpy as np

TprName, XtcName, OUTPUT = sys.argv[1], sys.argv[2], sys.argv[3]
top = mda.topology.TPRParser.TPRParser(TprName).parse()
u = mda.Universe(top, XtcName)
with open(OUTPUT,'w') as NDX:
    for i in range(u.residues.resnames.shape[0]):
        ids = u.residues[i].atoms.ids
        ats = " ".join(map(str, ids))
        NDX.write("[ bead%d ]\n%s\n" % (i+1, ats))

7. gromacs command

7.1. gmx editconf

Synopsis
gmx editconf [-f [<.gro/.g96/…>]] [-n [<.ndx>]] [-bf [<.dat>]] [-o [<.gro/.g96/…>]] [-mead [<.pqr>]] [-[no]w] [-[no]ndef] [-bt ] [-box ] [-angles ] [-d ] [-[no]c] [-center ] [-aligncenter ] [-align ] [-translate ] [-rotate ] [-[no]princ] [-scale ] [-density ] [-[no]pbc] [-resnr ] [-[no]grasp] [-rvdw ] [-[no]sig56] [-[no]vdwread] [-[no]atom] [-[no]legend] [-label ] [-[no]conect]
Note:

  1. -scale value 可以缩放坐标,如0.9,可以将所有坐标缩放至原来的0.9倍,盒子信息也会更新为原来的0.9倍。
  2. -density value 改变体系密度,应该也是通过scale坐标来实现的。
  3. -bt value 决定盒子类型triclinic, cubic,
  4. -box value 对于cubic的盒子类型只需要一个参数,因为盒子类型已经决定了其三边相等(nm单位)
  5. -d-box 都会将体系重置到盒子中心
    gmx_mpi editconf -f conf.gro -o shrink.gro -bt cubic -scale 0.968555900621118