仓库源文，站点原文

layout: post title: GROMACS团簇分析代码 categories:

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

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

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

    <div class="highlight"><pre style="line-height:125%"><span></span><span style="color: #00BB00; font-weight: bold">int</span> <span style="color: #00A000">gmx_mdcluster</span>(<span style="color: #00BB00; font-weight: bold">int</span> argc, <span style="color: #00BB00; font-weight: bold">char</span> <span style="color: #666666">*</span>argv[]);</pre></div>

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

    <div class="highlight"><pre style="line-height:125%"><span></span>registerModule(manager, <span style="color: #666666">&</span>gmx_mdcluster, <span style="color: #BB4444">"mdcluster"</span>, <span style="color: #BB4444">"md Cluster structures"</span>);</pre></div>

  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版本, 都正常通过.

使用方法

<div class="highlight"><pre style="line-height:125%"><span style="color:#A2F">gmx</span> mdcluster <span style="color:#666">-f</span> <span style="color:#666">-s</span> <span style="color:#666">-n</span></pre></div>

选项

输出文件

log文件

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

<div class="highlight"><pre style="line-height:125%"><span></span>Time: 7.4 Clusters: 3 Distance: 0.168967 to 1.00864 nm Average: 0.473927 Energy of the matrix: 2.95684 <span style="color: #008800; font-style: italic">#cluster | members(#res)</span> <span style="color: #666666">1</span> | <span style="color: #666666">1</span> <span style="color: #666666">4</span> <span style="color: #666666">5</span> <span style="color: #666666">7</span> 9 <span style="color: #666666">2</span> | <span style="color: #666666">2</span> <span style="color: #666666">6</span> 10 <span style="color: #666666">3</span> | <span style="color: #666666">3</span> 8 Time: 7.6 Clusters: 4 Distance: 0.171863 to 1.03121 nm Average: 0.462375 Energy of the matrix: 2.8525 <span style="color: #008800; font-style: italic">#cluster | members(#res)</span> <span style="color: #666666">1</span> | <span style="color: #666666">1</span> <span style="color: #666666">4</span> <span style="color: #666666">5</span> <span style="color: #666666">7</span> 9 <span style="color: #666666">2</span> | 2 <span style="color: #666666">3</span> | <span style="color: #666666">3</span> <span style="color: #666666">6</span> 10 <span style="color: #666666">4</span> | 8 Time: 7.8 Clusters: 6 Distance: 0.176229 to 1.07498 nm Average: 0.472193 Energy of the matrix: 3.18141 <span style="color: #008800; font-style: italic">#cluster | members(#res)</span> <span style="color: #666666">1</span> | <span style="color: #666666">1</span> <span style="color: #666666">7</span> 9 <span style="color: #666666">2</span> | 2 <span style="color: #666666">3</span> | <span style="color: #666666">3</span> 8 <span style="color: #666666">4</span> | 4 <span style="color: #666666">5</span> | 5 <span style="color: #666666">6</span> | <span style="color: #666666">6</span> 10</pre></div>

团簇数目

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

<div class="highlight"><pre style="line-height:125%"><span></span><span style="color: #008800; font-style: italic"># This file was created Fri Jul 6 09:24:17 2018</span> <span style="color: #008800; font-style: italic"># Created by:</span> <span style="color: #008800; font-style: italic"># :-) GROMACS - gmx mdcluster, 2016.5 (double precision) (-:</span> <span style="color: #008800; font-style: italic">#</span> <span style="color: #008800; font-style: italic"># Executable: C:\Users\Jicun\Downloads\GMX2016.5\bin\Release\gmx_d.exe</span> <span style="color: #008800; font-style: italic"># Data prefix: C:\Program Files\Gromacs</span> <span style="color: #008800; font-style: italic"># Working dir: C:\Users\Jicun\Downloads\GMX2016.5\bin\Release</span> <span style="color: #008800; font-style: italic"># Command line:</span> <span style="color: #008800; font-style: italic"># gmx_d mdcluster -g -no -cutoff .2 -xyz</span> <span style="color: #008800; font-style: italic"># gmx mdcluster is part of G R O M A C S:</span> <span style="color: #008800; font-style: italic">#</span> <span style="color: #008800; font-style: italic"># Great Red Owns Many ACres of Sand</span> <span style="color: #008800; font-style: italic">#</span> @ title <span style="color: #BB4444">"Cluster Numbers"</span> @ xaxis label <span style="color: #BB4444">"Time (ps)"</span> @ yaxis label <span style="color: #BB4444">"Cluster Number"</span> @TYPE xy @g0 <span style="color: #AA22FF">type</span> 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</pre></div>

团簇结构

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

值得注意的是:

<div class="highlight"><pre style="line-height:125%"><span></span>9 t: 20.000000 cl: <span style="color: #666666">1</span> 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 <span style="color: #008800; font-style: italic">#res:1</span> HW1 17.101772 16.680043 7.841409 <span style="color: #008800; font-style: italic">#res:1</span> HW2 16.426502 17.503887 6.765685 <span style="color: #008800; font-style: italic">#res:1</span> OW 13.016180 17.398965 3.690518 <span style="color: #008800; font-style: italic">#res:3</span> HW1 13.727098 17.004212 3.185545 <span style="color: #008800; font-style: italic">#res:3</span> HW2 13.171825 17.110682 4.589906 <span style="color: #008800; font-style: italic">#res:3</span> OW 13.643707 16.886055 6.493071 <span style="color: #008800; font-style: italic">#res:6</span> HW1 13.265058 16.279565 7.129487 <span style="color: #008800; font-style: italic">#res:6</span> HW2 14.567202 16.945082 6.737825 <span style="color: #008800; font-style: italic">#res:6</span> 9 t: 20.000000 cl: <span style="color: #666666">2</span> 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 <span style="color: #008800; font-style: italic">#res:2</span> HW1 12.713196 -0.143198 3.540757 <span style="color: #008800; font-style: italic">#res:2</span> HW2 11.454754 0.658140 3.283711 <span style="color: #008800; font-style: italic">#res:2</span> OW 11.333401 1.919177 0.772347 <span style="color: #008800; font-style: italic">#res:7</span> HW1 11.814307 1.534172 1.504967 <span style="color: #008800; font-style: italic">#res:7</span> HW2 11.917793 2.596026 0.430886 <span style="color: #008800; font-style: italic">#res:7</span> OW 14.576921 0.086829 1.378225 <span style="color: #008800; font-style: italic">#res:8</span> HW1 13.888679 0.222418 2.029508 <span style="color: #008800; font-style: italic">#res:8</span> HW2 14.172984 -0.473075 0.715220 <span style="color: #008800; font-style: italic">#res:8</span> 6 t: 20.000000 cl: <span style="color: #666666">3</span> 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 <span style="color: #008800; font-style: italic">#res:4</span> HW1 17.057798 17.103391 1.756840 <span style="color: #008800; font-style: italic">#res:4</span> HW2 16.319146 18.408840 1.961993 <span style="color: #008800; font-style: italic">#res:4</span> OW 16.936209 18.858645 4.993168 <span style="color: #008800; font-style: italic">#res:10</span> HW1 16.868745 19.803851 4.858015 <span style="color: #008800; font-style: italic">#res:10</span> HW2 16.985083 18.492818 4.109984 <span style="color: #008800; font-style: italic">#res:10</span> 3 t: 20.000000 cl: <span style="color: #666666">4</span> 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 <span style="color: #008800; font-style: italic">#res:5</span> HW1 12.790335 17.045542 1.170352 <span style="color: #008800; font-style: italic">#res:5</span> HW2 12.953310 17.324260 -0.308717 <span style="color: #008800; font-style: italic">#res:5</span> 3 t: 20.000000 cl: <span style="color: #666666">5</span> 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 <span style="color: #008800; font-style: italic">#res:9</span> HW1 15.229665 1.230862 3.144304 <span style="color: #008800; font-style: italic">#res:9</span> HW2 16.380787 1.424401 2.180275 <span style="color: #008800; font-style: italic">#res:9</span></pre></div>

未完成