GROMACS团簇分析代码

类别:    标签: gmx   阅读次数:   版权: (CC) BY-NC-SA

以前简单说过使用GROMACS进行团簇分析的方法, 是利用脚本组合mdmatcluster这两个工具完成的. 最近又需要进行团簇分析, 就把这两部分的代码整理了一下, 组合成一个新程序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.cppsrc\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

选项

输出文件

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格式, 标题行中还给出时刻, 团簇编号以及盒子大小的信息. 对团簇进行后处理的时候需要这些信息.

值得注意的是:

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

未完成

◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: GROMACS编写轨迹分析工具的示例代码(gmx2018.2)
后一篇: MSD算扩散系数的几种方法

访问人次(2015年7月 9日起): | 最后更新: 2024-05-27 02:08:45 UTC | 版权所有 © 2008 - 2024 Jerkwin