- 2018-07-11 20:41:53
以前简单说过使用GROMACS进行团簇分析的方法, 是利用脚本组合mdmat和cluster这两个工具完成的. 最近又需要进行团簇分析, 就把这两部分的代码整理了一下, 组合成一个新程序mdcluster, 并将其编译为GROMACS的内置模块, 这样就可以直接gmx mdcluster调用了, 更方便, 速度也更快了.
编写GROMACS分析程序涉及到的知识在前两篇文档中有说明, 但对于不熟悉C++的人来说, 还是很难看明白的. 所以我这里就具体地说下我的作法, 供参考.
- 在
src\gromacs\gmxana\gmx_ana.h中定义要添加的程序int gmx_mdcluster(int argc, char *argv[]);
- 在
src\programs\legacymodules.cpp中的// Modules from gmx_ana.h.部分注册模块, 这样才能在gmx主程序中调用registerModule(manager, &gmx_mdcluster, "mdcluster", "md Cluster structures");
-
将
src\gromacs\gmxana\gmx_mdmat.cpp和src\gromacs\gmxana\gmx_cluster.cpp两个文件合并, 命名为mdcluster.cpp -
整理
mdcluster.cpp的头文件, 修改主函数的名称为mdcluster -
试着进行编译. 遇到函数重定义的错误, 就修改函数名称. 直到编译通过, 执行
gmx mdcluster成功 -
根据需要编写需要的功能. 主要是弄明白
mdmat函数输出的距离矩阵, 并在cluster函数中调用它. - 根据需要, 添加其他需要的功能. 必要时可参考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:6
9
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:8
6
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:10
3
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:5
3
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格式可能更好, 虽然文件大些
- 可换用另外的团簇分析算法
