• 2018-07-11 20:41:53

以前简单说过使用GROMACS进行团簇分析的方法, 是利用脚本组合mdmat和cluster这两个工具完成的. 最近又需要进行团簇分析, 就把这两部分的代码整理了一下, 组合成一个新程序mdcluster, 并将其编译为GROMACS的内置模块, 这样就可以直接gmx mdcluster调用了, 更方便, 速度也更快了.

编写GROMACS分析程序涉及到的知识在前两篇文档中有说明, 但对于不熟悉C++的人来说, 还是很难看明白的. 所以我这里就具体地说下我的作法, 供参考.

1. 在src\gromacs\gmxana\gmx_ana.h中定义要添加的程序

   int gmx_mdcluster(int argc, char *argv[]);

2. 在src\programs\legacymodules.cpp中的// Modules from gmx_ana.h.部分注册模块, 这样才能在gmx主程序中调用

   registerModule(manager, &gmx_mdcluster, "mdcluster",                "md Cluster structures");

3. 将src\gromacs\gmxana\gmx_mdmat.cpp和src\gromacs\gmxana\gmx_cluster.cpp两个文件合并, 命名为mdcluster.cpp

4. 整理mdcluster.cpp的头文件, 修改主函数的名称为mdcluster

5. 试着进行编译. 遇到函数重定义的错误, 就修改函数名称. 直到编译通过, 执行gmx mdcluster成功

6. 根据需要编写需要的功能. 主要是弄明白mdmat函数输出的距离矩阵, 并在cluster函数中调用它.

7. 根据需要, 添加其他需要的功能. 必要时可参考GROMACS自带的其他分析程序.

步骤说起来并不复杂, 但实际做起来还是有点难度的, 虽然不需要你熟悉C++, 但至少需要要像我一样, 学过基本的C, 否则很多语句会看不懂, 也就无从修改了.

另外, 我的做法比较脏快, 好在也能解决问题, 所以我也就不再考虑更优雅的作法了.

代码

测试GROMACS 2016.5和2018版本, 都正常通过.

使用方法

gmx mdcluster -f -s -n

选项

• -f: 默认traj.trr, 要分析的轨迹文件

• -s: 默认topol.tpr, 运行输入文件

• -n: 默认index.ndx, 可选索引文件

• -g: 默认cluster.log, 输出文件, 包含每个时刻的团簇信息

• -num: 默认clust-num.xvg, 可选输出xvg文件, 每个时刻的团簇数目

• -xyz: 默认clust-xyz.pdb, 可选输出xyz格式的坐标文件, 列出每个团簇的坐标. 之所以扩展名为pdb, 是因为GROMACS不支持指定xyz格式的输出文件, 只能使用pdb代替了.

• 其他涉及团簇分析的选项没有改变, 参考cluster的文档.

输出文件

log文件

按时间顺序列出团簇数目, 距离矩阵的简单信息, 以及每个团簇包含的残基的编号.

Time: 7.4    Clusters: 3
Distance: 0.168967 to 1.00864 nm    Average: 0.473927
Energy of the matrix: 2.95684#cluster | members(#res)
       1 | 1 4 5 7 9       2 | 2 6 10       3 | 3 8

Time: 7.6    Clusters: 4
Distance: 0.171863 to 1.03121 nm    Average: 0.462375
Energy of the matrix: 2.8525#cluster | members(#res)
       1 | 1 4 5 7 9       2 | 2       3 | 3 6 10       4 | 8

Time: 7.8    Clusters: 6
Distance: 0.176229 to 1.07498 nm    Average: 0.472193
Energy of the matrix: 3.18141#cluster | members(#res)
       1 | 1 7 9       2 | 2       3 | 3 8       4 | 4       5 | 5       6 | 6 10

团簇数目

很简单的xvg文件, 无需多说

# This file was created Fri Jul  6 09:24:17 2018# Created by:#           :-) GROMACS - gmx mdcluster, 2016.5 (double precision) (-:## Executable:   C:\Users\Jicun\Downloads\GMX2016.5\bin\Release\gmx_d.exe# Data prefix:  C:\Program Files\Gromacs# Working dir:  C:\Users\Jicun\Downloads\GMX2016.5\bin\Release# Command line:#   gmx_d mdcluster -g -no -cutoff .2 -xyz# gmx mdcluster is part of G R O M A C S:## Great Red Owns Many ACres of Sand#@    title "Cluster Numbers"@    xaxis  label "Time (ps)"@    yaxis  label "Cluster Number"@TYPE xy
@g0 type bar
  0.000000        9
  0.200000       10
  0.400000       10
  0.600000        9
  0.800000        8
  1.000000        8
  1.200000        9
  1.400000        9
  1.600000        9
  1.800000        9

团簇结构

列出每个时刻下每个团簇中所有原子的名称及其xyz坐标, 是常用的xyz格式, 标题行中还给出时刻, 团簇编号以及盒子大小的信息. 对团簇进行后处理的时候需要这些信息.

值得注意的是:

• 这些坐标是原始轨迹中的坐标, 没有经过PBC处理, 所以有些构型看起来好像不正确.

• 不要使用VMD来查看xyz文件, 因为VMD要求每个构型的原子数目是固定的, 但团簇分析结果显然无法满足这个条件, 所以VMD显示出来会有错误. 我一般用JMOL来查看这种文件.

9
t: 20.000000    cl: 1    box: 1.959580 0.000000 0.000000 0.000000 1.959580 0.000000 0.000000 0.000000 1.959580
      OW  16.299245  17.166265   7.652279 #res:1
     HW1  17.101772  16.680043   7.841409 #res:1
     HW2  16.426502  17.503887   6.765685 #res:1
      OW  13.016180  17.398965   3.690518 #res:3
     HW1  13.727098  17.004212   3.185545 #res:3
     HW2  13.171825  17.110682   4.589906 #res:3
      OW  13.643707  16.886055   6.493071 #res:6
     HW1  13.265058  16.279565   7.129487 #res:6
     HW2  14.567202  16.945082   6.737825 #res:69
t: 20.000000    cl: 2    box: 1.959580 0.000000 0.000000 0.000000 1.959580 0.000000 0.000000 0.000000 1.959580
      OW  12.382308   0.610490   3.052188 #res:2
     HW1  12.713196  -0.143198   3.540757 #res:2
     HW2  11.454754   0.658140   3.283711 #res:2
      OW  11.333401   1.919177   0.772347 #res:7
     HW1  11.814307   1.534172   1.504967 #res:7
     HW2  11.917793   2.596026   0.430886 #res:7
      OW  14.576921   0.086829   1.378225 #res:8
     HW1  13.888679   0.222418   2.029508 #res:8
     HW2  14.172984  -0.473075   0.715220 #res:86
t: 20.000000    cl: 3    box: 1.959580 0.000000 0.000000 0.000000 1.959580 0.000000 0.000000 0.000000 1.959580
      OW  16.912391  17.797985   2.399200 #res:4
     HW1  17.057798  17.103391   1.756840 #res:4
     HW2  16.319146  18.408840   1.961993 #res:4
      OW  16.936209  18.858645   4.993168 #res:10
     HW1  16.868745  19.803851   4.858015 #res:10
     HW2  16.985083  18.492818   4.109984 #res:103
t: 20.000000    cl: 4    box: 1.959580 0.000000 0.000000 0.000000 1.959580 0.000000 0.000000 0.000000 1.959580
      OW  13.452642  17.216688   0.500806 #res:5
     HW1  12.790335  17.045542   1.170352 #res:5
     HW2  12.953310  17.324260  -0.308717 #res:53
t: 20.000000    cl: 5    box: 1.959580 0.000000 0.000000 0.000000 1.959580 0.000000 0.000000 0.000000 1.959580
      OW  16.102296   1.618148   3.075337 #res:9
     HW1  15.229665   1.230862   3.144304 #res:9
     HW2  16.380787   1.424401   2.180275 #res:9

未完成

• 团簇分析时采用的截断距离只有一个, 这对复杂的体系可能不合适. 可扩展一下, 对不同的原子使用不同的截断值.

• 输出构型使用pdb格式可能更好, 虽然文件大些

• 可换用另外的团簇分析算法

◆本文地址: https://jerkwin.github.io/2018/07/11/GROMACS团簇分析代码/, 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆