博文

分子片段分析

已有 3050 次阅读 2016-9-24 09:31 |系统分类:科研笔记

2016-09-23 20:15:59

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

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

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

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

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

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

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

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

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

Talk is cheap, show me the CODE

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

FrgMol.bsh
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	usage=">>>>>>>>>>>>>>>> FrgMol <<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>> Jicun LI <<<<<<<<<<<<<<<<>>>>>>>>>> 2016-09-22 16:19:51 <<<<<<<<<>> Usage: FrgMol <File.gjf>"[[$# -lt 1]]&&{echo"$usage"; exit; }File=${1%.gjf} awk -v File=$File ' BEGIN{Nblk=0; Ntot=0}NF==0&& Nblk<=2{ Nblk++; next }Nblk==2{while(NF>0){ getline; Ntot++ S[Ntot]=$1; X[Ntot]=$2; Y[Ntot]=$3; Z[Ntot]=$4} Ntot-- for(i=1; i<=Ntot; i++){ getline Imol[i]=0; Adj[i,i]=1 for(k=1; k<=(NF-1)/2; k++){ Adj[i,$(2k)]=1; Adj[$(2k),i]=1}}Nmol=0 for(i=1; i<=Ntot; i++){if(Imol[i]==0){ Nmol++; que[i]++ while(length(que)>0){for(v in que){ Imol[v]=Nmol; Natm[Nmol]++ for(j=1; j<=Ntot; j++){if(Adj[v,j]&& Imol[j]==0) que[j]++ } delete que[v]}}}}for(j=1; j<=Nmol; j++){ print Natm[j] print "Mol " j "/"Nmol Fout=File"~"j".gjf" print "# nnFrg "j"/"Nmol"nn0 1" >Fout for(i=1; i<=Ntot; i++){if(Imol[i]==j){printf"%4s %14.8f %14.8f %14.8fn", S[i], X[i], Y[i], Z[i]printf"%4s %14.8f %14.8f %14.8fn", S[i], X[i], Y[i], Z[i] > Fout }}}}'$File.gjf > $File~Frg.xyz

FrgMol.bsh

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

usage=">>>>>>>>>>>>>>>>    FrgMol      <<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>    Jicun LI    <<<<<<<<<<<<<<<<>>>>>>>>>>     2016-09-22 16:19:51     <<<<<<<<<>>   Usage: FrgMol <File.gjf>"[[$# -lt 1]]&&{echo"$usage"; exit; }File=${1%.gjf}

awk -v File=$File ' BEGIN{Nblk=0; Ntot=0}NF==0&& Nblk<=2{ Nblk++; next }Nblk==2{while(NF>0){
			getline; Ntot++
			S[Ntot]=$1; X[Ntot]=$2; Y[Ntot]=$3; Z[Ntot]=$4}
		Ntot--

		for(i=1; i<=Ntot; i++){
			getline
			Imol[i]=0; Adj[i,i]=1
			for(k=1; k<=(NF-1)/2; k++){ Adj[i,$(2*k)]=1; Adj[$(2*k),i]=1}}Nmol=0
		for(i=1; i<=Ntot; i++){if(Imol[i]==0){
				Nmol++; que[i]++
				while(length(que)>0){for(v in que){
						Imol[v]=Nmol; Natm[Nmol]++
						for(j=1; j<=Ntot; j++){if(Adj[v,j]&& Imol[j]==0) que[j]++
						}
						delete que[v]}}}}for(j=1; j<=Nmol; j++){
			print Natm[j]
			print "Mol " j "/"Nmol
			Fout=File"~"j".gjf"
			print "# nnFrg "j"/"Nmol"nn0 1" >Fout
			for(i=1; i<=Ntot; i++){if(Imol[i]==j){printf"%4s %14.8f %14.8f %14.8fn", S[i], X[i], Y[i], Z[i]printf"%4s %14.8f %14.8f %14.8fn", S[i], X[i], Y[i], Z[i] > Fout
				}}}}'$File.gjf > $File~Frg.xyz