给cp2k输出的pdb轨迹加残基名，以及把pdb轨迹转化为trr的脚本

ggdh

给cp2k输出的pdb轨迹加残基名的脚本：
更新：目前可以把atomname也一并转移过去了
cp2k输出的pdb轨迹没有残基名，残基号，chain等信息，分析起来有点不方便。于是根据ai的提示写了个脚本（ai胡编了一个不可行的方案，但是它的思路可以用），用法如下：
1. 用户需要准备一个带有残基名的pdb文件,REF.pdb，其中残基名，残基号，以及chain信息(pdb文件的第18-26列)会被使用。
    一般用packmol产生的pdb文件就已经带有这些信息了

2. 跑完cp2k后得到 sys-pos-1.pdb轨迹文件，这个文件中应该没有任何残基信息了。然后运行

1.    sh setres.sh REF.pdb sys-pos-1.pdb > sys-with-resname.pdb

复制代码

注意：
a. REF.pdb中的原子顺序必须和轨迹中的原子顺序一模一样。
b. xyz格式的轨迹文件不支持加入残基信息。

pdb轨迹转化为trr的脚本：
其实有残基名的轨迹pdb文件也不好直接分析所以setres.sh脚本的实际意义大概就是方便把pdb轨迹文件直接拖入vmd中看带有resname的轨迹
所以如果要做分析（使用gromacs的工具或者python包比如MDAnalysis）更好方法是把pdb轨迹转化成trr文件，这样可以用gromacs读取带残基的pdb文件作为top，然后读取trr文件后，直接对轨迹进行分析
论坛上的一个可行的方案(http://bbs.keinsci.com/thread-35685-1-1.html)是用vmd读取有残基名的pdb，然后载入后续轨迹信息，然后再保存轨迹
这个方案要涉及到鼠标操作很费时间，于是我把该方案改编成一个命令，也就是这里的pdb转化为trr的脚本。(需要安装了vmd）
用法

1. sh pdb2trr sys-pos-1.pdb

复制代码

这会产生一个sys-pos-1.trr
然后再gmx中使用比如

1. gmx_mpi rdf -f sys-pos-1.trr -s REF.pdb -n index.idx

复制代码

或者MDAnalysis中使用

1. u=mda.Universe("REF.pdb","sys-pos-1.trr")

复制代码

开始进行分析


使得MDAnalysis直接读取cp2k pdb轨迹的命令：
使用下面的命令把轨迹改一下，就能使得cp2k的pdb轨迹cp2k_traj.pdb 变成能够被MDAnalysis读取的mda_traj.pdb

1. awk '{if ($0 ~ /^CRYST1/ || $0 ~/^END/ || $0 ~/^ATOM/) print $0; else if ($0 ~ /REMARK/) print "MODEL " sub(",","",$3)}' cp2k_traj.pdb  > mda_traj.pdb

复制代码

主要是把REMARK line变成MODEL line 以及删除一些没用的行