一种在无机晶体结构中检索特定子结构的计算机方法

霍卫峰 李 乙,* 卢君然 于吉红 徐如人 李 晶

(1吉林大学无机合成与制备化学国家重点实验室,长春130012;2北京科技大学应用力学系,北京100083)

一种在无机晶体结构中检索特定子结构的计算机方法

霍卫峰1李 乙1,*卢君然1于吉红1徐如人1李 晶2

(1吉林大学无机合成与制备化学国家重点实验室,长春130012;2北京科技大学应用力学系,北京100083)

提出了一种针对无机晶体化合物的子结构检索方法.该方法以VF2子图同构算法为基础,针对无机晶体化合物的结构特点,采用了两种策略以提高子结构检索的效率:(1)引入晶体的对称性信息避免了在等价原子间进行的大量重复性计算;(2)采用结构编码预筛选可以有效地减少目标结构的数量.我们以在无机微孔分子筛数据库中进行子结构检索为例测试该方法的有效性.测试结果表明,该方法可以快速且准确地在分子筛数据库中检索包含特定子结构的记录.两种检索策略的引入大大降低了子结构检索的复杂度,检索速度可提高3-5个数量级.该方法通过Perl语言实现,具有较好的可移植性.

无机晶体材料;数据库;子结构检索;VF2算法;分子筛;构筑基元

1 引言

分子结构检索是利用计算机程序从数据库中提取分子结构信息的过程,是化学信息学的核心研究领域之一,1其中的子结构检索问题是研究最多的课题之一.简单而言,子结构检索就是利用计算机方法在分子结构数据库中寻找包含特定子结构的分子.目前,子结构检索主要应用于在有机分子结构中寻找特定的官能团,主要有ParMol、2MoFa、3FFSM、4gSpan、5Gaston、6SMSD7等方法.和有机分子一样,无机晶体材料的子结构(构筑基元)对无机晶体材料的结构、性质以及合成研究至关重要,8,9因此针对无机晶体材料的子结构检索同样具有非常重要的意义.

然而,和以独立分子为单位的有机晶体材料不同,无机晶体结构往往是三维方向无限延展的,没有像有机分子一样的特定独立单元,数以万计不同配位状态的原子相互连接形成一个巨大的拓扑网络结构.在如此巨大的体系中进行子结构检索是一个非常复杂的问题;而当检索范围扩大到由大量复杂无机结构组成的数据库时,子结构检索将变得难以实现.目前在化学领域中使用最广泛的晶体结构数据库是剑桥结构数据库(CSD),10数据库中包含剑桥晶体数据中心开发的子结构检索工具Conquest,11其原理和有机分子的检索类似,即把庞大的晶体结构分割为有限大小的晶体学不对称单元,进一步把这些不对称单元当作独立分子进行子结构检索.这种检索方法的局限性是显而易见的:首先,无机晶体的不对称单元在理论上是可以任意选取的,仅凭不对称单元无法代表整个无机骨架结构;更为重要的是,当待检索的子结构位于一个以上不对称单元的交界时,像Conquest这种只检索不对称单元内部的方法是无法识别不对称单元之间的关联的.因此以晶体学不对称单元或独立分子为基础的子结构检索方法并不适用于不具备独立单元的无机晶体结构.目前,对无机晶体结构进行子结构检索的唯一途径是将其骨架拓扑网络扩展到相当大的范围内(例如8个晶胞),而扩大的拓扑网络同时会使计算量以指数形式递增,当检索的对象是包含大量无机结构的数据库时,子结构检索将变得难以实现.

由此可见,无机晶体材料子结构检索的重要性和理论方法的局限性迫切要求新方法的诞生.基于这些问题,结合无机晶体材料的结构特点,我们提出了一种可以在无机晶体结构中有效地进行子结构检索的方法.

2 方法

将无机晶体结构中的中心原子看作节点,将节点之间的连接看作边,无机晶体结构就可以转化成图论中的图,子结构检索问题就转化成为图论中的子图同构问题.子图同构是非确定多项式时间完全问题(NPC),即没有方法可以在合理的时间内(多项式时间)解决它.12算法的时间消耗是随节点(中心原子)数的增加呈指数级增长的,可以想象,对于包含大量原子的无机晶体结构而言,传统的子图同构算法是难以实际应用的.

很多研究尝试将子图同构算法的复杂程度降至可以实际应用的程度,主要的算法包括Ullmann算法、13SD算法、14Nauty算法、15VF算法16和VF2算法17等.其中,VF2算法是在VF算法的基础上发展起来的,和VF算法相比,VF2算法成功地将空间复杂度从O(N2)降到了O(N),其中N为图中的节点数.特别是在规则图(每个节点都具有相同数目邻接节点的图)和大图(节点数大于600的图)的同构计算中,VF2算法具有比其他算法更高的性能.1,18-21因为无机晶体的拓扑网络结构绝大多数属于规则图和大图,所以我们采用VF2子图同构算法作为子结构检索方法的基础.

VF2算法是一种深度优先的回溯算法,算法的高效性得益于采用了状态空间表示(SSR)的方法和减少搜索空间的一些匹配规则.16

图1为检索图(Q)和目标图(T)示意图.如图1所示,假设要查询的图Q中有N1个节点,目标图T中有N2个节点,N1≤N2.M是Q和T之间的映射,包含点对(q,t),其中,q∈Q,t∈T,M⊂N1×N2.状态空间S是M的一个子集,用于存储图Q和图T的局部匹配.点对(qi,ti)是否可以加入S,取决于其是否满足匹配规则.在这里,我们规定点对(qi,ti)加入S的规则是:(1)qi与ti的类型要一致;(2)对于S中任一已有的点对(qi-1,ti-1),如果qi与qi-1相连,则ti和ti-1必须相连;(3)ti的邻接节点数要小于或等于qi的邻接节点数.

图1 检索图(Q)和目标图(T)示意图Fig.1 Diagrams of query graph(Q)and target graph(T)

状态空间S的初值S0为空.从查询图Q中选择任意一个未访问过的节点q1,在目标图T中寻找当前这一轮未被访问过的且满足匹配规则的对应节点.如果Q中的第一个节点q1没有在T中找到对应的节点,则Q不是T的子图;如果能找到这样的节点t1,就将点对(q1,t1)放入状态空间S中,然后从节点q1在Q中的后继节点集(与节点q相连的节点集合)中选择一个节点q2,再从节点t1在T中的后继节点集(与节点t1相连的节点集合)中寻找满足匹配规则的对应节点.如果成功找到对应节点t2,则将点对(q2, t2)放入到状态空间S中;如果找不到就回溯到上一轮,开始新的查找.当状态空间S中有N1个点对时,表明图Q是图T的子图,输出结果.

因为无机晶体结构中包含的节点数是非常庞大的,直接采用VF2算法对无机晶体结构进行检索仍然是不现实的.在此基础上,根据无机晶体材料的结构特点,我们在VF2算法的基础上做了两个重要的改进以提高检索效率.

(1)利用晶体结构的对称性.在进行匹配前,分别将待检子结构片段(查询图Q)和目标晶体结构片段(目标图T)中的晶体学独立或不对称的原子(节点)移至相应图的顶部.在进行子图匹配时,检验Q中是否有与T中晶体学独立节点相匹配的独立节点,如果没有找到满足匹配规则的节点,就可以认为待检结构不是目标结构的子结构,而并不需要遍历待查询的子结构及目标结构中的所有节点,从而大大提高计算的效率.

(2)采用编码预筛选方法.由于子图匹配算法的计算量非常大,因此如果在进行子图匹配之前能够预筛选掉大量确定无效的数据而仅保留少量有可能符合检索条件的数据,将极大提高子结构检索的效率.为此,我们设计了一个具有特定位数的二进制编码,编码中的每一位对应一种特定的构筑基元.以无机微孔分子筛晶体结构为例来说,我们采用了一个19位的二进制编码,编码中的每一位代表了19种在分子筛结构中最常见的构筑基元22(表1),编码中每一位的值1或0代表该结构中含有或不含有指定的构筑基元.在进行子图匹配计算之前,我们先通过计算机程序把待查询的子结构和目标结构分别转换成编码(如图2所示),然后对比编码中相应位上的值,如果待查询子结构编码中的某一位显示含有某个特定的构筑基元而目标结构编码显示不含有,说明待查询子结构中包含目标结构中所不包含的构筑基元,因此可以直接判定目标结构不包含待查询的子结构;如果待查询的子结构中所有的构筑基元在目标结构编码中都显示含有,则说明目标结构有可能包含待查询的子结构,需要进一步通过子图匹配计算以确定两者是否具有包含关系.对于仅包含少量目标结构的数据库而言,编码预筛选的优越性并不明显.随着数据库中目标结构数量的增加,编码预筛选的优越性将迅速显现出来.需要强调的是,目标结构的编码转换是在子结构检索之前完成并保存的,并不占用子图匹配计算的时间.

表1 19个分子筛结构编码中各位所代表的构筑基元Table 1 Prescreening code for 19 zeolite structure building units

图2 结构编码Fig.2 Structure encoding

3 结果与讨论

我们以在分子筛晶体结构中检索特定的子结构为例测试这种方法的效率.数据源采用了国际分子筛协会(IZA)官方网站提供的197种分子筛结构的cif文件.22从每个结构的cif文件出发,计算了8个晶胞内容所对应的邻接矩阵并合并在一起形成结构数据库.选取了一些已知的分子筛构筑基元作为待查询的子结构片段,在由197种分子筛晶体结构所组成的数据库中进行检索,并将检索结果与分子筛协会官方网站提供的手工检索结果相对比,以此检验我们这种检索方法的准确性.在此基础上,我们又任意定义了一些新的构筑基元,进一步考察这种结构检索方法在不同条件下的效率.测试的硬件环境是Intel Core2 Duo T9600 2.8 GHz,4 GB DDR2 800 MHz,操作系统为Windows 7 Ultimate SP1 32位,程序开发环境为ActivePerl 5.12.4.

表2 采用不同策略进行子结构检索所需时间的比较Table 2 Time needed for substructure search by different approaches

选取的已知分子筛构筑基元包括lov、nat、vsv、bre、bea等,22部分测试结果列在表2中.国际分子筛协会经过人工分析,在其官方数据库中给出了这些常见的构筑基元在197个分子筛结构中出现的情况.经过对比发现,采用我们的子结构检索方法得到的结果与国际分子筛协会人工统计的结果完全一致.说明我们的子结构检索方法在准确性上是可靠的.

此外,表2还列出了各种检索策略所需的CPU时间，包括单纯的VF2检索、引入对称性的VF2检索、采用编码预筛选的VF2检索以及同时引入对称性和编码预筛选的VF2检索.尽管这四种策略所得到的检索结果完全一致,但检索效率却相差较大.单纯的VF2算法耗时非常长,在由197个结构组成的数据库中进行子结构检索一般需要1 h以上的CPU时间.而在引入了对称性信息之后,检索消耗的时间会降低2-3个数量级.从理论上说,对称性的引入可以大大降低VF2算法所需遍历的节点个数,因此VF2算法的效率会得到指数级的提升.此外,编码预筛选能够在子图同构计算之前大幅度缩小目标结构的数量,同样也可以大幅度提高结构检索的效率.需要说明的是,待检的结构片段越特殊、越少见,编码预筛选提升的效率就越高,这个特性恰恰可以满足实际需求.因为在实际应用中,人们往往对新颖少见的子结构更关心,更需要在数据库中挖掘其结构特殊性.我们提出的采用编码预筛选的子图同构算法恰恰可以有效地解决这类问题.引入晶体结构对称性和编码预筛选这两套方案分别基于不同的原理,因此将这两种方法同时结合在一起,可以将子结构检索的效率进一步提升.从表2的最终结果可以看出,同时采用这两种改进方法的检索策略是效率最高的,对于简单的结构通常可以在1 s之内完成检索;对于相对复杂的结构,检索时间会略有增加,但仍在实际可接受的范围内.

4 结论

提出了一种针对无机晶体材料的子结构检索方法.该方法基于子图同构的VF2算法,利用晶体结构的对称性和编码预筛选策略降低结构检索的复杂度,可以在无机晶体数据库中快速准确地检索包含特定子结构的记录,为无机晶体材料的结构分析提供了重要的研究工具.本文以分子筛结构为例展示了这种结构检索方法的效率.除此之外,该方法同样可以用于其他无机晶体结构的检索,包括在假想分子筛结构数据库、23多配位态的开放骨架结构数据库、24甚至致密的矿物晶体结构数据库中进行的子结构检索.采用Perl语言实现了这种检索方法,其高度的可移植性使得这种检索方法可以在绝大多数计算机平台上得以实现.

(1) Rijnbeek,M.;Steinbeck,C.J.Cheminf.2009,1,17.

(2) Worlein,M.Extension and Parallelization of a Graph-Mining-Algorithm.Ph.D.Dissertation, Friedrich-Alexander-Universität Erlangen-Nürnberg,Germany, 2006.

(3) Meinl,T.;Borgelt,C.;Berthold,M.R.Discriminative Closed Fragment Mining and Perfect Extensions in MoFa.In Frontiers in Artificial Intelligence and Applications;Onaindia,E.,Staab, S.Eds.;IOS Press:Amsterdam,2004;Vol.109,pp 3-14.

(4) Huan,J.;Wang,W.;Prins,J.Efficient Mining of Frequent Subgraphs in the Presence of Isomorphism.In Proceedings of the 3rd IEEE International Conference on Data Mining, Melbourne,FL,USA,Nov 19-22,2003;Wu,X.D.,Tuzhilin, A.,Shavlik,J.Eds.;IEEE Computer Soc.:LosAlamitos,CA, 2003.

(5)Yan,X.F.;Han,J.W.gSpan:Graph-Based Substructure Pattern Mining.In Proceedings of 2002 IEEE International Conference on Data Mining,Maebashi City,Japan,Dec 9-12,2002; Kumark,V.Ed.;IEEE Computer Soc.:LosAlamitos,CA,2002.

(6) Nijssen,S.;Kok,J.N.Electronic Notes in Theoretical Computer Science 2005,127(1),77.

(7) Rahman,S.A.;Bashton,M.;Holliday,G.L.;Schrader,R.; Thornton,J.M.J.Cheminf.2009,1,12.

(8) Sastre,G.;Vidal-Moya,J.A.;Blasco,T.;Rius,J.;Jordá,J.L.; Navarro,M.T.;Rey,F.;Corma,A.Angew.Chem.Int.Edit. 2002,41,4722.

(9) Corma,A.;Rey,F.;Valencia,S.;Jordá,J.L.;Rius,J.Nature Mater.2003,2,493.

(10) Cambridge Structure Database.http://www.ccdc.cam.ac.uk/ products/csd/(accessedAug 30,2011).

(11) Bruno,I.J.;Cole,J.C.;Edgington,P.R.;Kessler,M.K.; MacRae,C.F.;McCabe,P.;Pearson,J.;Taylor,R.Acta Crystallogr.B:Struct.Sci.2002,58,389.

(12) Cook,S.A.The Complexity of Theorem-Proving Procedures.In Proceedings of the Third Annual ACM Symposium on the Theory of Computing,Ohio,USA,May 3-5,1971;Harrison, M.A.,Banerji,R.B.,Ullman,J.D.Eds.;ACM:New York, 1971.

(13) Ullmann,J.R.J.Assoc.Comput.Mach.1976,23,31.

(14) Schmidt,D.C.;Druffel,L.E.J.Assoc.Comput.Mach.1976, 23,433.

(15) McKay,B.D.Congressus Numerantium 1981,30,45.

(16) Cordella,L.P.;Foggia,P.;Sansone,C.;Vento,M.Performance Evaluation of the VF Graph MatchingAlgorithm.In Proceedings of the 10th International Conference on Image Analysis and Processing,Venice,Italy,Sept 27-29,1999; Roberto,G.,Cantoni,V.,Levialdi,S.Eds.;IEEE Computer Society Press:LosAlamitos,1999.

(17) Cordella,L.P.;Foggia,P.;Sansone,C.;Vento,M.IEEE Transactions on Pattern Analysis and Machine Intelligence 2004,26(10),1367.

(18) Su,Z.Q.;Liao,C.Z.;Xie,A.H.;Lu,X.P.;Shi,L.M. Computers and Applied Chemistry 2003,20(5),556.[苏振强,廖晨钟,谢爱华,鲁先平,石乐明.计算机与应用化学, 2003,20(5),556.]

(19) Li,X.;Song,T.T.;He,X.F.Computers and Applied Chemistry 2007,24(11),1551.[李 欣,宋婷婷,何险峰.计算机与应用化学,2007,24(11),1551.]

(20) Song,T.T.;He,X.F.;Wen,H.Computers and Applied Chemistry 2008,25(9),1152.[宋婷婷,何险峰,温 浩.计算机与应用化学,2008,25(9),1152.]

(21) Feng,H.J.;Wang,Y.;Zhou,J.L.;Haji,A.Computer Applications and Software 2010,27(10),117. [冯红君,汪 漪,周俊林,阿吉艾克拜尔·艾萨.计算机与应用软件, 2010,27(10),117.]

(22) The Database of Zeolite Structures.http://www.iza-structure.org (accessedAug 30,2011).

(23) Li,Y.;Yu,J.H.;Xu,R.R.Hypothetical Zeolite Database.http:// mezeopor.jlu.edu.cn/hypo/(accessedAug 30,2011).

(24) Li,Y.;Yu,J.H.;Xu,R.R.AlPO Database.http://mezeopor.jlu. edu.cn/alpo/(accessedAug 30,2011).

November 17,2011;Revised:December 27,2011;Published on Web:January 4,2012.∗

.Email:yili@jlu.edu.cn;Tel:+86-431-85168609.

A Computational Method for Specified Substructure Search in Inorganic Crystal Structures

HUO Wei-Feng1LI Yi1,*LU Jun-Ran1YU Ji-Hong1XU Ru-Ren1LI Jing2
(1State Key Laboratory of Inorganic Synthesis and Preparative Chemistry,Jilin University,Changchun 130012,P.R.China;2Department of Applied Mechanics,University of Science and Technology Beijing,Beijing 100083,P.R.China)

In this paper,a computational method for the substructure search in inorganic crystal structures is proposed.This method is based on the VF2 subgraph isomorphism algorithm.Furthermore, two additional approaches have been introduced into this method to improve the calculation efficiency of VF2:(1)introduction of crystal symmetry information with a view to avoiding redundant calculations among equivalent nodes(atoms);(2)a prescreening encoding treatment to enhance the calculation efficiency by greatly reducing the number of target structures.We tested the efficiency of this method by searching the zeolite crystal structure database from the International Zeolite Association for entries containing specified building units.The test results showed that this method could quickly and correctly retrieve all the entries containing the queried substructure in the zeolite structure database.The introduction of crystal symmetry information and the prescreening encoding treatment greatly reduce the complexity of substructure search. The search speed was significantly enhanced by at least 3-5 orders of magnitude.This method was developed using Perl programming language,ensuring that this method could be easily applied to various platforms.

Inorganic crystal structure;Database;Substructure search;VF2 algorithm;Zeolite; Building unit

10.3866/PKU.WHXB201201041

O641

The project was supported by the National Natural Science Foundation of China(21001049).

国家自然科学基金(21001049)资助项目