仓库源文，站点原文

layout: post title: 分子片段分析 categories:

在涉及固体和液体的计算中, 一个体系中通常会包含多个相同的分子. 这些分子加入时使用的方法不同, 其原子的排列次序可能很乱, 相互之间的对应次序也会不同. 而在做计算时, 我们习惯让这些分子按顺序排列, 属于同一个分子的原子依次排列, 而且不同分子间的原子一一对应, 这样在编程处理时会方便很多.

上面的论述可能有点抽象, 下面是一个例子.

这是一个由晶体文件创建的体系, 里面包含12个分子, 每个分子又包含65原子.

这些原子的排列顺序是乱的, 隶属于同一个分子的原子会分散到多处, 而不是集中在一起. 我们想要对这些原子的排列顺序进行整理, 使得首先是属于第一个分子的原子, 然后是属于第二个的, 第三个的..., 而且最好不同分子之间的原子可以一一对应好.

第一个目的相当于对整个体系进行片段分析, 根据原子间的连接性将原子划分到不同的分子中. 只要知道了原子间的连接性, 可以使用图论中的遍历算法来完成. 具体来说可采用广度优先的遍历方法(BFS). 这个算法还有其他的名字, 比如最大相连子图, 连接元件等, 是个经典算法. 详细的说明网上很多资料.

将原子划分为不同的分子后, 不同分子间的原子如何一一对应呢? 这个问题本质上属于广义的匹配问题. 二维图的匹配已经是NP难度的了, 再考虑三维结构的匹配就更加困难了. 虽然有一些算法可以解决, 但分子中的原子数多了以后就不适用了, 而且当分子间的差异过大时, 没有算法可以保证一定有效.

但当各个分子间的差异不大时, 我想到的方法有下面几种:

  1. 先根据连接性进行子图匹配, 列出所有可能的匹配方案, 然后再根据叠合后的RMSD来选择最佳匹配. 这种方法很严格, 但可能很耗时.
  2. 先对分子进行旋转, 使不同分子的取向一致, 然后根据距离最小规则进行原子匹配. 确定如何旋转时, 可使用某些特征量做为坐标轴, 如惯性主轴, 电四极矩主轴等, 也可手动选择不共线的三个特征原子, 以它们确定坐标轴. 前一方法对高对称性的分子可能会有问题, 而后一方法需要手动操作, 各有优缺点, 应根据实际情况选用.
  3. GaussView软件的Edit | Connection提供了一个基于连接信息的匹配功能. 如果分子没有对称性的话, 基本可以正确匹配, 但当分子具有对称性时, 会存在多种匹配可能, 给出的匹配可能有误.

总言之, 上面的这些方法, 无论哪种都不能保证结果一定合理, 但至少可以提供一个比较合理的初步匹配, 供你最后手动调整一下.

Talk is cheap, show me the CODE

下面的代码用于对*.gjf格式的文件进行片段分析, 这种文件在最后可以保留原子间的链接信息. 当然, 稍加修改就可以用于其他格式的文件, 如pdb, mol2等. 代码最后会将所有的分子片段输出到一个*~Frg.xyz文件, 用于检查, 同时还会将每个分子片段分别输出到一个*~#.gjf文件, 用来提供给GaussView以便使用其Edit | Connection功能.

<table class="highlighttable"><th colspan="2" style="text-align:left">FrgMol.bsh</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span style="color: #B8860B">usage</span><span style="color: #666666">=</span><span style="color: #BB4444">"\</span> <span style="color: #BB4444">>>>>>>>>>>>>>>>> FrgMol <<<<<<<<<<<<<<<<</span> <span style="color: #BB4444">>>>>>>>>>>>>>>>> Jicun LI <<<<<<<<<<<<<<<<</span> <span style="color: #BB4444">>>>>>>>>>> 2016-09-22 16:19:51 <<<<<<<<<</span> <span style="color: #BB4444">>> Usage: FrgMol <File.gjf>"</span> <span style="color: #666666">[[</span> <span style="color: #B8860B">$#</span> -lt <span style="color: #666666">1</span> <span style="color: #666666">]]</span> <span style="color: #666666">&&</span> <span style="color: #666666">{</span> <span style="color: #AA22FF">echo</span> <span style="color: #BB4444">"</span>$<span style="color: #BB4444">usage"</span>; exit; <span style="color: #666666">}</span> <span style="color: #B8860B">File</span><span style="color: #666666">=</span><span style="color: #BB6688; font-weight: bold">${</span><span style="color: #B8860B">1</span>%.gjf<span style="color: #BB6688; font-weight: bold">}</span> awk -v <span style="color: #B8860B">File</span><span style="color: #666666">=</span><span style="color: #B8860B">$F</span>ile <span style="color: #BB4444">'</span> BEGIN<span style="color: #666666">{</span> <span style="color: #B8860B">Nblk</span><span style="color: #666666">=</span>0; <span style="color: #B8860B">Ntot</span><span style="color: #666666">=0</span> <span style="color: #666666">}</span> <span style="color: #B8860B">NF</span><span style="color: #666666">==0</span> <span style="color: #666666">&&</span> Nblk<<span style="color: #666666">=2</span> <span style="color: #666666">{</span> Nblk++; next <span style="color: #666666">}</span> <span style="color: #B8860B">Nblk</span><span style="color: #666666">==2</span> <span style="color: #666666">{</span> <span style="color: #AA22FF; font-weight: bold">while</span><span style="color: #666666">(</span>NF>0<span style="color: #666666">)</span> <span style="color: #666666">{</span> getline; Ntot++ S<span style="color: #666666">[</span>Ntot<span style="color: #666666">]=</span><span style="color: #B8860B">$1</span>; X<span style="color: #666666">[</span>Ntot<span style="color: #666666">]=</span><span style="color: #B8860B">$2</span>; Y<span style="color: #666666">[</span>Ntot<span style="color: #666666">]=</span><span style="color: #B8860B">$3</span>; Z<span style="color: #666666">[</span>Ntot<span style="color: #666666">]=</span><span style="color: #B8860B">$4</span> <span style="color: #666666">}</span> Ntot-- <span style="color: #AA22FF; font-weight: bold">for</span><span style="color: #666666">(</span><span style="color: #B8860B">i</span><span style="color: #666666">=</span>1; i<<span style="color: #666666">=</span>Ntot; i++<span style="color: #666666">)</span> <span style="color: #666666">{</span> getline Imol<span style="color: #666666">[</span>i<span style="color: #666666">]=</span>0; Adj<span style="color: #666666">[</span>i,i<span style="color: #666666">]=</span>1 <span style="color: #AA22FF; font-weight: bold">for</span><span style="color: #666666">(</span><span style="color: #B8860B">k</span><span style="color: #666666">=</span>1; k<<span style="color: #666666">=(</span>NF-1<span style="color: #666666">)</span>/2; k++<span style="color: #666666">)</span> <span style="color: #666666">{</span> Adj<span style="color: #666666">[</span>i,<span style="color: #AA22FF; font-weight: bold">$(</span>2*k<span style="color: #AA22FF; font-weight: bold">)</span><span style="color: #666666">]=</span>1; Adj<span style="color: #666666">[</span><span style="color: #AA22FF; font-weight: bold">$(</span>2*k<span style="color: #AA22FF; font-weight: bold">)</span>,i<span style="color: #666666">]=1</span> <span style="color: #666666">}</span> <span style="color: #666666">}</span> <span style="color: #B8860B">Nmol</span><span style="color: #666666">=</span>0 <span style="color: #AA22FF; font-weight: bold">for</span><span style="color: #666666">(</span><span style="color: #B8860B">i</span><span style="color: #666666">=</span>1; i<<span style="color: #666666">=</span>Ntot; i++<span style="color: #666666">)</span> <span style="color: #666666">{</span> <span style="color: #AA22FF; font-weight: bold">if</span><span style="color: #666666">(</span>Imol<span style="color: #666666">[</span>i<span style="color: #666666">]==</span>0<span style="color: #666666">)</span> <span style="color: #666666">{</span> Nmol++; que<span style="color: #666666">[</span>i<span style="color: #666666">]</span>++ <span style="color: #AA22FF; font-weight: bold">while</span><span style="color: #666666">(</span>length<span style="color: #666666">(</span>que<span style="color: #666666">)</span>>0<span style="color: #666666">)</span> <span style="color: #666666">{</span> <span style="color: #AA22FF; font-weight: bold">for</span><span style="color: #666666">(</span>v in que<span style="color: #666666">)</span> <span style="color: #666666">{</span> Imol<span style="color: #666666">[</span>v<span style="color: #666666">]=</span>Nmol; Natm<span style="color: #666666">[</span>Nmol<span style="color: #666666">]</span>++ <span style="color: #AA22FF; font-weight: bold">for</span><span style="color: #666666">(</span><span style="color: #B8860B">j</span><span style="color: #666666">=</span>1; j<<span style="color: #666666">=</span>Ntot; j++<span style="color: #666666">)</span> <span style="color: #666666">{</span> <span style="color: #AA22FF; font-weight: bold">if</span><span style="color: #666666">(</span>Adj<span style="color: #666666">[</span>v,j<span style="color: #666666">]</span> <span style="color: #666666">&&</span> Imol<span style="color: #666666">[</span>j<span style="color: #666666">]==</span>0<span style="color: #666666">)</span> que<span style="color: #666666">[</span>j<span style="color: #666666">]</span>++ <span style="color: #666666">}</span> delete que<span style="color: #666666">[</span>v<span style="color: #666666">]</span> <span style="color: #666666">}</span> <span style="color: #666666">}</span> <span style="color: #666666">}</span> <span style="color: #666666">}</span> <span style="color: #AA22FF; font-weight: bold">for</span><span style="color: #666666">(</span><span style="color: #B8860B">j</span><span style="color: #666666">=</span>1; j<<span style="color: #666666">=</span>Nmol; j++<span style="color: #666666">)</span> <span style="color: #666666">{</span> print Natm<span style="color: #666666">[</span>j<span style="color: #666666">]</span> print <span style="color: #BB4444">"Mol "</span> j <span style="color: #BB4444">"/"</span>Nmol <span style="color: #B8860B">Fout</span><span style="color: #666666">=</span>File<span style="color: #BB4444">"~"</span>j<span style="color: #BB4444">".gjf"</span> print <span style="color: #BB4444">"# \n\nFrg "</span>j<span style="color: #BB4444">"/"</span>Nmol<span style="color: #BB4444">"\n\n0 1"</span> >Fout <span style="color: #AA22FF; font-weight: bold">for</span><span style="color: #666666">(</span><span style="color: #B8860B">i</span><span style="color: #666666">=</span>1; i<<span style="color: #666666">=</span>Ntot; i++<span style="color: #666666">)</span> <span style="color: #666666">{</span> <span style="color: #AA22FF; font-weight: bold">if</span><span style="color: #666666">(</span>Imol<span style="color: #666666">[</span>i<span style="color: #666666">]==</span>j<span style="color: #666666">)</span> <span style="color: #666666">{</span> <span style="color: #AA22FF">printf</span> <span style="color: #BB4444">"%4s %14.8f %14.8f %14.8f\n"</span>, S<span style="color: #666666">[</span>i<span style="color: #666666">]</span>, X<span style="color: #666666">[</span>i<span style="color: #666666">]</span>, Y<span style="color: #666666">[</span>i<span style="color: #666666">]</span>, Z<span style="color: #666666">[</span>i<span style="color: #666666">]</span> <span style="color: #AA22FF">printf</span> <span style="color: #BB4444">"%4s %14.8f %14.8f %14.8f\n"</span>, S<span style="color: #666666">[</span>i<span style="color: #666666">]</span>, X<span style="color: #666666">[</span>i<span style="color: #666666">]</span>, Y<span style="color: #666666">[</span>i<span style="color: #666666">]</span>, Z<span style="color: #666666">[</span>i<span style="color: #666666">]</span> > Fout <span style="color: #666666">}</span> <span style="color: #666666">}</span> <span style="color: #666666">}</span> <span style="color: #666666">}</span> <span style="color: #BB4444">'</span> <span style="color: #B8860B">$F</span>ile.gjf > <span style="color: #B8860B">$F</span>ile~Frg.xyz </pre></div> </td></tr></table>

如果需要测试的话, 点击这里下载脚本和测试文件吧.

GaussView连接匹配功能的使用

官网的说明见The Connection Editor, 下面以上面给出的测试文件示例一下.

运行bash FrgMol.bsh Mol.gjf会得到Mol~Frg.xyz和Mol~1.gjf, Mol~2.gjf, ...Mol~12.gjf.

打开GaussView, 点击File | Open..., 按住Shift键, 选中12个gjf文件, 注意最下面的选择.

打开所有这些分子后, 点击Y形图标将所有12个分子排列好

菜单Edit | Connection打开连接编辑器, 在第一个分子中随便选择一个原子, 然后点击Enable Autofixing, 后面所有分子就会自动调整顺序.

要保存的话, 点击OK, 退出连接编辑器, 然后分别保存每个分子, 保存后的结构中原子的顺序就是对应好的了. 如果分子很多的话, 手动做起来也还是麻烦.

还要注意的是, 在连接编辑器中也可以手动指定每个原子的对应原子, 作法是左键选中第一个分子中的一个原子, 在其他分子中使用右键指定对应原子. 这样做比手动调整坐标顺序方便一些.