仓库源文，站点原文

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

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

这是一个由晶体文件创建的体系, 里面包含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, 退出连接编辑器, 然后分别保存每个分子, 保存后的结构中原子的顺序就是对应好的了. 如果分子很多的话, 手动做起来也还是麻烦.

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